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:

# 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

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-12

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.

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.

from tidy3d.plugins.mode import ModeSolver
from tidy3d.plugins.mode.web import run as run_ms

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_solver = ModeSolver(
    simulation=sim_modesolver, plane=mode_plane, mode_spec=mode_spec, freqs=[freq0]
)
mode_data = run_ms(mode_solver)

05:03:54 UTC Mode solver created with                                           
             task_id='fdve-b82c2909-3821-4280-afd2-3512513cc3f0',               
             solver_id='mo-5ecccd28-3df1-4447-aa2d-0b8e11fb29da'.               

05:03:56 UTC Mode solver status: success                                        

f, ((ax1, ax2, ax3), (ax4, ax5, ax6)) = plt.subplots(2, 3, tight_layout=True, figsize=(10, 6))

mode_solver.plot_field("Ex", "abs", mode_index=0, ax=ax1)
mode_solver.plot_field("Ey", "abs", mode_index=0, ax=ax2)
mode_solver.plot_field("Ez", "abs", mode_index=0, ax=ax3)
mode_solver.plot_field("Ex", "abs", mode_index=1, ax=ax4)
mode_solver.plot_field("Ey", "abs", mode_index=1, ax=ax5)
mode_solver.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()

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.

mode_source_p = td.ModeSource(
    size=mode_plane.size,
    center=mode_plane.center,
    source_time=td.ContinuousWave(freq0=freqp, fwidth=freq0 / 10, amplitude=amp),
    mode_spec=td.ModeSpec(num_modes=2),
    mode_index=1,
    direction="+",
    num_freqs=1,
)

mode_source_s = td.ModeSource(
    size=mode_plane.size,
    center=mode_plane.center,
    source_time=td.ContinuousWave(freq0=freqs, fwidth=freq0 / 10, amplitude=amp),
    mode_spec=td.ModeSpec(num_modes=2),
    mode_index=1,
    direction="+",
    num_freqs=1,
)

Define monitors:

# field monitor for the first sideband
field_monitor = td.FieldMonitor(
    center=[0, 0, 0], size=[td.inf, td.inf, 0], freqs=[freq1], name="field"
)
# monitor the mode amps on the output waveguide
lambdas_measure = np.linspace(lambda_beg, lambda_end, 1001)
freqs_measure = td.C_0 / lambdas_measure[::-1]


mode_monitor = td.ModeMonitor(
    size=mode_plane.size,
    center=mode_plane.center,
    freqs=freqs_measure,
    mode_spec=td.ModeSpec(num_modes=2),
    name="mode",
)

mode_monitor = mode_monitor.copy(update=dict(center=[wg_insert_x, 0, 0]))

# flux monitor
flux_monitor = td.FluxMonitor(
    center=(3.5, 0, 0), size=mode_plane.size, name="fluxMon", freqs=freqs_measure
)

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.

# create simulation
sim = td.Simulation(
    normalize_index=None,
    center=[0, 0, 0],
    size=[x_span, y_span, 0],
    grid_spec=td.GridSpec.auto(min_steps_per_wvl=min_steps_per_wvl, wavelength=td.C_0 / freq0),
    structures=[waveguide],
    sources=[mode_source_p, mode_source_s],
    monitors=[field_monitor, mode_monitor, flux_monitor],
    run_time=run_time,
    boundary_spec=td.BoundarySpec(
        x=td.Boundary.pml(), y=td.Boundary.pml(), z=td.Boundary.periodic()
    ),
    medium=background,
    shutoff=0,
)

Visualize structure and sources.

# plot the two simulations
fig, ax = plt.subplots(1, 1, figsize=(6, 6))
sim.plot_eps(z=0.01, ax=ax)
plt.show()

05:03:59 UTC WARNING: An appropriate frequency could not be determined when     
             plotting the permittivity. The permittivity will be evaluated at   
             infinite frequency. Please supply a value for `freq` to plot at a  
             finite frequency.                                                  

