线型离子阱中钙离子库仑晶体结构和运动轨迹模拟*

张洪硕 周勇壮 沈咏 邹宏新

(国防科技大学,理学院量子信息研究所,长沙 410073)

离子阱中囚禁的离子在满足库仑耦合条件下,会形成库仑晶体,其结构分布和运动轨迹由离子阱的参数和离子种类决定.本文采用分子动力学模拟软件LAMMPS 和(py)LIon 程序包,仿真了40Ca+在线型离子阱中库仑晶体的形成过程,以及不同位置离子的微运动和宏运动轨迹.另外,本文还对混入少量同位素离子(44Ca+)和氢化钙离子(CaH+)后的40Ca+晶体结构进行了仿真,并对混入前后形成的库仑晶体结构变化进行对比和分析,期望能对离子阱实验中形成的暗离子进行识别和处理.

1 引 言

1953 年第一个四极质谱仪由德国科学家Paul等[1]提出,1959 年Fischer 实现了Paul 阱[2],同年Dehmelt 等实现了Penning 阱[3].当离子阱系统满足库仑耦合条件时,阱中囚禁的离子会形成库仑晶体,这是一种由经典电磁相互作用导致的相同带电粒子的固定排布,但在实际应用中又会受到量子力学的影响[4].

库仑晶体的研究和应用前景十分广泛.在非线性动力学方面,库仑晶体中存在着粒子的混沌行为,晶体结构稳定性与晶体频率和射频频率的比例相关[5].在热力学和统计物理方面,从一维离子串到二维锯齿形结构再到三维椭球结构相变过程中,出现的缺陷位置符合Kibble-Zurek 模型的预测,进而引出了对称性破缺问题的讨论,而单个杂质离子在离子串中的热运动涉及到的是Frenkel-Kontorova模型[6,7].在核物理方面,Penning 阱中形成的库仑晶体中存在杂质离子碰撞现象,这与致密恒星热核聚变类似[8].在腔量子电动力学方面,库仑晶体中的离子由于几乎静止,且离子彼此之间和离子与环境被较好的隔离,其内部电子态可以被高精度控制,因而可以作为单光子水平下原子和光相互作用的理想研究对象[9].在量子信息处理方面,当纠缠态下的离子数目达到边带冷却的囚禁极限时,需要保持在不同的离子阱中或相同离子阱中的不同囚禁位置,来满足不同区域特定离子之间的门操作等,因而库仑晶体可以应用于量子位日益增多的系统中[10].在量子模拟方面,库仑晶体作为量子比特比二进制的经典计算机有着天然的优势,同时也满足系统中离子稳定和与环境隔离的条件.在量子逻辑光谱学方面,通过激发双离子库仑晶体的振动模式,将光谱学离子的电子状态信息传递到逻辑离子,对逻辑离子进行激光冷却和状态检测,进而实现对光谱离子进行超分辨率光谱分析[11].

关于库仑晶体基本结构性质的研究,在实验方面,库仑晶体结构参数的测量主要涉及径向和轴向几何参数、宏运动和微运动频率以及离子阱参数与空间尺度关系等[12].通过调整离子阱参数,可以实现库仑晶体从一维至三维的相变,其中三维库仑晶体具有螺旋结构[13].多组分库仑晶体往往采用协同冷却的方式,由于经常选取相对质量相差较大的离子,因而具有明显的分层现象,如24Mg+,40Ca+等.也有选取较高电荷数离子的双组份晶体,如Ho14+和Be+,其排布方式与离子的荷质比和空间电荷效应相关[14,15].另外,调控离子阱参数可以精密控制库仑晶体,减小附加微运动等[12,16,17].在理论和仿真方面,现有研究主要涵盖了不同类型的离子阱中形成的库仑晶体结构,如四极阱、Penning阱、八极阱等[18].另外还有多组分库仑晶体的仿真,离子数量较多且荷质比相差较大时,库仑晶体有同心球壳的结构,而荷质比相差较小时分层不明显[19−22].在仿真严格的真空相变条件、分子动力学模拟上的优化等方面也有一些研究[23,24].在以往库仑晶体仿真方面,往往偏重于库仑晶体结晶后的状态,对不同位置和在冷却过程中离子运动轨迹研究相对较少.值得一提的是,在离子阱实验中经常遇到的暗离子问题,目前已开展了很多相关研究,如基于暗离子的数目变化的离子反应过程研究[25],对协同冷却作用下单个暗离子的成像研究[26],对激光冷却下与暗离子相关的碰撞系数研究等[27].本文针对多离子协同冷却的动力学问题,在库仑晶体中用暗离子替换了特定位置的离子,研究了离子晶体的动态演化过程.

