基于分子的磁性系统中的动态滞后现象：路径概率方法在正方晶格上混合自旋（3/2, 7/2）Blume-Capel伊辛系统中的应用

《Journal of Magnetism and Magnetic Materials》：Dynamic hysteresis in molecular-based magnetic systems: application of the path probability method to the mixed-spin (3/2, 7/2) Blume-Capel Ising system on a square lattice

　　穆斯塔法·根恰斯兰|贝扎·塔什
土耳其凯塞里市埃尔吉耶斯大学物理系，邮编38039

摘要
本研究利用路径概率方法，对二维正方晶格上具有混合自旋（3/2, 7/2）的Blume-Capel伊辛系统中的动态磁滞行为进行了全面研究。通过系统地改变温度、晶场强度、磁场频率以及动力学速率常数，我们分析了这些参数如何影响动态磁滞回线形态、矫顽场和剩磁。研究结果揭示了多种类型的磁滞回线，包括单环、双环和三环结构，这些差异取决于系统条件及交换相互作用（无论是铁磁还是反铁磁情形）。重要的是，该方法能够捕捉到传统模型中忽略的动力学效应，有助于理解系统在振荡磁场下的动态磁响应与磁滞行为。这些发现与先前理论研究结果以及快速处理磁性材料中的实验观测趋势具有定性一致性。本研究为设计具有可调动态特性的功能磁性材料提供了重要依据。

1. 引言
基于分子的磁性材料因其在自旋电子学、量子计算和高密度数据存储等先进技术中的应用潜力而备受关注[1]–[5]。这类材料的应用价值很大程度上取决于对其在时变外部刺激下动态磁响应的理解，尤其是动态磁滞回线的行为。这些以矫顽场和剩磁为特征的回线，是衡量磁性系统能量损失和存储能力的关键指标，因此在功能材料的设计与优化中起着至关重要的作用[6]–[7]。

为模拟伊辛型自旋系统的动力学行为，人们发展出了多种理论方法，包括平均场理论[8]–[13]、有效场理论[14]–[15]以及蒙特卡罗模拟[16]–[24]。这些方法在简单晶格结构和有限自旋配置的研究中取得了重要成果，但在处理非平衡现象的动力学方面往往存在不足，或者需要大量的计算资源，尤其是在面对持续外部驱动场的系统中。

具有不同自旋值和晶场参数的混合自旋Blume-Capel伊辛系统为探索复杂的磁现象提供了理想载体。然而，现有文献中针对此类系统——尤其是处于振荡外磁场作用下的（3/2, 7/2）Blume-Capel伊辛模型——的系统性研究十分有限。这一局限源于自旋状态增多带来的挑战、相应动态方程的复杂化以及庞大的数值计算需求。尽管有一些研究探讨了这类系统的某些方面，但往往侧重于静态性质或较为简单的动态场景。

例如，Bouda等人[25]使用蒙特卡罗方法研究了在单离子各向异性和外磁场作用下的混合自旋（3/2, 7/2）铁磁伊辛系统。他们的研究包括计算基态相图，并分析晶场常数如何影响特定温度下的磁滞曲线，发现了包括单环、双环、三环和五环在内的多种磁滞回线形态。类似地，Mohamad和Zuhair[26]在纵向磁场作用下，利用正方晶格上的混合自旋（3/2, 7/2）Blume-Capel伊辛系统确定了其复杂结构的基态相图。Chakravorty和Chowdhury[27]以及Chakravorty[28]通过蒙特卡罗模拟计算了该系统的临界指数和磁热性质。Espriella等人[29]、[30]则同时运用蒙特卡罗方法和平均场方法，研究了在外磁场作用下的正方晶格混合自旋（3/2, 7/2）铁磁伊辛系统，分析了晶场和外磁场对磁化强度及磁滞曲线的影响，同样观察到了多环结构。Elidrysy等人[31]、[32]采用蒙特卡罗方法研究了在外磁场作用下的正方晶格混合自旋（3/2, 7/2）亚铁磁伊辛系统，探讨了系统参数对磁化强度、比热、补偿温度及磁滞曲线的影响，还进一步研究了其磁热性质。Eraki等人[33]利用蒙特卡罗模拟研究了具有混合自旋（3/2, 7/2）伊辛结构的陶瓷釉面（4–8）核壳纳米结构，分析了核心尺寸对磁性能、磁化强度及磁化率的影响。Kake等人[34]则利用精确递归关系，研究了在纵向磁场作用下、贝特晶格上的混合自旋（3/2, 7/2）亚铁磁Blume-Capel伊辛系统的临界行为、补偿行为及磁滞行为，给出了通过分析不同平面上序参量温度依赖性得到的相图，并观察到了多环磁滞曲线现象。最后，Ghafel和Mohamad[35]采用平均场近似方法，研究了在外磁场作用下正方晶格和三角形晶格上的混合自旋（3/2, 7/2）Blume-Capel伊辛系统，发现该系统表现出超顺磁性。

