Theoretical basis for phase field modeling of polycrystalline grain growth using a spherical-Gaussian-based 5-D computational approach

Using a previously developed phase field modeling method, where interface energies are described by spherical gaussians that allow the modeling of complex anisotropies, a new phase field model was developed to model 5-D anisotropy in polycrystalline grain growth. We present the use of quaternions, assigned to individual grains as orientations and misorientations for grain boundaries, as a means of simulating the ongoing mesoscale changes during anisotropic polycrystalline grain growth. The full 5-D landscape is scanned in MATLAB, and the grain boundary (GB) energy of each grain boundary is calculated from the continuous function developed by Bulatov et al. MATLAB is then used to find all local minima in the GB energy which are stored for use in the phase field model. The methodology of including these minima in the phase field model involves using 2-D gaussian switches, which match the misorientation between grains with misorientations for the GB energy minima. Within a threshold range of the minima misorientation, the switch activates a spherical Gaussian to set the GB energy to the desired value creating in combination a full 5D GB energy space. This creates a GB energy that morphs in real time and space as the GB plane or grain orientations change. Implementation methods of the model are outlined for the Multiphysics Object Oriented Simulation Environment (MOOSE), where reduced order parameters still retain individual grain identification useful for individually assigned quaternions.

three degrees represent the misorientation between the lattices of the two grains, and the other two degrees are associated with the interface plane between the grains. A complex 5-D anisotropy is therefore involved in the relationship between material processing, structure, and properties.
An effective method of simulating evolutions at the mesoscale involves using phase field modeling (PFM). The application space of PFM is wide including solidification (Warren et al. 2003;Gránásy et al. 2004a), solid state phase transformations (Bair et al. 2017;Heo et al. 2019), and grain growth (Moelans et al. 2008a, b;Admal et al. 2019). In the case of grain growth, one PFM method uses non-conserved order parameters individually assigned to each grain represent continuous functions of that evolve with space and time (Chen and Yang 1994). Although highly efficient, PFM has limitations in modeling materials with complex anisotropy, directly represented by direction-dependent interface energies, instead of energies variation between different grain combinations (Hirschhorn et al. 2019;Greenquist et al. 2020;Bair et al. 2020).
In this work, we expand on the use of spherical gaussians in adding complex directional anisotropy to phase field models. This is based on work by Bair et al. (2020), which proved efficient in the case of comparative dendritic growth (Kobayashi 1993). This work outlines the implementation of spherical gaussians to simulate anisotropies in polycrystalline grain growth using PFM. The high flexibility of spherical gaussians allow the application of anisotropy in any desired directions, shaping the anisotropy into any possible complex form (Bair et al. 2020). 2-D gaussian switches, are used to turn on and off specific 3-D anisotropic parameters creating a GB energy that morphs in real time and space as the GB plane or grain orientations change.
The data used to generate the expanded PFM are GB energy minima. The 5-D continuous function derived by Bulatov et al. (2014) is too computationally expensive to implement directly into a PFM, so a MATLAB code is used to scan the full 5-D space of GB's by simulating all the theoretical GBs in a sample. The GB energy is then calculated using the 5-D continuous function developed by Bulatov et al. (2014). Once all the energies are stored, the energy array is scanned for all the local minima. These GB energy values, along with their associated GB geometry (GB plane orientation and misorientation between grains), are the basis for the spherical gaussian and expanded PFM. This approach assumes that the key driving force for anisotropic polycrystalline grain growth will be the minima in energies, and that while some aspects of the continuous function will be lost, the overall effects will be similar.

