FDTD学习笔记|Lumerical官方案例、FDTD时域有限差分法仿真学习（六）——等离子体超材料吸收器（Plasmonic metamaterial absorber）仿真|时域有限差分|算法|FDTD

一、超材料模拟方法
本节主要涉及模拟人造“原子”，例如金属线对和裂环，可用于创建不寻常的有效体积特性，例如负折射率。我们还将讨论有效参数提取的细节。有关创建大块电磁介质的示例，请参阅大块超材料示例。
1.1 模拟设置提示
为确保以最有效的方式设置超材料模拟，请考虑将以下一些技术应用于您的模拟。这些技巧对于在低频（相对于光频率）下运行的设备来说是最重要的。
1.1.1 完美的金属近似
对于低频下的 3D 对象，金属以“理想”方式表现（100% 反射，0% 吸收），使用完美电导体 (PEC) 材料模型，而不是采样 3D 数据（ Sampled 3D data）、导电 3D（Conductive 3D ）或其他材料模型。 PEC 模型是数值效率最高的选项，因为它不需要那么精细的模拟网格。
如果您的设备在金属不理想的频率下运行（即在可见范围内），则不应使用 PEC 模型。
这些材料模型在添加材料的时候能够看到，如图所示。

文章图片

1.1.2 薄层
与波长相比，一些超材料具有极薄的层。可以使用 2D 矩形或 2D 多边形对象来表示薄层，并使用导电材料模型指定材料。
如果使用 3D 对象来表示薄层，为了获得准确的模拟结果，至少有几个网格单元来解析层的厚度是很重要的。当层非常薄时，这需要极小的网格，这会显着增加模拟的总内存和时间要求。在这种情况下，可以在模拟中使用更厚的层。例如，如果该层实际上是 1/1000 波长厚，则将您的模拟设置为 1/10 或 1/50 的层厚。这允许您使用更大的网格尺寸，而不会显著损失精度。
1.1.3 检查低频模拟的默认设置
FDTD 中的默认设置是针对光频率设置的。如果您在低得多的频率下工作，则某些默认值将不正确。特别是，检查以下设置：
单位 - 在设置菜单中。您可能想要更改单位设置。例如，从 THz 到 GHz。
模拟时间 - 默认值 1000fs 对于低频模拟来说不够长。在源脉冲之后，场必须有足够的时间衰减回零。
网格大小 - 在大多数情况下，网格大小将默认为合理的值，即使对于非常低频率的模拟也是如此。但是，在某些情况下（例如从不同的模拟文件复制网格覆盖区域），最终可能会导致网格尺寸过大或过小。
1.1.4 使用periodic（或 Bloch）边界条件
大多数超材料都是周期性的，这意味着您只需要模拟一个晶胞。使用周期性（如果源垂直入射）或 Bloch（如果源倾斜）边界条件。当使用周期性或布洛赫边界条件时，没有理由模拟设备的多个周期。如果你用两个单元重新运行模拟，你最终会得到完全相同的答案，但运行时间至少是原来的两倍！
1.2 感兴趣的量
可以从超材料器件的 FDTD 模拟中直接测量许多感兴趣的量，包括
? \bullet ?整个结构的场增强
? \bullet ?透射和反射光谱
? \bullet ?散射和吸收截面
? \bullet ?圆二色性
此外，有效的材料特性通常也很重要。参数提取是可能的，但需要对模拟结果进行额外的后处理。有关计算的更多信息，请参见以下页面：
? \bullet ?S parameters
? \bullet ?effective material parameters: refractive index, impedance, permittivity and permeability（折射率、阻抗、介电常数和磁导率）
二、超材料 S 参数提取
这一部分就是第一章S parameters链接的内容，介绍如何计算超材料器件的超材料 S 参数。
2.1 了解 S 参数
S 参数描述 2×2 网络或传输线的行为（见下图）：

文章图片

