4. Coupled Wave Equations#

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()

We have seen in the previous notebooks how new frequency components can be generated by nonlinear optical processes. Those new optical fields can act as source fields for even further nonlinear interactions. In order to describe the coupled interactions, we need to solve the relevant wave equations. We’ll start with Maxwell’s equations. Here we go:

Maxwell’s equations with no free charges or curents can be written

\[\begin{split} \begin{align} \nabla \cdot \tilde{\mathbf{D}} &= \tilde{\rho} \\ \nabla \cdot \tilde{\mathbf{B}} &= 0 \\ \nabla \times \tilde{\mathbf{E}} &= -\frac{\partial \tilde{\mathbf{B}}}{\partial t} \\ \nabla \times \tilde{\mathbf{H}} &= \frac{\partial \tilde{\mathbf{D}}}{\partial t} \end{align} \end{split}\]

We assume that the material is nonmagnetic:

\[ \begin{align} \tilde{\mathbf{B}} = \mu_0 \tilde{\mathbf{H}} \end{align} \]

To allow for nonlinear interactions, we let \(\tilde{\mathbf{D}}\) and \(\tilde{\mathbf{E}}\) be related by

\[ \begin{align} \tilde{\mathbf{D}} = \epsilon_0 \tilde{\mathbf{E}} + \tilde{\mathbf{P}} \end{align} \]

where \(\tilde{\mathbf{P}}\) is the nonlinear polarization.

We now perform the usual trick of taking the curl of the curl-\(\tilde{\mathbf{E}}\) equation and substituting in the curl-\(\tilde{\mathbf{B}}\) equation:

\[\begin{split} \begin{align} \nabla \times \nabla \times \tilde{\mathbf{E}} &= -\frac{\partial}{\partial t} \nabla \times \tilde{\mathbf{B}} \\ & = -\mu_0\frac{\partial^2}{\partial t^2} \tilde{\mathbf{D}} \end{align} \end{split}\]

Substituting in the expression for \(\tilde{\mathbf{D}} \) we obtain

\[ \begin{align} \nabla \times \nabla \times \tilde{\mathbf{E}} + \epsilon_0\mu_0\frac{\partial^2}{\partial t^2} \tilde{\mathbf{E}} = - \mu_0\frac{\partial^2}{\partial t^2} \tilde{\mathbf{P}} \end{align} \]

This is the most general form of the wave equation, but it is somewhat difficult to work with because of the repated curl operations. We can use some vector calculus by remembering that \(\nabla \times \nabla \times \tilde{\mathbf{E}} = \nabla\left( \nabla \cdot \tilde{\mathbf{E}}\right) - \nabla^2 \tilde{\mathbf{E}} \) and noting that the \(\nabla\left( \nabla \cdot \tilde{\mathbf{E}}\right)\) term is generally very small. In these cases (the whole course!), we have the wave equation

\[ \begin{align} \nabla^2 \tilde{\mathbf{E}} - \epsilon_0\mu_0\frac{\partial^2}{\partial t^2} \tilde{\mathbf{E}} = \mu_0\frac{\partial^2}{\partial t^2} \tilde{\mathbf{P}} \end{align} \]

Let’s ignore for a minute the fact that \(\tilde{\mathbf{E}}\) and \(\tilde{\mathbf{P}}\) are vectorial fields, and write this wave equation down for the case of plane waves. We’ll also separate the polarization into its linear and nonlinear parts.

\[ \begin{align} \frac{\partial^2\tilde{E}}{\partial z^2} - \epsilon_0\mu_0\frac{\partial^2}{\partial t^2} \tilde{{E}} = \mu_0\frac{\partial^2}{\partial t^2} \left( \tilde{P}^{(1)}+ \tilde{P}^{(NL)} \right) \end{align} \]

We can now insert the electric field and polarization into these equations as we did in the last notebook. However, we add one thing: the propogating nature of the wave in the \(z\) direction. A single-frequency wave will then be

\[\begin{split} \begin{align} \tilde{E}(z,t) = E(z)e^{j(kz - \omega t)} + c.c. \\ \tilde{P}(z,t) = P(z)e^{j(kz - \omega t)} + c.c. \end{align} \end{split}\]

Note that \(E(z)\) is a function of \(z\), since the amplitude of the electric field will increase or decrease in space (this is what we’re solving fore!). But we have separated out the rapidly varying parts oscillating on time-scales of the optical wave.

Note that the representation of the polarization as a wave in space is different than the electric field, and can feel kind of funny. For any individual atom, the polarization oscillating in time makes sense. But there is not an atom at every point in space. So the polarization wave is a wave moving throught the atoms, which we assume are closely spaced enough that it’s fairly continuous. This is ok since the polarization is just a function of the electric field (which is continuos), so if we know what the field is, we also know what the polarization is.

Anyway, the complete collection of fields will can be represented as the sum of individual frequency components:

\[\begin{split} \begin{align} \tilde{E}(z,t) = \sum_n E_n(z)e^{j(k_nz - \omega_n t)} + c.c. \\ \tilde{P}(z,t) = \sum_n P_n(z)e^{j(k_nz - \omega_n t)} + c.c. \end{align} \end{split}\]

The wave equation requires that we take two derivatives in both time and space. Let’s look at the space derivative first:

\[ \begin{align} \frac{\partial^2\tilde{E}(z, t)}{\partial z^2} &= \sum_n -k_n^2 E_n(z)e^{j(k_nz - \omega_n t)} + 2jk_n\frac{\partial E_n(z)}{\partial z}e^{j(k_nz - \omega_n t)} + \frac{\partial^2 E_n(z)}{\partial z^2}e^{j(k_nz - \omega_n t)} \end{align} \]

It turns out the third term that involves a second spatial deriviative of \(E\) can be neglected when it is much smaller than the second term involving only a single derivative. This approximation is true whenever the field amplitude does not grow or shrink very much on the order of a wavelenght (distance of \(1/k\)).

The time derivatives have fewer complications, since the complete time dependence is contained in the exponential:

\[\begin{split} \begin{align} \frac{\partial^2\tilde{E}(z, t)}{\partial t^2} &= \sum_n -\omega_n^2 E_n(z)e^{j(k_nz - \omega_n t)} \\ \frac{\partial^2\tilde{P}(z, t)}{\partial t^2} &= \sum_n -\omega_n^2 P_n(z)e^{j(k_nz - \omega_n t)} \\ &= \sum_n -\omega_n^2 \left(P_n^{(1)}(z) + P_n^{(NL)}(z)\right)e^{j(k_nz - \omega_n t)} \\ &= \sum_n -\omega_n^2 \left(\epsilon_0\chi^{(1)}E_n(z) + P_n^{(NL)}(z)\right)e^{j(k_nz - \omega_n t)} \\ \end{align} \end{split}\]

Inserting these expressions back into the wave equation we obtain

\[ \begin{align} \sum_n -k_n^2 E_n(z)e^{j(k_nz - \omega_n t)} + 2jk_n\frac{\partial E_n(z)}{\partial z}e^{j(k_nz - \omega_n t)} + \epsilon_0\mu_0 \omega_n^2 E_n(z)e^{j(k_nz - \omega_n t)} = \mu_0\sum_n -\omega_n^2 \left(\epsilon_0\chi^{(1)}E_n(z) + P_n^{(NL)}(z)\right)e^{j(k_nz - \omega_n t)} \end{align} \]

We note that each frequency component must independently satisfy this wave equation, so we can drop the summations. Re-arranging a bit and dividing out all the exponentials we get

\[ \begin{align} 2jk_n\frac{\partial E_n(z)}{\partial z} + \left[-k_n^2 +\epsilon_0\mu_0 \omega_n^2 (\chi^{(1)}+1)\right]E_n(z) = - \omega_n^2\mu_0 P_n^{(NL)}(z) \end{align} \]

It turns out we have a choice for \(k_n\). We would like to choose it such that the field amplitude doesn’t change when there is no nonlinear polarization. In other words \(\frac{\partial E_n(z)}{\partial z} = 0\) when \(P_n^{(NL)}(z) = 0\). This requires that the second term on the left entirely vanish, or in other words that

\[\begin{split} \begin{align} k_n^2 &= \epsilon_0\mu_0 \omega_n^2 (\chi^{(1)}+1) \\ k_n &= \frac{\omega_n}{c}n_n \end{align} \end{split}\]

where we have used the relationships \(n_n = \sqrt{\chi^{(1)}+1}\) and \(\epsilon_0\mu_0 = 1/c^2\).

We finally can write the wave equation as

\[ \begin{align} \frac{\partial E_n(z)}{\partial z} &= j\frac{\omega_n}{2 \epsilon_0 c n_n} P_n^{(NL)}(z) \end{align} \]

This simple differential equation is the basis for nearly all coupled equations in nonlinear optics. Note, however that we have ignored any time dependence of the field amplitudes, and rather written them as \(E_n(z)\) and \(P_n(z)\). This means we can’t solve nonlinear optics problems with short pulses using this equation. We’ll tackle that later.

4.1. Second Harmonic Generation#

Let’s use this wave equation to look at second harmonic generation. Our first job is to find the nonlinear polarization terms \(P_n^{(NL)}(z)\) oscilatting at frequency \(\omega_n\). The electric field will have frequencies at \(\omega_1\) and \(\omega_2\):

\[ \begin{align} E(z,t) = E_1(z)e^{j(k_1z - \omega t)} + E_2(z)e^{j(k_2z - 2\omega t)} \end{align} \]

Likewise, nonlinear polarization will have fields oscillating at the same frequencies:

\[ \begin{align} P^{(NL)}(z,t) = P_1^{(NL)}(z)e^{j(k_1z - \omega t)} + P_2^{NL}(z)e^{j(k_2z - 2\omega t)} \end{align} \]

