基于符号回归方法探索磁性斯格明子结构近似解析式*

史猛 王伟伟 杜海峰

1) (中国科学院合肥物质科学研究院,强磁场科学中心,合肥 230031)

2) (中国科学技术大学,合肥 230026)

3) (安徽大学物质科学与信息技术研究院,合肥 230601)

1 引言

磁性斯格明子是一种特殊的拓扑磁结构[1],如图1(a)和图1(b)所示,存在于具有由反转对称性破缺诱导的Dzyaloshinskii-Moriya (DM)相互作用的手性铁磁材料或多层膜体系中[2–5].由于其独特的拓扑保护性、非易失性、小尺寸以及低能耗的驱动特性,磁性斯格明子在自旋电子器件领域具有广泛的应用前景,尤其是在不同条件下的静态与动态特性引起了广泛兴趣.磁性斯格明子的特性与其自身的结构密切相关,因此研究其在不同相互作用下的结构变化具有重要意义.

图1 磁性斯格明子结构 (a) 布洛赫型斯格明子;(b) 奈尔型斯格明子;(c)单磁畴壁Fig.1.Magnetic skyrmion structures: (a) Bloch-type skyrmion;(b) Néel-type skyrmion;(c) magnetic domain wall.

目前,一般的磁结构如磁畴壁(图1(c))已经有了精确且普适的结构解析式,这些解析式可以直观地展现不同相互作用对磁畴壁结构的影响.然而,对于磁性斯格明子而言,并不存在这样一般性的精确解析式.研究者们提出了一些近似公式来描述其基本结构.例如,对于较大半径的磁性斯格明子,Braun[6]基于显式一维畴壁结构提出了初步模拟方法,这在结构和动力学特性的检测中得到了广泛应用[7–10].在易轴磁各向异性模型的情况下,Komineas 等[11,12]针对不同半径的磁性斯格明子获得了渐进意义上精确的结构图像,甚至扩展到中等半径大小的磁性斯格明子[13].然而,上述结果仍然是基于经验性公式的近似推导.随着对磁性斯格明子的实验研究日益精细,上述公式逐渐难以满足解释实验现象的需求.因此,如何更高效地获取更精确、更普适的磁性斯格明子结构解析式,以及这些解析式的具体形式,都迫切需要研究者们的深入探讨.

当前,随着人工智能技术的迅猛发展,利用机器学习等人工智能技术从实验和计算模拟产生的数据中捕捉不同物理量之间的关系[14–17],构建能够合理解释体系变化的模型和公式,在物理研究领域得到了广泛应用[18–21].在众多机器学习方法中,符号回归方法是一种强大的数据驱动模型构建方法.与一般的参数回归方法不同,符号回归方法无需预先指定模型或公式的初始结构.它通过遗传规划等方法,能够同时搜索数学解析式的形式以及内部参数的数值,因此在很大程度上不依赖于既有知识,却能够生成一系列在理解物理体系方面具有参考价值的公式.对于那些变量之间联系模糊、整体变化理解不够深刻的物理体系,符号回归方法成为一项强有力的工具.

本文一方面提供了通过符号回归方法处理磁结构数据,获得更精确、更普适的描述磁性斯格明子的结构解析式的新思路;另一方面将符号回归方法实际运用到由微磁学模拟以及数值计算得到磁结构数据处理中,获得描述二维磁性斯格明子的结构近似解析式.首先概述研究所需的磁性斯格明子理论研究成果以及符号回归方法的整体框架构建;其次描述并讨论通过符号回归方法处理磁结构数据的结果;最后做出总结.

2 模型与方法

2.1 磁性斯格明子理论模型

本文考虑了具有交换相互作用、DM 相互作用、磁各向异性以及外部静磁场的二维铁磁性材料薄膜.微磁结构通过单位磁化矢量m=m(x,y)描述,即m2=1,其体系能量表达式为[13]

其中,A,D和K分别是交换相互作用常数;DM相互作用常数以及磁性各向异性常数,外部静磁场h=H/Ms,其中Ms为 饱和磁 化强 度,eDM为块体形式下的DM 相互作用单位矢量:(∂µm×m),µ=1,2.同时,引入特征长度[13,22]:

其中,lex为交换长度;lw为畴壁宽度;lD给出了强DM 相互作用的螺旋畴周期.在微磁模拟或数值计算中,为了简化问题和降低运算成本,通常将所研究的体系划分成一个个的网格且一个网格区域内的所有磁矩被简化为一个磁矩,而交换长度lex决定了相邻磁矩之间存在相互作用的最大距离,lw和lD则给出了一个具体磁结构的大小,因此为了计算结果符合实际,网格尺寸不能大于以上三个特征长度.由于在参数设置符合实际的情况下,三个特征长度中最小的通常是交换长度lex,因此本文使用lex作为长度的单位,并依此将磁各向异性常数和DM 相互作用常数无量纲化,且假设外部静磁场垂直于铁磁性材料薄膜表面:

从而体系能量表达式修正为

使用磁化矢量的球面角 (Θ,Φ),具有轴对称的磁性斯格明子解可以通过极坐标 (r,φ) 来表示,其中角Θ满足下列方程:

下文将使用这样的约定: 在斯格明子中心Θ(r=0)=π,在空间无穷大r→∞处,Θ=0.

2.2 符号回归方法

随着二十世纪末计算机技术和机器学习的迅猛发展,当时的研究者开始探索一种能够自动编写软件代码的编程解决方案,而符号回归方法是其中相对较新的算法结构之一.这一方法最初是由John Koza[23]基于他的导师John Holland 所提出的遗传算法改进而来.遗传算法基于达尔文的进化论,采用种群进化和优胜劣汰的方式来选择最优模型.然而,遗传算法主要用于参数优化,需要预先提供一个已知的公式形式,因此难以解决寻找未知公式以描述数据内在联系的问题.John Koza 创新地采用二叉树结构来表示公式形式(图2(a)),并通过遗传算法对公式形式和内部参数进行优化,从而找到最佳解.因此,符号回归方法无需预先设定公式形式,研究者只需提供一组数据,该方法便能从数据中总结出能够描述数据内在联系的表达式.

图2 符号回归算法中二叉树结构与遗传操作示意图 (a) 用二叉树表示公式形式;(b) 复制操作;(c) 交叉操作;(d) 变异操作Fig.2.Binary tree structure and schematic diagrams of genetic operations in symbolic regression: (a) Representing the formula in binary tree form;(b) copy operation;(c) crossover operation;(d) mutation operation.

符号回归方法搜索最优公式的大致流程如图3 所示,具体步骤如下:

图3 符号回归方法流程图Fig.3.Flowchart of the symbolic regression.

1)设定初始算术符号集,确定运行过程中的设置条件如结束条件等;

2)生成初始种群(即初始公式集);

3)按预测结果与数据结果之间的误差大小以及公式的复杂度来评估初始种群;

4)通过复制、交叉以及变异操作(图2(b)—(d))从原种群中生成新的种群;

5)重复步骤3)和步骤4),直至满足结束条件,并根据评估结果将生成的公式集进行排序.

采用SymbolicRegression.jl 包(SR.jl)[24]来搜索磁性斯格明子结构解析式,并用均方差来评估所得公式对于数据的符合度,其形式为

其中,N为样本数量;yi,o与yi,p分别为样本真实值与预测值.同时根据公式复杂度,即

式中,Nvar,Nconst与Nop分别为一个公式中的变量数量、常数数量与算符数量,E为所评估的公式.从低到高列出每个公式复杂度下符合度最高的公式与其对应的均方根误差大小,来对搜索到的公式进行筛选,即帕累托最优分析,从而避免过拟合的出现.

3 利用符号回归算法求解磁性斯格明子结构近似解析式

3.1 通过搜索一维单磁畴壁结构解析式验证符号回归算法

首先为了验证符号回归算法对于搜寻具体磁结构解析式的能力,使用一维单磁畴壁解析式:

可以得到磁畴壁的结构数据,并用符号回归算法基于该结构数据来搜索公式,结果显示,符号回归算法成功搜索到了(8)式.

为进一步增加挑战,采用JuMag 微磁学模拟程序,通过改变交换相互作用常数和磁各向异性常数,生成了一系列的一维单磁畴壁结构数据,以供符号回归算法进行公式搜索.然而,微磁学模拟得到的磁畴位置是未对齐的,通过中心校准算法与符号回归算法的耦合,运用符号回归算法进行公式搜索,结果显示符号回归算法成功地搜索出(8)式.图4 展示了符号回归算法搜索得到的(8)式与经过处理的一维单磁畴壁结构数据的拟合情况.

图4 不同A 与K 值下,(8)式拟合结果与一维单磁畴壁结构模拟数据比较图 (a) A=1×10–12 J/m,K=1×103 J/m3;(b) A=5×10–12 J/m,K=2×103 J/m3;(c) A=13×10–12 J/m,K=3×103 J/m3Fig.4.Comparison between the fitting results of Eq.(8) and simulation data of one-dimensional magnetic domain wall under various values of A and K: (a) A=1×10–12 J/m,K=1×103 J/m3;(b) A=5×10–12 J/m,K=2×103 J/m3;(c) A=13×10–12 J/m,K=3×103 J/m3.

3.2 使用符号回归算法搜索磁性斯格明子结构近似解析式

关于磁性斯格明子结构数据,通过Matlab 软件中的bvp4c 函数求解二阶常微分方程((5)式)来获取.边值条件为Θ(r=0)=π ;r→∞处,Θ=0,其中r→∞在实际计算中用r=5 来替代.同时,为了简化问题的复杂度并考虑到磁各向异性对于磁性斯格明子结构的影响较为特殊,因此在生成磁性斯格明子结构数据时,未将磁各向异性作用考虑在内,即κ=0.通过改变DM 相互作用以及外部静磁场的大小,生成了56 组二维磁性斯格明子结构数据,并分别运用符号回归算法来搜索解析式.根据结果以及帕累托最优分析(图5),总共发现了两个相对合适的磁性斯格明子结构近似解析式,这两个解析式的形式为

