专利名称：磁斯格明子晶格和磁涡旋晶格的量子力学模拟方法
专利号：201810526044.6
专利权人：刘照森

项目简介：本发明涉及计算物理和材料学领域，将量子理论用于模拟磁斯格明子晶格和磁涡旋晶格。包括：1、磁系统哈密顿量中的各自旋为量子力学算符；2、所有物理量都严格按照量子理论算出；3.采用自洽算法；4.模拟自高温逐渐降至低温度，以保证程序的正确收敛；5.采用周期性边界条件；6.考虑垂直外加磁场，使因DMI和海森堡交换作用产生的螺旋状条纹，变为斯格明子或涡旋晶格；7.在海森堡交换和DMI作用共存时，考虑Compass型各向异性作用，模拟出近邻涡旋方向相反的斯格明子晶格。克服了微磁学和蒙特‑卡洛二经典方法的不足，使得模拟十分灵活、便捷，尤其在量子尺度内，模拟出二经典方法无法算出的周期性磁结构。

磁斯格明子晶格和磁涡旋晶格的量子力学模拟方法
技术领域

本发明属于计算物理和计算材料学领域，尤其涉及磁斯格明子晶格和磁涡旋晶格的量子力学模拟方法。由于磁斯格明子和磁涡旋的大小常在纳米量级，且需要操控的电流甚小，所以它们被视为下一代磁存储和磁逻辑运算的；理想候选者，因而具有重要的理论和应用价值。

背景技术

数十年来，国内外研究者广泛采用蒙特-卡洛(Monte Carlo)[1]和微磁学(Micromagnetics)[2，3]两种数值方法模拟磁性材料的微观磁结构和研究其宏观物理性质。此二法在科学研究中取得的成果及其在计算物理中的地位是举世公认的，国内外研究者运用它们发表的文献数不胜数，在此无需引证。

然而，这二种方法都建立于经典物理之上：磁系统中的自旋或磁矩被当做长度不变的经典矢量。这显然会对磁系统的精确描述、计算速度和最终的计算结果产生不良的影响。

为了达到系统的平衡态，蒙特-卡洛法采用Metropolis算法，运算中的每一步都开始于旋转各自旋的空间取向，以期降低系统的总能量；在运行数万甚至数百万个循环后，对系统中每个自旋的矢量值求平均，以确定其大小和空间取向；然后再计算整个系统的磁化强度、磁化率、比热等宏观物理量。

微磁学模拟法通常把磁体分成许多网格，每个网格的磁矩设为Mi，并假设它们的大小MS在整个磁体中处处相等，Mi随时间的变化规律满足Landau-Lifshitz-Gilbert方程[2，3]。微磁学通常不考虑温度效应，即仅研究零温的特殊情况。为了克服这一缺陷，Skomski等人把磁矩的大小MS以及系统参量K1等作为与温度相关的量[4]。但是，如何确定MS(T)和K1(T)二函数又成了新的问题。此外，诸如纳米等有限系统，其中各处的磁矩大小显然不同。所以，微磁学的描述是不够完备的。

蒙特-卡洛和微磁学二方法的经典局限性可能大大影响收敛速度和模拟结果的正确性。可以发现许多文献中用蒙特-卡洛法算出的磁结构并不具应有的对称性。例如，国外学者考虑磁偶极矩之间的相互作用，用二法算出的圆型纳米盘上的磁结构都不对称[5，6]；还有，采用蒙特-卡洛根本无法模拟出非零温度下无限二维正方格子上的反铁磁斯格明子晶格[7，8]。

