光学实验室里折腾过激光器的人大概都绕不开两个名字Hermite-Gaussian光束HG和Laguerre-Gaussian光束LG。用Matlab把它们算出来、画出来不是玄学是很多课题的第一步——无论是设计光学谐振腔、做原子捕获还是搞轨道角动量复用这两个模式家族都会反复出现在你的代码里。HG和LG本质上都是傍轴波动方程在自由空间中的本征解区别只在于你选用什么坐标系去求解。矩形对称的谐振腔自然长出HG模柱对称的光纤或者具有旋转对称性的系统则更容易出现LG模。很多教材把它们列成两章看起来公式又多又长实际上在Matlab里实现并不复杂。只要抓住Hermite多项式、Laguerre多项式、高斯包络和相位因子这几个关键部件半小时就能写出一个能出图的仿真脚本。这篇博文我会按自己平时写代码的习惯来走一遍先讲清楚HG和LG的数学结构再给出可直接跑的Matlab代码然后是相位、强度、传播这些常踩的坑。内容面向光学方向的研究生和工程师也欢迎刚接触光束模拟的本科生直接抄作业。1. 动手之前先弄清楚HG和LG在数学上到底长什么样1.1 为什么这两个模式总是成对出现激光谐振腔里的稳定横模可以用一套正交完备的模式函数来展开。用直角坐标解出来的那套正交基就是Hermite-Gaussian模用柱坐标解出来的那套就是Laguerre-Gaussian模。HG的场分布沿x和y方向都是分离的所以它的空间图案是规则的矩形网格状光斑LG的场分布则带一个方位角相位因子强度呈环形中心往往有个暗点。两者的关系不是并列而是互转。一个HG模经过适当的柱透镜变换可以变成LG模这就是所谓的模式转换器原理。我最早写这个模拟脚本就是因为要研究这种转换前后的光场变化不能用实验反复调镜片只能先在Matlab里把两个模式的场分布都算出来。对于做光通信的人LG模更诱人——它的相位涡旋携带轨道角动量OAM不同拓扑荷对应不同复用信道。而对于做激光器设计的人HG模更实用因为矩形增益介质或半导体激光器里天然就容易出现低阶HG模。把两个模式放同一个代码框架里用起来非常顺。1.2 两个核心公式HG模在z0截面的标量场可以写成[ E_{n,m}(x,y)C_{nm}\cdot H_n\left(\frac{\sqrt{2}x}{w_0}\right)H_m\left(\frac{\sqrt{2}y}{w_0}\right)\exp\left(-\frac{x^2y^2}{w_0^2}\right) ]其中(n,m)是模式阶数(H_n)是n阶Hermite多项式(w_0)是束腰半径(C_{nm})是归一化常数。注意这里的Laguerre多项式用带参数形式场分布写作[ E_{p,l}(r,\phi)C_{pl}\left(\frac{\sqrt{2}r}{w_0}\right)^{|l|}L_p^{|l|}\left(\frac{2r^2}{w_0^2}\right)\exp\left(-\frac{r^2}{w_0^2}\right)\exp(-il\phi) ]这里的(p)是径向阶数(l)是方位角阶数也叫拓扑荷。(L_p^{|l|})是关联Laguerre多项式。HG里的(n,m)表示x方向、y方向各有多少个暗纹交点LG里的(p)表示径向暗环个数(l)则决定了中心相位奇点的拓扑荷大小和方向。所有变量里最容易写错的是Hermite多项式的自变量必须是(\sqrt{2}x/w_0)不是(x/w_0)也不是(2x/w_0)。如果你发现算出来的光斑零点位置和文献对不上九成是这里出了问题。2. 建模思路与参数选择这步决定了后续代码好不好改2.1 先想清楚你要模拟什么写仿真脚本最大的忌讳是一上来就写循环和画图。我的习惯是先问自己三个问题要算近场还是传播后的光场要看强度分布还是相位分布模式阶数是大是小如果只是验证HG和LG的经典图案直接在z0截面上计算就够不需要算传播。如果是研究光束在自由空间或透镜后的演化就得走角谱传播或菲涅尔衍射。如果模式阶数超过十几普通双精度计算也会遇到多项式数值溢出的问题那就要换更大范围的归一化或者干脆用符号计算先推导。在Matlab里我一般把网格、波长、束腰、模式阶数、传播距离全部放在脚本最前面方便反复调整。代码风格不需要多花哨变量名直观最重要。我自己的模板大致是这样的% 基础参数 lambda 1.064e-6; % 波长 1064nm w0 0.5e-3; % 束腰半径 0.5mm L 6e-3; % 空间网格范围 N 512; % 网格点数 z 0; % 传播距离0表示束腰处网格范围可能有人会写成10cm甚至1m这样其实非常浪费。光束在束腰附近的有效范围通常就是几个(w_0)。网格太小也不行高阶模的光斑扩展明显截断了会导致强度图案边缘畸变。N取512在大多数情况下都够既不会太慢也能让HG(2,1)的精细条纹看得清楚。2.2 这一步最容易被低估网格中心对齐用meshgrid生成坐标网格时如果不注意中心对齐光斑中心会偏出画面而且相位图会出现奇怪的斜条纹。我最常用的做法是x linspace(-L/2, L/2, N); [X, Y] meshgrid(x, x); R sqrt(X.^2 Y.^2); Phi atan2(Y, X);R和Phi几乎每个模式模拟都要用建议在初始化阶段就算好后面直接调用。需要留意的是meshgrid的两个输出X和Y是完整矩阵不是向量。第一次写的人容易用X.^2 Y.^2时忘记加号两边的点乘号结果得到的是一个标量或者报错矩阵维度不一致。为了后面复用我把主要参数整理成一张表写进注释里也方便读者对照自己的实验条件参数符号示例值说明波长(\lambda)1064 nm决定衍射尺度束腰半径(w_0)0.5 mm高斯光束最小光斑半径网格范围L6 mm一般取5~8倍(w_0)网格点数N512越大图案越平滑但计算越慢模式阶数n,m / p,l2,1 / 0,1根据研究目标设置3. Hermite-Gaussian光束的Matlab实现3.1 Hermite多项式不需要符号工具箱Matlab的符号数学工具箱里有现成的hermiteH函数但小规模网格上逐点调用符号函数速度让人头疼而且高维数组时内存开销很大。更实用的方法是用递推公式自己实现一个函数[ H_0(x)1,\quad H_1(x)2x,\quad H_{n1}(x)2xH_n(x)-2nH_{n-1}(x) ]写成Matlab函数就是function Hn hermite_poly(n, x) % 递推计算Hermite多项式 H_n(x) if n 0 Hn ones(size(x)); return; end if n 1 Hn 2 * x; return; end Hn_2 ones(size(x)); Hn_1 2 * x; for k 2:n Hn 2 * x .* Hn_1 - 2 * (k - 1) * Hn_2; Hn_2 Hn_1; Hn_1 Hn; end end用这个函数HG(n,m)的主代码可以写得很短n 2; m 1; Hx hermite_poly(n, sqrt(2) * X / w0); Hy hermite_poly(m, sqrt(2) * Y / w0); E_HG exp(-(X.^2 Y.^2) / w0^2) .* Hx .* Hy; I_HG abs(E_HG).^2; figure; pcolor(X * 1e3, Y * 1e3, I_HG); shading interp; colormap(hot); axis image; xlabel(x (mm)); ylabel(y (mm)); title(sprintf(HG(%d,%d) intensity, n, m));这段代码直接跑就能看到经典的瓣状图案。HG(2,1)的图案是6个亮斑x方向有两个暗纹将光斑分成三段y方向有一个暗纹分成两段组合出3乘2的亮斑网格。3.2 束腰、范围和归一化的关系如果不做定量分析只画强度分布归一化常数可以省。但如果你想计算能量占比、模式耦合系数或者把HG和LG放在同一个能量标尺下比较就必须带上归一化常数。完整的HG归一化常数是[ C_{nm}\sqrt{\frac{2}{\pi\cdot 2^{nm}\cdot n!\cdot m!}}\cdot \frac{1}{w_0} ]在Matlab里实现只需要一行Cnm sqrt(2 / (pi * 2^(nm) * factorial(n) * factorial(m))) / w0; E_HG Cnm * E_HG;记住这里有两个容易踩的点一是factorial在高阶数时会很快溢出二是归一化常数里的1/w0不能漏。如果你只做相对强度图漏掉系数无所谓但如果你把E_HG和另一个光束的场做干涉叠加系数不一致干涉条纹的对比度就是错的。我遇到过一位同学他把HG(5,5)的场算出来光斑中心数值大得离谱边缘又趋近于零。后来发现问题不是出在高斯包络而是Hermite多项式的高阶项已经远超双精度能表示的范围。对这种高阶情况最好改用归一化的Hermite函数或者把自变量缩放到(|x| \le w_0/\sqrt{2})附近再算。4. Laguerre-Gaussian光束的Matlab实现4.1 拉盖尔多项式同样用递推LG模的径向部分涉及关联Laguerre多项式递推关系比Hermite稍复杂但也很固定[ L_0^{m}(x)1,\quad L_1^{m}(x)m1-x ][ L_{k1}^{m}(x)\frac{(2k1m-x)L_k^{m}(x)-(km)L_{k-1}^{m}(x)}{k1} ]我写的是function Lp laguerre_poly(p, m, x) % 递推计算关联Laguerre多项式 L_p^m(x) if p 0 Lp ones(size(x)); return; end if p 1 Lp (m 1) - x; return; end L_km2 ones(size(x)); L_km1 (m 1) - x; for k 2:p Lp ((2*k - 1 m - x) .* L_km1 - (k - 1 m) .* L_km2) / k; L_km2 L_km1; L_km1 Lp; end endLG模的主代码比HG多一个方位角相位项p 0; l 1; rho2 2 * R.^2 / w0^2; Lpl laguerre_poly(p, abs(l), rho2); E_LG (sqrt(2) * R / w0).^abs(l) .* Lpl ... .* exp(-R.^2 / w0^2) ... .* exp(-1i * l * Phi); I_LG abs(E_LG).^2; figure; pcolor(X * 1e3, Y * 1e3, I_LG); shading interp; colormap(hot); axis image; xlabel(x (mm)); ylabel(y (mm)); title(sprintf(LG(%d,%d) intensity, p, l));这里abs(l)不要写错。对方位角阶数取绝对值是因为Laguerre多项式的上标必须是非负整数而相位因子里的l可以带正负号。正负拓扑荷的强度分布完全一样但相位涡旋方向相反。4.2 相位涡旋和轨道角动量如果是第一次看LG模式除了环形强度更值得关注的是相位分布。把相位画出来figure; pcolor(X * 1e3, Y * 1e3, angle(E_LG)); shading interp; colormap(hsv); axis image; xlabel(x (mm)); ylabel(y (mm)); title(sprintf(LG(%d,%d) phase, p, l));你会看到相位图从0到(2\pi)绕中心转了一圈中心处相位是不连续的。这就是相位涡旋对应每个光子携带(l\hbar)的轨道角动量。拓扑荷为1时相位旋转一圈拓扑荷为2时旋转两圈。强度中心始终是暗的因为相位奇点处场必须为零这也带来了LG模经典的甜甜圈形状。实验中大家经常用拓扑荷为1的LG模去做受激辐射损耗显微STED用不同拓扑荷的LG模做模式复用通信。这个模拟脚本虽然简单但它是理解这些应用的第一步。有了场分布后面算传播、算聚焦、算耦合都会顺很多。5. 从束腰截面到自由空间传播用角谱法把光束推出去5.1 角谱法的核心思想只算z0截面的场分布其实只是一个静态快照。很多场景下我们希望看到光束传播一段距离后的样子比如经过透镜聚焦、经过一段光纤或者单纯在空气中传播。自由空间传播最省事的办法是角谱法思路一句话就能说清楚把空间场做二维傅里叶变换得到不同方向平面波分量的振幅然后给每个分量乘一个传播相位因子再做逆傅里叶变换。角谱法在Matlab里只需要几行z_prop 0.02; % 传播距离 20mm dx L / N; fx (-N/2 : N/2-1) * (1/L); [FX, FY] meshgrid(fx, fx); % 角谱传递函数 H_AS exp(1i * 2 * pi * z_prop / lambda) ... .* exp(-1i * pi * lambda * z_prop * (FX.^2 FY.^2)); spectrum fftshift(fft2(E_HG)); spectrum_prop spectrum .* H_AS; E_prop ifft2(ifftshift(spectrum_prop)); figure; pcolor(X * 1e3, Y * 1e3, abs(E_prop).^2); shading interp; colormap(hot); axis image;这里的传递函数包含了两个指数项第一个是整体载波相位第二个是由横向空间频率决定的色散相位。角谱法的好处是单次快速傅里叶变换代价是网格采样率要满足(dx \lambda)否则频率混叠会把结果污染得很厉害。对可见光和近红外光这个条件在常规毫米级网格上很容易满足。5.2 传播模拟中最常出现的三个问题第一个问题是坐标轴方向。Matlab的fft2默认把零频放在矩阵左上角而meshgrid的FX、FY是零频在中心所以必须先做fftshift对齐。很多初学者漏掉这一步出来的光斑是斜的或者直接跑出画面。第二个问题是传播距离符号。角谱法里传播因子中的相位是(-i\pi\lambda z(f_x^2f_y^2))如果符号写反光束不但不发散反而会聚焦回去物理上完全讲不通。对比验证的办法很简单看传播到(z_R\pi w_0^2/\lambda)处的光斑宽度是否约为(\sqrt{2}w_0)。第三个问题是边缘截断。角谱法默认场在网格边缘外是周期的如果网格范围太小边缘的高斯尾巴被硬切掉逆变换后就会出现周期性的旁瓣条纹。我的经验是网格范围至少取6倍束腰高阶模可以放宽到8倍。6. 常见问题与排坑实录这些东西教材里一般不写6.1 一张表对照排查我把自己和周围人踩过的坑整理成一张速查表今后你写类似脚本时可以直接照着查现象可能原因解决方案光斑图案不对称或整体偏移网格中心不对齐meshgrid范围不是对称的用linspace(-L/2, L/2, N)HG零点位置和文献不符Hermite多项式自变量少乘了(\sqrt{2})检查sqrt(2)*X/w0LG中心不是暗斑方位角相位项写成了0或者l0检查exp(-1i*l*Phi)传播后图案出现周期性条纹网格范围太小导致边缘截断增大L至少6倍(w_0)传播结果异常发散角谱符号错误或fftshift顺序不对与理论w(z)对比验证高阶模数值溢出Hermite多项式递推溢出改用归一化多项式或缩小自变量范围相位图跳变无法辨认angle结果位于(-\pi)到(\pi)有跳变用unwrap或wrapToPi处理6.2 一个实用的验证技巧代码写完之后别急着画炫酷的彩色图先用一个最简单的物理场景验证逻辑对不对。我常用的验证方法是把HG(0,0)写跑通它应该退化成普通高斯光束把LG(0,0)写跑通它也应该退化成同一个高斯光束。如果这两个结果不一致说明你的常数或坐标定义有偏差。再验证LG的环形半径。对(p0, l1)的模式强度极大值出现在(r_{\max}w_0/\sqrt{2})附近。这个解析结果用来检查网格尺度和径向指数特别有效。我试过用0.5mm束腰算强度峰值正好落在0.353mm附近和解析值对上了基本可以放心继续算高阶模。还有一个画图层面的小技巧强度分布跨度可能很大直接pcolor会把暗纹理压掉。用imagesc(I)加colorbar或者对强度取对数坐标imagesc(10*log10(I))都能让纹理更清楚。尤其是高阶LG模径向明暗环之间的对比非常大线性色标几乎看不清暗环。6.3 从标量模式到混合模式扩展HG和LG不一定是孤立存在的。实际光路里经常出现的是几个模式的叠加比如一个HG(1,0)加上一个HG(0,1)相位关系不同合成光斑可能是一个斜向条纹也可能是一个旋转结构。在Matlab里做叠加非常方便E_mix E_HG E_LG; I_mix abs(E_mix).^2;叠加之前一定要保证两个场的网格、波长、束腰完全一致否则相位关系是乱的。叠加的结果对相对相位极其敏感稍微改一个符号中心暗点就可能变成亮斑。这就是为什么前面所有归一化问题都要提前理清楚。混合模式在实际中很重要比如柱透镜模式转换器产生的就是HG到LG的连续变换中间状态的每一帧都是一个混合模式。你还可以写一个循环脚本让两束分量之间的相位差从0扫到(\pi)把结果导出成GIF授课时放给学生看理解涡旋光束会直观很多。高阶扩展方面你还可以把Hermite和Laguerre多项式替换成Ince多项式得到Ince-Gaussian模式它在直角坐标和极坐标之间架了一座桥。不过那是另一个故事了先把手上的HG和LG跑稳后面的路自然就顺了。结尾这套代码在实际中还能怎么用我平时拿这套Matlab脚本做的事不只是画几张好看的图。研究生开题阶段用它预估实验光斑形状写论文时用它生成示意图做光纤耦合模拟时把输出场投影到本征模上算耦合效率。最近我在做一个涡旋光与原子相互作用的仿真也是从这段LG代码改出来的——把标量场换成多个带不同拓扑荷的场叠加配合原子能级方程就能初步验证轨道角动量的转移。如果你也想快速上手建议先把我上面给的HG和LG两段代码原样跑一遍然后把n,m,p,l这四个阶数来回改观察图案变化。弄清楚“零点和暗环怎么来的”比记忆任何公式都重要。后面无论换成矢量光束、部分相干光束还是脉冲光束底层框架都能复用。