17. Combining Dispersion and Nonlinearity in Optical Pulse Propagation#

Name:

Total Points: –/100 pts

### UNCOMMENT AND RUN THIS CELL IF USING GOOGLE COLAB
# !pip install ipympl -q
# from google.colab import output
# output.enable_custom_widget_manager()

In the last two notebooks we have explored how pulses propagate in nonlinear media. We first examined how dispersion affect pulse propogation, and found that we can model dispersion by transforming to the frequency domain and allowing each frequency component to accummulate a different phase. We then examined how nonlinearity affects pulse propogation by studying how the intensity-dependent refractive index leads to self-phase modulation. In this notebook, we will combine these two effects and look at how dipsersion and nonlinearity act simultaneously.

In the last notebook, we derived the nonlinear wave equation that includes a time-varying amplitude:

\[ \begin{align} \frac{\partial \tilde{E}}{\partial z} = \frac{1}{2}j\beta_2\frac{\partial^2 \tilde{E}}{\partial T^2} -j\frac{\omega}{2c\epsilon_0n_0}\left(1-\frac{j}{\omega}\frac{\partial}{\partial T} \right)\tilde{P}^{NL} \end{align} \]

It is very often the case that the second term on the right (the time derivative of the polarization) is negligible, so we’ll first consider a scenario in which is the case. This simplifies the pulse propagation equation to

\[ \begin{align} \frac{\partial \tilde{E}}{\partial z} = \frac{1}{2}j\beta_2\frac{\partial^2 \tilde{E}}{\partial T^2} -j\frac{\omega}{2c\epsilon_0n_0}\tilde{P}^{NL} \end{align} \]

Let’s start by modelling what happens when the nonlinear polarization results from the intensity-dependent refractive index. In this case we can write

\[ \begin{align} \tilde{P}^{NL} = 3\epsilon_0\chi^{(3)}|E|^2 E \end{align} \]

so that we have

\[\begin{split} \begin{align} \frac{\partial \tilde{E}}{\partial z} &= \frac{1}{2}j\beta_2\frac{\partial^2 \tilde{E}}{\partial T^2} -j\frac{\omega}{2cn_0}3\chi^{(3)}|E|^2 E \\ \frac{\partial \tilde{E}}{\partial z} &= \frac{1}{2}j\beta_2\frac{\partial^2 \tilde{E}}{\partial T^2} -j\gamma |E|^2E \\ \end{align} \end{split}\]

where \(\gamma = 3 \omega \chi^{(3)}/ 2cn_0 \).

This equation is known as the nonlinear Schrodinger equation. The first term on the right is responsible for dispersion and the second term on the right is responsible for the nonlinearity. A very common (and simple)technique for solving this equation numerically is the Split-Step Fourier Method. We will not go through the details of the implementation, but the idea is that we will step forward a distance \(h\) in two steps. In the first step we’ll take care of dispersion by transforming to the frequency domain and letting each frequency component acquire its own phase shift. In the second step, will transform back to the time domain and take a nonlinear step. In practice, the frequency domain step is itself usually divided into two half-steps, one before and one after the time step. This reduces numerical errors.

Here is a simple implementation of the Split-Step Fourier method for self-phase modulation in a dispersive fiber. We’ll start with an initial pulse of the form

\[ \begin{align} E_0 = A_0e^{-(1-jC)\frac{t^2}{2\tau^2}} \end{align} \]

This is the same as the pulse we used in the previous two notebooks with two modifications:

We have included the possibility of a linear chirp in the initial pulse, parameterized by the constant \(C\).
We have divided out the carrier wave \(e^{j\omega_0 t}\).

We’ll also be a little more careful to make our time and frequency arrays have a length of \(2^N\) for more numerically efficient FFT’s.

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

P0 = .00064*10  #input power in Watts
A0 = np.sqrt(P0) #amplitude

alpha = 0 #Fiber loss value in dB/km
alph = alpha/4.343 # attenuation coefficient satisfying Pout = Pin*e^(-alpha*z)

gamma = 0.003  # effective nonlinearity in units of rad/(W*m).  Typical value for silica fiber.
#Note n_2 = gamma*Aeff*lambda/(2*pi), where Aeff is the effective core area of the fiber.
#Also note that the nonlinear phase shift is phi = gamma*P*z