我们的研究采用了由Kikuchi[36]首次提出的路径概率方法，相较于以往主要依赖平均场理论、有效场理论或蒙特卡罗方法、且往往未充分考虑动力学效应的研究，该方法具有显著优势。这一创新方法使得我们能够更细致深入地探究动态参数如何影响磁滞回线的特征。因此，我们的工作为当前相关研究带来了重要进展，揭示了非平衡系统中微观动力学与宏观磁响应之间复杂的相互作用机制。路径概率方法提供了一个自洽的框架，可同时处理热力学过程和动力学过程，从而研究动态条件下时间依赖的占据概率及其演化规律[37]–[40]。该方法的独特优势在于能够引入不同的速率常数——这些参数可对应于诸如弛豫动力学或凝固动力学之类的物理现象，甚至可类比为快速淬火技术中的轮速参数[39]–[42]。正是这些特性，使得路径概率方法特别适合用于研究具有复杂相互作用的自旋系统。

近期研究已证实，路径概率方法在分析不同自旋配置和晶格结构下、处于振荡磁场作用下的混合自旋Blume-Capel伊辛系统的动态性质方面具有高效性[43]–[50]。在此基础上，本研究进一步拓展了路径概率方法的应用范围，旨在探究更高阶自旋系统中的复杂动态磁滞回线行为。此前已有研究主要聚焦于混合自旋（1, 7/2）和（3/2, 7/2）的Blume-Capel伊辛系统，并对其磁滞回线行为进行了系统分析[49]–[50]。在本研究中，我们将研究对象定为二维正方晶格上的混合自旋（3/2, 7/2）Blume-Capel伊辛系统，这是一种典型的二维参考几何结构（配位数z=4），常被用于研究伊辛模型中的非平衡磁滞与动态响应。在这一最简的二维参考结构中，对更高阶自旋系统的系统研究在理论文献中仍相对较少。此外，以往研究大多采用固定的交换相互作用参数，往往忽视了铁磁（正J值）和反铁磁（负J值）耦合情景对动态磁响应的影响，包括磁滞回线的形态和数量。在本研究中，我们并未试图完整描述系统的动态相变过程，而是集中分析动态磁滞行为，包括动态磁滞回线形态、矫顽场及剩磁。我们详细研究了这些量与关键系统参数之间的依赖关系，这些参数包括温度、晶场相互作用、外加振荡磁场的幅值与频率，以及动力学速率常数。因此，本研究的宗旨并非构建一个完整的DPT诊断框架，而是阐明所研究的混合自旋系统中动态磁滞特性的参数演变规律。

本文其余部分结构如下：第2节详细介绍二维正方晶格上混合自旋（3/2, 7/2）Blume-Capel伊辛系统的理论模型及路径概率方法的构建方式。第3节展示并解读数值计算结果，重点分析系统参数如何影响矫顽场、剩磁以及动态磁滞回线的形态。最后，第4节总结我们的主要发现，并提出未来的研究方向。

2. 方法与建模
我们研究的是由交替排列在正方晶格上的自旋3/2（σ）和自旋7/2（S）磁性离子构成的二维混合自旋Blume-Capel伊辛系统。此处选用正方晶格作为标准的二维参考几何结构，其配位数为z=4，便于在路径概率方法框架内清晰分析动力学效应；该结构并不等同于任何具体的实验系统。该晶格被划分为两个相互渗透的子晶格：子晶格A包含自旋σ=3/2的离子，子晶格B包含自旋S=7/2的离子。每个子晶格中的自旋都会与其在另一子晶格中的四个最近邻自旋发生相互作用。

