Magnetic APFC modeling and the influence of magneto-structural interactions on grain shrinkage

We derive the amplitude expansion for a phase-field-crystal (APFC) model that captures the basic physics of magneto-structural interactions. The symmetry breaking due to magnetization is demonstrated, and the characterization of the magnetic anisotropy for a BCC crystal is provided. This model enables a convenient coarse-grained description of crystalline structures, in particular when considering the features of the APFC model combined with numerical methods featuring inhomogeneous spatial resolution. This is shown by addressing the shrinkage of a spherical grain within a matrix, chosen as a prototypical system to demonstrate the influence of different magnetizations. These simulations serve as a proof of concept for the modeling of manipulation of dislocation networks and microstructures in ferromagnetic materials within the APFC model.


Introduction
External magnetic fields offer the possibility to manipulate microstructure in ferromagnetic materials. This has several potential applications in microstructural engineering [1]. However, a detailed understanding of the interactions between magnetic fields and solid-state matter transport is unavailable. This is due to the lack of a theoretical model which allows describing the fundamental physics of magneto-structural interactions in a multiscale approach, combining the dynamics of defects, dislocation networks, and grain boundaries with experimentally accessible microstructure evolution on diffusive time scales.
A promising approach in this direction is provided by a recently introduced phase-field-crystal (PFC) model that captures the basic physics of magnetocrystalline interactions [2,3]. The PFC approach [4,5] is a continuum method that has been shown in numerous publications to capture the essential physics of atomic-scale elastic and plastic effects that accompany diffusive phase transformations, such as solidification, dislocation kinetics, and solid-state precipitation, see [6] for a review. In [2] the PFC density is coupled with magnetization to generate a ferromagnetic solid below a critical temperature (the Curie temperature). The latter can depend on the local mean density. In [3] this PFC model is extended to multiferroic binary solid solutions and used to demonstrate the influence of magnetic fields on the growth of crystalline grains. The model, which is a system of evolution equations for the rescaled particle density field ϕ and an averaged magnetization m, is used in [7] together with efficient and scalable numerical algorithms [8] to study the role played by external magnetic fields on the evolution of defect structures and grain boundaries. The induced magnetic anisotropy has consequences also on larger scales. They result from the dependence of the bulk energy of single crystals on the direction of the magnetization w.r.t. to crystallographic orientations. Indeed, simulations show that the application of external magnetic fields affects the growth and coarsening of grains in microstructures, leading to the selection of preferential orientations. These results are in qualitative agreement with experiments, e.g., on Zn and Ti sheets [9], and classical grain growth simulations of Mullins type with an analytical magnetic driving force [10]. The additional driving force, due to the external magnetic field, also enhances the coarsening process in the simulations [7]. This effect has been observed experimentally, e.g., during annealing of FeCo under high steady magnetic fields [11]. It has also been shown that the grain boundary mobility is anisotropic with respect to the applied magnetic field. This kinetic effect leads to elongated grains. The simulations show that under the influence of a strong external magnetic field, the magnetization m can be assumed to be constant. This simplifies the considered model and allows the simulation of larger systems, thus enabling the analysis of long-time scaling behaviors and various geometrical and topological properties in grain growth. In [12] it is investigated how the scaling and characteristic geometric and topological properties change under the influence of strong magnetic fields. The results are in qualitative agreement with experiments for thin Zr sheets annealed under the influence of magnetic fields [13]. However, even with these fundamental contributions, simulation-driven applications in microstructural engineering are still out of reach. While the microscopic details are well resolved with the considered PFC model and also experimental relevant time scales can be reached, the required spatial resolution essentially restricts physically relevant simulation to two-dimensional settings or relatively small samples in three dimensions, see, e.g., [14,15].
Required for a true multiscale approach is a coarse-graining in space. The complex amplitude PFC (APFC) model originally introduced in [16,17] provides such a framework. The idea of the APFC approach is to model the amplitude of the density fluctuations instead of the density itself. This enables to reach larger spatial scales by still retaining essential microscopic effects, thus enabling mesoscale investigations of crystalline systems [18][19][20][21][22]. For a recent review of APFC models, we refer to [23].
We here extend the magnetic PFC model [2,3], or to be more precise, the simplified version used in [12], to the APFC framework. In Section 2 we formulate and explain the model. The symmetry breaking due to magnetization is discussed in Section 3. In Section 4 we report the main concept behind the numerical approach exploited to solve the coupled system of equations of the magnetic APFC model, and in Section 5 we demonstrate the applicability for a simple three-dimensional setting. We consider a spherical grain in a matrix and simulate its evolution under the influence of a strong magnetic field. The results are compared with PFC and APFC simulations without magnetic field [15,24]. Finally we draw conclusions in Section 7.

From magnetic PFC to magnetic APFC models
We consider here the limit of strong external magnetic fields. In this limit, the local magnetization is assumed to be constant in the crystal. The direction and the magnitude of the magnetization is solely defined by the external magnetic field. Then the magnetic PFC model in [2,3] in terms of its free energy functional, reduces to: where ϕ denotes the scaled particle density and m the magnetization. q 0 defines the lattice spacing at equilibrium, which we set here to 1 without loss of generality. Ω is the domain of integration. B x 0 , ∆B 0 , τ and v are parameters as introduced in reference [25]. Together with the average densityφ they define crystal structure and physical properties. Parametric studies and phase diagrams obtained by varying such parameters in classical formulations of the PFC model can be found in the literature [5,[25][26][27]. α 2 controls the strength of the magnetic interaction. The considered form α 2 ϕ(m∇) 2 ϕ provides the simplest possible term leading to symmetry breaking by the magnetization, m, which is assumed to be constant. The magnetization is scaled to unit length, |m| = 1 without loss of generality. Then, the evolution equation for conserved dynamics reads Note that this consists of a 6th order partial differential equation. We chose model parameters such to describe a BCC crystal at equilibrium, similarly to references [28,29]. They are reported in Table 1. The magnetic coupling is controlled by α 2 , which is specified below. In the APFC model the particle density ϕ is expressed in terms of the Fourier modes of a relaxed reference crystal with complex amplitudes A j (r) accounting for its deformation, where A * j (r) is the complex conjugate of A j (r) in order to describe the density wave. The reference crystal is then defined by a set of N reciprocal space vectors k j . The set {k j } explicitly entering the APFC model is typically chosen by considering the minimum number of modes describing the crystal, i.e. the shortest k j , with −k j accounted for through the complex conjugate in eq. (4) [23], see also Figure 1. A single crystal, corresponding to the reference structure without deformations, is then described by real and constant amplitudes. For deformed crystals, the set of reciprocal-space vectors, {k j }, is modified by a constant deformation D: With respect to the original description, eq. (4), this leads to A j (r) → A j (r)e ik j (D −1 −1)r . Thus, a constant deformation leads to spatial varying complex amplitudes. The gradient of the phases of the amplitudes, k j (D −1 − 1), define their periodicity. The same holds also for rotations. Further on, the explicit space dependence for ϕ(r) and A j (r) will be omitted.
To derive the corresponding magnetic APFC model, eq. (4) is substituted in the PFC energy including the effect of magnetization, eq. (1). If the amplitudes, A j , vary on larger length scales than the density wave, ϕ, then the fluctuations on small scales can be averaged [30]. This can be formally done by multi-scale analysis [28,31], or renormalization techniques [31,32]. All techniques show that only resonant terms contribute to the energy in this limit. That is, only terms, where the phases of the reference density wave, k j · r, cancel each other, remain: e.g. A j A * j = A j e ik j ·r A * j e −ik j ·r are the only quadratic resonant terms. Following these concepts, the magnetic APFC energy reads with . A convenient choice of the reference crystal is the one corresponding to an expansion as in eq. (4) which minimizes the free energy, eq. (5). It corresponds to |k j |= 1 with constant and real amplitude, whose values are dependent on model parameters [23]. f S is a polynomial function in {A j } and {A * j }. It depends on the τ vφ B x 0 ∆B 0 1 1 0 0.98 0.02 Table 1. Modeling parameters for magnetic PFC and APFC model [28,29].
reference crystal structure, see [23,28,29]. An example of f S is given in Section 3 for a BCC crystal. The energy reported above, eq. (5), can be split in three characteristic parts: encodes the dependence of the energy on amplitudes' phases through their gradients, and thus on lattice deformation with respect to the reference crystal. As such, it leads to the elastic properties of the model [5,20].
Thus, it does not change with deformation of the crystal. It contains the lattice , independent of the crystal structure. These terms determine the absolute value of the amplitudes.
iii) The magnetic coupling term, N j=1 α 2 A * j M 2 j A j , accounts for the magnetization m. It depends on the relative orientation of k j and m. Thus, as it will be also shown in the following, it breaks the rotational symmetry of the free energy, leads to magnetic anisotropy, and introduces magneto-striction.
Within the long wavelength limit [33], the approximation of the conserved dynamics for the particle density leads to a non-conserved dynamics for amplitudes. The evolution equations for each amplitude reads with Thus, the amplitude expansion leads to a 4th order partial differential equation for every amplitude. The number of coupled equations depends on the considered crystal structure, namely on the number of Fourier modes to represent the targeted lattice. Details on how to generally solve this equations numerically including defects and grain boundaries are reported in Section 4. However, in order to analyze the magnetostructural interaction, we first consider a defect free crystal with constant deformation. This leads to a simpler formulation used to analyse magnetic-induced anisotropies in the next section.

