For most waveguide that are uniform in the direction of propagation, such as fibers, ridge waveguide, it is easy to find the dispersion relations using Lumerical's FDE solver. For example, see Dispersion Analysis and Waveguides. For waveguides with periodic patterning in the direction of propagation, we can use the FDTD solver to calculate the dispersion relations and bandstructure, as shown in Waveguides 3D. In this section, we will use FDE and FDTD method to analyze the dispersion of waveguides that are uniform in the direction of propagation. For simple material systems, the FDE solver is the best solver for studying uniform waveguides, but it is not capable of simulating waveguides composed of some complex materials such as nonlinear materials. For these advanced material systems, the FDTD method must be used.
In this example, we start by demonstrating the simulation methodology for a simple uniform waveguide composed of a simple linear, isotropic material. Next, we repeat the calculation using a diagonal anisotropic material. Finally, we repeat the calculation again using a magneto-optical material (non-diagonal anisotropic material).
Simulation setup
A screenshot of the FDE yz view to show the waveguide cross section
Isotropic waveguide
The isotropic waveguide is composed of a rectangular core (2.6 μm wide in y and 2 μm thick in z, refractive index 1.5), on a substrate of refractive index 1.4.
Anisotropic waveguide
We now set the object “core” to have an anisotropic refractive index of “1.5;1.49;1.51”. This will break the degeneracy of the TE and TM modes.
Magneto optical waveguide
Magneto optical material has an anisotropic permittivity with off-diagonal components. When light passes through such material, the polarization plane will rotate, and the left- and right- rotated circular polarization will propagate in different speed. Such anisotropic permittivity can be diagonalized (please refer to anisotropic materials). We enable the grid attribute called “faraday_effect” and make sure that the object “core” has the "grid attribute name" property also set to faraday_attribute. The grid attribute faraday_effect applies a complex unitary transformation of that was setup with the following script:
# Define U
# EL = (-Ey + 1i*Ez)/sqrt(2);
# ER = (-Ey - 1i*Ez)/sqrt(2);
# E_circular = U * E_cartesian, therefore:
U = [ sqrt(2), 0, 0;
0, -1, 1i;
0, -1, -1i] / sqrt(2);
# Set U as transformation matrix
select("faraday_attribute");
set("U",U);
FDE Simulation approach
The FDE solver supports material anisotropy through the addition of matrix grid attributes. This example simulates a rib waveguide with a slightly anisotropic core. The slight anisotropy breaks the degeneracy of the first two modes (TE and TM). This example is about computing the dispersion relation of those modes, and the generated modes have elliptical (almost circular) polarization as expected. The MODE simulation (nr_waveguide2d.lms) and script file (nr_waveguide2d.lsf) that generates the dispersion relation (kx vs k0) for the anisotropic case.
Isotropic waveguide
Anisotropic waveguide
Magneto optical waveguide
FDTD Simulation approach
The waveguide has been drawn up in the simulation file nr_waveguide3d.fsp, which is a simple rib waveguide.
Suppose the propagation is along x, we can use a "x span" of 50 nm or so and a override mesh with dx = 50 nm to have only one cell in x. This means we have approximately lambda/dx=1.55microns/1.5/50nm ~ 21 points per wavelength in the direction of propagation. Such small mesh size makes the additional grid dispersion small (see details on grid dispersion in FDTD).
If we have no idea of the mode profiles and want to extract the full bandstructure information then a randomized dipole cloud should be used as a source, with a bandwidth that covers the full range of frequencies we want to consider (please refer to rectangular photonic crystal bandstructure example for more details). For the current example, we know that the waveguide mode we are looking for is a small perturbation to the waveguide system (for example for some nonlinear or anisotropy effects), then we can use a mode source that starts with the solution from Mode Source. This will converge much faster than a randomized dipole cloud. However, the user should be cautious because some modes may be completely eliminated by this approach, particularly if they have a different symmetry from the mode used to excite the system. This waveguide supports two almost degenerate TE and TM modes. With the Mode Source, we can select both the TE and TM modes and then track these modes in a wavelength sweep of 1.5μm to 1.6μm. Such sources will be approximately polarized at 45 degrees in the waveguide and will excite both modes equally. We cannot use only one mode source because the modes are orthogonal. We use the simulation band from 1.4 to 1.7 μm. The value of beta or effective index vs wavelength can then be calculated through sweeping kx. For a fast demonstration, we only sweep 5 points from kx = 5.6 e6 to 6.0 e6 rad/m.
For periodically patterned waveguides, such as the line-defect photonic crystal waveguide, we simulate one unit cell along its propagation direction. For any uniform waveguide, any length of the waveguide can be considered as "periodic" along its propagation direction. Therefore, for these waveguides, we can also simulate one "unit cell". In this case, the value of the period is "arbitrary". For the sake of simulation efficiency, we use one mesh along the propagation direction, eg., the span is dx if along x axis. This means that we have only one single FDTD Yee cell in this direction. The result is to push the edge of the Brillouin zone to k=pi/dx.
To reduce the grid dispersion, this mesh size dx should be typically less than lambda/neff/10 where neff is the largest effective index we want to study. In the case of index guided structures, we can set dx.
Isotropic waveguide
We first simulate the regular waveguide without anisotropy or nonlinear property. Open the nr_waveguide3d.lsf file, and run the lsf file to sweep different k ( set in the sweep group). Then the effective index can be calculated as k/k0 where k0 is the wavenumber in vacuum. The beta and the spectrum vs. wavelength are plotted as follows:
Propagation constant
Resonant spectrum
It clearly shows that there is only one resonant peak at each k due to the mode degeneration.
Anisotropic waveguide
We now set the object “core” to have an anisotropic refractive index of “1.5;1.49;1.51”. This will break the degeneracy of the TE and TM modes. After re-running the sweep and the script file, we see the following figures showing that the degeneracy is indeed broken. At each k, there are two resonate frequencies corresponding to the two modes.
Propagation constants of the two modes
Effctive index of two modes
Amplitude of Ey and Ez
Resonant spectrum
It is also useful to run the simulation with an intermediate value of kx = 5.9e6 rad/m. We can then plot the E field in the time monitor “waveguide center”, and particularly the values of Ey and Ez as shown above (bottom left). We can see that while we inject both modes, there is no coupling between them which makes sense because the two modes are orthogonal.
Magneto optical waveguide
Now we have a waveguide made of a material with left and right circularly polarized light which have indices of 1.49 and 1.51. This means that our initially injected light which is approximately 45 degree polarized should rotate as the simulation runs. Additionally, we can extract the neff vs lambda using the same approach for both of the resulting modes. These results are
Propagation constants of the two modes
Effctive index of two modes
Amplitude of Ey and Ez
Resonant spectrum
We can again run a single simulation with kx = 5.9e-6 m-1 and plot |Ey| and |Ez| from the time monitor waveguide_center. We indeed see that the polarization rotates with time (above, bottom left).Finally, if we want to extract the mode profiles for the 2 resulting modes, we need to add a frequency domain monitor (called "mode profiles") and ensure that the wavelengths chosen correspond to the peak values of the 2 modes. To find the exact frequency, we can have the time monitor apodized starting from 300fs to 900fs. Run the simulation and using the following script to extract the frequency:
t=getdata("waveguide center","t");
Ey=getdata("waveguide center","Ey");
f=linspace(190,197,100)*1e12;
eey=czt(pinch((Ey)),t,f*2*pi);
plot(c/f*1e6,abs(eey));
plot(f*1e-12,abs(eey));
n=findpeaks(abs(eey),2);
?c/f(n)*1e6;
For the profile monitor, we also need to use the same apodization, and record the profiles for the found wavelengths or frequencies. Please note that these mode profiles look very similar when considering |E|. Looking at a line plot of the phase through the center of the mode, we can see that the phase of Ey and Ez is advanced or retarded by 90 degrees, corresponding mainly to a right or left circularly polarized mode.
Similarly, you may simulate waveguide with nonlinear material.