Orientations and misorientations
A clear understanding of the representation of individual grains and grain boundaries is necessary before proceeding to the development of the phase field model for polycrystalline grain growth. Specific definitions of the ideas of orientations and misorientations are henceforth needed. Whereas orientations are defined as passive rotations which transfer any point in a crystal (grain) reference system to that of a specimen (sample) reference system, misorientations represent the passive rotations from one crystal (grain A) reference frame to another crystal (grain B) reference frame (Krakow et al. 2017). The passive rotations in both cases can be represented by various methods including Euler angle representation, rotation matrices, quaternions, and axis-angle pairs. Quaternions will be used for the representation of both orientations and misorientations in the grain growth PFM considering their computational efficiency and their simpler expressions. Modified versions of such expressions are presented as follows (Horn 1987): where q 0 is the scalar component; and q a the vector component of the unit quaternionq; ⇀ a = (a 1 , a 2 , a 3 ) is the unit vector of the axis of rotation; and θ is the angle of rotation.
Considering the case of a specific grain boundary between two grains A and B, the quaternion representation of the system is made up by representing the orientation of each grain with respect to the sample reference frame R as q A − R (orientation of grain A in reference frame of R), q B − R (orientation of grain B in reference frame of R), and the misorientation between grain A and grain B as q B − A (misorientation from grain B to grain A in reference frame of A). The misorientation q B − A can be determined via quaternion multiplication using either of the following expressions: or with: where q −1 A−R is the inverse (reciprocal); q A−R the conjugate; and ‖q A − R ‖ the norm of q A − R . Since q A − R is defined as a unit quaternion, ‖q A − R ‖ = 1, and q −1 A−R = q A−R . Two different expressions are defined for q B − A due to quaternion multiplication being non-commutative. The two expressions are not equivalent, but conjugate of each other; corresponding to rotations by same angles around axes vectors of opposite directions (Dorr et al. 2010). An arbitrary convention will need to be specified based on the comparative quaternions from the library of misorientations, to be used during the phase field simulation for polycrystalline grain growth. Further explanations of this aspect are presented in the phase field model development. A representation of the orientation and misorientation is given in Fig. 1. (1) q = (q 0 , q a ) (2) q 0 = cos θ 2 (3) q a = q a 1 , q a 2 , q a 3 = a 1 sin θ 2 , a 2 sin θ 2 , a 3 sin θ 2 Another representation of GB geometry is the axis-angle pair (Krakow et al. 2017). This form is useful because it allows for direct and intuitive interpretation from grain orientation and GB misorientation to their mathematical description. The axis-angle is represented by a vector in the direction of the axis of rotation with a magnitude equal to the rotation (misorientation) angle.

Building the library of GBs
The proposed driving force for the expanded PFM are the GB energy minima in the 5-D anisotropic space. Energies are calculated using the MATLAB function developed by Bulatov et al. (2014) and available in their supplemental data. The function requires two inputs in the form of rotation matrices that represent two separate grains in the space. A MATLAB simulation has been developed to generate the rotation matrices and feed them into the 5-D continuous function developed by Bulatov et al. (2014), calculate the GB energies, and output the energy minima along with their corresponding quaternions.
The simulation creates two separate grains (grains A and B) in the 5-D space. Grains A are generated by inputting their axis-angle pairs iteratively. Using axis elements ranging from zero to three and angles from zero to 180 degrees, 51,120 unique grain As are produced after normalizing the rotation axes and removing duplicates. To generate the set of grain Bs, each grain A is rotated about 284 axes by angles five through 180. In total, this generates a space of approximately 2.54 × 10 9 unique GBs. The energy for each boundary is calculated using Bulatov et al. 's continuous function (Bulatov et al. 2014), then all local minima are determined using MATLAB's local minima function with a minimum prominence value specified. When finding local minima, the prominence is the smaller of the two vertical distances measured between the local extreme and its two neighboring, opposing extrema. Specifying a minimum prominence introduces a requirement in the algorithm that local minima must be at least a certain amplitude from their neighboring maxima to be included in the list. Figure 2 demonstrates the usage of the minimum prominence parameter. Without a minimum prominence specified, the minimizing function found 38 local minima along the rotation axis. In this case, a minimum prominence set at 0.18 eliminated the noise and only included two (the lowest 5%) of the total local minima.
To verify that the code generates and correctly feeds real grain boundaries to the GB energy code, sample geometries were analyzed and verified with molecular dynamics data from the work by Olmsted et al. (2009). Eight of the minima found in the generated 5D space come close to coherent twin boundaries, and four of them are nearly exactly coherent twins. The discrepancy can be attributed to the fact that the space was not covered densely enough; increasing the misorientation angles in smaller increments should allow the generation of exact coherent twin boundaries. Despite the absence of perfect coherent twin boundaries, the existing space encompassed boundaries that gave an energy within 10 % of that reported by Olmstead et al. for coherent twins. It could also be seen that, as the boundaries generated moved closer towards a coherent twin, the energies moved closer to the expected value. The GB energy minima for a Σ3 misorientation were plotted and are shown in Fig. 3. There are 2 minima for each of 2 coherent twin boundaries, one boundary plane normal at [111] and one at 111 . By using the minimum prominence parameter when calculating the energy minima, we only account for deep minima. We consider that the effect when computing the gaussian will be highly more significant compared to if shallower minima are included. The cutoff for our minima ends up being 55% of the maximum energy along each rotation axis when using a minimum prominence of 0.18. In the future, a smaller minimum prominence requirement can be employed to include a larger set of minima if desired.

