应用于微波钟的十六极线型阱中汞离子囚禁数值及仿真研究*

王朝全，涂建辉，党文强，崔敬忠

(1.兰州空间技术物理研究所， 兰州 730000；2.真空技术与物理重点实验室， 兰州 730000)

0 引言

基于多极线型阱的汞离子微波钟由于不依赖于庞大、笨重的激光系统使得它重量轻、体积小，并且具有稳定度高、漂移率小等特点，非常适合作为新一代星载原子钟，可用于新一代导航系统、深空探测等项目[1-2]。1966年H.G.Dehmelt等人首次提出将离子阱技术应用到频标上的想法[3]，并阐述了离子阱微波钟的概念。世界上第一台离子阱微波频标是1981年法国LHA实验室的M.Jardion 和C.Audion 在双曲面Paul离子阱中研制成功的[4]。1989年，美国喷汽推进实验室(JPL)首次提出用线型阱代替双曲面阱的方案，开发出了检测区和共振区合一的单区式线型阱并实现了汞离子的微波频标[5-6]。理论分析表明，对于相同量级的二阶多普勒效应，线型Paul阱可以比双曲面Paul阱多囚禁20倍的离子[7]。1999年JPL又首次将多极线型阱应用到了汞离子微波频标中[8]。与四极阱相比，多级线型阱等效囚禁势相对较宽并且在很大范围内势场平缓，离子云密度较低，离子的相互作用较小，二级多普勒频移将更小[9]。目前国内的主要研究单位有中国科学院武汉物理与数学研究所、清华大学、航天五院504研究所等。中科院武汉物数所研究起步较早，并且对汞离子微波钟的研究较为全面[10-11]。2013年他们采用最新研制的延展型线型阱结构，使空间用的汞离子微波钟的物理部分的体积减小到与铯束管相近。同年底，在国内率先实现了线型阱囚禁汞离子微波钟原理样机的闭环锁定。2014年他们又成功完成了线型阱囚禁汞离子的控制和检测，并且观察到了Ramsey振荡谱线。

综上所述，为实现汞离子微波钟，线型阱中汞离子囚禁的研究是重要的。首先，线型离子阱几何结构，特别是杆电极半径的优化得到了研究[12-14]。研究表明，可根据实验测定给出优化参量，在一定程度上消除实际情况下离子囚禁中的非谐势以及由此引起的离子微运动效应。多极线型阱中囚禁势的对称性破缺已经得到了实验证实，并利用数值模拟得到了合理解释，这将进一步提高设计离子阱中囚禁势的合理性。对多极阱中稳定电气参数的研究是建立在离子阱几何结构确定的基础上的，稳定的电气参数对于离子囚禁稳定性和粒子运动轨迹比如微运动效应影响较大。对于任意确定的离子阱结构或尺寸，寻求与之匹配的稳定电气参数是必须要研究的。然而，对于多级阱而言，由于其囚禁方程的复杂性，无法得到稳定区的解析解，只能通过囚禁方程得到数值解[9,15]。尽管如此，对多极阱中电势场的数值分析以及近似赝势模型的研究已有相关报道[16]。另外，囚禁离子冷却问题是实现高性能微波钟必须考虑的，基于碰撞冷却的缓冲气体冷却由于其简单可行、不受离子种类和能级结构的限制在离子阱频标中得到了广泛研究和应用[17-19]。其中，缓冲气体冷却汞离子过程中引进的碰撞频移以及弛豫过程得到了尤为重要的研究[20-22]。

离子阱离子囚禁的仿真研究目前主要集中在离子阱微波钟[9,17]、离子阱质量分析器[15,23-24]以及基于离子阱的量子模拟等方面。研究方法主要有两种，一种是基于四阶龙科库塔方法的Matlab数值计算[25]，这种方法计算较为精确，但是计算量很高，计算速度很慢，对于模拟时间较长的离子囚禁并不适用；另一种是基于有限差分网格计算方法的离子光学仿真软件，比如利用Simion软件或者Comsol软件来模拟离子囚禁等[26-27]。研究的离子阱类型主要有双曲型Paul阱、线型Paul阱、矩形阱等。仿真研究在离子囚禁稳定电压、离子阱系统参数、离子冷却条件等方面均对系统的性能优化提供了有效指导[23-27]。尽管如此，目前在离子阱微波钟方面仿真研究的线型离子阱多为四极阱，其目标是为系统搭建提供设计支持，仿真研究仅作为系统性能研究的一个前期准备工作，而在离子阱质量分析器方面仿真研究的离子阱结构多为矩形阱,并不涉及多极阱。目前在离子阱微波钟系统构建层面一般还停留在利用四极阱或者有些涉及到十二极阱，因此其中的仿真研究作为一个前期工作也仅是针对四极阱和十二极阱开展。总之，对多极阱，尤其是十六极阱中离子囚禁问题的系统性仿真研究还很少有文献报道。

