尧图网站建设 
简介面向激光物理、光学工程及光通信方向的学生与研究者这份资源以MATLAB程序演示拉盖尔-高斯光束p01模态的模拟与可视化帮助理解LG光束的径向结构、轨道角动量及中心暗斑等概念相关特性在光学通信、量子光学与光镊领域有重要应用。压缩包体积小巧共2个文件包括.m脚本和.png结果图整体仅66KB脚本直接运行即可生成光束截面强度图无需额外配置。已有2190人学习适合本科高年级或研究生作为课堂辅助、课程设计及自学参考。代码依据拉盖尔多项式与高斯项构建复光场清晰呈现“甜甜圈”状的强度分布读者还可修改p、l参数观察不同模式变化为后续研究涡旋光束、光镊或量子光学应用提供可复现的算法基础。1. 从一张环状光斑说起LG01拉盖尔-高斯光束的MATLAB模拟能做什么在光学实验室里想看到LG01这种带有轨道角动量的涡旋光束通常需要空间光调制器、涡旋波片或计算全息图成本高、调光路也费时间。而用MATLAB模拟LG01你可以在几分钟内得到强度环、螺旋相位和干涉叉形条纹甚至在还没有搭建光学系统前就评估光束参数对面型的影响。这个任务对于学习物理光学、设计涡旋光实验或者做光通信复用仿真的工程师都非常常见。LG01是拉盖尔-高斯模式中p0、l1的特例它看起来像一个明亮的空心圆环实际上每个点在传播方向上都携带一份轨道角动量。下文会从数学表达式出发逐步给出可在MATLAB里直接运行的代码和参数设置并指出常见的数值陷阱。2. 拉盖尔-高斯光束的数学结构从复振幅公式到LG01的相位奇点2.1 LG_{p,l}模式的复振幅与拉盖尔多项式拉盖尔-高斯光束是傍轴波动方程在柱坐标下的一组完备解。它用两个整数来标记模式径向指数 p0、1、2…和角向指数 l也称拓扑荷可为正负整数。复振幅在 z 处可以写成$$E_{p,l}(r,\phi,z) A_{p,l} \frac{1}{w(z)} \left(\frac{\sqrt{2}r}{w(z)}\right)^{|l|} L_p^{|l|}\left(\frac{2r^2}{w(z)^2}\right) \exp\left(-\frac{r^2}{w(z)^2}\right) \exp\left(-i\frac{k r^2}{2R(z)}\right) \exp(-i l \phi) \exp(i (2p|l|1)\psi(z))$$其中 $A_{p,l} \sqrt{\frac{2p!}{\pi (p|l|)!}}$ 是归一化系数。$L_p^{|l|}$ 是关联拉盖尔多项式$w(z)$ 是光束半径$R(z)$ 是等相面曲率半径$\psi(z)$ 是Gouy相位。这些参数都与束腰 w0、波长 λ 和 z 相关$w(z)w_0\sqrt{1(z/z_R)^2}$$R(z)z(1(z_R/z)^2)$$z_R\pi w_0^2/\lambda$$\psi(z)\arctan(z/z_R)$。这个公式里最容易被忽视的是 $(\sqrt{2} r / w)^{|l|}$ 因子。当 l 不为零时它导致光强在 r→0 处趋于零从而形成空心。同时 $\exp(-i l \phi)$ 给出了绕轴旋转的螺旋相位这是轨道角动量的来源。在MATLAB中实现的关键不是自己推导拉盖尔多项式而是利用递推关系。关联拉盖尔多项式满足以下递推式$L_0^k(x)1$$L_1^k(x)1k-x$$(n1)L_{n1}^k(x) (2n1k-x)L_n^k(x) - (nk)L_{n-1}^k(x)$。这个递推式在数值上很稳定适合任意 p 和 k。对于 LG01 来说 p0多项式等于1所以代码可以简化很多。2.2 束腰、瑞利距离与Gouy相位参数间的耦合很多人在模拟时只关心强度分布但“LG01”里的物理参数会直接影响光斑尺寸和相位曲率。束腰 w0 决定瑞利距离 $z_R \pi w_0^2/\lambda$。如果 z 远大于 zR光束半径会近似按 $\lambda z/(\pi w_0)$ 扩散如果 z0则 $R(z)$ 无穷大等相面是平面。在仿真中如果你在束腰处生成了LG01则复振幅中的曲率项 $\exp(-i kr^2/(2R))$ 等于1可以省略。但这只对 z0 成立。Gouy相位项 $\exp(i(2p|l|1)\psi(z))$ 通常不会影响强度但会影响干涉或相干叠加。比如将 LG01 与 LG00 (p0,l0) 叠加时两者经历的Gouy相位不同传播过程中花瓣图样会绕轴旋转这是Gouy相位在干涉上的直接体现。因此在编写通用函数时我一般会把 z、zR、R 和 ψ 都作为可计算参数传进去而不是写死为束腰处的特例。这样同一个函数既能生成初始面也能生成传播后的光束后面做FFT角谱对比时也更方便。下列参数表给出了一个典型LG01仿真需要关注的量参数符号LG01推荐值物理影响波长λ632.8 nm决定瑞利距离和光斑尺度束腰半径w01 mm决定空心环大小和发散角传播距离z0束腰0时无球面相位便于验证拓扑荷l1决定螺旋相位阶数径向指数p0决定径向节线数p0单环2.3 LG01的强度为什么是空心的相位为什么是螺旋的当 p0, l1 时拉盖尔多项式部分等于1强度公式变为 $I(r) \propto (2r^2/w^2) \exp(-2r^2/w^2)$。这个函数在 $rw/\sqrt{2}$ 处达到最大值中心处为零因此形成暗核的空心环。环的半径约等于 $w/\sqrt{2}$即强度峰值位置注意不要和光束半径 w 混淆。这个峰值位置是实验上测量拓扑荷的重要参考。相位方面$\exp(-i l \phi)$ 中包含一个沿方位角从0到2π变化的螺旋面。l1 意味着相位绕轴一圈刚好变化 2π中心处相位不定、强度为零。若用 unwrap 相位图会看到从0到2π的连续渐变并且用相位图配色时能明显看到“旋转楼梯”状结构。这个螺旋相位的方向由 l 的正负决定l1 时是逆时针l-1 时顺时针但光斑形状完全相同。理解了这一点你就知道为什么生成 LG01 时不能只画强度必须同时检查相位图否则根本无法确认螺旋拓扑荷和旋向是否正确。在后面的代码中我会把强度图和相位图放在同一个 figure 里。3. 在MATLAB中搭建LG01光束的网格与场函数3.1 设计坐标网格物理尺寸与采样点数模拟LG01的第一步是建立一个二维网格。物理尺寸要保证光束完全落在网格内一般取网格边长 L 为光束半径的 4~6 倍。对于 w01 mm 的LG01峰值半径约 0.707 mm所以 L 取 5 mm 足够。采样点数 N 建议取 256 或 512既能让环的圆滑度足够又不会让矩阵过大。网格间距 dxL/N要满足采样定理尤其是后续要用FFT传播时dx 与波长、传播距离要满足关系。核心代码如下% LG01_laguerre_gaussian.m lambda 632.8e-9; % 波长[m]常用氦氖激光器红光 w0 1e-3; % 束腰半径[m] L 6e-3; % 网格边长[m] N 512; % 采样点数2的幂 x linspace(-L/2, L/2, N); y linspace(-L/2, L/2, N); % 注意meshgrid的默认输出是 X 按行变化Y 按列变化 % 与图像坐标一致直接使用即可。 [X, Y] meshgrid(x, y); [Phi, R] cart2pol(X, Y); % 极坐标Phi 是方位角meshgrid生成的 X 和 Y 是 N×N 矩阵。cart2pol返回的 Phi 范围是 [-π, π]而公式里的 φ 通常取 [0, 2π) 也没有问题因为相位差只在指数上但为了 view 时连续建议对 Phi 加一个 pi 或直接用 unwrap。R 是径向距离矩阵。参数选择上L 太小会截断光束导致边缘出现衍射纹L 太大则会浪费采样点。一个经验公式是 $L \geq 6w(z)$而 N 取 512 时内存约 4 MB per double matrix非常轻量。采样间隔 dx 需要小于某个极限特别是做角谱传播时dx 必须满足 $dx \leq \lambda z / N$ 可能不一定容易但这里不展开第五章会说明。3.2 实现拉盖尔多项式用循环还是符号计算对于任意 p 和 l推荐用递推循环计算广义拉盖尔多项式。下面这个函数返回 $2r^2/w^2$ 处对应的多项式的值。因为我们要在矩阵上计算所以函数接收矩阵 x 并返回同尺寸矩阵。function Lval laguerre_assoc(p, k, x) % 关联拉盖尔多项式递推 % p 为非负整数k 为标量 |l|x 可为矩阵 L0 ones(size(x)); % n0 if p 0 Lval L0; return; end L1 1 k - x; % n1 if p 1 Lval L1; return; end Ln_prev L0; Ln L1; for n 1:(p-1) Ln_next ((2*n 1 k - x) .* Ln - (n k) .* Ln_prev) / (n 1); Ln_prev Ln; Ln Ln_next; end Lval Ln; end代码里的递推关系与上一章公式一致。这里用ones(size(x))而不是标量1是为了让整块矩阵都能参与运算。由于表达式里只涉及矩阵乘法和除法运行速度非常快即便 p 取到两位数也不会有性能问题。符号计算symslaguerreL只适合推导不适合放进 512×512 的循环里所以不推荐。另外注意 MATLAB 的laguerreL函数在高版本中支持符号输入但它返回的是符号对象转换为 double 后再填充矩阵会损失性能。自己写递推函数是更可控的方案。3.3 生成LG01的复振幅并绘制强度与相位有了网格和拉盖尔函数生成LG01的复振幅就简单了。以下代码放在同一个脚本中% 物理参数 p 0; l 1; k 2*pi/lambda; z 0; % 束腰位置 zR pi*w0^2/lambda; wz w0 * sqrt(1 (z/zR)^2); % 曲率项在z0时为1省去直接乘以exp(0) psi atan(z/zR); % 极坐标的径向变量 rho sqrt(2) * R / wz; x_rho 2 * R.^2 / wz^2; % 拉盖尔多项式自变量 Lval laguerre_assoc(p, abs(l), x_rho); % 复振幅z0时省略曲率项 E sqrt(2*factorial(p)/(pi*factorial(pabs(l)))) ... / wz * rho.^abs(l) .* Lval .* exp(-R.^2/wz^2) ... .* exp(-1i * l * Phi); I abs(E).^2; I I / max(I(:)); % 归一化方便显示 % 相位 phase_E mod(angle(E), 2*pi); % 折叠到[0,2pi) figure(Color,white); subplot(1,2,1); imagesc(x*1e3, y*1e3, I); axis image; colormap(hot); colorbar; xlabel(x (mm)); ylabel(y (mm)); title(LG01 Intensity); subplot(1,2,2); imagesc(x*1e3, y*1e3, phase_E); axis image; colormap(jet); colorbar; xlabel(x (mm)); ylabel(y (mm)); title(LG01 Phase (wrapped));注意l是1所以rho.^abs(l)与rho等价。exp(-R.^2/wz^2)在束腰处是 $\exp(-r^2/w_0^2)$。相乘之后得到复场。phase_E是折叠后的相位看起来是从0到2π的调色板中心由于强度为零相位噪声会被放大但那是正常的。参数上factorial(pabs(l))对于 p0,l1 是 1所以系数简化为 $\sqrt{2/\pi}/w_0$。如果 z≠0还要补乘 $\exp(-i k R^2/(2R_z))$并将 Rz 计算为 $z(1(z_R/z)^2)$。这里省略是因为 z0 时 Rz 无穷大避免除零。3.4 关键参数表从波长到束腰的推荐取值下面给出常用参数范围特别是当你希望生成的光斑直接对应真实实验时参数符号实验室常用值仿真建议波长λ532 nm / 632.8 nm / 1064 nm对应实际激光器束腰半径w00.5~2 mm越小环越小但要知道分子网格边长L—6~8倍w(z)采样点数N—256~1024传播距离z0~数个zR初始面设为0最省事拓扑荷l1, 3, 5l的正负决定旋向径向指数p0, 1, 2增加环数上表里“仿真建议”列的 w0 越小空心环峰值位置就越靠近中心但如果 w0 小于几个网格间距中心会因采样不足而出现菱形畸变。因此模拟时要保证 $w_0/dx \geq 20$。如果 w01mmN512L6mm则 dx11.7μmw0/dx≈85足够。4. 验证LG01的核心特征干涉条纹、轨道角动量与模式分析4.1 平面波干涉用叉形条纹确认拓扑荷生成LG01后最直接的验证是让它与沿x方向倾斜的平面波干涉。干涉场强度为 $|E A \exp(i k x \sin\theta)|^2$其中A是参考光振幅θ是倾斜角。干涉图在涡旋中心会出现典型的叉形条纹——条纹分叉数等于拓扑荷的绝对值。下面的代码演示了这一过程% 参考平面波沿x方向传播并带一个横向相位梯度 theta 0.5e-3; % 倾斜角单位 rad E_ref exp(1i * k * X * theta); % 相当于叠加横向波矢 k*theta E_int E / max(abs(E(:))) 0.8 * E_ref; % 控制相对振幅 I_int abs(E_int).^2; figure; imagesc(x*1e3, y*1e3, I_int); axis image; colormap(gray); colorbar; xlabel(x (mm)); ylabel(y (mm)); title(Interference of LG01 with a tilted plane wave);这里的k*X*theta容易理解倾斜角 θ 使平面波具有横向相位变化。条纹间距为 λ/θ当 θ0.5mrad、λ633nm 时约 1.27mm在 6mm 网格内能容纳约4条亮纹足够看到叉形。实际调节 θ 时注意不要让条纹间距小于3个像素否则叉形看不清。叉形中间的暗线分叉方向取决于 l 的符号——这比强度图更可靠。4.2 数值计算轨道角动量密度拉盖尔-高斯模式的一个可测量特征是每光子携带 $l\hbar$ 的角动量。在数值上可以直接计算相位矩阵在圆周上的差分和。因为 l1 时绕中心一圈的总相位变化应为 $2\pi$。用MATLAB可以取一个半径 r0 的圆环上的相位值unwrap后做线性拟合获得斜率。方法对 l 较大时依然准确但要求环半径不能太小否则中心奇点附近的相位噪声会干扰。下面给出一个简单的相位斜率检验代码% 取半径接近强度峰值位置的圆环 r0 w0 / sqrt(2); % 峰值半径处 t linspace(0, 2*pi, 720); % 角度向量 cx 0; % 光束中心假设在原点 pts [cx r0*cos(t); cx r0*sin(t)]; % 2x720 % 双线性插值提取相位 phase_ring interp2(X, Y, phase_E, pts(1,:), pts(2,:), cubic); % 去除NaN中心奇点附近没有有效相位 phase_ring unwrap(phase_ring); % 线性拟合相位 slope * t offset coeff polyfit(t, phase_ring, 1); slope coeff(1); fprintf(Phase slope around ring: %.4f (should be %.4f)\n, slope, l);unwrap会自动补偿跳变得到的斜率等于 l。这里用interp2插值而不要直接取Phit的索引因为圆环不一定正好落在网格点上。slope应接近 1.0误差小于0.1可以认为模式正确。峰值半径 r0w0/√2 是LG01强度最大处信噪比最高。4.3 扫描拓扑荷l与径向指数p观察模式演变进一步可以让 l 或 p 在一个循环中变化验证强度环数量和相位奇点数的关系。一个典型脚本figure(Color,white); for l_val [1 2 3] % 将第3章的核心代码封装成函数lgBeam(p,l,w0,lambda,z,N,L) E_l lgBeam(0, l_val, w0, lambda, 0, N, L); subplot(2,3,l_val); imagesc(x*1e3, y*1e3, abs(E_l).^2); axis image; colormap(hot); title(sprintf(l %d, p 0, l_val)); end这里省略生成函数实际可以把第3章代码封装成lgBeam(p,l,w0,lambda,z,N,L)。l 绝对值越大空心环半径越大相位奇点强度保持不变但仍为零。p 增大时会出现 p1 个同心环并且环与环之间有相位节点径向相位反转。下表总结了几种常见 $LG_{p,l}$ 的可见特征模式pl强度环数相位奇点LG00001实心无LG01011空心1拓扑荷1LG02021空心半径更大1拓扑荷2LG11112同心环1拓扑荷1扫描 l 可以帮助你直观理解拓扑荷对环半径的影响环半径正比于 $\sqrt{|l|} \cdot w_0/\sqrt{2}$大约如此。对于 l3峰值位置约在 $\sqrt{3}\times 0.707 w_0 \approx 1.22 w_0$。这个关系在实验上往往被用来从测量强度环半径推算拓扑荷。5. 进阶LG01的传播仿真与模式叠加5.1 用FFT角谱法传播LG01光束实际仿真中除了观察束腰处的光斑你往往还需要知道 LG01 在传播几个瑞利距离后的样子。常见做法是用角谱衍射传播。在 MATLAB 中实现% 角谱传递函数 dx L / N; fx (-N/2 : N/2-1) / L; [FX, FY] meshgrid(fx, fx); % 角谱传递函数自由空间 H exp(1i * k * z * sqrt(1 - (lambda*FX).^2 - (lambda*FY).^2)); % 初始场E前面生成的LG01 E0 E; Ef ifft2(ifftshift(fftshift(fft2(E0)) .* H)); % 显示 I_prop abs(Ef).^2;这里有一个很关键的细节fft2的输出频率中心在数组第一个元素需要用fftshift移动到中心而ifft2之前再用ifftshift恢复。如果少一步光斑会偏移半格。这里的H在sqrt内可能为负倏逝波一般可以忽略也可以将负值置零。在自由空间传播、距离较小时没有大问题。传播距离 z 的单位与网格尺度一致。若 z 太大角谱中的高频相位变化过快会出现混叠。经验法则是保证 z 远小于 $N dx^2 / \lambda$。例如 N512, dx11.7μm, λ633nm则 $N dx^2 / \lambda \approx 0.111m$即传播距离超过 11cm 时就需要加大网格或增大 dx。你在模拟一个几十厘米甚至1米后的LG01时要重新设计网格尺寸。5.2 叠加LG01与LG00生成横向花瓣结构拉盖尔-高斯模式之间是正交的但叠加后会因多普勒相位差产生干涉图案。例如将 LG01 与 LG00 同轴叠加会得到一个沿方位角变化的光强分布花瓣数等于两个模式拓扑荷之差。实现起来只需要E_lg00 lgBeam(0, 0, w0, lambda, 0, N, L); % 假定函数已封装 E_sum E / max(abs(E(:))) 1.0 * E_lg00 / max(abs(E_lg00(:))); I_sum abs(E_sum).^2; imshow(I_sum, []);这里的相对振幅比会影响对比度。当两个模式幅度相等时花瓣达到最大对比度如果想模拟实验上常见的“不完全对准”可以给其中一个模式加很小的横向位移花瓣会变得左右不对称。这种叠加常被用于产生径向偏振类似结构的矢量光束或者用于测量模式之间的Gouy相位差。注意由于 LG00 的中心实心叠加图不再是空心而是被切割成两瓣或更多瓣。l 差为1时是两瓣像“花生”。这在LG01作为纠缠光子源调试中常被用来判断空间模式是否匹配。5.3 与实验对照相机积分效应与像素化实验上的CCD或CMOS相机记录的并不是精确的光强点而是每个像素上的积分强度。模拟中如果直接显示abs(E).^2会忽略像素的积分效应。要更接近实验结果可以用conv2对强度平面做一次“像素平均”。比如让模拟的每个网格对应一个像素可先计算局部均值I_sim abs(E).^2; fs 4; % 假设每个物理像素对应4×4个网格 kernel ones(fs, fs) / fs^2; I_camera conv2(I_sim, kernel, same);在实际对照实验时还需要考虑相机噪声、位深量化等。最省事的方法是先用imresize把强度图降到相机分辨率再显示。这样做不会改变结论但能让你在用MATLAB写报告或论文插图时视觉上更接近拍摄照片。6. 把LG01仿真写扎实的3个实用检查6.1 检查数值归一化与总能量守恒每次都验证一下总能量是否恒定。对初始场和传播后的场分别做sum(abs(E(:)).^2) * dx^2在角谱传播中能量应几乎不变。如果能量偏离超过1%多半是网格尺寸不够或传递函数在边缘被截断。归一化到峰值为1很容易但能量守恒更关键。6.2 避开中心奇点的采样错误LG01中心相位奇异相位上的值是未定义的。在angle计算后中心会出现随机数这是正常的。如果强行用unwrap从中心开始解卷会产生大幅振荡。避免方法是先遮蔽半径小于 0.1w0 的区域再求相位。而插值环提取相位时选择的半径不应小于峰值半径否则插值到奇点附近会被NaN污染。6.3 用单精度还是双精度生成初始场时用 double 即可但对于超大规模网格N2048以上或者在循环里做上千次传播可以考虑把场转成 single 以节省一半内存。注意fft2支持 single 输入但传递函数 H 里同时包含 k 和 z乘积精度在传播距离很大时可能不足所以传播时建议仍用 double。基本的原则是内存可控时全程 double兼顾速度和精度只有反复调参的大批量扫描时才用 single。本文还有配套的精品资源点击获取 拉盖尔-高斯光束LG01的MATLAB仿真：涡旋光场与轨道角动量分析
简介面向激光物理、光学工程及光通信方向的学生与研究者这份资源以MATLAB程序演示拉盖尔-高斯光束p01模态的模拟与可视化帮助理解LG光束的径向结构、轨道角动量及中心暗斑等概念相关特性在光学通信、量子光学与光镊领域有重要应用。压缩包体积小巧共2个文件包括.m脚本和.png结果图整体仅66KB脚本直接运行即可生成光束截面强度图无需额外配置。已有2190人学习适合本科高年级或研究生作为课堂辅助、课程设计及自学参考。代码依据拉盖尔多项式与高斯项构建复光场清晰呈现“甜甜圈”状的强度分布读者还可修改p、l参数观察不同模式变化为后续研究涡旋光束、光镊或量子光学应用提供可复现的算法基础。1. 从一张环状光斑说起LG01拉盖尔-高斯光束的MATLAB模拟能做什么在光学实验室里想看到LG01这种带有轨道角动量的涡旋光束通常需要空间光调制器、涡旋波片或计算全息图成本高、调光路也费时间。而用MATLAB模拟LG01你可以在几分钟内得到强度环、螺旋相位和干涉叉形条纹甚至在还没有搭建光学系统前就评估光束参数对面型的影响。这个任务对于学习物理光学、设计涡旋光实验或者做光通信复用仿真的工程师都非常常见。LG01是拉盖尔-高斯模式中p0、l1的特例它看起来像一个明亮的空心圆环实际上每个点在传播方向上都携带一份轨道角动量。下文会从数学表达式出发逐步给出可在MATLAB里直接运行的代码和参数设置并指出常见的数值陷阱。2. 拉盖尔-高斯光束的数学结构从复振幅公式到LG01的相位奇点2.1 LG_{p,l}模式的复振幅与拉盖尔多项式拉盖尔-高斯光束是傍轴波动方程在柱坐标下的一组完备解。它用两个整数来标记模式径向指数 p0、1、2…和角向指数 l也称拓扑荷可为正负整数。复振幅在 z 处可以写成$$E_{p,l}(r,\phi,z) A_{p,l} \frac{1}{w(z)} \left(\frac{\sqrt{2}r}{w(z)}\right)^{|l|} L_p^{|l|}\left(\frac{2r^2}{w(z)^2}\right) \exp\left(-\frac{r^2}{w(z)^2}\right) \exp\left(-i\frac{k r^2}{2R(z)}\right) \exp(-i l \phi) \exp(i (2p|l|1)\psi(z))$$其中 $A_{p,l} \sqrt{\frac{2p!}{\pi (p|l|)!}}$ 是归一化系数。$L_p^{|l|}$ 是关联拉盖尔多项式$w(z)$ 是光束半径$R(z)$ 是等相面曲率半径$\psi(z)$ 是Gouy相位。这些参数都与束腰 w0、波长 λ 和 z 相关$w(z)w_0\sqrt{1(z/z_R)^2}$$R(z)z(1(z_R/z)^2)$$z_R\pi w_0^2/\lambda$$\psi(z)\arctan(z/z_R)$。这个公式里最容易被忽视的是 $(\sqrt{2} r / w)^{|l|}$ 因子。当 l 不为零时它导致光强在 r→0 处趋于零从而形成空心。同时 $\exp(-i l \phi)$ 给出了绕轴旋转的螺旋相位这是轨道角动量的来源。在MATLAB中实现的关键不是自己推导拉盖尔多项式而是利用递推关系。关联拉盖尔多项式满足以下递推式$L_0^k(x)1$$L_1^k(x)1k-x$$(n1)L_{n1}^k(x) (2n1k-x)L_n^k(x) - (nk)L_{n-1}^k(x)$。这个递推式在数值上很稳定适合任意 p 和 k。对于 LG01 来说 p0多项式等于1所以代码可以简化很多。2.2 束腰、瑞利距离与Gouy相位参数间的耦合很多人在模拟时只关心强度分布但“LG01”里的物理参数会直接影响光斑尺寸和相位曲率。束腰 w0 决定瑞利距离 $z_R \pi w_0^2/\lambda$。如果 z 远大于 zR光束半径会近似按 $\lambda z/(\pi w_0)$ 扩散如果 z0则 $R(z)$ 无穷大等相面是平面。在仿真中如果你在束腰处生成了LG01则复振幅中的曲率项 $\exp(-i kr^2/(2R))$ 等于1可以省略。但这只对 z0 成立。Gouy相位项 $\exp(i(2p|l|1)\psi(z))$ 通常不会影响强度但会影响干涉或相干叠加。比如将 LG01 与 LG00 (p0,l0) 叠加时两者经历的Gouy相位不同传播过程中花瓣图样会绕轴旋转这是Gouy相位在干涉上的直接体现。因此在编写通用函数时我一般会把 z、zR、R 和 ψ 都作为可计算参数传进去而不是写死为束腰处的特例。这样同一个函数既能生成初始面也能生成传播后的光束后面做FFT角谱对比时也更方便。下列参数表给出了一个典型LG01仿真需要关注的量参数符号LG01推荐值物理影响波长λ632.8 nm决定瑞利距离和光斑尺度束腰半径w01 mm决定空心环大小和发散角传播距离z0束腰0时无球面相位便于验证拓扑荷l1决定螺旋相位阶数径向指数p0决定径向节线数p0单环2.3 LG01的强度为什么是空心的相位为什么是螺旋的当 p0, l1 时拉盖尔多项式部分等于1强度公式变为 $I(r) \propto (2r^2/w^2) \exp(-2r^2/w^2)$。这个函数在 $rw/\sqrt{2}$ 处达到最大值中心处为零因此形成暗核的空心环。环的半径约等于 $w/\sqrt{2}$即强度峰值位置注意不要和光束半径 w 混淆。这个峰值位置是实验上测量拓扑荷的重要参考。相位方面$\exp(-i l \phi)$ 中包含一个沿方位角从0到2π变化的螺旋面。l1 意味着相位绕轴一圈刚好变化 2π中心处相位不定、强度为零。若用 unwrap 相位图会看到从0到2π的连续渐变并且用相位图配色时能明显看到“旋转楼梯”状结构。这个螺旋相位的方向由 l 的正负决定l1 时是逆时针l-1 时顺时针但光斑形状完全相同。理解了这一点你就知道为什么生成 LG01 时不能只画强度必须同时检查相位图否则根本无法确认螺旋拓扑荷和旋向是否正确。在后面的代码中我会把强度图和相位图放在同一个 figure 里。3. 在MATLAB中搭建LG01光束的网格与场函数3.1 设计坐标网格物理尺寸与采样点数模拟LG01的第一步是建立一个二维网格。物理尺寸要保证光束完全落在网格内一般取网格边长 L 为光束半径的 4~6 倍。对于 w01 mm 的LG01峰值半径约 0.707 mm所以 L 取 5 mm 足够。采样点数 N 建议取 256 或 512既能让环的圆滑度足够又不会让矩阵过大。网格间距 dxL/N要满足采样定理尤其是后续要用FFT传播时dx 与波长、传播距离要满足关系。核心代码如下% LG01_laguerre_gaussian.m lambda 632.8e-9; % 波长[m]常用氦氖激光器红光 w0 1e-3; % 束腰半径[m] L 6e-3; % 网格边长[m] N 512; % 采样点数2的幂 x linspace(-L/2, L/2, N); y linspace(-L/2, L/2, N); % 注意meshgrid的默认输出是 X 按行变化Y 按列变化 % 与图像坐标一致直接使用即可。 [X, Y] meshgrid(x, y); [Phi, R] cart2pol(X, Y); % 极坐标Phi 是方位角meshgrid生成的 X 和 Y 是 N×N 矩阵。cart2pol返回的 Phi 范围是 [-π, π]而公式里的 φ 通常取 [0, 2π) 也没有问题因为相位差只在指数上但为了 view 时连续建议对 Phi 加一个 pi 或直接用 unwrap。R 是径向距离矩阵。参数选择上L 太小会截断光束导致边缘出现衍射纹L 太大则会浪费采样点。一个经验公式是 $L \geq 6w(z)$而 N 取 512 时内存约 4 MB per double matrix非常轻量。采样间隔 dx 需要小于某个极限特别是做角谱传播时dx 必须满足 $dx \leq \lambda z / N$ 可能不一定容易但这里不展开第五章会说明。3.2 实现拉盖尔多项式用循环还是符号计算对于任意 p 和 l推荐用递推循环计算广义拉盖尔多项式。下面这个函数返回 $2r^2/w^2$ 处对应的多项式的值。因为我们要在矩阵上计算所以函数接收矩阵 x 并返回同尺寸矩阵。function Lval laguerre_assoc(p, k, x) % 关联拉盖尔多项式递推 % p 为非负整数k 为标量 |l|x 可为矩阵 L0 ones(size(x)); % n0 if p 0 Lval L0; return; end L1 1 k - x; % n1 if p 1 Lval L1; return; end Ln_prev L0; Ln L1; for n 1:(p-1) Ln_next ((2*n 1 k - x) .* Ln - (n k) .* Ln_prev) / (n 1); Ln_prev Ln; Ln Ln_next; end Lval Ln; end代码里的递推关系与上一章公式一致。这里用ones(size(x))而不是标量1是为了让整块矩阵都能参与运算。由于表达式里只涉及矩阵乘法和除法运行速度非常快即便 p 取到两位数也不会有性能问题。符号计算symslaguerreL只适合推导不适合放进 512×512 的循环里所以不推荐。另外注意 MATLAB 的laguerreL函数在高版本中支持符号输入但它返回的是符号对象转换为 double 后再填充矩阵会损失性能。自己写递推函数是更可控的方案。3.3 生成LG01的复振幅并绘制强度与相位有了网格和拉盖尔函数生成LG01的复振幅就简单了。以下代码放在同一个脚本中% 物理参数 p 0; l 1; k 2*pi/lambda; z 0; % 束腰位置 zR pi*w0^2/lambda; wz w0 * sqrt(1 (z/zR)^2); % 曲率项在z0时为1省去直接乘以exp(0) psi atan(z/zR); % 极坐标的径向变量 rho sqrt(2) * R / wz; x_rho 2 * R.^2 / wz^2; % 拉盖尔多项式自变量 Lval laguerre_assoc(p, abs(l), x_rho); % 复振幅z0时省略曲率项 E sqrt(2*factorial(p)/(pi*factorial(pabs(l)))) ... / wz * rho.^abs(l) .* Lval .* exp(-R.^2/wz^2) ... .* exp(-1i * l * Phi); I abs(E).^2; I I / max(I(:)); % 归一化方便显示 % 相位 phase_E mod(angle(E), 2*pi); % 折叠到[0,2pi) figure(Color,white); subplot(1,2,1); imagesc(x*1e3, y*1e3, I); axis image; colormap(hot); colorbar; xlabel(x (mm)); ylabel(y (mm)); title(LG01 Intensity); subplot(1,2,2); imagesc(x*1e3, y*1e3, phase_E); axis image; colormap(jet); colorbar; xlabel(x (mm)); ylabel(y (mm)); title(LG01 Phase (wrapped));注意l是1所以rho.^abs(l)与rho等价。exp(-R.^2/wz^2)在束腰处是 $\exp(-r^2/w_0^2)$。相乘之后得到复场。phase_E是折叠后的相位看起来是从0到2π的调色板中心由于强度为零相位噪声会被放大但那是正常的。参数上factorial(pabs(l))对于 p0,l1 是 1所以系数简化为 $\sqrt{2/\pi}/w_0$。如果 z≠0还要补乘 $\exp(-i k R^2/(2R_z))$并将 Rz 计算为 $z(1(z_R/z)^2)$。这里省略是因为 z0 时 Rz 无穷大避免除零。3.4 关键参数表从波长到束腰的推荐取值下面给出常用参数范围特别是当你希望生成的光斑直接对应真实实验时参数符号实验室常用值仿真建议波长λ532 nm / 632.8 nm / 1064 nm对应实际激光器束腰半径w00.5~2 mm越小环越小但要知道分子网格边长L—6~8倍w(z)采样点数N—256~1024传播距离z0~数个zR初始面设为0最省事拓扑荷l1, 3, 5l的正负决定旋向径向指数p0, 1, 2增加环数上表里“仿真建议”列的 w0 越小空心环峰值位置就越靠近中心但如果 w0 小于几个网格间距中心会因采样不足而出现菱形畸变。因此模拟时要保证 $w_0/dx \geq 20$。如果 w01mmN512L6mm则 dx11.7μmw0/dx≈85足够。4. 验证LG01的核心特征干涉条纹、轨道角动量与模式分析4.1 平面波干涉用叉形条纹确认拓扑荷生成LG01后最直接的验证是让它与沿x方向倾斜的平面波干涉。干涉场强度为 $|E A \exp(i k x \sin\theta)|^2$其中A是参考光振幅θ是倾斜角。干涉图在涡旋中心会出现典型的叉形条纹——条纹分叉数等于拓扑荷的绝对值。下面的代码演示了这一过程% 参考平面波沿x方向传播并带一个横向相位梯度 theta 0.5e-3; % 倾斜角单位 rad E_ref exp(1i * k * X * theta); % 相当于叠加横向波矢 k*theta E_int E / max(abs(E(:))) 0.8 * E_ref; % 控制相对振幅 I_int abs(E_int).^2; figure; imagesc(x*1e3, y*1e3, I_int); axis image; colormap(gray); colorbar; xlabel(x (mm)); ylabel(y (mm)); title(Interference of LG01 with a tilted plane wave);这里的k*X*theta容易理解倾斜角 θ 使平面波具有横向相位变化。条纹间距为 λ/θ当 θ0.5mrad、λ633nm 时约 1.27mm在 6mm 网格内能容纳约4条亮纹足够看到叉形。实际调节 θ 时注意不要让条纹间距小于3个像素否则叉形看不清。叉形中间的暗线分叉方向取决于 l 的符号——这比强度图更可靠。4.2 数值计算轨道角动量密度拉盖尔-高斯模式的一个可测量特征是每光子携带 $l\hbar$ 的角动量。在数值上可以直接计算相位矩阵在圆周上的差分和。因为 l1 时绕中心一圈的总相位变化应为 $2\pi$。用MATLAB可以取一个半径 r0 的圆环上的相位值unwrap后做线性拟合获得斜率。方法对 l 较大时依然准确但要求环半径不能太小否则中心奇点附近的相位噪声会干扰。下面给出一个简单的相位斜率检验代码% 取半径接近强度峰值位置的圆环 r0 w0 / sqrt(2); % 峰值半径处 t linspace(0, 2*pi, 720); % 角度向量 cx 0; % 光束中心假设在原点 pts [cx r0*cos(t); cx r0*sin(t)]; % 2x720 % 双线性插值提取相位 phase_ring interp2(X, Y, phase_E, pts(1,:), pts(2,:), cubic); % 去除NaN中心奇点附近没有有效相位 phase_ring unwrap(phase_ring); % 线性拟合相位 slope * t offset coeff polyfit(t, phase_ring, 1); slope coeff(1); fprintf(Phase slope around ring: %.4f (should be %.4f)\n, slope, l);unwrap会自动补偿跳变得到的斜率等于 l。这里用interp2插值而不要直接取Phit的索引因为圆环不一定正好落在网格点上。slope应接近 1.0误差小于0.1可以认为模式正确。峰值半径 r0w0/√2 是LG01强度最大处信噪比最高。4.3 扫描拓扑荷l与径向指数p观察模式演变进一步可以让 l 或 p 在一个循环中变化验证强度环数量和相位奇点数的关系。一个典型脚本figure(Color,white); for l_val [1 2 3] % 将第3章的核心代码封装成函数lgBeam(p,l,w0,lambda,z,N,L) E_l lgBeam(0, l_val, w0, lambda, 0, N, L); subplot(2,3,l_val); imagesc(x*1e3, y*1e3, abs(E_l).^2); axis image; colormap(hot); title(sprintf(l %d, p 0, l_val)); end这里省略生成函数实际可以把第3章代码封装成lgBeam(p,l,w0,lambda,z,N,L)。l 绝对值越大空心环半径越大相位奇点强度保持不变但仍为零。p 增大时会出现 p1 个同心环并且环与环之间有相位节点径向相位反转。下表总结了几种常见 $LG_{p,l}$ 的可见特征模式pl强度环数相位奇点LG00001实心无LG01011空心1拓扑荷1LG02021空心半径更大1拓扑荷2LG11112同心环1拓扑荷1扫描 l 可以帮助你直观理解拓扑荷对环半径的影响环半径正比于 $\sqrt{|l|} \cdot w_0/\sqrt{2}$大约如此。对于 l3峰值位置约在 $\sqrt{3}\times 0.707 w_0 \approx 1.22 w_0$。这个关系在实验上往往被用来从测量强度环半径推算拓扑荷。5. 进阶LG01的传播仿真与模式叠加5.1 用FFT角谱法传播LG01光束实际仿真中除了观察束腰处的光斑你往往还需要知道 LG01 在传播几个瑞利距离后的样子。常见做法是用角谱衍射传播。在 MATLAB 中实现% 角谱传递函数 dx L / N; fx (-N/2 : N/2-1) / L; [FX, FY] meshgrid(fx, fx); % 角谱传递函数自由空间 H exp(1i * k * z * sqrt(1 - (lambda*FX).^2 - (lambda*FY).^2)); % 初始场E前面生成的LG01 E0 E; Ef ifft2(ifftshift(fftshift(fft2(E0)) .* H)); % 显示 I_prop abs(Ef).^2;这里有一个很关键的细节fft2的输出频率中心在数组第一个元素需要用fftshift移动到中心而ifft2之前再用ifftshift恢复。如果少一步光斑会偏移半格。这里的H在sqrt内可能为负倏逝波一般可以忽略也可以将负值置零。在自由空间传播、距离较小时没有大问题。传播距离 z 的单位与网格尺度一致。若 z 太大角谱中的高频相位变化过快会出现混叠。经验法则是保证 z 远小于 $N dx^2 / \lambda$。例如 N512, dx11.7μm, λ633nm则 $N dx^2 / \lambda \approx 0.111m$即传播距离超过 11cm 时就需要加大网格或增大 dx。你在模拟一个几十厘米甚至1米后的LG01时要重新设计网格尺寸。5.2 叠加LG01与LG00生成横向花瓣结构拉盖尔-高斯模式之间是正交的但叠加后会因多普勒相位差产生干涉图案。例如将 LG01 与 LG00 同轴叠加会得到一个沿方位角变化的光强分布花瓣数等于两个模式拓扑荷之差。实现起来只需要E_lg00 lgBeam(0, 0, w0, lambda, 0, N, L); % 假定函数已封装 E_sum E / max(abs(E(:))) 1.0 * E_lg00 / max(abs(E_lg00(:))); I_sum abs(E_sum).^2; imshow(I_sum, []);这里的相对振幅比会影响对比度。当两个模式幅度相等时花瓣达到最大对比度如果想模拟实验上常见的“不完全对准”可以给其中一个模式加很小的横向位移花瓣会变得左右不对称。这种叠加常被用于产生径向偏振类似结构的矢量光束或者用于测量模式之间的Gouy相位差。注意由于 LG00 的中心实心叠加图不再是空心而是被切割成两瓣或更多瓣。l 差为1时是两瓣像“花生”。这在LG01作为纠缠光子源调试中常被用来判断空间模式是否匹配。5.3 与实验对照相机积分效应与像素化实验上的CCD或CMOS相机记录的并不是精确的光强点而是每个像素上的积分强度。模拟中如果直接显示abs(E).^2会忽略像素的积分效应。要更接近实验结果可以用conv2对强度平面做一次“像素平均”。比如让模拟的每个网格对应一个像素可先计算局部均值I_sim abs(E).^2; fs 4; % 假设每个物理像素对应4×4个网格 kernel ones(fs, fs) / fs^2; I_camera conv2(I_sim, kernel, same);在实际对照实验时还需要考虑相机噪声、位深量化等。最省事的方法是先用imresize把强度图降到相机分辨率再显示。这样做不会改变结论但能让你在用MATLAB写报告或论文插图时视觉上更接近拍摄照片。6. 把LG01仿真写扎实的3个实用检查6.1 检查数值归一化与总能量守恒每次都验证一下总能量是否恒定。对初始场和传播后的场分别做sum(abs(E(:)).^2) * dx^2在角谱传播中能量应几乎不变。如果能量偏离超过1%多半是网格尺寸不够或传递函数在边缘被截断。归一化到峰值为1很容易但能量守恒更关键。6.2 避开中心奇点的采样错误LG01中心相位奇异相位上的值是未定义的。在angle计算后中心会出现随机数这是正常的。如果强行用unwrap从中心开始解卷会产生大幅振荡。避免方法是先遮蔽半径小于 0.1w0 的区域再求相位。而插值环提取相位时选择的半径不应小于峰值半径否则插值到奇点附近会被NaN污染。6.3 用单精度还是双精度生成初始场时用 double 即可但对于超大规模网格N2048以上或者在循环里做上千次传播可以考虑把场转成 single 以节省一半内存。注意fft2支持 single 输入但传递函数 H 里同时包含 k 和 z乘积精度在传播距离很大时可能不足所以传播时建议仍用 double。基本的原则是内存可控时全程 double兼顾速度和精度只有反复调参的大批量扫描时才用 single。本文还有配套的精品资源点击获取