3. Mueller矩阵与Stokes-Mueller formalism
好，咱们进入第三章。前面我们聊了Stokes矢量，它像个快照，能告诉你一束光当前的偏振态。但光学系统是会改变偏振态的——透镜会引入应力双折射，反射镜会改变s和p分量的相位，散射介质会退偏。怎么描述这种「输入→输出」的变换关系？
答案就是Mueller矩阵。我个人觉得，这是偏振测量里最优雅的数学工具之一。它把整个光学系统当成一个黑盒子，你输入一个Stokes矢量，它输出另一个Stokes矢量，中间全靠一个4×4矩阵搞定。
3.1 Mueller矩阵的定义
定义其实很简单。假设入射光的Stokes矢量为 Sin，出射光的Stokes矢量为 Sout，那么：
S_out = M · S_in展开写就是：
┌ S0_out ┐ ┌ m00 m01 m02 m03 ┐ ┌ S0_in ┐
│ S1_out │ = │ m10 m11 m12 m13 │ │ S1_in │
│ S2_out │ │ m20 m21 m22 m23 │ │ S2_in │
└ S3_out ┘ └ m30 m31 m32 m33 ┘ └ S3_in ┘每个元素mij都有物理意义。m00是总透过率，第一行描述系统对光强的响应，第一列描述系统对非偏振光的响应。剩下的3×3子矩阵，才是真正处理偏振信息的部分。
关键点：Mueller矩阵可以处理部分偏振光和退偏效应。这是它比Jones矩阵更强大的地方——Jones矩阵只能处理完全偏振光。
嗯，这里要注意：不是随便一个4×4矩阵都能当Mueller矩阵用。它必须满足物理可实现条件——输出光的偏振度不能大于1，总光强不能为负。我在项目里见过有人拿测量数据硬凑矩阵，结果算出来的输出光偏振度1.3，明显是错的。
3.2 典型光学元件的Mueller矩阵
记住几个常用的，工作中会反复遇到。
3.2.1 理想线偏振器
透光轴沿水平方向（0°）的理想线偏振器：
M_pol(0°) = ½ ┌ 1 1 0 0 ┐
│ 1 1 0 0 │
│ 0 0 0 0 │
└ 0 0 0 0 ┘如果透光轴在θ角度，需要做旋转。旋转矩阵长这样：
M_pol(θ) = R(-θ) · M_pol(0°) · R(θ)其中旋转矩阵R(θ)为：
R(θ) = ┌ 1 0 0 0 ┐
│ 0 cos(2θ) sin(2θ) 0 │
│ 0 -sin(2θ) cos(2θ) 0 │
└ 0 0 0 1 ┘为什么是2θ？因为偏振是二阶张量效应，旋转一次光学元件，偏振态变化两倍角度。我第一次推导时也卡在这儿，后来画了个Poincaré球才想明白。
3.2.2 波片（相位延迟器）
快轴沿水平方向、延迟量为δ的波片：
M_ret(δ, 0°) = ┌ 1 0 0 0 ┐
│ 0 1 0 0 │
│ 0 0 cos(δ) sin(δ) │
│ 0 0 -sin(δ) cos(δ) │
└ 0 0 0 1 ┘常用的两个特例：
- 四分之一波片（QWP）：δ = 90°，cos(δ)=0，sin(δ)=1
- 半波片（HWP）：δ = 180°，cos(δ)=-1，sin(δ)=0
实战技巧：我在搭建Mueller矩阵成像系统时，习惯先用QWP和线偏振器做系统标定。QWP能把线偏振光变成圆偏振光，反过来也能把圆偏振光变成线偏振光——这个特性在测量中特别有用。
3.2.3 反射镜
理想反射镜的Mueller矩阵：
M_mirror = ┌ 1 0 0 0 ┐
│ 0 1 0 0 │
│ 0 0 -1 0 │
└ 0 0 0 -1 ┘注意S2和S3分量反号。这意味着右旋圆偏振光经反射后会变成左旋圆偏振光。如果你在光学系统中用了奇数个反射镜，记得补偿这个符号变化。
3.2.4 退偏器
理想退偏器把所有偏振光变成非偏振光：
M_depol = ┌ 1 0 0 0 ┐
│ 0 0 0 0 │
│ 0 0 0 0 │
└ 0 0 0 0 ┘实际退偏器做不到这么完美，但原理就是这个意思。
3.3 Stokes矢量与Mueller矩阵的运算规则
这部分其实不复杂，但有几个坑要避开。
3.3.1 级联系统
光依次经过多个元件时，总Mueller矩阵是各元件矩阵的乘积，顺序从右往左：
S_out = M_N · ... · M_2 · M_1 · S_in光先经过M1，再经过M2，最后经过MN。写代码时最容易搞反顺序，我吃过这个亏——有一次仿真结果死活不对，查了两天才发现矩阵乘反了。
3.3.2 光强计算
出射光的总光强就是S0分量：
I_out = S0_out = m00·S0_in + m01·S1_in + m02·S2_in + m03·S3_in如果入射光是非偏振光（Sin = [1, 0, 0, 0]T），那么I_out = m00。所以m00也叫「非偏振光透过率」。
3.3.3 偏振度计算
从输出Stokes矢量可以算出输出光的偏振度：
DoP = sqrt(S1² + S2² + S3²) / S0如果DoP > 1，说明你的Mueller矩阵有问题，或者测量数据有噪声。我曾经用有噪声的数据算出一个DoP=1.08，后来发现是探测器饱和了。
避坑指南：我曾经在测量生物组织样本时，发现Mueller矩阵的m00元素异常大。排查了半天，原来是样本架有应力双折射，相当于额外串了一个波片。从那以后，我每次实验前都会先测一遍空载系统的Mueller矩阵做背景扣除。
3.4 知识结构总览
下面这张图把本章的核心逻辑串起来了：
3.5 典型元件Mueller矩阵速查表
| 元件 | 参数 | Mueller矩阵（简化形式） |
|---|---|---|
| 水平线偏振器 | 透光轴0° | ½·diag(1,1,0,0) + ½·off(1,1,0,0) |
| 四分之一波片 | 快轴0°, δ=90° | diag(1,1,0,1) + off(0,0,1,0) 调整符号 |
| 半波片 | 快轴0°, δ=180° | diag(1,1,-1,-1) |
| 理想反射镜 | 垂直入射 | diag(1,1,-1,-1) |
| 理想退偏器 | 完全退偏 | diag(1,0,0,0) |
我的习惯：做实验前，我会先在Python里把整个光路的Mueller矩阵链算一遍。用numpy做矩阵乘法，几分钟就能跑完。如果仿真结果和实验对不上，至少能排除数学错误。代码大概长这样：
import numpy as np
def polarizer(theta):
"""返回透光轴在theta的线偏振器Mueller矩阵"""
R = np.array([[1, 0, 0, 0],
[0, np.cos(2*theta), np.sin(2*theta), 0],
[0, -np.sin(2*theta), np.cos(2*theta), 0],
[0, 0, 0, 1]])
M0 = 0.5 * np.array([[1, 1, 0, 0],
[1, 1, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, 0]])
return R.T @ M0 @ R
def qwp(fast_axis):
"""返回快轴在fast_axis的四分之一波片Mueller矩阵"""
R = np.array([[1, 0, 0, 0],
[0, np.cos(2*fast_axis), np.sin(2*fast_axis), 0],
[0, -np.sin(2*fast_axis), np.cos(2*fast_axis), 0],
[0, 0, 0, 1]])
M0 = np.array([[1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 0, 1],
[0, 0, -1, 0]])
return R.T @ M0 @ R
# 示例：光经过QWP(45°)再经过水平偏振器
M_total = polarizer(0) @ qwp(np.pi/4)
S_in = np.array([1, 0, 0, 1]) # 右旋圆偏振光
S_out = M_total @ S_in
print(f"输出Stokes矢量: {S_out}")跑一下这段代码，你会发现右旋圆偏振光经过45°QWP变成水平线偏振光，再经过水平偏振器全部通过——这就是QWP的线偏振光→圆偏振光转换的逆过程。
好了，Mueller矩阵这块儿先到这儿。记住核心思想：一个4×4矩阵，搞定所有偏振变换。后面讲测量方法时，你会发现所有技术都是在想办法求解这个矩阵的16个元素。