Phase field crystal modeling of grain boundary migration: Mobility, energy and structural variability
相场晶体模型中的晶界迁移：迁移率、能量和结构可变性

Characterization of grain boundaries (GBs) and identification of their properties is central to advance the understanding and engineering of crystalline materials. Yet it remains a complex and largely unresolved challenge to fully characterize the effect GBs exert on the evolution of polycrystal microstructures. Most GB properties, such as GB energy, have traditionally been derived under assumptions of equilibrium, often only focusing on stationary GBs in their lowest-energy configurations. These simplified models do not, however, fully capture the complexity of real-world systems, where a vast range of alternative GB configurations may be found. This is due, in part, to the variations made possible by the five macroscopic degrees of freedom (DOF) defining a GB configuration, along with three microscopic DOF, representing relative translations of the adjoining grains within the displacement shift complete (DSC) lattice. Variations in this full set of DOF affect key GB properties such as GB energy [1]. Additionally, other microscopic DOF can be defined, such as the excess volume of the GB and the chemical composition of the GB. Studies that do account for GB energy variations, such as [2], tend to show better agreement with experiments, underscoring the importance of these considerations.
晶界（GBs）的表征及其特性的识别是推进对晶体材料理解和工程化的关键。然而，全面表征晶界对多晶材料微观结构演变的影响仍然是一个复杂且尚未完全解决的挑战。大多数晶界特性，如晶界能，传统上是在平衡假设下推导得出的，通常只关注处于最低能态的静止晶界。然而，这些简化模型并不能完全捕捉现实世界系统的复杂性，因为在现实世界中可能存在广泛的其他晶界构型。这部分是由于定义晶界构型的五个宏观自由度（DOF）以及三个微观自由度（代表相邻晶粒在位移移位完全（DSC）晶格中的相对平移）所允许的变化。这一整套自由度的变化会影响关键的晶界特性，如晶界能[1]。此外，还可以定义其他微观自由度，例如晶界的过量体积和晶界的化学成分。 对晶界能量变化进行考虑的研究，例如[2]，往往与实验结果更为吻合，这突显了这些考虑的重要性。

Experimental constraints further complicate the study of GB properties and dynamics, as conducting in-situ measurements at the atomic scale over meaningful timescales ($>1$ second [3], [4]) and sufficient spatial resolution remains challenging, especially in 3D. Although recent advancements, such as those presented in [4], have made it possible to observe GB dynamics over extended time scales and spatial range, the studies are usually limited to particular GB configurations under specific conditions.
实验限制进一步增加了对晶界性质和动力学的研究难度，因为在有意义的时标（ $>1$ 秒[3], [4]）内进行原子尺度的原位测量，并保持足够的空间分辨率仍然具有挑战性，尤其是在三维情况下。尽管像[4]中展示的近期进展使得在较长时间尺度和空间范围内观察晶界动力学成为可能，但这类研究通常仅限于在特定条件下特定的晶界构型。

Due to the experimental challenges, computational simulations, particularly those based on molecular dynamics (MD), have become essential for investigating GB properties and behaviors. Simulations are especially valuable when stationary GBs fail to reveal certain dynamic properties, such as GB migration velocity. But numerical simulations also have their limitations: MD studies are typically confined to short timescales and small domain sizes, often only capturing migration over a few picoseconds or nanoseconds [5], [6], which may not accurately represent the long-term dynamic processes occurring in real materials during GB migration. MD-based GB studies involving longer time spans are scarce, but examples can be found in [7], [8].
由于实验上的挑战，计算模拟，特别是基于分子动力学（MD）的模拟，已成为研究晶界（GB）性质和行为的重要手段。当静态晶界无法揭示某些动态特性时，如晶界迁移速度，模拟尤其有价值。但数值模拟也有其局限性：MD 研究通常局限于较短的时间尺度和较小的域尺寸，往往只能捕捉到几皮秒或纳秒内的迁移过程[5], [6]，这可能无法准确反映晶界迁移过程中真实材料中发生的长期动态过程。涉及更长时间跨度的基于 MD 的晶界研究较为罕见，但可以在[7], [8]中找到例子。