Symmetry breaking due to magnetization
In this section we analyze the proposed magneto-structural coupling in (A)PFC for a single crystal without defects and restrict deformations to be constant in space. In this case, the density, eq. (4), can be expanded according to deformed reciprocal space vectors, k j = D k k j . The deformation matrix D k is independent on k j and is connected to the spatial deformation of the crystal, D, by D k = (D −1 ) T . With these assumptions, the free energy, eq. (5), reduces to where |Ω| is the volume of the integration domain. We consider a BCC crystal (S = BCC), which is described by a set of N=6 k j 's, see Figure 1.
It defines the size of the unit cell of the crystal without any magnetic coupling. The second term, N j A j (m·k j ) 2 A j , couples local magnetization to the local structure in the crystal. The structure dependent part of the energy reads We recall that this term is independent on deformation. Due to symmetry considerations we only allow for deformations parallel and perpendicular to m. Then the deformation can be written as D = R TD−1 R withD a diagonal matrix.   magnetic coupling depends on (m · k j ) 2 and thus on the angle of m and the reciprocal space vectors describing the crystal, k j . Thus, the magnetic coupling acts differently on the single amplitudes. For m parallel to [1 1 1], there are two sets of equivalent amplitudes: those corresponding to {k 1 , k 2 , k 3 }, which have the same angle with respect to m and those corresponding to {k 4 , k 5 , k 6 } which are perpendicular to m. That is, the magnetic coupling breaks the symmetry of the energy w.r.t. the amplitudes. This symmetry breaking is shown in Figure 3. For vanishing coupling, α 2 = 0, all amplitudes are the same. Increasing the coupling and α 2 > 0 leads to increasing amplitudes {A 1 , A 2 , A 3 } and decreasing amplitudes {A 4 , A 5 , A 6 }. The model breaks for α 2 outside the shown range as some amplitudes vanish, and the energy minimizing state does not describe a BCC crystal anymore. The magnetic coupling for α 2 > 0 also leads to a compression of the crystal in the direction of m and a smaller contraction perpendicular to m. The total energy decreases with increasing magnetic coupling α 2 , see Figure 3.
The magnetic anisotropy induced by the proposed magnetic coupling can be shown by comparing the energy for different orientations of m. The highest energy for a fixed α 2 is always for m parallel to 1 0 0 , therefore referred to as the hard directions of magnetization in this model. The directions with the lowest energy, namely the easy directions, depend on the sign of α 2 . For α 2 > 0 they are 1 1 0 , while 1 1 1 are the easy directions otherwise.
The general dependency of the energy can be directly explained by the sign of the magnetic coupling. For α 2 < 0 the magnetic coupling term introduces a positive contribution in the energy and, thus, the total energy increases. However, the magnetic anisotropy depends on the relative influence of the magnetic coupling on the different amplitudes. In order to understand the induced magneto-striction, the influence of the magnetic coupling on {k j } has to be examined. The term (1 − (k j ) 2 ) 2 defines the k j in equilibrium without magnetic interaction. That is, this term has a minimum at |k j |= 1. The energy contribution of this term vanishes in this case. The magnetic coupling changes the situation. The minimum of (1 − (k j ) 2 ) 2 − α 2 (m · k j ) 2 depends on the relative orientation of k j and m. For α 2 > 0, the minimum is shifted towards larger k j , and the energy contribution becomes negative. Thus, amplitudes influenced by this term, namely A j for which m · k j = 0, become larger. The shift of the minimum towards larger k j leads to a compression of the crystal in this direction. This analysis is exact for a one-dimensional crystal (so namely for a single k j vector). For a structure as BCC it is more involved as the k j 's do not enter in the energy independently. They all can only transform with the same deformation, D k . Extending these computations to all orientations allows to compute the magnetic anisotropy. The corresponding polar plots for various α 2 are provided in Figure 4.

