# standard python imports
import matplotlib.pyplot as plt
import numpy as np
import tidy3d as td
# tidy3D import
import tidy3d.web as web
from numpy import random
# define geometry
wg_width = 0.25
wg_length = 2.5
wg_spacing = 0.5
buffer = 1.0
# compute quantities based on geometry parameters
x_span = 2 * wg_spacing + 2 * wg_length + 2 * buffer
y_span = wg_width + 2 * buffer
wg_insert_x = wg_length + wg_spacing
TIDY3D
LEARNING CENTER
Example Library
Generation of Kerr sideband¶
Note: the cost of running the entire notebook is larger than 1 FlexCredits.
Nonlinear materials are of great interest in many industries due to their unique capability to exhibit nontrivial phenomena, such as wave mixing and frequency generation. An interesting phenomenon that occurs in media with third-order nonlinearities is the Kerr sidebands, a form of four-wave mixing. Due to phase-matching conditions, these sidebands appear at frequencies offset by integer multiples of the frequency difference between the pump and signal.
Kerr sidebands can significantly enhance sensor resolution through all-optical signal processing. By combining a tunable laser source and a pump laser with a Fiber Bragg Grating (FBG) in a nonlinear fiber, frequency sidebands are generated, with each sideband’s power dependent on input intensities. Filtering specific sidebands narrows the FBG-reflected signal, improving wavelength shift detection. This technique enhances FBG-based temperature sensors and can be applied to other optical systems for increased resolution.
In this notebook, we use Tidy3D to perform a simulation of the generation of Kerr sideband in a waveguide. This example is based on the paper from Ole Krarup, Chams Baker, Liang Chen, and Xiaoyi Bao " Nonlinear resolution enhancement of an FBG based temperature sensor using the Kerr effect." Optics Express Vol. 28, Issue 26, pp. 39181-39188 (2020) Doi: https://doi.org/10.1364/OE.411179
This example was kindly created by Dr. Chenchen Wang, Postdoctoral Researcher at the University of Wisconsin–Madison.
For more examples, please refer to our learning center, where you can find tutorials such as Bistability in photonic crystal microcavities.
FDTD simulations can diverge due to various reasons. If you run into any simulation divergence issues, please follow the steps outlined in our troubleshooting guide to resolve it.
Theoretical base
In this section, we will introduce the theoretical framework of generating sidebands in a Kerr medium, as developed in the reference paper.
To generate Kerr sidebands, we inject laser light with two distinct angular frequencies into a Kerr medium. The two frequencies are a signal frequency \(\omega_s\) and a pump frequency \(\omega_p\), where \(\omega_s < \omega_p\). The Kerr medium is a waveguide made of a \(\chi^{(3)}\) material. The total electric field amplitude at the input of the fiber is given by:
\[A_{\text{in}} = \left[\sqrt{P_s} \exp(-0.5i\omega_s t) + \sqrt{P_p} \exp(0.5i\omega_p t)\right]\]
where \(P_s\) and \(P_p\) are the powers of the signal and pump fields, respectively. The angular frequency difference between the signal and pump is defined as \(\omega_d = \omega_p - \omega_s\). The input field’s power is:
\[|A_{\text{in}}|^2 = P_s + P_p + 2 \sqrt{P_s P_p} \cos(\omega_d t)\]
Neglecting the effects of dispersion, loss, and polarization, the evolution of this field is governed by the Non-linear Schrödinger Equation (NLSE):
\[\frac{dA}{dz} = i\gamma |A|^2 A\]
where \(A\) is the complex amplitude of the electric field, \(z\) is the propagation distance, and \(\gamma\) is the nonlinear coefficient of the medium. The term \(i\gamma |A|^2\) represents the third order nonlinear interaction of the electric field with the medium, where the field strength \(|A|^2\) is proportional to the intensity.
Solving this differential equation, the field at the output of the waveguide is given by:
\[A_{\text{out}} = A_{\text{in}} \exp\left[i\gamma L (P_s + P_p)\right] \exp\left[i \gamma L \cdot 2 \sqrt{P_s P_p} \cos(\omega_d t)\right]\]
Here, \(L\) is the length of the waveguide. The exponential term \(\exp[i\gamma L (P_s + P_p)]\) can be neglected as it does not affect the overall output power.
The term involving \(\cos(\omega_d t)\) can be expanded using the Jacobi-Anger expansion:
\[\exp[i M \cos(\Omega t)] = \sum_{n=-\infty}^{\infty} i^n J_n(M) \exp(in\Omega t)\]
where \(J_n(M)\) is the Bessel function of the first kind of order \(n\). Applying this expansion to the second exponential term:
\[\exp\left[i \gamma L \cdot 2 \sqrt{P_s P_p} \cos(\omega_d t)\right]\]
The output field becomes a sum of frequency sidebands, spaced by \(\omega_d\):
\[A_{\text{out}} = \sum_{n=-\infty}^{\infty} i^n \exp\left(i \omega_d (n + \frac{1}{2}) t \right) \left[i J_{n+1}\left( 2 \gamma L \sqrt{P_s P_p} \right) \sqrt{P_s} + J_n \left( 2 \gamma L \sqrt{P_s P_p} \right) \sqrt{P_p}\right]\]
This result indicates that the output field is a superposition of several frequency components, with the sidebands separated by \(\omega_d\).
The power of the \(n\)-th sideband is given by:
\[|A_n|^2 = J_{|n+1|}^2 \left( 2 \gamma L \sqrt{P_s P_p} \right) P_s + J_{|n|}^2 \left( 2 \gamma L \sqrt{P_s P_p} \right) P_p\]
Introducing the normalization \(x = \gamma L P_p\), \(y = \gamma L P_s\), we have:
\[z_n = x J_{n+1}^2 \left( 2 \sqrt{x y} \right) + y J_n^2 \left( 2 \sqrt{x y} \right)\]
Here we agree that \(z_0\) refers to \(y\), \(z_{-1}\) refers to \(x\) and the Indexes grow from low frequency to high frequency.
Under the condition \(0 < M < \sqrt{1 + n}\), the power in the \(n\)-th order sideband can be approximated as:
\[z_n \approx x^{|n+1| + 1} y^{|n+1|} \frac{1}{\left[ (|n+1|)! \right]^2} + x^{|n|} y^{|n|+1}\frac{1}{\left[ (|n|)! \right]^2}\]
This equation shows that the normalized output power is proportional to the normalized input power, raised to an integer exponent.
As an example, filtering out the \(n = 1\) sideband yields:
\[z_{1} \approx x^3y^2\frac{1}{4} + xy^2\]
Similarly for \(n = -2\) sideband:
\[z_{-2} \approx y^3x^2\frac{1}{4} + yx^2\]
This implies an symmetry between \(z_{1}\),\(z_{-2}\) and \(x\),\(y\) which can be verified later in simulation.
Initial setup
First we start defining the parameters for the simulation:
Define frequency:
# wavelength range of interest
lambda_beg = 0.5
lambda_end = 0.6
# define pulse parameters
freq_beg = td.C_0 / lambda_end
freq_end = td.C_0 / lambda_beg
freq0 = (freq_beg + freq_end) / 2
fwidth = (freq_end - freq0) / 1.5
freqd = 1e13
freqp = freq0 + 0.5 * freqd
freqs = freq0 - 0.5 * freqd
# frequency for the first sideband
freq1 = freqs + (freqp - freqs)
min_steps_per_wvl = 30
run_time = 5e-12Define Materials:
To define the \(\chi^{(3)}\) material, we will create a NonlinearSpec object with a NonlinearSusceptibility model. Since the underlying mechanism for Kerr sidebands is four-wave mixing, NonlinearSusceptibility is a suitable choice as it utilizes real electric fields.
The num_iters parameter can be used if the convergence is poor, although it can’t prevent a simulation with high nonlinearities from diverging, as we will discuss below.
n_bg = 1.0
n_solid = 1.5
background = td.Medium(permittivity=n_bg**2)
solid = td.Medium(permittivity=n_solid**2)
# define the nonlinear parameters
n_kerr_2 = 2e-8
kerr_chi3 = 4 * (n_solid**2) * td.constants.EPSILON_0 * td.constants.C_0 * n_kerr_2 / 3
amp = 400
chi3_model = td.NonlinearSpec(models=[td.NonlinearSusceptibility(chi3=kerr_chi3)], num_iters=10)
kerr_solid = td.Medium(permittivity=n_solid**2, nonlinear_spec=chi3_model)Define structures:
waveguide = td.Structure(
geometry=td.Box(
center=[0, 0, 0],
size=[td.inf, wg_width, 0.22],
),
medium=kerr_solid,
name="waveguide",
)Compute and visualize the waveguide modes.
td.config.simulation.use_local_subpixel = True
mode_plane = td.Box(
center=[-wg_insert_x, 0, 0],
size=[0, 1, td.inf],
)
sim_modesolver = td.Simulation(
center=[0, 0, 0],
size=[x_span, y_span, 3],
grid_spec=td.GridSpec.auto(min_steps_per_wvl=min_steps_per_wvl, wavelength=td.C_0 / freq0),
structures=[waveguide],
run_time=1e-12,
boundary_spec=td.BoundarySpec.all_sides(boundary=td.Periodic()),
medium=background,
)
mode_spec = td.ModeSpec(num_modes=2)
mode_simulation = td.ModeSimulation.from_simulation(
simulation=sim_modesolver, plane=mode_plane, mode_spec=mode_spec, freqs=[freq0]
)
mode_sim_data = mode_simulation.run_local()
mode_data = mode_sim_data.modesf, ((ax1, ax2, ax3), (ax4, ax5, ax6)) = plt.subplots(2, 3, tight_layout=True, figsize=(10, 6))
mode_sim_data.plot_field("Ex", "abs", mode_index=0, ax=ax1)
mode_sim_data.plot_field("Ey", "abs", mode_index=0, ax=ax2)
mode_sim_data.plot_field("Ez", "abs", mode_index=0, ax=ax3)
mode_sim_data.plot_field("Ex", "abs", mode_index=1, ax=ax4)
mode_sim_data.plot_field("Ey", "abs", mode_index=1, ax=ax5)
mode_sim_data.plot_field("Ez", "abs", mode_index=1, ax=ax6)
ax1.set_title("|Ex|: mode_index=0")
ax2.set_title("|Ey|: mode_index=0")
ax3.set_title("|Ez|: mode_index=0")
ax4.set_title("|Ex|: mode_index=1")
ax5.set_title("|Ey|: mode_index=1")
ax6.set_title("|Ez|: mode_index=1")
mode_data.to_dataframe()
TERMS & CONDITIONS
How the process works:
Try this real URL to see an example: