Σ3 grain boundary dynamics studied by atomistic spherical bicrystal modeling

In this study, we propose an atomistic spherical bicrystal model that employs the synthetic driving force method to investigate the dynamics of curved interfaces. The model allows tracking the migration and faceting of spherical grain boundaries. As a demonstration case, a simulation for a grain boundary with Σ3 misorientation is performed. A subsequent analysis of interface energy and stiffness is performed, whereby it is found that the dynamics of [111] and [1¯1¯2] segments in the Σ3 boundary plane fundamental zone follow energy and stiffness trends, while others, e.g., the [01¯1] segment, do not.


Introduction
Molecular statics/dynamics simulations are commonly used numerical methods for investigating various grain boundary (GB) properties, such as energy and mobility, as well as boundary migration and faceting.Conventional planar bicrystal models, i.e., periodic bicrystal systems containing flat GBs with specific misorientations, are frequently employed as a standard simulation method for such studies [1][2][3].These models exclude effects from free surfaces.However, most GBs in real materials are curved, and energy and mobility can vary significantly along the GBs, influencing the overall dynamics of the boundary.Additionally, anisotropy in the angle of the GB plane normal also plays a crucial role in the dynamic processes of curved boundaries, which is typically described by the GB stiffness parameter.The stiffness of a GB can be calculated by summing the interface energy and the second derivative of that energy with respect to the boundary plane (BP) orientation, where the stiffness tensor Г(n) is defined by Г(n)=γI+∇n∇nγ.Here γ represents the free energy, I is the identity tensor, and ∇n is the surface gradient with respect to n, the GB normal.It is worth noting that the second term depends on the tangential directions.Previous research [4] has demonstrated that the GB stiffness strongly depends on the inclination change of a boundary.Thus, the evolution of curved GBs in real materials is complex and involves integrated effects from energy, mobility, and stiffness.Studying these combined effects in conventional planar bicrystal models has certain limitations.
Here we employ a spherical bicrystal (SB) model for studying curved boundaries for a better understanding of grain boundary dynamics [5,6].The strategy is to create a spherical misoriented grain that is embedded in a periodic cubic cell, which allows for studies of spherical interfaces under specific conditions.One key advantage of the SB model is that a full subspace of BP inclinations associated with one specific crystallographic misorientation is accessible in a single simulation run, as opposed to conventional studies that explore only one GB type (i.e., one selected misorientation and one BP inclination).Furthermore, the evolution of BP inclinations on the sphere may be greatly affected by neighboring BP inclinations, enabling the study of interactions between different BP inclinations.This provides rich information about interface dynamics, making the SB model especially suitable for investigating effects of GB stiffness.For instance, previous work by Abdeljawad et al. [6,7] analyzed the impact of GB stiffness on equilibrium shapes and migration for Σ3-11 boundaries in nickel using SB models.In these works, the stiffness is demonstrated to play a critical role in related dynamical processes, as the stiffness contribution is larger and more anisotropic than that of the energy.
In this work, we set up a new simulation procedure where we combine the SB model with the synthetic driving force method, which is typically used in planar bicrystal models, to stimulate the inward migration of a curved GB.The motivation for this is to stimulate migration of all types of boundaries.In the classic SB models, boundaries such as Σ3 and Σ7, do migrate driven by the initial curvature.The curvature however works less efficiently for boundaries like Σ5.Our new approach of accelerating the process enables dynamics SB studies for all misorientations and is thus particularly useful for comparisons between boundary dynamics.The specific aim of this paper is however to study the Σ3 GB dynamics only.This is because the availability of extensive datasets about Σ3 boundaries, and the simulated behaviors, such as faceting, may provide valuable data for a better understanding of the mechanism governing the dynamical processes of curved GBs.The work will be the starting point for studying many more GB misorientations.Understanding the dynamics of GBs is crucial for advancing GB engineering to design superior materials.

