磁各向异性介质平面波求解:从张量本构到法拉第旋转
2026/9/18 0:19:21 网站建设 项目流程

简介:《磁各向异性介质中的平面电磁波》是一篇面向电磁理论、通信技术与电子器件开发领域专业读者的理论文献。全文基于Maxwell方程组与对称磁化张量,针对线性、均匀、透明且电各向同性的磁晶体介质,系统推导了平面单色波的完整方程组,详细梳理了电场、磁场、波矢量与能流方向的几何关系,包括H、B、e、s四矢量共面、E垂直于B、能流方向与波矢方向不一致等关键结论,并给出E与B振幅比同相速度之间的定量联系。作者进而引入电磁对偶原则,从电各向异性介质中的已有结果导出磁晶体中的“菲涅耳方程”,用以分析平面波的结构、传播方向及偏振特性,为各向异性介质中波动问题的研究提供了完整推导示范和可扩展的理论方法。资源包仅含1个PDF文件,大小139KB,内容精炼但逻辑完整,适合作为高等电磁学、光学或电动力学课程的补充阅读材料,也可供相关技术研发人员在建模与仿真时参考。该资源已有83人浏览学习,值得对磁化介质波动物理机制有深入需求的读者下载研读。

1. 磁各向异性介质中的平面电磁波:先推翻两个直觉

各向异性介质里的平面波,第一反直觉是:波法线方向 k 和能量传播方向并不重合,坡印廷矢量与等相位面法线会有一个夹角。第二个反直觉是:给定一个传播方向,通常存在两个(而不是一个)本征平面波解,各以不同的相速度传播,偏振互相正交。磁各向异性介质把这两件事都放大——介电张量和磁导率张量同时是满阵时,连"寻常波/非寻常波"这种分类都不再安全。搞天线、做隔离器、写 FDTD 后处理的人,最容易在这里栽跟头。这篇直接把从张量本构到色散方程、再到数值求解和偏振演化的链路走一遍,目标是拿到任意 εr 和 μr 张量时,能立刻算出折射率、偏振态和特征模式。

2. 从张量本构到波法线方程:色散矩阵怎么来的

2.1 时谐约定与张量本构

所有推导从频域麦克斯韦方程组出发。这里必须先把时间因子钉死:本文统一用 exp(-iωt)。约定不同,后面所有张量虚部的符号都要跟着翻,尤其是含磁光效应或铁氧体的非对称 μr 时,i 的符号错一个,法拉第旋转方向就反了。

假定介质无源、无空间色散,本构关系写为:

D = ε0 εr E,B = μ0 μr H

其中 εr 和 μr 都是 3×3 复张量。电各向异性常见于晶体,磁各向异性主要来自铁氧体、磁光薄膜以及部分人工电磁材料。磁各向异性的特点是 μr 不对称,典型铁氧体在 z 向静磁化下的张量为:

μr = [[μr, iκ, 0], [-iκ, μr, 0], [0, 0, μz]]

κ 不为零意味着左旋和右旋圆极化波感受到不同的磁导率,这是后面讲的法拉第旋转的根源。同样重要的是,εr 和 μr 都可能含损耗,即张量为复矩阵。

约定时谐因子ε'' 与损耗关系铁氧体非对角元
物理/光学常用exp(-iωt)ε'' > 0 表示吸收写 +iκ
工程电路常用exp(+jωt)ε'' < 0 表示吸收写 -jκ

提示:和商业软件或论文对比前,先确认对方的时间因子。这是各向异性介质仿真对不上数的最常见原因,不是公式错,是符号约定没对齐。

2.2 平面波假设与广义波动方程

考虑均匀平面波解:

E(r,t) = E0 exp(ik·r - iωt)

对空间梯度做替换 ∇ → ik,旋度变为叉乘。把两个旋度方程写出来:

k × E = ω μ0 μr H
k × H = -ω ε0 εr E

从第一个式子解出 H,代入第二个,消去磁场,得到电场满足的代数方程:

