TE-like mode is not smooth when changing the geometry

[vorobus]

Hi, I am new to MPB and learning the basics and found some unclear behavior of TE mode frequency upon sweeping the geometry of the object with sharp edges.

I am trying to follow the algorithm from the paper

combined with this paper:

And want to simulate a bandgap of such a structure as a function of hole size:

To do this, I created the unit cell of interest (3D) with Prism as a waveguide and cylinder geometric object and ask for TE modes.

Searching for them, I found that when I gradually change the R of the holes, the mode frequency of the second even mode is having instabilities. Looking up the problem in the manual and documentation I found that this could come from the discontinuity of the epsilon, and could be tackled via the mesh_size parameter.

I wonder, what is the good practice with this? How the mesh_size parameter actually affects the obtained epsilon? Despite that, I made the mesh size 7,11, or 15, and it seems not removing the problem - I am still getting the jumps in the frequency when changing the R slowly. What could be the solution to this problem?

Where the X-axis is the R of the holes. Thank you in advance!

[stevengj]

If you are getting a 10% jump in the frequency (from 1.3 to 1.4) as you change the radius slightly, then it is too big to be a discretization effect. You must have some other problem.

e.g. a bug in how you are defining the geometry, or in how you are identifying the gap (you know that this structure only has a gap in the guided modes, right?), or something like that.

[vorobus]

Thank you, @stevengj!

Indeed, a 10% jump forced me to shift my direction of thoughts towards some bug search,

Just to clarify:

1. I am defining the geometry with just two objects - ellipsoid and prism. I will attach the code for that below.
2. I identify the mode frequency for the gap at the position of k_x = 0.5, at the end of the irreducible zone (and the modes are guided, i.e. they are below the light cone)

2D cylinder modes I first tried to see what happens if I consider just cylindrical glass waveguide and sweep the R of a cylinder and ask for the frequency of the modes at particular k = pi/a As expected I got a smooth picture: But I noticed that it could be a frustration between the modes if one solves only for 2 modes, but 2-nd and 3-rd are crossing.

I tried then to rerun my initial example of 3D waveguide photonic crystal consisting of the dielectric prism and elliptical vacuum hole in it, and solving for n_bands = 2 and n_bands = 3, and here is what I see:

When I solve for n_bands = 3, run_te solver

And here is where I use n_bands = 2, run_te solver (same resolution, and same mesh_size)

So, it is maybe related to the fact that when I ask only for 2 modes, the choice of the second mode is frustrating for the solver? It would be interesting to hear your opinion on that, Thank you, and with best regards, Vadim

Code for triangular waveguide nanobeam.py file:

  # Global parameters:
# # height of the prism
# h = 0.15         # units of um
# # width of the beam
# w = 0.3         # units of um
import meep as mp
from meep import mpb
import numpy as np
from matplotlib import pyplot as plt
from matplotlib import font_manager as fm
import matplotlib.colors as mcolors
import shutil as sh
fontpath = '/home/yy3/.fonts/times.ttf'
prop = fm.FontProperties(fname=fontpath, weight = 'ultralight')
plt.rc('text',usetex=False)
plt.rc('xtick', labelsize=8)
plt.rc('ytick', labelsize=8)
plt.rc('axes', labelsize=8)
# width as measured in inkscape
width = 3.487
height = width / 1.618
#%matplotlib inline
import os
from datetime import datetime
import scipy
from scipy import interpolate
from scipy import optimize
import numpy as np
import pandas as pd
import pickle
import time
import h5py

class analyze_geometry:
    c = 300.0 #um*THz speed of light
    def __init__(self, a = 1.0, h = 0.2, w = 0.4,res = 64, sc_y = 2, r = None, sc_z = 2, n_Host = 1.0, n_SiC = 2.6, k_point_x = None, plot = True):
        """
        h: height of the beam # um
        w: full top width of the beam
        sc_.. : width of the computational cell (in a)
        k_point_x: in units of $2\pi/a$
        """
        #recalibration to unit cell.
        if r is None:
            fiber = False
        else:
            fiber = True
        self.a = a # um
        self.h = h
        self.w = w
        self.fundmode = 'even'
        self.resolution = res # cells per a
        self.parity = mp.ODD_Y
        #self.h = h/self.a          # units of "a"
        #self.w = w/self.a          # units of "a"
        r = r

        self.sc_y = sc_y  # supercell width
        self.sc_z = sc_z  # supercell height
        self.Air = mp.Medium(index=n_Host)
        self.host_n = n_Host
        self.Glass = mp.Medium(index=1.45)

        self.geometry_lattice = mp.Lattice(size=mp.Vector3(0,sc_y,sc_z))
        self.SiC = mp.Medium(index=n_SiC)

        #geometry_lattice = mp.Lattice(size=mp.Vector3(period,4,4))
        self.vertices = [mp.Vector3(0,w/2,h/2), mp.Vector3(0,-w/2,h/2), mp.Vector3(0,0,-h/2)]
        self.geometry = [mp.Prism(vertices = self.vertices, axis = mp.Vector3(1,0,0),height = mp.inf, material=self.SiC)]
        if fiber:
            self.geometry.append(mp.Cylinder(center=mp.Vector3(0,0,r+h/2),axis = mp.Vector3(1,0,0),  radius=r, height=mp.inf, material=self.Glass))
        self.num_k = 20
        self.num_bands = 8

        self.k_min = 0.0
        self.k_max = 4*a
        self.k_points_all = mp.interpolate(self.num_k, [mp.Vector3(self.k_min), mp.Vector3(self.k_max)])
        self.k_points = [mp.Vector3(k_point_x,0,0)]
