Use of symmetries prevents to get results with harminv

[elaot]

Dear all,

I am currently working on the analyze of a photonic cavity that is composed of a waveguide of rectangular section, pierced with holes. At both extremities of the waveguide, the holes have the same dimensions and are equally spaced and they form Bragg mirrors ; between these mirrors, a defect is progressively introduced in the periodicity and the size of the holes, allowing the creation of a resonant cavity.

The waveguide is place onto a LiNbO3 crystal. As you can see on the figure, there are two mirror planes, respectively perpendicular to the x and the y axis. Nevertheless, when I place a Gaussian source at the center of the waveguide, and I indicate the symmetries to meep, the harminv functions I do compute are returned empty. To the contrary, I obtain coherent results when I do not indicate the symmetries.

The source is along Ey and I indicate the phase -1 for the Y mirror plane.

I would really appreciate any help.

You can find my code below :

`_
#!/usr/bin/env python3
import numpy as np
import meep as mp

Function that is building the desired device

def build_device(width, thickness,holes_length_0, holes_width_0, oblong_y,dpml):

# The different media
Si = mp.Medium(index=3.44)
vacuum = mp.Medium(index=1)
substrate = mp.Medium(index=2.21)

well=np.array([0.31808414, 0.31772006, 0.31669799, 0.31512321, 0.31310101,
       0.31073665, 0.30813542, 0.30540259, 0.30264344, 0.29996324,
       0.29746729, 0.29526085, 0.2934492 , 0.29213761, 0.29143138,
       0.29143138, 0.29213761, 0.2934492 , 0.29526085, 0.29746729,
       0.29996324, 0.30264344, 0.30540259, 0.30813542, 0.31073665,
       0.31310101, 0.31512321, 0.31669799, 0.31772006, 0.31808414])
holes_x =np.array([-4.25926992, -3.94154986, -3.62485187, -3.30972865, -2.99662765,
       -2.685891  , -2.37775558, -2.072353  , -1.76970956, -1.46974632,
       -1.17227903, -0.87701818, -0.58356899, -0.29143138,  0.        ,
        0.29143138,  0.58356899,  0.87701818,  1.17227903,  1.46974632,
        1.76970956,  2.072353  ,  2.37775558,  2.685891  ,  2.99662765,
        3.30972865,  3.62485187,  3.94154986,  4.25926992]) 
#holes_x contain the position of the center of holes in the defect

#the values composing the arrays well and holes_x result from a previous optimization that is not the subject of this 
#simulation

holes_length = holes_length_0 
holes_width = holes_width_0*(well/a)**(oblong_y)

holes_y = 0
# The beam itself : it is centered on (0,0,th/2), and extends infinitely in the x direction
blk = mp.Block(size=mp.Vector3(mp.inf, width, thickness),
                   center=mp.Vector3(0, 0, thickness/2),
                   material=Si)
   
geometry = [blk]

# Now we add the holes. We generate them as meep ellipsoids of
# inifinite z dimension

pad = 2                  # padding extension around the nanocavity


sx =( 2*a*n_mirror #Bragg mirrors at both sides of the central cavity
     + (holes_x[-1] + holes_length/2) - (holes_x[0] - holes_length/2) # size of the defect region that is forming the cavity
     + pad*2 + dpml*2) 

sy = width + pad*2 + dpml*2
sz = thickness + pad*2 + dpml*2

# First the defect holes
for x, w in zip(holes_x, holes_width):
    #zip() function creates an iterator that will aggregate elements from two or more iterables
    hole = mp.Ellipsoid(size=mp.Vector3(holes_length, w, mp.inf),  
                        center=mp.Vector3(x, holes_y, 0),
                        material=vacuum)
    geometry.append(hole)
# Then the mirror holes : before and after the defect
topmirror_x = holes_x[-1] + a
bottommirror_x = holes_x[0] - a


for m in range(n_mirror):
    # top mirror
    hole = mp.Ellipsoid(size=mp.Vector3(holes_length, holes_width[-1], mp.inf),
                        center=mp.Vector3((topmirror_x + m*a), holes_y, 0),
                        material=vacuum)
    geometry.append(hole)
    # bottom mirror
    hole = mp.Ellipsoid(size=mp.Vector3(holes_length, holes_width[0], mp.inf),
                        center=mp.Vector3((bottommirror_x - m*a), holes_y, 0),
                        material=vacuum)
    geometry.append(hole)

# The substrate, we need to keep it last, so that the holes are not punched through it
subs = mp.Block(size=mp.Vector3(mp.inf, mp.inf, sz/2),
                material=substrate)
subs += mp.Vector3(0, 0, -sz/4) # place the substrate in the bottom half
geometry.append(subs)


return(geometry,sx,sy,sz)

%% PARAMETER

a = 318.084145e-3 # Lattice parameter in µm
width = 594e-3 # Beam width in the y axis direction in µm
thickness = 250e-3 # Total beam thickness in m
holes_length_0 = 218.130412e-3 # Hole dimension in the x axis direction in µm
holes_width_0 = 158e-3 # Hole dimension in the y axis direction in µm
oblong_y = 2.23184998 # Oblongness factor. Positive values makes the holes stretch along the y axis
ndef = 14 # Number of holes in the defect region. The total number of defect holes is 2*p.ndef+1 because of symmetry
n_mirror = 12 # Number of holes in the mirro region.
pos=0.29143138/2 #middle betwteen the central hole and its neighbour

%% BUILDING OF THE DEVICE

resolution = 40
sym_x=True
sym_y=True

dpml = 1 # PML thickness for bounding conditions

geometry,sx,sy,sz=build_device(width, thickness, holes_length_0, holes_width_0, oblong_y,dpml)

cell = mp.Vector3(sx, sy, sz)
pml_layers = [mp.PML(dpml)]

Define the symmetries

symx = mp.Mirror(mp.X, phase=1)
symy = mp.Mirror(mp.Y, phase=-1)
sym=[]

if sym_x==True:
sym.append(symx)
if sym_y==True:
sym.append(symy)

%% SIMULATION

f1 = 1/1.520
f2 = 1/1.540
f_centre = (f1+f2)/2
f_width = f1-f2

harm0=mp.Harminv(mp.Ey, mp.Vector3(), f_centre, f_width)
harm1 = mp.Harminv(mp.Ey, mp.Vector3(pos, 0, thickness/2), f_centre, f_width)

src = [mp.Source(mp.GaussianSource(f_centre, fwidth=f_width), mp.Ey,
mp.Vector3(0,0,0))] #The center of the waveguide is at(0,0,thickness/2)

sim = mp.Simulation(cell_size=cell,
geometry=geometry,
boundary_layers=pml_layers,
sources=src,
symmetries=sym,
resolution=resolution)

sim.init_sim()

sim.run(mp.after_sources(harm0,harm1), until_after_sources=500)

print("Result of harminv0 computation \n",harm0.modes)
print("Result of harminv1 computation \n",harm1.modes)
_`

Sincerely,

Erell