Numerical approach and setup
We solve eqs. (6) numerically using the finite element method (FEM). At variance with other common approaches exploited in the literature, such as the pseudo-spectral Fourier method, it allows for an optimized spatial discretization for amplitudes that vary unevenly in space. Indeed, they feature significant variations at dislocations, requiring a resolution approaching the one needed by the classical PFC model, but vary slowly for rotated and strained crystals and are constant for relaxed crystalline domains. We consider in particular the approach reported in references [29,34] and adapt it to the novel coupling with magnetic fields proposed in this work. To discretize eqs. (6) we expand the fourth-order operator and rewrite eqs. (6) as systems of second order equations with G j ({A j }) := ∂g S ({A j })/∂A * j the nonlinear terms. We consider a semi-implicit time discretization evaluating linear terms implicitly and providing an approximation for the implicit evaluation of non linear terms. In particular, these terms are linearized for every amplitude via a one-iteration Newton method, such that where A n j and A n+1 j denote A j at timestep n and n + 1, respectively. The considered semi-implicit discretisation then reads where ∆t is the timestep. For the investigations reported in this work, a proper numerical convergence is achieved for ∆t = 1. In order to use standard solvers and real basis functions for the FEM implementation, the above system of equations has to be split further into the real and imaginary parts of the amplitudes. This leads to a system of four coupled second-order equations for each amplitude which can be solved by exploiting linear elements. See references [29,34] for further details. The weak form of these systems of equations and their FEM discretization are implemented in the parallel and adaptive finite element toolbox AMDiS [35,36]. Following reference [34] mesh refinement is set according to the local variation of the amplitudes. In particular, the spatially-dependent wavelength of their oscillation is computed and properly resolved by a fixed number of mesh elements per period of the oscillation. In addition, a minimum refinement is set for regions where the crystal does not exhibit any deformation or rotation, and a fine mesh is set at dislocations (see also figure 5).
We examine the influence of magnetization on the shrinkage of an initially spherical grain rotated with respect to the surrounding matrix, see figure 5, in a BCC crystal. The grain has a radius of 30π and is rotated by 10 • about the [1 0 0] direction. That is, the rotation axis is set as the z-axis in our frame. From a macroscopic point of view, the normal velocity of the grain boundary v can be described by [37,38] where M is the mobility, γ is the grain boundary energy, κ the local curvature and ∆f the difference of the bulk energies of matrix and grain, and neglecting additional elasticity effects [39,40]. In our case ∆f results from the magnetization, as the magnetic properties are different in the grain and the matrix. Without ∆f , eq. (14) reduces to mean curvature flow and describes isotropic grain shrinkage due to minimization of the grain boundary energy. As already shown by PFC [15] and APFC [24] simulations without magnetization, while recovering the general scaling described by eq. (14), the shrinkage is anisotropic. We will therefore directly compare our result for the grain shrinkage under the influence of magnetization with these previous findings. The crystal structure in the matrix surrounding the circular embedded grain is set as the reference structure. Thus, the amplitudes therein are constant and real. Within the uniformly rotated grain, the amplitudes vary regularly, still featuring a constant A 2 . This quantity varies significantly at dislocations only, and can be therefore exploited for their identification [24,29]. Indeed, we recall that dislocations are described in the APFC model as zeros of some of the complex amplitudes, namely the ones which have singular phases [23]. Far away from the dislocation core, they smoothly connect to their bulk values. A diffuse representation of the dislocation core over a few lattice spacing is then achieved [41], consistent with non-singular continuous descriptions [42,43]. Dislocation networks can then be reconstructed as regions where A 2 < A 2 with A 2 a selected threshold. For visualization purposes, the extension of these regions across the dislocations can be chosen arbitrarily by selecting a value of A 2 within the range of A 2 , provided that it ensures enough resolution [24,29]. Here, we set A 2 = 0.073., corresponding to ∼ 2 lattice spacing. Figure 5 illustrates the setup for the simulation of the embedded rotated grain and the resulting small-angle grain boundary, together with the features of the amplitudes, the reconstructed dislocation network after relaxation of the initial condition, and the inhomogeneous mesh refinement used for simulations.

