Atoms in external electric fields
Contents
1. Atoms in external electric fields#
Name:
Total Points: –/100 pts
Extra Credit: –/60 pts
1.1. Atomic polarizability#
A neutral atom contains an equal number of positive and negative charges. The positve charges are protons in the nucleus, and the negative charges are electrons orbiting the nucleus. In the absence of any external electromagnetic forces, the electrons orbit the nucleus in such a way that the center of charge of the electrons is at the same location as the center of charge of the nucleus. In reality, quantum mechanics tells us that electrons and protons do not have fixed positions, but consist of quantum mechanical wavefunctions. Protons, because of their higher mass, are much more localized than electrons, which form a “cloud” around the nucleus.
When an external field is applied, the electrons and protons are pushed in opposite directions, leading to a dipole moment:
The dipole moment can be expressed mathematically as \begin{align} \mathbf{p}=qr \end{align}
where \(q\) is the electric charge and \(r\) is the shifted position of the electron wave function relative to the nucleus. An extremely important quantity in nonlinear optics is the dipole moment per unit volume, or polarization \(\mathbf{P}\) of the material. If a material has a high polarization, it has many atoms with dipole moments close to one another.
In linear optics, the relative positions of the nucleus and the electron \(r\) varies linearly with the strength of the electric field. If you apply a field twice as strong, the center of the electronic wavefunction moves twice as far away from the nucleus. This leads to a linear relationship between the dipole moment and the applied electric field:
\begin{align} \mathbf{p}=\alpha\mathbf{E}_{ext} \end{align}
where \(\alpha\) is know as the linear polarizability of the atom. In nonlinear optics, we study situations in which the relative positions of the nucleus and the electrons vary nonlinearly with the strength of the applied field.
1.2. Magnitude of displacement: Linear Approximation#
So how large is the displacement of the nucleus relative to the electric cloud when an electric field is applied? Let’s make an approximation to get an estimate. If we assume a spherically uniform distribution of total electron charge \(q\) over a radius of \(R\), the charge density will be equal to
This electron charge distribution will create an electric field inside its interior equal to (using Gauss’s law):
where \(r\) is the position relative to the center of the sphere. Notice that this electric field is linear in \(r\). That means that the nucleus will feel a restoring force proportional to its displacement, just like a mass on a spring.
Suppose that as a result of the external field, the nucleus moves by an amount \(d\) with respect to the center of the electron cloud. The nucleus will feel a restoring force from the electric force from the cloud equal to
In equlibrium, this force \(\mathbf{F}_{cloud}\) is equal to the force applied by the external field \(q\mathbf{E}_{ext}\). So we have
and therefore
and the electronic polarizability is
Real atoms of course do not have perfectly spherical charge distributions, but this approximation yields polarizabilities that are quite close to measured values for many atoms. Using this approximation, let’s calculate the linear polarizability of an oxygen atom and get an order of magnitude estimate of the displacement of the electron cloud when an external field is applied:
### UNCOMMENT AND RUN THIS CELL IF USING GOOGLE COLAB
# !pip install ipympl -q
# from google.colab import output
# output.enable_custom_widget_manager()
#Import python packages we'll use
import numpy as np #import NumPy package. Did you know Travis Oliphant, former BYU ECE faculty, wrote NumPy?
from scipy.constants import * #get physical constants like pi and c
import matplotlib as mpl #plotting package
import matplotlib.pyplot as plt
#mpl.rc('text', usetex = True)
%matplotlib widget
#mpl.rcParams['figure.figsize'] = (10,5)
#Calculate electronic polarizibility using spherical electron shell approximation.
E = 1e4 # [V/m] external electric field magnitude
Z_Ox = 8 # oxygen atomic number
R_Ox = 100e-12 # effective atomic radius of Ox atom. Actual size I looked up is 48 pm.
d = 4*pi*epsilon_0*R_Ox**3/(e*Z_Ox)*E
alpha_calc = 4*pi*epsilon_0*R_Ox**3 #calculated polarizability
alpha_calc_cgs = R_Ox**3*1e6*1e24 #change to cgs units (most common unit in tables)
alpha_Ox_measured_cgs = 1.31 # 1.31 A^3 in cgs units (value I looked up for Oxygen)
print('Calculated polarizability: {:4.2f}, Measured polarizability: {:4.2f}'.format(alpha_calc_cgs, alpha_Ox_measured_cgs ))
print('Displacement with external electric field of E = {:.2e} V/m is d = {:.2e} m'.format(E, d))
Calculated polarizability: 1.00, Measured polarizability: 1.31
Displacement with external electric field of E = 1.00e+04 V/m is d = 8.68e-19 m
So with an applied field of \(10^4\) V/m in our approximate linear model, the electron cloud (uniform sphere!) is only displaced about \(10^{-18}\) m. Compared to the atomic radius of approximately \(10^{-10}\) m, this dispacement is very small. Like ridiculously small, which is why early researchers (pre-1950’s) were not very optimistic about nonlinear optics. The invention of the laser changed all that…but more on that later.
1.3. Exercises#
(30 pts) Calculate the electronic polarizibility of Hydrogen using the spherical, uniform electron cloud approximation. Compare with measured values. You’ll have to look them up using Google and/or the database of your choice.
Unlike our approximate model, real atoms don’t have a perfectly spherical charge distribution. Hence, we’d expect some nonlinear displacement if we apply a high enough external field. Let’s do a calculation using the actual shape of the electron cloud for the simplest atom: Hydrogen. It turns out that hydrogen gives us a pretty good estimate to how many materials behave.
1.4. Magnitude of displacement: Hydrogen Atom#
The actual radial charge distribution \(\rho(r)\) of an electron is proportional to the square of the quantum mechanical wavefunction describing the charge. For an electron in the Hydrogen 1s state, the wavefunction \(\psi(r)\) and its square \(|\psi(r)|^2\) are:
where \(e\) is the electron charge and \(a_0 = 53.92\) pm is the Bohr radius.
The probability \(P(r)\) of finding the electron in a shell of thickness \(dr\) at a distance \(r\) from the nucleus is given by
Let’s plot this:
#set up equations for quantum mechanical electron distribution
a0 = physical_constants["Bohr radius"][0] #get Bohr radius from scipy
r0 = 1e-5*a0 #set min radial coordinate. A little more than 0 to avoid dividing by zero
rf = 6*a0 #set max radial coordinate
N=1000 #number of points in r array
dr = (rf - r0)/(N) #dr step size
r = np.arange(r0,rf,dr) #radial coordinate array
Z = 1 # atomic number
psi = 1/np.sqrt(pi)*(Z/a0)**(3/2)*np.exp(-Z*np.abs(r)/a0) #electron density wavefunction
P = 4*pi*r**2*np.abs(psi)**2*dr # radial probability distribution
# plot quantum mechanical electron distribution
fig, ax = plt.subplots()
ax.plot(r/1e-12,psi**2/np.max(psi**2),label=r'$|\psi(r)|^2$', color='blue') #plot electron density as a function of radial coordinate
ax.plot(r/1e-12,P/np.max(P), label = r'$P(r)dr$', color='red') #plot radial distribution function
ax.vlines(a0/1e-12, 0, 1.1, colors='k', linestyles='dashed', linewidth=1, label=r'$a_0$ = {:.2f} pm'.format(a0/1e-12))
plt.title(r'Hydrogen 1s radial electron density')
plt.xlabel('Distance [pm]')
plt.legend()
plt.show()
Notice that the most likely place to find the electron is the Bohr radius \(a_0\) at about 53 pm from the nucleus. If we like, we can visualize the charge and probability distribution in 2D as an electron cloud:
# 2D visualization of actual quantum mechanical electron charge distribution
a0=1 #let's use distance units of a_0 for the plot
r0 = -3*a0 #set min radial coordinate
rf = 3*a0 #set max radial coordinate
# function to calculate electron wavefunction
def psi(x, y):
r = np.sqrt(x**2 + y**2)
return 1/np.sqrt(pi)*(Z/a0)**(3/2)*np.exp(-Z*r/a0)
#function to calculate radial probability distribution
def rad_dist(x, y):
r = np.sqrt(x**2 + y**2)
return 4*pi*r**2*np.abs(psi(x,y))**2
x = np.linspace(r0, rf, 100)
y = np.linspace(r0, rf, 100)
X, Y = np.meshgrid(x, y)
Zpsi = psi(X, Y)
ZD = rad_dist(X, Y)
cmap1 = plt.cm.Blues
cmap2 = plt.cm.Reds
fig, [ax1,ax2] = plt.subplots(nrows=1, ncols=2)
impsi = ax1.imshow(Zpsi, interpolation='bilinear', cmap=cmap1,
origin='lower', extent=[r0,rf,r0,rf],
vmax=abs(Zpsi).max(), vmin=0)
ax1.set_xlabel(r'distance in units of $a_0$')
ax1.set_title(r'charge distribution $\rho(r)$')
imD = ax2.imshow(ZD, interpolation='bilinear', cmap=cmap2,
origin='lower', extent=[r0,rf,r0,rf],
vmax=abs(ZD).max(), vmin=0)
ax2.set_xlabel(r'distance in units of $a_0$')
ax2.set_title(r'radial probability distribution $4\pi r^2 \rho(r)$')
plt.show()
1.5. Exercises#
(20 pts) What fraction of the electron charge is within a Bohr radius in Hydrogen? Does this mean our approximate model is really bad? Why or why not?
(20 pts) What fraction of the electron charge is inside the nucleus? Assume a nuclear radius is the diameter of a proton. Is physics broken if we just showed that the electron has some probability of being found inside the nucleus?
In order to calculate the restoring force, and hence potential energy of the electron as a function of position, we need to find the electric field resulting from this charge distribution. To do this we again use Gauss’s law, just like in the spherical approximation:
Since the 1s state in Hydrogen is radially symmetric, we may find the electric field at any distance \(r\):
That integral wasn’t too bad, but we should note that we have Python at our disposal (the whole point of doing this course in Python!), so we can just let Python do the work and numerically integrate. Let’s do that and compare to the analytic expressions for the enclosed charge and electric field:
# Write down expressions for enclosed charge as a function of radius and resulting electric field.
a0 = physical_constants["Bohr radius"][0] #get Bohr radius from scipy
r0 = 1e-6*a0 #set min radial coordinate
rf = 5*a0 #set max radial coordinate
N=1000 #number of points in r array
dr = (rf - r0)/(N) #dr step size
r = np.arange(r0,rf,dr) #radial coordinate array
Z = 1 # atomic number
psi = 1/np.sqrt(pi)*(Z/a0)**(3/2)*np.exp(-Z*np.abs(r)/a0) #electron density wavefunction
P = 4*pi*r**2*np.abs(psi)**2 # radial distribution function
#Numerical calcuations
qenc = e*np.cumsum(P)*dr #Calculate charge enclosed at radius r using a cumulative sum. This is a simple way to integrate in Python, but is not quite as accurate.
A = 4*pi*r**2 #Area of shell at radius r
E = qenc/epsilon_0/A #Electric field at radius r
#Analytic expressions
qenc_analytic = e*(-np.exp(-(2*r)/a0)*(a0**2 + 2*a0*r + 2*r**2)/a0**2)+e
E_analytic = qenc_analytic/epsilon_0/A
# plot electric field and enclosed charge
fig,ax = plt.subplots()
ax.plot(r/1e-12,qenc/np.max(qenc), label = r'$Q_{enc}(r)$ numeric', color='blue') #numerically calculated charge enclosed
ax.plot(r/1e-12,qenc_analytic/np.max(qenc_analytic), label = r'$Q_{enc}(r)$ analytic', color='k', linestyle = '--') #analytic charge enclosed
ax.plot(r/1e-12,E/np.max(E), label = r'$E(r)$ numeric', color='red') #plot Electric Field
ax.plot(r/1e-12,E_analytic/np.max(E_analytic), label = r'$E(r)$ analytic', color='k',linestyle = '--') #plot Electric Field
ax.vlines(a0/1e-12, 0, 1.1, colors='k', linestyles='dashed', linewidth=1, label=r'$a_0$ = {:.2f} pm'.format(a0/1e-12))
plt.title(r'Hydrogen 1s enclosed charge and resulting E-field')
plt.xlabel('Distance [pm]')
plt.legend()
plt.show()
The numeric and analytic expressions match very well, except perhaps right near the origin. Let’s do some slightly fancier integration using built-in Python integration tools from SciPy to see if we can do better. While we’re at it, let’s also do some fancy interactive plotting. This will be a useful Python skill for visualizing things later in the course. For this example we’ll vary the atomic number Z in the interactive plot to see how things change.
# Use Scipy integrate function to calculate electric field. Also demonstrate interactive plotting.
import ipywidgets as ipw
from scipy import integrate
def ode(rho, r, Z):
drho = 4*pi*r**2*np.abs(1/np.sqrt(pi)*(Z/a0)**(3/2)*np.exp(-Z*r/a0))**2
return drho
def update(Z = 1):
"""
update solution
"""
sol = e*integrate.odeint(ode, rho0, r, args = (Z,))
Esol = sol[:,0]/epsilon_0/A
line0.set_ydata(sol/e)
line1.set_ydata(Esol/np.max(Esol))
rho0 = 0
#r = np.logspace(np.log10(r0), np.log10(rf),N) #radial coordinate array
#A = 4*pi*r**2
dummy = np.zeros_like(r)
fig = plt.figure()
line0, = plt.plot(r/1e-12, dummy, label = "Enclosed charge")
line1, = plt.plot(r/1e-12, dummy, label = "E Numerical")
line3, = plt.plot(r/1e-12, E_analytic/np.max(E_analytic), linestyle = '--',label = "E Analytic, Z=1")
plt.grid()
plt.ylim(-0.1, 1.1)
plt.xlabel("$r$ [pm]")
#plt.ylabel("Enclosed Charge [e]")
plt.legend()
ipw.interact(update, Z = (1, 10, 1));
1.6. Exercises#
(30 pts) Add to the interactive plot so that it displays the value of the enclosed charge and the field at the Bohr radius.
Now we’ve plotted the electric field as a function of distance from the nuclues. It is very nonlinear! But as long as the electron doesn’t move too far from the origin, the restoring force will be approximately linear.
Let’s compare this to the linear approximation above in which we assumed a spherical charge distribution.
Note: We’re using the words “linear” and “nonlinear” a lot to get you used to what is actually nonlinear in nonlinear optics (the polarizibility!).
To summarize so far, we have calculated the electric field inside a Hydrogen atom using a linear approximation and the full quantum mechanical wavefunction:
E_approximate = 1/(4*pi*epsilon_0)*e*r/a0**3
fig,ax = plt.subplots()
ax.plot(r/1e-12,E_analytic, label = r'$E(r)$ actual', color='k') #plot Electric Field
ax.plot(r/1e-12,E_approximate, label = r'$E(r)$ approximate', color='r') #plot Electric Field
ax.vlines(a0/1e-12, 0, 1e12, colors='k', linestyle='dashed', linewidth=1, label=r'$a_0$ = {:.2f} pm'.format(a0/1e-12))
plt.xlabel('Distance [pm]')
plt.legend()
plt.xlim(0, 1.1*a0/1e-12)
plt.ylim(0, 0.5e12)
plt.xlabel('Distance [pm]')
plt.ylabel(r'$E(r)$ [V/m]')
plt.show()
That’s close but not perfect. If we do a power series expansion of the equation for the analytic field, we obtain
Looking at the linear term in the expansion, we find that the linear polarizibility of Hydrogen is
So our solid sphere model for the electric charge was off by a factor of \(4/3\) for the linear polarization, but still pretty close!
Let’s plot the correct value so that we have a satisfying pretty picture with linear and nonlinear curves:
E_taylor = 1/(3*pi*epsilon_0)*e*r/a0**3
fig,ax = plt.subplots()
ax.plot(r/1e-12,E_analytic, label = r'$E(r)$ actual', color='k') #plot Electric Field
ax.plot(r/1e-12,E_taylor, label = r'$E(r)$ power series', color='r') #plot Electric Field
ax.vlines(a0/1e-12, 0, 1e12, colors='k', linestyle='dashed', linewidth=1, label=r'$a_0$ = {:.2f} pm'.format(a0/1e-12))
plt.xlabel('Distance [pm]')
plt.legend()
plt.xlim(0, 1.1*a0/1e-12)
plt.ylim(0, 0.5e12)
plt.xlabel('Distance [pm]')
plt.ylabel(r'$E(r)$ [V/m]')
plt.show()
Ahh…that’s better. Note that the linear approximate devatiates significantly after about 10 pm from the orgin, where the electric field is approximately \(10^{11}\) V/m. This means we need a non-resonant external field of around \(10^{11}\) V/m to observe nonlinear effects! Fortunately, the oscillating electric fields in lasers can be tuned to be resonant with the oscillating motion of the proton-electron oscillating dipole moment, and hence not require such ridiculously large fields.
To further cement the analogy to mass-spring systems, let’s also plot the potential energy \(eV(r)\), just like it’s done in mass-spring systems. Recall that the potential \(V(r)\) is simply the line integral of the electric field from some point with zero potential energy (infinity) to the point in question.
r0 = 1e-6*a0 #set min radial coordinate
rf = 100*a0 #set max radial coordinate. Set really big so that potential is zero at this point.
N=50000 #number of points in r array
dr = (rf - r0)/(N) #dr step size
r = np.arange(r0,rf,dr) #radial coordinate array
def ode(V, r):
qenc_analytic = e*(-np.exp(-(2*r)/a0)*(a0**2 + 2*a0*r + 2*r**2)/a0**2)+e
A = 4*pi*r**2
E_analytic = qenc_analytic/epsilon_0/A
dV = E_analytic
return dV
V0 = 0
sol = integrate.odeint(ode, V0, np.flip(r)) #note we integrate from infinity to zero by flipping r
U_taylor = 0.25*1/(3*pi*epsilon_0)*e*r**2/a0**3 + np.min(sol)/2
fig = plt.figure()
plt.plot(np.flip(r)/1e-12, sol/2, color = 'k',label = "Actual")
plt.plot(-np.flip(r)/1e-12, sol/2, color = 'k')
plt.plot(r/1e-12,U_taylor, color = 'r', label = 'Linear Approximation')
plt.plot(-r/1e-12,U_taylor, color = 'r')
plt.vlines(a0/1e-12, -13.4, 0, colors='k', linestyle='dashed', linewidth=1, label=r'$a_0$ = {:.2f} pm'.format(a0/1e-12))
plt.grid()
plt.ylim(-13.5, 0)
plt.xlim(-2*a0/1e-12, 2*a0/1e-12)
plt.xlabel("$r$ [pm]")
plt.ylabel("$U(r)$ [eV]")
plt.legend()
plt.show()
The the potential energy for the linear approximation is a perfect parabola, just like in a mass-spring system. The potential energy for the actual system is approximately the same for small displacements, but much lower for larger displacements. Notice also that we recover the ground state energy of a Hydrogen atom of -13.6 eV.
1.7. Extra Credit Opportunity#
(20 pts) Can you write code that animates a “ball” in these potential wells, sloshing back an forth?
(40 pts, Much harder) For those who have taken quantum optics, can you solve the Schrodinger equation in the energy basis and plot the eigenfunctions for a linear and nonlinear harmonic oscillator with these potentials? Ignore radial coordinates for now and just assume a 1D potential. What is the spacing between the eigenvalues (energies) in each case?
1.8. Closing thought: what does polarizability have to do with optics?#
So we now know that in general the atomic polarization can be nonlinear. But what does inducing a dipole moment with an external electric field have to to with optics? A couple of things are worth pointing out:
First, for most of the nonlinear optics processes we will study, the external electric field is an optical field (often a laser). This optical field oscillates at optical frequencies (100’s of THz), and hence the relative position of the electron cloud and nucleus can oscillate at optical frequencies.
Second, dipoles produce electric fields themselves.
1.9. Bonus Material (Not Required)#
Remeber how we expanded the electric field in a power series? Let’s extend that a bit:
We might naturally wonder whether we could use this power series expansion can help us find the correct terms for 2nd, 3rd, and higher order nonlinearities. It turns out the answer is no. Why?
We have only expanded the field for positive values of \(r\), and so our Taylor expansion is not valid for negative values of \(r\). So all higher terms except the linear term are being “influenced” by these negative values. Say what? We can think about it this way: in writing the electric field as a function of \(r\), we have implicitly neglected the fact that the electric field always points radially outwards from the origin. So any line plot through the origin of the field projected onto one direction should be an odd function. The Taylor series expansion of that field should therefore only contain odd powers of \(r\), but our Taylor expansion clearly has even terms, so it can’t be used to find the nonlinear forces acting on an electron oscillating about the origin.
1.10. Mega-Bonus Opportunity#
Find an odd representation of the electric field of a Hydrogen atom that can be expanded in a Taylor series that accurately gives the nonliear polarization. I have not seen this in the literature, so it might be worth a peer-review publication if you can figure it out. Here is the beginning of an attempt I’ve made:
We know the positive values are correct, so we can numerically make an odd function that is correct simply by negating the field for negative values of \(r\). The strategy is then to find a power series expansion that converges within some radius and/or an analytic function for that field. Here’s a numerical approximation:
qenc_analytic = e*(-np.exp(-(2*r)/a0)*(a0**2 + 2*a0*r + 2*r**2)/a0**2)+e
A = 4*pi*r**2
E_analytic = qenc_analytic/epsilon_0/A
rodd = np.concatenate((-r[-1:0:-1], r))/a0 #Create an r-array that has negative and positive values
Eodd = np.concatenate((-E_analytic[-1:0:-1], E_analytic))/np.max(E_analytic) # Make E negative for negative r
#Plot the electric field
fig = plt.figure()
plt.plot(rodd/1e-12, Eodd)
plt.show()
# use SciPy's interpolate functions to fund an approximate Taylor expansion.
# Look for an expansion that has very small even terms
from scipy.interpolate import approximate_taylor_polynomial
from scipy.interpolate import interp1d
Ef = interp1d(rodd,Eodd) # Generate an interpolation function
fig = plt.figure()
plt.plot(rodd,Ef(rodd), color='k')
for degree in [1, 3, 5, 21]:
E_taylor = approximate_taylor_polynomial(Ef, 0, degree, 3, order=degree) #the scale parameter is very sensitive
E_taylor_trunc = np.poly1d(np.asarray(E_taylor)[-4::1])
plt.plot(rodd, E_taylor(rodd), linestyle='--', label=f"degree={degree}")
#plt.plot(rodd, E_taylor_trunc(rodd), linestyle='--',label=f"degree={degree} truncated")
print(np.polynomial.Polynomial(E_taylor_trunc))
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.0, shadow=True)
plt.tight_layout()
plt.axis([-5, 5, -1.1, 1.1])
plt.plot()
plt.show()
print(np.polynomial.Polynomial(E_taylor))
0.10218431402849848 - 3.907900746391579e-06·x¹
-0.03174901797013865 + 8.684223880561781e-07·x¹ + 0.3879228704925821·x² -
3.9079007460769075e-06·x³
-0.1551749040265325 + 2.3157930346684452e-06·x¹ + 0.8086852999295011·x² -
3.907900746650285e-06·x³
-3.993071446870283 + 3.4736894464082284e-05·x¹ + 2.709041470977424·x² -
3.907900746291206e-06·x³
3.04024358526746e-07 - 1.068377540892196e-10·x¹ -
1.651293752788933e-05·x² + 4.80770066653099e-09·x³ +
0.00039154695619807167·x⁴ - 9.248817263961262e-08·x⁵ -
0.00531950514117728·x⁶ + 9.930312114009993e-07·x⁷ +
0.04569349057083049·x⁸ - 6.5167691158593555e-06·x⁹ -
0.2584510748810722·x¹⁰ + 2.6910429478332887e-05·x¹¹ +
0.9732519948178003·x¹² - 6.937846480162269e-05·x¹³ -
2.4193834511552472·x¹⁴ + 0.00010704107698999987·x¹⁵ +
3.889342305399692·x¹⁶ - 9.031591908076342e-05·x¹⁷ -
3.993071446870283·x¹⁸ + 3.4736894464082284e-05·x¹⁹ +
2.709041470977424·x²⁰ - 3.907900746291206e-06·x²¹
#Look at just the first few terms.
E_taylor_trunc = np.polynomial.Polynomial(np.asarray(E_taylor)[-4::1])
print(np.polynomial.Polynomial(E_taylor_trunc))
-3.993071446870283 + 3.4736894464082284e-05·x¹ + 2.709041470977424·x² -
3.907900746291206e-06·x³