#         self.ms_all = mpb.ModeSolver(geometry_lattice=self.geometry_lattice,
#                             geometry=self.geometry,
#                             k_points=self.k_points_all,
#                             resolution=self.resolution,
#                             num_bands=self.num_bands)
        self.case_dir = os.path.join(os.getcwd(), "Nanobeam_data", datetime.now().strftime('%Y-%m-%d_%H-%M-%S'))

    def run(self,te, tm, *kwargs):

        self.ms_all = mpb.ModeSolver(geometry_lattice=self.geometry_lattice,
                                     default_material=mp.Medium(index=self.host_n),
                            geometry=self.geometry,
                            k_points=self.k_points_all,
                            resolution=self.resolution,
                            num_bands=self.num_bands)

        if te and tm:
            self.ms_all.run()
        elif te and not tm:
            self.ms_all.run_te()
        elif tm and not te:
            self.ms_all.run_tm()
        else:
            print('Error: nothing to run')

    def find_dk(self):
        """
        Finds delta from light cones...
        :return:
        """
        self.lambd = 0.916  # um V1 defect
        self.KS = np.linspace(self.k_min, self.k_max, len(self.k_points_all))
        self.dk = []
        for i in range(self.num_bands):  # number of bands
            x = self.KS / self.a
            y = self.ms_all.all_freqs.T[i] * self.c / self.a
            f = scipy.interpolate.interp1d(x,y,3)
            f0 = self.c / self.lambd
            g = lambda arg: abs(f(arg) - f0)
            res = scipy.optimize.minimize_scalar(fun=g, bounds=[self.k_min/self.a, self.k_max/self.a], method='Bounded', tol=1e-3)
            k_i = res.x
            KS_0 = f0*self.host_n / self.c  #: self.c / self.a / self.host_n * self.KS == f0
            self.dk.append(k_i - KS_0)
        #plt.plot(self.KS / self.a, self.c / self.a / self.host_n * self.KS, 'k--')
        #plt.hlines(xmin=0, xmax=4, y=self.c / self.lambd, colors='g', ls='-.')
        return self.dk

    def get_dispersion_curve(self, save= False, label_plot = None):

        f = plt.figure(dpi=150)
        self.lambd = 0.916 #um V1 defect
        self.KS = np.linspace(self.k_min,self.k_max,len(self.k_points_all))
        for i in range(self.num_bands):  # number of bands
            # plt.plot(self.KS/self.a,self.ms_all.all_freqs.T[i]*self.c/self.a, '-', lw = 0.5)
            plt.plot(self.KS/self.a, self.ms_all.all_freqs.T[i]/self.a, '-', lw=0.5)

        plt.plot(self.KS/self.a, self.host_n * self.KS / self.a, 'k', lw=0.8)
        # plt.plot(self.KS/self.a,-self.c/self.a/self.host_n*self.KS+self.c/self.a,'k',lw=0.8)
        # plt.hlines(xmin=0,xmax=4,y=self.c/self.lambd,colors='g', ls='-.')
        # plt.text(0, 360, r'$\lambda_{vac} = '+ str(1e3*self.lambd) + r' nm$',color = 'g',fontproperties = prop)
        # plt.text(0, 550, r'w = '+str(round(self.w,3))+r'$\mu m$   h = '+str(self.h)+r' $\mu m$',fontproperties = prop)
        # plt.text(0, 460, r'n0 = ' + str(round(self.host_n,3)), fontproperties=prop)

        plt.xticks(fontproperties=prop)
        plt.yticks(fontproperties=prop)
        plt.xlabel(r'$k_x a/2 \pi$', fontproperties=prop)
        plt.ylabel(r'Frequency $\omega a/2\pi c$', fontproperties=prop)
        #plt.hlines(xmin=0,xmax=4,y=self.c/self.lambd,colors='g', ls='-.')
        #plt.text(0, 360, r'$\lambda_{vac} = '+ str(1e3*self.lambd) + r' nm$',color = 'g',fontproperties = prop)
        #plt.text(0, 550, r'w = '+str(round(self.w,3))+r'$\mu m$   h = '+str(self.h)+r' $\mu m$',fontproperties = prop)
        #plt.text(0, 460, r'n0 = ' + str(round(self.host_n,3)), fontproperties=prop)

        #plt.xticks(fontproperties = prop)
        #plt.yticks(fontproperties = prop)
        #plt.xlabel(r'$k_x, 2\pi/\mu m$', fontproperties = prop)
        #plt.ylabel('Frequency c/a, THz',fontproperties = prop)
        #plt.xlim(0,2.0)
        #plt.ylim(0,800)
        f.set_size_inches(width, height)
        f.subplots_adjust(left=.20, bottom=.23, right=0.95, top=.98)
        plt.show()
        if save:
            self.save(f)
        plt.close()

    def estimate_mode_at_freq(self, freq=1, L = None, mode = 'abs', component = None):
        # estimate field profile at all bands.
        # modes: abs, real, imag
        # components: ex,ey,ez,hx,hy,hz,all,e,h
        self.all_fields = []
        self.all_fields_h = []
        self.fluxes = []
        self.mode_k_points = []
        for band in range(self.num_bands):

            if L is None:
                freq = freq
            else:
                freq = 1.0/self.a/L
            # find the P vector
            for parity in ['even', 'odd']:            
                k_point, pointing_x_real = self.find_k(band,freq, parity = parity)
                self.fluxes.append(pointing_x_real)
                #print(k_point)
                self.mode_k_points.append(k_point)

                time.sleep(0.1)

                fig, ax = plt.subplots()

                ax = plt.axes()
                shape = pointing_x_real.shape
                y = np.linspace(0,self.sc_y, shape[0])
                z = np.linspace(0,self.sc_y, shape[1])
                im = ax.pcolormesh(y,z,pointing_x_real, cmap = 'hot',linewidth=0,rasterized=True)
                im.set_edgecolor('face')
                ax.set_aspect("equal")
                ax.set_title('Pointing vector x real k ='+str(round(k_point[0],2)),fontproperties = prop)
                ax.set_xlabel(r'y coordinate, $\mu m$', fontproperties=prop)
                ax.set_ylabel(r'z coordinate, $\mu m$', fontproperties=prop)
                fig.set_size_inches(width, height*1.6)
                fig.subplots_adjust(left=.12, bottom=.22, right=0.98, top=.80)
                plt.colorbar(im)
                fig.savefig(os.path.join(self.case_dir, 'P_x'+str(band)+'_'+parity+'_'+'.pdf'))
                plt.close()

            if component is not None:
                self.ms = mpb.ModeSolver(geometry_lattice=self.geometry_lattice,
                                geometry=self.geometry,
                                k_points=[mp.Vector3(k_point[0],0,0)],
                                resolution=self.resolution,
                                num_bands=self.num_bands)
                self.ms.run()

                for fi in ['e','h']:
                    if fi in component:
                        for bi, base_vector in enumerate(["x","y","z"]):
                            f,ax = plt.subplots()

                            if fi == 'e':
                                fields = self.ms.get_efield(band+1,bloch_phase=True)
                            elif fi == 'h':
                                fields = self.ms.get_hfield(band+1,bloch_phase=True)

                            self.all_fields.append(fields)
                            d = np.real(fields[:,:,0,bi]).T

                            absmax = max([abs(d.min()),abs(d.max())])

                            norm = mcolors.DivergingNorm(vmin=-absmax, vmax = absmax, vcenter=0)
                            #plt.imshow(d, cmap=plt.cm.RdBu, norm=norm)
                            im = ax.pcolormesh(y,z,d, cmap = "RdBu",rasterized=True, norm = norm)
                            plt.colorbar(im)

                            plt.xticks(fontproperties=prop)
                            plt.yticks(fontproperties=prop)
                            #eps = self.ms.get_epsilon()
                            #ax.contour(eps[::-1,:].T, cmap='binary')

                            ax.set_xlabel(r'y coordinate, $\mu m$', fontproperties=prop)
                            ax.set_ylabel(r'z coordinate, $\mu m$', fontproperties=prop)
                            ax.set_aspect("equal")

                            f.set_size_inches(width, width)
                            ax.set_title('k = ' +str(k_point), fontproperties = prop)
                            f.subplots_adjust(left=.12, bottom=.22, right=0.91, top=.82)
                            f.savefig(os.path.join(self.case_dir,'band-'+str(band)+fi+'-'+base_vector+'-field'+'.pdf'))
                            plt.close()                        

    def find_k(self,band, freq, parity = None):

        f_mode = freq#1/0.86   # frequency corresponding to 0.86 um                                                                                                             
        band_min = 1+band
        band_max = 1+band
        kdir = mp.Vector3(1)
        tol = 1e-6
        kmag_guess = f_mode*2.0
        kmag_min = f_mode*0.1
        kmag_max = f_mode*4.0
        parity_def = 'odd'
        parity = parity_def if parity is None else parity
        self.parity = {'even':mp.EVEN_Y, 'odd':mp.ODD_Y}[parity]
        k = self.ms_all.find_k(self.parity, f_mode, band_min, band_max, kdir, tol, kmag_guess,
                        kmag_min, kmag_max,mpb.output_poynting_x)#, mpb.display_group_velocities)
        if self.parity == mp.ODD_Y:
            files = "flux.v.k01.b01.x.yodd.h5"
        elif self.parity == mp.EVEN_Y:
            files = "flux.v.k01.b01.x.yeven.h5"
        with h5py.File(files) as f:
            d = np.array(f['x.r']).T
        return k,d

    def save(self, fig):

        variables_to_save = ["h",
                            "w",
                            "a",
                            "resolution",
                            "k_points",
                            "ms_all.all_freqs",
                            "resolution",
                            "fundmode"]

        d = {}
        for var in variables_to_save:
            var2 = var.split('.')
            if len(var2) == 1:
                d[var] = self.__dict__[var]
            else:
                d[var] = self.__dict__[var2[0]].__dict__[var2[1]]
        try:
            os.makedirs(self.case_dir)
        except Exception as e:
            print(e)

        # 1. save the figure and the data
        folder_name = self.case_dir
        pickle.dump(d,open(os.path.join(folder_name, "variables.p"),'wb'))
        fig.savefig(os.path.join(folder_name, "plot.pdf"))
        sh.copy2(os.path.join("nanobeam.py"),os.path.join(folder_name, "_nanobeam.py"))
        #%notebook -e _backup.ipynb
        sh.copy2(os.path.join('Modes_of_the_waveguide_SiC.ipynb'),os.path.join(folder_name, "_backup.ipynb"))

