Propagation Family：从 Maxwell 方程到 RS、ASM、Fresnel 与 Fraunhofer

导论：传播模型不是彼此独立的公式

自由空间衍射传播算法看起来很多，但它们并不是彼此独立的公式。把近似发生的位置标出来以后，整个 Propagation Family 实际上非常清楚：

$$
\boxed{
\text{Maxwell}
\rightarrow
\text{Vector Helmholtz}
\rightarrow
\text{Scalar Helmholtz}
}
$$

在 scalar Helmholtz 这一层，有两种等价的传播表述：

$$
\boxed{
\text{Rayleigh-Sommerfeld}
\quad\Longleftrightarrow\quad
\text{Angular Spectrum}
}
$$

然后再施加近轴近似：

$$
\boxed{
\text{Scalar Helmholtz}
\xrightarrow{\text{paraxial}}
\text{Fresnel}
\xrightarrow{\text{far field}}
\text{Fraunhofer}
}
$$

而 Direct Integration、FFT-DI、Shifted Fresnel 并不是新的物理理论：

$$
\text{RS}
\rightarrow
[
\begin{cases}
\text{Direct Integration}\
\text{FFT-DI}
\end{cases}
]
$$

$$
\text{Fresnel}
\rightarrow
\text{Shifted Fresnel}.
$$

这篇文章重点不在“公式长什么样”，而在每一步究竟丢掉了什么物理信息，以及什么时候这个近似会失效。

一、物理起点：从 Maxwell 到 Scalar Helmholtz

起点：Maxwell 方程

考虑线性、均匀、各向同性、无自由电荷、无自由电流介质。

设时间依赖采用

$$
e^{-i\omega t}.
$$

频域 Maxwell 方程为：

$$
\nabla\times\mathbf E = i\omega\mu\mathbf H,
$$

$$
\nabla\times\mathbf H = -i\omega\epsilon\mathbf E,
$$

同时：

$$
\nabla\cdot(\epsilon\mathbf E)=0,
\qquad
\nabla\cdot(\mu\mathbf H)=0.
$$

如果介质均匀，即

$$
\epsilon=\text{const},
\qquad
\mu=\text{const},
$$

那么：

$$
\nabla\cdot\mathbf E=0.
$$

对第一式取旋度：

$$
\nabla\times(\nabla\times\mathbf E)=i\omega\mu\nabla\times\mathbf H.
$$

代入

$$
\nabla\times\mathbf H=-i\omega\epsilon\mathbf E,
$$

得到：

$$
\nabla\times(\nabla\times\mathbf E)=\omega^2\mu\epsilon\mathbf E.
$$

利用恒等式：

$$
\nabla\times(\nabla\times\mathbf E)=\nabla(\nabla\cdot\mathbf E)-\nabla^2\mathbf E,
$$

以及

$$
\nabla\cdot\mathbf E=0,
$$

得到：

$$
-\nabla^2\mathbf E=\omega^2\mu\epsilon\mathbf E.
$$

因此：

$$
\boxed{\nabla^2\mathbf E+k^2\mathbf E=0}
$$

其中

$$
k=\omega\sqrt{\mu\epsilon}=\frac{2\pi n}{\lambda_0}.
$$

这就是 vector Helmholtz equation。磁场同样满足：

$$
\nabla^2\mathbf H+k^2\mathbf H=0.
$$

到目前为止，还没有进入 Fresnel、Fraunhofer，也没有做 paraxial approximation。

第一个真正重要的近似：Vector → Scalar

因为介质均匀，电场的每个 Cartesian component 都满足：

$$
\nabla^2 E_x+k^2E_x=0,
$$

$$
\nabla^2 E_y+k^2E_y=0,
$$

$$
\nabla^2 E_z+k^2E_z=0.
$$

于是我们把其中一个分量抽象为 $U(x,y,z)$，得到 scalar Helmholtz equation：

$$
\boxed{\nabla^2 U+k^2U=0.}
$$

但是这里要区分两个概念。

在均匀介质内部，“每个场分量满足 Helmholtz equation”本身并不是近似。真正的 scalar approximation 是：

我们不再同时求解完整的 $\mathbf E,\mathbf H$，而只用一个 scalar complex field $U$ 描述传播。

这样做意味着忽略或者弱化：

$$
\text{polarization coupling},
$$

$$
E_z,
$$

以及复杂边界引起的 vector effects。

因此 scalar diffraction 通常适用于：

$$
\text{feature size}\gg\lambda,
$$

并且不是极高 NA、强聚焦、亚波长结构、复杂介质界面的问题。

