使用原子尺度上的随机LLG方程与3TM 耦合来描述由超快飞秒激光脉冲如何激发Fe 薄膜产生太赫兹波的过程。

在介观尺度上，当铁磁材料系统的温度为0 K时，内部的磁矩倾向于围绕外磁场的方向进动，进动过程遵循LLG方程（1）：

其中$\gamma $是旋磁比，$\alpha $是阻尼系数，$m$是三维空间中的单位磁化强度，$ H_{\mathrm{eff}} $是等效外场，如式（2）：

其中$ E $是磁性材料体系的总能量，有效场$ H\mathrm{_{eff}} $中的第一项描述了相互作用的影响，第二项$ H\mathrm{_{therm}} $描述了热效应对磁化动力学的影响, ${\mu _0}$是真空磁导率。热扰动场$ H_{\mathrm{therm}} $如式（3）所定义：

其中$ k_{\mathrm{B}} $是玻尔兹曼常数, T是开尔文温度，ΔV是网格单元的体积, $ \eta $是三维空间中的矢量，其各个分量遵循标准正态分布。

铁磁材料中电子、晶格和自旋三个子系统之间的相互作用可用三温度模型(3TM)来解释^[13]. 假设每个子系统都处于热平衡状态，并且其温度独立。声子系统和自旋系统的热演化将受到热电子的影响，三个系统之间的热力学过程可用微分方程（4），（5）和（6）来描述：

其中T_e, T_l和T_s, 分别表示电子系统、声子系统和自旋系统的温度。相应系统的热容由C_e, C_l和C_s表征，G_el, G_es和G_sl 是三个系统之间的交换系数, P(t)代表激光源产生的热通量。

基于微磁学模拟方法^[14]模拟了飞秒激光诱导Fe薄膜产生太赫兹波的过程。相关参数为：交换相互作用系数A为1.9545 × 10⁻¹¹ Am，面内单轴磁晶各向异性常数为4.4675 × 10⁻⁴ A/m，饱和磁化强度M_s为1.71 × 10⁶ A/m，旋磁比γ为2.21 × 10⁵ Hz/(A/m)，阻尼系数α为0.01。脉冲激光能量为2 Jm⁻²，环境温度为300 K, C_e = 1.4 × 10⁶ Jm⁻³K⁻¹, C_l = 8.8 × 10⁶ Jm⁻³K⁻¹, C_s = 4 × 10⁵ Jm⁻³K⁻¹, G_el = 8 × 10¹⁷ Wm⁻³K⁻¹, G_es = 1.8 × 10¹⁸ Wm⁻³K⁻¹, G_sl = 0.8 × 10¹⁷ Wm⁻³K⁻¹。

自旋波是磁体中的一种元激发，由于电子自旋之间的交换作用而产生，表现为波的形式。在纳米尺度上，物质的铁磁性受短程的交换相互作用比长程的磁偶极相互作用强得多。因此，在接下来的模拟过程中忽略了磁偶极相互作用对磁矩的影响，偶极场仅在对外辐射的过程中考虑。

考虑若干磁矩排成一列，相邻磁矩之间通过交换相互作用联系在一起。假设磁矩沿$ x $方向排列，磁矩的初始磁化方向为$ {\textit{z}} $方向，同时存在一个沿着$ {\textit{z}} $方向的单轴磁晶各向异性。此时，等效场可以写成$ H_{\mathrm{eff}}\text{ = }A\nabla^{\text{2}}m\text{ + }Km_{\text{z}}e_{\text{z}}\text{ + }H $。由于相邻的磁矩之间存在耦合，磁矩的进动会存在相位差，从而体现出波动性质，因此除了需要考虑进动频率外还要考虑沿一维结构传播的波矢，可以得到式（7）：

这里假定了交变的外加磁场相对于内场很小，可以将等效外场拆分为静态和动态两部分。同理，假设磁矩的进动幅度非常小，做同样的线性化处理。其中$ \omega $为磁场分量的角频率，同时也是磁矩进动的角频率。

由式（7）分量的形式和式（1）可以得到式（8），其中$ {{\boldsymbol{m}}}_{0} $和 $ {{\boldsymbol{h}}}_{0} $均是沿着${\textit{z}} $方向的矢量：

类似于久期方程，要想使磁矩的变化量有解，需要使行列式（9）为 0：

最终得到交换自旋波的色散关系，即频谱，如式（10）：

考虑一个沿着$x$方向长度为 500 nm的准一维模型，初始磁矩沿${\textit{z}} $方向饱和，在最左端沿$y$方向施加一个持续半个周期的交变外磁场（磁场成正弦函数变化，周期为 2 GHz），外磁场峰值为 10 mT。系统最左端的磁矩受到激励后，被激励点的磁矩将首先发生扰动，随后进动过程将沿$x$方向传递，生成如图1所示的沿着$x$方向传播的自旋波。自旋波的传递随着距离逐步衰减，即磁矩相对于平衡位置的偏转将沿着传播方向逐步减小。