my_analysis = analyze_geometry()

The instructions from jupyter notebook to reproduce the figures:

from importlib import reload
import nanobeam as nb; reload(nb)
from matplotlib import pyplot as plt
%matplotlib inline
import os
os.getcwd()
import numpy as np
l_te = 0.916 #nm
theta = 45 #deg
k = 10/7.0
res = 32
sc = 1.2317642832343687*0.9482403860729134
a = sc*0.319*l_te;print('a',a)
w = sc*0.58*l_te*(k)**0.5;print('w',w, w/a)
h = sc*(w/2)/np.tan(np.deg2rad(theta));print('h',h, h/a)
r = sc*0.087*l_te/a;print('r',r, r*a)
r1_max = sc*(k)**0.5*0.145*l_te/a;print('r1_max',r1_max, r1_max*a)
r1_min = sc*0.087*l_te/a;print('r1_min',r1_min, r1_min*a)

w*h*0.5

0.95/(a/res)
r1s = np.linspace(r1_min,r1_max,7)
for ms in [3,7]:
    plt.figure(figsize = (5,3), dpi=100)
    my_analysis_list2 = []
    for r1_I in r1s:
        time.sleep(1)
        print(r1_I)
        my_analysis_list2.append(nb.analyze_geometry(a=a, mesh_size = ms, 
                                                     num_k = 5, num_bands = 2, 
                                                     res = 15,n_SiC = 2.6, 
                                                     sc_x = 1.0, sc_y = 5.0, 
                                                     sc_z = 5.0, h=h/a, w = w/a, 
                                                     r = 2*r, r1 = 2*r1_I))
        ma = my_analysis_list2[-1]
        ma.run(tm=False,te=True,run=True)
        ma.get_dispersion_curve(save=True,label = 'te')
    lowhigh_tm = []
    for ma in my_analysis_list2:
        freqs_tm = ma.ms_te.all_freqs/ma.a
        lowhigh_tm.append(freqs_tm[:][-1])
    lowhigh_tm = np.array(lowhigh_tm)
    plt.title('Mesh_size:' + str(ms) + '_TE bandgap')
    plt.plot(r1s,lowhigh_tm.T[0],'.-')
    plt.plot(r1s,lowhigh_tm.T[1], '.-')
    #plt.plot(r1s,lowhigh_tm.T[2], '.-')
    plt.ylabel(r'frequency, $1/\lambda$')
    plt.xlabel('Radius, [a]')
    plt.hlines(xmin = r1s[0], xmax = r1s[-1], y = 1/0.916, color= 'r', ls = '--', label = '916 nm')
    plt.legend()

    lowhigh_tm.T[0][-1]*0.916