Today I was asked to calculate the Franck-Condon factors (FC) for nitrogen using the harmonic approximation. As I won’t go into physical details, here’s a two-sentence description: FC factors give a hint on relative intensities of radiation from a molecule and vary from 0 (no radiation) to 1 (quite high radiation). In general, they have rather little impact as selection rules are more important.
The bad news: the math is tricky. FC factors are the integral over the whole space of the product of the two transition states’ wave functions. The good news: as theoretical physics is the part of reality that can be approximated by an harmonic oscillator, math gets easier. Now the wave functions are real (no complex conjugation necessary) and only consist of Hermite polynomials with some prefactors and exponential sugar. For that case, python does the calculations (source code below). There are a few things to keep in mind:
- deltaX is the offset of the transitions target potential, so x=0 means the first potential’s parabola tip.
- The boundaries for the numerical integration are important. If you change the parameters in the order of magnitude, you might have to adjust the range. In general it should be from minus infinity to plus infinity. However, the wave functions’ product is zero for most of R. Hence, restricting is necessary. Obviously, in case the results depend on your integration boundaries, the boundaries are wrong 😉
#!/usr/bin/env python
# calculates the Franck-Condon-Factors using harmonic approximation for electronic transitions
# configuration
deltaX = 0.2e-10 # Distance between the two potentials in SI units
omega1 = 4.33e14 # angular frequency of the lower potential in SI units
omega2 = 4.14e14 # angular frequency of the upper potential in SI units
m = 1.66e-27*14 # vibrating masses
# -- code --
from scipy import integrate
from numpy.polynomial.hermite import hermval
import numpy
import math
hbar = 1.05e-34
# Hermite polynomial of n-th grade at x
def Hermite(n, x):
coeff = [0]*(n)
coeff.append(1)
return hermval(x, coeff)
# wavefunction of the harmonic oscillator
def psi(n, x, omega):
a = m*omega/hbar
return pow(a/math.pi, 0.25)/math.sqrt(pow(2, n)*math.factorial(n))* \
Hermite(n, math.sqrt(a)*x)*math.exp(-0.5*a*x*x)
# overlap product for transition n1 -> n2 at x0, where x0=(x=0) for the lower potential
def transition(x0, n1, n2):
return psi(n1, x0, omega1)*psi(n2, x0-deltaX, omega2)
#for i in range(-500, 500):
# print psi(2, float(i)/10e11, omega1)
print 'Franck-Condon-Factors for any transitions n1[0:9] to n2[0:9] with a value of more than 10e-3'
print '{0:3}{1:3}{2}'.format('n1', 'n2', 'FC factor')
for n1 in range(0, 10):
for n2 in range(0, 10):
# the integration bounds are to be determined by hand usually, they go from
# -numpy.inf to numpy.inf - biggest influence is the mass m
fc = integrate.quad(transition, -5e-9, 5e-9, args=(n1, n2))
fc = fc[0]*fc[0]
if fc > 10e-3:
print '{0:^3}{1:^3}{2}'.format(n1, n2, fc)