15. Dispersion in Optical Pulses#

Name:

Total Points: –/100 pts

In this notebook, you will explore how to numerically propagate a pulse including dispersion. We’ll also be getting familiar with taking Fourier transforms in Python.

First, let’s write some code that gives us a Gaussian pulse in the time domain:

import numpy as np
import matplotlib.pyplot as plt
import ipywidgets as ipw
from scipy.constants import *
%matplotlib widget

## Generate Gaussian Pulse in the time domain

lam0 = 1.0 # wavelength in microns 
f0 = c*1e-9*1e6/lam0 # cycles per ns
omega0 = 2*pi*f0  #angular frequency in rad/ns
T0 = 1/f0 # period in ns
tau = 2*T0 #width of gaussian pulse in time (ns)

fs = 5*f0

t = np.arange(-100*T0, 100*T0, fs)
N=1000

t = np.linspace(-100*T0, 100*T0, N)  #time array

A = 1  #peak field amplitude
E_pulse = A*np.exp(-t**2/(tau**2)) 
Et = E_pulse*np.exp(1j*omega0*t)
It = np.abs(Et)**2

#z = np.linspace(0, 20*lam, 1000)

fig, ax = plt.subplots()
ax.plot(t, It)

plt.show()

Ok, that looks good. Now let’s take the Fourier transform of that pulse to see what frequency components it contains. This can be a little tricky the first time you do it, since numerical FFT’s have a bit of subtlety. We’ll compare our result to an analytic formula in New to make sure we’ve got it right.

## Take Fourier Transform to see pulse in Frequency Domain

from scipy import fftpack
Ef = fftpack.fft(Et)/fs
freqs = fftpack.fftfreq(len(Ef))*fs

# Analytic solution for comparison
Efa = np.sqrt(pi)*A*tau*np.exp(-(2*pi*(freqs - f0))**2*tau**2/4)


fig, ax = plt.subplots()
ax.stem(freqs, np.abs(Ef)**2, use_line_collection=True)
ax.plot(freqs, np.abs(Efa)**2, color = 'red')
plt.show()

Notice that that Scipy’s FFT algorithm normalizes the FFT such that the area under the transformed function is preserve. We’d like to normalize to the intensity (square of the field), so we work that out briefly:

FFT normalizes such that \begin{align} \int |E|df = 1 \end{align}

We want \begin{align} A^2\int |E|^2df = 1 \end{align} \begin{align} A^2 = \frac{1}{\int |E|^2df} \end{align}

## Normalization
df = freqs[1]-freqs[0]
Afnorm = sum(np.abs(Ef)**2)*df
Efnorm = Ef/np.sqrt(Afnorm)
print(sum(np.abs(Efnorm)**2)*df)

fig, ax = plt.subplots()
ax.stem(freqs, np.abs(Efnorm)**2)
ax.plot(freqs, np.abs(Efa/np.sqrt(Afnorm))**2, color = 'red')
plt.show()
1.0000000000000009

Ok, now for the interesting part. We will assign each frequency component it’s own k-vector, allow it to propagate, and the take the inverse FFT to see what the pulse looks like in the time domain:

# Calculate the frequency-dependent refractive index

# slice containing nonzero frequency components
index0 = np.argmax(Efnorm)
indexi = int(index0-N/10)
indexf = int(index0+N/10)
freqslice = freqs[indexi:indexf]
lamslice = c/freqslice*1e-3
Efnormslice = Efnorm[indexi:indexf]


# Sellmeir coefficients for sapphire extraordinary wave
B1 = 1.5039759
C1 = 5.48041129e-3
B2 = 0.55069141
C2 = 1.47994281e-2
B3 = 6.5927379
C3 = 402.89514
B = [B1, B2, B3]
C = [C1, C2, C3]


eps = (1+(lamslice)**2*B[0]/((lamslice)**2-C[0]) 
      +(lamslice)**2*B[1]/((lamslice)**2-C[1])
      +(lamslice)**2*B[2]/((lamslice)**2-C[2]))

n = np.sqrt(eps)
k = n**2*pi*freqslice/c


fig, ax = plt.subplots()
ax.plot(lamslice, n)
plt.show()
# Propagate pulse
z = 0
tslice = t[indexi:indexf]
V = Efnormslice*np.exp(1j*k*z)
Etz = fftpack.ifft(V)

def update(z = 0):
    V = Efnormslice*np.exp(1j*k*z)
    Etz = fftpack.ifft(V)
    line0.set_ydata(np.abs(Etz)) #conversion fraction

fig, ax = plt.subplots()
dummy_I = np.zeros(V.shape)
line0, = ax.plot(-tslice, dummy_I, label = 'I2', color = 'blue')
#ax.plot(N, p2)
ax.set_ylabel('Intensity')
ax.set_xlabel('time')
ax.set_ylim(-1e-4, 1.5e-3)


ipw.interact(update, z = (0, 20000000, 100))

#plt.plot(z,I2)
plt.show()

You can see that the pulse widens and distorts.

15.1. Exercise#

  1. (30 pts) Choose another material (not Sapphire) and calculate the frequency-dependent index of refraction for that material using Sellmeir coefficients.

  2. (60 pts) Simulate the propgation of a short pulse in this medium. Try a few different pulse widths to see how the pulse spreads. How many meters does it take for a 100 ps pulse to double in width?

  3. (10 pts) Using your simulation, calculate the group velocity of the pulse.