Process of Generating Terahertz Waves in Fe Thin Film Simulated by Using Micromagnetic Simulation
-
摘要:
太赫兹辐射的独特特性使其成为一种多功能工具，被应用于成像、光谱学、通信、材料科学和生物医学等各个领域。通过基于Landau Lifshitz Gilbert方程的微磁学模拟研究了微磁学领域中的一维自旋波的色散关系、磁畴壁和热效应下的磁矩统计分布，在此基础上又结合三温度模型模拟了飞秒激光作用于Fe 薄膜后产生太赫兹波的过程。模拟产生的太赫兹波的电场强度的分量Ey在0.3 ps时达到了最大值, 由辐射电场的频谱可以看出辐射的电场强度的Ey分量的频谱连续性较好，没有突变或间断。
Abstract:The unique characteristics of terahertz radiation make it a multifunctional tool that can be applied in various fields such as imaging, spectroscopy, communication, materials science, and biomedical applications. In this paper, the dispersion relationship of one-dimensional spin waves, magnetic domain walls, and statistical distribution of magnetic moments under thermal effects in the field of micromagnetics are studied through micromagnetic simulations based on the Landau-Lifshitz-Gilbert equation. Based on these, the three-temperature model is combined to simulate the process of terahertz waves generated by femtosecond laser irradiation on Fe thin film. The component of the electric field intensity of the simulated terahertz wave reached its maximum value at 0.3 ps. From the spectrum of the radiated electric field, it can be seen that the spectral continuity in the y-direction of the radiated electric field intensity is good, without any abrupt changes or interruptions.
-
Keywords:
- Femtosecond laser /
- Micromagnetism /
- Terahertz /
- LLG equation
-
太赫兹 (terahertz, THz) 波是属于微波和红外辐射之间的电磁波谱，频率范围约为 0.1 至 10 THz（波长约为 30 μm 至 3 mm）。 该区域通常被称为“太赫兹间隙”，因为与电磁频谱的其他部分相比，生成、操纵和检测太赫兹辐射历来具有挑战性。然而，最近技术的进步使得研究者们对研究太赫兹波的兴趣大增，并在利用太赫兹辐射进行各种应用方面取得了重大进展[1]。太赫兹辐射的独特特性使其成为了一种多功能工具，可应用于成像、光谱学、通信、材料科学和生物医学等各个领域，并且持续的研究和技术进步在不断扩大太赫兹技术的能力和潜在应用[1-2]。
设计和生产先进的太赫兹光源的传统方法存在一些如成本、带宽等方面的缺点[3]。然而，基于与铁磁材料相关的超快自旋动力学产生太赫兹辐射的方法具有低成本、宽频带、高效率、易于集成等优点。因此，为了得到优良的太赫兹源，很多研究者聚焦于研究与铁磁材料相关的自旋太赫兹源[1]。如研究人员使用飞秒脉冲激光激发Fe、Co、Ni单层薄膜[4]、铁磁性/非铁磁性异质结构[5-6]和其他材料异质结结构[7-8]以激发太赫兹波。由此看出，用飞秒脉冲激光激发薄膜材料是太赫兹波产生的一个非常重要的手段。为了更深刻地理解脉冲激光诱导太赫兹波的产生机制，研究人员还建立了重要的理论模型和计算模拟[9-11]来解释其原因。对于飞秒脉冲激光作用于磁性薄膜后，磁性薄膜磁化动力学的理论研究大多是基于唯象模型。所报道的方法包括3TM(three temperature model)、M3TM(microcosmic three temperature model)、使用传统LLG(Landau-Lifshitz-Gilbert)方程的微磁学模拟方法、使用随机LLG方程的微磁学模拟方法和使用LLB(Landau-Lifshitz-Bloch)方程的微磁学模拟方法等。
由于传统LLG方程中的等效磁场和Gilbert阻尼项α不会改变磁化强度的振幅，振幅是常数。因此，当温度接近或高于居里温度时，LLG方程不能用于描述超快退磁和磁化过程的演变。它只能用于研究温度为0 K时的磁化过程。而LLB方程是解决此类问题最通用的方法，因为LLB方程只使用热力学自洽宏观方程，LLB方程式中的横向和纵向阻尼参数只与系统和热储层之间的耦合参数以及自旋系统的伪温度（由于自旋系统仍处于非平衡状态，因此无法定义其真实温度）和居里温度有关。由于热效应，磁化强度的大小不再是常数，而是随时间变化。LLB方程不仅考虑了LLG方程的横向Gilbert阻尼效应，还考虑了纵向阻尼效应。
另外，如果把一个形式化的朗之万随机场加到LLG方程中的等效磁场中，即随机LLG方程。使用随机LLG方程可用于研究居里温度附近的磁化动力学，如系统在很短的时间内完成退磁和磁化恢复的过程。将随机LLG方程应用于原子（即原子尺度上的网格划分）尺度，即用等效交换磁场表示海森堡短程相互交换耦合作用，也可以用来研究线性不超过几十纳米的系统的磁化动力学。
尽管研究人员用唯象理论在解释实验结果方面取得了一些成功，但对铁磁薄膜和铁磁/非铁磁异质结构中产生太赫兹波方面的微观机制仍存在争议，还有许多问题值得研究[12]，因此，作者采用微磁学模拟方法，研究了微磁学领域中的一维自旋波的色散关系、磁畴壁和热效应下的磁矩统计分布，在此基础上结合三温度模型（3TM）和原子尺度的随机LLG方程，研究了飞秒激光作用于Fe 薄膜后产生太赫兹波的过程。期望通过模拟为实验提供理论支撑，促进高性能太赫兹器件的发展。
1. 模拟方法
使用原子尺度上的随机LLG方程与3TM 耦合来描述由超快飞秒激光脉冲如何激发Fe 薄膜产生太赫兹波的过程。
1.1 LLG方程和有效磁场
在介观尺度上，当铁磁材料系统的温度为0 K时，内部的磁矩倾向于围绕外磁场的方向进动，进动过程遵循LLG方程（1）：
$$ \frac{\partial\boldsymbol{m}}{\partial t}=-\gamma\boldsymbol{m}\times\boldsymbol{H}\mathrm{_{eff}}+\alpha\boldsymbol{m}\times\frac{\partial\boldsymbol{m}}{\partial t} $$ (1) 其中$\gamma $是旋磁比，$\alpha $是阻尼系数，$m$是三维空间中的单位磁化强度，$ H_{\mathrm{eff}} $是等效外场，如式（2）：
$$ \begin{array}{c}{{\boldsymbol{H}}}_{{\mathrm{eff}}}=-\dfrac{1}{{\mu }_{0}{M}_{{\mathrm{s}}}}\dfrac{\delta E}{\delta {\boldsymbol{m}}}\end{array}\text{ + }{{\boldsymbol{H}}}_{{\mathrm{therm}}} $$ (2) 其中$ E $是磁性材料体系的总能量，有效场$ H\mathrm{_{eff}} $中的第一项描述了相互作用的影响，第二项$ H\mathrm{_{therm}} $描述了热效应对磁化动力学的影响, ${\mu _0}$是真空磁导率。热扰动场$ H_{\mathrm{therm}} $如式（3）所定义：
$$ {{\boldsymbol{H}}_{{\mathrm{therm}}}} = \eta \sqrt {\frac{{2{k_{\mathrm{B}}}T\alpha }}{{\gamma {\mu _0}{M_{\mathrm{s}}}\Delta V\Delta t}}} $$ (3) 其中$ k_{\mathrm{B}} $是玻尔兹曼常数, T是开尔文温度，ΔV是网格单元的体积, $ \eta $是三维空间中的矢量，其各个分量遵循标准正态分布。
1.2 三温度模型
铁磁材料中电子、晶格和自旋三个子系统之间的相互作用可用三温度模型(3TM)来解释[13]. 假设每个子系统都处于热平衡状态，并且其温度独立。声子系统和自旋系统的热演化将受到热电子的影响，三个系统之间的热力学过程可用微分方程（4），（5）和（6）来描述：
$$ {C}_{{\mathrm{e}}}\frac{{\mathrm{d}}{T}_{{\mathrm{e}}}}{{\mathrm{d}}t}=-{G}_{{\mathrm{el}}}\left({T}_{{\mathrm{e}}}-{T}_{{\mathrm{l}}}\right)-{G}_{{\mathrm{es}}}\left({T}_{{\mathrm{e}}}-{T}_{{\mathrm{s}}}\right)+P\left(t\right) $$ (4) $$ {C}_{{\mathrm{l}}}\frac{{\mathrm{d}}{T}_{{\mathrm{l}}}}{{\mathrm{d}}t}={G}_{{\mathrm{el}}}\left({T}_{{\mathrm{e}}}-{T}_{{\mathrm{l}}}\right)-{G}_{{\mathrm{sl}}}\left({T}_{{\mathrm{l}}}-{T}_{{\mathrm{s}}}\right) $$ (5) $$ {C_{\mathrm{s}}}\frac{{{\mathrm{d}}{T_{\mathrm{s}}}}}{{{\mathrm{d}}t}} = {G_{{\mathrm{es}}}}\left( {{T_{\mathrm{e}}} - {T_{\mathrm{s}}}} \right) - {G_{{\mathrm{sl}}}}\left( {{T_{\mathrm{l}}} - {T_{\mathrm{s}}}} \right) $$ (6) 其中Te, Tl和Ts, 分别表示电子系统、声子系统和自旋系统的温度。相应系统的热容由Ce, Cl和Cs表征，Gel, Ges和Gsl 是三个系统之间的交换系数, P(t)代表激光源产生的热通量。
1.3 模拟参数
基于微磁学模拟方法[14]模拟了飞秒激光诱导Fe薄膜产生太赫兹波的过程。相关参数为：交换相互作用系数A为1.9545 × 10−11 Am，面内单轴磁晶各向异性常数为4.4675 × 10−4 A/m，饱和磁化强度Ms为1.71 × 106 A/m，旋磁比γ为2.21 × 105 Hz/(A/m)，阻尼系数α为0.01。脉冲激光能量为2 Jm−2，环境温度为300 K, Ce = 1.4 × 106 Jm−3K−1, Cl = 8.8 × 106 Jm−3K−1, Cs = 4 × 105 Jm−3K−1, Gel = 8 × 1017 Wm−3K−1, Ges = 1.8 × 1018 Wm−3K−1, Gsl = 0.8 × 1017 Wm−3K−1。
2. 微磁学现象模拟和Fe薄膜中THz信号模拟结果
2.1 一维自旋波及其色散关系
自旋波是磁体中的一种元激发，由于电子自旋之间的交换作用而产生，表现为波的形式。在纳米尺度上，物质的铁磁性受短程的交换相互作用比长程的磁偶极相互作用强得多。因此，在接下来的模拟过程中忽略了磁偶极相互作用对磁矩的影响，偶极场仅在对外辐射的过程中考虑。
2.1.1 一维自旋波的解析解
考虑若干磁矩排成一列，相邻磁矩之间通过交换相互作用联系在一起。假设磁矩沿$ x $方向排列，磁矩的初始磁化方向为$ {\textit{z}} $方向，同时存在一个沿着$ {\textit{z}} $方向的单轴磁晶各向异性。此时，等效场可以写成$ H_{\mathrm{eff}}\text{ = }A\nabla^{\text{2}}m\text{ + }Km_{\text{z}}e_{\text{z}}\text{ + }H $。由于相邻的磁矩之间存在耦合，磁矩的进动会存在相位差，从而体现出波动性质，因此除了需要考虑进动频率外还要考虑沿一维结构传播的波矢，可以得到式（7）：
$$ \left\{ \begin{array}{l}{\boldsymbol{m}} = {{\boldsymbol{m}}}_{0}+ \delta {\boldsymbol{m}}{e}^{i\left(\omega t+kx\right)}\\ {H}_{{\mathrm{eff}}} = {{\boldsymbol{h}}}_{0}+ \delta {\boldsymbol{h}}{e}^{i\omega t}\end{array} \right.$$ (7) 这里假定了交变的外加磁场相对于内场很小，可以将等效外场拆分为静态和动态两部分。同理，假设磁矩的进动幅度非常小，做同样的线性化处理。其中$ \omega $为磁场分量的角频率，同时也是磁矩进动的角频率。
由式（7）分量的形式和式（1）可以得到式（8），其中$ {{\boldsymbol{m}}}_{0} $和 $ {{\boldsymbol{h}}}_{0} $均是沿着${\textit{z}} $方向的矢量：
$$ \left\{ \begin{gathered} i\omega \delta {m_x} + \gamma K\delta {m_y} + \gamma A{k^2}\delta {m_y} = \gamma \delta {h_{\textit{z}} } \\ i\omega \delta {m_y} - \gamma K\delta {m_x} - \gamma A{k^2}\delta {m_x} = - \gamma \delta {h_x} \\ i\omega \delta {m_{\textit{z}} } = 0 \\ \end{gathered} \right. $$ (8) 类似于久期方程，要想使磁矩的变化量有解，需要使行列式（9）为 0：
$$ \left| {\begin{array}{*{20}{c}} { - \gamma \left( {A{k^2} + K} \right)}&{i\omega } \\ { - i\omega }&{ - \gamma \left( {A{k^2} + K} \right)} \end{array}} \right| = 0 $$ (9) 最终得到交换自旋波的色散关系，即频谱，如式（10）：
$$ \omega = \gamma (A{k^2} + K) $$ (10) 2.1.2 一维自旋波的数值模拟结果
考虑一个沿着$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的参数相关。
![]()
图 2 磁化强度 y 分量模拟结果。（a）单位磁化强度 y 分量变化情况，（b）自旋波色散关系Figure 2. Simulation results of y component of magnetization. (a) Variation of y component of unit magnetization, (b) spin wave dispersion relationship可以看出，交换自旋波的色散关系呈抛物线，激发频率越高，所激发出的自旋波的波长越短。与此同时，由于磁晶各向异性的存在，自旋波存在频率阈值，阈值为$\gamma K$, 即只有当外场激发的频率大于某个阈值时才能够激发出自旋波，模拟结果与解析结果符合很好。
2.2 磁畴结构
在实际磁性体系中，交换相互作用使相邻磁矩平行排列从而造成自发磁化，磁晶各向异性使得所有磁矩都保持在易轴方向，然而仅靠这两个相互作用并不能直接导致磁畴的产生。大尺度上丰富的磁畴结构往往起源于磁偶极相互作用，对于某一磁化单元，它将在周围产生与自身磁化方向相反的磁场，促使周围的单元沿相反方向磁化，为了使包括磁偶极相互作用在内的总能量达到最低，磁畴的出现才成为必要。
2.2.1 一维磁畴壁
磁偶极相互作用主要在较大尺度上导致磁畴壁的产生，但磁畴壁的具体形貌在小尺度上主要受交换相互作用的影响，下面将讨论不考虑磁偶极相互作用，仅考虑交换相互作用和单轴磁晶各向异性时磁畴壁的磁化强度分布。首先将直角坐标系下的 LLG 方程（1）应用到球坐标系。在球坐标系中，由于各个单元的磁化强度大小不变，变量减少为$\theta $和 $\phi $，单位磁化强度矢量${\boldsymbol{m}}$为：
$$ {\boldsymbol{m}} = \left( {\begin{array}{*{20}{c}} {\sin \theta \,\, \cos \phi } \\ {\sin \theta \,\, \sin \phi } \\ {\cos\theta } \end{array}} \right) $$ (11) 为简单起见，这里只考虑一维情况，此时体系的总能量在球坐标系中为式（12）或式（13）。
$$ {E_{{\mathrm{tot}}}} = \frac{{{\mu _0}{{\boldsymbol{M}}_{\mathrm{s}}}}}{2}\smallint \left[ {A{{\left( {\frac{{\partial \theta }}{{\partial x}}} \right)}^2} + A{\text{si}}{{\text{n}}^2}\theta {{\left( {\frac{{\partial \phi }}{{\partial x}}} \right)}^2} + K{\text{si}}{{\text{n}}^2}\theta } \right]{\mathrm{d}}x $$ (12) $$ = \smallint L\left( {x;\theta \left( x \right),\theta '\left( x \right);\phi '\left( x \right)} \right){\mathrm{d}}x $$ (13) 要使整个系统的能量达到最小，泛函式（13）需取最小值。系统总能量的泛函包含多个一元函数，根据基本变分原理，求系统能量最小时等价于求解如下欧拉-拉格朗日方程组式（14）：
$$ \left\{\begin{array}{*{20}{c}}\dfrac{\partial L}{\partial\phi}-\dfrac{\mathrm{d}}{\mathrm{d}x}\dfrac{\partial L}{\partial\phi'}=0 \\ \dfrac{\partial L}{\partial\theta}-\dfrac{\mathrm{d}}{\mathrm{d}x}\dfrac{\partial L}{\partial\theta'}=0\end{array}\right. $$ (14) 通过求解系统的拉格朗日量对$\phi $的变分可以得到：
$$ {\text{2}}A{\sin ^2}\theta \frac{{\partial \phi }}{{\partial x}} = C $$ (15) 要使上式对于任意$\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）。
$$ \frac{\partial\left(K\sin^2\theta\right)}{\partial\theta}-\frac{\mathrm{d}}{\mathrm{d}x}\frac{2A\theta'\partial\theta'}{\partial\theta'}=0 $$ (16) 进一步可化简为式（17）。
$$ \frac{{\partial \left( {K{{\sin }^2}\theta } \right)}}{{\partial \theta }} = A{\left( {\frac{{\partial \theta }}{{\partial x}}} \right)^2} $$ (17) 由于这里主要考虑磁畴壁的磁化强度分布，这是一个稳态结构，故对于磁化强度分布而言只与空间位置有关，与时间无关，式（17）中的偏微分可以转化为全微分，两端同乘$ \mathrm{d}\theta\mathord{\left/\vphantom{d\theta\mathrm{d}x}\right.}\mathrm{d}x $，再求原函数可以得到式（18），
$$ x = \sqrt {\frac{A}{K}} {\text{ln tan}}\frac{\theta }{2} $$ (18) 进行整理可以得到磁畴壁磁化强度分布关于位置的表达式（19），
$$ \theta = 2{\text{arc}}\,\, \tan \theta \left( {{e^{\tfrac{{x - {x_0}}}{\delta }}}} \right) $$ (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）。
![]()
图 3 不同$ \phi $对应的磁畴壁结构。（a）奈尔壁结构图， （b）布洛赫壁结构图Figure 3. Magnetic domain walls corresponding to different $ \phi $. (a) Neel wall structure, (b) Bloch wall structure2.2.2 二维磁畴壁
对于二维情况，采用与对拉格朗日量式（12）求解变分方程组的方法, 此时泛函变分为式（20）。
$$ \left\{ {\begin{array}{*{20}{c}} {\dfrac{{\partial L}}{{\partial \phi }} - \dfrac{\partial }{{\partial x}}\dfrac{{\partial L}}{{\partial {\phi _x}}} - \dfrac{\partial }{{\partial y}}\dfrac{{\partial L}}{{\partial {\phi _y}}} = 0} \\ {\dfrac{{\partial L}}{{\partial \theta }} - \dfrac{\partial }{{\partial x}}\dfrac{{\partial L}}{{\partial {\theta _x}}} - \dfrac{\partial }{{\partial y}}\dfrac{{\partial L}}{{\partial {\theta _y}}} = 0} \end{array}} \right. $$ (20) 简单起见，在这里不对式（20）进行求解，只研究二维磁畴壁的形貌特征。数值模拟的结构如图4所示，图4中大圆柱面内的小圆盘是半径为 100 nm、厚度为 8 nm 的Fe 薄膜，薄膜的厚度很薄，近似认为是一个二维模型，其受到交换相互作用、$ {\textit{z}} $方向的单轴磁晶各项异性和磁偶极相互作用。 Fe 薄膜外部的灰色部分是真空，真空的设置可以确保薄膜中的磁偶极作用计算准确。
![]()
图 4 二维磁畴壁模拟结构示意图Figure 4. Schematic diagram of two-dimensional magnetic domain wall simulation structure初始条件为薄膜内的各单元磁矩随机分布，通过求解麦克斯韦方程组和LLG 方程，得到了稳态下薄膜内部的磁矩分布情况。如图5所示，颜色表示$ {\textit{z}} $方向上的磁矩，$ x $和 $ y $方向上的磁矩用黑色箭头表示。在磁偶极相互作用和单轴磁晶各向异性的相互竞争下，薄膜上的磁矩形成了涡旋状结构，这类涡旋结构广泛存在于圆盘状的磁性样品中。由于圆盘具有较高的对称性，圆盘内的磁矩倾向于顺时针或者逆时针排列，在这种构型下体系的退磁场能被尽可能地降到最低。由于磁偶极相互作用或者面内各向异性的作用，圆盘周边的磁矩倾向于排列于面内，然而这种排列方式会产生一个奇点，即中心位置。该点的磁矩由于对称性的原因只能垂直于薄膜表面向上或向下分布，涡旋核与圆盘边缘之间的磁矩由于交换相互作用而缓慢变化。
2.3 热平衡状态下的磁矩分布
接下来研究了温度场对磁矩运动的影响，并通过统计与理论上的磁矩分布对比。根据统计力学，单个粒子和一个无限大热库所组成的体系可以看作是一个巨正则系综。在系统达到平衡态后，系统处在某个能量状态$ E $的概率将服从玻尔兹曼分布式(21)。
$$ P\left( E \right) \propto {e^{ - \beta 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所示。
![]()
图 6 热效应下的磁矩变化和统计。（a）磁矩 z 分量随时间的变化，（b）磁矩统计分布Figure 6. Magnetic moments variation and statistics under thermal effect. (a) Variation of magnetic moment z component with time, (b) statistical distribution of magnetic moment对于同一个磁矩，时间是没有影响的，可以认为不同时间的磁矩指向是完全独立的。可以看出磁矩分布近似成指数分布，这里需要注意的是，磁矩为1的时候能量最小。为了更进一步研究磁矩统计分布情况，对式（21）两边同时取对数，得到$ \ln P = - \beta E + C $。通过取对数，可以更清晰地对比出理论与模拟之间的关系。根据模拟结果统计出 $ \ln P $与理论值$ - \beta E + C $相对比，得到图7，图7中蓝色圆点为模拟结果统计，红色直线为理论值。
由图7可以看出在能量较低时统计结果和理论符合得非常好，但随着能量增大理论值与统计值出现了较小的偏差，这主要是由两个原因造成的：一是磁矩被激发到高能级的概率较低，导致高能量时数据稀少。二是能量越高越偏离平衡态，会导致系统的能量分布偏离玻尔兹曼分布。
2.4 飞秒脉冲激励下的太赫兹辐射
前面我们用微磁学模拟方法模拟了一维自旋波的色散关系、磁畴壁和热平衡状态下的磁矩分布。自旋波和磁畴壁是飞秒激光作用于Fe薄膜内部后，在其内可能存在的微磁学现象，热效应的等效场更是连接了三温度模型和LLG方程，这为模拟飞秒激光与Fe薄膜的作用提供了支持。为了简化模型，只考虑在非磁性金属层上生长的磁性Fe薄膜。主要考虑飞秒激光作用在Fe薄膜上产生磁矩扰动的过程，从而产生纯自旋流。在模拟过程中，忽略了飞秒激光对可能产生自旋极化电流并辐射电磁波的铁磁材料的电子分布密度的影响。由于Fe薄膜（10 nm）的线性度远小于太赫兹波长（约10−4 m），因此整个Fe薄膜可以近似为点源，相当于Fe薄膜的平均磁矩辐射。所以可以使用随机LLG方程模拟Fe薄膜中的磁动力学过程，然后通过经典的磁偶极子辐射理论获得了飞秒激光诱导的太赫兹辐射。
图8为模拟得到的太赫兹辐射图谱。图8（a）为Fe薄膜辐射的太赫兹波的电场强度Ey，这与实验结果一致[3-4,15]。由图8（a）可知，电场强度Ey在0.3 ps 时达到了最大值。通过对太赫兹辐射的电场强度的时域关系图8（a）进行傅里叶变换，得到了辐射电场的频谱图8（b）。可以看出辐射的电场强度的$ y $方向的频谱连续性较好，这意味着辐射的太赫兹波具有一个良好信号的特征，没有突变或间断。它能够提供更多的信息、更高的数据传输速率和更好的调制灵活性，可以用作高传输效率的通讯媒介。
![]()
图 8 飞秒激光作用于Fe 薄膜上后辐射的电场强度的y分量及其频谱。（a）电场强度的y分量，（b）电场强度的y分量的频谱Figure 8. The y component and spectrum of electric field intensity radiated by femtosecond laser acting on Fe thin film. (a) The y component of the electric field intensity, (b) spectrum of the y component of the electric field intensity3. 结论
通过基于 LLG 方程的微磁学模拟研究了微磁学领域的一维自旋波的色散关系、磁畴壁和热效应下的磁矩统计分布，并将模拟结果与理论进行对比验证，证实了微磁学模拟的可行性和准确性。这些微磁学常见现象也是飞秒激光作用Fe薄膜后其内部可能存在的微磁学现象，热效应的等效场更是连接了三温度模型和LLG方程，这些模拟过程和结果为模拟飞秒激光与 Fe 薄膜的作用打下了基础。由模拟的飞秒激光与Fe薄膜的作用结果可以看出辐射的电场强度的$ y $方向的频谱连续性较好，具有一个良好信号的特征，能够提供更多的信息、更高的数据传输速率和更好的调制灵活性，可以用作高传输效率的通讯媒介。
-
![]()
图 2 磁化强度 y 分量模拟结果。（a）单位磁化强度 y 分量变化情况，（b）自旋波色散关系
Figure 2. Simulation results of y component of magnetization. (a) Variation of y component of unit magnetization, (b) spin wave dispersion relationship
![]()
图 3 不同$ \phi $对应的磁畴壁结构。（a）奈尔壁结构图， （b）布洛赫壁结构图
Figure 3. Magnetic domain walls corresponding to different $ \phi $. (a) Neel wall structure, (b) Bloch wall structure
![]()
图 4 二维磁畴壁模拟结构示意图
Figure 4. Schematic diagram of two-dimensional magnetic domain wall simulation structure
![]()
图 6 热效应下的磁矩变化和统计。（a）磁矩 z 分量随时间的变化，（b）磁矩统计分布
Figure 6. Magnetic moments variation and statistics under thermal effect. (a) Variation of magnetic moment z component with time, (b) statistical distribution of magnetic moment
![]()
图 8 飞秒激光作用于Fe 薄膜上后辐射的电场强度的y分量及其频谱。（a）电场强度的y分量，（b）电场强度的y分量的频谱
Figure 8. The y component and spectrum of electric field intensity radiated by femtosecond laser acting on Fe thin film. (a) The y component of the electric field intensity, (b) spectrum of the y component of the electric field intensity
-
[1] Pettine J, Padmanabhan P, Sirica N, et al. Ultrafast terahertz emission from emerging symmetry-broken materials[J]. Light-Science & Applications,2023,12(1):133
[2] Leitenstorfer A, Moskalenko A S, Kampfrath T, et al. The 2023 terahertz science and technology roadmap[J]. Journal of Physics D-Applied Physics,2023,56(22):223001 doi: 10.1088/1361-6463/acbe4c
[3] 张顺浓, 朱伟骅, 李炬赓, 等. 铁磁异质结构中的超快自旋流调制实现相干太赫兹辐射[J]. 物理学报,2018,67(19):197202 (in Chinese) doi: 10.7498/aps.67.20181178 Zhang S N, Zhu W H, Li J G, et al. Coherent terahertz radiation via ultrafast manipulation of spin currents in ferromagnetic heterostructures[J]. Acta Physica Sinica,2018,67(19):197202 doi: 10.7498/aps.67.20181178
[4] Huang L, Lee S H, Kim S D, et al. Universal field-tunable terahertz emission by ultrafast photoinduced demagnetization in Fe, Ni, and Co ferromagnetic films[J]. Scientific Reports,2020,10(1):15843 doi: 10.1038/s41598-020-72855-1
[5] Wu Y, Elyasi M, Qiu X P, et al. High-performance Thz emitters based on ferromagnetic/nonmagnetic heterostructures[J]. Advanced Materials,2017,29(4):1603031 doi: 10.1002/adma.201603031
[6] Li G, Medapalli R, Mikhaylovskiy R V, et al. Thz emission from Co/Pt bilayers with varied roughness, crystal structure, and interface intermixing[J]. Physical Review Materials,2019,3(8):084415 doi: 10.1103/PhysRevMaterials.3.084415
[7] Salikhov R, Ilyakov I, Körber L, et al. Coupling of terahertz light with nanometre-wavelength magnon modes via spin–orbit torque[J]. Nature Physics,2023,19(4):529−535 doi: 10.1038/s41567-022-01908-1
[8] Rongione E, Gueckstock O, Mattern M, et al. Emission of coherent Thz magnons in an antiferromagnetic insulator triggered by ultrafast spin-phonon interactions[J]. Nature Communications,2023,14(1):1818 doi: 10.1038/s41467-023-37509-6
[9] Yang Y S, Dal Forno S, Battiato M. Modeling spintronic terahertz emitters as a function of spin generation and diffusion geometry[J]. Physical Review B,2023,107(14):144407 doi: 10.1103/PhysRevB.107.144407
[10] Zhang W T, Maldonado P, Jin Z M, et al. Ultrafast terahertz magnetometry[J]. Nature Communications,2020,11(1):4247 doi: 10.1038/s41467-020-17935-6
[11] Liu L W, Ren Y F, Robert W, et al. Terahertz wave radiation simulation in the Fe thin film[J]. Journal of Physics: Condensed. Matter,2024,36:305902 doi: 10.1088/1361-648X/ad4222
[12] 褚欣博, 金钻明, 吴旭, 等. 铁磁异质结的远红外脉冲辐射及其光热调控研究[J]. 物理学报,2023,72:157801 (in Chinese) doi: 10.7498/aps.72.20230543 Chu X B, Jin Z M, Wu X, et al. Pulsed far-infrared generation in ferromagnetic heterostructue controlled by photo-thermal effect[J]. Acta Physica Sinica,2023,72:157801 doi: 10.7498/aps.72.20230543
[13] Beaurepaire E, Merle J C, Daunois A, et al. Ultrafast spin dynamics in ferromagnetic nickel[J]. Physical review letters,1996,76(22):4250−4253 doi: 10.1103/PhysRevLett.76.4250
[14] Zhang J, Yu W, Chen X, et al. A frequency-domain micromagnetic simulation module based on comsol multiphysics[J]. AIP Advances,2023,13(5):055108 doi: 10.1063/5.0143262
[15] Zhang S, Jin Z, Zhu Z, et al. Bursts of efficient terahertz radiation with saturation effect from metal-based ferromagnetic heterostructures[J]. Journal of Physics D: Applied Physics,2018,51(3):034001 doi: 10.1088/1361-6463/aa9e43
下载:
计量
- 文章访问数: 34
- HTML全文浏览量: 0
- PDF下载量: 6
