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:

Define 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.

Define structures:

Compute and visualize the waveguide modes.

Here we choose TM mode, then we create two mode sources $P_s$,$P_p$.

To ensure spectral overlap between the signal and pump sources, we will define the time profile of the sources as ContinuousSource, which simulates a continuous wave (CW) source.

Note that a GaussianPulse time profile is preferred in most cases, as the fields decay at the end of the simulation, making normalization in frequency domain meaningful. In contrast, for a ContinuousSource, the fields do not decay, so frequency domain normalization is not meaningful. Nevertheless, the temporal mean power of the source is equal to the square of the amplitude parameter.

Also, note that the fwidth argument controls the ramping of the CW amplitude rather than the spectral width, since a CW source is nearly monochromatic.

Define monitors:

Define simulation.

Note that we will set the shutoff argument as 0, as we are using CW sources so the fields will not decay. For the same reason, we can also disregard the warning messages in the log.

Visualize structure and sources.

Run simulation¶

Next, we can plot the fields at the frequency of the first sideband. As there is no source injecting power in this frequency, the fields are generated due to the nonlinear process.

Now, we can visualize the transmittance spectrum, where the two sidebands are visible at frequencies equally spaced from the pump and signal pulse.

Amplitude analysis¶

We can now vary the amplitude parameter and observe the power at the n = 1 sideband (at 565 Thz). Since the amplitudes of the pump and signal are identical, we expect the power of this band to vary approximately with $\text{amplitude}^5$.

It can be observed that the power in the n = 1 band follows a power law until amplitude values around 500. Beyond this point, the results deviate from the exponential trend until the simulation diverges for amplitudes greater than 700.

This occurs because the simulations remain stable only for small nonlinearities. The nonlinear permittivity should be smaller than the linear value. Therefore:

$ \epsilon_0 \chi^{(3)} E^3 \ll \epsilon_0 \chi^{(1)} E$

$ \chi^{(3)} E^2 \ll n_0^2-1$

Using the Poynting theorem, and the relation $\chi^{(3)} = (4/3)n_0^2 \epsilon_0 c n_2$, we have:

$\frac{8}{3}\frac{n_0 n_2}{n_0^2 - 1} I \ll 1$

Since the Intensity ($I$) is proportional to $\text{amplitude}^2$, the simulation can quickly become unstable for high amplitude values.

Additionally, we should mention that the $\chi^{(3)}$ formalism is only perturbative and may not be accurate for extremely intense electric fields.

Further verification¶

In the previous derivation, we obtained the following conclusions about $z_{-2}$,$z_{1}$:

$$ z_{-2} \approx x^2 \left( y + \frac{y^3}{4} \right) $$

$$ z_{1} \approx y^2 \left( x + \frac{x^3}{4} \right) $$

To verify this symmetry, We set the amplitude of the two input signals to half of the original amplitude and observe the change of the sideband signal.

Visualize the output spectrum of two further simulations.

From the figure above, we can observe the symmetry for $z_{-2}$, $z_{1}$ and $x$, $y$, which is consistent with the theoretical derivation.