The second order nonlinear polarization is produced by \(\chi^{(2)}\) and has the form:

\[ \begin{align} P^{(NL)}(z, t) = \epsilon_0 \chi^{(2)}E^2(z,t). \end{align} \]

Let’s take the square of the electric field using SymPy:

E1, E2, omega, k1, k2, z, t = symbols('E1, E2, omega,k_1, k_2, z, t') #create symbols for the input electric fields.  We'll ignore epsilon_0 and chi for now.
Esq = (E1*exp(I*(k1*z - omega*t))+ functions.conjugate(E1)*exp(-I*(k1*z - omega*t)) + E2*exp(I*(k2*z - 2*omega*t))+ functions.conjugate(E2)*exp(-I*(k2*z - 2*omega*t)))**2
Esq.expand()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[2], line 1
----> 1 E1, E2, omega, k1, k2, z, t = symbols('E1, E2, omega,k_1, k_2, z, t') #create symbols for the input electric fields.  We'll ignore epsilon_0 and chi for now.
      2 Esq = (E1*exp(I*(k1*z - omega*t))+ functions.conjugate(E1)*exp(-I*(k1*z - omega*t)) + E2*exp(I*(k2*z - 2*omega*t))+ functions.conjugate(E2)*exp(-I*(k2*z - 2*omega*t)))**2
      3 Esq.expand()

NameError: name 'symbols' is not defined

We are looking for terms that oscillate at \(\omega\) and \(2\omega\). We see that

\[\begin{split} \begin{align} P_2^{(NL)}(z)e^{j(k_2z-2\omega t)} = \epsilon_0\chi^{(2)}E_1^2 e^{j(2k_1z - 2\omega t)} \\ P_1^{(NL)}(z)e^{j(k_1z-\omega t)} = 2\epsilon_0\chi^{(2)}E_2E_1^* e^{j\left[(k_2-k_1)z - \omega t\right]} \end{align} \end{split}\]

and we can solve for the nonlinear polarizations:

\[\begin{split} \begin{align} P_2^{(NL)}= \epsilon_0\chi^{(2)}E_1^2 e^{-j\Delta k} \\ P_1^{(NL)}(z) = 2\epsilon_0\chi^{(2)}E_2E_1^* e^{j\Delta k} \end{align} \end{split}\]

where \(\Delta k = k_2 - k_1\).

We can now plug these nonlinear polarizations into Equation 20 to obtain the coupled differential equations that describe second harmonic generation:

\[\begin{split} \begin{align} \frac{\partial E_1(z)}{\partial z} &= j\frac{\omega}{4 c n_1}\chi^{(2)}E_2E_1^* e^{j\Delta k z} \\ \frac{\partial E_2(z)}{\partial z} &= j\frac{2\omega}{2 c n_2}\chi^{(2)}E_1^2 e^{-j\Delta k z} \end{align} \end{split}\]

Let’s solve these equations. The first and simplest solutions is to make the assumption that the amount of second harmonic generation is so small that it steals hardly any energy from the fundamental field. In this case, we set \(E_1(z)\) to be a constant and simply integrate the equation for \(\frac{\partial E_2(z)}{\partial z} \).

#init_printing(use_unicode=False,wrap_line=False)
dk, omega, n2, c, chi2, z = symbols('\Delta_k, \omega, n_2, c, \chi^{(2)}, z')
shg = integrate(I*omega/c/n2*chi2*E1**2*exp(I*dk*z),(z, 0, z))
simplify(shg)

Let’s plot this result

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

def update(Lcoh = 10e-3*np.pi):
    dk = np.pi/Lcoh
    E2 = E1**2*chi2*omega*(np.exp(1j*dk*z)-1)/(dk*c*n2) #second harmonic electric field
    I1 = c*epsilon_0**n1*np.abs(E1)**2 #Intensity of input beam
    I2 = c*epsilon_0**n2*np.abs(E2)**2  #Intensity of second harmonic
    line0.set_ydata(I2/I1) #conversion fraction

z = np.linspace(0, 10e-2, 1000) # propagation distance (10 cm crystal)
E1 = 1e4
chi2 = 1e-10
omega = 500e12
n1 = 1
n2 = 1

fig, ax = plt.subplots()
dummy_I = np.zeros_like(z)
line0, = ax.plot(z, dummy_I, label = 'I2', color = 'blue')
ax.set_ylabel('SHG Fraction')
ax.set_xlabel('z')
ax.set_ylim(0, 2e-2)


ipw.interact(update, Lcoh = (1e-4, 20e-2, 5e-4))

plt.show()

4.2. Exercises#

  1. (80 pts) Use python to solve for the second harmonic generation intensity without assuming that \(E_1(z)\) is constant. Plot the fundamental and second harmonic intensities in the same plot, and animate the result with a slider bar for the coherence length. Use the same parameters as the example above.

  2. (20 pts) What is the physical meaning of the wave-vector mismatch \(\Delta k\) and the coherence length \(L_{coh}\)?