Run simulation

sim_data = web.run(sim, task_name="kerr_sideband", path="data/simulation_data.hdf5", verbose=True)

05:04:00 UTC Created task 'kerr_sideband' with resource_id                      
             'fdve-f312f378-63bd-499e-ad8d-b7017548eca3' and task_type 'FDTD'.  

             View task using web UI at                                          
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-f312f378-63b
             d-499e-ad8d-b7017548eca3'.                                         

             Task folder: 'default'.                                            

05:04:01 UTC Estimated FlexCredit cost: 0.101. This assumes the FDTD solver runs
             for the full simulation time; if early shutoff is reached, the     
             billed cost can be lower. Use 'web.real_cost(task_id)' to get the  
             billed FlexCredit cost after a simulation run.                     

05:04:02 UTC status = success                                                   

05:04:04 UTC Loading results from data/simulation_data.hdf5                     

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.

# visualize the fields
fig, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True, figsize=(15, 6))

ax1 = sim_data.plot_field("field", "Hz", val="real", z=0, ax=ax1)
ax2 = sim_data.plot_field("field", "S", val="abs", z=0, ax=ax2)

plt.show()

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

transmission_amps = sim_data["mode"].amps.sel(mode_index=1, direction="+") ** 2
f, ax = plt.subplots(figsize=(10, 5))
transmission_amps.abs.plot.line(x="f", ax=ax, label="absolute value")
ax.legend()
ax.set_title("Flux of the fundamental TM mode (forward)")
ax.set_ylim(0, None)
ax.set_xlim(freqs_measure[0], freqs_measure[-1])
ax.set_ylabel("Power (a.u.)")
plt.show()

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\).

Simulations = {}
Amplitudes = [50, 100, 150, 200, 250, 300, 350, 400, 450, 500, 550, 600, 700]
for amplitude in Amplitudes:
    # creating the sources with the new amplitude value
    s1, s2 = sim.sources
    st1 = s1.source_time.updated_copy(amplitude=amplitude)
    st2 = s2.source_time.updated_copy(amplitude=amplitude)

    sim2 = sim.updated_copy(
        sources=[s1.updated_copy(source_time=st1), s2.updated_copy(source_time=st2)]
    )

    Simulations[str(amplitude)] = sim2

batch = web.Batch(simulations=Simulations, verbose=True)
results = batch.run(path_dir="batch_simulations")

05:04:06 UTC Got 1 simulation from cache.                                       

05:04:14 UTC Started working on Batch containing 13 tasks.                      

05:04:26 UTC Maximum FlexCredit cost: 1.213 for the whole batch.                

             Use 'Batch.real_cost()' to get the billed FlexCredit cost after    
             completion.                                                        

05:04:32 UTC Batch complete.                                                    

Amps = []
Power1 = []

for amplitude in Amplitudes:
    sim_data = results[f"{amplitude}"]

    # recording the power for the band n = 1
    band1 = sim.sources[0].source_time.freq0 + (
        sim.sources[0].source_time.freq0 - sim.sources[1].source_time.freq0
    )
    bm = (sim_data["fluxMon"].flux.f > band1 * 0.99) & (sim_data["fluxMon"].flux.f < band1 * 1.01)
    max = sim_data["fluxMon"].flux[bm].max()

    Amps.append(amplitude)
    Power1.append(max)

# fitting the data with a polynomial function
from scipy.optimize import curve_fit

func = lambda X, a: a * X**5
Y = np.array(Power1) * 10**19
res, err = curve_fit(func, Amps[:9], Y[:9], p0=[1e-9], bounds=([0], [1]))

fig, ax = plt.subplots()
ax.plot(Amps, Y, "o")
ax.plot(Amps, func(np.array(Amps), *(res)))

ax.set_xlabel("Amplitude value")
ax.set_ylabel("First sideband power (a.u.)")

