Tag Archives: Quantum

Atomic Orbitals

Prerequiresites: Quantum Mechanics course

Electrons around a nucleus. Do they look like little well behaved planets orbiting a sun?

NOPE!

We get spread out blobs in special little patterns called orbitals. Here, we will look at their shapes and properties a bit. Today we will look at graphs in 1D and 2D, but the next post, Atomic Orbitals Pt. 2, uses a fancy, but slightly unstable plotting package, GLVisualize to generate some 3D plots.

The Hamiltonian for our problem is:

begin{equation}
{cal H}Psi(mathbf{x}) =left[ -frac{hbar}{2 m} nabla^2 – frac{Z e^2}{4 pi epsilon_0 r}right]Psi(mathbf{x}) = E Psi(mathbf{x})
end{equation}
with
begin{equation}
nabla^2= frac{1}{r^2}frac{partial}{partial r} left(
r^2 frac{partial}{partial r}
right)+
frac{1}{r^2 sin theta} frac{partial}{partial theta} left(
sin theta frac{partial}{partial theta}
right)+
frac{1}{r^2 sin^2 theta} frac{partial^2}{partial phi^2}
end{equation}

To solve this problem, we begin by guessing a solution with separated radial and angular variables,
begin{equation}
Psi(mathbf{x}) = R(r) Theta ( theta,phi)
end{equation}

begin{equation}
frac{E r^2 R(r)}{2r R^{prime}(r) + r^2 R^{prime prime}(r)}=
frac{left( frac{1}{sin theta} frac{partial}{partial theta} left(
sin theta frac{partial Theta(theta,phi)}{partial theta}
right)+
frac{1}{sin^2 theta} frac{partial^2 Theta(theta,phi)}{partial phi^2}right) }{Theta( theta, phi)}
=C
end{equation}

Instead of going into the precise mechanisms of solving those two separate equations here, trust for now that they follow standard special functions, the associated Legendre polynomial and the generalized Laguerre polynomial. Try a standard quantum mechanics textbook for more information about this.

begin{equation}
Y^m_l(θ,ϕ) = (-1)^m e^{i m phi} P^m_l (cos(θ))
end{equation}
where $P^m_l (cos (theta))$ is the associated Legendre polynomial.

begin{equation}
R^{n,l} ( rho ) = rho ^l e^{- rho /2} L^{2 l+1}_{n-l-1} ( rho )
end{equation}
where $L^{2 l+1}_{n-l-1}(rho)$ is the generalized Laguerre polynomial.

begin{equation}
rho=frac{2r}{n a_0}
end{equation}

begin{equation}
N=sqrt{left(frac{2}{n}right)^3 frac{(n-l-1)}{2n(n+l)!}}
end{equation}

#Pkg.update();
#Pkg.add("GSL");
#Pkg.add("PyPlot");
using GSL;    #GSL holds the special functions
using PyPlot;

Cell to Evaluate

What’s below is a bunch of definitions that makes our calculations easier later on. Here I utilize the GNU scientific library, GSL imported above, to calculate the special functions.

Programming Tip!

Even though it’s not necessary, specifying the type of inputs to a function through m::Int helps prevent improper inputs and allows the compiler to perform additional optimizations. Julia also implements abstract types, so we don’t have to specify the exact type of Int. Real allows a numerical, non-complex type.

Type Greek characters in Jupyter notebooks via LaTeX syntax, e.g. alpha+tab

The function Orbital throws DomainError() when l or m do not obey their bounds. Julia supports a wide variety of easy to use error messages.

a0=1; #for convenience, or 5.2917721092(17)×10−11 m

# The unitless radial coordinate
ρ(r,n)=2r/(n*a0);

#The θ dependence
function Pmlh(m::Int,l::Int,θ::Real)
    return (-1.0)^m *sf_legendre_Plm(l,m,cos(θ));
end

#The θ and ϕ dependence
function Yml(m::Int,l::Int,θ::Real,ϕ::Real)
    return  (-1.0)^m*sf_legendre_Plm(l,m,cos(θ))*e^(im*m*ϕ)
end

#The Radial dependence
function R(n::Int,l::Int,ρ::Real)
    if isapprox(ρ,0)
        ρ=.01
    end
     return sf_laguerre_n(n-l-1,2*l+1,ρ)*e^(-ρ/2)*ρ^l
end

#A normalization: This is dependent on the choice of polynomial representation
function norm(n::Int,l::Int)
    return sqrt((2/n)^3 * factorial(n-l-1)/(2n*factorial(n+l)))
end

#Generates an Orbital function of (r,θ,ϕ) for a specified n,l,m.
function Orbital(n::Int,l::Int,m::Int)
    if l>n    # we make sure l and m are within proper bounds
        throw(DomainError())
    end
    if abs(m)>l
        throw(DomainError())
    end
    psi(ρ,θ,ϕ)=norm(n, l)*R(n,l,ρ)*Yml(m,l,θ,ϕ);
    return psi