Spherical gaussian anisotropy
The development of the 5-D anisotropy parameter is based on a previous work using spherical gaussians to create complex anisotropies in interface energies (Bair et al. 2020). The study involved using spherical Gaussians with multiple differing minima in any combination of directions to properly represent the anisotropy for simple phase field modeling in 3-D interfaces. It provides a high flexibility in modeling the anisotropy as shown in Fig. 4. The gaussian equation is as follows: with where G u is the gaussian in direction u; n is the total number of gaussians to be added; a u is the amplitude of the gaussian; λ u is a parameter controlling the width (sharpness) of the gaussian; μ u is the unit vector to the center of the gaussian; and υ is the normalized gradient vector of the order parameter η.
A comparative investigation of dendritic growth to the work by Kobayashi (1993) showed the high adaptivity of the methodology. Further analysis revealed that adding anisotropy to the bulk or the interface yielded major differences in the growth process; hence the need to limit it to the interface. This interpretation will therefore be considered in the broadly developed grain growth model. Additional considerations will also be given to the sharpness parameter λ u considering the results indicating its significant effect on dendrite morphology in materials with low anisotropy, which may require further investigations for the cases of polycrystalline grain growth.

Gaussian switch
As stated, the previously developed method is used as base to create a 5-D anisotropy parameter to be applied at the interface for the phase field modeling of polycrystalline grain growth. The 5-D energy is a combination of a 3-D bulk GB misorientation parameters and 2-D GB plane. The methodology involves using 2-D gaussian switches to turn on and off specific 3-D anisotropic parameters, creating a GB energy that morphs in real time and space as the GB plane or grain orientations change. The switch provides a comparison of any misorientation, from a laboratory sample, between two grains q B − A to misorientations form the library q B − Alib of misorientations generated using the work by Bulatov et al. (2014). The library misorientations are misorientations taken by assuming that only some threshold value for low GB energy is needed to produce accurate overall grain growth (e.g., lowest 50% of GB minima). The threshold is to be adjusted as phase field simulations of experimental data are ongoing. The directions, μ uB − Alib , in which the gaussian must be added or subtracted, are provided for each library misorientation q B − Alib . Whenever there is significant relation between q B − Alib and q B − A , the switch will act as a scaling factor for the gaussian anisotropy; whereas in the case of no relation, a base value is assumed.
Considering any misorientation q B−A = q B−A 0 , q B−A 1 , q B−A 2 , q B−A 3 and any library misorientation q B−A lib = q B−A 0 lib , q B−A 1lib , q B−A 2 lib , q B−A 3lib , the switch function is as follows: where Ƨ is the switch function; θ B − A is the angle of rotation for misorientation q B − A ; θ B−A lib is the angle of rotation for library misorientation q B−A lib ; θ ∆θ -Angle variance, is the angle difference of both misorientations; β is the acceptable range of angle from θ B−A lib ; − → V B−A is the unit vector of the axis of rotation for misorientation q B − A ; − → V B−A lib is the unit vector of the axis of rotation for library misorientation q B−A lib ; θ ∆V -Axis variance, is the angle between both unit vectors; α is the acceptable range of angle from − → V B−A lib .
Each grain is assigned a bulk orientation, or the orientations can be directly taken from experimental data. As previously discussed, quaternions, defined in the sample reference (11) frame R, are used for orientation representation. The misorientation between grain A and grain B, q B − A , is used to represent the on-off switch. The on-off switch equals 1 for a specific library misorientation q B − Alib (q B − A = q B − Alib ). Transitioning values (0 < x < 1) will be generated for nearby, somewhat equivalent misorientations (q B − A ≈ q B − Alib ), and 0 elsewhere for completely different misorientations (q B − A ≠ q B − Alib ). The effect of the switch with varying parameters is presented in Fig. 5.
The transitioning values are generated considering that no strict rotational effects exist in real life. A library misorientation q B − Alib representing a 30 o rotation around the (1,0,0) axis will not be very far off energetically from a laboratory misorientation q B − A representing 29.5 o rotation around the same axis. A laboratory misorientation rotational unit vector, separated by a 0.5 o to a library misorientation rotational unit vector, will also not be that far off in effect. The effect might be dampened, but not necessarily null as can be seen by looking at any continuous GB energy representation. Figure 6 shows the variances in possible comparison outcomes between q B − A and q B−A lib ; inducing transitioning values (0 < x < 1).

