Ce_YIG一维磁光晶体透射与法拉第旋转计算：从介电张量到4×4传输矩阵-平芜编程栈

简介：这份资源面向光学、磁光材料与物理仿真方向的学习者和研究人员，围绕掺铈钇铁石榴石（Ce:YIG）这一典型一维磁光晶体，聚焦其透射、反射与法拉第旋转效应的数值模拟。压缩包内共1个文件，为MATLAB脚本（.m格式），整体约1KB，体积轻量，便于直接运行与二次修改。脚本可用于计算不同磁场、波长与晶体厚度条件下光的偏振旋转角及透射、反射系数，帮助理解磁光隔离器、磁光调制器等器件的物理机理。目前已有438人学习下载，说明该方向具备一定关注度。对于需要快速搭建Ce:YIG磁光效应仿真框架、验证法拉第旋转理论公式或开展课程设计、科研预研的读者，这份代码可作为可复用的计算起点，节省从零编写数值模型的时间。

1. Ce_YIG 一维磁光晶体：透射与法拉第旋转到底在算什么

一块厚度几百微米的 YIG 晶体，外加一个可调磁场，就能让穿过它的线偏振光偏振面转几十度，同时透射率还能保持在 80% 以上——这件事在光通信和激光系统里被反复利用，但真正动手算的时候，很多人卡在同一个地方：透射谱和法拉第旋转角到底怎么从材料参数推出来，边界条件怎么设，磁场方向怎么定。Ce_YIG 是在 YIG 里掺铈，目的是把法拉第旋转角在 1550 nm 附近拉高一个量级，代价是吸收边红移、损耗上升。一维磁光晶体的透射和法拉第旋转计算，本质上是解一层或多层各向异性介质中的传播问题，核心变量是介电张量的非对角元。适合做磁光隔离器、环形器、磁场传感的从业者，也适合想从材料参数直接算器件性能、不想只靠仿真软件黑匣子出结果的人。下面从物理量定义一路写到可复现的计算流程和参数边界。

2. 从介电张量到透射矩阵：一维磁光传播的物理骨架

2.1 磁光效应的微观来源与介电张量形式

YIG 的磁光效应来自磁化后电子自旋轨道耦合导致的介电张量非对角元。在饱和磁化、磁化方向沿 z 轴时，介电张量写成：

ε = ε0 * [ ε1 -i*ε2 0 i*ε2 ε1 0 0 0 ε3 ]

其中 ε1 是普通介电常数，ε2 是磁光耦合项，正比于磁化强度 M，ε3 在立方晶系中通常等于 ε1。法拉第旋转角 θ_F 在弱吸收近似下正比于 ε2，而 ε2 又正比于 M，所以外加磁场通过改变 M 来调旋转角。Ce 掺杂的作用是增强自旋轨道耦合，把 ε2 在 1550 nm 附近抬高，代价是 ε1 的虚部（吸收）也增大。

这里有一个容易翻车的地方：很多教材直接给 θ_F = V·B·L，V 是费尔德常数，但这个公式只在远离吸收带、且磁化饱和时成立。Ce_YIG 在 1550 nm 附近吸收不可忽略，必须回到介电张量解 Maxwell 方程，否则算出来的旋转角会偏大 20% 以上。

2.2 一维传播的 4×4 传输矩阵怎么建

一维意味着只考虑沿 z 方向传播、层状结构。对每一层，把电场写成四个分量：Ex、Ey 及其对应的磁场分量。代入 Maxwell 方程后得到本征值问题，解出四个本征模式（两个前向、两个后向），每个模式有各自的传播常数和偏振态。层内传播用对角矩阵表示，层间界面用边界条件匹配。

具体步骤：

对每层材料，由介电张量构造 4×4 矩阵 Δ，解其特征值得到传播常数 k1~k4 和特征向量。
构造层内传播矩阵 P = diag(exp(-ik1d), exp(-ik2d), exp(-ik3d), exp(-ik4d))，d 是层厚。
构造界面矩阵 D，把电场和磁场切向分量从一层映射到下一层。
总传输矩阵 M = D1^-1 * D2 * P2 * D2^-1 * D3 * ... 按层序乘起来。
由 M 和入射/出射半空间的边界条件解出反射系数 r 和透射系数 t。

