Frequency-dependent correction factor

[Din-Y]

In issue #2503 I asked about phase in simple model consists of two homogeneous materials with different refractive indexes. I have a complex geometrical model (unit cell with two pillars) and planewave spread through the first part of homogeneous silicon material of substrate and the second part of two silicon pillars. In cross section of the second part the most of the volume is occupied by air. I confused about calculation of frequency-dependent correction factor exp(ikxL). Could you tell me please how to calculate the phase correction factor in heterogeneous layer? Thanks.

[stevengj]

It sounds like you still want exp(ikₓL), where kₓ(ω) is the propagation constant of the mode in the heterogeneous (constant-cross-section) layer. (Which can be computed by MPB or any other mode solver.)

[Din-Y]

My model has homogeneous layer of substrate and homogeneous layer of air above heterogeneous layer. The middle heterogeneous layer of two pillars and air. I applied MPB solver to find wavevector magnitudes in the heterogeneous layer via ModeSolver.find_k. Then I calculate total phase correction factor as

cmath.exp(1j * 2 * math.pi  * (kx_substrate * L_substrate + kx_pillars * L_pillars + kx_air * L_air))

for all homogeneous layers wavevector calculated as

kx = math.cos(theta) * n * fcen

n is refractive index, fcen is frequency. I guess that angle of wave incidence always is zero.
My results doesn’t converged with data obtained by other methods. Could you tell me please where I make mistakes?

Also I attached script of ModeSolver.find_k

import meep as mp
from meep import mpb
import numpy as np
from matplotlib import pyplot as plt
# %%

Si = mp.Medium(index=3.44) 
Air = mp.Medium(index=1)

wl_max = 12                                         # max wavelength 
wl_min = 8                                          # min wavelength 
f_min = 1/wl_max                                     # min frequency
f_max = 1/wl_min                                     # max frequency

resolution = 20

gp = 8                                              # period of meta_atoms
h = mp.inf

#define parameters of the first pillar
pillar_one_size = mp.Vector3(h, 2.4, 1.6)
#pillar_one_center = mp.Vector3(3.75, 1.27, 1.27)
pillar_one_center = mp.Vector3(0, 1.27, 1.27)

#define parameters of the second pillar
pillar_two_size = mp.Vector3(h, 1.4, 2)
#pillar_two_center = mp.Vector3(3.75, -1.55, -1.55)
pillar_two_center = mp.Vector3(0, -1.55, -1.55)

#define rotation of pillars
rotation_x = mp.Vector3(x=1),
rotation_y = mp.Vector3(y=1).rotate(mp.Vector3(x=1),-45.*np.pi/180)
rotation_z = mp.Vector3(z=1).rotate(mp.Vector3(x=1),-45.*np.pi/180)

geometry =[mp.Block(
    material=mp.Medium(index=3.44),
    size = pillar_one_size,
    center = pillar_one_center,
    e1=rotation_x,
    e2=rotation_y,
    e3=rotation_z
    ),
    mp.Block(
    material=mp.Medium(index=3.44),
    size=pillar_two_size,   
    center=pillar_two_center,
    e1=rotation_x,
    e2=rotation_y,
    e3=rotation_z
    )
    ]

k = mp.Vector3(0, 0, 0)                       #the boundaries are Bloch-periodic

sy = gp                                             # cell size on Y direction
sz = gp                                             # cell size on Z direction
geometry_lattice = mp.Lattice(size=mp.Vector3(0, sy, sz))

# %%

all_num_modes = []
k_vectors =[]
nfreqs = 40

frequencies = np.linspace(f_min, f_max, nfreqs)
for i in frequencies:
    t_nmode = np.floor((i-k.y)*gp)-np.ceil((-i-k.y)*gp)
    if t_nmode == 0:
        t_nmode += 1
    all_num_modes.append(t_nmode)


for j in range(len(all_num_modes)):
    ms = mpb.ModeSolver(
        geometry_lattice = geometry_lattice,
        geometry         = geometry,
        resolution       = resolution,
        k_points          = k,
        num_bands        = int(all_num_modes[j])
        )

for j in range(len(all_num_modes)):
    fcen = frequencies[j]
    num_modes = int(all_num_modes[j])
    print(fcen)
    print(num_modes)

    k = ms.find_k(
    mp.ODD_Y + mp.EVEN_Z, # mp.Ey
    fcen, # omega
    1, # band_min
    num_modes, # band_max
    mp.Vector3(1,0,0), # korig_and_kdir
    1e-4, # tol
    fcen * 3.44, # kmag_guess 
    fcen * 0.1, # kmag_min  
    fcen * 4, # kmag_max
    )

    k_vectors.append(k)

[Din-Y]

I can be mistaken, but a complex unit cell of grating requires to consider some effects and/or interactions between pillars and so on, to get correct value of phase via get_eigenmode_coefficients. Could you tell me please where I can make a mistake or do you know another approach to get phase?