One particularly difficult dynamical property to characterize is the GB mobility, a key parameter in the evolution of polycrystalline systems. The mobility can be identified by noting that the GB migration velocity is commonly expressed as
表征晶界迁移率这一特别困难的动力学特性，而迁移率是多晶系统演化中的一个关键参数。可以通过注意到晶界迁移速度通常表示为...来识别迁移率。(1)$v=MP$where $M$ is the mobility parameter, usually assumed to follow an Arrhenius-like temperature dependence, and $P$ is the driving pressure acting on the GB. Some models propose alternative formulations, such as $v=M{P}^n$, where $n$ is an exponent varying between 1 and 3 [9].
其中 $M$ 是迁移率参数，通常假定遵循类似阿伦尼乌斯型的温度依赖关系，而 $P$ 是作用在晶界上的驱动力。一些模型提出了替代公式，例如 $v=M{P}^n$ ，其中 $n$ 是在 1 和 3 之间变化的指数 [9]。

Although the temperature dependence of $M$ is well established to follow an Arrhenius-like behavior, as demonstrated in [10], [11], significant jumps in mobility can still occur due to thermally induced GB roughening at high temperatures, as discussed in [5], [12]. Temperature effects are not, however, the primary focus of the present study which is conducted at constant, non-dimensional, temperature. The influence of other DOF on mobility, such as GB orientation, remains less well understood. Orientation effects on GB mobility have been observed experimentally in [10] and by simulations in [5], [13], [14]. Moreover, GB mobility has been shown to also depend on the magnitude of the driving pressure [15], [16], [17].
尽管 $M$ 的温度依赖性已被证实遵循类似阿伦尼乌斯型行为，如 [10] 和 [11] 所示，但由于高温下的热诱导晶界粗糙化，迁移率仍可能发生显著跳跃，如 [5] 和 [12] 所讨论。然而，温度效应并非本研究的重点，该研究是在恒定的无量纲温度下进行的。其他自由度对迁移率的影响，例如晶界取向，仍了解较少。晶界取向对迁移率的影响已在实验 [10] 和模拟 [5]、[13]、[14] 中观察到。此外，晶界迁移率已被证明还取决于驱动力的大小 [15]、[16]、[17]。

Additional complexity arise as the driving pressure can stem from multiple origins, including, for example, stored deformation energy, chemical driving forces or magnetic fields [3]. In many cases, however, capillarity-driven migration is assumed, governed by the GB mean curvature. For this commonly adopted curvature-driven growth model, the driving pressure appears as
随着驱动力可能源自多个来源，包括储存的变形能、化学驱动力或磁场[3]，复杂性会增加。然而，在许多情况下，通常假设毛细作用驱动的迁移，由界面平均曲率控制。对于这种通常采用的曲率驱动生长模型，驱动力表现为(2)$P=\mathit{\kappa }:\mathit{\Gamma }$where $\mathit{\kappa }$ is the GB curvature tensor and $\mathit{\Gamma }$ the GB stiffness tensor. Recent experimental and computational findings suggest, however, that this model may overlook critical phenomena, such as GB migration stagnation via disruptive atomic jumps [7], energy dissipation via GB replacement [18], chemical composition at the GB [19] or the influence of topological constraints like GB junctions [20]. Other factors, such as shear-coupled GB migration, may also influence GB mobility [14], [21], [22], [23]. Shear coupling relates the GB migration velocity, $v_{\perp }$, to the shear velocity, $v_{\parallel }$, through the coupling factor
其中 $\mathit{\kappa }$ 是界面曲率张量， $\mathit{\Gamma }$ 是界面刚度张量。然而，最近的实验和计算发现表明，该模型可能忽略了关键现象，例如界面迁移停滞通过破坏性原子跳跃[7]、通过界面替代的能量耗散[18]、界面处的化学成分[19]或拓扑约束（如界面结点）的影响[20]。其他因素，如剪切耦合的界面迁移，也可能影响界面迁移率[14]、[21]、[22]、[23]。剪切耦合将界面迁移速度 $v_{\perp }$ 与剪切速度 $v_{\parallel }$ 通过耦合因子联系起来(3)$\beta =\frac{v_{\parallel }}{v_{\perp }}$