LAMMPS(large-scale atomic/molecular massively parallel simulator)是由美国桑迪亚国家实验室联合天普大学开发的基于C++的开源经典分子动力学模拟软件包,已被广泛应用于物理学、化学、材料科学、生物学等多个领域[28].由于LAMMPS 的建模需要自行建模或者第三方软件导入,因而Bentine 等[29]开发了一款基于Python 的针对线型离子阱中库仑晶体仿真的外部程序包(py)LIon,本文工作便是在此基础上完成的.

本文采用分子动力学模拟软件LAMMPS 和(py)LIon 程序包,基于线型离子阱囚禁钙离子开展了库仑晶体结构和运动轨迹的模拟仿真.对单组分的钙离子库仑晶体形成过程中的动能和温度变化进行分析,也对多组分荷质比相差比较大的库仑晶体结构进行仿真.钙离子阱中可能混入同位素离子(44Ca+)和氢化钙离子(CaH+),分析其对库仑晶体结构的影响,结果表明,在一定条件下有可能通过对离子阱中库仑晶体成像实现对暗离子种类的识别,这对基于离子阱的相关应用和实验研究具有一定的指导意义.

2 线型离子阱模型

关于线型离子阱的物理结构,主要是由加载高频交变电压的四极杆和两端的端盖直流电压组成,另外在实验中往往会在四极杆的多个位置施加额外的直流修正电极,囚禁场电势的表达式为(1)式,其中Urf为射频电压、Ωrf为射频电压的频率、2r0为离子阱的内径、Uend为端盖电压、2z0为有效极杆的长度(端盖电压之间的距离)、κr为径向几何参数、κz为轴向几何参数,x和y为径向方向,z为轴向方向.其径向和轴向宏运动的频率表达式如(2)和(3)式:

在线型离子阱中形成库仑晶体还需要施加激光冷却、郎之万热浴等方式进一步限制离子的运动,如(4)式所示:

离子在满足等离子体耦合参数Γ≥ 175 的情况时,在离子阱中可以形成库仑晶体,其中Q是粒子电荷,a是与粒子密度相关的Wigner-Seitz 半径,kB是玻尔兹曼常数,T是粒子温度[30].

本文通过设定离子阱的参数并生成LAMMPS输入文件,传输给LAMMPS 运行得到输出结果,而后在Python 中对仿真结果进行进一步地分析处理[29].在(py)LIon 的主程序中可将单线程修改为多线程,这在运行多离子体系时可以提升一定的运行速度.(py)LIon 工具包将激光冷却对离子的作用力Fl等价为与速度成正比的反向阻尼力(常量阻尼参数β),如(5)式,大部分仿真中一般小于βmax≈π2ℏ/λ2[22,29,31].而这在仿真中计算冷却极限时并不够严谨,因而根据参考文献[32]中多普勒激光冷却的半经典理论,在原工具包中加入新的激光冷却函数如(6)—(8)式和自发辐射反向冲量函数如(9)和(10)式,其中Δ为激光失谐量,Γ为自然线宽,Ω为拉比振荡频率,k为波矢,ρee为激发态概率[32].由于Fr的方向是随机的,在仿真中可以通过随机的球坐标中方位角φ和仰角θ的来给定方向.另外,还存在着离子吸收光子的离散过程对冷却力的影响,形式与Fr相同,同样有着加热作用.[32]在冷却激光弱饱和(共振饱和参数s0=2Ω2/Γ2→0),且激光失谐量Δ为自然线宽Γ的1/2 的状态下,当加热功率与冷却功率相等时,可以得到多普勒冷却极限,如(11)式,其中 ℏ 为约化普朗克常数[32,33].

仿真得到每个时间步长下的离子坐标和瞬时速度,进而生成离子的三维运动轨迹.由于库仑晶体中离子的宏运动动能贡献了主要的总动能,因而可以计算多个微运动周期的动能平均值来衡量库仑晶体的总动能和温度随时间的变化.[22]具体由(12)和(13)式给出,其中J为一个RF 周期(微运动周期)所包含的时间步长数量,α取x,y,z三个方向,N为总离子数,〈〉为多个RF 周期下的速度平均值,Esec为宏运动动能,Tsec为宏运动温度[22].另外动能与温度的关系如(14)式:

3 仿真结果与分析

3.1 钙离子库仑晶体

