Negative values in Purcell Enhancement of Cylinder metasurface.

[Fengkai-Wei]

I follow the example 'Planar Cavity with Lossless Metallic Walls' and tried to transfer it into my metasurface to simulate the Purcell Enhancement. Here is the general description of my Simulation Project.

Purcell Factor Simulation - Slice and Position sweep

1.Geometry:

1. Cylinder(Metastructure):

   • Height: 0.44 um
   • Center: (0,0,0.2) um
   • Radius: 0.40 um
   • Material: Silicon

2. Box1(Substrate):

   • Length = 1000 nm
   • Width = 1000 nm
   • Thickness = 500 nm
   • Center = (0,0,-250) nm
   • Material: Sapphire (Al2O3)

3. Box2(Reference Planar Surface)

   This geometric structure replaces Cylinder(Metasurface) in reference computation.

   • Length = 1000 nm
   • Width = 1000 nm
   • Thickness = 0.44 um
   • Center = (0,0,0.20) um
   • Material: Silicon

2.Sources:

All Dipoles are:

1. Position: Inside the center of cylindrical resonator.

   Sweeped Parameter:

   • Diople height: From 0 to 0.4 um, step size = 0.02
   • Dipole radial position: From 0 to 0.40 step size = 0.02

2. Magnitute = 1

3. Center frequency = 1/1.55, frequency width = 0.2*center frequency

4. Type: GaussianSource

   • Electric Dipole(ED_x, ED_y, ED_z)
   • Magnetic Dipole(MD_x, MD_y, MD_z)

3.Simulation Info:

1. FDTD size:(1,1,3.2) um
   • PML: z direction, thickness 0.1 um
2. Spaial resolution: 30 pixels/um
3. End Condition: 1000 Time Unit
4. Boundary Conditions:
   • x: peridoic
   • y: peridoic
   • z: perfect matched layer(PML), with thickness = 0.1 um
5. Materials to be used:
   • Crystalline Silicon(c-Si), defined by $\epsilon_r = 11.68$
   • Sapphire(Al2O3), defined by $\epsilon_r = 10.0$
6. Interested Frequency: 1/1.55

But the Purcell Factors I got are sometimes negative. Like in output lines:

on time step 57541 (time=959.017), 0.00247497 s/step 
Meep progress: 987.1333333333333/1000.0 = 98.7% done in 152.1s, 2.0s to go 
on time step 59256 (time=987.6), 0.00233249 s/step 
ldos0:, 0.6451612903225806, -25.815967412706385 
run 0 finished at t = 1000.0 (60000 timesteps)

-25 in ldos0 is obviously unphysical. I plotted the Purcell Enhancement sweep:

All Red Ns indicates either the ldos of metasurface or ldos of planar surface is negative. And both values are negative may also lead to positive Enhancement Factor...

And I was doubting whether the following factors may cause this weird performance:

1. Resolution = 30 is too low for a proper simulation in my model.
2. The resonator (Cylinder) is so small that the dipole may too close the interface of resonator.
3. Time = 1000 is not adequate for running, It should be replaced my examining the source decay.
4. Possibly my Boundary Condition is not correct... :-(
5. Other factors that I haven't been aware of

For your reference, here is my python script:

import meep as mp
import numpy as np
import pandas
import matplotlib.pyplot as plt
import math
from IPython.display import clear_output
from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets
from mpl_toolkits.mplot3d import axes3d
%matplotlib ipympl

pml_thickness = 0.1
resolution = 30
end_time = 1000
#end_condition = 
cell_length = 1
cell_width = 1
cell_height = 3 + 2*pml_thickness
cell_size = (cell_length,
            cell_width,
            cell_height)

pml_layers = [mp.PML(pml_thickness, direction=mp.Z)]

Sapphire_subs = mp.Block(mp.Vector3(mp.inf,mp.inf,0.5),
                center=mp.Vector3(0,0,-0.25),
                material=mp.Medium(epsilon = 10.0))


fcen = 1/1.55
fwidth = 0.2*fcen

Height = 0.44 #TBD
Radius = 0.44 #TBD
r_step = 0.02
h_step = 0.02
r_list = np.arange(0,Radius,r_step)
h_list = np.arange(0,Height,h_step)
print(r_list.shape,h_list.shape)

Si_ref = mp.Block(mp.Vector3(mp.inf,mp.inf,Height),
                              center=mp.Vector3(0,0,Height*0.5),
                              material=mp.Medium(epsilon = 11.68))

Si_meta = mp.Cylinder(radius = Radius,
                      height = Height,
                      center=mp.Vector3(0,0,Height*0.5),
                      material=mp.Medium(epsilon = 11.68))


class dipole(object):
    size = mp.Vector3(0,0,0)
    
    def __init__(self, dp_center, dp_type):
        self.dp_center = dp_center
        self.dp_type = dp_type
        self.dp_fcen = fcen
        self.dp_fwidth = fwidth
    def update_center(self,new_center):
        self.dp_center = new_center
    def update_type(self,new_type):
        self.dp_type = new_type
    def gen(self):
        dipole_out = mp.Source(mp.GaussianSource(self.dp_fcen, 
                                                 self.dp_fwidth),
                               component=self.dp_type,
                               center=self.dp_center,
                               size = self.size)
        return dipole_out


idx_h = 0

for h in h_list:
    idx_r = 0
    for r in r_list:
        dipole_list = [dipole(dp_center=mp.Vector3(0,r,h),dp_type = mp.Ex).gen(),
                       dipole(dp_center=mp.Vector3(0,r,h),dp_type = mp.Ey).gen(),
                       dipole(dp_center=mp.Vector3(0,r,h),dp_type = mp.Ez).gen(),
                       dipole(dp_center=mp.Vector3(0,r,h),dp_type = mp.Hx).gen(),
                       dipole(dp_center=mp.Vector3(0,r,h),dp_type = mp.Hy).gen(),
                       dipole(dp_center=mp.Vector3(0,r,h),dp_type = mp.Hz).gen(),]
        
        sim_Si = mp.Simulation(cell_size=cell_size,
                                    boundary_layers=pml_layers,
                                    geometry=[Si_meta, Sapphire_subs],
                                    sources=dipole_list,
                                    resolution=resolution)
        

        sim_Si_ref = mp.Simulation(cell_size=cell_size,
                                    boundary_layers=pml_layers,
                                    geometry=[Si_ref, Sapphire_subs],
                                    sources=dipole_list,
                                    resolution=resolution)
        
        sim_Si.run(mp.dft_ldos(fcen, 0, 1),
                            until = end_time)
        print('Reference Planar Surface is simulated.\n')
        
        sim_Si_ref.run(mp.dft_ldos(fcen,0,1),
                        until = end_time)
        print('Metasurface is simulated.\n')




        print('---------------------------')
        print(f'LDOS meta:{sim_Si.ldos_data[0]} \n')
        print('---------------------------')
        print(f'LDOS ref:{sim_Si_ref.ldos_data[0]} \n')
        print('---------------------------')
        print(idx_h,idx_r)
        print(Purcell_Enhancement.shape)
        Purcell_Enhancement[idx_h,idx_r] = sim_Si.ldos_data[0]/sim_Si_ref.ldos_data[0]
        print('Purcell Factor is simulated and stored.\n')
        idx_r += 1

    idx_h += 1

np.save("temp_Fp", Purcell_Enhancement)