2 端口网络的 S 参数
对于两个给定的输入信号 a1 和 a2，输出 b1 和 b2 可以计算为

文章图片

2.2 从 FDTD 仿真中测量 S 参数
S 参数是复振幅反射和传输系数（与功率反射和传输系数相反）。例如，对于a1入射，S11 是反射系数，S21 是透射系数；对于a2入射， S22 是反射系数，S12透射系数。

文章图片

源和监视器的相位补偿
认识到 a、b 系数与电场成正比，S 参数可以从上面的定义中计算出来，例如， S 11 = b 1 / a 1 = E r / E i S11=b1/a1=E_r/E_i S11=b1/a1=Er?/Ei?和 S 21 = b 2 / a 1 = E t / E i S21=b2/a1=E_t/E_i S21=b2/a1=Et?/Ei?，其中， E i E_i Ei?， E r E_r Er?， E t E_t Et?是入射、反射和透射电场。这种技术使得获取给定设备的 S 参数变得非常简单，但它确实做出了许多假设。确保所有假设对您的模拟都有效，这一点很重要。
? \bullet ?要测量近场中的 S 参数（如在 S 参数分析对象中所做的那样），透射场和反射场必须像测量点的简单平面波一样传播。如果存在渐逝场（evanescent fields），则测量将不正确。为避免渐逝场，请确保监视器离结构足够远。
? \bullet ?在您的模拟中，源和监视器始终与超材料表面有一定距离（例如，监视器必须离结构足够远以避免渐逝场）。当场从源传播到超材料，从超材料传播到监视器时，这些场将积累额外的相位。 S 参数旨在仅表征超材料，没有这个额外的传播阶段，因此我们必须补偿这个额外的阶段。假设入射空间中的波数为k i k_i ki?，传输空间中的波数为k t k_t kt?，Monitor T 中来自源 S 的额外相位为k i r s + k t r t k_ir_s+k_tr_t ki?rs?+kt?rt?，其中r s r_s rs?和 r t r_t rt?是距离（因此它们是正的）。对于反射监视器 R，来自源 S 的额外相位是k i r s + k i r r k_ir_s+k_ir_r ki?rs?+ki?rr?。
基于上述讨论，对象库中提供了一个现成的光栅 S 参数分析组。在大多数情况下，您只需在 Analysis–Variables 中设置板厚度、板中心位置和源位置沿 x 轴，以及背景折射率。相位补偿在脚本中完成。请注意，由于反射波和透射波必须像平面波一样传播，我们使用点式频域监视器（oint-type frequency-domain monitor）来记录电场分量。分析组还包含两个 2D 表面监视器，用于测量功率传输并检查场实际上是否像单个平面波一样传播。
该分析组根据 DFT 监视器为透射和反射收集的近场数据的光栅投影计算超材料的 S 参数（幅度和相位）。分析组隔离特定的输出光栅阶数，可由用户选择（参见下面的分析属性描述）。定义了 S 参数，以便对于平面结构，它们与根据参考文献 [1] 中的约定定义的菲涅耳系数相匹配（见下图）。此约定与 stackrt 脚本命令的行为一致。因此，如果超材料上方和下方的介质具有不同的折射率， ∣ S 12 ∣ 2 |S_{12}|^2 ∣S12?∣2和 ∣ S 21 ∣ 2 |S_{21}|^2 ∣S21?∣2可以大于 1。分析组可以处理斜入射以及输入和输出的 s 和 p 偏振光。
位置：

文章图片

stackrt函数：

文章图片

文章图片

该分析组有几个设置属性，在对象内部的设置脚本中有详细描述。这里值得一提的是：
? \bullet ?“起始波长”和“终止波长”，指定光源的波长范围。如果起始波长和终止波长不相同，可以在全局监视器属性中设置使用的频点数。
? \bullet ?“source_type”，用于在两种类型的源之间进行选择：1 用于 Bloch/periodic，2 用于 BFAST。仅将“source_type”= 2 用于斜入射的宽带模拟（参见 Source - BFAST）。在本例中，我们将仅考虑法向入射，因此我们设置“source_type”=1。

文章图片

此外，设置以下分析属性很重要（更多信息请参阅“s_params_adv”的分析脚本中的信息）：
（没找到上面这个脚本在哪儿，也许是指分析组里的script吧）
? \bullet ?“超材料中心（metamaterial center）”和“超材料跨度（metamaterial span）”指定您认为是全局坐标中超材料的区域。由于源的位置和计算远场的半球，此信息用于减去 S 参数的附加相位。这类似于此处描述的相位校正过程。
? \bullet ?“target_grating_order_out”指定将考虑计算 S 参数的特定输出光栅顺序。
? \bullet ?“suppress_warnings”禁用 (suppress_warnings = 1) 或启用 (supress_warnings = 0) 带有警告消息的弹出窗口（参见下面的描述）。如果生成警告，“警告”结果将为 1。

文章图片

如果出现以下情况，将生成警告：
? \bullet ?光源的极化既不是 S 也不是 P。
? \bullet ?未找到某些频率点的所需顺序。当不支持目标光栅顺序时会发生这种情况。
分析组生成两个 S 参数结果：“S”和“S_polarization”。第一个是主输入偏振和相同输出偏振的 S 参数（例如，如果源是 s 偏振，它将仅计算 s 偏振透射和反射光的结果）。第二个是主要输入极化和两个输出极化的 S 参数，它仅在我们期望极化被超材料旋转时才有用。在这里考虑的简单情况下，没有偏振旋转，因此我们可以简单地使用“S”结果及其属性：“S11_Gn”和“S21_Gn”（介质 1 中的源）和“S22_Gn”和“S12_Gn”（介质2中的源).
三、有效的超材料体积特性
这一章是第一章effective material parameters链接的内容，介绍如何使用材料的 S 参数计算超材料的有效散装材料属性。
3.1 有效材料属性（来自 S 参数）
我们使用 Smith 等人描述的技术。从 S 参数测量中提取有效的材料参数。特别是，以下技术假设设备在前向和后向传播中表现对称。如果设备不具有这种对称性，则必须对分析进行适当修改，如下所述。
从参考论文[4]方程式 (9) 中，有效折射率 neff 计算为：

文章图片

有效阻抗为

文章图片

一旦获得有效折射率和有效阻抗，就很容易得到有效介电常数和有效磁导率，如下所示：

文章图片

注意：重要的是要记住参数提取并非易事。明确地确定 n 和 z（因此确定 epsilon 和 mu）是该领域的挑战之一，并且仍然是一个活跃的研究领域。一个问题是上述函数是多值的。选择错误的根会导致错误的结果。上面的计算是在 S 参数分析对象中实现的，在某些情况下有效，但肯定不是在所有情况下都有效。有效参数 - Smith 示例中提供的实现应被视为参数提取工作的起点，而不是适用于所有情况的稳健分析。有关替代提取方法，请参阅本页顶部提供的 Szabó 参考（没找到他说的这个在哪儿）。
3.2 非对称（或非均匀）情况下的有效材料属性
本节将简要概述计算有效材料属性的方法，该方法适用于设备在前向和后向源传播方向上表现不同的情况。换言之，当散射参数S11不等于S22时。如果超材料的一侧有衬底，则可能会出现这种情况。在这种非对称情况下，以下等式仍然适用：

文章图片

我们可以使用以下等式来获得 n 和 z：

文章图片

T 矩阵的值也可以从 S 参数计算：

文章图片

Smith 等人也讨论了上述方法。与对称情况一样，由于选择方程的正确根存在挑战，参数提取并非易事。
四、等离子体超材料吸收器
在本例中，我们将根据参考文献 [6] 计算等离子体超材料吸收体的吸收特性。这些器件的吸收机制主要受局域电磁共振的激发支配，并且高度依赖于顶部金属层的几何形状和介电层的厚度。