Simulation details and analysis method
A spherical Σ3 GBs is investigated in the present work.The atomic interactions of the aluminum atoms are modeled in LAMMPS [8], using the Mishin-Farkas embedded-atom method (EAM) interatomic potential [9].The potential has been validated previously and found to be in excellent agreement with physical properties, including the generalized stacking fault energy curve [10].The visualizations were generated using Ovito [11].
Figure 1 shows the 400 Å edge length cubic simulation cell filled with face-centered cubic lattice atoms.The cell is periodic in all directions, containing ~3.86 million atoms.The recrystallizing grain (blue) is cube-oriented, consistent with the experimental coordinate system (x-y-z in figure 1), i.e., the [100], [010], [001] axes are along the experimental x, y, and z axes.To initiate the simulation, the system is equilibrated at 600 K, with zero pressure in all directions.The energyconserving orientational (ECO) driving force in its recently developed efficient implementation [12] is utilized to create a synthetic energy difference of 0.01 eV between the recrystallizing grain and the matrix.This drives the spherical interface to migrate towards the center grain.
In this study, we utilize a strategy proposed by Abdeljawad et al. [6] to determine the stiffness of the Σ3 GBs.Specifically, we employ a spherical harmonic fit with {2nθ, 6mφ} modes where θ and φ are polar and azimuthal angles respectively to represent the energy distribution BP inclination spaces.Considering that Σ3 boundaries on the sphere are symmetrical around the [111] axis, we rotate the coordinate system so that the [111] axis is along the z-axis to facilitate the harmonic fit. Figure 2 illustrates the Σ3 BP fundamental zone (BPFZ) [13] and the spherical coordinate system utilized for GB energy fitting.Per the definition of GB stiffness, the calculated stiffness value should be much more sensitive to the employed fitting equation due to the contribution from the second derivative term.So in the present case of aluminum, two harmonic equations based on {2θ, 4θ, 6φ, 12φ} modes and their cross terms are used in the calculation.The GB energy values are subsequently defined using the following equations (based on these equations the stiffness can be calculated) γgb = α0+α1+cos(2θ)+α2cos(4θ)+α3cos(6φ)+α4cos(2θ)cos(6φ)+α5cos(4θ)cos(6φ) (1) γgb = α0+α1+cos(2θ)+α2cos(4θ)+α3cos(6φ)+α4cos(2θ)cos(6φ)+α5cos(4θ)cos(6φ)+α6cos(12φ)

Spherical Σ3 interface migration
Figure 3a-c depict the configurations of the spherical Σ3 interface at initial, 10 ps, and 50 ps at 600 K, respectively.The GB atoms are identified according to the ECO parameters [14].As most of the initial spherical interface inclinations deviate from the Gibbs-Wulff equilibrium shape, the embedded grain undergoes anisotropic shrinking during annealing, exhibiting a strong dependence on the boundary plane inclination.The spherical interface gradually evolves towards a stable configuration with several faceted segments, and the top and bottom poles ([111] segments) hardly move, as evidenced by the atomistic plots.
To provide a clearer picture of the interface shape changes, we present cross-section plots in figures 3d-f.The [111] view cross-section (figure 3d) reveals that the initial circular trace on this cross-section completely transforms into a hexagram shape during the first 6 ps.The [01 ̅ 1] segments exhibit the highest mobility, while [1 ̅ 1 ̅ 2] segments remain almost stationary.Faceting occurs until most visible interfacial segments on this cross-section become perpendicular to [1 ̅ 1 ̅ 2] (or its equivalents).These newly formed facets are stable and move little thereafter.Other stationary segments mostly belong to

Discussion of factors impacting the evolution
Here we discuss the various factors that could impact the evolution of the initially spherical Σ3 GB.We consider the possible effects of energy, mobility, and stiffness.
Figure 4a and b show the fitted Σ3 GB energy on the spherical interface based on equations 1 and 2, respectively.The scattered atomistic energy data derived from Olmsted et al. [1] are plotted in the BPFZ.A comparison between the two sets of data is presented in figure 4c and d, demonstrating good fits with high fitting coefficients (R-squared) of 0.972 and 0.978, respectively.The slightly improved fit achieved by using equation 2 is attributed to the addition of {12φ} mode, which on the other hand introduces more fluctuations along the φ dimension.It is worth noting however that the three lowest energy facets, [111], [01 ̅ 1], and [1 ̅ 1 ̅ 2], are the points with the highest error because of the nature of the functions which impose sinusoidal periodicity rather than discontinuous first derivative cusps, which would make it difficult to calculate the stiffness.
By examining the interfacial configurations depicted in figure 3 and the Σ3 GB energy map shown in figure 4, a correlation can be established between them.Specifically, we observe that the final facets perpendicular to [1 ̅ 1 ̅ 2] and [111], which are the most stable during annealing, correspond to two local energy minima of the free energy function.However, the distinct difference in evolution between [1 ̅ 1 ̅ 2]-[111] segments and [01 ̅ 1]-[111] segments, as seen in figure 3, cannot be attributed to energy alone, as the energy distribution trends are similar (see figure 4e and f).[15].Thus, assuming similar aluminum behaviors, one could rationalize the final grain configuration based on mobility alone, which may be a factor because of the use of the ECO driving force.However, in the absence of a driving force, Abdeljawad et al. showed that the structure evolved first into a hexagram shape and then later into irregular facets, which have {112} character [6].
The GB stiffness tensor can easily be calculated based on the aforementioned energy fits.As negative stiffness is known to cause interfacial instabilities [6], the distribution of stiffness can significantly impact the dynamics of the spherical interface.Figure 5 plots the calculated principal stiffness values for aluminum Σ3 GBs in the polar θ and azimuthal φ directions.The plot indicates that negative stiffness values are prevalent across a wide range of Σ3 inclinations.The accuracy of the stiffness values depends on the fitting equation, as the formula contains second derivative terms that are highly sensitive to the equation used.In our case, both fitting equations (equations 1 and 2) satisfactorily represent the GB energy and give constant stiffness values along θ (figure 5a and b).However, the introduction of {12φ} mode in equation 2 induces greater fluctuations in stiffness distribution along φ (figure 5d) compared to equation 1. Notably, equation 1 predicts negative stiffness for the [01 ̅ 1] segment, while equation 2 gives a positive value.The negative stiffness based on equation 1 aligns more consistently with the [01 ̅ 1] segment dynamics shown in figure 3, but since both equations provide a reasonable fit to the energy, it is not clear how best to interpret the stiffness in the present case.Additional focus may need to be placed on obtaining an energy function that more accurately captures the cusps of the energy, which might then lead to more accurate calculations of the stiffness in order to better understand its role in the evolution of the structure.It should also be recognized that the spherical bicrystal model simulates a specific subset of neighbor correlations, which correspond to the limited curved segments of the initial sphere.The discretization of the segments is thus influenced by the size of the initial sphere.Nevertheless, this type of simulation allows analysis of the grain boundary plane condition during the whole range from the initial state to the steady-state, which would be of interest to do in the future with better stiffness and mobility data and can not be done based on experimental data due to lack of temporal resolution.