end

#We will calculate is spherical coordinates, but plot in Cartesian, so we need this array conversion
function SphtoCart(r::Array,θ::Array,ϕ::Array)
    x=r.*sin(θ).*cos(ϕ);
    y=r.*sin(θ).*sin(ϕ);
    z=r.*cos(θ);
    return x,y,z;
end

function CarttoSph(x::Array,y::Array,z::Array)
    r=sqrt(x.^2+y.^2+z.^2);
    θ=acos(z./r);
    ϕ=atan(y./x);
    return r,θ,ϕ;
end

"Defined Helper Functions"

Parameters

Grid parameters:
You might need to change rmax to be able to view higher $n$ orbitals.

Remember that
begin{equation}
0<n ;;;;; ;;;; 0 leq l < n ;;;;; ;;;; -l leq m leq l
;;;;; ;;;; n,l,m in {cal Z}
end{equation}

# Grid Parameters
rmin=.05
rmax=10
Nr=100 #Sampling frequency
Nθ=100
Nϕ=100

# Choose which Orbital to look at
n=3;
l=1;
m=0;
"Defined parameters"
#Linear Array of spherical coordinates
r=collect(linspace(rmin,rmax,Nr));
ϕ=collect(linspace(0,2π,Nθ));
θ=collect(linspace(0,π,Nϕ));
#3D arrays of spherical coordinates, in order r,θ,ϕ
ra=repeat(r,outer=[1,Nθ,Nϕ]);
θa=repeat(transpose(θ),outer=[Nr,1,Nϕ]);
ϕa=repeat(reshape(ϕ,1,1,Nϕ),outer=[Nr,Nθ,1]);

x,y,z=SphtoCart(ra,θa,ϕa);

Though I could create a wrapped up function with Orbital(n,l,m) and evaluate that at each point, the below evaluation takes advantage of the separability of the solution with respect to spherical dimensions. The special functions, especially for higher modes, take time to calculate, and the fewer calls to GSL, the faster the code will run. Therefore, this implementation copies over radial and angular responses.

Ψ=zeros(Float64,Nr,Nϕ,Nθ)
θd=Int64(round(Nθ/2))  ## gives approximately the equator.  Will be useful later

p1=Pmlh(m,l,θ[1]);
p2=exp(im*m*ϕ[1]);
for i in 1:Nr
    Ψ[i,1,1]=norm(n,l)*R(n,l,ρ(r[i],n))*p1*p2;
end

for j in 1:Nθ
    Ψ[:,j,1]=Ψ[:,1,1]*Pmlh(m,l,θ[j])/p1;
end

for k in 1:Nϕ
    Ψ[:,:,k]=Ψ[:,:,1]*exp(im*m*ϕ[k])/p2;
end
pygui(false)
xlabel("θ")
ylabel("Ψ")
title("Wavefunction for n= $n ,l= $l ,m= $m ")

annotate("l= $l Angular Node",
xy=[π/2;0],
xytext=[π/2+.1;.02],
xycoords="data",
arrowprops=Dict("facecolor"=>"black"))

plot(θ,zeros(θ))
plot(θ,reshape(Ψ[50,:,1],100)) #reshape makes Ψ 1D

2p Angle Slice

A slice along the θ plane showing an angular node for the 2p orbital.

pygui(false)
xlabel("r")
ylabel("Ψ")
title("Wavefunction for n= $n ,l= $l ,m= $m ")

plot(r,zeros(r))
plot(r,reshape(Ψ[:,50,1],100)) #reshape makes Ψ 1D

3p Radial Slice

A slice along the radial plane showing a radial node in the 3p orbital.

#rap=squeeze(ra[:,:,50],3)
#θap=squeeze(θa[:,:,50],3)
#ϕap=squeeze(ϕa[:,:,50],3)
#Ψp=squeeze(Ψ[:,:,50],3)
rap=ra[:,:,50]
θap=θa[:,:,50]
ϕap=ϕa[:,:,50]
Ψp=Ψ[:,:,50]
xp,yp,zp=SphtoCart(rap,θap,ϕap);
pygui(false)
xlabel("x")
ylabel("z")
title("ϕ-slice of Ψ for n=$n, l=$l, m=$m")
pcolor(xp[:,:],zp[:,:],Ψp[:,:],cmap="coolwarm")
colorbar()

3p in 2d

Slice of a 3p orbital in the x and z plane.

3dz2 in 2d

Slice of a 3dz2 orbital in the x and z plane.

Don’t forget to check out Atomic Orbitals Pt. 2!

Quantum Harmonic Oscillator

Prerequiresites: Quantum Mechanics course

Slinkies. They started out as toys. I still have one to play with on my desk.

