涡旋光束的夫琅禾费衍射仿真（含部分源代码）

北京大学微电子学与固体电子学硕士

涡旋光束

涡旋是自然界中普遍存在的现象，如海洋环流、台风以及光涡旋。光学奇点是光场空间中具有不确定值或突变的某个点，光学涡旋同样具有这样的相位奇点。作为光学涡旋的代表，涡旋光束（Vortex beam)是研究光学中涡旋特性的有效工具。涡旋光束不仅具有由于螺旋相位结构而产生的轨道角动量，而且具有与偏振相关的自旋角动量，其螺旋相位因子可写成 exp(il\theta) ,这种涡旋光束可承载 lh/2\pi 的轨道角动量 ^[1] ，涡旋光束的光场表达式为

E(r,\varphi,z)=E_{0}(r,z)exp(-ikz)exp(il\theta),(1)

实验中常见的涡旋光束包括拉盖尔-高斯光束（Laguerre-Gaussian beams，LG beams），贝塞尔光束（Bessel beams）以及完美涡旋光束，此外，还有一些携带轨道角动量的涡旋光束如Mathieu光束，Ince-Gaussian光束以及矢量涡旋光束等。这里，我们以拉盖尔-高斯光束为例，对涡旋光束的夫琅禾费（夫琅和费）衍射进行仿真。拉盖尔-高斯光束的光场表达式为

E(x_{0},y_{0})=E_{0}\left( \frac{\sqrt{x_{0}^{2}+y_{0}^{2}}}{w_{0}} \right)^{\left| l \right|}exp\left[ \frac{-\left( x_{0}^{2}+y_{0}^{2} \right)}{w_{0}^{2}} \right]exp(il\theta),(2)

其中， E_{0} 为光场振幅， w_{0} 为光束的束腰半径， l 表示光束的拓扑荷数。这里，我们展示了 l=3\sim5 的拉盖尔-高斯光束的光强和相位分布如图1所示

夫琅禾费衍射公式为

U(x,y)=\frac{1}{iz\lambda}exp[ik\left( z+\frac{x^{2}+y^{2}}{2z} \right)]\int\int U_{0}\left( x_{1}+y_{1} \right)exp\left[ ik\frac{-\left( xx_{1}+yy_{1} \right)}{z} \right]dx_{1}dy_{1}，(3)

这里，我们主要仿真矩孔衍射（也称为方形孔径光阑衍射），弓形孔衍射（也称为弓形孔径光阑衍射），阵列光阑衍射。

涡旋光束的矩孔衍射

矩孔光阑的透过率函数可表示如下

rect\left( \frac{x-x_{0}}{a} \right)rect\left( \frac{y-y_{0}}{b} \right)=1,\left| \frac{x-x_{0}}{a}\right|\leq1/2,\left| \frac{y-y_{0}}{b}\right|\leq1/2，\left( 4 \right)

其中 a 为 x 方向上的边长， b 为 y 方向上的边长， \left( x_{0},y_{0} \right) 是矩阵光阑的中心。对正方形光阑来说， a=b 。这里，正方形光阑的示意图如图2所示。

由式子（2），我们不难发现其积分项对光束衍射特性的影响占主导作用。基于式子（2-4），这里仿真了不同拓扑荷数的拉盖尔-高斯光束的矩孔衍射如图3所示，取值： \lambda=632.8 nm , a=b=0.2mm , z=100mm , w_{0}=0.5mm 。可以看到随着涡旋光束拓扑荷数的增加，衍射亮斑的数目也在增加，而且具有旋转拖尾的特征。

涡旋光束的弓形孔衍射

我们知道圆形孔径光阑可定义如下

circ(r)=circ(\sqrt{x^{2}+y^{2}})=1,r=\sqrt{x^{2}+y^{2}}\leq R

这里，弓形孔的透过率函数可以看作圆孔光阑与矩孔光阑的差值，弓形孔径光阑的示意图如图4所示。

基于式子（2-4），不同拓扑荷数的拉盖尔-高斯光束的弓形孔衍射仿真如图5所示。取值： \lambda=632.8nm , w_{0}=0.5mm , a=0.6mm , z=200mm ， R=\sqrt{2}a/2 。可以看到随着拓扑荷数的增大，衍射级数增多，受到光阑影响，光强分布基本趋于圆形分布，而且同样具有旋转特性。

涡旋光束的阵列光阑衍射

这里我们考虑一个由圆孔构成的六边形阵列光阑，如图6所示。

不同拓扑荷数的拉盖尔-高斯光束的弓形孔衍射仿真如图7所示。取值： z=600nm ， R=0.05mm ，相邻圆孔圆心间距 d=0.2mm ，其它取值如弓形孔那节取值。可以看到采用六边形阵列光阑，衍射亮斑数量与拓扑荷数之间存在复杂关系。

部分Matlab源代码

clc
clear
close all
%------仿真变量------%
N = 200;            %取点数
lambda = 632.8e-6;  %波长 [nm->mm]
k = 2*pi/lambda;    %波数
w0 = 0.5;           %束腰半径 [mm]
a = 0.2;            %方形孔/矩孔边长 [mm]
b = a;              %方形孔/矩孔边长 [mm]
E0 = 1;             %光场振幅
l = 3;              %拓扑荷数
z = 100;            %传播距离 [mm]
%------程序变量------%
x = linspace(-2,2,N); %直角坐标
y = linspace(-2,2,N);
[x,y]=meshgrid(x,y);
x0 = linspace(-2,2,N); %直角坐标
y0 = linspace(-2,2,N);
[x0,y0]=meshgrid(x0,y0);
dx = x0(1,2)-x0(1,1);
dy = y0(2)-y0(1);
r = sqrt(x0.^2+y0.^2);  %柱坐标
phi = cart2pol(x0,y0);
%------初始光场------%
E1 = E0.*(sqrt(x0.^2+y0.^2)./w0).^abs(l).*exp(-(x0.^2+y0.^2)./(w0.^2)).*exp(1i.*l.*phi);