ORIGINAL ARTICLE: Jing B. Wang, Paul C. Abbott, Jim F. Williams; Visualizing atomic orbitals. Comput. Phys. 1 January 1996; 10 (1): 69–82. DOI: https://doi.org/10.1063/1.168561
Article Abstract: The Mathematica software package provides a range of tools for working with atomic systems. The symbolic tools include orthogonal polynomials and Clebsch-Gordan coefficients whilst the graphical capabilities cover polar plots; spherical plots; density plots; contour plots in two and three dimensions; and animation. These tools are applied to the manipulation and visualization of atomic orbitals.

INTRODUCTION

The probability distribution of the electron charge cloud in an atom, namely the atomic orbital, often reveals crucial information in atomic physics. Visual representation of such orbitals is invaluable to researchers working in this field and also helps teaching atomic structures to physics or chemistry students. The atomic orbital concepts were discussed in detail by Berry and the cited papers therein, including several visualization devices. However, the newly developed symbolic package Mathematica provides a wider range of graphical tools (including polar plots, spherical plots, density plots, contour plots in two and three dimensions, and animation), which can be applied to help visualise the radial and angular behaviour of multi-dimensional systems.
Mathematica also has a wide range of built-in special functions appropriate for the analysis of quantum systems. Orthogonal polynomials — Legendre, Laguerre, Chebyshev, Hermite and Jacobi (corresponding to the Mathematica functions
P
n
(x)
,
a
L
n
(x)
,
U
n
(x)
,
T
n
(x)
,
H
n
(x)
, and
(a,b)
P
n
(x)
) — arise when solving a wide class of eigenvalue and boundary value problems. The spherical harmonics,
m
Y
l
(θ,ϕ)
, are eigenfunctions of the angular part of the Laplacian operator (which is the quantum mechanical angular momentum operator).
Clebsch-Gordan coefficients
arise in the study of angular momenta in quantum mechanics, and in other applications of the rotation group. These functions, taken together with the orthogonal polynomials and graphical tools, provide a powerful environment for the study of atomic systems. It is worth noting that these methods can be readily utilized in the studies of other multi-dimensional systems.
As this article is not intended to be a tutorial on Mathematica, commands and syntax will not usually be explained. However, it is hoped that by showing a well-chosen set of examples, even readers unfamiliar with Mathematica will be able to follow the examples. The graphical output is placed immediately after its corresponding input command as it appears in Mathematica notebooks. Arbitrary units are used in these graphs.

1
. NOTATION AND PREREQUISITES

The quantity
Ω(r)r=

ψ
n,l,m
(r)r
pure basis state
ψ
n,l,m
(r),

ψ
n,l
(r)r
coherently mixed state
ψ
n,l
(r)=
∑
m
c
m
ψ
n,l,m
(r),
where
ψ
n,l,m
(r)≡nlm=
R
n,l
(r)
m
Y
l
(θ,ϕ),r=(r,θ,ϕ),r=
2
r
rΩ,Ω=sin(θ)θϕ,
represents the probability of finding an electron in the volume element
r
. The spherical harmonics,
m
Y
l
(θ,ϕ)
, give the angular distribution of the atomic orbitals and
R
n,l
(r)
is the corresponding radial distribution.
In the following, the general spherical harmonic will be denoted using the Mathematica shorthand
m
Y
l
(θ,ϕ)
. For example,
1
Y
2
(θ,ϕ)
-
1
2
15
2π
ϕ

sin(θ)cos(θ)
Although it is a built-in operation, complex conjugation, required to yield real-valued charge densities, is more easily achieved using a replacement rule:
In[]:=
*
exp_
:=exp/.Complex[a_,b_]Complex[a,-b]
For example,
*
1
Y
2
(θ,ϕ)
Out[]=
-
1
2
15
2π
-ϕ