Rubber bands What was once something useful, is now a wonderful projectile weapon.

Swings I still love them, but people seem to not make them in adult sizes for some reason.

A person’s perception of these objects starts to change as they enter their first physics class. Even in that beginning class of classical mechanics, the problems are filled with harmonic oscillators, like slinkies, rubber bands, or swings, which exert a force proportional to their displacement

and therefore a quadratic potential

This is all extremely fun and useful in the classical regime, but we add Quantum Mechanics to the mix, and LOW AND BEHOLD! we have one of the few exactly solvable models in Quantum Mechanics. Moreso, this solution demonstrates some extremely important properties of quantum mechanical systems.

The Hamiltonian

The Solution

Today, I just intend to present the form of the solution, calculate this equation numerically, and visualize the results. If you wish to know how the equation is derived, you can find a standard quantum mechanics textbook, or stay tuned till I manage to write it up.

Physicists’ Hermite Polynomials

Note: These are not the same as the “probabilists’ Hermite Polynomial”. The two functions differ by scaling factors.

Physicists’ Hermite polynomials are defined as eigenfunctions for the differential equation

I leave it as an exercise to the reader to:

  • demonstrate orthogonality with respect to the measure $e^{-x^2}$, i.e.

  • demonstrate completeness. This means we can describe every function by a linear combination of Hermite polynomials, provided it is suitably well behaved.

Though a formula exists or calculating a function at $n$ directly, the most efficient method at low $n$ for calculating polynomials relies on recurrence relationships. These recurrence relationships will also be quite handy if you ever need to show orthogonality, or expectation values.

Recurrence Relations

#Pkg.update();
#Pkg.add("PyPlot");
#Pkg.update()
#Pkg.add("Roots")
using Roots;
using PyPlot;

Programming Tip!

Since Hermite polynomials are generated recursively, I wanted to generate and save all the functions up to a designated value at once. In order to do so, I created an array, whose values are anonymous functions.

function GenerateHermite(n)
    Hermite=Function[]

    push!(Hermite,x->1);
    push!(Hermite,x->2*x);

    for ni in 3:n
        push!(Hermite,x->2.*x.*Hermite[ni-1](x).-2.*n.*Hermite[ni-2](x))
    end
    return Hermite
end

Let’s generate some Hermite polynomials and look at them.
Make sure you don’t call a Hermite you haven’t generated yet!

Hermite=GenerateHermite(5)

Programming Tip!

Since the Hermite polynomials, and the wavefunctions after them, are composed of anonymous functions, we need to use map(f,x) in order to map the function f onto the array x. Otherwise our polynomials only work on single values.

x=collect(-2:.01:2);
for j in 1:5
    plot(x,map(Hermite[j],x),label="H_$j (x)")
end
legend()
ylim(-50,50)

Hermite Polynomials

First few Hermite Polynomials

# Lets make our life easy and set all units to 1
m=1
ω=1
ħ=1

#Finally, we define Ψ
Ψ(n,x)=1/sqrt(factorial(n)*2^n)*(m*ω/(ħ*π))^(1/4)*exp(-m*ω*x^2/(2*ħ))*Hermite[n](sqrt(m*ω/ħ)*x)

Finding Zeros

The eigenvalue maps to the number of zeros in the wavefunction. Below, I use Julia’s roots package to indentify roots on the interval from -3 to 3.

zeds=Array{Array{Float64}}(1)
zeds[1]=[] #ground state has no zeros
for j in 2:4
    push!(zeds,fzeros(y->Ψ(j,y),-3,3))
end
# AHHHHH! So Much code!
# Don't worry; it's all just plotting
x=collect(-3:.01:3)  #Set some good axes

for j in 1:4    #how many do you want to view?
    plot(x,map(y->Ψ(j,y),x)+j-1,label="| $j >")
    plot(x,(j-1)*ones(x),color="black")
    scatter(zeds[j],(j-1)*ones(zeds[j]),marker="o",s=40)
end
plot(x,.5*m*ω^2*x.^2,linestyle="--",label="Potential")

scatter([],[],marker="o",s=40,label="Zeros")
xlabel("x")
ylabel("Ψ+n")
title("Eigenstates of a Harmonic Osscilator")
legend()
xlim(-3,3);
ylim(-.5,4.5);

Example Result

Eigenstates

Eigenstates of the Quantum Harmonic Osscillator

More to come

This barely scratched the surface into the richness that can be seen in the quantum harmonic oscillator. Here, we have developed a way for calculating the functions, and visualized the results. Stay tuned to hear about ground state energy, ladder operators, and atomic trapping.

Atomic Orbitals Pt. 2

Prerequiresites: Quantum Mechanics course

If you haven’t read it already, check out Atomic Orbitals Pt. 1. Today, we try and make some prettier pictures. GLVisualize is quite a beautiful package, but not entirely the easiest to use at this point with some not so consistent documentation.