使用参考文献[34]的参数,采用激光冷却方式模拟了34 个40Ca+形成的库仑晶体.具体参数为Uec=1 V,Urf=252.4 V,Ωrf=2π×5.634 MHz,r0=3.5 mm,z0=5.0 mm,κz=0.17[34].针对阱中的离子施加一束在x,y和z轴上均有相同分量的激光.根据阱参数定义时间步长(timestep)为8.875×10–9s,总时间为107个时间步长,1 个微运动周期包含20 个时间步长,而当采样间隔为微运动周期的倍数时,可以得到仅包含宏运动的轨迹.绘制了位于顶部(靠近晶体的旋转对称轴)和中部偏下(偏离旋转对称轴)的两个离子的运动轨迹,为方便之后的描述,分别将其记为离子1 和离子2.

图1(a)所示为库仑晶体形成稳定后,在10 个微运动周期下的平均位置,离子1 和离子2 分别用红色的“×”和“+”来标记.图1(b)是两个离子加载了微运动的宏运动轨迹在径向平面上的投影,总时间为20 个微运动周期(1 个宏运动周期,即离子往复运动的一个周期),采样间隔为1 个时间步长,由浅至深的灰度表示轨迹随时间的演化,左侧为离子1 的运动轨迹,右侧为偏离z轴的离子2 轨迹.图1(c)是最后100 个微运动周期下,两个离子到旋转对称轴的距离随时间的变化.由于离子1 的运动轨迹经过旋转对称轴,因而其距离变化的曲线(红色)相对平滑,而离子2 的曲线(蓝色)有明显的锯齿.从图中可以得出离子1 的宏运动轨迹更加对称,离子2 的微运动幅度要强于离子1.

图1 激光冷却下不同位置离子的微运动和宏运动轨迹比较 (a)34 个40Ca+组成的库仑晶体在10 个微运动周期下的平均位置,“×”和“+”所在位置为选取的离子1 和2;(b)两个离子在径向平面上加载了微运动的宏运动轨迹,灰度代表时间演化;(c)两个离子与旋转对称轴的距离随时间的演化,红色平滑曲线是离子1 的变化曲线,蓝色曲线是离子2 的变化曲线Fig.1.Micromotion and secular motion of ions at different positions under laser cooling: (a) The average position of Coulomb crystals composed of 34 40Ca+ in 10 micromotion cycles,the positions of × and+are the selected particles 1 and 2;(b) the secular motion trajectory loaded with micromotion on the radial plane of the two particles,the gray scale represents the time evolution;(c)the evolution of the distance between the particle and axial rotational axis of symmetry evolving over time.The red smooth curve is for particle 1,and the blue curve is for particle 2.

为了进一步分析,图2 是在不同的演化时间下离子1 和2 与旋转对称轴的距离变化,8×106个时间步长下的变化曲线由黄色上三角标记,107个时间步长由红色“×”标记,1.5×107个时间步长由蓝色圆形标记,2×107个时间步长由绿色下三角标记.可以观察到其宏运动的振幅随着时间演化会不断增大,这说明在库仑晶体形成一段时间后,整体温度会缓慢升高,同时说明库仑晶体存在寿命，这也可能是由于算法误差积累导致的.

图2 激光冷却下离子1 和离子2 的径向微运动和宏运动对比（1000 个时间步长）.不同颜色表示不同演化时间后的运动轨迹：0.8×107 个时间步长(黄色上三角),107 个时间步长(红色“×”),1.5×107 个时间步长(蓝色圆形),2×107 个时间步长(绿色下三角)Fig.2.Comparison of radial micromotion and secular motion of ion 1 and ion 2 over 1000 timesteps under lasercooling.Different colors represent the trajectory of motion after different evolution times: 8×106 timesteps (yellow),107 timesteps (red),1.5×107 timesteps (blue),2×107 timesteps(green).

由于库仑晶体早在106个时间步长之前就已经大致形成,因此本文绘制了不同时刻下库仑晶体的温度变化如图3,最终冷却温度可以降至3.02 mK,时间范围是0—4×105(时间步长),采样间隔为1000 个时间步长,(8)式中用于平均动能的微运动周期数量为10,其中在初始段存在温度上升的阶段，这是由于初始条件并非稳定状态而造成的.

图3 激光冷却下34 个40Ca+形成的库仑晶体温度随时间的变化曲线,时间范围是0—0.4×106 个时间步长Fig.3.Temperature of Coulomb crystals formed by 34 40Ca+under laser cooling,the time range is 0 —0.4 × 106 timesteps.

值得一提的是,本文模拟了单个离子在线型离子阱中的冷却过程,使用了上文提到的新加入到工具包中的激光冷却函数和加热函数,其中40Ca+的397 nm 激光对应能级跃迁的自然线宽Γ约为22.4 MHz,在激光弱饱和的状态下,得到其冷却极限为0.536 mK,其冷却过程的温度和速率变化如图4 所示.

