Catalog
  1. 1. Meep
    1. 1.1. 介绍
    2. 1.2. 波导模式求解
    3. 1.3. 求解耦合效率
    4. 1.4. 收敛性测试
    5. 1.5. 与Lumerical和FDTD对比
Meep教程(3):Meep的模式投影方法计算偶极子耦合到光波导的效率

Meep

介绍

因为实验室需要,我的研究生生涯相当一部分就是在计算这个耦合效率,偶极子放在波导附近时,其一部分能量会耦合到波导中去,偶极子可以是任意的荧光分子、量子点等,我们希望提高这个耦合效率,其定义为耦合到波导部分的功率比上总的辐射功率

其中$p_{guide}$是偶极子耦合到波导的功率,需要在传播方向上将总场与波导模式场进行投影得到,而$p_{total}$则是偶极子总的辐射功率,需要计算包裹偶极子的一个封闭面上的坡印廷矢量的积分。这是最一般最有效的方法,其实还有基于模式分析的半解析方法来计算,我这里就不在详细讨论,本笔记目的只是为了计算出正确的耦合效率。

波导模式求解

这次的结构要复杂一些,如下图所示

我们还是需要先求解一下对应的波导模式,和上篇笔记的做法一样,只是结构更加复杂了,设置结构的代码如下

1
import meep as mp
2
from meep import mpb
3
import numpy as np
4
from matplotlib import pyplot as plt
5
6
omega = 1/0.78
7
# index of material
8
index_SiN=2.02
9
index_SiO2=1.46
10
index_PVA=1.5
11
index_An=1.8
12
# The size of the guide
13
width=0.65
14
height=0.25
15
16
height_An=0.2
17
height_PVA=0.03
18
19
# The size of simulation domain
20
width_cal=2.5
21
length_cal=2.5
22
height_cal=2.5
23
24
# Thinkness of PML
25
t_pml=0.4
26
27
# The size of box to obatin the dipole power
28
size_x=0.3
29
size_y=0.3
30
size_z=0.3
31
32
# The position of dipole
33
pos_x=0
34
pos_y=0
35
pos_z=height+0.005
36
37
resolution=50
38
numx=int(width_cal*resolution)
39
numy=int(height_cal*resolution)
40
x_list=np.linspace(-width_cal/2,width_cal/2,numx)
41
y_list=np.linspace(-height_cal/2,height_cal/2,numy)
42
x_grid,y_grid=np.meshgrid(x_list,y_list)
43
# %%
44
cell_size = mp.Vector3(length_cal,width_cal,height_cal)
45
46
pml_layers = [mp.Absorber(thickness=t_pml,direction=mp.Y),
47
                   mp.PML(thickness=t_pml,direction=mp.Z),
48
                   mp.PML(thickness=t_pml,direction=mp.X)]
49
50
# Y dipole 
51
source = [mp.Source(mp.ContinuousSource(frequency= omega), 
52
            component = mp.Ex, center=mp.Vector3(0, 0, 0.255) )]
53
54
geometry = [mp.Block(material=mp.Medium(epsilon=index_SiN**2), 
55
                     size=mp.Vector3(mp.inf,width, height),
56
                     center=mp.Vector3(0, 0, height/2)),
57
58
            mp.Block(material=mp.Medium(epsilon=index_SiO2**2),
59
                     size=mp.Vector3(mp.inf,mp.inf, height_cal/2),
60
                     center=mp.Vector3(0, 0, -height_cal/4)),
61
62
            mp.Block(material=mp.Medium(epsilon=index_An**2),
63
                     size=mp.Vector3(mp.inf,mp.inf, height_An),
64
                     center=mp.Vector3(0, 0, height+height_An/2)),
65
            mp.Block(material=mp.Medium(epsilon=index_PVA**2),
66
                     size=mp.Vector3(mp.inf,mp.inf, height_PVA),
67
                     center=mp.Vector3(0, 0, height+height_An+height_PVA/2))
68
                     ]

然后是计算波导模式和绘图

1
geometry_lattice = mp.Lattice(size=mp.Vector3(0, width_cal, height_cal))
2
num_modes = 1
3
ms = mpb.ModeSolver(
4
    geometry_lattice = geometry_lattice,
5
    geometry         = geometry,
6
    resolution       = resolution,
7
    num_bands        = num_modes)