In addition to its dependence on GB orientation, shear coupling has also been shown to vary with both the origin of the GB driving pressure — like a jump in chemical potential or originating from a difference in stress state — and its magnitude [23].
除取决于晶界取向外，剪切耦合还显示出与晶界驱动压力的来源（如化学势的跃迁或源自应力状态差异）及其幅值的变化有关[23]。

Most GB-related research relies on assuming stationary, minimum-energy GB configurations, which may fail to capture essential mechanics of dynamical GB behavior. Notably, [5] expanded the study of GB mobility by using MD to assess various GBs. The study is limited, however, by the inherently short timescales of MD. Other simulation-based GB investigations have examined highly curved boundaries, such as in [24], which may not be representative of real-world, less curved, flat, or stepped GBs. GB studies also commonly focus on symmetrical tilt GBs, as is also the case in this study, which might not be representative of all possible GB migration behaviors and as mixed GBs, for example, might have distinctly different properties, as illustrated in [25].
大多数与晶界相关的研究依赖于假设静止的、最小能量的晶界构型，这可能无法捕捉到动态晶界行为的基本力学。值得注意的是，[5]通过使用分子动力学来评估各种晶界，扩展了对晶界迁移率的研究。然而，该研究受限于分子动力学的固有短时间尺度。其他基于模拟的晶界研究考察了高度弯曲的边界，例如在[24]中，这些可能无法代表现实世界中较不弯曲、平坦或阶梯状的晶界。晶界研究通常也关注对称倾转晶界，本研究的情形也是如此，这可能无法代表所有可能的晶界迁移行为，因为混合晶界（例如）可能具有显著不同的特性，如[25]所示。

This study aims to address some of the common assumptions and limitations in simulation-based GB characterization by investigating the migration dynamics of flat symmetrical tilt GBs using the phase field crystal (PFC) method. The PFC method offers atomic spatial resolution over diffusional timescales, making it well suited for capturing long-term GB evolution, beyond the capabilities of conventional MD simulations.
本研究旨在通过使用相场晶体（PFC）方法研究平面对称倾转晶界的迁移动力学，来探讨基于模拟的晶界表征中的一些常见假设和局限性。PFC 方法在扩散时间尺度上提供原子级空间分辨率，使其非常适合捕捉长期晶界演化，这超出了传统分子动力学模拟的能力。

The structure of this paper is organized as follows: Section 2 introduces the PFC model employed in this study and Section 3 describes the simulation model setup and data extraction procedures. The results are presented and discussed in Section 4, followed by conclusions and some final remarks in Section 5.
本文的结构安排如下：第二节介绍了本研究中使用的 PFC 模型，第三节描述了模拟模型设置和数据提取程序。第四节展示并讨论了结果，第五节为结论和最终评论。