本文在1.1节中首先给出了在多极阱中无缓冲气体条件下的离子囚禁理论，然后基于Matlab数值计算研究了十六极阱中汞离子囚禁的稳定电压和离子运动轨迹；在1.2节中首先给出了缓冲气体对离子冷却的理论，然后基于Matlab数值计算研究了缓冲气体的种类、压强和温度等因素对离子冷却的影响。在2.1节中基于Simion软件仿真研究了十六极阱中不同电气和几何参数下汞离子囚禁的稳定电压和离子运动轨迹；在2.2节中基于Simion软件仿真研究了缓冲气体的种类、压强和温度等因素对离子冷却的影响。在第3节中给出了本文研究的总结和讨论。

1 数值分析

本节首先从理论上根据势阱中囚禁势满足的拉普拉斯方程和牛顿第二定律推导出囚禁离子运动方程以及根据分子动力学理论推导出缓冲气体离子冷却方程。然后利用Matlab分别数值求解以上方程，根据所得数值解研究无缓冲气体和有缓冲气体两种情况下离子囚禁问题。

1.1 简单条件下(无缓冲气体)离子囚禁

本小节考虑无缓冲气体时离子囚禁情况。首先根据势阱中囚禁势满足的拉普拉斯方程和牛顿第二定律得出囚禁离子运动方程，然后依据该方程利用Matlab数值求解得出不同参数下该方程的解，从而研究无缓冲气体时离子囚禁稳定区及离子运动特征。

1.1.1囚禁离子运动方程

求解离子阱势函数有多种方法，如复变函数、无穷级数展开、赝势模型等[3,12,20,28-30]。在无穷级数展开方法中高阶项对计算势阱的修正将直接影响汞离子的微运动。尽管如此，实验证明当电极半径和阱半径的比率为某个值时高阶项修正可以被补偿[12]。本文中由于十六极线型阱的囚禁空间是电荷无源的，因此可以利用求解拉布拉斯方程的方法来推导十六极阱的势函数，进而由牛顿第二定律得出十六极阱中离子的运动方程。在极坐标系下，假设势函数为Φ(r,θ)，则对应的拉布拉斯方程可以表达为[9]

(1)

利用分离变量法可以得到方程式(1)的通解为

(2)

式(2)中，A,B,C,D为常数。对于十六极阱，若杆电极所加电压为Φ0=U0-V0cos (Ωt)，其中U0为直流电压，V0为射频电压，则式(1)的边界条件为

(3)

式(3)中，r0为阱半径，即阱中心到杆电极表面的距离。将式(3)边界条件带入方程式(1)，得到囚禁空间势场函数为

(4)

(5)

(6)

1.1.2第一稳定区及离子运动

首先利用Matlab编程，采用四阶龙科库塔方法对离子运动方程式(6)进行数值求解，得出汞离子囚禁稳定区并对在十六极阱中的运动特征进行理论研究。此时，方程式(6)中的m=199×1.67×10-27kg为汞离子质量，qe=1.602×10-19，C为汞离子所带电量。另外，设置十六极阱中心(原点)到杆电极的距离r0=5.0 mm。首先研究十六极阱中离子囚禁的第一稳定区(所加射频场圆频率以Ω=3.14 MHz为例)，设置的稳定条件为囚禁时间大于2 ms。通过数值计算得到结果如图1所示。

图1 十六极阱中第一稳定区