5-D parameter
Using the spherical gaussian anisotropy equation Eq. (9) and the developed gaussian switch Eq. (11), the 5-D anisotropy parameter between two grains A and B in a laboratory sample of reference frame R can be generated as follows: with: where E AB is the 5-D anisotropy parameter; E base is the base value of the gradient energy coefficient to which the gaussian is added; G lib is the gaussian for a library misorientation;; a ulib is the amplitude of the gaussian; λ ulib is a parameter controlling the width (sharpness) of the gaussian; μ uB − Alib is the unit vector to the center of the gaussian; µ uB−A lib R is μ uB − Alib in the reference frame of R; q A − R is the orientation of grain A in reference frame of R; q B − R is orientation of grain B in reference frame of R; m is the total number of library misorientations; n is the total number of gaussians to be added for a specific library misorientation; v AB is the normalized gradient vector of the order parameters that specifies the orientation of the normal to the grain boundary between grains A and B, and is perpendicular to the contour at which the phase field profiles intersect.
where p is the total number of non-conserved order parameters; E ij is the 5-D anisotropy parameter presented in Eq. (18) for two grains of orientation i and j; E(θ, v) is the overall gradient energy coefficient dependent on misorientation between adjacent grains θ and the inclination of the grain boundaries v; η i and η j are non-conserved order parameters. At a diffuse grain boundary with changes for η i and η j , E(θ, v) = E ij . This formulation accounts for each grain contributing to the grain boundary; hence it accounts for both double junctions and triple junctions where there is a smooth transition. The proposed free energy functional, based works by Chen et al. (1994) and Moelans et al. (2008a), is as followed: with where F is the total free energy of the system; η i and η j are the non-conserved order parameters; m is a driving force parameter; f bulk is the bulk free energy; φ, ω, and γ are phenomenological parameters. φ and ω are usually equal to 1, but variations may occur. Moelans et al. (2008a) provide means of calculating the GB energy and GB width. Modified versions are as follows: where σ ij is the misorientation and inclination dependent specific GB energy between grains with orientations i and j; l ij is a measure for the width of the diffuse grain boundary region; f bulk, interf is the interface bulk free energy dependent on γ ij ; g ij is a numerically calculated function values dependent on γ ij , which a graph is provided for determination. γ ij is to be selected as to minimize the variations of the grain boundary widths in an acceptable range to be determined after further simulations. If necessary, it can be given anisotropy similar to E ij to produce more uniform interface widths. Moelans et al. (2008a) propose 1.5 as an acceptable value; however, variations may still occur. The time evolution of the non-conserved order parameter is controlled using the Allen-Cahn equation: where L is the constant grain boundary mobility.
The functional derivative in Eq. (27) can be expanded as follows: The derivative of E(θ, v) with respect to ∇η i can further be determined using Eqs. (18) and (22). As previously noted, E(θ, v) = E ij at diffuse grain boundary. The formulation of Eq. (22) just provides a continuous function for the phase field modeling considering all types of junctions. The summation term " p i=1 p j>1 η 2 i η 2 j " is first evaluated at each quadrature point. At locations where there are no interaction between grains ( p i=1 p j>1 η 2 i η 2 j = 0 , E(θ, v) = 0 to avoid any undefined division by zero. The sample derivative of E ij with respect to the gradient of the order parameter in a direction, d, can be approximated as followed: hence: Further simplifying Eq. (28) yields the following equation: In a future study, phase field simulations will be performed using Eq. (31) in the Multiphysics Object Oriented Simulation Environment (MOOSE) to prove the theory proposed in this work. The necessary physics are already available in MOOSE. Spherical gaussians kernels and material class developed in previous works (Bair et al. 2020) will be modified to develop a new material class and kernel in MOOSE to simulate Eq. (31) for polycrystalline grain growth.
The model was specifically developed for a finite element differential solver application through MOOSE. Investigation of different forms of differential solver (e.g., semi-implicit Fourier spectral method) would be reserved for later works when successful application of the finite-element scheme yields accurate enough results or when determined that simpler modifications may be more advantageous than opting for a new method.
It is noteworthy that this method of adding anisotropy could be modified to fit into the orientation phase field models (Pusztai et al. 2005(Pusztai et al. , 2008Gránásy et al. 2004b;Korbuly et al. 2017a, b;Henry et al. 2012). Pusztai et al. already proposed the possibility of including the full 5D GB energy using the 4 values of quaternions for each grain to make 4 orientation field variables in a simulation (Pusztai et al. 2005). This method could be an attractive alternative as it would reduce the total number of variables necessary to model the system. However, it would require significant modification of the form of anisotropic equations used in this work. Both should be developed further and tested to compare the benefits of each.

Conclusions
A theoretical phase field model for polycrystalline gain growth was developed using a spherical gaussian method to add 5-D anisotropy to the grain boundary energy. The validity, performance, and possible improvements of the model will need to be determined by performing simulations using MOOSE, with comparison to experimental studies.
With a determined minima set acquired by verified methodology, one could apply the gaussian method to specify anisotropy. Bulatov et al. (2014) continuous function for grain boundary energy for nickel and different materials is used in this sense by creating numerous boundaries that incorporate the whole of the 5D space using quaternions, calculating the energies, and finally finding the minima set while considering the cut-off energy percentage by using MATLAB's minimum prominence parameter. From there, the gaussian model with the on/off switch would be able to look at any misorientation and determine whether to subtract the energy in case of a minima (% subtraction with the continuous switch) or keep the maximum energy as default value. Similar methods as those described in this work could be applied to add anisotropy to GB mobility in PFM. Several studies have shown significant changes in mobility with temperature and crystallography, which will likely have a strong effect on the morphology of the microstructure. However, due to the current lack of a continuous function for GB mobility, this addition was beyond the scope of the current study.