tau = 50e-12  #initial pulse width
C = -2  #input chirp parameter 
b2 = -20e-27  #  2nd order dispersion in units of s^2/m
Ld = (tau**2)/np.absolute(b2)  # Dispersion length in meters (see equation 6.16 in New)

N = 10  # power of 2 of time array length
l = 2**N # number of elements in time array
t = np.linspace(-50*tau, 50*tau, l)  # time array
dt = t[1]-t[0]
E0 = A0*np.exp(-((1+1j*(-C))/2.0)*(t/tau)**2)  # initial optical pulse



#find FWHM length
fwhml = np.nonzero(np.absolute(E0) > np.absolute(np.max(np.real(E0))/2.0))  
fwhml = len(fwhml[0])


# frequency array    
domega = 2*pi/(dt*l)  # frequency discretization
omega = domega*np.arange(-l/2, l/2, 1) # frequency array
omega = np.fft.fftshift(omega) #shift frequency array to zero frequency
    
spectrum = np.fft.fft(E0) # fourier transform of time-domain pulse


fig, (ax1, ax2) = plt.subplots(2)
ax1.set_title('Initial Pulse')
ax1.plot(t, np.real(E0))
ax1.set_xlabel('Time [s]')
ax1.set_ylabel('Amplitude')

fig.subplots_adjust(hspace=.5)

ax2.plot(omega, np.abs(spectrum))
ax2.set_xlabel('Frequency [rad/s]')
ax2.set_ylabel('Amplitude')


plt.show()

#Split-step Fourier Method  (See Wikipedia Article on Split-Step Fourier for details)  
def split_step_fourier(spectrum, z, h):  
    for step in np.arange(h, z+1, h):
        spectrum = spectrum*np.exp(-alph*(h/2.0)+1j*b2/2.0*omega**2*(h/2.0)) # dispersive half step 
        f = np.fft.ifft(spectrum) # transform back to time domain
        f = f*np.exp(1j*gamma*np.abs(f)**2*h) # nonlinear propogation full step 
        spectrum = np.fft.fft(f) # transform back to frequency domain
        spectrum = spectrum*np.exp(-alph*(h/2.0)+1j*b2/2.0*omega**2*(h/2.0)) #dispersive half step

    f = np.fft.ifft(spectrum) #Fourier transform final pulse back to time domain
    return f, spectrum

h = Ld/50 # step size
spectrum = np.fft.fft(E0) # fourier transform of initial time-domain pulse

def update(z = 0):
    g, s = split_step_fourier(spectrum, z, h)
    line0.set_ydata(np.real(g))
    line1.set_ydata(np.abs(g))
    line2.set_ydata(np.abs(s))

fig, (ax1, ax2) = plt.subplots(2)
dummy_f = np.zeros_like(t)
line0, = ax1.plot(t, dummy_f, label = 'Time-Domain', color = 'red')
line1, = ax1.plot(t, dummy_f, label = 'Time-Domain', color = 'black')
line2, = ax2.plot(omega, dummy_f, label = 'Spectrum', color = 'blue')

ax1.set_xlabel('Time [s]')
ax1.set_ylabel('Amplitude')
ax1.set_ylim(-(1.5*A0), 1.5*A0)

fig.subplots_adjust(hspace=.5)

ax2.set_xlabel('Frequency [rad/s]')
ax2.set_ylabel('Amplitude')
ax2.set_ylim(-0.05, 1.5)


ipw.interact(update, z = (0, 10*Ld, 0.5*Ld))

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

17.1. Exercises#

(40 pts) Track the FWHM of the pulse in the time domain at every h-step. Hint: I’ve included above some code to calculate the length of the FWHM time array for the initial pulse.
(30 pts) Calculate the pulse broadening ratio at every step. This is the ratio of the FWHM to the original pulse. Question: why does the pulse get narrower at first when you start with a negative chirp?
(30 pts) Explore what happens when you increase the input pulse power. How does the pulse broadening and frequency content change? Why? At what power does the nonlinearity start to have as much effect as the dispersion?

Nonlinear Optics

Combining Dispersion and Nonlinearity in Optical Pulse Propagation

Contents

17. Combining Dispersion and Nonlinearity in Optical Pulse Propagation#

17.1. Exercises#