k × (μr⁻¹ (k × E)) + k0² εr E = 0

其中 k0 = ω/c。这个式子和各向同性介质的 k×(k×E) + k0² εE = 0 形式相似,但 μr⁻¹ 夹在两次叉乘中间,导致方程不再能简单地化成标量 k² 的关系。用折射率矢量 n = k/k0 代替 k,整理成:

n × (μr⁻¹ (n × E)) + εr E = 0

这是后面所有数值工作的出发点。

2.3 波法线方程:把叉乘写成矩阵

为了程序化处理,把叉乘运算转成矩阵乘法。令单位波法线方向为 n̂,构造斜对称矩阵:

N = [[0, -n̂z, n̂y], [n̂z, 0, -n̂x], [-n̂y, n̂x, 0]]

这个矩阵满足 N·E = n̂ × E。注意 n 的模长并非 1,我们写成 n = √λ n̂,其中 λ = n² 就是待求的折射率平方。代入广义波动方程后:

λ N μr⁻¹ N E + εr E = 0
(λ Q + εr) E = 0,Q = N μr⁻¹ N

方程有非零解的条件是行列式为零:

det(λ Q + εr) = 0

这就是波法线方程,也叫广义 Fresnel 方程。对各向同性介质,Q 的两个横向本征值都是 -1,方程退化出 λ = ε,重根;对单轴晶体,退化成寻常波 λ = εo 和非常波 λ = εoεe / (εo sin²θ + εe cos²θ),θ 是 k 与光轴夹角。磁各向异性情况下 λ 多项式仍只有两个有限根,第三个根对应无穷大,物理上是纵向静电场解,需要丢弃。

3. 用 Python 求解色散方程:特征值法替代三次多项式求根

3.1 为什么不用行列式展开

det(λ Q + εr) 是 λ 的三次多项式,但 Q 是奇异矩阵(秩最大为 2),所以 λ³ 项系数为零,实际是二次多项式。有人会先展开系数再用 np.roots 求根,但面对 3×3 复数张量时,系数展开很容易在数值上损失精度,尤其当 εr 或 μr 接近奇异时。更稳的做法是把行列式方程改写成广义特征值问题:

det(εr + λ Q) = 0 ⇔ det(-εr - λ Q) = 0

这正是 scipy.linalg.eig(a, b) 的标准形式:求 det(a - λ b) = 0 的根。取 a = -εr,b = Q 即可。

3.2 最小可运行代码

import numpy as np from scipy.linalg import eig def plane_wave_modes(eps_r, mu_r, n_hat): """求解磁各向异性介质中的平面波本征模。 参数: eps_r : 3x3 相对介电张量,复数数组 mu_r : 3x3 相对磁导率张量,复数数组 n_hat : 波法线方向矢量,会被自动归一化 返回: modes : 列表,每个元素是 (λ, 电场偏振向量) λ 为折射率平方,偏振向量已归一化 """ n_hat = np.asarray(n_hat, dtype=float) n_hat = n_hat / np.linalg.norm(n_hat) nx, ny, nz = n_hat N = np.array([ [0.0, -nz, ny], [nz, 0.0, -nx], [-ny, nx, 0.0] ]) Q = N @ np.linalg.inv(mu_r) @ N # 求 det(-eps_r - λ Q) = 0 的广义特征值 lam, V = eig(-eps_r, Q) modes = [] for lam_val, evec in zip(lam, V.T): # 无穷大特征值对应纵向静电解,直接丢弃 if not np.isfinite(lam_val): continue if abs(lam_val) < 1e-12: continue evec = np.asarray(evec).reshape(3) evec = evec / np.linalg.norm(evec) modes.append((lam_val, evec)) return modes

这段代码把整个色散问题压缩成了不到 20 行。核心是 Q 的构造:N μr⁻¹ N 把波法线方向的两次叉乘和磁导率逆矩阵揉在一起,任何方向的斜入射都自动处理,不需要针对特殊方向写分支。广义特征值问题的妙处在于允许 Q 奇异,scipy 会返回 inf 特征值,对应非物理解,过滤掉即可。