To add this package:

Pkg.add("GLVisualize")

and test with:

Pkg.test("GLVisualize")

But, other steps may be necessary to get the package working. On a Mac, I was required to install the Homebrew.jl package.

#Pkg.update();
#Pkg.add("GSL");
using GSL;
using GLVisualize;
a0=1; #for convenience, or 5.2917721092(17)×10−11 m

# The unitless radial coordinate
ρ(r,n)=2r/(n*a0);

#The θ dependence
function Pmlh(m::Int,l::Int,θ::Real)
    return (-1.0)^m *sf_legendre_Plm(l,m,cos(θ));
end

#The θ and ϕ dependence
function Yml(m::Int,l::Int,θ::Real,ϕ::Real)
    return  (-1.0)^m*sf_legendre_Plm(l,abs(m),cos(θ))*e^(im*m*ϕ)
end

#The Radial dependence
function R(n::Int,l::Int,ρ::Real)
    if isapprox(ρ,0)
        ρ=.001
    end
     return sf_laguerre_n(n-l-1,2*l+1,ρ)*e^(-ρ/2)*ρ^l
end

#A normalization: This is dependent on the choice of polynomial representation
function norm(n::Int,l::Int)
    return sqrt((2/n)^3 * factorial(n-l-1)/(2n*factorial(n+l)))
end

#Generates an Orbital Funtion of (r,θ,ϕ) for a specificied n,l,m.
function Orbital(n::Int,l::Int,m::Int)
    if (l>n)    # we make sure l and m are within proper bounds
        throw(DomainError())
    end
    if abs(m)>l
       throw(DomainError())
    end
    psi(ρ,θ,ϕ)=norm(n, l)*R(n,l,ρ)*Yml(m,l,θ,ϕ);
    return psi
end

#We will calculate is spherical coordinates, but plot in cartesian, so we need this array conversion
function SphtoCart(r::Array,θ::Array,ϕ::Array)
    x=r.*sin(θ).*cos(ϕ);
    y=r.*sin(θ).*sin(ϕ);
    z=r.*cos(θ);
    return x,y,z;
end

function CarttoSph(x::Array,y::Array,z::Array)
    r=sqrt(x.^2+y.^2+z.^2);
    θ=acos(z./r);
    ϕ=atan(y./x);
    return r,θ,ϕ;
end

"Defined Helper Functions"

Here, create a square cube, and convert those positions over to spherical coordinates.

range=-10:.5:10
x=collect(range);
y=collect(range);
z=collect(range);
N=length(x);
xa=repeat(x,outer=[1,N,N]);
ya=repeat(transpose(y),outer=[N,1,N]);
za=repeat(reshape(z,1,1,N),outer=[N,N,1]);
println("created x,y,z")

r,θ, ϕ=CarttoSph(xa,ya,za);
println("created r,θ,ϕ")
Ψ=Orbital(3,2,-1)
Ψp=Orbital(3,1,0)
Ψv = zeros(Float32,N,N,N);
ϕv = zeros(Float32,N,N,N);
for nn in 1:N
    for jj in 1:N
        for kk in 1:N
            val=Ψ(ρ(r[nn,jj,kk],2),θ[nn,jj,kk],ϕ[nn,jj,kk]);
            #val+=Ψp(ρ(r[nn,jj,kk],2),θ[nn,jj,kk],ϕ[nn,jj,kk]);
            Ψv[nn,jj,kk]=convert(Float32,abs(val));
            ϕv[nn,jj,kk]=convert(Float32,angle(val));
        end
    end
end

mid=round(Int,(N-1)/2+1);
Ψv[mid,mid,:]=Ψv[mid+1,mid+1,:]; # the one at the center diverges
Ψv=(Ψv-minimum(Ψv))/(maximum(Ψv)-minimum(Ψv) );
w,r = glscreen()

robj=visualize(Ψv)

#choose this one for surfaces of constant of intensity
view(visualize(robj[:intensities],:iso))

#choose this for a block of 3D density
#view(visualize(Ψv))
r()

2p Orbital

2p

2p Orbital block showing the density of the wavefunction.

2p

2p Orbital shown via isosurface.

3d orbitals

3d0

3dz2 Orbital shown via isosurface. This corresponds to $n=3$, $l=2$, $m=0$.

3dm1

A 3d Orbital shown via isosurface. This corresponds to $n=3$, $l=2$, $m=-1$. This is not one of the canonical images, but instead an $m$ shape.

3dxy

3dxy (x2-y2) orbital shown in density. This is the sum of an $m=-2$ and $m=2$ state, for $n=3,l=2$.

3p

In order to get this 3p surface image to come out correctly, I used the square root of the values instead in order to be able to see the much fainter outer lobe.

3p

3p surface plot.

3p

3p space plot.