8
# ms.init_sim()
9
# %%
10
E = []
11
store_fields = lambda ms_temp, mode: E.append(ms_temp.get_efield(which_band=mode,bloch_phase=False))    
12
13
mode=1
14
k = ms.find_k(
15
    mp.NO_PARITY, # p
16
    omega, # omega
17
    mode, # band_min
18
    mode, # band_max
19
    mp.Vector3(1,0,0), # korig_and_kdir
20
    1e-4, # tol
21
    omega * 2, # kmag_guess
22
    omega * 0.1, # kmag_min
23
    omega * 4, # kmag_max
24
    store_fields # band_funcs
25
    )
26
# %%
27
neff=k[0]/omega
28
eps = ms.get_epsilon()
29
eps_arr=np.transpose(np.array(eps))
30
31
# %%
32
plt.figure(figsize=(12,4))
33
# Plot the E fields
34
mode=0
35
Ex=np.squeeze(E[mode][:,:,0,0]).transpose()
36
Ey=np.squeeze(E[mode][:,:,0,1]).transpose()
37
Ez=np.squeeze(E[mode][:,:,0,2]).transpose()
38
39
print('Current band: {}'.format(mode+1))
40
plt.subplot(1,3,1)
41
# plt.pcolormesh(x_grid,y_grid,eps.transpose(), cmap='binary')
42
plt.pcolormesh(x_grid,y_grid,np.abs(Ex), cmap='jet', alpha=0.9)
43
plt.axis('off')
44
plt.title('Ex')
45
plt.subplot(1,3,2)
46
# plt.pcolormesh(x_grid,y_grid,eps.transpose(), cmap='binary')
47
plt.pcolormesh(x_grid,y_grid,np.abs(Ey), cmap='jet', alpha=0.9)
48
plt.axis('off')
49
plt.title('Ey')
50
plt.subplot(1,3,3)
51
# plt.pcolormesh(x_grid,y_grid,eps.transpose(), cmap='binary')
52
plt.pcolormesh(x_grid,y_grid,np.abs(Ez), cmap='jet', alpha=0.9)
53
plt.axis('off')
54
plt.title('Ez')
55
plt.tight_layout()
56
plt.show()
57
plt.savefig('SimulatedField_onChip.png')

其模场分布如下

其折射率平方的分布绘制出来,便于我们进一步检查结构是否准确

1
fig1=plt.figure()
2
plt.pcolormesh(x_grid,y_grid,eps_arr,  cmap='jet')
3
plt.colorbar()
4
plt.savefig('StructuresonChip.png')
5
plt.show()

其分布如下

因此我们的结构设置是没有问题的。

求解耦合效率

我们需要设置一个封闭的面包裹住偶极子,将表面的功率密度积分,得到总功率,还需要设置两个模式投影,其实只需要一个,因为结构是对称的

1
sim = mp.Simulation(resolution=resolution,
2
3
 cell_size=cell_size,
4
5
 boundary_layers=pml_layers,
6
7
 geometry=geometry,
8
9
 sources=source,
10
11
 eps_averaging= False,
12
13
 Courant = 0.5)
14
15
 # 模式投影面设置
16
17
box_1 = sim.add_mode_monitor(omega, 0, 1,
18
19
mp.FluxRegion(center=mp.Vector3(0.2, 0, 0),
20
21
 size=mp.Vector3(0,width_cal,height_cal))) 
22
23
box_2 = sim.add_mode_monitor(omega, 0, 1,
24
25
mp.FluxRegion(center=mp.Vector3(-0.2, 0, 0),
26
27
 size=mp.Vector3(0,width_cal,height_cal))) 
28
29
  
30
31
# 总积分面设置
32
33
box_x1 = sim.add_flux(omega, 0, 1,
34
35
mp.FluxRegion(center=mp.Vector3(pos_x-size_x/2,pos_y,pos_z),size=mp.Vector3(0,size_y,size_z)))
36
37
box_x2 = sim.add_flux(omega, 0, 1,
38
39
mp.FluxRegion(center=mp.Vector3(pos_x+size_x/2,pos_y,pos_z),size=mp.Vector3(0,size_y,size_z)))
40
41
  