为克服上述二方法的经典局限性，本人近年来研发了一种新的量子模拟方法。因其使用自洽算法(Self-Consistent Algorithm)，故简称为SCA方法。在此量子模型中，系统哈密顿量中的自旋或总角动量是量子力学算符，而不是经典矢量；任意温度下各种物理量，如磁化强度、系统总能量、总自由能、热容量等都严格按照量子理论算出。至今，笔者用此量子方法模拟了纳米颗粒、纳米线、纳米管、纳米盘等[9-16]，且模拟结果与实验和理论相符[10，12]。如考虑海森堡交换作用(Heisenberg exchange interaction)与Dzyaloshinsky-Moriya作用(DMI)共存的磁性纳米盘，模拟出了铁磁(FM)和反铁磁(AFM)的涡旋。

在文献[12，13]中，本人利用量子SCA法模拟出纳米盘上的铁磁和反铁磁涡旋。由于纳米盘的有限尺度，故仅算出至多两个磁涡旋。文献[13]题目中使用反铁磁斯格明子一词，但是实际仅仅算出两个反铁磁涡旋，并非周期性的晶格(crystal)。2017年6月25日，本人在第三届全国凝聚态物理会议上的报告题目中也仅涉及纳米盘上的多个磁涡旋而已[16]。

磁斯格明子存在于反演对称性破缺的磁性表面或多层金属的界面上，它通常因Dzyaloshinsky-Moriya作用而产生。磁斯格明子的直径在几个纳米到100纳米之间。它们微小的尺度，及易被微弱电流驱动的优点，被看作下一代理想的磁存储和磁逻辑单元。这些年来，国内外学者都运用微磁学和蒙特-卡洛法，对斯格明子进行理论研究。然而，在多层磁材料的界面上形成的直径仅仅几个纳米的斯格明子[17-20]才具有极高的应用价值。因微磁学基于连续模型，它尚能描述大尺度的斯格明子，但在此微小的尺度内，连续模型不再适用，所以微磁学的描述不再精确和可靠[21]。另一方面，因微观尺度上显著的量子效应，蒙特-卡洛也变得不理想，如用它无法模拟出正方晶格上的反铁磁斯格明子晶格[8]。因此，模拟和研究纳米尺度大小的磁斯格明子和磁涡旋，迫切需要运用量子理论写出的新软件。本发明就是为此目的提出和展开的。


发明内容

本发明的目的，一方面是要克服蒙特-卡洛和微磁学二方法的经典局限，另一方面则是要大大拓展本人量子SCA方法的应用，尤其是模拟和研究纳米尺度的斯格明子晶格和涡旋晶格，从而为这一具有重要应用前景的热点问题提供理论方法。

量子力学理论的应用，使得本发明能够细致、正确地描述任意温度下和外磁场中系统磁结构，特别是磁斯格明子晶格和磁涡旋晶格的变化规律，达到为实验和现代技术应用服务的目的。

根据量子力学理论，磁系统的哈密顿量可写为



其中，和表示位于第i个格点和第j个格点处的自旋；Jij和分别表示和之间的海森堡交换作用和DMI作用的强度；KA表示单轴垂直各向异性作用强度；表示第i个格点处磁系统表面的单位法向矢量；μB表示波尔磁子；gJ为朗德因子；为外加磁场强度。式(1)中第一项表示Heisenberg交换作用，第二项表示DMI作用，第三项表示单轴垂直各项异性，最后一项为系统在外加磁场中的能量。DMI存在于反演对称性破缺的表面或界面上，它具有手性，是形成磁斯格明子和磁涡旋的主要原因。如果(表示第i个自旋到第j个自旋的矢量)，则斯格明子为涡旋状的Bloch型；如(表示垂直于表面/界面的单位矢量)，则斯格明子为Néel型。

研究中根据需要，还常常计入磁偶极矩之间的长程相互作用：



本人在研究中发现，对于有限尺度的铁磁纳米盘，这种磁偶极矩之间的长程相互作用也可诱发磁涡旋。

为了克服蒙特-卡洛和微磁学的经典局限，根据量子理论，本人在过去的研究工作已采取了以下技术方案：

(1)引入量子理论，磁系统哈密顿量中的自旋不再是经典矢量，而是量子力学算符；