文章图片

4.1 模拟设置
等离子体吸收器由位于薄介电层和高反射厚金属层顶部的亚波长金属贴片的周期性阵列组成。对于光学模拟，我们只需要模拟单个晶胞。文件plasmonic_absorber.fsp 包含如图 1 所示的结构。设备的设计参数可以在对象树中“model”对象的“model”选项卡中定义。

文章图片

MIM等离子超材料吸收器的晶胞
对称和反对称边界条件都用于减少模拟时间。两个功率监视器，“R”和“T”，用于计算设备的透射和反射（以及吸收）光谱。为了与参考文献 6 中的模拟一致，将一种称为“silver”的Plasma (Drude) 材料添加到材料数据库中，并分配给顶部和底部金属层，并将 1.75 的折射率用于介电层。
4.1.1 对称边界的设置
来自官网的介绍：

文章图片

文章图片

文章图片

也就是说：
当使用对称边界时，该边界仅应用于最小边界（min）,大边界（max）应设置为没有对称性的任何边界，通常是PML。

文章图片

唯一的例外是对周期结构使用对称界面时，在这种情况下，在周期方向上同时设置最小和最大边界以使用对称或反对称边界。需勾选allow symmetry on all boundaries.

文章图片

4.1.2 相关参数设置
这里有关结构的大小、4.1.3的材料模型参数的设置，在参考文献6中均有体现，如下图的截图所示。

文章图片

由于这里仿的是一个周期，所以论文中a=250nm，被用在了FDTD区域的设置参数上。

文章图片

边界条件设置：
由于是周期性结构，这里勾选了allow symmetry on all boundries，将x,y的min和max都设置为了对称或反对称条件。

文章图片

mesh
设置了mesh的大小，使其刚好括住金属铁片，其z方向上向下再延伸一部分到衬底底部。

文章图片

这里dx,dy,dz设置为0.0025，相当于将x,y方向上分成200份，进一步细化。

文章图片

4.1.3 Plasma (Drude)材料

文章图片

以下为示例中添加的称为“silver”的Plasma (Drude)材料

文章图片

4.1.4 光源
光源范围设置为0.5-1.5um

文章图片

4.1.5 监视器
【FDTD学习笔记|Lumerical官方案例、FDTD时域有限差分法仿真学习（六）——等离子体超材料吸收器（Plasmonic metamaterial absorber）】R和T这两个监视器都很大，搞不懂为啥这么设置，反正也只记录仿真区域内的数据。

文章图片

折射率监视器

文章图片

因为用了对称条件，所以这里也仅显示x上0-0.125um的折射率分布。
也是仅记录FDTD仿真区域内的数据，所以z上范围是-2到2，不是geometry里的数据-2.4到2.4

文章图片

时间监视器
由于只是记录场的衰减，所以没必要记录那么多没用的数据。

文章图片

xz_profile监视器
看场强分布的

文章图片

文章图片

profile和power monitor的区别：profile网格化计算，需要插值计算，有一点误差，求算一个严格指定区域的场分布；power位置变化不超过半个网格，以减少插值带来的误差，计算透射率和功率。
但是我记得以往的示例中，一般用profile的监视器，这里可能是为了保证精度吧。

文章图片

4.2 结果
模拟 plasmonic_absorber.fsp 完成后，运行plasmonic_absorber_RTA.lsf 以创建如下所示的图像。

文章图片

R / T / A 光谱（即来自文献[6]的图2）
文献中的图：

文章图片

4.2.1 plasmonic_absorber_RTA.lsf

