Bistability in photonic crystal microcavities

Introduction

Here we take advantage of Tidy3D’s nonlinear capability to demonstrate optical bistability. We base this example on the paper from Yanik et al. titled High-contrast all-optical bistable switching in photonic crystal microcavities. In this paper, a photonic crystal waveguide coupled with a point-defect cavity with Kerr nonlinearity achieves high-contrast transmission between two stable states.

For more simulation examples, please visit our examples page. If you are new to the finite-difference time-domain (FDTD) method, we highly recommend going through our FDTD101 tutorials. 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.

Setup

We first import the packages we’ll need.

# standard python imports
import gdstk
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import scipy

# tidy3D import
import tidy3d as td
import tidy3d.web as web

The waveguide-cavity system we examine is given by a square lattice (with lattice constant a) of high dielectric rods (\(n\) = 3.5) of radius \(.2a\). The waveguide is given by removing a row of these rods, and the cavity is given by a point defect in the crystal - instead of a rod, the defect is an ellipse with major and minor axes of \(a\) and \(2a\), respectively. Here we will set \(a\) = 1.

To build this in Tidy3D, we take advantage of the “square_cylinder_array” method in the common photonic crystal structures page in the Tidy3D learning center.

Since we are dealing with a photonic crystal with discrete translational symmetry, we can get spurious reflections if we simply use a PML boundary condition. Thus, instead we use the custom absorbing condition given in Mekis et al.’s Absorbing boundary conditions for FDTD simulations of photonic crystal waveguides. This is a \(k\)-matched distributed Bragg reflector waveguide consisting of a periodic array of alternating periodic slabs, with the waveguide given by a single line defect of a single slab with larger thickness. Since this waveguide has continuous translational symmetry, we can now feed it into a standard PML boundary. For simplicity, in defining the parameters for this waveguide, we use the parameters given in Mekis et al.

def square_cylinder_array(
    x0,
    y0,
    z0,
    R,
    hole_spacing_x,
    hole_spacing_y,
    n_x,
    n_y,
    height,
    medium,
    reference_plane="bottom",
    sidewall_angle=0,
    axis=2,
):
    # parameters
    # ------------------------------------------------------------
    # x0: x coordinate of center of the array (um)
    # y0: y coordinate of center of the array (um)
    # z0: z coordinate of center of the array (um)
    # R: radius of the circular holes (um)
    # hole_spacing_x: distance between centers of holes in x direction (um)
    # hole_spacing_y: distance between centers of holes in y direction (um)
    # n_x: number of cylinders in x direction
    # n_y: number of cylinders in y direction
    # height: height of array
    # medium: medium of the cylinders
    # reference_plane
    # sidewall_angle: angle slant of cylinders
    # axis

    cylinder_group = []

    start_x, start_y = x0 + hole_spacing_x * (1 - n_x) / 2, y0 + hole_spacing_y * (1 - n_y) / 2
    for i in range(0, n_x):
        for j in range(0, n_y):
            c = td.Cylinder(
                axis=axis,
                sidewall_angle=sidewall_angle,
                reference_plane=reference_plane,
                radius=R,
                center=(start_x + i * hole_spacing_x, start_y + j * hole_spacing_y, z0),
                length=height,
            )
            cylinder_group.append(c)
    structure = td.Structure(geometry=td.GeometryGroup(geometries=cylinder_group), medium=medium)
    return structure


def DBR(
    x0,
    y0,
    z0,
    length,
    thickness,
    spacing,
    num_layers,
    medium,
    direction="+",
):
    # parameters
    # ------------------------------------------------------------
    # x0: x coordinate of center of the array (um)
    # y0: y coordinate of center of the array (um)
    # z0: z coordinate of center of the array (um)
    # length: Length of slabs
    # thickness: thickness of slabs (um)
    # spacing: space between each slab (um)
    # num_layers: number of slabs
    # medium: medium of the slabs
    # direction: does it extend right or left from the center?

    slab_group = []
    orientation = 1
    if direction != "+":
        orientation = -1
    start_x, start_y = x0 + orientation * length / 2, y0 - num_layers / 2 * a
    for i in range(0, num_layers + 1):
        if i != num_layers // 2:
            s = td.Box(center=(start_x, start_y + i * a, 0), size=(length, thickness, td.inf))
            slab_group.append(s)
    structure = td.Structure(geometry=td.GeometryGroup(geometries=slab_group), medium=medium)
    return structure