(2)采用自洽算法，模拟中不必每一步都计算和比较系统的总能量，计算会自发收敛于平衡态；

(3)所有物理量都严格按照量子理论计算。

为模拟磁系统中周期的斯格明子/磁涡旋晶格，还采取了以下技术方案：

(4)采用周期性边界条件，而不是自由边界条件，使得模拟无穷大系统成为可能；

(5)必要时考虑垂直的外加磁场。

这里必须指出采用周期性边界条件的重要性。实验中观测到，斯格明子的直径可小到10nm。所以，宏观上看来很小的磁性薄膜内部存在数以万计以上的周期分布的斯格明子或磁涡旋。然而，计算机的内存和速度是有限的，要使得模拟可行，就必须采用周期性边界条件。否则，如采用自由边界条件，如在垂直强磁场的作用下，有限尺度的磁性薄膜上可能仅仅出现一个磁涡旋或斯格明子。此外，为使周期性边界条件得到满足，需要在预计算中，初步确定磁结构空间周期的近似量值。

DMI是诱发斯格明子的重要原因。但是，如果此作用甚强，与Heisenberg交换作用竞争的结果却产生螺旋状(helical)的基态磁结构。研究表明，施加适量强度的垂直外磁场，有助于产生周期的斯格明子晶格。

(6)在海森堡交换和DM作用共存的情况下，再考虑Compass型的各向异性作用。此时系统的哈密顿量为：



这里不存在外加磁场的作用，<ij>表示相互作用限于相邻自旋之间，第三项表示单自旋各向异性，最后一项即为Compass型各向异性。如使用参量J＝1K，，D＝1.5K，A1＝0.5K和A2＝-2K，可模拟出周期的斯格明子阵列，且相邻斯格明子涡旋方向相反。

(7)由于系统中各种相互作用的激烈竞争，磁过程变得十分复杂，要模拟出正确的磁结构，需要在系统达到平衡态后，再运行数百次，对结果求平均。如此便可消除热扰动造成的偏离，算得的磁结构更加可靠，且具有完好的对称性。

在介绍了以上各点之后，我们在附图一中绘出在某一强度固定的垂直磁场中，模拟系统磁结构(如斯格明子)随温度变化的程序设计流程图。现在说明如下：

第一步：考虑量子化的系统哈密顿量。

第二步：模拟从高温T＝T0开始，初始化磁结构，使得所有自旋取向在空间随机分布。

第三步：采用周期性边界条件。

第四步：根据量子力学公式，逐一计算每个自旋的热平均值。

第五步：计算前后二次循环算出的自旋或系统总能量的差值。如果某个相对差值大于给定小量δ，则返回第三步。

第六步：计算和输出系统微观磁结构，及其它宏观物理量。

第七步：降低系统温度，即令T＝T-ΔT。

第八步：判断是否到达最低温度(T≤Tf)，否则返回第三步。

第九步：模拟结束。

在以上步骤，可看到权利要求中的各点：

C1、根据量子理论使磁系统哈密顿量中的自旋为量子力学算符；

C2、模拟始于高温，产生三个随机数分别赋予自旋的三个分量，模仿初始态的自旋取向随机分布；

C3、采用周期性边界条件以模拟无穷大系统的周期性磁结构；

C4、采用自洽算法，量子理论的运用使得运算能自发地收敛于系统的平衡态；

C5、判断前后两次循环算出的自旋差是否小于δ，如是，则满足收敛条件，输出磁结构(如周期分布的斯格明子晶格或涡旋晶格)；如否，则返回步骤C3。

C6、模拟自高温的顺磁态开始，逐渐降低温度，以保证计算收敛于正确的平衡态；

C7、判断逐渐降低的温度是否为预设的最低温度，如是，则模拟结束，如否，则返回步骤C3。

模拟中采用了前述的技术方案：