图2 给出了一维结构在受到外界交变磁场激励后的模拟结果。在 500 nm 的一维结构上，单位磁化强度的$y$ 分量变化情况如图2（a）所示，横坐标为$x$方向上的不同位置，纵坐标为激励后的时间，颜色代表强度，可以清晰地看出自旋波的传递过程。$x$方向的磁矩变化与$y$方向变化趋势一致，只存在一定的相位延迟，这里不再给出$x$方向分量的变化情况。通过对图2（a）进行傅里叶变换，得到交换自旋波的色散关系$\omega \sim k$。图2（b）中的背景即是对图2（a）的数据进行傅里叶变换后的结果，图2（b）中红色的线为理论值。对比发现，大多数自旋波的模式与理论值符合很好，但也存在一些其他模式的自旋波，这与激励的外场和Fe的参数相关。

可以看出，交换自旋波的色散关系呈抛物线，激发频率越高，所激发出的自旋波的波长越短。与此同时，由于磁晶各向异性的存在，自旋波存在频率阈值，阈值为$\gamma K$, 即只有当外场激发的频率大于某个阈值时才能够激发出自旋波，模拟结果与解析结果符合很好。

在实际磁性体系中，交换相互作用使相邻磁矩平行排列从而造成自发磁化，磁晶各向异性使得所有磁矩都保持在易轴方向，然而仅靠这两个相互作用并不能直接导致磁畴的产生。大尺度上丰富的磁畴结构往往起源于磁偶极相互作用，对于某一磁化单元，它将在周围产生与自身磁化方向相反的磁场，促使周围的单元沿相反方向磁化，为了使包括磁偶极相互作用在内的总能量达到最低，磁畴的出现才成为必要。

磁偶极相互作用主要在较大尺度上导致磁畴壁的产生，但磁畴壁的具体形貌在小尺度上主要受交换相互作用的影响，下面将讨论不考虑磁偶极相互作用，仅考虑交换相互作用和单轴磁晶各向异性时磁畴壁的磁化强度分布。首先将直角坐标系下的 LLG 方程（1）应用到球坐标系。在球坐标系中，由于各个单元的磁化强度大小不变，变量减少为$\theta $和 $\phi $，单位磁化强度矢量${\boldsymbol{m}}$为：

为简单起见，这里只考虑一维情况，此时体系的总能量在球坐标系中为式（12）或式（13）。

要使整个系统的能量达到最小，泛函式（13）需取最小值。系统总能量的泛函包含多个一元函数，根据基本变分原理，求系统能量最小时等价于求解如下欧拉-拉格朗日方程组式（14）：

通过求解系统的拉格朗日量对$\phi $的变分可以得到：

要使上式对于任意$\theta $均成立，$C$应等于0，即$ {{\partial \phi } \mathord{\left/ {\vphantom {{\partial \phi } {\partial x}}} \right. } {\partial x}} = 0 $，此时$\phi $为任意常数。将$ {{\partial \phi } \mathord{\left/ {\vphantom {{\partial \phi } {\partial x}}} \right. } {\partial x}} = 0 $带入对$\theta $的变分式中得到式（16）。

进一步可化简为式（17）。

由于这里主要考虑磁畴壁的磁化强度分布，这是一个稳态结构，故对于磁化强度分布而言只与空间位置有关，与时间无关，式（17）中的偏微分可以转化为全微分，两端同乘$ \mathrm{d}\theta\mathord{\left/\vphantom{d\theta\mathrm{d}x}\right.}\mathrm{d}x $，再求原函数可以得到式（18），

进行整理可以得到磁畴壁磁化强度分布关于位置的表达式（19），

显然磁畴壁的宽度与参数$ \delta $紧密相关，$ \delta = \sqrt {{A \mathord{\left/ {\vphantom {A K}} \right. } K}} $，体现了磁畴壁是不同能量之间相互竞争的结果，该结果对于一维磁畴壁是普遍适用的。

最终得到不同$ \phi $对应的磁畴壁结构类型，$ \phi = 0 $时对应奈尔壁，如图3（a）, $ \phi = {\pi \mathord{\left/ {\vphantom {\pi 2}} \right. } 2} $时对应布洛赫壁，如图3（b）。

对于二维情况，采用与对拉格朗日量式（12）求解变分方程组的方法, 此时泛函变分为式（20）。

简单起见，在这里不对式（20）进行求解，只研究二维磁畴壁的形貌特征。数值模拟的结构如图4所示，图4中大圆柱面内的小圆盘是半径为 100 nm、厚度为 8 nm 的Fe 薄膜，薄膜的厚度很薄，近似认为是一个二维模型，其受到交换相互作用、$ {\textit{z}} $方向的单轴磁晶各项异性和磁偶极相互作用。 Fe 薄膜外部的灰色部分是真空，真空的设置可以确保薄膜中的磁偶极作用计算准确。