调用方式很简单:

eps_r = np.diag([2.0, 2.0, 3.0]) # 单轴晶体,光轴在 z mu_r = np.eye(3) # 非磁性 theta = np.deg2rad(45.0) modes = plane_wave_modes(eps_r, mu_r, [np.sin(theta), 0, np.cos(theta)]) for lam, evec in modes: print(f"n^2 = {lam.real:.6f} {lam.imag:+.6f}j") print(f"E = ({evec[0]:.4f}, {evec[1]:.4f}, {evec[2]:.4f})")

期望输出两组根:一组接近 2.0,另一组接近 2.4。前者是寻常波,偏振垂直于 k 与光轴构成的平面;后者是非常波,偏振在该平面内。如果输出与预期不符,先检查 μr 是否为单位阵,再看 n̂ 方向有没有归一化。

3.3 与解析公式对标:参数怎么调才对

单轴介质为数值代码提供了绝佳的验证基准。非常波解析公式为:

λe(θ) = εo εe / (εo sin²θ + εe cos²θ)

扫几个角度对比:

θ (度)λ_numericalλ_analytic误差
02.0000002.0000000
302.1818182.181818< 1e-12
452.4000002.400000< 1e-12
602.6666672.666667< 1e-12
903.0000003.0000000

θ = 0 时两个模式简并,都看到 εo;θ = 90° 时非常波看到 εe。中间角度验证了插值行为。误差全部来自浮点舍入,说明广义特征值法在这个问题上精度足够。

提示:如果你的 εr 或 μr 元素量级差异超过 1e6,先在代码里做归一化。比如把 εr 除以 max(abs(εr)),对应 λ 结果再乘回去。否则广义特征值求解器可能报收敛警告。

4. 特征波与偏振演化:从双折射到法拉第旋转

4.1 从特征值拿回偏振向量

第 3 章的代码已经返回了电场偏振向量:特征向量 V 的每一列就是对应 λ 的模式。原理上,如果 A 是广义特征问题的解,那么 (λ Q + εr) 的零空间向量就是该模式的电场方向。对无损耗介质,两个模式偏振严格正交;有损耗时仍近似正交,但会出现微小的椭圆度。

拿到 E 之后,磁场 H 也能算出来:

H = (1/(ωμ0)) μr⁻¹ (k × E)

这个式子在做能量计算时必须要用。坡印廷矢量平均值为:

⟨S⟩ = 0.5 Re(E × H*)

它的方向就是能量传播方向。把 ⟨S⟩ 和 k 放在一起看,两者夹角就是能流偏转角,这在各向异性介质中可能达到几十度。做天线罩或透镜设计时,这个偏转角直接决定出射波束指向,不能忽略。

4.2 双折射相位差与偏振片设计

设波沿 z 轴传播,两个本征模式折射率分别为 n1 和 n2,板厚为 L。入射场的两个偏振分量分别获得相位延迟:

Δφ = (n2 - n1) k0 L

这个公式是波片设计的基础。Δφ = π 是半波片,可以把线偏振旋转 2α(α 为偏振方向与快轴的夹角);Δφ = π/2 是四分之一波片,把线偏振变成椭圆偏振。若要频率扫描特性,把 k0 = ω/c 代入,Δφ 随频率线性变化,这就是色散型波片的工作原理。

4.3 磁光效应下的圆偏振分裂

各向异性介质中特别值得单独看的是铁氧体。设波沿磁化方向(z 向)传播,介电张量各向同性,μr 取 2.1 节的形式。数值求解后会发现两个本征模式不再是线偏振,而是左旋和右旋圆偏振,且折射率不同:

n±² = ε(μr ± κ)

两个圆偏振的传播常数差导致线偏振入射波在传播过程中偏振面连续旋转,旋转角为:

θF = 0.5 k0 L (n₋ - n₊)