图4 单个40Ca+在线型离子阱中的冷却过程(a)温度随时间的变化(b)速度随时间的变化,蓝色、绿色和红色曲线分别对应在x,y 和z 方向上的速度平均值变化Fig.4.The cooling process for a single 40Ca+ in a linear ion trap: (a) The change in temperature over time;(b) the change in speed over time.The blue,green,and red curves correspond to changes in the average velocity in the x,y,and z direction respectively.

3.2 多组分库仑晶体

在之前的研究中有不少是通过协同冷却的方式对多组份库仑晶体进行的研究,这在离子的荷质比相差比较大的情况下,库仑晶体才会有比较明显的分层排布[22].本文仿真了4 种荷质比的离子,Q∶M分别为1∶40,1∶80,1∶120 和1∶160,离子阱参数设定与之前相同,初始每种离子各125 个,同样采用协同冷却的方式,对荷质比为1∶40 的离子施加激光冷却,其他离子在库仑相互作用下也会得到冷却,演化时间为3.5×106个时间步长.图5 为4 种不同荷质比的离子形成的库仑晶体,用不同的颜色形状来区分.图中可以明显看到其分层现象,越靠近中心其荷质比越大,在三维结构上有类似圆柱壳的结构,而在径向平面投影上有近似圆形的排布方式.本文进一步计算了径向平面投影上的类圆排布结构的平均半径和间隔,如表1 数据.由计算可知,最外层荷质比为1∶160 的离子,没有严格按照1∶1 的间隔方式排布,可能的原因是最外层离子数目不够充足,虽然形成了类似圆柱壳结构的排布方式,但是轴向方向上部分离子更加靠近中心部分,导致整体结构向外凸,在径向平面上壳壁厚度增大,在轴向平面上圆柱高度变短.

表1 4 种荷质比离子组成的类圆排布结构信息Table 1. Information on the circle-like arrangement of four charge-to-mass ratio ions.

图5 4 种不同荷质比的离子在库仑晶体中的排布.离子种类分别为1∶40(蓝色圆形),1∶80(红色下三角),1∶120(黄色“×”),1∶160(绿色上三角)(a)为三维图片;(b)为轴向平面投影;(c)为径向平面投影Fig.5.The arrangement of ions with four different charge-to-mass ratios in the Coulomb crystal.Ion types are 1∶40 (blue circle),1∶80 (red lower triangle),1∶120 (yellow ×),1∶160 (green upper triangle) separately: (a) Three-dimensional image;(b) axial plane image projection;(c) radial plane image projection.

3.3 Ca+库仑晶体中混入相关杂质对晶体结构的影响

自然界钙的同位素中占比最多的前两位是40Ca 和44Ca,丰度分别为96.941%和2.086%[35].在钙离子阱实验中可能会引入同位素44Ca+和氢化离子CaH+,由于荷质比相差不大,其分层结构不太明显.如图6 所示,首先仿真了125 个40Ca+、8 个CaH+和8 个44Ca+形成的库仑晶体.在径向平面投影上CaH+和44Ca+均位于库仑晶体外侧,而44Ca+位于更靠近外侧.

图6 多组分库仑晶体结构,包含36 个40Ca+(蓝色圆形),8 个CaH+(红色下三角形)和8 个44Ca+(黄 色“×”)(a)三维图片;(b)轴向平面投影;(c)径向平面投影Fig.6.Multicomponent Coulomb crystal structures,36 40Ca+ (blue circle),8 CaH+ (red upper triangle) and 8 44Ca+ (yellow ×):(a) Three-dimensional image;(b) axial plane image projection;(c)radial plane image projection.

由于同位素离子和氢化离子形成的暗离子无法直接观测,可间接分析其对周围离子排布的影响.本文仿真了8 个40Ca+离子形成一维离子串的库仑晶体,如图7 所示.在3 种离子组合的情况下进行了仿真,第1 组为8 个40Ca+,第2 组为7 个40Ca+和1 个CaH+,第3 组 为7 个40Ca+和1 个44Ca+.从轴向角度来看,离子串轴向排布间隔整体上没有太大的变化,而从径向角度来看,44Ca+对周围40Ca+离子的排布影响较大,造成了纳米量级上的位置变动(最大为1.70 nm),但目前的实验手段无法观测到.而对于CaH+的混入,在成像上只会观察到出现了暗离子,但基本对40Ca+的排布基本无影响.