Summary
Previous numerical studies of grain boundary (GB) dynamics have mainly focused on atomistic planar boundary models.However, the limited simulation scale of these models has resulted in less attention being paid to the study of curved GBs.To address this gap, we employ an approach that utilizes the spherical boundary model coupled with the synthetic driving force method to study the dynamics of curved GBs.This approach enables simulations that take into account the effects of interfacial energy, stiffness, and mobility.As a first step into this new area of research, the current provides an energy and stiffness analysis of curved Σ3 GBs.The dynamics of [111] and [1 ̅ 1 ̅ 2] boundary segments are found to follow the energy and stiffness trends, while others, e.g., the [01 ̅ 1] segment, do not.We acknowledge that GB mobility has not been considered and could be another key factor.Effects hereof will be explored in future works.

Figure 4 .
Figure 4.The fitted energy as a function of the BP inclination for Σ3 GBs using harmonic equations 1 (a) and 2 (b).(c, d) Comparison of GB energies between previous atomistic data and the present fitted data corresponding to (a) and (b).(e) Energies of [1 ̅ 1 ̅ 2]-[111] orientations and (f) Energies of [01 ̅ 1]-[111] orientations.One could also examine the mobility trends of Σ3 inclinations as investigated in other FCC materials.Priedeman et al. showed that in Ni, the [111] facet is immobile, the [1 ̅ 1 ̅ 2] facet migrates slowly while the [01 ̅ 1] facet migrates very quickly[15].Thus, assuming similar aluminum behaviors, one could rationalize the final grain configuration based on mobility alone, which may be a factor because of the use of the ECO driving force.However, in the absence of a driving force, Abdeljawad et al. showed that the structure evolved first into a hexagram shape and then later into irregular facets, which have {112} character[6].The GB stiffness tensor can easily be calculated based on the aforementioned energy fits.As negative stiffness is known to cause interfacial instabilities[6], the distribution of stiffness can significantly impact the dynamics of the spherical interface.Figure5plots the calculated principal stiffness values for aluminum Σ3 GBs in the polar θ and azimuthal φ directions.The plot indicates that negative stiffness values are prevalent across a wide range of Σ3 inclinations.The accuracy of the stiffness values depends on the fitting equation, as the formula contains second derivative terms that are highly sensitive to the equation used.In our case, both fitting equations (equations 1 and 2) satisfactorily represent the GB energy and give constant stiffness values along θ (figure5a and b).However, the introduction of {12φ} mode in equation 2 induces greater fluctuations in stiffness distribution along φ (figure5d) compared to equation 1. Notably, equation 1 predicts negative stiffness for the [01 ̅ 1] segment, while equation 2 gives a positive value.The negative stiffness based on equation 1 aligns more consistently with the [01 ̅ 1] segment dynamics shown in figure3, but since both equations provide a reasonable fit to the energy, it is not clear how best to interpret the stiffness in the present case.Additional focus may need to be placed on obtaining an energy function that more accurately captures the cusps of the energy, which might then lead to more accurate calculations of the stiffness in order to better understand its role in the evolution of the structure.It should also be recognized that the spherical bicrystal model simulates a specific subset of neighbor correlations, which correspond to the limited curved segments of the initial sphere.The discretization of the segments is thus influenced by the size of the initial sphere.Nevertheless, this type of simulation allows analysis of the grain boundary plane condition during the whole range from the