Correlated Rotational Alignment Spectroscopy

Visualizing Rotational Wave Packets (superposition of spherical harmonics)

Attached below is a Python script to display spherical harmonics and spherical harmonic superpositions. Spherical harmonics are the angular wavefunctions (“eigenfunctions”) used to describe rotational states of molecules and the angular properties of electrons in atoms (i.e., “atomic orbitals”). A single eigenfunction of the Schrödinger equation is time-independent but dynamics can be described by a superposition (sum) of eigenfunctions. E.g., the square of a sum of spherical harmonics predicts the angular orientation of a molecule as function of time.

To use the script, you need to have a functional version of Python running on your computer. The header describes how to set up the correct environment using Anaconda Python. Have fun and comment if you found any issues or want to share your experience.

# -------------------
# This Python code plots a superposition of spherical harmonics (angular
# wave functions), that correspond to the energy eigenstates of a rotating 
# molecule or to the angular wavefunction of the electron in the hydrogen atom.
# The square of the angular angular wavefunction describes the probability to 
# observe the particle at the corresponding angle.
# Note that each eigenfunction is a time-independent solution to the angular 
# Schrödinger equation. "Add" additional eigenfunctions with different J and m 
# quantum numbers s to create a state superposition ("wave packet"), which 
# describes time-dependent properties, such as pseudo-classical molecular rotation.
# To run Mayavi in Conda, it is necessary to create an environment for 
# Python 3.6 (or lower) including numpy, scipy and pyqtgraph. It is neccessary 
# to install VTK and Mayavi via PIP ("conda install" does not work).
# I use Anaconda Python [download from:]
# In an administrator Conda shell, run:
# > conda create -n mayavi python=3.6 pyqt=5
# > conda activate mayavi
# > conda install numpy scipy pyqtgraph ipython
# > pip install vtk
# > pip install mayavi
# It is required to activate the (mayavi) environment (line 2) each time when 
# opening a new shell. Return to your (base) shell with "conda activate base".
# To run the script, I suggest you open an ipython shell and copy/paste the 
# code in three batches (marked [1],[2], and [3]). Note that only wavefunctions 
# with quantum numbers m = -J...J are allowed. In a conda shell, type:
# > conda activate mayavi
# > ipython
# and start copy/pasting the code from below.
# Written by Thomas Schultz, Apr. 2021
# Feel free to use and modify the code according to the Creative Commons license 
# [see:]

# [1a] Import mathematical functions and define wave functions ---#
from numpy import sin, cos, pi, exp, mgrid
from scipy.special import sph_harm

# The rotational wavefunctions are the Spherical Harmonics. 
# For explicit ways to calculate the angular wavefunctions, refer to a textbook.
def TiWf (J,m,theta,phi):               # Time-independent angular wave function  
    return sph_harm(m, J, theta, phi)

def TdWf(k,t,phase):                    # Time dependent wave function. WF(r,t) = TiWf(r) * TdWf(t)
    return np.exp(-1j*k*(t+phase))      # k is the wavevector (angular frequency), == T(J,B)/2*pi

# To look at time-dependent wavefunctions, we need the state energies / terms
B = 3.27e9                              # Rotational constant for CS2 (in Hz) 
def EJ(J): return J*(J+1)               # Energy of a J state in a rigid linear rotor
def T(J,B): return B*J*(J+1)            # Energy term of J state (in freq. units)

# [1b] Plot a time-independent wave function             ---#
from mayavi import mlab
# Spatial coordinates for plotting:
# Create a grid over a sphere defined by r, phi, theta
# More grid points look better, but take more time to calculate
r = 0.3         # radius
grid = 80j
phi, theta = mgrid[0:pi:grid, 0:2*pi:2*grid]

# Spherical to Cartesian coordinate transformation
x = r * sin(phi) * cos(theta)
y = r * sin(phi) * sin(theta)
z = r * cos(phi)    

SqWF = abs(TiWf(2,1,theta,phi))**2      # Wave function psi-squared for state J=2, m=1
plot = mlab.mesh(SqWF *x, SqWF *y, SqWF *z, colormap='RdBu', opacity=1., transparent=True)

# [2a] Initialize a list of wavefunctions to plot a state superposition ---#