透射率 T = |t|^2 * (n_out/n_in)，法拉第旋转角 θ_F = 0.5 * atan(2Re(t_xyt_yy* + t_xxt_yx) / (|t_xx|^2 + |t_yy|^2 - |t_xy|^2 - |t_yx|^2))，椭圆率由虚部给出。

注意：这里 t 是 2×2 琼斯矩阵，不是标量。只取 t_xx 算透射率、忽略交叉项，是新手最常见的错误，会导致旋转角算出来恒为零。

2.3 材料参数从哪来：Ce_YIG 的 ε1 和 ε2 取值

Ce_YIG 没有像 Si 或 SiO2 那样通用的数据库。常见做法是：ε1 的实部由折射率 n ≈ 2.2~2.4（1550 nm）反推，虚部由吸收系数 α 换算，ε2 由法拉第旋转角实验值反推。如果手头没有自己的椭偏或磁光克尔谱数据，可以用文献里 Ce_YIG 的典型值起步：ε1 ≈ 5.3 + i*0.02，ε2 ≈ 0.01~0.03，具体随 Ce 浓度和退火条件变化很大。

我一般会先固定 ε1 实部，扫 ε2 从 0.005 到 0.05，看透射率和旋转角怎么变，再拿实验值去卡。这样比一上来就拟合所有参数快得多，也能看出哪个参数是主导。

3. 用 Python 跑通 Ce_YIG 透射谱与法拉第旋转角的最小实现

3.1 环境准备与依赖

只需要 numpy 和 matplotlib，不需要 COMSOL 或 Lumerical。Python 3.9 以上即可。

pip install numpy matplotlib

不依赖任何磁光专用库，因为 4×4 传输矩阵本身只有几十行，自己写反而可控。用现成软件包的问题是参数含义不透明，出了异常值不知道是物理还是设置问题。

3.2 构造介电张量与 4×4 传输矩阵

import numpy as np def eps_tensor(eps1, eps2): """构造磁化沿z轴的介电张量""" eps = np.array([ [eps1, -1j*eps2, 0], [1j*eps2, eps1, 0], [0, 0, eps1] ], dtype=complex) return eps def layer_matrices(eps, k0, d): """返回单层的传播矩阵P和界面矩阵D""" # 构造Delta矩阵，解本征值 # 这里用简化形式：对正入射、磁化沿z，本征模式为左右圆偏振 n_r = np.sqrt(eps[0,0] + eps[0,1]*1j) # 右旋 n_l = np.sqrt(eps[0,0] - eps[0,1]*1j) # 左旋 k_r = k0 * n_r k_l = k0 * n_l P = np.diag([np.exp(-1j*k_r*d), np.exp(-1j*k_l*d), np.exp(1j*k_r*d), np.exp(1j*k_l*d)]) # 界面矩阵D把圆偏振基映射到线偏振基 D = np.array([ [1, 1, 0, 0], [1j, -1j, 0, 0], [0, 0, 1, 1], [0, 0, 1j, -1j] ], dtype=complex) / np.sqrt(2) return P, D

逻辑说明：对正入射、磁化沿传播方向的情况，本征模式是左右圆偏振，折射率分别为 sqrt(ε1 ± ε2)。这个简化在偏离正入射超过 10 度后误差迅速增大，那时必须回到完整 4×4 本征值求解。参数 d 是层厚，单位与 k0 一致，k0 = 2π/λ。

3.3 计算透射率和法拉第旋转角的主循环

def calculate(lam, eps1, eps2, d): k0 = 2*np.pi / lam eps = eps_tensor(eps1, eps2) P, D = layer_matrices(eps, k0, d) # 单层：M = D^-1 * P * D M = np.linalg.inv(D) @ P @ D # 半空间边界：入射n0=1，出射n0=1 # 提取琼斯矩阵t（简化：取M的左上2x2块） t = M[:2, :2] T = np.abs(t[0,0])**2 + np.abs(t[1,0])**2 # 法拉第旋转角 theta = 0.5 * np.angle(t[0,0] + 1j*t[1,0]) - 0.5 * np.angle(t[0,0] - 1j*t[1,0]) return T, np.degrees(theta) # 扫波长 lams = np.linspace(1500, 1600, 200) T_list, th_list = [], [] for lam in lams: T, th = calculate(lam*1e-9, 5.3+0.02j, 0.02, 500e-6) T_list.append(T) th_list.append(th)