Simulation results
To study the influence of the magnetization on grain shrinkage, three configurations are considered: a simulation without magnetization, α 2 = 0, (M0), corresponding to [15,24], and two settings with magnetization, one that hinders grain shrinkage (M1) and one that enhances grain shrinkage (M2). We consider a coupling strength, α 2 = 0.1, and only vary the direction of magnetization, m, which is chosen according to the magnetic anisotropy to enhance or hinder grain boundary motion. For m parallel to [1 1 1], the energy in the matrix would be larger than in the grain. This results from the grain rotation by θ about [1 0 0] of the matrix crystal. The specific setup can be seen in Figure 5. This energy difference would allow us to consider the effect of the anisotropy induced by the magnetization. However, as the differences in the magnetic anisotropy for this configuration are relatively low (see also Figure 6), we consider a different orientation for m. We rotate m about the [1 0 0] direction by β = π/8 + θ for (M1) and β = −π/8 − θ/2 for (M2), respectively. For positive β (M1), the bulk energy in the grain is lower than in the matrix. Thus, ∆f in eq. (14) becomes positive, which reduces the velocity of the grain boundary. Rotation of m in the other direction (M2) leads to lower bulk energy in the matrix than in the grain. Thus, ∆f becomes negative, which provides an additional driving force in the direction of grain shrinkage and thus enhances grain shrinkage. Figure 7 shows the shrinkage of the grain for the three cases (M0), (M1), and (M2). We consider the area of the grain boundary, computed as in [24], as well as lengths l 0 andl = 1 2 (l 1 + l 2 ), as defined in Figure 8, to account for the anisotropic shrinkage. As expected in the setting of enhanced grain shrinkage (M2), the grain shrinks faster than in the case of hindered grain shrinkage (M1). However, the fastest shrinkage is obtained for the reference case (M0). This behavior has already been observed in [7], and it can be explained by the dependency of M γ, namely the product of the mobility and the grain boundary energy, on m. The anisotropy of the shrinkage is influenced by m but remains within the same order for all three cases. To further highlight details of the grain shrinkage discussed above, dislocation networks corresponding to the open symbols in Figure 7(left) are shown in Figure 8. They correspond to the same grain boundary area and indicate the deformation of the dislocation network by the magnetization. Figure 9 shows the dislocation networks during their shrinkage over time. Representative stages are selected, corresponding to the full symbols in Figure 7.
They also indicate the magnetization-dependent deformation of the dislocation network. Besides the global shapes, the dislocation networks also show a qualitatively different

