1 引言

衍射是波普遍具有的性质。因为衍射效应的存在，所以即便不考虑成像系统的误差也不可能完全分辨极小物点。1987年，Durnin等[提出存在径向光强分布满足贝塞尔函数的无衍射光束，即贝塞尔光束。贝塞尔光束是亥姆霍兹方程的精确解，传播距离只影响其相位分布，而不影响横截面光强，是无衍射光束。贝塞尔光束的性质使得其可以用于超分辨成像、光子操作、微纳加工等领域。其等相面呈圆锥形，所以只要赋予平面波圆锥形的等相面即可产生贝塞尔光束。早期的研究中采用轴锥透镜来实现，除了轴锥透镜还可以用空间光调制器或者亚波长圆环[来产生锥形等相位面。但是传统方法面临着光学系统复杂、透镜加工困难等问题。近年来超材料和超表面的提出，为产生贝塞尔光束提供了新的思路。

超材料[是一种将特定形状的结构按照一定规律排列而成的三维材料。根据单元结构形状以及排列规律的不同可以得到不同的等效电磁特性(例如介电常数、磁导率)。使用超材料，人们实现了负折射[、电磁隐身[等传统材料难以实现的功能。但是超材料存在三维加工和金属损耗等问题，影响了其进一步的应用研究。超表面[。当圆偏振光入射到几何相位超表面，与结构相互作用后的透射光除了含有与入射光相同偏振态的分量外，还包含与入射光偏振态正交的光束，即交叉极化光束。交叉极化光相对入射光的相位突变与结构的指向角呈2倍关系[。

通常几何相位超表面采用离散结构来实现，离散结构的旋转角与几何相位有着一一对应关系。应用几何相位超表面可以设计涡旋光束产生器件[、全息器件[、贝塞尔光束产生器件[等有着复杂相位分布的二维光学器件。但这些器件将相位进行离散化，离散的相位会使得器件的性能降低，导致较窄的低频工作带宽等。2015年，中国科学院光电技术研究所罗先刚课题组提出了悬链线型超表面，可以实现连续的相位调控[。正如上文介绍的，悬链线狭缝由悬链线垂直平移构成，所以狭缝的两端将十分狭窄，不利于加工。本文提出用等宽悬链线狭缝构成的几何相位超表面来产生贝塞尔光束，相对于传统悬链线狭缝，等宽悬链线狭缝更易于加工。

在超表面的设计以及仿真中，通常用电磁仿真软件CST进行建模仿真和分析结果，并将数据保存，然后用Matlab绘制图形，分析效果。在设计V型天线、矩形孔径，以及悬链线狭缝等超表面结构时，可以直接在CST中进行建模。但是对于法向等宽悬链线以及其他不规则形状的结构时，CST建模较为复杂。另一方面，在涉及到需要对亚波长单元阵列的排列方式进行优化时，则需要多次在CST里重复建模，然后在Matlab里分析数据、优化结构，这将影响工作效率。在以往的研究中，有文章提出Matlab与CST联合仿真，但是因为软件版本、工作目的不同、以及CST繁多的底层命令而难以复现并应用到自己的工作之中。

本文在设计等宽悬链线及其阵列构成的超表面时，采用Matlab调用CST Microwave Studio进行建模，仿真。得到了可产生贝塞尔光束的等宽悬链线型器件。联合仿真的方法可用于设计更加复杂的结构，并结合Matlab中的优化算法实现一体化电磁结构仿真设计。

2 结构设计

2.1 等宽悬链线

两端固定的链在重力作用下形成的曲线是一种悬链线。1670年，约翰·雅可比、惠更斯等人研究了悬链线的数学性质。其表达式为

$

y = \frac{\varLambda }{\pi }{\rm{ln}}\left( {\left| {{\rm{sec}}\left( {\frac{{\pi x}}{A}} \right)} \right|} \right),

$

(1)

其中：Λ表示悬链线的水平跨度。由式可知，y在x=0.5Λ处为无穷大。所以在设计结构时，将x的取值范围限制在区间(-0.5Λ+dx，0.5Λ-dx)内。在悬链线型超表面中，因为其切向角与位置x成正比，从而可以实现线性的、连续的相位调控。如Λ+dx, 0.5Λ-dx)内，绘制N个悬链线上的点。对于其中一个点的坐标(xn, yn)，其切角为

图 1 悬链线狭缝结构示意图。(a)普通悬链线狭缝；(b)等宽悬链线狭缝 Fig. 1 Diagrams of catenary slits. (a) Ordinary catenary slit; (b) Equal-width catenary slit

$

\theta = {\rm{arctan}}\left( {\frac{{{y_{n + 1}} - {y_n}}}{{{x_{n + 1}} - {x_n}}}} \right){\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} ,

$

(2)

则与其法向垂直距离为w的点的坐标(xn', yn')为

$

\left\{ {\begin{array}{*{20}{l}}

{x_n^\prime = {x_n} - w \cdot {\rm{cos}}\theta }\\

{y_n^\prime = {y_n} + w \cdot {\rm{sin}}\theta }

\end{array}} \right.。$

(3)

绘制出的等宽悬链线如

2.2 基于等宽悬链线阵列的贝塞尔光束产生器

贝塞尔光束是亥姆霍兹方程的精确解，其表达式可以写作[

$

\begin{array}{*{20}{l}}

{E(x, y, z, t) = \frac{1}{{2\pi }}{\rm{exp}}[{\rm{i}}({k_z}z - \omega t)]}\\

{{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot \int_0^{2\pi } {{\rm{exp}}} [{\rm{i}}{k_{\rm{r}}}(x{\kern 1pt} {\rm{cos}}{\kern 1pt} {\kern 1pt} \varphi + y{\kern 1pt} {\rm{sin}}{\kern 1pt} \varphi )]{\rm{d}}\varphi , }

\end{array}

$

(4)

式中：ω是光的角频率，kr是横向波矢量分量，kz是波矢的纵向分量。由式(4)可知，传播距离z只影响贝塞尔光束的相位项，横截面光强在不同距离处保持不变，即为无衍射光束。上式的积分项表明，贝塞尔光束可以看作许多平面波的叠加，具有干涉场的性质。并且各平面子波与z轴有相同的夹角，即等相面为圆锥形。这就是平面波附加圆锥形相位可以得到贝塞尔光束的原因。

垂直于传播方向的贝塞尔光束电场分布表达式：

$

E(r, \varphi ) = {{\rm{J}}_l}{\rm{exp}}({\rm{i}}{k_{\rm{r}}} + {\rm{i}}l\varphi + {k_z}z),

$

(5)

式中：Jl表示第l阶贝塞尔方程，φ表示横截面的方向角。上式中指数项决定了贝塞尔光束的相位分布，对于零阶(l=0)贝塞尔光束，其相位分布如

图 2 等宽悬链线阵列结构。(a)零阶贝塞尔光束相位分布图；(b)等宽悬链线零阶贝塞尔光束产生器件 Fig. 2 Equal-width catenary array structure. (a) Phase distribution of zero-order Bessel beam generator; (b) Zero-order Bessel beam generator constructed by equal-width catenary arrays

$

\varPhi = - 2\sigma {\kern 1pt} {\rm{arctan}}\left( {\frac{{{\rm{d}}y}}{{{\rm{d}}x}}} \right) = - 2\sigma \frac{\pi }{\varLambda }x,

$

(6)

其中σ=1表示入射圆偏振光为右旋，σ=-1表示入射圆偏振光为左旋。arctan(dy/dx)表示悬链线切线和x轴的夹角。将悬链线按极坐标排列，可以实现贝塞尔光束半径方向和切向的相位延迟。悬链线端点轨迹为[

$

r = \frac{{(l - 2)\varphi + (2m + 1)\pi }}{{{k_{\rm{r}}}}},

$

(7)

其中m为正整数。使用上式，得到如

3 Matlab与CST联合仿真的具体方法

对于

CST建模仿真有两种方式，一种是使用交互界面直接建模仿真，另一种是通过脚本语言Visual Basic for Applications (VBA)操作。使用交互界面进行的所有操作，均会在CST的历史命令里写入相应的VBA语句。所以，只要使用Matlab将VBA命令语句写入CST历史命令就可以实现建模仿真等一系列操作。Matlab调用CST建模仿真流程如

图 3 Matlab调用CST仿真流程图 Fig. 3 Flowchart of CST simulation called by Matlab

Matlab调用CST建模仿真的部分语句如

表 1 Matlab和CST联合仿真的基本语句

Table 1 Basic codes of Matlab and CST co-simulation

Function

Code

Calling CST

cst=actxserver(‘CSTStuio.application’)

Creating new MWS file

mws=cst.invoke(‘NewMWS’);

Saving MWS file

mws.invoke(‘saveas’, ‘Bessl.cst’);

Saving MWS file

mws.invoke(‘AddToHistory’, VBA);

4 仿真过程及结果

采用Matlab与CST联合进行建模仿真。首先在Matlab中将‘CSTStudio.application’作为参数传入‘axtxserver’函数，返回CST的COM服务器。使用‘invoke’方法新建微波仿真文件，建立仿真模型，以及设置仿真参数。通过Matlab调用CST，在金膜上建立p=2 μm，半径方向的波矢为kr=2π/p。背景材料设置为‘normal’，边界条件为‘open’，使用CST中的时域求解器进行求解，所用的入射光波长为λ=632.8 nm。

图 4 CST中的等宽悬链线阵列 Fig. 4 Equal-width catenary arrays in CST

y=0时，x-z平面光强分布图。从仿真结果可以看到明暗交替的条纹，并且中心为亮条纹，符合零阶贝塞尔光束的光强分布特点。x-y平面的光强分布和相位分布图。其中x-y平面光场相位随着距离发生变化，而光强基本保持不变。

图 5 CST仿真得到的x-z平面Bessel光束光强分布 Fig. 5 Intensity distribution of Bessel beam on x-z plane obtained by CST

图 6 x-y平面贝塞尔光束光强和相位分布。(a)~(d)距离结构一定距离处，x-y平面贝塞尔光束相位分布；(e)~(h)与相位对应的光强分布图；(i) ~ (l)图(e) ~ (h)虚线所示的光强截线图 Fig. 6 Intensity and phase distributions of Bessel beam on x-y plane. (a)~(d) Represent the phase distributions on x-y plane at certain distances from the structure; (e)~(h) Denote the intensity distributions corresponding to the phase distributions (a)~(d), respectively; (i)~(l) are the cross views of the intensity shown by the dashed lines in figures (e)~(h), respectively

5 结论

本文介绍了Matlab调用CST Microwave Studio完成联合仿真的具体方法，以及需要使用的核心函数。通过Matlab强大的矩阵处理与数值建模能力，调用CST建立了等宽悬链线狭缝。介绍了设计等宽悬链线的方法，并给出了等宽悬链线狭缝可以产生的几何相位分布。按照产生贝塞尔光束所需的相位梯度，用等宽悬链线阵列设计了超表面贝塞尔光束产生器件，并得到了零阶贝塞尔光束。除了零阶贝塞尔光束，应用等宽悬链线作为几何相位超表面同样的单元结构，还可以设计器件来产生高阶贝塞尔光束、涡旋光束以及聚焦透镜等。联合仿真的方法可以用于设计复杂的结构，也可以结合Matlab强大的数值优化能力，实现一体化结构设计与优化。另外，也可以用同样的思路实现Matlab与HFSS、FEKO等常用的电磁软件的联合仿真。