Define Structure

We now carve out our waveguide and introduce or ellipse point defect. While the end goal of this notebook is to take advantage of the ellipse’s nonlinearity, we must first find the resonant frequency of the cavity when it behaves linearly - thus we set the material as a linear material for now.

We also want to compute \(\kappa\), a nonlinear feedback parameter crucial to our simulations - more on this later.

n_rods = 3.5
rod = td.Medium(permittivity=n_rods**2)
n_air = 1
air = td.Medium(permittivity=n_air)

# Kerr coefficient
n2 = 1.5e-17  # m^2/W, from paper
n2 *= 1e12  # convert to um^2/W

a = 1
radius = 0.2 * a
block_rows = 21
block_cols = 17

# ######################## SIMPLE CRYSTAL ##########################
block = square_cylinder_array(
    0, 0, -a, radius, a, a, block_cols, block_rows, td.inf, rod, reference_plane="middle"
)
# ######################## WAVEGUIDE WITHOUT CAVITY ##########################
waveguide = td.Structure(
    geometry=td.Box(
        center=[0, 0, 0],
        size=[td.inf, a, td.inf],
    ),
    medium=air,
    name="waveguide",
)
# ########################## WAVEGUIDE WITH CAVITY ###########################
cavity_air = td.Structure(
    geometry=td.Box(
        center=[0, -3 * a, 0],
        size=[a, a, td.inf],
    ),
    medium=air,
    name="cavity air",
)
air_block = td.Structure(
    geometry=td.Box(
        center=[-22 * a, 0, 0],
        size=[30 * a - a / 2, td.inf, td.inf],
    ),
    medium=air,
    name="air block",
)

# ################## DISTRIBUTED BRAGG REFLECTOR BOUNDARY ###################
thickness = 0.25 * a
n_DBR = 10.2
dbr_med = td.Medium(permittivity=n_DBR**2)
dbr_boundary = DBR(block_cols * a / 2, 0, 0, 20, thickness, a, block_rows - 1, dbr_med)

# ################################# CAVITY ##################################
major = a
minor = 0.2 * a
cavity_gdstk = gdstk.ellipse((0, -3 * a), (minor / 2, major / 2))
cavity_cell = gdstk.Cell("Ellipse")
cavity_cell.add(cavity_gdstk)
cavity = td.PolySlab.from_gds(
    cavity_cell,
    gds_layer=0,
    gds_dtype=0,
    axis=2,
    slab_bounds=(-td.inf, td.inf),
    reference_plane="bottom",
)[0]

cavity_structure = td.Structure(
    geometry=cavity,
    medium=rod,
)

