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.

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.

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.

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

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:

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

Visualize Simulation¶

Run Simulation¶

Calculate Linear Resonance¶

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

Compute Resonance Frequency and Q¶

We import Tidy3D's resonance finder

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

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$:

We can now define the characteristic power of the cavity:

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

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

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:

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}$.

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.

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:

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

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.

Define Nonlinear Simulation¶

Analytic Solution¶

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

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:

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.