はじめに
メタサーフェイス等の検討を行うときに、RCWAといった電場計算の手法を用いて透過率や遅延位相を計算することがある。この手法は他の代表的な電場計算手法と比較しても計算リソースが乏しくても計算を行うことが可能である。
一方で、RCWAの計算を行おうとすると意外とハードルは高い。有料ソフトであるzamaxやLumericalといった光学計算ソフトウェアや有料プログラミング言語のmatlabにあるreticoloといったモジュールが中心である。これらの計算は最適化がされ、計算速度も速いがやはり無料で計算を行いたいといった希望もある。
Pythonの計算packageにもS4(参考：RCWA(FMM)のフリーライブラリーS4の紹介)が存在するが、依存関係が複雑で少し試したが使い勝手がわかりづらかった。
この背景にはLuaといった言語で書かれたらしくて、githubまでいってソースコードを読みたい自分としてはハードルの高さを否定できなかった。
また、Pythonであることから計算速度に限界があり、今後拡張したいと思っても計算速度が非常に重要であるRCWAとは相性が悪い。
そこで今回は計算速度が速く、無料で使えるプログラミング言語であるJuliaのモジュール：RigorousCoupledWaveAnalysisを使って計算を行っていく。
基本的な動作確認
まずはお決まりの方法でjuliaのmoduleをインストールする
using Pkg
Pkg.add(RigorousCoupledWaveAnalysis)RigorousCoupledWaveAnalysis.jlは近年の研究にも使われており、[John2021]はこのパッケージの開発者であり、この論文内でその実用例が示されている。
計算速度に関してもmatlabコードを提示することは不可能だが、reticoloと同等か少し遅い程度の計算速度となっている。まずはsampleコードにある[Zhijie2018]のfig2を再現する。
このコードはgithubに上がっているが、一部間違いがあったので正しいコードを次に示す。
"""
Zhijie Ma, Yi Li, Yang Li, Yandong Gong, Stefan A. Maier, and Minghui Hong,
"All-dielectric planar chiral metasurface with gradient geometric phase,"
Opt. Express 26, 6067-6078 (2018)
"""
#Figure 2 A
using RigorousCoupledWaveAnalysis
#required materials
ge=InterpolPerm(RigorousCoupledWaveAnalysis.ge_nunley) #Ge from interpolated measured values
ox=ModelPerm(RigorousCoupledWaveAnalysis.sio2_malitson) #SiO2 from dispersion formula
air=ConstantPerm(1.0) #superstrate material is air
#parameters of structure and kgrid
N=6
wls=1300:5:2000.0 #wavelength axis
a=900.0 #cell size
lmid=615/a #length of main arm
larm=410/a #length of side arms
w=205/a #width of arms
use_cuda=true
function zshapegeo(lmid,larm,w) #parametrized z-shaped unit cell
cent=Rectangle(w,lmid) #center arm
top=Shift(Rectangle(larm-w,w),.5larm,lmid/2-w/2) #top arm
bottom=Shift(Rectangle(larm-w,w),-.5larm,-lmid/2+w/2) #bottom arm
return Combination([cent,top,bottom]) #combine them
end
#forward propagation
geo=zshapegeo(lmid,larm,w) #call the method
act=PatternedLayer(500,[air,ge],[geo]) #the active layer is air with Ge structures, 500 nm thick
mdl=RCWAModel([act],air,ox) #build the model, with superstrate and substrate material
Rrf=zeros(length(wls)) #Forward rcp reflectivity
Trf=zeros(length(wls)) #Forward rcp transmissivity
Rlf=zeros(length(wls)) #Forward lcp reflectivity
Tlf=zeros(length(wls)) #Forward lcp transmissivity
@time for i in eachindex(wls) #iterate over all wavelengths
λ=wls[i] #get wavelength from array
grd=rcwagrid(N,N,a,a,1E-5,0,λ,air,use_cuda) #build a reciprocal space grid
ste,stm=rcwasource(grd,1) #define source
Rlf[i],Tlf[i]=etm_reftra(sqrt(.5)*(stm+1im*ste),mdl,grd,λ) #lcp propagation
Rrf[i],Trf[i]=etm_reftra(sqrt(.5)*(1im*stm+ste),mdl,grd,λ) #rcp propagation
end
#now invert it
function zshapegeo(lmid,larm,w) #parametrized z-shaped unit cell
cent=Rectangle(w,lmid) #center arm
top=Shift(Rectangle(larm-w,w),.5larm,-lmid/2+w/2) #top arm
bottom=Shift(Rectangle(larm-w,w),-.5larm,lmid/2-w/2) #bottom arm
return Combination([cent,top,bottom]) #combine them
end
#backward propagation, flip device
geo=zshapegeo(lmid,larm,w) #inverted structure
act=PatternedLayer(500,[air,ge],[geo])
mdl=RCWAModel([act],ox,air) #model is now inverted
Rrb=zeros(length(wls))#Backward rcp reflectivity
Trb=zeros(length(wls))#Backward rcp transmissivity
Rlb=zeros(length(wls))#Backward lcp reflectivity
Tlb=zeros(length(wls))#Backward lcp transmissivity
@time for i in eachindex(wls)
λ=wls[i]
grd=rcwagrid(N,N,a,a,1E-5,0,λ,ox) #<- source codeではairだったがoxが正しい
ste,stm=rcwasource(grd,1)
Rlb[i],Tlb[i]=etm_reftra(sqrt(.5)*(stm+1im*ste),mdl,grd,λ)
Rrb[i],Trb[i]=etm_reftra(sqrt(.5)*(1im*stm+ste),mdl,grd,λ)
end
### plot(筆者追記)
using PyPlot
plt = PyPlot
plt.figure(figsize=(8,6))
plt.plot(wls, Tlf, label="LCP, forward", color="blue")
plt.plot(wls, Trf, label="RCP, forward", color="red")
plt.plot(wls, Tlb, label="LCP, backward", color="blue", linestyle="--")
plt.plot(wls, Trb, label="RCP, backward", color="red", linestyle="--")
plt.xlim(1300, 2000)
plt.ylim(0,1)
plt.xlabel("Wavelength (nm)")
plt.legend()
display(plt.gcf())このコードは最後のrcwagridの部分がairとなっていたが、正しくはoxであった。
その結果が以上の図となる。当たり前ではあるが非常に一致した値となっていることが確認できる。
コードの解説と使用例
[ma2024]の論文の結果を参考に実際に透過率の計算を行ってみる。
この論文によるとメタアトムのパラメータは以下のようにまとめられる。
| 項目 | 数値・仕様 |
|---|---|
| 材料（メタアトム） | アモルファスシリコン（α-Si） |
| 基板材料 | 石英（quartz, n ~ 1.45） |
| 波長 | 1064 nm |
| メタアトム周期 | P = 500 nm |
| メタアトム高さ | H = 800 nm |
| メタアトム半径 | 80 nm ～ 147 nm（2π位相範囲で選定） |
また、図より形状は円柱状であることがわかっている。
これらの数値を使って計算を行ってみると以下のようなコードとなる。
using RigorousCoupledWaveAnalysis
a_Si = InterpolPerm(RigorousCoupledWaveAnalysis.si_schinke) # シリコンの屈折率
quartz = ConstantPerm(1.45^2) # 石英の屈折率
air = ConstantPerm(1.0)
p = 500.0 #メタアトムのピッチ
H = 800.0 #メタアトムの高さ
λ = 1064.0 #波長
radius = 50:1:147 #半径設定
N = 7 #フーリエ次数
T_te = zeros(length(radius))
T_tm = zeros(length(radius))
# グリッドを定義
# Nはフーリエ次数、pはピッチ、λは波長、airは空気の屈折率
grd = rcwagrid(N, N, p, p, 1E-5, 0, λ, air)
#入射波を定義
#TE偏光とTM偏光の光源を生成
ste, stm = rcwasource(grd)
for (i, r) in enumerate(radius)
geo = Circle(2*r/p) #ここが重要、pで規格化する必要あり、また直径を代入
#厚さ、屈折率、パターンを指定してパターン化レイヤーを作成
pattern_layer = PatternedLayer(H, [air, a_Si], [geo])
# モデルを定義
# airは上層、quartzは下層
model = RCWAModel([pattern_layer], air, quartz)
# 反射・透過係数を取得
# 入射光, モデル, グリッド, 波長を代入
R_te, T_te[i] = etm_reftra(ste, model, grd, λ)
R_tm, T_tm[i] = etm_reftra(stm, model, grd, λ)
println("Progress: $i/$(length(radius))")
end
# 平均透過率
T_avg = (T_te + T_tm) / 2
ここでフーリエ次数には注意が必要で、20とかにすると一気に計算が重くなる。
また、メタアトムの屈折率は複素屈折率を入れないとうまくいかなかった。
各変数についてわかりやすいようにコメントも残した。
これを可視化してみる。
using PyPlot
plt = PyPlot
plt.figure()
plt.plot(radius, T_avg, label="Transmittance", color="black", marker="o")
plt.xlabel("Radius (nm)")
plt.ylabel("Transmittance")
plt.xlim(minimum(radius), maximum(radius))
plt.ylim(0, 1.0)
plt.grid(true, linestyle="--")
plt.legend()
display(plt.gcf())
こちらも再現性高く計算できたことが確認できた。
今後について
今後はphaseの計算と斜入射の計算を行っていきたい。
また、意外とmatlabと比較しても計算がとても重いこともなく、早かったので今後とも使っていきたい。
コメント