The Difference Between a 2D Source and 3D Source For a 3D Periodic Waveguide

[rsmkim]

Hello everyone,

I have a question about the usage of a 2D vs 3D source for a 3D periodic structure in Meep. I understand that sources in Meep are current sources, but I have noticed that this source can take on a volume. I have created a 3D simulation that takes a single unit of a periodic silicon waveguide and excites it using an Eigenmode Source. I have been trying to match it to the mode profile I have made in COMSOL, but the fields are always a bit off.

For such a case, would it be better to have a 2D Eigenmode Source at the start of the period or have a 3D Eigenmode Source across a whole period of the waveguide and what are the differences of the two situations?

To me, it makes more sense to have a 3D source so that it can create the fields at each part of the mode profile for a 3D structure.
Thank you!

Here is my code:

import meep as mp
import numpy as np
import matplotlib.pyplot as plt
import math

c = 300000000
UM = 1000000

f = 234 #In THz
fsrc = f*THz/c / UM

wl = 1/fsrc

F = 0.7

Lx = 0.6
Ly = 0.25
Lz = 0.3

sx = Lx + 3*wl
sy = Ly + 3*wl
sz = Lz

bnum = 1    # band number of eigenmode

resolution = 40 # pixels/μm

cell = mp.Vector3(sx,sy,sz) # These sim bounds include the PMLs

pml_layers = [
                mp.PML(thickness=wl/2, direction=mp.X),
                mp.PML(thickness=wl/2, direction=mp.Y)
              ]

geometry = [
    mp.Block(
        size=mp.Vector3(Lx, Ly, Lz*F),
        center=mp.Vector3(),
        material=mp.Medium(epsilon=12.11),)
    ]

eig_kpoint = mp.Vector3(0, 0, 0.45/Lz) #in units of 2*pi/(unit length)
k_point = mp.Vector3(0, 0, 0.45) 

sources = [mp.EigenModeSource(src=mp.ContinuousSource(fsrc),
                            center=mp.Vector3(0,0,0),
                            size=mp.Vector3(Lx +  2*wl, Ly +  2*wl, Lz),
                            direction=mp.Z,
                            eig_kpoint=eig_kpoint,
                            eig_band=bnum,
                            eig_parity=mp.NO_PARITY,
                            eig_match_freq=False,
                            eig_lattice_center=mp.Vector3(0, 0, 0),
                            eig_lattice_size=mp.Vector3( 2*wl + Lx,  2*wl + Ly, Lz))] # This has to be larger than the continous source.

sim = mp.Simulation(cell_size=cell,
                    resolution=resolution,
                    boundary_layers=pml_layers,
                    sources=sources,
                    geometry=geometry,
                    k_point = k_point,
                    #ensure_periodicity=True,
                    )



sc_fr = mp.FluxRegion(
    center=mp.Vector3(0, Ly/2 + wl, 0), size=mp.Vector3(sx, 0, Lz)
)
sc = sim.add_flux(fsrc, 1, 1, sc_fr)

sim.run(until=100/fsrc)

# Plotting

eps_data = sim.get_array(center=mp.Vector3(), size=cell, component=mp.Dielectric)
ez_data = sim.get_array(center=mp.Vector3(), size=cell, component=mp.Ez)

eps_data = eps_data.real
ez_data = ez_data.real

idx_x = math.floor(ez_data.shape[0]/2)
idx_y = math.floor(ez_data.shape[1]/2 + 2)
idx_z = math.floor(ez_data.shape[2]/2)

plt.figure()
plt.plot((sy/2, -sy/2), (sx/2 - wl/2, sx/2- wl/2), color="white", linewidth=1.0) 
plt.plot((sy/2 - wl/2, sy/2 - wl/2), (sx/2, -sx/2), color="white", linewidth=1.0) 
plt.plot((-sy/2, sy/2), (-sx/2 + wl/2, -sx/2 + wl/2), color="white", linewidth=1.0) 
plt.plot((-sy/2 + wl/2, -sy/2 + wl/2), (sx/2, -sx/2), color="white", linewidth=1.0) 
plt.plot((Ly/2, -Ly/2), (Lx/2, Lx/2), color="white", linewidth=1.0) 
plt.plot((Ly/2, -Ly/2), (-Lx/2, -Lx/2), color="white", linewidth=1.0) 
plt.plot((-Ly/2, -Ly/2), (Lx/2, -Lx/2), color="white", linewidth=1.0) 
plt.plot((Ly/2, Ly/2), (Lx/2, -Lx/2), color="white", linewidth=1.0) 
plt.imshow(ez_data[:,:,idx_z], extent=[-sy/2, sy/2, -sx/2, sx/2])
plt.colorbar()
plt.savefig('Field xy Plane')