The PFC formulation is based on energy minimization, where the total normalized free energy, $F$, is given by the sum of the ideal gas energy, $F_{id}$, and the excess energy, $F_{exc}$, following
PFC 公式基于能量最小化，其中总归一化自由能 $F$ 由理想气体能量 $F_{id}$ 和过剩能量 $F_{exc}$ 之和给出，遵循(4)$\frac{F\left[n\left(\mathit{r}\right)\right]}{k_{b}T{\rho }_0}=F_{id}\left[n\left(\mathit{r}\right)\right]+F_{exc}\left[n\left(\mathit{r}\right)\right]$which is a functional of the density field $n$, itself a function of the spatial coordinates $\mathit{r}(x,y,z)$. Furthermore, $k_b$ is the Boltzmann constant, $T$ the absolute temperature and ${\rho }_0$ the reference density. Following [26], the ideal gas energy is expressed as
这是一个关于密度场 $n$ 的泛函，而密度场本身是空间坐标 $\mathit{r}(x,y,z)$ 的函数。此外， $k_b$ 是玻尔兹曼常数， $T$ 是绝对温度， ${\rho }_0$ 是参考密度。根据[26]，理想气体能量表示为(5)$F_{id}\left[n\left(\mathit{r}\right)\right]=\int d\mathit{r}\left(\frac{1}{2}n{\left(\mathit{r}\right)}^2-\frac{1}{6}n{\left(\mathit{r}\right)}^3+\frac{1}{12}n{\left(\mathit{r}\right)}^4\right)$whereas the excess energy is taken as
而多余的能量被当作(6)$F_{exc}\left[n\left(\mathit{r}\right)\right]=-\frac{1}{2}\iint d\mathit{r}d{\mathit{r}}^{\prime }n\left(\mathit{r}\right)n\left({\mathit{r}}^{\prime }\right)C_2\left(|{\mathit{r}}^{\prime }-\mathit{r}|\right)=-\frac{1}{2}\int d\mathit{r}n\left(\mathit{r}\right)(C_2\ast n)\left(\mathit{r}\right)$where $C_2$ is the two-point correlation function and $\ast$ denotes a convolution. As FCC and BCC lattices are considered in this work, the restriction to a two-point correlation function is appropriate. Higher-order correlations could be used, however, for more complex crystal structures, as done for diamond cubic lattices in [27]. Following the structural PFC (XPFC) formalism introduced in [26], the Fourier transform of $C_2$, denoted as ${\widetilde{C_2}}_{,k}$, is given by the loci of a sum of Gaussian peaks:
其中 $C_2$ 是两点相关函数， $\ast$ 表示卷积。由于本工作考虑了面心立方和体心立方晶格，限制为两点相关函数是恰当的。然而，对于更复杂的晶体结构，可以使用高阶相关，正如在[27]中对金刚石立方晶格所做的那样。遵循在[26]中引入的结构性相场晶体（XPFC）形式主义， $C_2$ 的傅里叶变换，记作 ${\widetilde{C_2}}_{,k}$ ，由一系列高斯峰值的和的轨迹给出：(7)${\widetilde{C_2}}_{,k}\left(\mathit{k}\right)=\underset{i}{max}\left[c_{i}exp\left(-\frac{{(\left|\mathit{k}\right|-k_i)}^2}{2{\alpha }_i^2}\right)\right]$where $\mathit{k}$ is the reciprocal vector, $c_i$ is the amplitude of the $i$’th symmetry plane, $k_i$ is its wave number and ${\alpha }_i$ is the width of the Gaussian peak. To stabilize an FCC phase, two peaks located at $k_1=2\pi \sqrt{3}$ and $k_2=4\pi$ are used. Peaks at $k_1=2\pi \sqrt{2}$ and $k_2=4\pi$ are used to stabilize a BCC phase. The amplitudes $c_i$ are set to unity and the widths to ${\alpha }_i=2$ for all peaks. The widths were set according to [1], [28] to facilitate direct comparison with those previous studies. It should be noted, however, that the widths of the peaks affect properties such as the width of liquid–solid interfaces [26] and elastic properties [29], [30]. Simulations also show that decreasing the Gaussian width has some effect on the results. For instance, lowering the value of ${\alpha }_i$ to unity, a common standard choice for XPFC simulations, results in a smaller variation of the GB energy. As the GB energy presented in [1], [28] appears to be in good agreement with MD simulations, setting ${\alpha }_i=2$ is considered a reasonable choice, and though values of GB velocity and energy may change quantitatively, their qualitative variations remain similar.
其中 $\mathit{k}$ 是倒易矢量， $c_i$ 是第 $i$ 个对称平面的振幅， $k_i$ 是其波数， ${\alpha }_i$ 是高斯峰的宽度。为了稳定面心立方相，使用了位于 $k_1=2\pi \sqrt{3}$ 和 $k_2=4\pi$ 的两个峰。为了稳定体心立方相，使用了位于 $k_1=2\pi \sqrt{2}$ 和 $k_2=4\pi$ 的两个峰。振幅 $c_i$ 被设置为 1，宽度被设置为 ${\alpha }_i=2$ 。宽度是根据 [1]、[28] 设置的，以便与之前的研究直接比较。然而，需要注意的是，峰的宽度会影响诸如液固界面宽度 [26] 和弹性性质 [29]、[30] 等性质。模拟也表明，减小高斯宽度对结果有一定影响。例如，将 ${\alpha }_i$ 的值降低到 1（这是 XPFC 模拟的常用标准选择），会导致 GB 能量变化较小。由于 [1]、[28] 中给出的 GB 能量与分子动力学模拟结果吻合较好，因此认为设置 ${\alpha }_i=2$ 是一个合理的选择，尽管 GB 速度和能量的数值可能发生变化，但它们的定性变化仍然相似。

