简介：面向磁光材料与光学仿真研究者的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 为什么能顶替纯 YIG
2.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⁻¹) |
|---|---|---|---|---|
| YIG | 1064 | 2.35 | +280 | 5~10 |
| YIG | 1550 | 2.22 | +217 | 2~5 |
| Ce:YIG | 1064 | 2.40 | -3500 | 80~120 |
| Ce:YIG | 1550 | 2.25 | -1300 | 20~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 增大后带隙变硬、带边变陡，增强倍数上升，但透射峰下降、旁瓣变密，稍微偏离设计波长旋转角就会大幅衰减。常见做法是先固定 N=3，做厚度初步扫描，确定工作点后再加层数验证。项目标题中“透射和法拉第旋转”与“透射反射”并列，本意就是提醒仿真时 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, n0=1.0, ns=1.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, n0=1.0, ns=1.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, N=3) theta, T, R = faraday_response(layers) print(f"theta_F = {theta:.3f} deg, T = {T:.4f}, R = {R:.4f}")这里每一步的关键点：
scalar_tmm用的是电场和磁场边界条件推导的 2×2 传输矩阵，相位因子phi = k0 * n * d直接带进cos/sin，白磷损耗只需把折射率写成复值。- 透射系数公式中
n0在分母上，入射介质为空气时n0=1会让透射系数直接可比。 - 线偏振 x 入射时，透射场的 x 分量是 t₊ 和 t₋ 的平均，y 分量来自两者的差，这对应圆偏振叠加回线性偏振的变换。
theta_far用复模比计算，同时适用于存在磁圆二色性、出射光带有椭圆度的情形；如果只取相位差的一半，会丢失吸收引起的误差。
3.3 参数表：怎么从目标波长倒推层厚
| 参数 | 默认值 | 作用 | 调大后的效果 |
|---|---|---|---|
| d_ce | 150 nm | Ce:YIG 单层厚度，决定腔镜反射率 | 反射增强，透射下降，θ_F 峰红移 |
| d_yig | 320 nm | YIG 间隔层厚度，决定纵模位置 | 透射峰和 θ_F 峰都向长波移动 |
| N | 3 | 周期数 | 增强倍数上升，带宽变窄，旁瓣变多 |
| eps_xy | 由 θ_F 反推 | 磁光耦合强度 | θ_F 整体抬升，但共振峰位置几乎不动 |
初值确定可以按 λ/(4n) 估：Ce:YIG 在 1550 nm 的四分之一波长约 172 nm，YIG 约 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, N=3) 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_ce=140, d_yig=330 附近，theta 约 12~18°，T 约 0.7~0.8这个双重循环会算 6×16=96 个结构，每个结构跑两个圆偏振传输矩阵，在普通笔记本上秒级完成。扫完后要做两件事：
- 看记录里 R 是否合理，如果 R 大于 0.4，说明腔镜反射过强，后续增加周期数意义不大；
- 把 FOM 最大组合的 d_ce 固定，再对 N 从 2 到 6 扫一遍，观察增强倍数是线性增长还是饱和。
周期数扫描常见结果是 N 从 2 到 4 增幅明显，N=5 之后 θ_F 还会涨，但带宽收窄到几十纳米，且对加工误差变得极敏感。用于实际器件设计建议 N=3 或 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 却找不到原因：空气到衬底界面折射率不匹配时，T+R 不等于 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, N=3) res = verify_nonreciprocity(layers) print(f"fw={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 结构从“算出旋转角”到“镀出的膜能对上”就算闭环了。
本文还有配套的精品资源，点击获取