t = 0; deltaT = 5e-12                   # Initial time (0), time-increment (1 ps)
J = [0,]                                # List of J-values
m = [0,]                                # List of m-values
Amp = [1.,]                             # List of amplitudes
phase = np.array([0.,])                 # List of initial phases
showReal = [True,]                      # Show the real or imaginary part
showImag = [True,]                      # ... of the wavefunction
TiWfList = []                           # Initialize list of time-independent wavefunctions (don't recalculate for each step)
kList = []                              # Initialize list of k-vectors for propagation
for _J,_m,_a in zip(J,m,Amp):           # Calculate TiWf for each state in the list. (recalculate when J, m, Amp change)
    TiWfList.append(_a* TiWf(_J,_m,theta,phi))  # Time-independent wavefunctions WF(r)
    kList.append(T(_J,B))               # k-vectors for WF propagation: WF(r,t) = WF(r)*exp(-ikt)

# [2b] Animated plot: Update TiWf*TdWf for incrementing time values ---#

def anim():
    global TiWfList, kList, phase, t, deltaT, x,y,z
    while True:
        if showReal[0] and showImag[0]:
            sumWF = TiWfList[0]*TdWf(kList[0],t,phase[0])
          if showReal[0]:
            sumWF = (TiWfList[0]*TdWf(kList[0],t,phase[0])).real +(0+0j)
          if showImag[0]:
            sumWF = (TiWfList[0]*TdWf(kList[0],t,phase[0])).imag +(0+0j)
          if not (showReal[0]) and not (showImag[0]):
            sumWF = (TiWfList[0]*TdWf(kList[0],t,phase[0])).imag *(0+0j)
        if len(kList)>1: 
          for i in range(len(kList)-1):
            if (showReal[i+1] and showImag[i+1]):
              sumWF += TiWfList[i+1]*TdWf(kList[i+1],t,phase[i+1])
            elif showReal[i+1]:
              sumWF += ((TiWfList[i+1]*TdWf(kList[i+1],t,phase[i+1])).real +(0+0j))
            elif showImag[i+1]:
              sumWF += ((TiWfList[i+1]*TdWf(kList[i+1],t,phase[i+1])).imag +(0+0j))
        plot.mlab_source.x = np.absolute(sumWF)**2 *x       # replace x data
        plot.mlab_source.y = np.absolute(sumWF)**2 *y       # replace y data
        plot.mlab_source.z = np.absolute(sumWF)**2 *z       # replace z data
        t += deltaT                           # increment time 

# [3a] Create a user interface for easy control of the plot. ---#
#--- I use a PyQt ParameterTree as user interface            ---#
import pyqtgraph as pg
from pyqtgraph.Qt import QtCore, QtGui, QtWidgets
import pyqtgraph.parametertree.parameterTypes as pTypes
from pyqtgraph.parametertree import Parameter, ParameterTree, ParameterItem, registerParameterType
qt_app = QtWidgets.QApplication(sys.argv)