To allow for a differential pressure across flat GBs, a directional correlation function is added to the model. This has the effect of energetically favoring a specific crystal orientation, which could be likened with the introduction of a stored energy jump across the GB and is similar to how a GB driving pressure is achieved in some MD studies, for example in [25], [31]. Following [32], the directional correlation function is defined as
为了在平面晶界间产生压差，模型中添加了一个方向相关性函数。这会在能量上优先选择特定的晶体取向，类似于在晶界间引入了储存能的跃迁，并且类似于某些分子动力学研究中实现晶界驱动压力的方式，例如在[25]、[31]中。根据[32]，方向相关性函数定义为(8)${\widetilde{H_2}}_{,k}\left(\mathit{k}\right)=H_0\underset{i}{max}\left[\sum _{i}exp\left(-\frac{{\Vert \mathit{k}-{\mathit{k}}_i\Vert }^2}{2{\beta }_i^2}\right)\right]$where $H_0$ is the amplitude, ${\mathit{k}}_i$ the center position and ${\beta }_i$ the width of a Gaussian peak. The amplitude $H_0$ is set equal for all peaks, with a value chosen to get a correct energy difference between two grains in a bicrystal, as discussed further in Section 3.2. The reciprocal vectors ${\mathit{k}}_i$ are set to coincide with the $〈111〉$ reciprocal lattice vectors for FCC lattices and the $〈110〉$ reciprocal lattice vectors for BCC lattices, with orientations equal to the target crystal orientation. Finally, the peak widths are set to ${\beta }_i=0.6$ for all vectors ${\mathit{k}}_i$. This choice of $\beta$ defines relatively sharp peaks, which is set to minimize interactions with neighboring crystals that only differ by small values of disorientation—the smallest possible misorientation between two crystals, accounting for switching symmetries.
其中 $H_0$ 表示高斯峰的振幅， ${\mathit{k}}_i$ 表示中心位置， ${\beta }_i$ 表示高斯峰的宽度。振幅 $H_0$ 被设置为所有峰相等，其值的选择是为了在双晶体的两个晶粒之间得到正确的能量差，这将在第 3.2 节中进一步讨论。倒易矢量 ${\mathit{k}}_i$ 被设置为与面心立方（FCC）晶格的 $〈111〉$ 倒易矢量重合，以及与体心立方（BCC）晶格的 $〈110〉$ 倒易矢量重合，其方向等于目标晶体取向。最后，所有矢量 ${\mathit{k}}_i$ 的峰宽被设置为 ${\beta }_i=0.6$ 。这种选择 $\beta$ 定义了相对尖锐的峰，这是为了最小化与仅具有微小取向差值的相邻晶体的相互作用——两个晶体之间可能的最小取向差值，考虑了转换对称性。