该系统处于随时间变化的纵向振荡磁场中，磁场强度可表示为H=H0sinωt，其中H0为磁场幅值，ω为角频率。该系统的时变哈密顿量如下：
(1) ?=?J∑σiASjB?D∑iσiA2+∑jSjB2?H∑iσiA+∑jSjB。
其中，J表示不同子晶格上相邻自旋之间的双线性交换相互作用，其符号决定了磁耦合的性质：J>0表示铁磁耦合，而J<0则表示反铁磁耦合。第二项代表单离子各向异性，即晶场参数D，其正负及大小决定了非零自旋状态的能量高低。最后一项则描述了自旋与时变磁场之间的塞曼相互作用。

为研究系统的动力学演化过程，我们采用了由Kikuchi[36]最初提出的路径概率方法。该方法为研究系统在动态外场作用下的微观构型随时间变化的概率分布提供了强大工具。与静态或准平衡方法不同，路径概率方法能够同时处理动力学效应和热力学效应。

我们定义子晶格A的状态概率为XiA，子晶格B的状态概率为XjB，它们分别表示自旋处于相应状态的时间依赖概率。对于自旋3/2的子晶格A，共有4种可能的自旋状态（σ=±3/2, ±1/2）；而对于自旋7/2的子晶格B，则有8种可能的自旋状态（S=±7/2, ±5/2, ±3/2, ±1/2）。因此，整个系统由12个状态概率和10个相应的平均序参量来描述，这样的设定保证了在归一化约束下路径概率方程的封闭性。

相关的动态序参量包括：
• 对于A层：平均磁化强度mA≡<σia>，四极矩qA≡<σia2>，八极矩rA≡<σia3>；
• 对于B层：平均磁化强度mB≡，四极矩qB≡，八极矩rB≡，十六极矩oB≡，三十二极矩pB≡，四十六极矩tB≡，十八极矩uB≡。这些平均序参量可通过状态概率按以下公式表示：
(2a) mA=32X1A?X4A+12X2A?X3A
(2b) qA=94X1A+X4A+14X2A+X3A
(2c) rA=278X1A?X4A+18X2A+X3A
(3a) mB=72X1B+52X2B+32X3B+12X4B?12X5B?32X6B?52X7B?72X8B
(3b) qB=494X1B+254X2B+94X3B+14X4B+14X5B+94X6B+254X7B+494X8B
(3c) rB=3438X1B+1258X2B+278X3B+18X4B?18X5B?278X6B?1258X7B?3438X8B
(3d) oB=240116X1B+62516X2B+8116X3B+116X4B+116X5B+8116X6B+62516X7B+240116X8B
(3e) pB=1680732X1B+312532X2B+24332X3B+132X4B?132X5B?24332X6B?312532X7B?1680732X8B
(3f) tB=11764964X1B+1562564X2B+72964X3B+164X4B+164X5B+72964X6B+1562564X7B+11764964X8B
(3g) uB=823543128X1B+78125128X2B+2187128X3B+1128X4B?1128X5B?2187128X6B?78125128X7B?823543128X8B

在此建模体系中，所有高阶矩均定义为来自各子晶格占据概率的单点矩，因此目前路径概率方法近似中并未明确考虑自旋间的点间关联。不过，由于状态概率及相关量对于推导动态序参量至关重要，因此可以通过方程(2)–(3)并结合归一化条件，将它们用平均序参量表示出来。对于A层：
(4a) X1A=?124mA+14qA+16rA?116
(4b) X2A=18mA?14qA?12rA+916
(4c) X3A=98mA?14qA+12rA+916
(4d) X4A=124mA+14qA?16rA?116
对于B层则有以下表达式：
(5a) –X1B=‐57168mB+25923040qB+3711520rB?71152oB?1576pB+11440tB+15040uB?52048
(5b) X2B=495120mB?4994608qB?49911520rB+591152oB+592880pB?1288tB?1720uB+492048
(5c) –X3B=‐2453072mB+12992560qB+4331280rB?15128oB?564pB+1160tB+1240uB?2452048
(5d) X4B=12251024mB?18914608qB?18912304rB+831152oB+83576pB?1288tB?1144uB+12252048
(5e) –X5B=‐12251024mB?18914608qB+18912304rB+831152oB?83576pB?1288tB+1144uB+12252048
(5f) X6B=2453072mB+12992560qB?4331280rB?15128oB+564pB+1160tB?1240uB?2452048
(5g) –X7B=‐495120mB?4994608qB+49911520rB+591152oB?592880pB?1288tB+1720uB+492048
(5h) X8B=57168mB+25923040qB?3711520rB?71152oB+1576pB+11440tB?15040uB?52048

