第一次在Matlab里看到涡旋光束的相位图时我盯着屏幕愣了好一会儿——整个相位面像被拧成了螺旋从中心一点向外一圈圈转出去。那个中心点相位是个奇点光强正好是零所以光斑中间是个标准的小圆洞。也就是这一眼让我把课本上“携带轨道角动量的光”这句抽象的描述彻底看明白了。这篇内容主要围绕“Matlab仿真常见涡旋光束”展开我会从涡旋光束的物理图像讲起把拉盖尔-高斯光束的复振幅构造、干涉图样验证、自由空间传播仿真、高阶涡旋与叠加态等内容逐步串起来最后附上我反复踩过的坑。无论你是光学方向的本科生、研究生还是刚开始接触光场仿真的工程师只要能用Matlab写几行矩阵运算就可以照着这篇文章造出属于自己的涡旋光束。1. 螺旋相位与轨道角动量先弄懂涡旋光束在“旋”什么1.1 指数因子 exp(ilφ) 的物理直觉涡旋光束最核心的东西其实就是一个相位因子E ∝ exp(ilφ)这里φ是光场横截面上的方位角从0到2π绕一圈l是一个整数叫拓扑荷数。这个因子意味着当你在横截面上沿着某个半径绕中心走一整圈时相位并不是回到原值而是变化了2πl。相位变化了整数个2π所以电场本身仍然连续但中心那一点的相位无法定义是一个奇点。为了让波动方程成立中心振幅必须强迫为零这就是为什么涡旋光束的光强图中间总有一个暗洞。用生活化的方式理解普通高斯光束的等相位面像个平整的球冠光往前传播时波前是光滑的涡旋光束的等相位面则像被“拧”过的弹簧面拧了几圈拓扑荷数就是几。拓扑荷的正负代表旋转方向好比左旋螺纹和右旋螺纹。我第一次仿真l1的涡旋光束时觉得这个“拧”看不太出来直到画出相位图才意识到整个相位分布就是一把“旋转楼梯”。相位从-π渐变到π中间夹着一条相位跳变线。很多资料把这个跳变线画成从中心向外的一条“切缝”实际上它对涡旋结构至关重要。1.2 涡旋光束家族为什么拉盖尔-高斯束最常用涡旋光束不是一个单一的光束而是一类带有螺旋相位的光束的总称。常见的有拉盖尔-高斯光束Laguerre-GaussianLG柱坐标系下傍轴波动方程的本征解携带确定轨道角动量实验上最容易用螺旋相位板或空间光调制器产生。贝塞尔-高斯光束Bessel-Gaussian具有无衍射特性中心亮斑或暗斑可以在较长距离内保持不扩散适合光镊和成像。艾里涡旋光束、马丢涡旋光束等特殊场景下使用。但如果只做基础仿真我强烈建议从LG束入手原因有三个。第一它的复振幅表达式是解析的只需几条Matlab语句就能构造第二它在自由空间传播时保持涡旋结构很适合用角谱法观察演化第三它直接对应轨道角动量(OAM)的本征态后续不管做OAM复用、涡旋光通信还是超表面设计都是绕不开的基础。需要区分的是圆偏振光对应的自旋角动量SAM和涡旋光束对应的轨道角动量OAM是两种不同的自由度。自旋角动量来自电场矢量的旋转轨道角动量来自波前相位结构的旋转。一个线偏振涡旋光束也能携带OAM只是没有SAM。这个区别在仿真中不会直接影响光强和相位但如果做的是矢量光束或紧聚焦仿真就要特别注意。2. 仿真起点二维网格怎么搭才不出“鬼影”2.1 网格点数、物理尺寸与采样间隔的取舍Matlab仿真光场本质上就是把连续光场离散成二维矩阵。矩阵的每个元素对应空间上的一个采样点所以第一步必须是搭好坐标系。这一步看似简单实际上决定了后面所有结果的可靠性。我常用的参数组合是网格点数 N 512 或 1024网格物理尺寸 L 3 mm 到 5 mm束腰半径 w0 0.4 mm 左右波长 lambda 632.8 nm氦氖激光器采样间隔 dx L / N。如果L3 mmN512则dx约5.86 μm。对可见光波段来说这个采样间隔足够描述束腰在百微米量级的高斯包络。经验法则是网格尺寸至少要覆盖束腰的4到5倍否则高斯包络在边界被截断仿真的光强图会出现一圈不自然的环而dx至少要小于光场最小特征尺度的一半。为什么要纠结采样间隔这关系到后面的傅里叶变换。角谱法传播中要用到fft2而fft2输出的频率范围由dx决定最大可表示空间频率是1/(2dx)大于这个频率的成分会混叠回低频形成所谓的“鬼影”或折叠条纹。如果网格太小最初看不出问题一旦做传播仿真高频混叠非常明显。2.2 用 meshgrid 构造径向和极角坐标Matlab里最常用的网格构造方式如下lambda 632.8e-9; % 波长 N 512; % 网格点数 L 3e-3; % 网格尺寸 w0 0.4e-3; % 束腰半径 l 1; % 拓扑荷数 p 0; % 径向指数 dx L / N; x (-N/2 : N/2-1) * dx; [X, Y] meshgrid(x, x); [Phi, R] cart2pol(X, Y);meshgrid生成的是二维网格坐标X和Y然后用cart2pol直接得到径向坐标R和极角Phi。这一步非常方便不需要手动写 atan2。有一个细节需要提醒x数组我习惯用(-N/2 : N/2-1)而不是linspace(-L/2, L/2, N)。两者很接近但前者保证相邻点间距精确为dx后者采样点之间间距略有差异在做fft2时容易引入额外的数值误差。这个差异很小但既然能做对就别留隐患。2.3 一个容易忽略的问题原点处的 0^0构造涡旋光束时我们会遇到(R.^abs(l))这样的项。当l0时R0处的0^0在Matlab里恰好返回1所以普通高斯光束不会出问题。但当l≠0时(R.^abs(l))在原点为0这正好是涡旋光束需要的中心振幅为零。麻烦往往出在负指数或除法运算上。有些人为了构造高阶涡旋会写(R.^(-abs(l)))或1./(R.^abs(l))这在原点会得到Inf。物理上涡旋光束的振幅不会发散所以这种写法是错的。正确做法是保留(R.^abs(l))在原点为零的特性后面再乘上高斯包络exp(-R.^2/w0^2)这样原点处振幅就是0乘1干净利落。3. 拉盖尔-高斯光束复振幅构造与可视化3.1 焦平面复振幅公式拉盖尔-高斯光束在z0处的复振幅可以写成E_lp(r, φ) A * (√2 r / w0)^|l| * L_p^|l|(2r²/w0²) * exp(-r²/w0²) * exp(ilφ)其中A是归一化常数如果只看相对强度可以省去L_p^|l|是关联拉盖尔多项式exp(-r²/w0²)是高斯包络exp(ilφ)是螺旋相位因子。很多教程为了省事只写p0的情形这时L_0^|l|恒为1公式简化成E_l0(r, φ) ∝ (√2 r / w0)^|l| * exp(-r²/w0²) * exp(ilφ)我做仿真时p0的模式用得最多因为实验中最常见的涡旋光就是p0的LG束。但为了让代码更通用建议把拉盖尔多项式也一并实现这样以后想要多环结构就能直接改参数。3.2 拉盖尔多项式的Matlab实现Matlab的符号数学工具箱里有laguerreL函数可以直接调用。但很多情况下并没有装这个工具箱而且符号计算在高分辨率网格上非常慢。我习惯用递推关系手写速度快而且可控L_0^k(x) 1L_1^k(x) 1 k - xL_{n1}^k(x) [(2n1k-x) * L_n^k(x) - (nk) * L_{n-1}^k(x)] / (n1)对应的Matlab代码function L assoc_laguerre(n, k, x) % 递推计算关联拉盖尔多项式 L_n^k(x) if n 0 L ones(size(x)); return; end L_prev ones(size(x)); % L_0 L_curr 1 k - x; % L_1 if n 1 L L_curr; return; end for m 1 : n-1 L_next ((2*m 1 k - x) .* L_curr ... - (m k) .* L_prev) / (m 1); L_prev L_curr; L_curr L_next; end L L_curr; end使用时只需要在构造光场的代码里写x 2 * R.^2 / w0^2; Lp assoc_laguerre(p, abs(l), x); E (sqrt(2)*R/w0).^abs(l) .* Lp .* exp(-R.^2/w0^2) .* exp(1i*l*Phi);需要注意的是递推公式中的x不能是二维矩阵吗完全可以Matlab的数组运算天然支持逐元素计算。这也是为什么Matlab非常适合做光场仿真整个横截面光场一次矩阵运算就出来了。3.3 光强与相位图的正确画法光强图很简单I abs(E).^2然后用imagesc显示。我常用自定义的颜色映射比如hot或parula因为默认的jet虽然好看但色带不单调容易让人误读强度差异。相位图稍微讲究一点。Matlab的angle函数返回的是[-π, π]范围内的包裹相位直接画出来会出现密密麻麻的跳变线这是正常现象不是bug。对于涡旋光束来说这些跳变线恰好就是从中心延伸到边界的那条线非常有特征。一个常用技巧是相位图用mod(angle(E), 2*pi)把范围转换到[0, 2π]再画螺旋结构更直观跳变线从-π到π的转换变成从2π到0的转换视觉上干净许多。figure(Position, [100 100 1200 500]); subplot(1,2,1); imagesc(x*1e3, x*1e3, abs(E).^2); axis image; colormap(gca, hot); title(光强分布); xlabel(x (mm)); ylabel(y (mm)); subplot(1,2,2); imagesc(x*1e3, x*1e3, mod(angle(E), 2*pi)); axis image; colormap(gca, hsv); title(相位分布); xlabel(x (mm)); ylabel(y (mm));画出图之后你应该能看到一个清晰的暗环或暗盘以及相位图上的螺旋结构。中心点的相位奇点就在螺旋的中心处。4. 干涉图样仿真用叉形条纹给涡旋光束“验明正身”4.1 平面波干涉叉形条纹的产生原理光强图上的暗斑可以承认涡旋光束但还不能完全确认它就是涡旋——任何中空光强都可能长这样。真正的“指纹”在干涉图里。把涡旋光束和一个倾斜的平面波叠加E_total E_l A * exp(i * kx * X)其中kx k * sinθ ≈ kθ是参考光的横向波矢。干涉图样的光强为I |E_l|² A² A * E_l * exp(-i kx X) A * conj(E_l) * exp(i kx X)交叉项携带了涡旋光束的螺旋相位因子exp(ilφ)所以干涉条纹会出现分叉结构。涡旋相位的等相位线绕中心一圈变化2πl但平面波的条纹是平行线两者叠加后在相位奇点处条纹被迫“断开再重新接上”产生分叉分叉的条数正好等于拓扑荷数l。这一段原理必须讲清楚因为很多人在仿真时看到分叉条纹觉得好看就完事了实际上分叉数量和方向才是关键。4.2 Matlab干涉图样仿真代码% 参考光参数 theta 0.008; % 倾斜角单位弧度大约0.46度 kx 2*pi/lambda * theta; % 横向波矢 E_ref exp(1i * kx * X); % 平面波参考光 % 干涉 E_interf E E_ref; I_interf abs(E_interf).^2; figure; imagesc(x*1e3, x*1e3, I_interf); axis image; colormap(gca, gray); title([l , num2str(l), 的叉形干涉条纹]);运行之后你会看到明暗相间的竖直条纹在中心附近出现了叉形结构。当l1时条纹中间出现一个分叉像字母Y倒过来当l2时分叉条数变成两条l的绝对值越大分叉越多。我一开始做干涉仿真时犯过一个错误参考光的倾角设得太小结果条纹间距比整个光斑还大分叉特征完全看不清。后来总结的经验是参考光的空间频率应该能让一个周期内装下至少4到5个条纹也就是让干涉条纹间距δ 2π/kx ≈ λ/θ大约为光斑尺寸的1/5到1/10。把θ设成0.005到0.02之间一般都有不错的效果。叉形条纹的朝向也能反映l的正负。沿着参考光倾斜的方向观察l0时叉形开口朝某个方向l0时朝另一个方向。这个特征实验上用来快速判断拓扑荷符号。4.3 用球面波干涉得到螺旋条纹除了平面波还可以用一个球面波作为参考光。球面波在傍轴近似下可以写成E_ref exp(1i * k * (X.^2 Y.^2) / (2 * d))其中d是参考点源到观察面的距离。球面波和涡旋光束干涉后由于球面波本身带有二次相位等相位线变成圆弧叠加后会产生螺旋条纹从中心往外看就像蜗牛壳一样。这种干涉图样在实验中也常见尤其在用马赫-曾德尔干涉仪时如果一臂是平面波一臂是涡旋光看到的就是叉形条纹如果用球面波看到的就是螺旋条纹。螺旋条纹的旋转方向和拓扑荷符号直接相关仿真时可以通过改变l的符号观察条纹绕向反转。这个操作很直观也能加深对拓扑荷方向的理解。5. 角谱法模拟自由空间传播看涡旋光束如何演变5.1 角谱传递函数的推导与近轴近似光场仿真最有意思的部分不是只看焦平面而是看它传一段距离之后长什么样。涡旋光束传几米之后会不会散了中心暗斑还在不在这两个问题都可以用角谱法回答。角谱法的思路很简单把初始光场E_in做二维傅里叶变换得到它在空间频率域的角谱自由空间传播只是给每个平面波分量叠加一个相位延迟再逆傅里叶变换就得到传播后的光场E_out。频域传递函数可以写成H(fx, fy) exp(i * 2π * z * sqrt(1/λ² - fx² - fy²))这是精确形式。在傍轴近似下也就是fx² fy²远小于1/λ²时可以做泰勒展开得到常用形式H(fx, fy) ≈ exp(i * k * z) * exp(-i * π * λ * z * (fx² fy²))其中k 2π/λz是传播距离。Matlab代码function Eout angular_spectrum(Ein, L, lambda, z) % 角谱法自由空间传播近轴近似 N size(Ein, 1); dx L / N; fx (-N/2 : N/2-1) / L; [FX, FY] meshgrid(fx, fx); k 2 * pi / lambda; H exp(1i * k * z) .* exp(-1i * pi * lambda * z * (FX.^2 FY.^2)); H ifftshift(H); Eout ifft2(fft2(Ein) .* H); end使用这个函数前一定要先自检把z设成0输出必须等于输入。如果不等说明频域网格构造或ifftshift用错了。频域网格fx (-N/2 : N/2-1) / L这一句是整个传播函数的关键。L是网格物理尺寸所以频率间隔是1/L最大频率是N/(2L) 1/(2dx)对应奈奎斯特频率。5.2 传播距离、采样条件与参数选择角谱法看起来简单但用起来有一堆限制。传播距离z过大时H矩阵中高频分量的相位旋转非常快导致数值上出现混叠。判断标准可以用菲涅耳衍射的采样条件对于光场横向范围L和波长λ最大可传播距离大致满足z_max ≈ N * dx² / λ以N512、dx5.86 μm、λ632.8 nm为例z_max ≈ 512 * (5.86e-6)² / 632.8e-9 ≈ 27.8 mm。也就是说上面这套参数只能传播不到3厘米。如果你要模拟光束传输20 cm甚至1 m就必须增大网格尺寸或减小网格点数。更好的做法是加一个缩放网格的传播算法比如两步角谱法或快速傅里叶分步法否则远距离传播一定会发散或产生折叠伪影。很多朋友说“仿真发散”其实绝大多数情况不是物理发散而是数值欠采样导致的频率混叠。光传着传着就从边界绕回来了看起来像从四个方向往中心钻那就是欠采样无疑了。5.3 传播过程中涡旋结构的稳定性用合适的参数跑一次传播你会发现涡旋光束在自由空间传播时中心暗斑会一直保持。虽然暗斑的尺寸会随着光束衍射而增大光强的空心轮廓也慢慢放大但相位奇点和拓扑荷数不会消失。这就是涡旋光束最迷人的地方——拓扑荷是拓扑保护的只要没有扰动破坏相位奇点它就一直在。我常用一个循环来观察不同距离处的光强和相位zlist [0, 10e-3, 20e-3, 30e-3]; figure; for m 1:4 Ez angular_spectrum(E, L, lambda, zlist(m)); subplot(2, 4, m); imagesc(x*1e3, x*1e3, abs(Ez).^2); axis image; colormap(gca, hot); title([z , num2str(zlist(m)*1e3), mm]); subplot(2, 4, m4); imagesc(x*1e3, x*1e3, mod(angle(Ez), 2*pi)); axis image; colormap(gca, hsv); end这个方法可以直观地看到涡旋光束的衍射行为。我建议大家都跑一遍这个程序它对建立“光束如何传播”的空间感帮助非常大。6. 高阶涡旋、叠加态与更多花样6.1 增大拓扑荷数 l光强环会变大吗同样的束腰w0下增大拓扑荷数l会发生什么从复振幅公式可以看到(√2 r / w0)^|l|这一项让中心区域的光强更早趋于零所以暗核半径会变大。具体来说p0的LG束光强极大值位置大约在r_peak w0 * sqrt(|l| / 2)所以l1时亮环半径约0.707w0l2时约w0l3时约1.225w0。这个公式可以用来快速估算高阶涡旋的光斑大小。相位图上的变化更明显l越大螺旋臂越多相位跳变线也越多。画出来的干涉图样叉形条纹的分叉数也会按l值增加。6.2 共轴叠加与涡旋劈裂涡旋光束叠加是OAM复用和模式分析的基础。把两个拓扑荷不同的涡旋光直接复振幅相加l1 1; l2 2; E1 (sqrt(2)*R/w0).^abs(l1) .* exp(-R.^2/w0^2) .* exp(1i*l1*Phi); E2 (sqrt(2)*R/w0).^abs(l2) .* exp(-R.^2/w0^2) .* exp(1i*l2*Phi); E_sum E1 E2;共轴叠加后的光强会出现方位角上的花瓣结构。比如l11和l2-1的叠加因为交叉项中有cos(2φ)光强呈对称的双瓣或四瓣结构中间可能还会出现额外的涡旋位错。相位图上可以看到多个相位奇点对这些奇点成对出现或消失是涡旋光场中非常有趣的动力学现象。两个涡旋如果横向错开一定距离再叠加光强图上可以看到两个暗核相位图上对应两个相位奇点。随着横向位移增加暗核也会分开得更明显。这个仿真对理解涡旋光场的拓扑结构和奇点演化学非常有用。6.3 径向指数 p 与多环结构拉盖尔-高斯光束除了拓扑荷l还有径向指数p。p0是单环p1、p2时会出现多个同心亮环环与环之间是径向暗线。p越大环数越多光场分布越复杂。用前面写的assoc_laguerre递推函数把p设成1或2就能轻松观察多环结构。这类光束在粒子囚禁和超分辨成像实验中有实际用途模型本身也非常适合入门练习。7. 仿真踩坑记录这些错误我几乎每次都遇到7.1 相位图显示“毛玻璃”效果第一次画相位图时如果直接用imagesc(angle(E))你会看到很多细碎的跳变线像毛玻璃一样完全看不出螺旋结构。这不代表光场有问题而是因为angle返回的相位被包裹在[-π, π]里螺旋相位每绕一圈就跳变一次。解决方法是用mod(angle(E), 2*pi)把范围改到[0, 2π]视觉上跳变从一条密线变成一条干净的“色带边界”。如果连这条边界都不想看到可以做相位解包裹但Matlab自带的unwrap只能沿一个方向解包对二维螺旋效果一般。我通常在演示代码里直接用mod简单且足够直观。7.2 传播后图样紊乱或从边界“绕回来”前面提到过传播距离超出采样条件后角谱法结果会出现混叠。判断标准是输出光场边缘是否出现高强度的重复结构。如果有说明频率域采样不足。处理办法有三个增大网格点数N不要觉得512就够了传播距离长时用1024甚至2048很正常减小最大传播距离改用缩放网格算法比如chirp z变换或两步角谱法。最简单粗暴的方法是直接减小dx、增大L并同时增大N但这样内存占用会明显上升毕竟三维矩阵在Matlab里并不便宜。我实际操作中会先用512点跑通逻辑再根据传播距离调大N。7.3 中心暗斑没有出现或光强溢出如果构造涡旋光束后光强图中心并没有暗斑最常见的原因是l设成了0只是普通高斯光束。另一个可能是(R.^abs(l))这一项在原点为0但整个网格分辨率太低暗核半径小于一个像素看起来就像没有暗斑。解决办法是增大N或增大l的值让暗核至少覆盖2到3个像素。光强溢出一般是因为没有归一化。构造E时exp(-R.^2/w0^2)把幅度限制在1以内所以光强峰值最多是1不太容易溢出。但如果你把多个模式叠加比如l从-3到3共7个模式全加起来峰值可能达到7画图时白花花一片。这时需要做归一化比如E E / max(abs(E(:)))或者把光强图用caxis限制显示范围。7.4 干涉条纹没分叉或对比度太低干涉条纹没分叉通常不是涡旋光的问题而是参考光设置不对。参考光横向波矢kx太小条纹间距太大整个光斑内只有一两根条纹分叉特征显示不出来。把倾斜角θ调大一些或者直接减小参考光波长与光斑尺寸的比值都能改善。对比度低则往往是因为涡旋光和参考光的振幅相差太大。比如涡旋光中心暗边缘亮参考光振幅取1叠加后交叉项的调制深度被淹没。我把参考光振幅调成与涡旋光峰值振幅同一量级干涉条纹立刻清晰许多。7.5 常见问题汇总表现象可能原因解决方法光强图中心没有暗斑l0或网格分辨率不足检查l值增大N或l相位图一片毛玻璃包裹相位显示方式不当用mod(angle(E), 2*pi)传播后图样折叠混叠传播距离或采样间隔不匹配减小z增大N用缩放网格法干涉条纹少见分叉参考光倾角太小增大θ到0.005~0.02干涉条纹对比度低两束光振幅差异过大调整参考光振幅光强分布不对称网格中心不是偶数索引导致偏移用(-N/2:N/2-1)构造坐标出现NaNR0处除零或负指数幂避免1./R.^abs(l)的写法结尾的一些实际操作体会这套仿真流程我自己反复用过很多次从最基本的拉盖尔-高斯光束到干涉图样再到传播演化每次跑通都有新的理解。最想提醒你的是不要一口气追求复杂功能先把l1、p0的单涡旋在z0处的光强和相位搞清楚再逐步加干涉、加传播、加高阶项。涡旋光束的物理图像很大程度上是靠这一步步的仿真建立起来的。我个人比较推荐的一个练习是用Matlab生成一个l3的涡旋光束仿真它与平面波干涉的叉形条纹然后把图存下来对比实验室里空间光调制器产生的实验结果。你会发现最重要的差异往往不在叉形条纹本身而在背景噪声、光斑对称性和条纹衬比度这些细节会促使你回头修正仿真参数也让你对光路调校有更具体的感知。如果你已经跑通了这些基础内容后面可以试着仿真分数阶涡旋光束拓扑荷是非整数、偏振涡旋矢量光束、或者用角谱法模拟涡旋光束通过薄透镜的聚焦过程。每往前迈一步都会有新的坑和新的乐趣。