初始条件为薄膜内的各单元磁矩随机分布，通过求解麦克斯韦方程组和LLG 方程，得到了稳态下薄膜内部的磁矩分布情况。如图5所示，颜色表示$ {\textit{z}} $方向上的磁矩，$ x $和 $ y $方向上的磁矩用黑色箭头表示。在磁偶极相互作用和单轴磁晶各向异性的相互竞争下，薄膜上的磁矩形成了涡旋状结构，这类涡旋结构广泛存在于圆盘状的磁性样品中。由于圆盘具有较高的对称性，圆盘内的磁矩倾向于顺时针或者逆时针排列，在这种构型下体系的退磁场能被尽可能地降到最低。由于磁偶极相互作用或者面内各向异性的作用，圆盘周边的磁矩倾向于排列于面内，然而这种排列方式会产生一个奇点，即中心位置。该点的磁矩由于对称性的原因只能垂直于薄膜表面向上或向下分布，涡旋核与圆盘边缘之间的磁矩由于交换相互作用而缓慢变化。

接下来研究了温度场对磁矩运动的影响，并通过统计与理论上的磁矩分布对比。根据统计力学，单个粒子和一个无限大热库所组成的体系可以看作是一个巨正则系综。在系统达到平衡态后，系统处在某个能量状态$ E $的概率将服从玻尔兹曼分布式(21)。

其中$ \beta=1\mathord{\left/\vphantom{1k_BT}\right.}k_{\mathrm{B}}T $。假设磁矩沿$ {\textit{z}} $方向磁化且仅受外磁场作用，那么体系的塞曼能为$ E_Z=(1-m_{\text{z}})\mu_0M\mathrm{_S}VH $。为了测试模拟方法的有效性，首先在三维空间中计算了5000步时间步长，通过并行计算同时计算14个单元，共计70000组数据，并对磁矩的能量进行统计，结果如图6所示。

对于同一个磁矩，时间是没有影响的，可以认为不同时间的磁矩指向是完全独立的。可以看出磁矩分布近似成指数分布，这里需要注意的是，磁矩为1的时候能量最小。为了更进一步研究磁矩统计分布情况，对式（21）两边同时取对数，得到$ \ln P = - \beta E + C $。通过取对数，可以更清晰地对比出理论与模拟之间的关系。根据模拟结果统计出 $ \ln P $与理论值$ - \beta E + C $相对比，得到图7，图7中蓝色圆点为模拟结果统计，红色直线为理论值。

由图7可以看出在能量较低时统计结果和理论符合得非常好，但随着能量增大理论值与统计值出现了较小的偏差，这主要是由两个原因造成的：一是磁矩被激发到高能级的概率较低，导致高能量时数据稀少。二是能量越高越偏离平衡态，会导致系统的能量分布偏离玻尔兹曼分布。

前面我们用微磁学模拟方法模拟了一维自旋波的色散关系、磁畴壁和热平衡状态下的磁矩分布。自旋波和磁畴壁是飞秒激光作用于Fe薄膜内部后，在其内可能存在的微磁学现象，热效应的等效场更是连接了三温度模型和LLG方程，这为模拟飞秒激光与Fe薄膜的作用提供了支持。为了简化模型，只考虑在非磁性金属层上生长的磁性Fe薄膜。主要考虑飞秒激光作用在Fe薄膜上产生磁矩扰动的过程，从而产生纯自旋流。在模拟过程中，忽略了飞秒激光对可能产生自旋极化电流并辐射电磁波的铁磁材料的电子分布密度的影响。由于Fe薄膜（10 nm）的线性度远小于太赫兹波长（约10⁻⁴ m），因此整个Fe薄膜可以近似为点源，相当于Fe薄膜的平均磁矩辐射。所以可以使用随机LLG方程模拟Fe薄膜中的磁动力学过程，然后通过经典的磁偶极子辐射理论获得了飞秒激光诱导的太赫兹辐射。

图8为模拟得到的太赫兹辐射图谱。图8（a）为Fe薄膜辐射的太赫兹波的电场强度E_y，这与实验结果一致^[3-4,15]。由图8（a）可知，电场强度E_y在0.3 ps 时达到了最大值。通过对太赫兹辐射的电场强度的时域关系图8（a）进行傅里叶变换，得到了辐射电场的频谱图8（b）。可以看出辐射的电场强度的$ y $方向的频谱连续性较好，这意味着辐射的太赫兹波具有一个良好信号的特征，没有突变或间断。它能够提供更多的信息、更高的数据传输速率和更好的调制灵活性，可以用作高传输效率的通讯媒介。

通过基于 LLG 方程的微磁学模拟研究了微磁学领域的一维自旋波的色散关系、磁畴壁和热效应下的磁矩统计分布，并将模拟结果与理论进行对比验证，证实了微磁学模拟的可行性和准确性。这些微磁学常见现象也是飞秒激光作用Fe薄膜后其内部可能存在的微磁学现象，热效应的等效场更是连接了三温度模型和LLG方程，这些模拟过程和结果为模拟飞秒激光与 Fe 薄膜的作用打下了基础。由模拟的飞秒激光与Fe薄膜的作用结果可以看出辐射的电场强度的$ y $方向的频谱连续性较好，具有一个良好信号的特征，能够提供更多的信息、更高的数据传输速率和更好的调制灵活性，可以用作高传输效率的通讯媒介。