从图1中可以发现，比较于四极阱和双曲面型阱中稳定区的特征[9]，十六极阱的数值计算结果与四极阱中第一稳定区的解析解有着相同的形状特征，其稳定区都是关于直流电压U0=0 V(即a=0)对称，而与双曲面型阱中第一稳定区解析解的形状有明显不同。这主要是由于两种离子阱的结构不同造成的。双曲面型阱中电极结构的对称性以及所形成势场分布的对称性在径向(r方向)和轴向(z方向)两个方向上并不相同。而在四极或十六极的线型阱中电极结构对称性以及势场分布对称性在x轴和y轴两个方向上是相同的。另外，在图1中可以看到，这里数值计算得到的结果虽然可以给出第一稳定区的基本轮廓和特征，但并不能得到完美的第一稳定区(即在稳定区轮廓内有稳定点遗漏)，这是因为在十六极阱中稳定区对参数的依赖更加敏感，这种敏感性在非解析解的情况下例如数值计算或者软件模拟时表现的更为突出。图1中稳定区轮廓内的漏点正是对数值计算中的时间步长、计算精度等因素敏感性的体现。因此，由于这种较强的参数敏感性，实际中将很难达到数值计算与解析解的完美符合。除此之外，十六极阱的这种参数敏感性在实际的物理设计中也将表现出来，例如在实验中十六极阱对实际物理参数设计的偏差容忍度将比四极阱时小得多或者说十六极阱将对物理参数的要求更加苛刻。出现这种情况的原因从方程层面讲是由于在四极阱中离子运动方程在x和y方向是独立的，而在十六极阱中离子运动方程在x和y方向并不是独立的，而是相互影响的。从物理角度讲是由于四极阱中仅存在四极矩势场，而在十六极阱中除了四极矩以外还将激发多极矩势场。设计中由于有阱形缺陷和囚禁场不纯等因素将导致高阶非谐势的存在，非谐势对囚禁的离子具有加热效应，囚禁离子会从外界(主要是囚禁场)吸收能量使得离子动能增加，温度上升，当离子的动能增加到大于囚禁场的势阱深度时，离子就会被抛出阱外。特别是多极阱中的高阶非谐势对离子的加热效应更加明显，且越远离阱中心非谐势越强，这是十六极阱对参数依赖更加敏感的物理原因(高阶非谐势会使离子运动加速、温度升高，这将导致一方面离子无法被长期稳定囚禁，另一方面由于离子被加速，会增大多普勒频移，而使光谱展宽)。

在十六极阱中的稳定区内接着研究汞离子的运动轨迹。为了方便这里选择Ω=3.14 MHz时稳定区内的一对稳态电压U0=5 V，V0=400 V情况，将设定条件带入方程式(6)，借助Matlab软件，利用四阶龙科库塔数值迭代方法求解得到离子的运动轨迹，如图2所示。本文数值求解时的时间步长设定为10-8s(注：步长太大会降低方程的求解精度，步长太小会产生随机误差)。

图2 第一稳定区内离子稳定运动

从图2中可以看到，在2 ms内离子在x和y方向上的运动均是稳定的，离子的运动仍然是宏运动和微运动的叠加。另外，由于十六极阱中的势场分布是一个多势阱震荡，而四极阱中是单势阱震荡，所以与四极阱相比，十六极阱中随着时间离子宏运动的振幅并不相同，出现的幅值较多，正如这里数值计算的结果即图2所示。其实，如果考虑多离子囚禁时这也表现为两种阱中的粒子数分布不同，十六极阱中离子为环状分布，而四极阱中离子为球状分布。另外，由于在十六极阱中囚禁势场相对较宽且在很大范围内势场较为平缓，因此在十六极阱中离子微运动的影响将大大减小，如图2所示。其实，对于线型阱来说极数越高离子的微运动就越不明显。

1.2 复杂条件下(有缓冲气体)离子囚禁

本小节考虑缓冲气体存在时离子囚禁情况。首先根据分子动力学理论建立缓冲气体对离子的冷却模型，推导出缓冲气体冷却方程，然后依据该方程，利用Matlab数值求解得出不同参数下该方程的解，从而研究缓冲气体对离子的冷却作用。

1.2.1缓冲气体冷却方程