如果研究

$$
\text{high-NA focusing},
$$

$$
\text{nanophotonics},
$$

$$
\text{metasurfaces},
$$

$$
\text{subwavelength structures},
$$

通常就不能只靠 scalar propagation，而需要 FDTD、FEM、RCWA、vector Debye integral 等 Maxwell solver。

所以 Propagation Family 的第一条边界是：

$$
\boxed{
\text{Maxwell}
\xrightarrow{\text{scalar approximation}}
\text{Scalar Helmholtz}
}
$$

后面的 RS、ASM、Fresnel、Fraunhofer 全都在这条边界之后。

二、Scalar Helmholtz 的两种等价传播表述

Scalar Helmholtz 怎么产生传播公式？

现在的问题变成：已知 $U(x’,y’,0)$，怎样求 $U(x,y,z)$？

Scalar Helmholtz：

$$
(\nabla^2+k^2)U=0
$$

对应的自由空间 Green’s function 是：

$$
\boxed{G(\mathbf r)=\frac{e^{ikr}}{4\pi r}}
$$

因为它满足：

$$
(\nabla^2+k^2)G=-\delta(\mathbf r).
$$

Green’s function 可以理解成：

一个点源产生什么场。

所以只要知道边界面上的场，就可以通过 Green’s theorem，把整个传播问题写成边界积分。这一步最终得到 Rayleigh–Sommerfeld diffraction integral。

Rayleigh–Sommerfeld：空间域中的 scalar propagation

设：

$$
r=\sqrt{(x-x’)^2+(y-y’)^2+z^2}.
$$

Rayleigh–Sommerfeld 第一类衍射积分可以写成：

$$
U(x,y,z)=\iint U(x’,y’,0),h_{\rm RS}(x-x’,y-y’,z),dx’dy’,
$$

其中：

$$
\boxed{
h_{\rm RS}=
\frac{1}{2\pi}
\frac{e^{ikr}}{r}
\frac{z}{r}
\left(\frac{1}{r}-ik\right)
}
$$

也可以理解成 Green’s function 对传播方向法向量的导数。

这里最重要的是：

$$
r=\sqrt{z^2+\Delta x^2+\Delta y^2}
$$

被完整保留下来了。没有做：

$$
r\approx z+\frac{\Delta x^2+\Delta y^2}{2z}.
$$

因此 RS 并不是 paraxial diffraction formula。只要 scalar Helmholtz model 本身成立，RS 能够处理相对大的传播角度。

所以：

$$
\boxed{\text{RS 不要求 Fresnel 的近轴条件}}
$$

但它仍然不是 Maxwell-exact。它的物理边界仍然是：

$$
\boxed{\text{scalar diffraction}}
$$

而不是 vector electromagnetic propagation。

三、Angular Spectrum：频域中的 Scalar Helmholtz 解

Angular Spectrum：从另一个方向直接解 Helmholtz

现在回到：

$$
\nabla^2U+k^2U=0.
$$

不要用 Green’s function，而是在 transverse plane 对 $x,y$ 做 Fourier transform：

$$
U(x,y,z)\leftrightarrow\tilde U(k_x,k_y,z).
$$

因为：

$$
\frac{\partial^2}{\partial x^2}\rightarrow-k_x^2,
$$

$$
\frac{\partial^2}{\partial y^2}\rightarrow-k_y^2,
$$

Helmholtz equation 变成：

$$
\frac{\partial^2\tilde U}{\partial z^2}+(k^2-k_x^2-k_y^2)\tilde U=0.
$$

定义：

$$
\boxed{k_z=\sqrt{k^2-k_x^2-k_y^2}}
$$

向 $+z$ 传播的解为：

$$
\tilde U(k_x,k_y,z)=\tilde U(k_x,k_y,0)e^{ik_zz}.
$$

因此：

$$
\boxed{
U(x,y,z)=
\mathcal F^{-1}
\left[
\mathcal F{U_0}
e^{iz\sqrt{k^2-k_x^2-k_y^2}}
\right].
}
$$

这就是 Angular Spectrum Method。

ASM 还能自然看到 evanescent waves

如果：

$$
k_x^2+k_y^2<k^2,
$$

那么 $k_z\in\mathbb R$，对应 propagating plane wave：

$$
e^{ik_zz}.
$$

如果：

$$
k_x^2+k_y^2>k^2,
$$

则：

$$
k_z=i\alpha.
$$

传播项变成：

$$
e^{-\alpha z}.
$$

这就是 evanescent wave。因此 ASM 比 Fresnel 更清楚地暴露了：