图7 由8 个40Ca+形成的离子串在两种混入杂质的情况下排布方式的对比,两种情况分别为混入1 个氢化离子(CaH+)或者1 个同位素离子(44Ca+)Fig.7.Comparison of the arrangement of the ion string formed by pure 40Ca+ and ion string containing 1 hydride ion (CaH+) or 1 isotopic ion (44Ca+).

如图8 所示,本文仿真了16 个40Ca+离子形成二维平面结构的库仑晶体,同样分为3 种情况.图8(a)是16 个40Ca+形成的库仑晶体,红色标记“×”为另外两种情况下替换40Ca+的位置,绿色实线是选取在图8(b)中进行对比的部分.图8(b)中第1 组为16 个40Ca+,用蓝色点和连线标记,第2 组为15 个40Ca+和1 个CaH+,用黄色点和连线标记,第3 组为15 个40Ca+和1 个44Ca+用绿色点和连线标记.另外,CaH+用红色 “+” 标记,44Ca+用红色“×”标记,横坐标为离子在二维库仑晶体平面的位置.同样从轴向角度看,在整个库仑晶体的轴向间隔没有太大的变化,而径向角度来看,44Ca+对离子串的排布影响较大,对附近40Ca+排布的最大影响为0.453 µm,这在目前的实验手段是可以观察到的,而CaH+的影响相对较小,最大为0.154 µm.

图8 由16 个40Ca+形成的平面库仑晶体在两种混入杂质的情况下排布方式的对比.两种情况分别为混入1 个氢化离子(CaH+)和1 个同位素离子(44Ca+)(a)三维图片;(b)轴向平面投影的对比Fig.8.Comparison of the arrangement of the planar Coulomb crystal formed by 16 pure 40Ca+ and Coulomb crystal containing 1 hydride ion (CaH+) or 1 isotopic ion (44Ca+): (a) three-dimensional image;(b) axial plane image projection.

多个离子的库仑晶体呈现出类似DNA 的螺旋结构[13].本文仿真了36 个40Ca+,分别混入1 个CaH+或1 个44Ca+所形成的库仑晶体,如图9 所示,按照与图7 相同的方式进行了标识,区别是选取对比的部分是被替换的离子所在的螺旋线.在径向方向来看,通过计算替换离子对所在螺旋线上所有40Ca+造成位置变化，得到44Ca+造成最大的影响为0.314 µm,而CaH+最大影响为0.246 µm.值得一提的是,通过遍历所有离子的位置变化,暗离子对周围离子排布影响最大的离子并不在螺旋线上,而是在轴向上最靠近暗离子的位置,其位置变化为0.905 µm.

图9 由36 个40Ca+形成的三维库仑晶体在两种混入杂质的情况下排布方式的对比.两种情况分别为混入1 个氢化离子(CaH+)和1 个同位素离子(44Ca+)(a)三维图片;(b)螺旋结构的对比Fig.9.Comparison of the arrangement of the three-dimensional Coulomb crystal formed by 36 40Ca+ and Coulomb crystals containing 1 hydride ion (CaH+) or 1 isotopic ion (44Ca+): (a) Three-dimensional image;(b) comparison of helical structure.

4 结 论

本文使用分子动力学模拟软件,仿真了线型离子阱中形成的40Ca+库仑晶体的相关内容,主要为冷却过程中不同位置离子的运动轨迹、动能和温度随时间变化的分析,等数目的多组分库仑晶体的仿真以及关于40Ca+库仑晶体中引入杂质后的排布结构变化等3 个方面.

对于单组分库仑晶体,靠近旋转对称轴的离子宏运动轨迹更加对称,幅值也更小.而对于等数目的多组分库仑晶体,在荷质比相差较大且离子总数较多的情况下,虽然有明显的分层排布方式,但最外层离子的排布方式并不严格按照圆柱壳的方式排布,这对于实验中多组分离子库仑晶体成像也具有一定的指导意义.在40Ca+库仑晶体中引入杂质的排布结构方面,从一维到三维结构的都进行了仿真.在一维离子串中混入的暗离子对周围离子的影响较小,而在二维平面和三维螺旋结构中的暗离子对周围离子造成的排布位置变化较为明显,在目前实验精度下是可观测的,通过对比分析或许可以用于区分暗离子的种类,进一步提高离子阱的性能.

总的来说,当下通过离子阱搭建的量子系统占据着现代物理的重要地位,库仑晶体虽然由经典物理条件约束,但在量子系统中或许有着不可估量的应用前景,因而本文关于库仑晶体在3 个方面开展的仿真工作对今后的实验和仿真具有一定的参考价值,尤其是在离子阱中暗离子的识别等方面.