在离子阱中当缓冲气体分子与囚禁离子达到热平衡时，离子和缓冲气体分子碰撞能量的损失等于离子从电场中得到的能量，单位时间内动量改变的大小等于离子受到的电场力大小。此时离子的速度vm=kmE，其中km是离子迁移率，E为离子所在位置的电场强度。热平衡状态下离子受到缓冲气体的粘滞阻力Fm与受到的电场力Fq为一对平衡力，因此离子受到的粘滞力可以表示为[32]

(7)

(8)

1.2.2缓冲气体冷却

这里仍然采用四阶龙科库塔方法进行数值求解，对缓冲气体存在情况下的离子运动进行分析。首先研究不同种类的缓冲气体对离子冷却的作用。这里分析了Ω=3.14 MHz，U0=-5 V，V0=180 V即稳态点时缓冲气体压强为4.0×10-5Torr，缓冲气体温度为0.5 K的氦气、氖气和氩气对所囚禁汞离子的冷却情况。假设缓冲气体与汞离子的相互作用温度为800 K，则数值计算结果如图3所示。图3(a)，(b)和(c)分别为氦气、氖气和氩气的结果，实线表示离子x方向运动，虚线表示离子y方向运动。

图3 不同种类的缓冲气体对离子冷却的作用

从图3可以看到以上3类缓冲气体对离子均有冷却作用，而且在离子运动的各个方向均有冷却，特别是在相同温度和压强下氦气、氖气和氩气的冷却作用依次增强。这主要是由两方面决定的，一方面3类气体中氦气的分子半径最小，所以与汞离子的碰撞截面较小，而氩气分子半径最大，因此碰撞截面也最大。碰撞截面越大越有利于对离子的冷却。然而，另一方面氦气分子的质量最小，而氩气分子的质量最大，根据缓冲气体与离子碰撞过程中动量守恒可知，相较于囚禁离子的质量，缓冲气体分子质量越小冷却效果越好。因此图3中缓冲气体对离子的冷却作用其实是以上两方面因素的一个综合表现。此外，不同于四极阱的情况，十六极阱中离子被缓冲气体冷却后的稳态并不总是在阱中心的原点处，这是因为十六极阱中存在多个势阱并且绕着原点(阱中心)呈环状分布，而四极阱中只存在一个势阱且其中心在原点处，因此在四极阱中被缓冲气体冷却后的离子更趋于稳定在原点附近。从图3中还可以看到，离子被不同种类缓冲气体冷却后的稳态位置并不相同，这主要是因为一方面，在十六极阱中存在多个小势阱绕原点分布，当离子被冷却后的动能仍大于各个小势阱深度的时候，离子并不能被某个小势阱囚禁，而仅仅是被整个离子阱区域的整体势阱囚禁，使得离子无法逃脱到阱外，仍然被囚禁于阱中，此时离子将会在整个离子阱区域内的各个势阱之间穿梭震荡，并且震荡过程中不断被缓冲气体冷却减速，因此离子运动速度越来越小，在不同小势阱之间穿梭振荡的周期越来越长，如图3(a)中虚线所示。若离子在各个小势阱之间穿梭震荡过程中，被缓冲气体冷却减速到其动能小于各个小势阱深度的时候，离子将被囚禁于某个小势阱中，进而被缓冲气体进一步减速最终使得离子在某个小势阱中达到稳态。另一方面3种缓冲气体对离子的冷却能力不同，因此离子在各个小势阱之间穿梭震荡过程中被冷却减速到动能小于各个小势阱深度的位置和时间也不同，所以离子最终被囚禁的稳态位置也不同。

图4 缓冲气体在不同压强下对离子的冷却作用

图4中实线表示离子x方向运动，虚线表示离子y方向运动。从图4中看到，由于缓冲气体的存在，离子的震荡很快被削弱并稳定下来，并且随着缓冲气体压强的增大冷却速度和冷却效果愈加明显。另外，从图中可以看到，当无缓冲气体时离子的运动速度不会被减小，最终运动状态为动态稳定即离子在多个势阱中周期性穿梭，最终达到运动轨迹周期性震荡稳定，此时离子稳定运动的空间范围较大。然而当冲入缓冲气体时离子的运动速度很快被减小，最终运动状态将趋于静态稳定即离子最终被囚禁在某一个小势阱中，此时离子稳定运动的空间范围较小。因此，在十六极阱中缓冲气体起着减小粒子运动速度和压缩离子空间运动范围作用。