本发明的技术方案：所述磁系统哈密顿量中的自旋为量子力学算符，而非经典矢量

本发明的技术方案：所述量子力学模拟方法在计算的过程中所有物理量都严格地按照量子理论算出。

本发明的技术方案：所述量子力学模拟方法中因量子理论的引入，采用自洽算法运算便能自发地收敛于系统的正确平衡态。

在蒙特-卡洛模拟中，每一步都需要人为地旋转各自旋，计算因旋转引起的能量变化，按照一定的几率确定此旋转操作是否被接受，以期系统收敛于最终的平衡态。然而，因量子理论的引入，本发明采用自洽算法，运算能自发地收敛于系统的正确平衡态，使得算法更为便捷、灵活。

本发明的技术方案：在海森堡作用和DM作用共存的系统中，再考虑Compass型各向异性作用，模拟出近邻涡旋方向相反的斯格明子晶格。

本发明的技术方案：模拟自高温的顺磁态开始，然后逐渐降低至温度；因在高温，各自旋的幅值甚小，大的热能足以使得量子系统克服势垒，隧穿至能量的最小点，到达系统的平衡态；模拟中逐渐降低温度，便能保证每一温度下都实现正确的收敛。

本发明的技术方案：采用周期性边界条件，以小的系统模拟大的、以致无穷的系统成为可能。模拟微小尺度的磁系统，通常采用自由边界条件；模拟尺度较大的磁系统，则需要运用周期边界条件，使得采用小的系统模拟大的、以致无穷的系统成为可能；否则，边界的影响太大，就模拟不出周期性的磁结构。

本发明的特征是：在磁系统中微观结构和宏观物理性质随温度及外加磁场的变化而变化。而微磁学通常仅仅研究零温的特殊情况。

本发明的特征是：所述量子力学模拟方法模拟出的系统微观磁结构通常具有严格的周期性和完美的对称性；相较于经典模型，当磁条纹或者斯格明子晶格的周期为数个纳米时，量子模型更能细致、精确地描述磁系统。

本发明的有益效果是：由于引入量子理论，本发明克服了微磁学和蒙特-卡洛经典方法的不足，且程序会自发收敛于系统平衡态，使得模拟十分灵活、便捷和快速。

特别是：微磁学仅模拟零温的情况，本发明可模拟任意温度下的情况；微磁学把样品中所有点的磁矩大小看作相等，本发明则可细致地算出材料中各点处磁矩的不同幅值及方向。

此外，蒙特-卡洛法无法模拟出无限二维正方结构系统中的反铁磁斯格明子阵列[7，8]，但利用本发明的方法，却能很容易地算出。

用本法模拟，在临界温度附近需要较多次循环达到如10-6的相对精度；但是运算到较低温区，往往仅经过数次循环，即可达到上述精度。例如最近模拟56×56个自旋构成的二维网格，从高温到低温模拟30温度点，用笔记本ThinkPad T430S，仅需要约1.15个小时。

附图说明


图1是模拟程序设计流程图。

图2是本发明实施例提供的Bloch型(a)铁磁、和(b)反铁磁斯格明子阵列。

图3是本发明实施例提供的Néel型的(a)铁磁、和(b)反铁磁斯格明子阵列。

具体实施方式

首先，要模拟二维正方晶格上的磁斯格明子或涡旋阵列，需采用式(1)给出的力量子哈密顿量。其中，如Jij＞0，系统为铁磁耦合；当Jij＜0时，系统为反铁磁耦合。如采用(矢量沿着二自旋连线的方向)，可模拟Bloch型斯格明子；如用(矢量在二维平面内，且与二自旋的连线垂直)，则可模拟Néel型斯格明子。

以下给出本发明的四个应用实例。模拟中都仅考虑最近邻自旋间的相互作用，且大小处处相等，即Jij＝J，Dij＝D。

实施例一：Bloch型的铁磁斯格明子阵列(Ferromagnetic Skymion Crystal of Bloch-Type)