其中,c表示任意常数.在不同参数下,搜索到的两个解析式在磁性斯格明子结构数据的符合度方面互有胜负.为了进一步深入研究这两个解析式谁更符合磁结构数据的条件,对所获得到的数据做了综合分析.发现在特定情况下,即λ2/h≤0.16 时,(10b)式符合更好;而在另一情况下,即λ2/h>0.16时,(10a)式更为吻合.具体结果如图6 所示.

接着,基于参数的设定对搜索的公式进行初步分析.首先,两个解析式共同的特点是它们都包含双曲正切函数,且这一函数位于解析式计算的最后一步,这一点确定了磁性斯格明子的基本结构,与磁性材料中的交换相互作用相关联;这两个解析式的不同点是双曲正切函数中包含的另外一种一元函数.在DM 相互作用与外部静磁场相比较小时,该一元函数为对数函数,该函数在拐点处的变化较为急剧,表明虽然DM 相互作用导致了磁性斯格明子的生成,但由于外部静磁场是主要因素,DM相互作用的影响受到压制,从而导致磁性斯格明子尺寸较小.Romming 等[7]的研究中也展示了斯格明子尺寸随着垂直外磁场的增大会逐渐减小的现象.而在DM 相互作用与外部静磁场相比较大时,该一元函数为平方根函数,其在拐点处的变化较为平缓,这表明随着DM 相互作用的增大,其影响范围扩大,逐渐占据主导地位,因此磁性斯格明子尺寸增大.在Wu 等[25]的研究中,表明了斯格明子的尺寸与A和D的比值大小有关,随着A/D值的减小,斯格明子尺寸逐渐增大,即意味着D的增大会导致斯格明子尺寸的增大,这和本文的结果是一致的.

随后,再次运用符号回归算法,针对先前搜索出的两个解析式中的常数项关于λ和h来进行公式搜索,从而探寻两个解析式中的常数项与这两者之间的关系.将得到的关系式分别代入(10a)式与(10b)式中,其结果为

其中,

接下来,选取了4 个λ2/h的关键数值所对应的磁性斯格明子结构数据,分别代入(11a)式与(11b)式,如图7 所示,拟合结果仍然相当良好.因此,将DM 相互作用、外部静磁场以及磁性斯格明子结构近似解析式紧密联系起来,得到了一个更为普适的结论.

图7 不同 λ2/h 大小下 (11a)式或(11b)式的拟合情况 (a) 0.01;(b) 0.16;(c) 0.167;(d) 0.9Fig.7.Fitting results of Eq.(11a) or Eq.(11b) under various λ2/h values: (a) 0.01;(b) 0.16;(c) 0.167;(d) 0.9.

最后为了检验所得(11a)式与(11b)式的泛化能力,本文又另外生成6 组数据进行拟合,其结果如图8 所示.(11a)式与(11b)式仍然拟合得较好,说明符号回归算法搜索到的这两个公式具有较强的泛化能力,这进一步证实了符号回归方法在搜索磁性斯格明子结构近似解析式方面的强大能力.

图8 不同 λ2/h 大小下,(11a)式或(11b)式的拟合情况 (a) 0.009;(b) 0.015;(c) 0.139;(d) 0.192;(e) 0.227;(f) 0.769Fig.8.Fitting results of Eq.(11a) or (11b) under various λ2/h values: (a) 0.009;(b) 0.015;(c) 0.139;(d) 0.192;(e) 0.227;(f) 0.769.

4 结论

本文首先通过搜索一维单磁畴壁结构解析式,初步验证了符号回归算法的应用潜力,随后通过求解二阶常微分方程,获取不同DM 相互作用和外部静磁场大小下的磁性斯格明子结构数据.借助符号回归算法,对这些数据进行公式搜索,得到了两个适合的结构近似解析式: 在λ2/h≤0.16 时,(8b)式拟合效果较佳;在λ2/h>0.16 时,(8a)式表现更优.同时,对解析式中函数的构成进行了分析: 当DM 相互作用与外部静磁场相比较小时,磁性斯格明子结构与对数函数相联系;在DM 相互作用与外部静磁场相比较大时,磁性斯格明子结构与平方根函数相联系.不过,无论何种情况,双曲正切函数始终决定了磁斯格明子的基本结构.最后,再次运用符号回归算法,将上述两个解析式中的常数项与DM 相互作用以及外部静磁场相关联,以使上述两个解析式具有普适性.本文提出了使用符号回归算法研究磁性斯格明子结构的新方法,同时通过符号回归算法初步探索了磁性斯格明子结构近似解析式,验证了符号回归算法在该研究领域中的强大能力.