简介:本资源是一套面向光学工程初学者与科研人员的MATLAB实践工具包,聚焦涡旋光束生成原理与超表面相位调控建模,解决轨道角动量光场仿真、螺旋相位设计及纳米结构光学响应模拟等核心问题。压缩包共10个文件(9个.m源码+1个.docx说明文档),总大小仅18KB,轻量紧凑;其中.m文件涵盖涡旋相位分布构建(含不同拓扑荷l设置)、傅里叶域传播模拟、超表面等效相位调制模块及菲涅耳波带片辅助验证等功能,.docx文档则系统梳理了物理模型、代码逻辑与关键参数含义。已有880人学习下载,适合在无硬件实验条件时快速开展原理验证与算法调试。用户可直接运行各脚本观察涡旋光强度/相位图、理解l值对光束暗核与绕数的影响,并基于代码框架拓展超表面单元设计或与实际器件响应数据对接,是连接理论光学与计算仿真的实用入门载体。 做超表面的同学应该都经历过这么一个场景:课题任务书里写着“设计并仿真一个涡旋光束生成超表面”,配图是一张漂亮的甜甜圈光强图。你兴冲冲打开MATLAB,却发现不知道该从哪一行代码开始写。是先算相位?还是先设计单元结构?相位分布和超表面单元的旋转角到底怎么对应?远场仿真怎么看结果?这些环节看起来各自独立,实际是一条完整的技术链路,任何一个环节断裂,后面的仿真结论都站不住。
标题里的“B-涡旋光matlab_phasechange_涡旋光_涡旋光matlab_超表面_matlab超表面”,我理解核心任务是在MATLAB环境下完成涡旋光的相位变化(phase change)分析与超表面器件的仿真验证。这类工作常见于OAM光通信、微粒操控、超分辨成像等领域,近几年论文发得非常多。这篇文章就把我平时做这套东西的完整流程摊开讲:从涡旋光相位怎么生成,到相位怎么离散化、怎么映射到超表面单元,再到远场怎么验证结果、常见的坑有哪些。目标是让拿到类似题目的同学能直接照着跑出一版能用的仿真。
1. 涡旋光的相位内核:为什么相位变化是它的灵魂
1.1 从相位因子exp(ilθ)理解涡旋光的“身份标签”
涡旋光和普通平面波、高斯光束最大的区别,不在振幅,而在相位。它的电场表达式可以写成:
E(r, θ, z) = E0(r, z) · exp(ilθ)
这里面的 θ 是方位角,l 是一个整数,叫做拓扑荷(topological charge),这个 l 就直接决定了涡旋光的“身份标签”。l=0 的时候,相位项等于1,光场退化成普通的高斯光束;l≠0 的时候,波前会绕着传播轴扭转,形成螺旋状的相位面。
我在给课题组新生讲这个的时候,常用一个生活化的类比:普通光束的波前像一张平整的纸,你从侧面看它,所有位置的相位都一样;涡旋光的波前像是把这张纸沿着中心线剪开,然后其中一边往上拧、另一边往下拧,拧了 l 圈,最后首尾相接的地方就会出现一个不连续点。这个不连续点就是相位奇点(phase singularity),对应光强分布中一个中心暗核——理论计算中,奇点处的光强严格为零。
这个相位奇点极其关键。它意味着涡旋光在中心位置的光强是零,于是整个光束的强度分布呈现“甜甜圈”形,而且甜甜圈半径会随 |l| 增大而变大。做超表面设计时,你所有关于“怎么产生涡旋光”的努力,本质上都是在空间上构造这个螺旋相位分布,让光通过超表面之后,每个位置都带上对应的相位偏移。
1.2 相位变化(phase change)与相位差:调控的底层逻辑
标题里的“phasechange”,在超表面仿真的语境下有两层含义。一层是物理上的“相位变化量”——光经过超表面某个位置时,透射或反射光的相位相对入射光发生了多少改变;另一层是设计上的“相位分布”——为了生成目标光场,超表面上每个位置需要提供多少相位增量。
这两层是因果关系:你先明确目标光场需要的相位分布(也就是第二层),再通过结构参数去实现每个位置的相位变化量(也就是第一层)。所以整个MATLAB仿真的第一步,不是急着去画超表面单元,而是先把目标相位分布算出来。
涡旋光的相位变化量,数学上是一个绕拓扑荷 l 的线性角向变化。它的梯度方向沿角向切向,大小等于 l/r,注意这个大小随半径增大而减小。这意味着靠近中心的相位变化非常剧烈,中心点本身是奇点,相位不确定。做数值仿真的时候,这个特征会引起一系列问题,我在第4章会专门讲。
2. MATLAB生成螺旋相位:从网格坐标到相位图
2.1 第一步:构建坐标网格
MATLAB里做这类光学仿真的第一步都是建立网格。网格密度直接决定仿真精度,网格太少,相位奇点附近的相位变化会失真;网格太多,后面的FFT传播计算会非常慢。我一般用一个经验值:模拟区域边长取 100μm 左右,网格数取 512×512,对于可见光波段的超表面仿真,这个分辨率可以兼顾精度和速度。
主代码就这几行:
% 涡旋光螺旋相位生成 clear; close all; % 参数定义 L = 100e-6; % 物理区域边长 100 um N = 512; % 网格数 512x512 lambda = 633e-9; % 设计波长 633 nm l = 1; % 拓扑荷,可改为 2, 3, -1 等 % 坐标网格 x = linspace(-L/2, L/2, N); y = linspace(-L/2, L/2, N); [X, Y] = meshgrid(x, y); % 直角坐标转极坐标 [Theta, R] = cart2pol(X, Y); % 螺旋相位 phase_vortex = l * Theta;这里用meshgrid生成二维网格,cart2pol把直角坐标转成极坐标,得到方位角矩阵 Theta 和径向距离矩阵 R。phase_vortex = l * Theta就是涡旋光的相位分布核心表达式,它让每个点的相位按照方位角线性变化,绕一圈相位累计变化 2πl。
2.2 相位卷绕(phase wrapping)的处理
直接运行上面的代码,然后imagesc(phase_vortex)看相位图,你会发现一个很有意思的现象:相位图是一圈一圈的条纹,颜色从蓝到红突然跳变。这个跳变不是错误,而是相位周期性造成的。相位本身以 2π 为周期,计算机会把 -π 和 +π 判定为同一个相位值,但在分布图上,它们会显示成急剧的颜色跳跃。
处理方式有两种。一种是保留这种“包裹相位”(wrapped phase),画图时用phase_vortex = mod(phase_vortex + pi, 2*pi) - pi或者直接mod(phase_vortex, 2*pi)把它限制在固定区间内;另一种是用unwrap解包裹。对于涡旋光仿真,我强烈建议你用包裹相位做后续计算,不要轻易解包裹——因为涡旋光中心是真实的相位奇点,解包裹算法会沿着奇点产生一条人为的相位断裂线,反而把问题搞复杂。
% 相位包裹到 [0, 2*pi) 区间,便于可视化 phase_wrapped = mod(phase_vortex, 2*pi); figure; imagesc(x*1e6, y*1e6, phase_wrapped); axis xy; axis equal; colormap(hsv); colorbar; xlabel('x (um)'); ylabel('y (um)'); title(['拓扑荷 l = ', num2str(l), ' 的螺旋相位分布']);用hsv色彩映射画相位图最合适,因为它首尾相连,颜色环的跳变位置正好对应相位的 2π 周期边界。你会看到,l=1 的时候相位图有一圈渐变色带,l=2 的时候有两圈,非常直观。
2.3 复振幅构造与振幅调制
相位只是半边信息,实际仿真还需要构造完整的复振幅场。无振幅调制的情况下,涡旋光的复振幅可以写成:
E_vortex = exp(1i * phase_vortex);如果你想要更接近实际激光器输出的拉盖尔-高斯模,可以加上径向振幅项:
w0 = 20e-6; % 束腰半径 amplitude = (sqrt(2)*R/w0).^abs(l) .* exp(-R.^2/w0^2); E_vortex_lg = amplitude .* exp(1i * phase_vortex);很多初学者容易忽略振幅调制,直接拿纯相位分布去做远场仿真,结果也能看到甜甜圈,但暗核周围的次级环结构会偏亮,和真实实验对不上。建议一开始就把振幅项加上,仿真结果更接近实验。
3. 从连续相位到超表面单元:几何相位的映射方案
3.1 为什么超表面实现的是“离散相位”
理论上的涡旋光相位是连续变化的,但真实的超表面单元是有限的亚波长结构,你不可能让每个位置都完美实现任意相位的相位变化。通常的做法是把 0~2π 的连续相位范围量化成 N 个台阶,每个台阶对应一组结构参数。
这个离散化过程和把一个连续的圆弧用多边形逼近是一个道理。台阶数越多,逼近越精确,但同样面积的超表面上就需要更多种不同尺寸的单元,设计加工难度随之上升。实际设计中,N=4、8、16 是比较常见的取值,我以前用 8 台阶设计过一个涡旋光生成器,远场中心暗核已经很干净了。
% 相位离散化:将连续相位量化为 Nstep 个台阶 Nstep = 8; phase_quantized = round(phase_vortex / (2*pi/Nstep)) * (2*pi/Nstep); figure; imagesc(x*1e6, y*1e6, mod(phase_quantized, 2*pi)); axis xy; axis equal; colormap(hsv); colorbar; title([num2str(Nstep), '台阶离散化的涡旋光相位']);3.2 两种主流调控机制:传播相位与几何相位
把离散相位落实到具体超表面单元上,有两条路线。
第一条是传播相位(propagation phase)。通过改变纳米柱的尺寸(直径或边长),让x偏振和y偏振的光在柱内传播时积累不同的相位延迟。这种方法的优势是偏振无关性好、工作带宽较宽,缺点是每个离散相位台阶都需要对应一组独立的柱尺寸,需要做大量的参数扫描,设计周期长。
第二条是几何相位(geometric phase),也叫PB相位(Pancharatnam-Berry phase)。这里有一个极其优美的物理关系:各向异性的纳米柱绕自身中心旋转一个角度 α 之后,入射左旋圆偏振光(LCP)透过它后会变成右旋圆偏振光(RCP),同时携带一个大小为 2α 的附加相位。你不需要改变柱的任何尺寸,只需要旋转柱的方向,就能实现相位变化,而且这个相位在理论上是无损耗的。
我个人的经验是,对于涡旋光生成这类只需要相位调控、不需要偏振转换的器件,PB相位是性价比最高的方案。因为它的相位-转角关系是精确线性的,不需要查表,不依赖结构扫描精度,MATLAB代码写起来也简洁得多。
3.3 PB相位超表面:旋转角分布的MATLAB实现
用PB相位实现涡旋光,核心关系就一句话:要让透射光在某一点获得相位 φ(x,y),需要把该点的纳米柱旋转角度 α = φ(x,y)/2。
把前面生成的离散相位映射成转角矩阵,代码非常简洁:
% PB相位超表面:单元旋转角 alpha = phase / 2 alpha_map = mod(phase_quantized / 2, pi); figure; imagesc(x*1e6, y*1e6, alpha_map/pi*180); axis xy; axis equal; colormap(parula); colorbar; title('超表面单元旋转角分布 (度)');这一张图就是后续画版图的直接依据。你只要在GDS文件里,把每个单位格的纳米椭圆柱旋转对应的角度,整个超表面就设计完成了。我在实际项目中,会再把这个旋转角分布导出成合适的坐标格式,然后批量生成GDS版图。MATLAB的gdsii工具箱可以干这件事,如果你的课题还涉及加工流片,这一步会帮你省去大量手工画图的时间。
3.4 单元结构选型与参数扫描说明
PB相位方案中的单元结构,推荐用椭圆纳米柱或者矩形纳米柱。因为各向异性是产生几何相位的前提,圆对称结构旋转了跟没旋转一样,产生不了相位差。
单元周期 p 建议取设计波长的 1/2 到 1/3,满足亚波长条件,避免高阶衍射。比如 633nm 波长,周期取 300nm 左右比较合适。柱的长轴直径 a 和短轴直径 b 则需要保证足够大的双折射效应,一般做参数扫描来确定。扫描的时候用电磁仿真软件(如FDTD或COMSOL)计算不同长宽比下的透射效率,目标是找到一组 a、b,使得正交偏振光之间的相位差最大、同时透射率尽量高。
这一环节虽然不在MATLAB里完成,但MATLAB可以帮你处理电磁仿真软件导出的数据,比如读取扫描结果、挑选最优参数组合、把离散相位表拟合成结构参数表。做好这个衔接,整个设计流程会流畅很多。
4. 远场仿真验证:怎么看懂你的涡旋光
4.1 夫琅禾费衍射与FFT计算远场
超表面设计完,不能光看相位图就说“设计好了”,必须验证它真的能产生涡旋光。验证方法分近场和远场两种。近场可以直接看透射光的振幅和相位分布,但这只能说明入射光被调制了,不能证明结果是涡旋光。真正有说服力的验证是远场光强分布——看看光束传一段距离之后,能不能形成一个中心暗核的甜甜圈。
远场计算最常用的方法是夫琅禾费衍射。道理不复杂:入射光经过超表面后形成近场分布 E(x, y),传播到远场观察面上的复振幅分布,等于 E 的傅里叶变换(忽略常系数)。MATLAB里用fft2就能算出远场分布。
% 近场复振幅,取涡旋相位对应的相位图 E_near = exp(1i * phase_vortex); % 简化为纯相位调制 % 远场:傅里叶变换计算夫琅禾费衍射场 E_far = fftshift(fft2(fftshift(E_near))); I_far = abs(E_far).^2; % 由于超表面是周期结构,通常还需要乘上一个单元因子。 % 若只有一个单元周期,则直接看整体远场。 % 画远场光强分布 figure; imagesc(I_far); colormap(hot); axis image; axis off; title('远场光强分布 (甜甜圈)');代码里做了两次fftshift,这个细节很关键。fft2的输出低频分量在四个角上,不搬移的话,甜甜圈会显示成四个半圆拼起来的样子,看着非常困惑。用fftshift把零频移到中心,才能看到完整的环形分布。
4.2 甜甜圈、暗核与拓扑荷的判断
看到环形光强图只是第一步,你还得验证拓扑荷对不对。根据涡旋光的性质,中心暗核的半径随 |l| 增大而增大,所以单靠一张远场图,你可以粗略判断拓扑荷的大小,但无法精确判断。
更可靠的验证方法是干涉。把涡旋光和同一波长的高斯平面波叠加,干涉条纹会出现旋臂状分叉结构,分叉的数目恰好等于 |l|。MATLAB里实现干涉就是两束光的复振幅相加,代码也不长:
% 涡旋光与高斯平面波的干涉图 E_gaussian = ones(size(E_near)); % 理想平面波近似 E_interference = E_near + E_gaussian; I_interference = abs(E_interference).^2; figure; imagesc(I_interference); colormap(gray); axis image; axis off; title(['涡旋光与平面波干涉图, l = ', num2str(l)]);我处理过的很多案例里,拓扑荷方向也很容易出错。l 的符号决定了涡旋光的旋向——正值是右旋(沿传播方向看相位顺时针增加),负值是左旋。MATLAB的atan2生成的方位角 θ 是逆时针增加的,所以直接l * Theta得到的相位对正 l 是逆时针旋转。如果你设计的超表面加工出来旋向不对,检查这个符号问题,通常就能找到原因。
4.3 仿真中的边界效应与观察窗选择
还有一个很多人踩过的坑:FFT远场结果中,中心暗核附近会出现规则排列的条纹或亮斑,这不是涡旋光本身的特征,而是有限孔径的边界效应。因为你的超表面尺寸是有限的,矩形孔径的傅里叶变换会叠加上一个二维sinc函数旁瓣,导致远场中出现十字状闪亮条纹。
减少边界效应有三个办法:
- 给近场振幅加一个软边界的超高斯窗函数,让振幅在边界处平滑衰减到零;
- 让模拟区域远大于超表面区域,四周留一圈“暗区”;
- 观察远场时,只取中央的一小块区域分析,避开旁瓣。
我实际做仿真的时候,习惯在超表面外围加一圈2~3个单元格不产生涡旋相位的“边框”,这样既能抑制边界效应,也能让器件在实验里更好对准。
5. 相位离散化的工程代价:量化误差与杂散衍射
5.1 量化台阶数对远场纯度的直接影响
连续相位被量化成 N 台阶之后,远场中必然会出现额外的杂散衍射级。道理可以这么理解:量化后的相位分布相当于“理想涡旋相位”加上一个“量化误差相位”,这个误差是周期性的,会在远场产生高阶衍射峰。
我做过一个量化台阶数的参数扫描,结论很明确:台阶数少于4时,远场会出现明显的杂散亮点,甜甜圈形状被污染;达到8台阶之后,中心暗核和环状结构已经非常干净;16台阶以上,肉眼很难再看出差别。所以不是台阶数越多越好,够用就行。8台阶是精度和加工成本之间的一个均衡点,我的项目里基本都是这个值。
5.2 中心区域相位奇点的特别处理
上一章我提到涡旋光中心有一个相位奇点,这里有个非常隐蔽的问题:量化之后,奇点附近的相位台阶会出现“螺旋阶梯”结构,这些台阶是正常的,但相邻台阶之间的相位不连续会产生强烈的相位梯度。实际加工出来的结构,不可能真的实现无穷大的相位梯度,这部分区域的纳米柱通常是圆形的、不具备各向异性,导致奇点附近的衍射效率降低,暗核边缘出现亮环。
解决这个问题的方法是把奇点附近一定半径内的中心区域切掉(振幅设为零),或者用更细的网格局部加密。这个半径一般取 2~3 个单元格间距就够。中心挖掉的区域虽然损失了一部分能量,但换来的是远场暗核位置更干净,对于最终应用来说,这通常是划算的。
5.3 单元相位覆盖不全带来的设计偏差
实际设计超表面单元库时,你可能没法在0到2π范围内做到完全覆盖。比如用传播相位方案,材料吸收会导致某些尺寸的柱透射率太低;用PB相位方案虽然相位全覆盖,但交叉偏振转换效率会随波长漂移。
遇到这种情况,最直接的做法是把未覆盖的相位区间重新映射到最近的已有相位值上。但要注意,这个映射过程本身会再次引入系统误差,导致远场中心暗核不对称。所以在单元库设计阶段,我会多做一步:把期望相位值附近可用的结构参数全部列出来,从中选择透射率最高的一组,而不是只看相位接近的。相位差一点点没关系,效率损失是不可逆的。
6. 从单功能到复合功能:相位叠加与扩展思路
6.1 涡旋光+聚焦:单一超表面同时完成两项任务
超表面最大的优势之一,就是可以把多种功能需要的相位分布“叠加”在同一片表面上。只要每个功能对应的相位是独立的,总相位就是它们的和(或差)。
最常见的扩展是涡旋光生成加聚焦。让透射光在产生螺旋相位的同时,还汇聚到一个焦点上。聚焦项相位是双曲面的球面波相位:
% 涡旋 + 聚焦 复合相位 f = 50e-6; % 设计焦距 k = 2*pi/lambda; phase_focus = -k .* (sqrt(X.^2 + Y.^2 + f^2) - f); % 复合相位:叠加后包裹 phase_combo = mod(phase_vortex + phase_focus, 2*pi);这一步做完,你就能设计一个“在焦平面上产生涡旋光”的超表面,这也是很多光镊应用的基础配置。叠加的时候建议用严格的球面相位公式而不是傍轴近似,因为超表面的数值孔径可能很大,傍轴公式会引入明显像差。
6.2 多涡旋复用:阵列分光与拓扑荷叠加
再进一步,可以做一个涡旋光阵列。比如把一个超表面划分成多个子区域,每个子区域生成不同拓扑荷的涡旋光,远场就会出现多个暗核,每个暗核对应一个OAM通道。这在光通信的模分复用里非常有用。
具体做法有两种:
- 划分扇区,不同角度范围分配不同拓扑荷;
- 叠加不同方向的空间频率偏折相位,让不同拓扑荷的涡旋光向不同方向传播。
第二种做法在MATLAB里实现,就是在总相位里加一个光栅项:
% 向 x 方向偏折的光栅相位 kx = 2*pi / 10e-6; % 光栅周期 10 um phase_blazed = kx * X; % 比如 l=0 和 l=2 两个通道 phase_ch0 = mod(phase_blazed, 2*pi); phase_ch2 = mod(2*Theta + phase_blazed, 2*pi);远场仿真之后,你会看到两个不同方向的涡旋光斑。调整光栅周期可以控制通道间距,这在做OAM阵列器件时非常实用。
6.3 一点点个人建议
做了几个涡旋光超表面项目之后,我最深的感受是:这套仿真的技术壁垒并不高,真正的难点在于“确认每一步的结果都对”。很多人第一次跑出甜甜圈就以为自己成功了,实际上可能只是看到了FFT的环形旁瓣,或者相位奇点没落在网格上造成的伪暗核。
建议拿到题目后,先复现一篇论文里的经典结果。把相位生成、离散化、单元映射、远场仿真这条链路完整跑通一次,中间每一步都停下来对比一下已知结论,比如l=1的远场光强半径应该等于什么量级、干涉条纹应该有几条分叉,全部对上之后再开始做自己的设计。这个“慢就是快”的过程,是我能给出的最实在的建议。等你跑通了一遍,后面换波长、换结构、换拓扑荷,都只是改几个参数的事。
本文还有配套的精品资源,点击获取