$$
\boxed{
\text{propagating spectrum}
+
\text{evanescent spectrum}
}
$$

的结构。

不过对于普通毫米/厘米尺度自由空间光学，evanescent components 通常传播很短距离就衰减掉。

RS 和 ASM 为什么其实是同一个物理模型？

RS：

$$
U_z=U_0*h_{\rm RS}
$$

是 spatial-domain representation。

ASM：

$$
\tilde U_z=\tilde U_0H_{\rm ASM}
$$

是 spatial-frequency representation。

二者本质上都是 scalar Helmholtz equation + outgoing-wave boundary condition 的解。因此可以理解为：

$$
\boxed{\text{RS}\overset{\mathcal F}{\Longleftrightarrow}\text{ASM}}
$$

一个是 Green-function / impulse-response view，一个是 plane-wave / transfer-function view。

所以在连续理论、正确边界条件和充分数值精度下：

$$
\boxed{U_{\rm RS}=U_{\rm ASM}.}
$$

如果程序结果差很多，通常应该先怀疑 sampling，而不是怀疑两种物理理论冲突。

四、统一传输函数模型：空域卷积与频域乘积

RS 与 ASM 分别给出了 scalar Helmholtz 传播的空域和频域表述。现在可以把两者放入同一个线性平移不变系统框架：传播形式本身是统一的，模型之间真正变化的是传播核或传输函数。

空域中的传播核

若点 $(x’,y’)$ 对输出点 $(x,y)$ 的影响只取决于坐标差，则传播可写成：

$$
\boxed{U_z=U_0*h_z}
$$

即：

$$
U(x,y,z)=\iint U_0(x’,y’)h_z(x-x’,y-y’)dx’dy’.
$$

其中 $h_z$ 是传播系统的 impulse response。对于 Rayleigh–Sommerfeld，传播核就是前面得到的 $h_{\rm RS}$。

傅里叶域中的传输函数

根据 convolution theorem：

$$
\mathcal F{U_0*h_z}=\mathcal F{U_0}\mathcal F{h_z}.
$$

定义传输函数：

$$
\boxed{H_z=\mathcal F{h_z}},
$$

传播便可以写成：

$$
\boxed{\tilde U_z=\tilde U_0H_z},
$$

以及：

$$
\boxed{U_z=\mathcal F^{-1}\left[\mathcal F{U_0}H_z\right]}.
$$

所以同一个传播过程具有两种等价表述：

$$
\boxed{\text{空域卷积}\quad\Longleftrightarrow\quad\text{傅里叶域乘积}}.
$$

RS 的 $h_{\rm RS}$ 与 ASM 的

$$
H_{\rm ASM}=\exp\left[iz\sqrt{k^2-k_x^2-k_y^2}\right]
$$

互为 Fourier transform。RS 和 ASM 因而不是两套冲突的理论，而是同一个 scalar Helmholtz 传播模型的空域与频域表示。

Direct Integration：直接计算空域卷积

将 RS 的连续卷积离散化：

$$
\boxed{
U[m,n]\approx\sum_{p,q}U_0[p,q]h_{\rm RS}(x_m-x_p,y_n-y_q)\Delta x\Delta y
}
$$

就得到 Direct Integration。它没有增加新的物理近似，只是在对 RS integral 做 numerical quadrature。

对于 $N\times N$ 的输入和输出，每个 $N^2$ 输出点都要累加 $N^2$ 个输入点，因此计算量约为：

$$
O(N^4).
$$

DI 很慢，但没有 FFT 自带的 periodicity、circular convolution 和 frequency-grid mapping 问题，因此通常适合作为 numerical reference。

FFT-DI：在频域计算同一个 RS 卷积

对 RS kernel 使用 convolution theorem：

$$
\boxed{
U_z=\mathcal F^{-1}\left[\mathcal F(U_0)\mathcal F(h_{\rm RS})\right]
}
$$

就得到 FFT-DI，计算量降低为：

$$
O(N^2\log N).
$$

因此 DI 和 FFT-DI 不是前后递进的物理近似，而是同一个 Rayleigh–Sommerfeld 模型的两种计算方法：

$$
\boxed{
\text{Rayleigh–Sommerfeld}
\longrightarrow
\begin{cases}
\text{DI：空域直接积分}\
\text{FFT-DI：频域快速卷积}
\end{cases}
}
$$

FFT-DI 引入的是数值假设，不是物理假设

DFT 隐含周期边界：

$$
u[n+N]=u[n].
$$

因此直接计算：