Discussion
The classical classical APFC model has been proved powerful in describing crystalline structure and defect networks at coarse-grained and mesoscale lengthscales while exploring diffusive timescale similarly to the PFC model. In [24] it has been shown that the dislocation network and its (overdamped) dynamics predicted by APFC corresponds to the one achieved by PFC [15]. More importantly, static properties, symmetries and energetics, are qualitatively as expected by other atomistic models [44] and experiments [45]. It is worth mentioning that, although minimal formulation have been mostly considered, PFC and APFC formulations may be adapted to fit quantitatively grain boundary properties of specific materials [46,47]. In this work, we showed that the influence of strong magnetic fields on the dynamics of grain boundaries could be included in the APFC model. Depending on the orientation of magnetization, the grain growth is hindered or enhanced as predicted by continuum models, eq. (14). The dislocation network is deformed due to the deformation of the crystal and matrix due to magneto-striction. Although the qualitative character of the reported results, extended parameterizations may be devised and extended model may be considered together with the coupling with magnetic field here proposed.
Many relevant magnetic materials are alloys where material properties are dependent on the solute concentration. This give rise to variation of magnetic properties due to variation in solute concentration. Especially at grain boundaries segregation may lead to change of magnetic properties [48][49][50] and, vice-versa, segregation is influenced by magnetization [51,52]. PFC and APFC models for binary alloys and ordered crystals have been already reported in literature [21,25,28,[53][54][55]. It has been demonstrated that these models can track grain boundary segregation [56] and precipitation [53], as well as dislocation dynamics accounting for Cottrell atmospheres [21]. These models are formulated either by treating the density of every species as a density field or by reformulating the model in terms of a mean mass density and a concentration field. However, in both approaches, the magneto-structural interaction can easily be in-cooperated the same way as proposed in eq. (1). The magnetization is coupled to the gradient of density the density field, α 2 ϕ(m∇) 2 ϕ, leading to magnetostriction and magnetic anisotropy if ϕ is a density field representing a crystal. Also, the coupling constant α 2 may be set as dependent on local concentration to account for magneto-structural interaction due to solute concentration at dislocations and grain boundaries. Additional terms which are isotropic but depend on magnetization and concentration may be introduced to account for local variations in the crystal, e.g., m 2 ϕ 2 [2,3]. However, it is worth noting that this ansatz is not expected to be sufficient for an accurate, quantitative description of any crystalline material. Indeed, magnetic anisotropy and magnetostriction are controlled by a single parameter α 2 , and, for BCC, the hard direction in this model is always [1 1 1]. Thus, higher-order terms for comprehensive modeling of magnetic coupling are needed as proposed in [3]. Further research is ongoing to develop a minimal magnetic coupling model that can be adapted to a broader class of materials within the APFC model. Furthermore, to account for energies due to variation in magnetization, a formulation beyond the assumption of constant magnetization [12] is necessary. To this goal, the free energy has then to be complemented by a Landau-Lifshitz-Ginzburg type energy for the magnetization field as done in [2,3].
In most of the extensions and perspectives mentioned above, the magneto-structural properties in the crystal can still be explained by an analysis similar to what has been reported in this work. In particular, the simplest magneto-structural coupling as studied here is the basic building block to introduce anisotropic magnetic interaction in the modeling of binary materials with PFC and, thus, APFC. This paves the way for the development of APFC models including magneto-structural interaction for alloys and, therefore, towards quantitative application of the APFC model to the modeling of specific materials.

Conclusion
We have given a proof of concept for the simulation of the interactions between magnetic fields and dislocation networks in a three-dimensional setting with the APFC model. We derived a magnetic APFC model and adapted an existing, state-of-the-art numerical approach to solve the corresponding equations. As a result, the basic physics of magneto-structural interactions are accounted for in a multiscale approach. We have demonstrated that the considered framework accounts for magnetic anisotropy and provided an example of its influence on the shrinkage of a rotated spherical grain within an unrotated matrix. This consists of a minimal setting that may be extended toward the modeling and simulation of defect networks and microstructure evolution on diffusive time scales and at relatively large length scales. The proposed approach is indeed the basic building block to introduce magneto-structural coupling in PFC models of complex crystals and alloys. However, various open problems remain. One central issue is an appropriate parametrization of the material parameters for specific materials. While this has been done for the bulk modulus and the grain boundary energy, e.g., for Fe [57], the magnetic anisotropy in the model has to be tuned to resemble known functional forms for the magnetic anisotropy of materials of interest.