模拟中选用30×30的正方格子，每个格点上有一个S＝1的自旋。为了用它模拟无穷二维系统，采用了周期性边界条件。此处不考虑垂直各向异性的作用，再设J＝1K，令D/J＝1.02733，以使斯格明子之间的周期距离λ＝10。

从模拟结果可以看到：在没有外场的情况下，此二维系统的自发基态为沿着[110]方向的条状结构；逐渐增大垂直外加磁场当着0.11Tesla≤B≤0.27Tesla时，便在正方格子内产生正六角密排结构的斯格明子晶格；在弱磁场中，此阵列的空间周期为10，与理论相符；逐渐增大外加磁场，正方格子内的斯格明子数目，从最初的18个逐渐减少到8个，而且，整个晶格形状随之旋转或扭曲；而当B＝0.27Tesla，整个阵列又恢复到正六角密排的结构。

图二(a)绘出B＝0.11Tesla，T/J＝0.1K时产生的六角密排结构的铁磁斯格明子阵列，它具有完美的几何对称性，空间周期为10，与理论和实验的结果完全相符[22]。

外磁场沿着z方向，每个斯格明子中心区域的自旋磁矩的z分量却沿着-z方向，而其周边区域的自旋磁矩的z分量取向相反，所以是严格的铁磁斯格明子。此外，每个斯格明子都是顺时针的，却被四个逆时针涡旋围绕着，从而降低系统的总能量，实现稳定的磁结构。

实施例二：Bloch型反铁磁斯格明子阵列(Antiferromagnetic Skymion Crystal of Bloch-Type)

因此例研究反铁磁情况，故取J＝-1K，D＝1K(于是D/J＝-1)，KA＝0[8]。采用56×56的正方格子，每个格点上有一个5＝1的自旋。利用周期边界条件，以模拟无穷大二维系统。模拟中发现，当着垂直外加磁场强度满足3.9Tesla≤B≤4.1Tesla时，在T/J＜1.6的温区内能够形成反铁磁斯格明子阵列；稍微增大或减弱外加磁场，反铁磁斯格明子图晶格便消失，而被反铁磁结构取代。

图二(b)示出在B＝4Tesla，T/J＝0.1时模拟出的反铁磁斯格明子晶格。它具有完美的对称性，空间周期λ＝7。每个斯格明子中心区域自旋的z分量平均值最小；而沿着[110]或者[-110]方向，二相邻斯格明子中心连线中点处的自旋的z分量值最大；相邻自旋磁矩的xy投影取向相反。这些都赋予斯格明子晶格的“反铁磁”特征。

实例一中B1＝0.11Tesla，实例二中B2＝4Tesla，且二者D/J的比值相近，但是B2/B1≈36.4。因此，要观测到反铁磁斯格明子阵列，需要施加36倍的强磁场。此外，反铁磁磁斯格明子晶格仅存在于很窄的B区间，因而缺乏很好的稳定性。这些都使得实验观测困难了许多。

有必要指出，为了模拟正方格子上的反铁磁斯格明子，R.Keesman等人令J＝-1，D/|J|＝1，H/|J|＝4，仅在8×8的有限大小的方格子上模拟出单个斯格明子，却模拟不出无限正方系统上的反铁磁斯格明子晶格[8]。然而Tretiakov等人的理论研究表明，此种斯格明子存在于非零温度，在外加磁场中会更加稳定[23]。因此，此实例说明了本发明的重要意义，即当着斯格明子很小时，经典理论的描述已不再精确、可靠，需要运用量子理论和方法。

实施例三：Néel型铁磁斯格明子晶格(Ferromagnetic Skymion Crystal of Néel-Type)