The directional correlation function, Eq. (8), can be set to both increase the energy of one crystal orientation by having a negative amplitude, or to reduce the energy of the same crystal by having a positive amplitude. This gives the option to define the pressure as either pushing or pulling the GB. For this work we found the positive energy contribution to give the most stable simulation conditions for higher pressures. This choice have some effect on certain simulated cases, as also noted for MD in [5], where different modes of migration are observed. These differences are, however, seen as a second-order phenomenon and characterization of the differences is considered out of scope for the current study.
方向相关性函数，公式(8)，可以设置为通过具有负振幅来增加某一晶体取向的能量，或通过具有正振幅来减少同一晶体的能量。这提供了定义压力是推动还是拉动晶界的选项。对于这项工作，我们发现正能量贡献为高压提供了最稳定的模拟条件。这种选择对某些模拟案例有一定影响，正如在 MD[5]中提到的，其中观察到不同的迁移模式。然而，这些差异被视为二级现象，对差异的表征被认为超出了当前研究的范围。

With the addition of the directional correlation function to Eq. (4), the free energy appears as
在将方向相关性函数添加到公式(4)后，自由能表现为(9)$\frac{F\left[n\right]}{k_{b}T{\rho }_0}=\int d\mathit{r}\left(\frac{1}{2}{n}^2-\frac{1}{6}{n}^3+\frac{1}{12}{n}^4-\frac{1}{2}n\left[n\ast (C_2+H_2)\right]\right)$where the explicit dependence on the coordinates $\mathit{r}$ has been omitted for brevity and where $H_2$ is the reverse Fourier transform of ${\widetilde{H_2}}_{,k}$.
其中为了简洁起见省略了对坐标 $\mathit{r}$ 的显式依赖，而 $H_2$ 是 ${\widetilde{H_2}}_{,k}$ 的逆傅里叶变换。

To study flat GBs, a simple bicrystal geometry is employed, as illustrated in Fig. 1. The bicrystal comprises two crystals, $G_1$ and $G_2$, which have different orientations relative to each other. Due to the periodic boundary conditions, which are required by the spectral method, two identical GBs are modeled. As an alternative to having two GBs, a single GB domain could have been used, as done in the bicrystal PFC model in [1], and as also commonly employed in MD simulations [5], [23]. While it is relatively straightforward to include free surfaces in MD simulations, incorporating them in PFC using the methodology in [1] remains to be tested for the study of migrating GBs. Nevertheless, a single GB domain may yield different results with respect to shear-coupled GB migration, which may alter the GB migration rate in some cases. This issue is less prominent in the two-GB geometry, as the shear-coupled migration of both boundaries effectively cancels out, restricting the translation of the bulk crystals. To better capture a constrained GB, the two-GB setup is adopted.
为研究平面晶界，采用了一种简单的双晶几何结构，如图 1 所示。该双晶由两个取向不同的晶体 $G_1$ 和 $G_2$ 组成。由于光谱方法需要周期性边界条件，因此模拟了两个相同的晶界。作为使用两个晶界的替代方案，也可以使用单个晶界区域，正如文献[1]中双晶相场晶格模型所做的那样，并且这也是分子动力学模拟中常用的方法[5],[23]。虽然将自由表面包含在分子动力学模拟中相对直接，但使用文献[1]中的方法在相场晶格模型中包含它们，对于研究迁移的晶界仍需进行测试。然而，单个晶界区域可能会在剪切耦合晶界迁移方面产生不同的结果，这可能会在某些情况下改变晶界迁移速率。在双晶界几何结构中，这个问题不太明显，因为两个边界的剪切耦合迁移会相互抵消，从而限制体相晶体的平移。为了更好地捕捉受约束的晶界，采用了双晶界设置。