这就是法拉第旋转。与自然双折射不同,法拉第旋转是非互易的——波反向传播时旋转方向不还原,而是叠加。这个性质被用在隔离器和环行器里。

写几行验证代码,观察旋转角随 κ 的变化:

eps_val = 5.0 mu_val = 0.8 + 0.0j kappa_val = 0.4 mu_r = np.array([ [mu_val, 1j * kappa_val, 0], [-1j * kappa_val, mu_val, 0], [0, 0, mu_val] ]) modes = plane_wave_modes(eps_val * np.eye(3), mu_r, [0, 0, 1]) print("两个模式的折射率平方:") for lam, evec in modes: print(f" n^2 = {lam.real:.4f} {lam.imag:+.4f}j")

注意非对角元填 1jκ 和 -1jκ,符号对应 exp(-iωt) 约定。两个 λ 的差正比于 κ,旋转角可以直接从折射率实部差算出。如果发现两个 λ 相等,说明 κ 被设成了 0,或者 n̂ 方向与磁化方向不平行。

5. 落地校验:符号约定、分支选取与对比技巧

5.1 折射率开方时的分支选取

色散方程解出来的是 λ = n²。实际工程中需要的是折射率 n = √λ,这里有一个必须处理的复变函数分支问题。选错分支的后果是场随距离指数增长,看起来像介质在放大信号,实际是数值假象。

电磁波因子 exp(ik·r) = exp(ik0 n z),无源无增益介质的物理要求是场沿传播方向衰减或不增长。这意味着 n 的虚部必须满足 Im(n) ≥ 0。具体分支判断可以按下面这张表来:

λ 的形态n 的取法物理含义
λ > 0 实数n = +√λ正常传播模
λ < 0 实数n = i√λ
λ 复数,Im λ > 0选 Im n > 0 的那一支有损耗介质
λ 复数,Im λ < 0选 Im n > 0 的那一支增益介质需额外判断

提示:开方后务必验证 Im n ≥ 0。如果程序里出现 Im n < 0 的模式,先怀疑分支选错,再怀疑 εr 或 μr 的虚部符号与时间因子不匹配。

5.2 与商业仿真软件对比的注意点

把这段代码的结果和 CST、HFSS 或 FDTD 结果对比时,有三个常见坑。

第一个坑是时间因子。商业软件内部大多用 exp(jωt),而本文代码用 exp(-iωt)。这导致所有非对称张量的非对角元符号相反,法拉第旋转方向、旋磁效应的旋向全部翻过来。对比之前,把 μr 复共轭或者把非对角元取负号再跑一遍。

第二个坑是折射率的定义域。软件里扫频结果给出的是某个模式传播常数 β 随频率的变化,但 β 对应的是哪个模式,需要根据偏振向量去匹配,而不是只看折射率大小。尤其在模式交叉频率附近,偏振状态迅速变化,按数值大小排序会接错支。

第三个坑是能流方向。FDTD 看的是场在网格里的实际传播,而平面波本征分析给出的是波法线方向 k。各向异性介质中两者不一致,所以直接对比"波束出射角"没有意义,要先从模式偏振和 εr、μr 算出坡印廷方向再对比。

5.3 一套快速的合理性自检流程

拿到新的一组张量参数时,我一般按下面三步检查结果。第一步做各向同性退化测试,把 εr 和 μr 都设成单位阵的倍数,确认两个模式简并且 λ = εμ。第二步扫传播方向,对无损耗介质,λ 应该始终为实数,虚部为零或接近 1e-12 量级,出现明显虚部说明某个张量参数出现非物理的耗散项。第三步做互易性抽查,把传播方向取反,重新求解,得到的两组折射率应当互换,即对称关系成立;只有非互易介质才允许破坏这种对称。

这套流程我每次换新材料参数都会跑一遍。磁各向异性介质的计算本身不难,难的是在符号约定、分支选取和模式配对这些环节上保持清醒。把这几个校验步骤固化成脚本,比任何一次理论推导都更能保证结果可靠。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询