模拟图三(a)和图二(a)使用了相同的参数，即J＝1K，D/J＝1.02733，T/J＝0.1，垂直外加磁场B＝0.11Tesla。用于模拟的正方格子大小为30×30，每个格点上的自旋S＝1，采用周期边界条件以研究无穷大的二维平面系统。有趣的是，尽管二种斯格明子分别为Bloch型和Néel型，但二晶格的空间周期都为λ＝10，30×30的方格上也出现18个斯格明子，也形成正六角密排结构。在每个斯格明子的中心区域内，各自旋的z分量都沿着-z方向，与外磁场方向相反；其周围区域的各自旋z分量与外场方向相同，所以都是典型的斯格明子。与图二(a)不同之处在于，斯格明子内的自旋的xy投影都收敛于(指向)其中心。

实施例四：N6el型反铁磁斯格明子晶格(Antiferromagnetic Skymion Crystal of Néel-Type)

模拟图三(b)和图二(b)使用了相同的参数，即J＝-1K，D/|J|＝1，B＝4Tesla。尽管二图中的斯格明子分别为Bloch型和Néel型，但二阵列空间的周期都为λ＝7，这些斯格明子也都形成正方结构。同样，图三(b)中每个斯格明子中心区域自旋的z分量的平均值最小；而沿着[110]或者[-110]方向，二相邻斯格明子中心连线中点处的自旋的z分量值最大；相邻自旋的xy投影取向相反。这些都使得斯格明子阵列具有“反铁磁”的特征。与图二(b)不同之处在于，除边缘区域外，各斯格明子内的自旋的xy投影都沿着径向反向排列。

有必要再次指出，R.Keesman等人令J＝-1，D/|J|＝1，H/|J|＝4，仅在8×8的有限大小的方格子上模拟出单个反铁磁斯格明子，却未能模拟出无限正方结构系统上的反铁磁斯格明子晶格[8]。然而Tretiakov等人的理论研究表明，此种斯格明子存在于非零温度，在外加磁场中会更加稳定[23]。因此，此实例再次说明经典方法的局限性，以及本发明的重要意义。

总的说来，图二(a)与图三(a)具有对偶性，图二(b)也与图三(b)具有对偶性，表明自然规律的完美对称性，也从另一方面说明我们的量子模拟方法和计算结果的正确性。

随着温度和外加磁场强度的变化，斯格明子及其晶格的形状、周期也随之变化。因篇幅所限，这里仅列出四张典型的代表图。更为详细的描述，将在未来发表的论文中给出。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

附注

1.关于量子力学算符

根据量子力学理论，在笛卡尔直角坐标表象中动量三个分量的算符



要求角动量的算符，则只需把此式中动量的三个分量换成算符即可。

同样，动能算符可由动量算符求得，



而能量算符则为



单个粒子的定态哈密顿量



其中表示势函数，其定态薛定谔方程为：



另一方面，如果选取一组正交、归一的完备基矢量，则力学量(动量、角动量、哈密顿量)的算符可以用厄米矩阵表示，它们的本征值为实数。

2.关于力学量平均值的计算

根据量子力学，在状态ψ下，力学量A的平均值



如果采用矩阵表示，则这里Ψ+是波函数Ψ的转置共轭矩阵。例如自旋的三个分量的平均值，就可如此求出。

如考虑温度效应，则需要先计算粒子在温度T下的配分函数



这里εn为单自旋哈密顿量的能量本征值，kB为波尔兹曼常量。状态n对力学量贡献的权重为：

wn＝exp(-εn/KBT)/z(T)。

3.关于斯格明子晶格的模拟

用有限大小的格子模拟无限大系统，需要采用周期性边界条件。

如果DM作用甚弱，则诱发简单的铁磁或反铁磁结构；如DM作用较强，如不考虑外加磁场，则模拟出螺旋状的条纹。对于铁磁交换作用，垂直外加磁场需满足



才会算出斯格明子晶格。但是，对于正方结构的二维材料，要模拟出反铁磁斯格明子晶格，必须采用本发明的量子方法才能成功。