sin(θ)cos(θ)
We make use of the Notation package to define custom input and output notations.
In[]:=
<<Notation`
In[]:=
Symbolize[
′
j
]

2
. ANGULAR DISTRIBUTION OF ATOMIC ORBITALS

1
. Pure states
ψ
nlm
(r)

For pure basis states,
ψ
n,l,m
(r)
, the angular distribution:
ang
Ω
l,m
(θ)=
2

m
Y
l
(θ,ϕ)
,
is axially symmetric with respect to the azimuthal angle
ϕ
. In this case, PolarPlot can provide a complete view of the atomic orbitals.
After defining the angular distribution function:
In[]:=
Ω
l_,m_
(θ_):=
m
Y
l
(θ,ϕ)
*
m
Y
l
(θ,ϕ)
individual charge densities can be immediately examined, e.g.,
Ω
3,2
(θ)
Out[]=
105
4
sin
(θ)
2
cos
(θ)
32π
Here are plots of the atomic orbitals for the states
l=1,2,3
,
m=0,…,l
.
GraphicsColumn@TablePolarPlotEvaluateTable
Ω
l,m
(θ),{m,0,l},{θ,0,2π},PlotRangeAll,{l,1,3}
Out[]=
In[]:=
ManipulatePolarPlotEvaluate
Ω
l,m
(θ),{θ,0,2π},AspectRatioAutomatic,PlotRange{{-0.8,0.8},{-0.3,0.3}},{l,Range[4]},{{m,0},Range[0,l],ControlTypeSetterBar},SaveDefinitions->True

2
. Mixed states
∑
m
c
m
ψ
nlm
(r)

For coherently mixed states
ψ
nl
(r)=
∑
m
c
m
ψ
nlm
(r)=r
R
nl
(r)
∑
m
c
m
m
Y
l
(θ,ϕ),
the angular probability distribution is
ang
Ω
l
(θ,ϕ)=
2

∑
m
c
m
m
Y
l
(θ,ϕ)
.
In this case, the charge cloud is no longer axially symmetric and its distribution is a function of both
θ
and
ϕ
. Three-dimensional plots are thus required to depict the structure of the coherently mixed atomic orbitals.
In general, the expression
ω
l
(θ,ϕ)=
l
∑
m=-l,2
a
m

b
m

m
Y
l
(θ,ϕ),
where
a
m
(
b
m
) are the amplitude ratio (relative phase) of each state, can be evaluated in Mathematica using:
In[]:=
ω
l_
[c_][θ_,ϕ_]:=
ω
l
[c][θ,ϕ]=
l
∑
i=0
c
〚i+1,1〛

c
〚i+1,2〛

l-2i
Y
l
(θ,ϕ)
The argument to
c
needs to be a list of
{
a
m
,
b
m
}
pairs. The general angular charge distribution then reads:
In[]:=
Ω
l_
[c_][θ_,ϕ_]:=Factor[ComplexExpand[Expand[
ω
l
(c)(θ,ϕ)
*
ω
l
(c)(θ,ϕ)
]]]
where the core of this expression is the product of
ω
l
(c)(θ,ϕ)
with its complex conjugate,
ω
l
(c)
*
(θ,ϕ)
. The other operations, Expand, ComplexExpand, and Factor, are used to simplify the expression to a human readable form.
As an example, for the general
P
-state wavefunction can be obtained from the mixed state consisting of
l=1
,
m=±1
:
ψ
P
(a_,b_)=
Ω
1

1
0
a
b
[θ,ϕ]
Out[]=
3
2
sin
(θ)(
2
a
-2acos(b-2ϕ)+1)
8π
and then plot the orbitals for the states with parameters
{a=0.1,b=0}
:
In[]:=
Manipulate[SphericalPlot3D[Evaluate[
ψ
P
(a,b)],{θ,0,π},{ϕ,0,2π},BoxRatios1,PlotRange0.5,MeshFalse,AxesFalse,BoxedFalse,SphericalRegionTrue],{{a,0.8},0,1},{{b,π/2},0,π},SaveDefinitions->True]
Note that the full range of ϕ was not shown to highlight the structure of the orbital.
Choosing the parameters arbitrarily as
one can produce a three-dimensional graph:
We clear the values assigned in this section.
Due to the interaction between the intrinsic magnetic dipole moment of the electron (spin) and the internal magnetic field of the atom (orbital angular momentum), the atomic charge density evolves with time. Animation provides views of such dynamic effects. Taking the spin-orbit (L-S) coupling into account, we obtained the time-dependent charge density function for a mixed state:
where
implemented in Mathematica as
are the so-called state multipoles and the corresponding time evolution factor is
implemented in Mathematica as
The time-dependent charge density function becomes
This can be plotted using the commands,
This command displays a sequence of graphics objects in succession and thus simulates the time evolution of the orbitals. Alternatively, here is an array of animation frames:
We clear the values assigned in this section.
For the hydrogen atom the radial wavefunction in atomic units is
where
The radial wavefunction is easily defined:
The radial distribution function,
then becomes
We have so far plotted separately the angular and the radial distribution functions of various atomic orbitals. However, these pictures are incomplete on their own and one needs to imaginatively assemble the angular and the radial plots in order to obtain a complete picture of the charge distribution in three-dimensional space.
Mathematica's DensityPlot is then adequate for showing the detailed distribution of the two-variable function, where the lighter region indicates higher charge density:
Defining a set of rules to transform from spherical polar to Cartesian coordinates:
the distribution in Cartesian coordinates is easily obtained:
Finally, for a general coherent mixture of several magnetic states
the total charge density distribution is a function of three variables (r,θ,ϕ):
can be plotted using ContourPlot3D. Here are the surfaces for the contour levels 0.0003, 0.0015, and 0.003.
and for the molecular charge density of the bonding,
and antibonding orbitals,
To see these orbitals in three dimensions, we also plot their contour surface at an arbitrary level of 0.005:
This article has demonstrated the capabilities of Mathematica for analyzing and visualizing atomic systems. The interplay between symbolic and graphical methods yields a tool that is more powerful than the sum of its parts and one that can be superior to special purpose visualization tools because of its flexibility.
Even for readers that have never used a computer algebra system, an examination of some of the code fragments presented above — say those implementing the radial wavefunctions for hydrogenic atoms or the triple integral of spherical harmonics — illustrates that implementing some quite complicated formulae in Mathematica is straightforward and transparent and, hence, less likely to contain or conceal an error. One of the goals of object oriented programming is to implement high-level structures. All of the computations in this article could, of course, be done in FORTRAN or C (or C++). However, in general, the physicist or chemist does not want to write packages for computing orthogonal polynomials, spherical harmonics, Clebsch-Gordan coefficients or basic graphical tools. Of course, there exists libraries of public domain and commercial software which, when linked together, can provide equivalent capabilities. But, since all these capabilities, along with general symbolic manipulation, are built-in to Mathematica it can be seen to be an ideal tool for a wide range of analysis of atomic and molecular systems.

REREFENCES

1. R.S. Berry, J. Chem. Ed. 43, 283 (1966)
2. S. Wolfram, Mathematica: A System for Doing Mathematics by Computer (Addison-Wesley, Reading, 1991)
3. A.R. Edmonds, Angular Momentum in Quantum Mechanics (Princeton, New Jersey, 1974)
4. B.H. Bransden and C.J. Joachain, Physics of Atoms and Molecules (Longman, UK, 1989)
5. K. Blum, Density Matrix Theory and Applications (Plenum: New York, 1981)
6. M. Abramowitz and I.A. Stegun, Handbook of Mathematical Functions (Dover, New York, 1972)
7. J.F. Shriner and W.J. Thompson, Computers in Physics, 7, 144 (1993)

CITE THIS NOTEBOOK

Visualizing atomic orbitals​
by Jingbo Wang, Paul C. Abbott, and Jim F. Williams
Wolfram Community, STAFF PICKS, October 10, 2023
​https://community.wolfram.com/groups/-/m/t/3031372