$$
IFFT{FFT(U)FFT(h)}
$$

得到的是 circular convolution，而物理上的 RS 是 linear convolution。为了避免 wrap-around，必须进行 zero padding，典型条件为：

$$
N_{\rm FFT}\ge N_U+N_h-1.
$$

所以 FFT-DI 的误差边界来自 sampling、finite window、zero padding 和 quadrature，而不是新的 diffraction approximation。

五、Fresnel：从空域和频域引入近轴近似

频域推导：近似 ASM 的传输函数

现在从 ASM 开始最清楚。

ASM：

$$
k_z=\sqrt{k^2-k_x^2-k_y^2}.
$$

提取 $k$：

$$
k_z=k\sqrt{1-\frac{k_x^2+k_y^2}{k^2}}.
$$

如果：

$$
\boxed{k_x^2+k_y^2\ll k^2}
$$

利用：

$$
\sqrt{1-\epsilon}\approx1-\frac{\epsilon}{2},
$$

得到：

$$
k_z\approx k-\frac{k_x^2+k_y^2}{2k}.
$$

因此 propagation transfer function $e^{ik_zz}$ 变成：

$$
\boxed{
H_F=e^{ikz}\exp\left[-i\frac{z}{2k}(k_x^2+k_y^2)\right].
}
$$

这就是 Fresnel propagation。

所以 Fresnel approximation 真正表达的是：

$$
\boxed{k_\perp\ll k}
$$

也就是：

$$
\boxed{\sin\theta\ll1.}
$$

这就是 paraxial approximation。

空域推导：近似 RS 的传播距离

RS 中：

$$
r=\sqrt{z^2+\rho^2},
$$

其中：

$$
\rho^2=(x-x’)^2+(y-y’)^2.
$$

写成：

$$
r=z\sqrt{1+\frac{\rho^2}{z^2}}.
$$

展开：

$$
r=z+\frac{\rho^2}{2z}-\frac{\rho^4}{8z^3}+\cdots
$$

Fresnel 只保留：

$$
\boxed{r\approx z+\frac{\rho^2}{2z}.}
$$

但是传播相位是 $kr$。所以真正判断误差大小的不是单纯 $\rho\ll z$，而是被丢掉的相位：

$$
\boxed{\Delta\phi\approx\frac{k\rho^4}{8z^3}.}
$$

需要：

$$
\boxed{\frac{k\rho_{\max}^4}{8z^3}\ll1.}
$$

这个条件比一句“$z$ 足够远”更准确。

Fresnel 的空域卷积形式

经过这个近似：

$$
r\approx z+\frac{x^2+y^2}{2z}.
$$

同时 amplitude slowly varying：

$$
\frac{1}{r}\approx\frac{1}{z},
$$

以及：

$$
\frac{z}{r}\approx1.
$$

于是 RS kernel 简化成：

$$
\boxed{
h_F(x,y,z)=
\frac{e^{ikz}}{i\lambda z}
\exp\left[i\frac{k}{2z}(x^2+y^2)\right].
}
$$

因此：

$$
\boxed{U_z=U_0*h_F.}
$$

这是 Fresnel diffraction 的 convolution form。

Fresnel 并不等于“近场”

这也是一个很常见的误解。

Fresnel 的真正条件是：

$$
\boxed{\text{paraxial}}
$$

而不是简单 $z$ 小。

一个系统即使传播距离很长，如果角度仍然很小，Fresnel 仍然可以非常准确。相反，如果存在大角度、高空间频率，即使距离不算特别短，$k_\perp/k$ 比较大，Fresnel 也可能失效。

所以更合理的理解是：

$$
\text{ASM/RS}\xrightarrow{\text{small angle}}\text{Fresnel}.
$$

六、Fraunhofer：从空域和频域走向远场

空域推导：忽略输入面的二次相位

Fresnel integral 可以写成：

$$
U(x,y,z)=
\frac{e^{ikz}}{i\lambda z}
e^{i\frac{k}{2z}(x^2+y^2)}
$$

$$
\times
\iint U_0(x’,y’)
e^{i\frac{k}{2z}(x’^2+y’^2)}
e^{-i\frac{2\pi}{\lambda z}(xx’+yy’)}
dx’dy’.
$$

Fraunhofer approximation 再进一步认为，在整个输入 aperture 上：

$$
e^{i\frac{k}{2z}(x’^2+y’^2)}\approx1.
$$

假设 aperture 最大半径为 $a$，则最大附加相位：

$$
\Delta\phi_{\rm aperture}\approx\frac{ka^2}{2z}.
$$