接下来，我们可以利用路径概率方法推导出平均序参量的动态方程。Kikuchi[36]提出的路径概率方法的核心方程描述了状态变量的变化率，其表达式为：
(6) dXidt=∑i≠jxji?xij。
其中xji表示系统从状态i转变为状态j的路径概率函数，而详细平衡条件要求xji=xij。在本研究中，我们将采用Kikuchi的第二种方法[36]，其表达式为：
(7) xij=kijZ?1Xiexp?βN?E?Xj。
在此公式中，β=1/kBT，kB为玻尔兹曼常数，在后续计算中我们统一取kB=1.0。kij为速率常数，且满足kij=kji。相比其他常见方法，路径概率方法的一个显著优势在于它具有三个速率常数，我们将其分别记为k1、k2和k3，它们实际上就是各种耦合常数。这些速率常数对于研究系统的动态性质至关重要。在Blume–Capel伊辛模型的固定晶格框架下，k1对应于同号（保持取向）的单点自旋转变，k2代表纯粹的自旋反转过程（+m ????m），而k3则描述了反转与磁化强度变化同时发生的转变过程。表I展示了A层和B层的速率常数。k1表示自旋粒子的平移运动，k2表示在特定位置上自旋粒子的旋转运动，而k3则表示自旋的平移与旋转同时发生。几何平均值k3可表示为k3=k1k2。层AX1A的构型为(3/2)X2A(1/2)X3A(?1/2)X4A(?3/2)X1A(3/2)，对应的参数关系为k1k3k2；X2A为(1/2)k1k2k3；X3A为(?1/2)k3k2k1；X4A为(?3/2)k2k3k1。层BX1B的构型为(7/2)X2B(5/2)X3B(3/2)X4B(1/2)X5B (?1/2)X6B (?3/2)X7B (?5/2)X8B (?7/2)X1B (7/2)，其参数关系分别为k1k1k1k3k3k3k2、X2B为k1k1k1k3k3k2k3、X3B为k1k1k1k3k2k3k3、X4B为k1k1k1k2k3k3k3、X5B为k3k3k3k2k1k1k1、X6B为k3k3k2k3k1k1k1、X7B为k3k2k3k3k1k1k1、X8B为k2k3k3k3k1k1k1。A层和B层的配分函数Z分别由下式给出：ZA=∑i=14exp?βN?E?XiA，ZB=∑j=18exp?βN?E?XjB，其中N为自旋粒子数，E为系统的内能，可通过公式(1)–(3)用XjA和XjB表示，即EN=?J{32X1A?X4A+12X2A?X3A72X1B+52X2B+32X3B+12X4B?12X5B?32X6B?52X7B?72X8B}?D{94X1A+X4A+14X2A+X3A+494X1B+254X2B+94X3B+14X4B+14X5B+94X6B+254X7B+494X8B}?H{32X1A?X4A+12X2A?X3A+72X1B+52X2B+32X3B+12X4B?12X5B?32X6B?52X7B?72X8B}。在混合自旋(3/2, 7/2)的Blume-Capel伊辛系统中，两个序参量mA和mB用于描述系统的动态磁性质。为确保方程数量与未知数数量相等，我们还定义并计算了其他序参量。为便于阐述，此处仅展示mA和mB的动态方程，其余方程见附录。这些动态方程可通过公式(4)–(9)得出。针对A层，平均磁化的动态方程为(10a)ΩdmAdξ={16cosha13k1+15k2?28k3sinhb+14k1?12k2?26k3coshbmA?sinha2?k2+k3sinhb+k1?k2coshbqA+23cosha?k1?3k2+4k3sinhb+2?k1+k3coshbrA+14sinha4k1?6k2+10k3sinhb+5k1+3k2+8k3coshb}/4kcoshacoshb；(11a)ΩdmBdξ={[?1225512k2+337126880k1+1436353760k3eb4cosha4+212918960k1?8699717920k3+245512k2eb3cosha3+?49512k2?7022510752k3+248755376k1eb2cosha2+264113840k1?682577680k3+5512k2eb1cosha1]mB+[1693311520k3?18912304k2?37395760k1eb4sinha4+55192304k3?37395760k1+38971280k2eb3sinha3+?37395760k1+22171280k3?24952304k2eb2sinha2+11332304k3+181311520k2?37395760k1eb1sinha1]qB+[12891152k3+18911152k2?301576k1eb4cosha4+1299640k2?301192k1+69071920k3eb3cosha3+4991152k2+279128k3?1505579k1eb2cosha2+2595760k2+213295760k3?2107579k1eb1cosha1]rB+[83576k2+19288k1?121576k3eb4sinha1+19288k1+367576k3?4564k2eb3sinha3+?3764k3+19288k1+295576k2eb2sinha2+11576k3+19288k1?49576k2eb1sinha1]oB+[2931440k3+61720k1?83288k2eb4cosha4+1532k2+61240k1?347480k3eb3cosha3+?732k3+61144k1?59288k2eb2cosha2+427720k1?8991440k3+7288k2eb1cosha1]pB+[?1360k1+7720k3?1144k2eb4sinha3+?1360k1?5144k3+380k2eb3sinha3+380k3?1360k1?5144k2eb2sinha2+?1144k3?1360k1+77202eb1sinha1]tB+[?5504k3?1252k1+7720k2eb4cosha4+?140k2+31840k3?184k1eb3cosha3+1168k3?5252k1+172k2eb2cosha2+?136k1+11360k3?1360k2eb1cosha1]uB+12251024k2?3631024k3+81512k1eb4sinha4+593512k1+26211024k3?7351024k2eb3sinha3+1105512k1+26651024k3+2451024k2eb2sinha2+1617512k1+39691024k3?351024k2eb1sinha1}/(2kcosha1eb1+cosha2eb2+cosha3eb3+cosha4eb4)，其中a=JmB+h/T，b=d/T，c=JmA+h/T，a1=72c，a2=52c，a3=32c，a4=12c，b1=494b，b2=254b，b3=94b，b4=14b，h=H/z，d=D/z，T=βz?1，z=4，Ω=ω/k，k=k1/k2，k3=k1k2，ξ=ωt。由于系统的动态磁性质主要由mA和mB决定，因此我们仅研究这两个量以分析系统的动态磁滞回线。A层和B层的动态磁滞回线面积定义为[47](12)A=?∮mA,BtdH=H0ω∮mA,Btsinωtdt。总磁化强度则定义为mT=mA+mB/2。由此可导出10个包含10个未知数的非线性微分方程，即公式(10a)、(11a)以及公式(A.10b)、(A.10c)、(A.11b)、(A.11c)、(A.11d)、(A.11e)、(A.11f)、(A.11g)。这些耦合的动态方程（PPM动态方程）是通过计算机代数程序MAPLE推导并简化得到的，而所有的数值时间积分则是通过我们自行编写的C++代码完成的。这些PPM动态方程（一组10个耦合的非线性常微分方程）是采用固定步长的预测-校正法进行数值积分的。最初几步使用标准的四阶龙格-库塔算法，之后的时间步长则通过亚当斯-莫尔顿预测-校正算法来确定。时间步长设定为Δt?=?11/112（无量纲单位），相当于振荡场相位变量的一个驱动周期（2π）内有64个均匀的时间步长。对于每组参数，系统会演化N_per?=?301个场周期；前N_tr?=?300个周期被视为瞬态阶段而被忽略，最终的磁滞量（回线面积、矫顽场和剩磁）则是从最后一个周期中计算得出的。通过减小时间步长以及延长忽略的瞬态阶段，我们验证了数值稳定的性，结果显示，在绘图精度范围内，磁滞回线并无变化。这种数值积分策略（亚当斯-莫尔顿预测-校正法结合罗姆伯格积分法计算回线面积）在我们之前的PPM动态磁滞研究中也有应用[43]、[44]、[45]、[46]、[47]、[48]、[49]、[50]。本研究的成果与讨论将在下一节中给出。3. 数值结果与讨论在本节中，我们将展示并分析通过PPM方法针对正交晶格上的混合自旋(3/2, 7/2) BCIS体系在振荡外磁场作用下的耦合非线性微分方程组所获得的数值结果。为避免过多图表给读者带来困扰，这里仅呈现最值得关注的结果。我们将把这些结果分为两类：一类是J?>?0的情况，另一类是J?5.85时，就不再观察到磁滞回线。这一趋势与关于动态自旋（1/2, 1）伊辛系统的理论研究[83]，以及利用熔融纺丝方法进行的实验研究[69]、[73]、[84]、[85]的结果一致。它也与Bai等人关于冷却速率对Fe53Nd37Al10合金磁性能影响的 research结果高度契合[86]。如前所述，熔融纺丝过程中的k2值对应着轮速，而轮速决定了材料的冷却速率。在高k2值下观察到宽大的磁滞回线，说明形成了硬磁体。这一结果与熔融纺丝工艺的实验结果一致，那些实验表明提高轮速能够提升合金的磁硬度[41]、[42]、[71]、[74]、[87]、[88]。下载：下载高分辨率图像（770KB）下载：下载完整尺寸图像图9。（彩色在线版）(a) 在固定h0=1.50、ω=6.00π、T=1.0、k1=1.0的前提下，不同速率常数（即快速凝固技术中的轮速）k2值下mA、mB和mT的DMHLB曲线。(b) 磁场矫顽力和剩磁强度随k2的变化规律。从某个特定的k2值开始直到k2=5.85，环状区域的面积基本保持不变，这意味着磁体的硬度也保持稳定。在Ni-Cr-Mo合金[41]以及通过熔融纺丝工艺制备的其他几种合金[69]、[73]、[74]、[84]、[85]、[86]中也有类似的实验观察结果。此外，mA、mB和mT的磁滞回线最初呈椭圆形。随着k2值的增大，环状区域先扩大，然后变为S形，随后又逐渐变窄。而且在某个特定的k2值之前，S形区域的面积还会进一步增大，接近矩形形态。这些结果与使用熔融纺丝技术制备的Nd-Fe-B合金[88]以及SmCoFeCu/NiZr7.5纳米复合材料[84]的实验研究结果相符。总之，如果我们考虑所有合理的k2值而非特定的数值范围，就会发现，与图9(a)的情况类似，但此时在较低的k2值下，较大的动态磁滞回线面积最初会减小到零，随后又随着k2的增大而再次上升。这一现象在图9(b)中有体现，该图展示了随k2变化的动态RF值和RM值。将图9(a)与理论计算结果进行比较后可以发现，它与之前使用PPM方法针对各种混合自旋伊辛系统所得到的DMHLB曲线[43]、[44]、[45]、[46]、[47]、[48]、[49]、[50]存在差异，而那些理论计算结果与该领域的其他理论和实验研究结果十分吻合。我们本文得到的这些新颖结果与现有实验研究结果之间的差异，可能源于熔融纺丝工艺中所用磁性合金中各组分浓度的不同。另外，也可能是因为无法实现更高的轮速，而轮速直接决定了合金带的厚度。4. 结论在本文中，我们运用路径概率法分析了在振荡磁场作用下，方形晶格上混合自旋（3/2, 7/2）Blume-Capel伊辛系统的动态磁滞特性。我们发现，该系统会展现出多种多样的磁滞行为，这些行为受温度、晶体场常数、振荡磁场的频率以及各种速率常数的影响，尤其是那些与自旋弛豫过程以及快速凝固过程中的冷却速率等物理过程相关的因素。该模型成功再现了实验中观察到的多重环结构、动态矫顽力以及剩磁行为，尤其是在铁磁性和反铁磁性相互作用存在的条件下。这些结果不仅证明了PPM方法在描述非平衡动力学过程中的有效性，还凸显了动力学参数在调控磁响应方面的作用。未来的研究可以将这些探索拓展到三维系统以及其他晶格结构中，从而进一步深入理解复杂的自旋动力学现象。作者贡献声明Mustafa Gen?aslan：撰写——审稿与编辑、指导、项目管理、方法设计、概念构思。Beyza Ta?：软件应用、方法设计、实验研究。未引用的参考文献[33]、[52]、[54]