
简介面向磁光材料与光学仿真研究者的Ce:YIG晶体磁光效应数值模拟资源包围绕掺铈钇铁石榴石Ce:YIG的磁光特性重点演示一维磁光晶体透射、反射和法拉第旋转的计算方法可作为物理、光学工程及光通信方向学生与工程师的理论验证工具和入门示例。Ce:YIG具有较高的磁化强度和较低的磁损耗是研究磁光效应的理想材料本脚本正好将其用于一维磁光结构分析。压缩包仅含1个MATLAB脚本.m文件大小约1KB轻量精简无需复杂环境即可直接运行也便于根据研究需要修改材料参数或结构尺寸。已有437人学习下载。脚本以Ce:YIG晶体为对象将磁光介质中的偏振演化、界面透反射及法拉第旋转角统合在一套数值过程中帮助使用者直观理解磁场对光传播路径和偏振状态的影响也为磁光隔离器、磁光调制器等器件的快速评估提供参考。1. Ce_YIG_Ce 在算什么磁光晶体的透射与法拉第旋转先把这个标题拆开Ce_YIG_Ce 指的是 Ce:YIG / YIG / Ce:YIG 这样一组对称三层外加一层层重复成的一维磁光光子晶体。做集成光隔离器或磁光调制器的人看到这串字第一反应通常是“法拉第旋转到底能增强多少透射还剩下多少”。纯 YIG 单晶在 1550 nm 附近比法拉第旋转只有 220 deg/cm 上下要转 45° 得堆到接近 2 mm波导上根本放不下Ce:YIG 能把旋转系数抬高一个量级但吸收也跟着上来。把 Ce:YIG 当作高旋转、高损耗的“磁光镜”把 YIG 当作透明相位层叠成 Ce_YIG_Ce 腔就能借共振把等效相互作用长度拉长几倍同时保持透射率不过分恶化。下面这套做法从材料参数出发用传输矩阵法把透射率、反射率和法拉第旋转谱一起算出来。2. 磁光材料选型Ce:YIG 为什么能顶替纯 YIG2.1 磁光材料靠什么转偏振介电张量非对角项法拉第旋转的根源是磁化介质对左旋、右旋圆偏振光给出不同的折射率。光沿 z 方向传播、外加磁场也沿 z 时磁光材料的介电张量可以写成ε_xx 与 ε_yy 相等代表普通介电响应ε_xy -ε_yx i·g其中 g 是磁光参量正负号取决于磁化方向其余分量为零。这样线偏振光进入介质后会被分解成两个圆偏振本征模各自按折射率 n_ 和 n_- 传播。每走 1 cm 产生的旋转角也就是比法拉第旋转直接等于θ_F (π / λ) · (n_ - n_-)实际写仿真脚本时几乎所有文献给的材料参数都是“折射率 n 比法拉第旋转 θ_F”而不是直接给 g。这时候要做一步反推由 n± sqrt(ε_xx ± ε_xy) 和一阶展开 n_ - n_- ≈ g / n 的关系用 θ_F 反算非对角元 ε_xy。import numpy as np LAM 1550e-9 # 设计波长单位 m def eps_from_theta(n, theta_deg_cm): # 输入折射率 n比法拉第旋转 deg/cm带符号 # 输出介电张量对角元 eps_xx 和非对角元 eps_xy theta_rad_m np.deg2rad(theta_deg_cm) * 100.0 # deg/cm - rad/m dn theta_rad_m * LAM / np.pi # n_plus - n_minus eps_xx n ** 2 eps_xy dn * n # 一阶近似 g ≈ n * Δn return eps_xx, eps_xy # Ce:YIG 常用参数示例符号取负表示与纯 YIG 旋转方向相反 print(eps_from_theta(2.25, -1300))这个一阶近似在 g 远小于 ε_xx 时成立。强磁光材料如 Ce:YIG 在 1064 nm 附近 g 与 ε_xx 的比例能达到百分之几误差仍在可接受范围进入可见光波段或遇到磁光谐振则应该直接用复折射率 n_± 构造每层矩阵避开近似。2.2 Ce 掺杂如何改变 YIG纯 YIG 是亚铁磁石榴石Fe³⁺ 在四面体和八面体位的电荷转移跃迁产生本征磁光响应但跃迁能量离近红外较远所以 1550 nm 处旋转系数只有 200 deg/cm 上下。掺入 Ce³⁺ 后Ce³⁺ 的 4f→5d 电偶极跃迁能量刚好覆盖近红外波段而且自旋轨道耦合很强会显著增强非对角介电张量元。宏观表现是 θ_F 的绝对值跳一个数量级且符号通常取负也就是与纯 YIG 的旋转方向相反。材料参数取不到实物时工程上一般按下表给初值材料波长 (nm)折射率 nθ_F (deg/cm)吸收系数 (cm⁻¹)YIG10642.352805~10YIG15502.222172~5Ce:YIG10642.40-350080~120Ce:YIG15502.25-130020~40注意 Ce:YIG 的掺杂浓度、薄膜应力、退火条件都会移动 n 和 θ_F实际镀膜后通常要用椭偏仪加磁光测量重新定标。仿真阶段按上表运行够用但要保证结构结论的稳健性最后一章会给出一个自检技巧。2.3 Ce_YIG_Ce 腔结构让增强不只靠材料单层 Ce:YIG 也有旋转但薄膜厚度通常在百纳米级单程旋转角不到 1°做不出器件级效果。Ce_YIG_Ce 的做法是把旋转和反射耦合起来两层 Ce:YIG 既是法拉第旋转源又是高折射率差形成的部分反射镜中间的 YIG 间隔层决定腔的纵模位置。共振波长附近光在腔内往返多次每次经过 Ce:YIG 层都累积一次旋转等效长度 L_eff 大于物理厚度。这就是一维磁光光子晶体的基本增强机制。周期数 N 也不是越多越好。N 增大后带隙变硬、带边变陡增强倍数上升但透射峰下降、旁瓣变密稍微偏离设计波长旋转角就会大幅衰减。常见做法是先固定 N3做厚度初步扫描确定工作点后再加层数验证。项目标题中“透射和法拉第旋转”与“透射反射”并列本意就是提醒仿真时 T 和 R 必须同时监控因为增强最大的波长往往位于禁带边缘驻波尾部那个位置反射率已经起来了。3. 用传输矩阵法算一维磁光晶体的透射率与法拉第旋转角3.1 为什么用圆偏振本征模而不是 TE/TM 展开普通多层膜仿真用 TE/TM 展开因为各向同性介质中 TE 和 TM 是独立本征模。磁光材料在法拉第配置下线偏振不再是对角化状态本征模变成左旋、右旋圆偏振。若强行用 TE/TM 递推每经过一层磁光层就要做一次偏振基变换矩阵一行内就会混杂偏振耦合项公式复杂且容易错。圆偏振基的优势在于垂直入射且磁化方向沿光传播方向时每个圆偏振本征模在任意界面和任意层内都是独立传播的。换句话说同一个一维多层结构只要分别用 n_ 和 n_- 跑两遍标量传输矩阵就得到两个复透射系数 t₊ 和 t₋之后再做线性偏振基变换就能提取旋转角。下面代码采用的正是这个思路。3.2 可运行的最小 Python 实现从透射谱到法拉第旋转import numpy as np LAM 1550e-9 K0 2 * np.pi / LAM # ---------- 材料参数 ---------- # YIG: 折射率 2.22, 比法拉第旋转 217 deg/cm, 厚度 320 nm # Ce:YIG: 折射率 2.25, 比法拉第旋转 -1300 deg/cm, 厚度 150 nm def scalar_tmm(n_list, d_list, n01.0, ns1.0): 标量传输矩阵跑一遍圆偏振本征模。 n_list: 各层复折射率按入射到出射排序 d_list: 各层厚度单位 m n0, ns: 入射和出射介质折射率 返回复透射系数 t 和复反射系数 r Z [1.0 / n for n in n_list] # 相对阻抗 M np.eye(2) for j in range(len(n_list)): phi K0 * n_list[j] * d_list[j] P np.array([ [np.cos(phi), 1j * Z[j] * np.sin(phi)], [1j * np.sin(phi) / Z[j], np.cos(phi)] ]) M M P # 由总矩阵求 t 和 r t (2.0 / n0) / ((M[0, 1] M[1, 1] / ns) / n0 M[0, 0] M[1, 0] / ns) r t * (M[0, 0] M[0, 1] / ns) - 1.0 return t, r def faraday_response(layers, n01.0, ns1.0): layers: [(n_plus, n_minus, thickness), ...] n_plus_list [np.sqrt(l[0]) for l in layers] n_minus_list [np.sqrt(l[1]) for l in layers] d_list [l[2] for l in layers] t_p, r_p scalar_tmm(n_plus_list, d_list, n0, ns) t_m, r_m scalar_tmm(n_minus_list, d_list, n0, ns) # 线偏振 x 入射转回线性基得到透射场分量 Ex 0.5 * (t_p t_m) Ey 0.5j * (t_m - t_p) chi Ey / Ex theta_far 0.5 * np.arctan2(2 * np.real(chi), 1 - np.abs(chi) ** 2) T 0.5 * (np.abs(t_p) ** 2 np.abs(t_m) ** 2) R 0.5 * (np.abs(r_p) ** 2 np.abs(r_m) ** 2) return np.rad2deg(theta_far), T, R # ---------- 组装 Ce_YIG_Ce 周期结构 ---------- def build_layers(n_ce, n_yig, d_ce, d_yig, N): layers [] for _ in range(N): eps_xx_ce, eps_xy_ce eps_from_theta(n_ce, -1300) eps_xx_yig, eps_xy_yig eps_from_theta(n_yig, 217) # 按入射方向Ce:YIG - YIG - Ce:YIG layers.append((eps_xx_ce eps_xy_ce, eps_xx_ce - eps_xy_ce, d_ce)) layers.append((eps_xx_yig eps_xy_yig, eps_xx_yig - eps_xy_yig, d_yig)) layers.append((eps_xx_ce eps_xy_ce, eps_xx_ce - eps_xy_ce, d_ce)) return layers layers build_layers(2.25, 2.22, 150e-9, 320e-9, N3) theta, T, R faraday_response(layers) print(ftheta_F {theta:.3f} deg, T {T:.4f}, R {R:.4f})这里每一步的关键点scalar_tmm用的是电场和磁场边界条件推导的 2×2 传输矩阵相位因子phi k0 * n * d直接带进cos/sin白磷损耗只需把折射率写成复值。透射系数公式中n0在分母上入射介质为空气时n01会让透射系数直接可比。线偏振 x 入射时透射场的 x 分量是 t₊ 和 t₋ 的平均y 分量来自两者的差这对应圆偏振叠加回线性偏振的变换。theta_far用复模比计算同时适用于存在磁圆二色性、出射光带有椭圆度的情形如果只取相位差的一半会丢失吸收引起的误差。3.3 参数表怎么从目标波长倒推层厚参数默认值作用调大后的效果d_ce150 nmCe:YIG 单层厚度决定腔镜反射率反射增强透射下降θ_F 峰红移d_yig320 nmYIG 间隔层厚度决定纵模位置透射峰和 θ_F 峰都向长波移动N3周期数增强倍数上升带宽变窄旁瓣变多eps_xy由 θ_F 反推磁光耦合强度θ_F 整体抬升但共振峰位置几乎不动初值确定可以按 λ/(4n) 估Ce:YIG 在 1550 nm 的四分之一波长约 172 nmYIG 约 175 nm但实际扫描时一般从比四分之一波长薄 10~20 nm 起步原因是腔共振所需相位与单层驻波不完全一致。后面章节会给出扫描方法把两个厚度同时放开找最优工作点。4. 透射和反射谱扫描把共振峰挪到目标波长4.1 先看频谱透射率峰和 θ_F 峰的位置不同以上面默认参数仿真会看到透射峰出现在约 1560 nm 附近而 θ_F 的峰值在约 1540 nm两者相差 20 nm 左右。这不是误差而是一维磁光晶体带边态的特性透射峰对应腔模中心光在腔内驻波场分布最强θ_F 峰值偏向禁带边缘因为边缘态里两个圆偏振态在反射叠加时相位差积累更明显但该波长反射率更高、透射更低。所以设计目标波长不能只看 T 峰值也不能只看 θ_F 峰值。实际器件通常在目标波长选一个“综合优值”最大的点优值定义为FOM θ_F * sqrt(T)其中 θ_F 取绝对值。sqrt(T) 的作用是对透射损耗加权因为旋转角再大透射只有 0.1 也没有使用价值。4.2 做参数扫描层厚、周期数与波长一起迭代下面这段扫描代码以 1550 nm 为目标波长把 d_ce 和 d_yig 按 10 nm 步长扫一圈记录 FOM 最大的组合。best_fom 0 best_params {} for d_ce_nm in range(120, 181, 10): for d_yig_nm in range(250, 401, 10): layers build_layers(2.25, 2.22, d_ce_nm * 1e-9, d_yig_nm * 1e-9, N3) theta, T, R faraday_response(layers) fom abs(theta) * np.sqrt(T) if fom best_fom: best_fom fom best_params {d_ce_nm: d_ce_nm, d_yig_nm: d_yig_nm, theta_deg: theta, T: T, R: R} print(best_fom, best_params) # 典型结果: d_ce140, d_yig330 附近theta 约 12~18°T 约 0.7~0.8这个双重循环会算 6×1696 个结构每个结构跑两个圆偏振传输矩阵在普通笔记本上秒级完成。扫完后要做两件事看记录里 R 是否合理如果 R 大于 0.4说明腔镜反射过强后续增加周期数意义不大把 FOM 最大组合的 d_ce 固定再对 N 从 2 到 6 扫一遍观察增强倍数是线性增长还是饱和。周期数扫描常见结果是 N 从 2 到 4 增幅明显N5 之后 θ_F 还会涨但带宽收窄到几十纳米且对加工误差变得极敏感。用于实际器件设计建议 N3 或 4。4.3 反射谱在参数反校中的作用反射谱对厚度误差比透射谱更敏感因为 R 的带边位置由每一层的光学厚度共同决定。镀膜样品实测时通常先拿反射谱对着仿真找带边如果实验带边整体右移说明实际 n·d 比设计值大可能是折射率偏高或膜厚偏厚带边变糊则多半是吸收比预期大或界面粗糙散射。这时候不要去调 ε_xy先调 n 或 d 把带边对齐再对比透射峰的深度最后才动磁光参数。这个顺序能避免“用错的厚度配出对的曲线”。4.4 三个容易翻车的坑符号混乱Ce:YIG 的 θ_F 与 YIG 相反如果两层材料都用正号会互相抵消计算出的 θ_F 突然掉到 1° 以下。检查办法是单独算一层 Ce:YIG 的旋转角看符号与输入是否一致。忽略吸收导致伪共振近红外波段纯 YIG 吸收极小但 Ce:YIG 吸收不可忽略至少要在折射率里加 0.005~0.01 的虚部。否则扫描出的高透射峰在实际中不存在。反射功率与透射功率加和不等于 1 却找不到原因空气到衬底界面折射率不匹配时TR 不等于 1 是正常的剩下的是背向散射或衬底吸收。需要检查scalar_tmm里 n0 和 ns 是否设置正确不要一味怀疑程序。5. 一个收尾技巧用正反双向入射差分标定纯法拉第旋转磁光多层里经常同时存在线性双折射和法拉第旋转前者来自膜层应力或斜入射后者来自磁化强度。实验上区分两者的标准办法是互换入射和出射方向线性双折射是互易的正反方向结果一样法拉第旋转是非互易的反向入射时旋转方向相反。仿真里也可以用同样的差分技巧来验证数值实现是否正确。def verify_nonreciprocity(layers): theta_fw, T_fw, R_fw faraday_response(layers) theta_bw, T_bw, R_bw faraday_response(layers[::-1]) diff 0.5 * (theta_fw - theta_bw) # 非互易部分 even 0.5 * (theta_fw theta_bw) # 互易部分应接近 0 return theta_fw, theta_bw, diff, even layers build_layers(2.25, 2.22, 150e-9, 320e-9, N3) res verify_nonreciprocity(layers) print(ffw{res[0]:.3f}, bw{res[1]:.3f}, diff{res[2]:.3f}, even{res[3]:.3f})如果even的绝对值小于 0.01°说明代码里不含线性双折射所有旋转都来自磁光项如果even明显不为零检查是不是 eps_xy 符号在反向结构时没跟着翻转或者层排序写错。这个检查对 Ce_YIG_Ce 这种对称结构特别合适因为正反结构形状完全相同even理论上必须为零任何数值偏差都提示程序有 bug。另一个实际用途是反推材料参数实测时分别测正反向法拉第旋转用上面公式分离出纯非互易旋转再和仿真谱对比。如果仿真中diff和实测相差超过 10%优先怀疑 Ce:YIG 层的 θ_F 取值因为 YIG 层参数相对稳定。这样把仿真、数值自检和实验标定连起来Ce_YIG_Ce 结构从“算出旋转角”到“镀出的膜能对上”就算闭环了。本文还有配套的精品资源点击获取