要求：

$$
\boxed{\frac{ka^2}{2z}\ll1.}
$$

代入：

$$
k=\frac{2\pi}{\lambda},
$$

得到：

$$
\boxed{z\gg\frac{\pi a^2}{\lambda}.}
$$

不同教材根据允许的 phase error，会给出略有不同的常数，因此工程中经常简写成：

$$
\boxed{z\gg\frac{a^2}{\lambda}.}
$$

频域结果：远场是输入场的 Fourier transform

去掉输入 quadratic phase：

$$
U(x,y,z)\approx
\frac{e^{ikz}}{i\lambda z}
e^{i\frac{k}{2z}(x^2+y^2)}
$$

$$
\times
\iint U_0(x’,y’)e^{-i2\pi(f_xx’+f_yy’)}dx’dy’,
$$

其中：

$$
\boxed{f_x=\frac{x}{\lambda z},\qquad f_y=\frac{y}{\lambda z}.}
$$

因此：

$$
\boxed{
U(x,y,z)\propto
\mathcal F{U_0}
\left(\frac{x}{\lambda z},\frac{y}{\lambda z}\right).
}
$$

所以 Fraunhofer pattern 就是 aperture 的 Fourier transform。

Lens：在有限距离实现 Fourier transform

Fraunhofer condition：

$$
z\gg\frac{a^2}{\lambda}
$$

针对的是：

$$
\boxed{\text{free-space propagation}}.
$$

如果加入 thin lens：

$$
t_L(x,y)=\exp\left[-i\frac{k}{2f}(x^2+y^2)\right],
$$

lens 的 quadratic phase 可以和传播中的 quadratic phase 抵消。结果是在后焦面 $z=f$ 直接得到输入场的 Fourier transform。

因此实验室里不需要真的传播几十米甚至几百米才能观察 far-field diffraction。Lens 实际上实现了：

$$
\boxed{\text{finite-distance Fourier transform}.}
$$

七、总结：模型谱系、适用边界与选择方法

最终 Propagation Family

现在可以把整个关系重新画成：

$$
\boxed{\text{Maxwell Equations}}
$$

$$
\downarrow
$$

均匀、线性、各向同性、无源、单色场

$$
\boxed{\text{Vector Helmholtz}}
$$

$$
\downarrow\quad\text{scalar approximation}
$$

$$
\boxed{\text{Scalar Helmholtz}}
$$

从这里分成两个等价表示：

$$
[
\begin{array}{ccc}
\text{Spatial domain} && \text{Spatial-frequency domain}\[4pt]
\downarrow && \downarrow\
\text{Rayleigh-Sommerfeld} && \text{Angular Spectrum}\
\downarrow &&\
\text{DI / FFT\text{-}DI} &&
\end{array}
]
$$

然后施加：

$$
k_\perp\ll k
$$

得到：

$$
\boxed{\text{Fresnel}}
$$

数值实现可以有：

$$
\text{one-step},\quad\text{two-step},\quad\text{multi-step},\quad\text{shifted Fresnel}.
$$

继续增加 far-field 条件：

$$
\frac{ka^2}{2z}\ll1
$$

得到：

$$
\boxed{\text{Fraunhofer}}.
$$

各模型的“前提—丢失信息—失效边界”

最重要的一条判断链

以后遇到任何 propagation 问题，不应该第一反应问：

“用 ASM 还是 Fresnel？”

更合理的是先判断你在哪一层。

如果 polarization、high NA、subwavelength structure 已经重要：

$$
\boxed{\text{Scalar model 本身就不成立}}
$$

应该回 Maxwell solver。

如果 scalar optics 成立，但存在明显大角度传播：

$$
\boxed{\text{RS / ASM}}
$$

更适合。

如果：

$$
k_\perp\ll k,
$$

也就是传播主要沿 optical axis：

$$
\boxed{\text{Fresnel}}
$$

可以成立。

如果进一步：

$$
\frac{ka^2}{2z}\ll1,
$$

则：

$$
\boxed{\text{Fraunhofer}.}
$$

而 DI、FFT-DI、Shifted Fresnel、one-step Fresnel 这些问题，已经不是新的物理层级，而是在确定物理模型之后：

$$
\boxed{\text{How do we compute it numerically?}}
$$

这也是 waveprop 后面真正开始变得有价值的地方：从这里开始，问题不再是“公式会不会”，而是如何保证离散计算真的逼近上面的连续物理模型。

许可协议

本文采用 署名-非商业性使用-相同方式共享 4.0 国际 许可协议，转载请注明出处。

分享文章