42
43
box_y1 = sim.add_flux(omega, 0, 1,
44
45
mp.FluxRegion(center=mp.Vector3(pos_x,pos_y-size_y/2,pos_z),size=mp.Vector3(size_x,0,size_z)))
46
47
box_y2 = sim.add_flux(omega, 0, 1,
48
49
mp.FluxRegion(center=mp.Vector3(pos_x,pos_y+size_y/2,pos_z),size=mp.Vector3(size_x,0,size_z)))
50
51
  
52
53
box_z1 = sim.add_flux(omega, 0, 1,
54
55
mp.FluxRegion(center=mp.Vector3(pos_x,pos_y,pos_z-size_z/2),size=mp.Vector3(size_x,size_y,0)))
56
57
box_z2 = sim.add_flux(omega, 0, 1,
58
59
mp.FluxRegion(center=mp.Vector3(pos_x,pos_y,pos_z+size_z/2),size=mp.Vector3(size_x,size_y,0)))
60
61
  
62
63
# %%

然后就可以计算然后输出结果了

1
# sim.run(until_after_sources=
2
3
# mp.stop_when_fields_decayed(20, mp.Ey, pt, 1e-3))
4
5
# sim.run(until_after_sources=20)
6
7
sim.run(until=20)
8
9
# %%
10
11
x1 = mp.get_fluxes(box_x1)
12
13
x2 = mp.get_fluxes(box_x2)
14
15
y1 = mp.get_fluxes(box_y1)
16
17
y2 = mp.get_fluxes(box_y2)
18
19
z1 = mp.get_fluxes(box_z1)
20
21
z2 = mp.get_fluxes(box_z2)
22
23
  
24
25
mode1 = sim.get_eigenmode_coefficients(box_1,
26
27
 [1], eig_parity=mp.NO_PARITY)
28
29
mode2 = sim.get_eigenmode_coefficients(box_2,
30
31
 [1], eig_parity=mp.NO_PARITY)
32
33
  
34
35
# %%
36
37
ptotal = -x1[0]+x2[0]-y1[0]+y2[0]-z1[0]+z2[0]
38
39
pmode_1 = abs(mode1.alpha[0,0,0])**2
40
41
pmode_2 = abs(mode2.alpha[0,0,1])**2
42
43
effic_1=pmode_1/ptotal
44
45
effic_2=pmode_2/ptotal
46
47
neff1 = mode1.kdom[0].x/omega
48
49
neff2 = mode2.kdom[0].x/omega
50
51
print(neff1)
52
53
print(neff2)
54
55
print("efficency to position direction:, {} ".format(effic_1))
56
57
print("efficency to negative direction:, {} ".format(effic_2))

收敛性测试

为了检查收敛性,我还改变了不同的分辨率(resolution),看结果是否准确,
代码地址
其结果如下

可以看见其实只需要50的分辨率就可以了,后面结果变化的很小。

与Lumerical和FDTD对比

另外我又扫描了偶极子的竖直高度和Lumerical FDTD以及COMSOL的解析结果对比,代码地址

其结果如下


可以看见Meep的Pguide的结果与COMSOL理论符合的非常好,与LumeircalFDTDSlutions不那么好,这与LumeircalFDTDSlutions归一化不好有关,而耦合效率吧则是完全对上的。

以上就是计算的全过程了。如果对你有帮助记得点赞~

Author: Knifelee
Link: https://knifelees3.github.io/2021/07/11/C_%E6%95%99%E7%A8%8B_%E7%94%A8Meep%E8%AE%A1%E7%AE%97%E5%81%B6%E6%9E%81%E5%AD%90%E8%80%A6%E5%90%88%E5%88%B0%E5%85%89%E6%B3%A2%E5%AF%BC%E6%A8%A1%E5%BC%8F%E7%9A%84%E6%95%88%E7%8E%87/
Copyright Notice: All articles in this blog are licensed under CC BY-NC-SA 4.0 unless stating additionally.
Donate
  • 微信
  • 支付寶

Comment