plt.figure()
plt.plot((F*Lz/2, -F*Lz/2), (Ly/2, Ly/2), color="white", linewidth=1.0) 
plt.plot((F*Lz/2, -F*Lz/2), (-Ly/2, -Ly/2), color="white", linewidth=1.0) 
plt.plot((F*Lz/2, F*Lz/2), (Ly/2, -Ly/2), color="white", linewidth=1.0) 
plt.plot((-F*Lz/2, -F*Lz/2), (Ly/2, -Ly/2), color="white", linewidth=1.0) 
plt.plot((-Lz/2, Lz/2), (sy/2 - wl/2, sy/2 - wl/2), color="white", linewidth=1.0) 
plt.plot((-Lz/2, Lz/2), (-sy/2 + wl/2, -sy/2 + wl/2), color="white", linewidth=1.0) 
plt.imshow(ez_data[idx_x,:,:], extent=[-sz/2, sz/2, -sy/2, sy/2])
plt.colorbar()
plt.savefig('Field yz Plane')

plt.figure()
plt.plot((F*Lz/2, -F*Lz/2), (Lx/2, Lx/2), color="white", linewidth=1.0) 
plt.plot((F*Lz/2, -F*Lz/2), (-Lx/2, -Lx/2), color="white", linewidth=1.0) 
plt.plot((F*Lz/2, F*Lz/2), (Lx/2, -Lx/2), color="white", linewidth=1.0) 
plt.plot((-F*Lz/2, -F*Lz/2), (Lx/2, -Lx/2), color="white", linewidth=1.0) 
plt.plot((-Lz/2, Lz/2), (sx/2 - wl/2, sx/2 - wl/2), color="white", linewidth=1.0) 
plt.plot((-Lz/2, Lz/2), (-sx/2 + wl/2, -sx/2 + wl/2), color="white", linewidth=1.0) 
plt.imshow(ez_data[:,idx_y,:], extent=[-sz/2, sz/2, -sx/2, sx/2])
plt.colorbar()
plt.savefig('Field xz Plane')

[stevengj]

It's the difference between:

• 2d source: the mode of a cross-section as if that cross-section were infinitely extended
• 3d source: the Bloch mode of the periodic unit cell. This can be very different if you actually have a periodic waveguide, and is the right source to launch a Bloch mode in that waveguide. e.g. for waveguides like those in chapter 7 of our textbook. (However, the wavevector k means something different in this case because it is modulo 2π/a where a is the period.)
  

(I'm not sure which mode you are solving for in Comsol, the mode of the cross-section or the Bloch mode?)

I'm not sure what you are expecting your simulation to show, because you are trying to fix both $k_z$ and the frequency. But what if the mode at that $k_z$ doesn't match the frequency fsrc? Shouldn't you look at the frequency of the mode, if you want to see the mode? (With eig_match_freq=False, you are computing the mode at the specified $k_z$, which is probably at a completely different frequency.) Also, if you want to just look at a mode, which is what it seems like since you are just looking at one unit cell, it would be easier to use MPB.

PS. Your flux region makes no sense to me. Your mode is propagating in the z direction, but you are computing flux in the y direction?

[rsmkim]

Hi, thank you for the response.

In COMSOL, I have been solving for the bloch mode. This gives a pretty good estimate for the frequency and wavevector of the mode in MEEP. (There is a slight difference because of different meshing and different solvers, but they are usually about 2-4 THz off at a given wavevector.)

To clarify the weird inconsistencies of my code:

• This has been my test code where I have been messing around with the different settings. This is why it looks like I am both setting the wavevector and the frequency. I have the eig_match_freq set to false in this current example because I was playing around with setting the wavevector instead.
• The flux region was also something I was playing around with, I should have just deleted it in the post here for clarity. I was trying to look if there was scattering from when I was using a 2D source.

I have been using MEEP rather than MPB because I ultimately want to use the adjoint solver to inverse design a periodic waveguide where the fields are concentrated above the center of the waveguide. I haven't looked too deep yet at the adjoint solver because I am focusing on learning the basics of MEEP first. Is the adjoint solver able to optimize using the fields calculated from MPB?

Edit:
I've looked moreinto MEEP's adjoint solver and it does not seem I can use MPB for the adjoint solver.