下面作者进一步研究了不同缓冲气体温度对离子冷却的影响，这里分别以氦气和氖气为例(压强均为4.0×10-5Torr)。氦气情况时数值计算结果如图5所示，图5(a)和(b)分别为温度0.8 K和0.5 K时的结果；氖气情况时数值计算结果如图6所示，图6(a)和(b)分别为温度1.5 K和0.5 K时的结果。在图5和图6中实线表示离子x方向运动，虚线表示离子y方向运动。

图5 氦气在不同缓冲气体温度下对离子冷却

从图5(a)中可以看到当缓冲气体温度比较高时离子受到缓冲气体的冷却作用使得在不同小势阱间穿梭的周期变长。尽管如此，离子仍然通过在不同的小势阱间震荡穿梭达到稳态，此时并没有被囚禁在单一势阱中。然而随着温度进一步降低，离子将很快被囚禁在单一小势阱中如图5(b)所示。从图6中可以看到缓冲气体的温度越低对离子的冷却作用越强，因为低温度的缓冲气体与囚禁离子发生碰撞，当两者之间达到热平衡的稳态后离子的温度将会被降低，离子运动的速度将会被减小。因此，温度越低的缓冲气体越能更快地减慢离子的运动速度从而达到稳态。

图6 氖气在不同缓冲气体温度下对离子冷却

2 仿真模拟

在这一部分中作者利用仿真模拟的方法研究了十六极阱中汞离子囚禁的相关问题。研究手段主要是借助于Simion电磁仿真软件以及LUA编程对十六极阱中离子囚禁的过程进行建模并对模拟结果进行分析。

2.1 第一稳定区及离子运动模拟

首先模拟十六极阱中离子囚禁的第一稳定区以及不同条件下第一稳定区的变化。下面分别给出了阱的不同电气参数(以电极射频频率为例)和不同几何参数(以阱半径为例)情况下对第一稳定区的模拟。模拟结果如下：当r0=5.0 mm时不同射频频率下第一稳定区模拟如图7所示，图7(a)和(b)分别为Ω=3.14 MHz和Ω=2.44 MHz时的结果；当Ω=1.74 MHz时不同阱半径下第一稳定区模拟如图8所示，图8(a)和(b)分别为r0=5.0 mm和r0=8.5 mm时的结果。

图7 不同电气参数时第一稳定区比较

从图7和图8中可以看到，随着频率和阱半径的减小，第一稳定区将变小，并且直流和射频电压的幅度都将减小，但是仍保持稳定区的形状不变，其稳定区特征与第1节中数值分析所得结果一致。尽管如此，在同一个稳定区内，随着射频电压变大稳定区的模拟逐渐开始出现稳定漏点，导致稳定区轮廓逐渐变得不再清晰，甚至出现了失真，这也是十六极阱中离子囚禁参数敏感性在仿真软件模拟精度上的体现。

图8 不同几何参数时第一稳定区比较

然后作者模拟了十六极阱中汞离子的稳定运动以及非稳定运动。这里选择Ω=3.14 MHz时第一稳定区内的稳定点U0=-7 V，V0=1 240 V的情况，其模拟结果如图9所示，图9(a)为x方向运动结果,图9(b)为y方向运动结果。另外，选择稳定区外的一点U0=-140 V，V0=1 240 V作为非稳定点，模拟此时离子的运动轨迹，其模拟结果如图10所示，实线表示离子x方向运动，虚线表示离子y方向运动。

图9 十六极阱中汞离子的稳定运动模拟

图10 十六极阱中汞离子的非稳定运动模拟

从图9中可以看到，这里模拟的离子稳定运动的轨迹特征与第1节中数值分析所得结果一致，这里不再赘述。从图10中可以看到离子非稳定运动经历了3个阶段即随着时间演化离子逐渐从存在微运动状态，到微运动消失状态，最终到离子逃离状态。此外，十六极阱中离子运动在x和y两个方向上的稳定性是关联的，一旦在某个方向上不稳定，很快将导致在另一个方向上也不稳定，两个方向上的不稳定几乎同时出现，如图10所示，原因是十六极阱中离子在x和y两个方向上的运动方程是不独立的，这与四极阱情况不同。

