6. Focused Beams#

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

So far, we have been studying nonlinear optics processsing under the assumption that all electric fields are plane waves. In reality, laser beams are not plane waves and furthermore are usually focused to achieve maximum intensity. In this notebook, we’ll examine some of the consequences of nonlinear optics experiements using focused beams.

The first important property of a laser beam is that it is more intense in the center than on the edges. This means that nonlinear conversion efficiency will be higher in the center. As a result, converted beams are narrower than the input beam.

Let’s go back to the place in our derivation of the wave equation wave equation where we assumed plane waves. After starting with Maxwell’s equations and taking taking the curl of the curl-\(\tilde{\mathbf{E}}\) equation and substituting in the curl-\(\tilde{\mathbf{B}}\) equation, we arrived at

\[ \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} \]

and then we assumed plane waves. If we do not assume plane waves at this point, the spatial part of the field will have transverse variations in addition to a periodic \(z\) dependence.

We can represent each frequency component of the electric field and polarization as

\[\begin{split} \begin{align} \tilde{\mathbf{E}}_n(\mathbf{r},t) = \mathbf{E}_n(\mathbf{r})e^{j(k_nz - \omega_n t)} + c.c. \\ \tilde{\mathbf{P}}_n(\mathbf{r},t) = \mathbf{P}_n(\mathbf{r})e^{j(k'_nz - \omega_n t)} + c.c. \\ \end{align} \end{split}\]

Notice that this is the same definition as before except we have allowed \(\mathbf{E}_n\) and \(\mathbf{P}_n\) to be spatially varying in the transverse dimensions (\(x, y\)) in addition to \(z\). This means that the Laplace operator \(\nabla^2\) will be \(\partial/dx^2 + \partial/dy^2 + \partial/dz^2 = \nabla_T^2 + \partial/dz^2\)

With this small change, we proceed to derive the wave equation the same as before. After taking spatial and temporal derivatives and making the slowly varying amplitude approximation, we arrive at

\[ \begin{align} 2jk_n\frac{\partial \mathbf{E}_n}{\partial z} + \nabla_T^2 \mathbf{E}_n = - \frac{\omega_n^2}{\epsilon_0 c^2 }\mathbf{P}_n \end{align} \]

This is known as the paraxial wave equation. Compare this to the plane wave version:

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

Notice that we have the extra transverse Laplacian and as a result have to keep the factor of \(2jk_n\) with the the \({\partial \mathbf{E}_n}/{\partial z}\) term. In the case of free propagation (e.g. \(\mathbf{P}_n = 0\)), a solution to the paraxial wave equation is a beam with a Gaussian transverse profile:

\[ \begin{align} E(r,z) = E_0 \frac{w_0}{w(z)}e^{-r^2/w(z)^2}e^{jkr^2/2R(z)e^{j\Phi(z)}} \end{align} \]

where

\[ \begin{align} w(z) = w_0\left[ 1+ \left( \lambda z / \pi w_0^2 \right) ^2 \right]^{1/2} \end{align} \]

is the \(1/e\) radius of the beam,

\[ \begin{align} R(z) = z\left[1 + \left( \pi w_0^2 / \lambda z \right)^2 \right] \end{align} \]

is the radius of curvature of the wavefront, and

\[ \begin{align} \Phi(z) = -\arctan \left( \lambda z / \pi w_0^2 \right) \end{align} is the phase. \]

Let’s plot this to get a feel feel for the Gaussian beam

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

n = 1
lam = 1.55e-6
w0 = 1e-3

zi = -5
zf = 5
Nz = 100
z = np.linspace(zi,zf,Nz)

ri = -5*w0
rf = 5*w0
Nr = 100
r = np.linspace(ri, rf, Nr)

w = w0*(1 + (lam*z/(pi*w0**2))**2)**(1/2)
R = z*(1+ (pi*w0**2)**2/(lam*z))
Phi = -np.arctan(lam*z/(pi*w0**2)**2)
E = w0/w*np.exp(-r**2/w**2)

fig, (ax1, ax2, ax3) = plt.subplots(3, 1)
ax1.plot(r,E) # plot slice of E at origin
ax2.plot(z, w, 'b') # plot beam width as a function of z
ax2.plot(z, -w, 'b') # plot negative also to make it look like a beam
ax3.plot(z, Phi) #plot phase
plt.show()

def Efield(r, z):
    w = w0*(1 + (lam*z/(pi*w0**2))**2)**(1/2)
    R = z*(1+ (pi*w0**2)**2/(lam*z))
    Phi = -np.arctan(lam*z/(pi*w0**2)**2)
    E = w0/w*np.exp(-r**2/w**2)
    return w0/w*np.exp(-r**2/w**2)


R1, Z1 = np.meshgrid(r, z)
ZE = Efield(R1, Z1)



cmap1 = plt.cm.Reds
fig, ax = plt.subplots()
imE = ax.imshow(ZE, cmap=cmap1, origin='lower', extent=[1000*ri,1000*rf,zi,zf])
ax.set_xlabel(r'$r$ [mm]')
ax.set_ylabel(r'$z$ [m]')
ax.set_title(r'Gaussian Beam')


plt.show()

6.1. Exercises#

  1. (20 pts) Animate the above plot so that you can change the beam radius \(w_0\) and the wavelength \(\lambda\).

To solve the paraxial equation when the nonlinear polarization is not included, we must solve the partial differential equation that includes \(r\) and \(z\). This usually involves setting up a mesh grid and defining the partial differential equations at each point, and then iterating until a solution is found with a desired accuracy.

Let’s get a feel for how these solutions look by looking at harmonic generation of the \(qth\) harmonic. In this case, we have \(\omega_q = q\omega\). The polarization in this case is:

\[ \begin{align} {P}_n = \epsilon_0 \chi^{(q)}E_1^{q}e^{j\Delta k z} \end{align} \]

where \(\Delta k = qk_1 - k_q\). Inserting this into the wave equation, we have

\[ \begin{align} 2jkq \frac{\partial E_q}{\partial z} + \nabla^2_T E_q = \frac{\omega_q^2}{c^2}\chi^{(q)}E_1^{q}e^{j\Delta k z} \end{align} \]

where

\[ \begin{align} E_1(r,z) = E_1 \frac{w_0}{w(z)}e^{-r^2/w(z)^2}e^{jkr^2/2R(z)e^{j\Phi(z)}} \end{align} \]

At this point it is usually more convenient to switch to a more compact form by letting

\[\begin{split} \begin{align} \zeta &= 2z/b \\ b &= 2\pi w_0^2/\lambda = kw_0^2 \end{align} \end{split}\]

so that

\[ \begin{align} E_1(r,z) = \frac{E_1}{1+j\zeta}e^{-qr^2/w_0^2(1+j\zeta)} \end{align} \]

We can guess that the solution to this equation will look like a Gaussian at each frequency. Following Boyd and others, we use the trial solution

\[ \begin{align} E_q(r,z) = \frac{E_q(z)}{1+j\zeta}e^{-qr^2/w_0(1+j\zeta)} \end{align} \]

This obeys the paraxial wave equation so long as (1) We work in the constant pump regime and (2) \(E_q(z)\) obeys the ordinary differential equation

\[ \begin{align} \frac{d E_q}{dz} = \frac{jq\omega}{2 n q c} \chi^{(q)}E_1^q \frac{e^{j\Delta k z}}{(1+j\zeta)^{(q-1)}} \end{align} \]

6.2. Exercises#

Solve Eq 8 numerically using lithium niobate as the nonlinear medium. Explore a range of beam radii, field intensities, and phase matching values \(\Delta k\). Make plots similar to the 2D mesh plot above.

  1. (40 pts) Compare your solution to the analytic solution in Boyd pg 113-114. Keep working at it until they match.

  2. (20 pts) Verify that the beam decreases in radius for higher harmonics

  3. (20 pts) What do you notice when \(\Delta k\) = 0? Explain the Guoy phase shift and why no frequency conversion takes place?