逻辑说明：t 矩阵的左上 2×2 块对应线偏振基下的琼斯矩阵。透射率取第一列（x 偏振入射）的模方和。旋转角用左右旋相位差的一半来算，这是正入射下的标准做法。参数 500e-6 是 500 微米厚，Ce_YIG 典型器件厚度在 200~1000 微米之间。

3.4 结果解读与参数敏感性

跑出来的透射谱在 1550 nm 附近如果出现明显干涉条纹，说明厚度和折射率匹配，条纹间距 Δλ ≈ λ²/(2nd)。旋转角随 ε2 线性增大，但透射率随 ε2 增大而下降，因为吸收项被放大。实际设计要在旋转角和插入损耗之间取折中，Ce_YIG 的 ε2 通常选在 0.015~0.025 之间。

提示：如果旋转角算出来是负的，检查 ε2 的符号和磁场方向定义是否一致。符号约定不统一是磁光计算里最常见的玄学问题。

4. 透射反射联算时最容易翻车的五个地方

4.1 现象：透射率大于 1

原因：出射半空间折射率没有归一化，或者 t 矩阵取了错误的块。解决：确认 T = |t|² * Re(n_out)/Re(n_in)，正入射下 n_in = n_out = 1 时退化为 |t|²。

4.2 现象：旋转角随厚度线性增大但实验不线性

原因：厚度超过吸收长度后，多次反射和吸收导致有效旋转饱和。解决：在传输矩阵里保留后向模式，不要只用前向近似。吸收长度 1/α 在 Ce_YIG 里约 1~5 mm，厚度接近这个量级时必须算全矩阵。

4.3 现象：反射率算出来和实验差一个数量级

原因：界面矩阵 D 用了错误的本征模式排序，导致前后向模式混淆。解决：检查 D 的列顺序是否和 P 的对角元顺序一致，前向两个、后向两个，不能交叉。

4.4 现象：改变磁场方向后旋转角不变

原因：介电张量里 ε2 的符号没有随磁化方向翻转。解决：磁化反向时 ε2 → -ε2，法拉第旋转角随之反号，这是法拉第效应的非互易性来源，代码里要显式处理。

4.5 现象：波长扫描出现非物理尖峰

原因：k0*d 在某些波长下使矩阵接近奇异，数值求逆不稳定。解决：改用解线性方程组而不是显式求逆，或者把波长步长减小到 0.1 nm 以下。

5. 进阶：用群论约化参数并做实验对标

5.1 用对称性减少独立参数

Ce_YIG 是立方晶系，磁化沿 [111] 和沿 [100] 时介电张量的非对角元结构不同。沿 [111] 磁化时，ε2 的有效值要乘一个方向因子。如果做的是磁场角度依赖实验，这一步不能省，否则拟合出来的 ε2 会随角度漂移，看起来像材料不均匀，其实是坐标没转对。

def rotate_eps(eps, theta, phi): """把介电张量从磁化坐标系转到实验室坐标系""" # theta, phi 是磁化方向球坐标 # 构造旋转矩阵R，返回 R @ eps @ R.T ...

参数 theta 是磁化与 z 轴夹角，phi 是方位角。对 [111] 方向，theta ≈ 54.7 度，phi = 45 度。旋转后非对角元不再只有 ε2 一个独立量，会出现 ε_xy 和 ε_xz 同时非零。

5.2 实验对标：椭偏仪和磁光克尔谱怎么对

椭偏仪给的是 ψ 和 Δ，对应反射系数比 r_p/r_s。磁光克尔谱给的是克尔旋转角和椭圆率。把计算出的反射矩阵转成 ψ、Δ 和克尔角，和实验曲线叠在一起看。如果透射谱对得上但克尔谱对不上，问题多半在界面层——Ce_YIG 表面容易形成非磁性的死层，厚度几纳米到几十纳米，对透射影响小但对反射影响大。

我一般会在模型里加一层 5~20 nm 的界面层，ε2 = 0，ε1 取体材料值，然后看克尔谱能不能压下去。这个死层参数没有通用值，必须用自己的样品去卡。