2.2 缓冲气体冷却模拟

图11 不同种类的缓冲气体对离子冷却

接着作者模拟了缓冲气体在不同压强和温度下对囚禁离子的冷却，缓冲气体以氦气为例。当缓冲气体温度Tgas=0.1 K时不同压强下缓冲气体对离子冷却的模拟结果如图12所示，其中L1线、L2线和L3线分别相应于缓冲气体压强为Pgas=3.0×10-5Torr，Pgas=7.0×10-5Torr和Pgas=2.0×10-4Torr的情况。与之对应的离子三维运动模拟如图13所示。模拟发现，当缓冲气体压强较小时，离子运动很快变得不再稳定进而逃离囚禁区域，如图12中L1线及图13(a)所示。随着压强增大离子将被囚禁在了十六极阱势场中的多个子势阱中，各个子势阱绕十六极阱中心呈环状分布，离子将在各个子势阱间穿梭达到稳态，如图12中L2线及图13(b)所示。如果进一步增大压强，离子将进一步被冷却，直至离子仅被囚禁在了十六极阱势场中的单一子势阱中，如图12中L3线及图13(c)所示。

图12 缓冲气体在不同压强下对离子冷却

图13 缓冲气体在不同压强下对离子冷却的三维运动

当缓冲气体压强Pgas=5.0×10-5Torr时不同温度下缓冲气体对离子冷却的模拟结果如图14所示，其中L1线、L2线和L3线分别相应于缓冲气体温度为Tgas=0.2 K，Tgas=0.05 K和Tgas=0.01 K的情况。与之对应的离子三维运动模拟如图15所示。从图14和15中可以看到，随着缓冲气体温度的降低，离子的运动出现了和缓冲气体压强增大时类似的特征，此处不再赘述。综上所述可知，此处对十六极阱中不同缓冲气体压强和温度条件下的离子运动模拟得到了和第1节数值分析时相同的结论。

图14 缓冲气体在不同温度下对离子冷却

图15 缓冲气体在不同温度下对离子冷却的三维运动

3 结语

作者分别从数值分析和仿真模拟两个角度研究了无缓冲气体情况下十六极阱中汞离子囚禁的第一稳定区、离子运动轨迹等，以及复杂条件下缓冲气体种类、缓冲气体压强和缓冲气体温度等因素对离子冷却的影响。通过笔者的研究得到了十六极阱中汞离子囚禁的某些动力学以及运动学特征，同时得到了十六极阱的几何参数、电气参数以及环境参数的设计要点。总结得到了与四极阱情况不同的一些特征，例如十六极阱离子囚禁中离子运动对设计参数特别敏感、离子阱内势场平坦粒子数密度低等。理论分析表明，与四极阱相比产生以上不同的根本原因是在多极阱中势场分布不同，其直接原因是在多极阱中离子运动方程高阶项的存在以及两个方向上运动方程的非独立性。总之，通过对十六极阱中离子囚禁技术的研究使作者掌握了多极阱中离子囚禁的多方面特征，例如十六极阱离子囚禁与四极阱有相同的第一稳定区特征、随着频率和阱半径的减小第一稳定区将变小、离子的微运动相对较弱且离子稳态震荡幅值较多、离子运动在x和y两个方向上的稳定性是关联的，两个方向上的不稳定几乎同时出现、在相同温度和压强下氦气、氖气和氩气对离子的冷却作用依次增强、不同种类缓冲气体或不同压强下离子最终被囚禁的稳态位置也不同、缓冲气体压强越大温度越低对离子的冷却能力越强等。优化离子囚禁的稳定电压直接关系到离子能否被稳定囚禁或者囚禁离子的寿命，优化离子阱中离子的分布可以提高光抽运和信号检测效率或者减小离子与微波作用时的二阶多普勒频移，优化离子运动轨迹可以减小由于离子微运动造成的加热效应，优化缓冲气体冷却条件可以在实现对离子有效冷却的前提下尽可能减少带来的额外碰撞频移。总之，通过研究作者可以从离子囚禁系统的优化设计角度提高微波钟的性能指标，最终为基于多极线型阱的汞离子微波钟实现提供理想的离子囚禁系统。