# List of parameters for each state
rotParameters = [
        {'name': 'amplitude', 'type': 'float', 'value': 1, 'limits': (0, 1), 'default': 1.},
        {'name': 'J', 'type': 'int', 'value': 0, 'limits': (0, 100)},
        {'name': 'm', 'type': 'int', 'value': 0, 'limits': (0, 100)},
        {'name': 'phase', 'type': 'float', 'value': 0., 'default': 0.},
        {'name': 'Show real**2', 'type': 'bool', 'value': True, 'tip': "add abs(real(phi)) to plot"}, 
        {'name': 'Show imag**2', 'type': 'bool', 'value': True, 'tip': "add abs(imag(phi)) to plot"}, 

# This class allows the user to add new parameters into the parameter tree
class ScalableGroup(pTypes.GroupParameter):
    """ Extendable group of parameters for parameter tree."""
    def __init__(self, **opts):
        opts['type'] = 'group'
        opts['addText'] = "Add"
        opts['addList'] = ['State']
        pTypes.GroupParameter.__init__(self, **opts)
        self.counter = 1
        self.addChild(dict(name="State %d" %(self.counter), children=rotParameters, removable=True, renamable=True))
    def addNew(self, typ):
        self.addChild(dict(name="State %d" %(self.counter+1), children=rotParameters, removable=True, renamable=True))
        self.counter += 1

# List of global property parameters 
parameters = [ 
    {'name': 'Global Properties', 'type': 'group', 'expanded': True, 'children': [
        {'name': 'Radius', 'type': 'float', 'value': 0.3, 'limits': (0.001, 10), 'default': 0.3, 'tip': "radius of image"},
#        {'name': 'Animate', 'type': 'bool', 'value': False, 'tip': "animate plot"},
        {'name': 'Speed (ps)', 'type': 'float', 'value': 10, 'limits': (0.01, 200), 'default': 0.5, 'tip': "time step for animation"},
    ScalableGroup(name="Rotational State Properties"),  # Call the ScalableGroup to add rotational parameters to the tree

# If anything changes in the parameter tree, update the wavefunction parameters
def change(parameterTree, changes):
    """ This function must handle all user requests by interpreting the 
        parameter origin, name, and value"""
    global B, Amp, J, m, phase, TiWfList, kList, showReal, showImag, theta, phi, deltaT, x,y,z
    for parameterName, change, data in changes:
        path = parameterTree.childPath(parameterName)
        print(f'path: {path}')
        if path is not None:
            childName = '.'.join(path)
            childName =
        if change == "childAdded": 
            print("== child added ==")            # New state added: Now update all wavefunctions
            Amp, J, m, phase, TiWfList, kList, showReal, showImag = updateWaveFunctions(parameterTree)  
        if change == "value":                     # Value change: Update only the affected wavefunction
          if path[0] == 'Rotational State Properties':
            stateNo = int(path[1].split()[1])-1
            parameter = path[2]
            print(f'state: {stateNo};   parameter: {parameter} = {data}')
            if parameter == 'amplitude': Amp[stateNo] = float(data)
            if parameter == 'J': J[stateNo] = int(data)
            if parameter == 'm': m[stateNo] = int(data)
            if parameter == 'phase': phase[stateNo] = float(data)
            if parameter == 'Show real**2': showReal[stateNo] = float(data)
            if parameter == 'Show imag**2': showImag[stateNo] = float(data)
            TiWfList[stateNo] = Amp[stateNo]* TiWf(J[stateNo],m[stateNo],theta,phi)
            kList[stateNo] = T(J[stateNo],B)
          elif path[0] == 'Global Properties':    # Global property change: update the affected parameters
            parameter = path[1]
            if parameter == 'Speed (ps)': deltaT = data*1e-12; 
            if parameter == 'Radius': 
              r = data
              x = r * sin(phi) * cos(theta)
              y = r * sin(phi) * sin(theta)
              z = r * cos(phi)

# This recalculates the whole list of wavefunctions. 
def updateWaveFunctions(p):
    global B, theta, phi
    Amp,J,m,phase,TiWfList,kList,showReal,showImag = [],[],[],[],[],[],[],[]
    values = p.getValues()["Rotational State Properties"][1]  # Table values for States
    for state in list(values):
        data = values[state][1]
        _a = float(data['amplitude'][0])
        _j = data['J'][0]
        _m = data['m'][0]
        _p = float(data['phase'][0]) 
        _real = data['Show real**2'][0]
        _imag = data['Show imag**2'][0]
        TiWfList.append(_a* TiWf(_j,_m,theta,phi))
    return Amp,J,m,phase,TiWfList,kList,showReal,showImag

#Now we actually create a window and show the parameter tree
win = QtGui.QWidget()
layout = QtGui.QGridLayout()
p = Parameter.create(name='params', type='group', children=parameters)
tree = ParameterTree()
tree.setParameters(p, showTop=False)
tree.setWindowTitle('Angular states')

CRASY: Observing Rotation in the Time Domain

Below, I posted a short video explaining how we observe rotation in the time domain. Please note that the video does not explain the quantum mechanical principles required to describe molecular rotation, but gives a pseudo-classical picture. This picture does not explain why molecular rotation is quantized or how we can relate rotational frequencies to molecular structure.

CRASY Data Analysis

Here is a guide for CRASY data analysis, written for a UNIST lab course: PC_lab_course_CRASY_Schultz.pdf. You can find an installation guide for the data analysis software on page 6 and a guide for CRASY data analysis on the following pages. If you want to learn more about the CRASY experiment, navigate to the landing page.

Required Software

You will require the following software:

(1) Python (I suggest to install the “Miniconda” package.)

(2) The crasyPlots program.

If you have trouble installing or running the software, watch the walk-through video below.


Data Analysis

The data is stored in compressed format (zip files). Here are a few data sets from 2020:
2-ns scan of carbon disulfide: 2020_CS2_Mar31_15.25_2ns
20-ns scan of pyridine: May25_17.31_20ns
100-ns scan of pyridine: 2020_pyridine_Mar26_16.00_100ns
500-ns scan of furan: 2020_furan_Jun02_10.45_500ns

Just follow the steps described in the guide. The pyridine (mass 79 u) and furan (mass 68 u) data also contains a small signal for CS2 (mass 76 u), which allows an easy calibration of the mass spectrum and also gives a nice frequency spectrum. If you have difficulties with the data analysis steps, you can watch the walk-through video below.

Listen to CRASY data

Our CRASY experiment measures rotational quantum waves of molecules. Usually we look at our data, but for the long traces we now measure this becomes a bit tedious and we really really have to zoom in to see anything. The example below is from a 17 nanosecond scan, but now we measure microsecond scans to get ever more precise data of our molecular structures.

How else could we explore our data?

Listen to it of course! So I slowed our data by a factor 50’000 to bring the Gigahertz frequencies into the audible range. Now you can hear our molecules rotate. The first example is our 300 nanosecond scan of the carbon disulfide rotational quantum-wave-packet. The molecule is linear and has a very simple harmonic spectrum.

So what do we actually hear? Like most things in the quantum world of atoms and molecules, the rotational motion is quantized and occurs only with discrete frequencies. We take snapshots of the molecular orientation at different points in time (in a mass spectrometer). Just like a rotating loudspeaker sounds louder when turned towards you, our signals are stronger when the molecules are oriented in a particular direction.

The hissing noise in the beginning is the signal from individual measurements. As we fill in more data points, the discrete molecular frequencies start to emerge. The spectrum is not quite harmonic because the molecules distort under the centripetal force of the rotation. This leads to the oscillating beating pattern at the end of the sound file.

To analyze the spectrum of frequencies, we can Fourier-transform the time-domain data and we get a spectrum as shown below. This data was published in PNAS.


Next, let’s listen to a 100 ns scan of 1,3-butadiene. This molecule is an asymmetric rotor and offers more complex harmonies. I’ll plot the corresponding frequency spectrum soon.

Finally the rather lengthy 1 microsecond scan of benzene. This molecule is a ‘symmetric rotor’ and has a much simpler spectrum than butadiene.

Below I plot some graphs from this 1 microsecond data-set. This scan was the longest time-domain measurement ever performed by us. The resolution is directly proportional to the scan length (Heisenberg’s time-energy uncertainty) and we obtained order-of-magnitude better resolution than any preceding Raman spectrum. The data is published in PCCP.

We actually observe our molecules in a mass spectrometer,  so first we look at the mass spectrum to identify the molecules. Mass 78 is the mass of benzene, so in the second row we plot only this signal amplitude, but now as function of delay time. This second trace nicely illustrates the trick we use to get to large delays without too much effort: we skip most point and only measure data at a random subset of delays. This explain the grainy noise in the beginning of all the sound files: We don’t have enough information to piece together the waves from the few measured data points.  The bottom graph shows the Fourier-transformed data from the middle graph. We can see the nice harmonic (evenly spaced) spectrum of benzene.

From the measured rotational frequencies, we determine the molecular shapes. But that story is better left for another post.


Scientific writing Style-Guide

The following are some resources that may help you with your scientific writing. Scientific writing can be a challenge even for native writers, so don’t despair. Practice makes perfect and you have to write to learn to write.

If you have to write a report, please use a template and submit your report with double line spacing for easier correction. I propose the template from JACS for MS Word or Latex.

To master structure and style, please refer to the relevant Style Guides. Here is a short 1-page cheat sheet to help you getting started: Writing_in_english_(1-page_for_koreans). Also look at the style guide of the American Institute of Physics for more extensive guidelines on the structure and style of scientific reports. You can find an excerpt from the AIP style guide at AIP_Style_4thed_extract or you can download the full guide from here. The ACS offers an exhaustive style guide at

Use dictionaries to find the right words (e.g.: (US-English), (British-English)). Use a spell checker to avoid unnecessary spelling mistakes.

If you are a non-native speaker, please make an extra effort to check your articles. Here is a very short primer about article placement in the English language:

  • “a” (also: “one”): Indefinite article. Always use when referring to one object that the reader does not yet know about.
  • “the” (also: “that”, “this”): Definite article. Always use when referring to one or multiple objects that the reader already knows about.
  • “” (no article): Plural indefinite article for objects that the reader did not yet know about, or for general statements about objects.

E.g.: I teach in a University. (You didn’t know about it yet.) The University is new. (I talk about the same University, so you already know about it.) I like a University. (Now I talk about a different University that you don’t know about.) Universities are cool. (I make a general statement about all Universities.)

And finally, read English books or English newspapers to improve your English style and vocabulary. There are lots of exciting English texts out there and reading is the easiest way to learn!

Here is some specific feedback for the 2016 Lab course reports: report_feedback.

IR Correlation Table and other Lab Course Material

  • The introduction slides for the lab course can be found on the UNIST Blackboard system.
  • Vibrational absorption frequencies for common chemical groups are listed in the IR_correlation_table.
  • Relevant vibrational Raman reference spectra [1] are summarized in this document: Literature Raman spectra.
  • You can ask questions to the teaching assistants 인호 ( and Begum Ozer (

[1] SDBSWeb : (National Institute of Advanced Industrial Science and Technology, Aug. 9, 2016)