freqs = np.linspace(1.08e14, 1.1e14, 1000)
freq0 = freqs[len(freqs) // 2]
fwidth = freqs[-1] - freqs[0]

run_time = 10000 / freq0

# simulation parameters
x_span = 27 * a
y_span = 17 * a

min_steps_per_wvl = 30
gridSpecUniform = td.GridSpec.uniform(dl=a / 24)

pulse = td.GaussianPulse(freq0=freq0, fwidth=fwidth)
source_distance = 9
gaussianBeam = td.GaussianBeam(
    center=[-source_distance * a, 0, 0],
    size=[0, td.inf, td.inf],
    source_time=pulse,
    direction="+",
    pol_angle=np.pi / 2,
    waist_radius=2.0 * a,
)

Things to Compute in the Linear Case

We wish to compute the resonant frequency \(\omega_{res}\) and the characteristic power of the cavity, given by

\(P_0=1/[\kappa Q^2 n_2 /c]\)

where \(Q\) is the cavity quality factor, \(n_2\) is the Kerr coefficient we’ll use, and \(\kappa\) is a nonlinear feedback parameter, a measure of the efficiency of the feedback system. We then use this to compute the cavity decay rate \(\gamma=\omega_{res}/(2Q)\).

To compute \(\omega_{res}\) and \(Q\), we use Tidy3D’s resonance finder.

In computing \(\kappa\), we use the definition given by Soljacic et al. in Optimal bistable switching in nonlinear photonic crystals:

\(\kappa=(\frac{c}{\omega_{res}})^2\frac{\int d^2r[|E(r)\cdot E(r)|^2 + 2|E(r)\cdot E^*(r)|^2]n^2(r)n_2(r)}{[\int d^2r|E(r)|^2n^2(r)]^2n_2}\)

Define Monitors and Simulation

We now input our monitors - a field monitor and permittivity monitor over the extent of the cavity’s mode to compute \(\kappa\) (field_monitor_cavity and permittivity_monitor_cavity), a field time monitor to compute our resonant frequency and cavity quality factor (field_time_monitor_wg), and a flux monitor at the end of the waveguide to find transmission vs frequency (flux_monitor).

cavity_monitor_size = 7
field_monitor_cavity = td.FieldTimeMonitor(
    fields=["Ez"],
    center=[0, -3 * a, 0],
    size=[cavity_monitor_size * a, cavity_monitor_size * a, 0],
    start=run_time * 8 / 10,
    stop=run_time * 8 / 10,
    name="field cavity",
)
permittivity_monitor_cavity = td.PermittivityMonitor(
    center=[0, -3 * a, 0],
    size=[cavity_monitor_size * a, cavity_monitor_size * a, 0],
    freqs=freqs[0],
    name="permittivity cavity",
)

flux_distance = 4.5

field_time_monitor_wg = td.FieldTimeMonitor(
    fields=["Ez"],
    center=[flux_distance * a, 0, 0],
    size=(0, 0, 0),
    start=run_time
    * 8
    / 10,  # time to start monitoring after source has decayed, units of 1/frequency bandwidth
    name="field time wg",
)
flux_monitor = td.FluxMonitor(
    center=[flux_distance * a, 0, 0],
    size=[0, td.inf, td.inf],
    freqs=freqs,
    normal_dir="+",
    name="flux",
)

x_boundary = td.Boundary(minus=td.PML(num_layers=80), plus=td.PML(num_layers=80))

sim_linear = td.Simulation(
    center=(0, -a, 0),
    size=[x_span, y_span, 0],
    medium=air,
    grid_spec=gridSpecUniform,
    structures=[block, waveguide, cavity_air, cavity_structure, dbr_boundary, air_block],
    sources=[gaussianBeam],
    monitors=[
        field_monitor_cavity,
        permittivity_monitor_cavity,
        field_time_monitor_wg,
        flux_monitor,
    ],
    run_time=run_time,
    boundary_spec=td.BoundarySpec(x=x_boundary, y=td.Boundary.periodic(), z=td.Boundary.periodic()),
    shutoff=False,
)

Visualize Simulation

sim_linear.plot(z=0)
plt.show()

Run Simulation

sim_linear_data = web.run(sim_linear, task_name="bistable_pc_cavity")

03:10:26 UTC Created task 'bistable_pc_cavity' with resource_id
             'fdve-d808890e-ad62-485f-b6ad-bf229be9543e' and task_type 'FDTD'.

             View task using web UI at
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-d808890e-ad6
             2-485f-b6ad-bf229be9543e'.

             Task folder: 'default'.

03:10:28 UTC Estimated FlexCredit cost: 0.097. 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.

03:10:29 UTC status = success

03:10:31 UTC Loading results from simulation_data.hdf5

Calculate Linear Resonance

We plot the flux through the waveguide vs frequency to verify that the the transmission vanishes on cavity resonance.

flx_cavity = sim_linear_data["flux"].flux
plt.plot(freqs, flx_cavity, label="Output Flux")
plt.xlabel("Frequency")
plt.ylabel("Transmission Flux")
plt.legend()
plt.show()

Compute Resonance Frequency and Q

We import Tidy3D’s resonance finder

from tidy3d.plugins.resonance import ResonanceFinder

resfinder = ResonanceFinder(freq_window=(freqs[0], freqs[-1]))

res_data = resfinder.run_scalar_field_time(signal=sim_linear_data.monitor_data["field time wg"].Ez)
res_data = res_data.sortby("error")
res_data.to_dataframe()

As you can see, the resonant frequency here agrees with the transmission spectrum shown above.

Q = res_data.Q[0].item()
fres = res_data.freq[0].item()
w_res = fres * 2 * np.pi
gamma = w_res / 2 / Q

Computing \(\kappa\)

We take advantage of our permittivity monitor data to use our \(n(r),n_2(r)\) functions needed in the integrals of the above definition of \(\kappa\). Note that, due to symmetry, the only nonzero field component is \(E_z\). To get \(n_2(r)\), we mask the permittivity monitor outside of the ellipse, and matshow to verify that this window only contains the nonlinear ellipse, and does not cut it off. We then compute the integrals and solve for \(\kappa\):

n_r = np.where(
    sim_linear_data.monitor_data["permittivity cavity"].eps_zz[1:, 1:, 0, 0] != 1.0 + 0.0j,
    n_rods,
    1,
)
n2_r = np.zeros(n_r.shape)
x, y = n_r.shape

# compute length between entries in monitor data
entry_length_x, entry_length_y = cavity_monitor_size * a / x, cavity_monitor_size * a / y
# compute number of entries that span a/2
wx, wy = int(a / entry_length_x / 2), int(3 * a / entry_length_y / 4)
# coordinates in entries of the monitor
rod_x, rod_y = x // 2, y // 2  # + wy//2

# isolate rod within window for n_2(r)
n2_r[rod_x - wx : rod_x + wx, rod_y - wy : rod_y + wy] = n_r[
    rod_x - wx : rod_x + wx, rod_y - wy : rod_y + wy
]
n2_r = np.where(n2_r > 1, n2, 0)

# get field components at resonant frequency
Ez = sim_linear_data.monitor_data["field cavity"].Ez[:, :, 0, 0]

integrand_num = (np.abs(Ez * Ez) ** 2 + 2 * np.abs(Ez * np.conjugate(Ez)) ** 2) * n_r**2 * n2_r
integrand_den = np.abs(Ez) ** 2 * n_r**2

int_num = integrand_num.integrate(coord=("x", "y")).item()
int_den = integrand_den.integrate(coord=("x", "y")).item()

kappa = (td.C_0 / w_res) ** 2 * int_num / ((int_den) ** 2 * n2)
print("kappa = ", kappa)

fig, ax = plt.subplots(1, 2)
ax[0].imshow(n2_r.T)
ax[0].set_title("$n_2(r)$")
ax[1].imshow(np.abs(Ez).T)
ax[1].set_title("$|E_Z(r)|$")
plt.show()

kappa =  0.16265699460554087

We can now define the characteristic power of the cavity:

P_0 = td.C_0 / (kappa * Q**2 * w_res * n2)

Waveguide normalization

We need to normalize the input power of our nonlinear simulation by the transmission intrinsic to the waveguide (without the cavity). The power frequency of our nonlinear source must be detuned by at least \(\delta=\sqrt{3}\), where

\(\delta=(\omega_{res}-\omega_{input})/\gamma\),

in order to be considered bistable. In Yanik et al., the frequency is detuned by \(2\sqrt{3}\).

detuning = 2 * np.sqrt(3)
detuned_freq = (w_res - detuning * gamma) / (2 * np.pi)
print(detuned_freq)

109024213322480.95

Define Waveguide-Only simulation for Normalization

We now update the previous simulation to only contain the waveguide in order to properly gauge the waveguide power input:

We take advantage of Tidy3D’s CustomSourceTime method to use this power function as a source with our desired detuned frequency:

sim_waveguide = sim_linear.updated_copy(
    monitors=[flux_monitor], structures=[block, waveguide, dbr_boundary, air_block]
)

sim_waveguide.plot(z=0)
plt.show()

sim_waveguide_data = web.run(sim_waveguide, task_name="bistable_pc")

03:10:35 UTC Created task 'bistable_pc' with resource_id
             'fdve-68e7d4c2-c417-4544-b838-9ee96aaacdfe' and task_type 'FDTD'.

             View task using web UI at
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-68e7d4c2-c41
             7-4544-b838-9ee96aaacdfe'.

             Task folder: 'default'.

03:10:36 UTC Estimated FlexCredit cost: 0.097. 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.

03:10:37 UTC status = success

03:10:38 UTC Loading results from simulation_data.hdf5

flx_waveguide = sim_waveguide_data["flux"].flux
detuned_transmission = sim_waveguide_data["flux"].flux.sel(f=detuned_freq, method="nearest").item()
print("Transmission at detuned frequency: ", detuned_transmission)

Transmission at detuned frequency:  0.4398850202560425

Define Input Power

We define the source. We first use scipy’s erf function to ramp up input power from 0 to \(3.95P_0\), hold this power, then add a pulse with peak power of \(20.85P_0\) to perform the switch, then go back down to \(3.95P_0\). The frequency of this power curve is detuned by \(\delta=2\sqrt{3}\).

# ramp up to CW
muramp = 3 / gamma  # central time for ramp up
sigmaramp = 1 / gamma  # time width for ramp up

PCW = 3.95 * P_0
# switching pulse
mupulse = 12.5 / gamma
sigmapulse = 1.5 / gamma
Ppeak = 20.85 * P_0


def power_curve(t):
    t /= gamma
    sig = np.sqrt(PCW) * (1 + scipy.special.erf((t - muramp) / (sigmaramp * np.sqrt(2)))) / 2
    sig += (np.sqrt(Ppeak) - np.sqrt(PCW)) * np.exp(-(((t - mupulse) / sigmapulse) ** 2) / 2)
    return sig


nt = 2000
tmax = 19 / gamma
times = np.linspace(0, tmax, nt)
pulse_run_time = tmax

fig, ax = plt.subplots(1, 1, figsize=(6, 5))
ax.set_xlabel(r"t$\gamma$")
ax.set_ylabel("$P/P_0$")

plt.plot(times * gamma, np.abs(power_curve(times * gamma)) ** 2 / P_0)
plt.plot(times * gamma, np.abs(power_curve(times * gamma)) ** 2 / P_0 / detuned_transmission)
plt.show()

We now create our source time with power amplitude given by this power curve, normalized with the transmission through the waveguide waveguide at our detuned frequency.

Normalization Note:

Notice that we also normalize by a factor of two. This is because, in Yanik et al., peak power is given, whereas the source power gives the cycle average power, which is 1/2 the peak power, so we need to multiply by 2 to get the peak power.

source_time = td.CustomSourceTime.from_values(
    freq0=detuned_freq,
    fwidth=fwidth,
    values=power_curve(times * gamma) / np.sqrt(2 * detuned_transmission),
    dt=times[1] - times[0],
)
nonlinear_source = gaussianBeam.updated_copy(
    center=[-source_distance * a, 0, 0], source_time=source_time
)
fig, ax = plt.subplots(1, 1, figsize=(5, 4))
source_time.plot(np.linspace(0, pulse_run_time, 20000), ax=ax)
plt.show()

Nonlinear Simulation

Now we add nonlinearity to the cavity, and input a pulse to switch between bistable states. The bistable states will have the same input power at the same (detuned from \(\omega_{res}\)) frequency, but will have very different transmissions through the waveguide.

To add nonlinearity, we compute \(\chi^{(3)}\) according to the formula:

\(\chi^{(3)}=\frac{n_2nRe(n)}{283}\)

and add this to the ellipse geometry we’ve already defined:

chi3 = n2 * n_rods * np.real(n_rods) / 283

cavityMedium = td.Medium(
    permittivity=n_rods**2,
    nonlinear_spec=td.NonlinearSpec(models=[td.NonlinearSusceptibility(chi3=chi3)], num_iters=10),
)

cavity_structure_nonlinear = td.Structure(
    geometry=cavity,
    medium=cavityMedium,
)

Define Further Monitors

We wish to examine the field profiles before and after the pulse at the different stable states, so we add FieldTimeMonitors at these times.

field_before_pulse = td.FieldTimeMonitor(
    center=[0, 0, 0],
    size=[td.inf, td.inf, 0],
    name="field before pulse",
    start=6 / gamma,
    stop=6 / gamma,
    fields=["Ez"],
)
field_after_pulse = td.FieldTimeMonitor(
    center=[0, 0, 0],
    size=[td.inf, td.inf, 0],
    name="field after pulse",
    start=19 / gamma,
    stop=19 / gamma,
    fields=["Ez"],
)
flux_t = td.FluxTimeMonitor(
    center=[flux_distance * a, 0, 0],
    size=[0, td.inf, td.inf],
    name="flux t",
    normal_dir="+",
)

03:10:39 UTC WARNING: The monitor 'interval' field was left as its default      
             value, which will set it to 1 internally. A value of 1 means that  
             the data will be sampled at every time step, which may potentially 
             produce more data than desired, depending on the use case. To      
             reduce data storage, one may downsample the data by setting        
             'interval > 1' or by choosing alternative 'start' and 'stop' values
             for the time sampling. If you intended to use the highest          
             resolution time sampling, you may suppress this warning by         
             explicitly setting 'interval=1' in the monitor.                    

Define Nonlinear Simulation

sim_nonlinear = sim_linear.updated_copy(
    run_time=pulse_run_time,
    medium=air,
    sources=[nonlinear_source],
    monitors=[field_before_pulse, field_after_pulse, flux_t],
    structures=[block, waveguide, cavity_air, cavity_structure_nonlinear, air_block, dbr_boundary],
)

             WARNING: 'normalize_index' 0 is a source with 'CustomSourceTime'   
             time dependence. Normalizing frequency-domain monitors by this     
             source is only meaningful if field decay occurs.                   

sim_nonlinear.plot(z=0)
plt.show()

sim_nonlinear_data = web.run(sim_nonlinear, task_name="bistable_pc")

03:10:40 UTC WARNING: 'normalize_index' 0 is a source with 'CustomSourceTime'   
             time dependence. Normalizing frequency-domain monitors by this     
             source is only meaningful if field decay occurs.                   

             WARNING: 'normalize_index' 0 is a source with 'CustomSourceTime'   
             time dependence. Normalizing frequency-domain monitors by this     
             source is only meaningful if field decay occurs.                   

             WARNING: 'normalize_index' 0 is a source with 'CustomSourceTime'   
             time dependence. Normalizing frequency-domain monitors by this     
             source is only meaningful if field decay occurs.                   

             WARNING: 'normalize_index' 0 is a source with 'CustomSourceTime'   
             time dependence. Normalizing frequency-domain monitors by this     
             source is only meaningful if field decay occurs.                   

             WARNING: 'normalize_index' 0 is a source with 'CustomSourceTime'   
             time dependence. Normalizing frequency-domain monitors by this     
             source is only meaningful if field decay occurs.                   

             WARNING: Simulation has 2.05e+06 time steps. The 'run_time' may be 
             unnecessarily large, unless there are very long-lived resonances.  

             Created task 'bistable_pc' with resource_id
             'fdve-5ad2939b-4c29-4e31-9ee5-90c75c80dffb' and task_type 'FDTD'.

             View task using web UI at
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-5ad2939b-4c2
             9-4e31-9ee5-90c75c80dffb'.

             Task folder: 'default'.

03:10:42 UTC Estimated FlexCredit cost: 0.321. 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.

03:10:43 UTC status = success

03:10:46 UTC Loading results from simulation_data.hdf5

             WARNING: 'normalize_index' 0 is a source with 'CustomSourceTime'   
             time dependence. Normalizing frequency-domain monitors by this     
             source is only meaningful if field decay occurs.                   

             WARNING: Warning messages were found in the solver log. For more   
             information, check 'SimulationData.log' or use                     
             'web.download_log(task_id)'.                                       

sim_nonlinear_data.plot_field(
    field_monitor_name="field before pulse", field_name="Ez", vmin=-10, vmax=10
)
plt.show()

sim_nonlinear_data.plot_field(
    field_monitor_name="field after pulse", field_name="Ez", vmin=-10, vmax=10
)
plt.show()

Analytic Solution

The analytic transmission we should expect to get is given by the formula in Yanik et al.:

\(\frac{dS_{ref}}{dt}=i\omega_{res}(1-\frac{1}{2Q}\frac{|S_{ref}|^2}{P_0})S_{ref}-\gamma S_{ref}-\gamma S_{in}\)

where \(P_{in}=|S_{in}|^2\) and \(P_{ref}=|S_{ref}|^2\) are the input power and reflected power, respectively.

Transmission

The transmission flux oscillates according to the source. In the figure from Yanik et al., the output results are the peak power with each optical period. We thus define a function that finds the peak power with each optical period:

def rolling_max(sig, freq, dt):
    interval = int(1 / (dt * freq))
    N = len(sig) // interval
    times = np.linspace(0, (len(sig) - interval) * dt, N - 1)
    sig_max = np.zeros(N - 1)
    for i, t in enumerate(times):
        start = int(t / dt)
        stop = start + interval
        sig_max[i] = np.max(sig[start:stop])
    return times, sig_max

We plot the analytically calculated transmission against the source. This transmission matches Figure 4 from Yanik et. al. Errors in our calculations are likely due to reflections in the boundary (parameters in the DBR boundary condition), meshing size, and cavity characteristic power calculation. Differences between Tidy3D calculations and the theoretical curve may also be due to the fact that the theoretical curve does not take higher harmonics due to nonlinearity into account.

times_nonlinear, nonlinear_transmission = rolling_max(
    sim_nonlinear_data["flux t"].flux.data, detuned_freq, sim_nonlinear.dt
)


# integrate a complex velocity f(x, t) using forward euler
def integrate(f, num_steps, T):
    dt = T / num_steps
    integral = np.zeros(num_steps, dtype=complex)
    time = 0
    for time_step in range(1, num_steps):
        x0 = integral[time_step - 1]
        integral[time_step] = x0 + dt * f(x0, time)
        time += dt
    return integral


interp_values = np.abs(source_time.amp_time(times_nonlinear)) * np.sqrt(2 * detuned_transmission)

power_input = scipy.interpolate.interp1d(times_nonlinear, interp_values)


def f(x, t):
    return 1j * gamma * (detuning - np.abs(x) ** 2 / P_0) * x - gamma * x - gamma * power_input(t)


nt = len(times_nonlinear)
tmax = times_nonlinear[-1]
dt = tmax / nt
res = integrate(f, num_steps=nt, T=tmax)

Sin = power_input(times_nonlinear)
Sref = res

Pin = np.abs(Sin) ** 2
Pout = np.abs(Sin + Sref) ** 2

fig, ax = plt.subplots(1, 1, figsize=(6, 5))
ax.set_xlabel(r"t$\gamma$")
ax.set_ylabel("$P/P_0$")

plt.plot(times_nonlinear * gamma, Pin / P_0, label="Input power")
plt.plot(times_nonlinear * gamma, Pout / P_0, label="Analytic")
plt.plot(times_nonlinear * gamma, nonlinear_transmission / P_0, label="Tidy3D")

ax.legend()
plt.show()