# Reflection/transmission/absorption plot f = getdata('R','f'); T = -transmission('T'); R = transmission('R'); ?A = 1-R-T; # 计算吸收 plot(c*1e6/f,R,T,A,'wavelength (um)','R/T/A'); legend('R','T','A'); # |H| at resonance H = getresult('xz_profile','H'); H2 = pinch(H.H2); # 获取|H|^2 x = H.x; # 得到x的坐标 z = H.z; lambda = H.lambda; ?fi = find(A,max(A)); # find the frequency point where the absorption becomes maximum ?lambda = lambda(fi); H2 = pinch(H2,3,fi); image(x*1e6,z*1e6,H2,"x (um)","z (um)","|H|^2 at "+num2str(round(lambda*1e9))+" (nm)");

计算吸收公式

FDTD学习笔记|Lumerical官方案例、FDTD时域有限差分法仿真学习（六）——等离子体超材料吸收器（Plasmonic metamaterial absorber）

文章图片

H.H2 表示获取|H|^2

文章图片

H2 = pinch(H2,3,fi); 这一句的解释
H2有三个维度的信息，通过函数size(H2)可发现，其为31×64×200，31是x坐标的个数，64是y的，200是之前设置监视器的200个频率点，即0.5到1.5um之间200个波长的强度信息。这里表示把第三个维度也就是波长删除，但保留需要的那个波长的维度，也就相当于只取了目标波长，在x,y平面记录的强度信息。
这句操作后的size结果为31×64。

文章图片

共振时的磁场强度
运行后，可见其结果与官网相同，如下图。

文章图片

对于这些类型的超材料吸收体，当正确模式被激发时，共振的光谱位置不会随着入射角的变化而变化。这是红外传感应用非常理想的特性，因为检测到的光通常来自许多不同的角度。为了验证所需的吸收光谱对于非法向入射角是否稳健，可以利用宽带固定角源技术 (BFAST) 来模拟一系列入射角的吸收光谱。脚本 plasmonic_absorber_angle_sweep.lsf 使用 BFAST 扫描宽带平面波源的入射角。从下图中可以看出，吸收光谱与照射角度的依赖性很小，正如预期的那样。

文章图片

吸收作为波长和 TE 偏振入射角的函数（即来自文献[6]的图3）
4.2.2 plasmonic_absorber_angle_sweep.lsf

theta = linspace(0,60,7); # 7个角度数据 A = matrix(length(theta),200); for (i = 1:length(theta)){switchtolayout; setnamed('source1','plane wave type','BFAST'); # enable BFAST setnamed('source1','angle theta',theta(i)); run; T = -transmission('T'); R = transmission('R'); A(i,1:200) = 1-R-T; ?i; }# revert to original switchtolayout; setnamed('source1','plane wave type','Bloch/periodic'); # disable BFAST setnamed('source1','angle theta',0); f = getdata('R','f'); image(c*1e6/f,theta,transpose(A),'wavelength (um)','theta','A'); # transpose矩阵转置f = getdata('R','f'); image(theta,c*1e6/f,(A),'wavelength (um)','theta','A');

transpose

文章图片

参考文献
[1]D. R. Smith et al., “Electromagnetic parameter retrieval from inhomogeneous metamaterials”, Phys Rev E 71, 036617 (2005)
[2]Szabó, Z., G.-H. Park, R. Hedge, and E.-P. Li, “A unique extraction of metamaterial parameters based on Kramers-Kronig relationship,” IEEE MTT, vol. 58, no. 10, 2010, pp. 2646-2653
[3]P. Yeh, “Optical Waves in Layered Media”, Wiley-InterScience, chap. 3, 2005.
[4]D. R. Smith et al., “Electromagnetic parameter retrieval from inhomogeneous metamaterials”, Phys Rev E 71, 036617 (2005)
[5]Szabó, Z., G.-H. Park, R. Hedge, and E.-P. Li, “A unique extraction of metamaterial parameters based on kramers-kronig relationship,” IEEE MTT, vol. 58, no. 10, 2010, pp. 2646-2653
[6] J. Hao et al., “Nearly total absorption of light and heat generation by plasmonic metamaterials,” Phys. Rev. B 83, 165107 (2011).