ax.set_ylim(-0.1, 1.1 * Y.max())

plt.show()

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:

$ _0 ^{(3)} E^3 _0 ^{(1)} E$

$ ^{(3)} E^2 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.

mode_source_p = mode_source_p.copy(
    update=dict(source_time=td.ContinuousWave(freq0=freqp, fwidth=freq0 / 10, amplitude=amp))
)
mode_source_s = mode_source_s.copy(
    update=dict(source_time=td.ContinuousWave(freq0=freqs, fwidth=freq0 / 10, amplitude=amp / 2))
)

sim_1 = td.Simulation(
    normalize_index=None,
    center=[0, 0, 0],
    size=[x_span, y_span, 0],
    grid_spec=td.GridSpec.auto(min_steps_per_wvl=min_steps_per_wvl, wavelength=td.C_0 / freq0),
    structures=[waveguide],
    sources=[mode_source_p, mode_source_s],
    monitors=[field_monitor, mode_monitor],
    run_time=run_time,
    boundary_spec=td.BoundarySpec(
        x=td.Boundary.pml(), y=td.Boundary.pml(), z=td.Boundary.periodic()
    ),
    medium=background,
    shutoff=0,
)

mode_source_p = mode_source_p.copy(
    update=dict(source_time=td.ContinuousWave(freq0=freqp, fwidth=freq0 / 10, amplitude=amp / 2))
)
mode_source_s = mode_source_s.copy(
    update=dict(source_time=td.ContinuousWave(freq0=freqs, fwidth=freq0 / 10, amplitude=amp))
)
sim_2 = td.Simulation(
    normalize_index=None,
    center=[0, 0, 0],
    size=[x_span, y_span, 0],
    grid_spec=td.GridSpec.auto(min_steps_per_wvl=min_steps_per_wvl, wavelength=td.C_0 / freq0),
    structures=[waveguide],
    sources=[mode_source_p, mode_source_s],
    monitors=[field_monitor, mode_monitor],
    run_time=run_time,
    boundary_spec=td.BoundarySpec(
        x=td.Boundary.pml(), y=td.Boundary.pml(), z=td.Boundary.periodic()
    ),
    medium=background,
    shutoff=0,
)

Simulations2 = {"1": sim_1, "2": sim_2}

batch2 = web.Batch(simulations=Simulations2, verbose=True)
results2 = batch2.run(path_dir="batch_simulations")

sim_data_1 = results2["1"]
sim_data_2 = results2["2"]

05:05:50 UTC Started working on Batch containing 2 tasks.                       

05:05:52 UTC Maximum FlexCredit cost: 0.202 for the whole batch.                

             Use 'Batch.real_cost()' to get the billed FlexCredit cost after    
             completion.                                                        

05:05:53 UTC Batch complete.                                                    

transmission_amps_1 = sim_data_1["mode"].amps.sel(mode_index=1, direction="+") ** 2
transmission_amps_2 = sim_data_2["mode"].amps.sel(mode_index=1, direction="+") ** 2
fig, (ax1, ax2) = plt.subplots(nrows=2, figsize=(10, 10))


transmission_amps_1.abs.plot.line(x="f", ax=ax1, label="absolute value (Top)")
ax1.legend()
ax1.set_title("Fundamental TM mode flux with half signal")
ax1.set_ylim(0, None)
ax1.set_xlim(freqs_measure[0], freqs_measure[-1])
ax1.set_ylabel("Power (a.u.)")

# Bottom plot1
transmission_amps_2.abs.plot.line(x="f", ax=ax2, label="absolute value (Bottom)")
ax2.legend()
ax2.set_title("Fundamental TM mode flux with half pump")
ax2.set_ylim(0, None)
ax2.set_xlim(freqs_measure[0], freqs_measure[-1])
ax2.set_ylabel("Power (a.u.)")

plt.tight_layout()
plt.show()