Multi-phase field modeling and simulation of magnetically driven grain boundary migration in SmCo polycrystals

There is a growing demand for magnetic materials straight forward in wind turbines and electric motors. Their functional properties depend critically on their microstructure and thus on the microstructure evolution during sintering or heat treatment. Field-assisted selective grain growth allows to optimize the microstructure. However, the simultaneous modeling of the structural and magnetic degrees of freedom on the micrometer scale is not possible with most simulation packages. Therefore, we extend the open-source software project OpenPhase and implement the micromagnetic equations needed to treat both degrees of freedom in the framework of the multi-phase field method. We apply our model to the field-assisted grain growth in Sm2Co17 polycrystal films.


Introduction
The mechanic and magnetic properties of polycrystalline materials are crucially dependent on their microstructure, * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
particularly the texture of their grains [1][2][3]. The remanent magnetization and the coercive field can be optimized by the size and alignment of grains [3,4]. However, the change of these elements of microstructure under heat treatment and sintering is not fully understood.
For nonferromagnetic metals, an external magnetic field has been reported to influence solidification and grain growth by selective nucleation and growth of grains with favorable alignment relative to the field [5][6][7]. For example, an external field influences the coarsening of Zn-Al alloys due to the anisotropy of the magnetic susceptibility χ [8]. Even stronger coupling effects between grain growth and an external magnetic field are possible in ferromagnetic (FM) materials with magnetocrystalline anisotropy. For example, it has been shown that an external field influences the grain growth of FM Alnico 8 alloys for a field strength of 0.4 MA m −1 whereas a field strength of 25.5 MA m −1 is needed to change the orientation distribution of the grains in the Zn-Al alloys [8,9]. It is therefore of importance to model the time evolution of the microstructure in an external magnetic field.
A model for field-induced magnetization in bismuth has been developed by Mullins [10]. For two grains (α and β) with a difference of susceptibilities along and perpendicular to, e.g. a trigonal axis (∆χ), the different angles between this axis and the applied field (θ α and θ β ) result in a magnetic driving force on the in-between boundary: where µ 0 is the vacuum permeability. This model was applied by Molodov et al to the grain boundary motion in bismuth bicrystal [11] and has been applied successfully to other nonferromagnetic metals [8,[12][13][14].
In recent years, the model has been introduced in a mesoscopic phase field model and successfully reproduces the field assisted coarsening in a large scaled polycrystal [15]. For FM polycrystal in an external magnetic field, a solution based on a phase field model is provided by Feng et al [16,17]. However, so far a uniform magnetization field has been used in contrast to the spatial distribution of magnetization in larger grains [18] and multi-phase systems cannot be addressed.
In the present work, we set up a general model for polycrystalline materials with and without an FM domain structure. We combine the multi-phase field (MPF) method and the time-dependent Ginzburg-Landau (TDGL) method to treat the evolution of the grain structure and magnetization on an equal footing [19,20]. A similar strategy has already been used successfully to study the martensitic transformation in magnetic shape-memory alloys [21][22][23].
We use the SmCo as a test case. As a prototypical material for high-performance permanent magnets, it combines large spontaneous magnetization, magnetocrystalline anisotropy, and high Curie temperature (up to 1000 K) [24][25][26]. We show that a moderate field (<1T) influences the time-evolution of grains during coarsening.

MPF model
The MPF method as developed by the last author and coworkers is a powerful approach to moving boundary problems in general [19]. It allows us to address the solidification and coarsening of metals [27,28] and is compatible with transport methods in bulk phases, for example, the Lattice Boltzmann method applied to microfluidity [29,30]. Each grain α is represented by one scalar field ϕ α . The value of ϕ is limited to the interval [0, 1]. The field can be understood as the volume fraction of this grain in a given unit volume at position x. Then the summation of phase fields on this point obeys the constraint: The equation of motion of N phase fields is given by the phase field equation: where α and β are the indices of phase fields, M ϕ αβ is the mobility of the interface between the phase fields α and β. Note that we use M ϕ for the interface mobility while M for the magnetization. The energy functional F is defined as the volume integral of all energy densities and typically consists of bulk and interface terms, which can be written as: where f local int (ϕ(x)) is the local interface energy at the point x in the whole volume Ω, which has two contributions: the gradient energy term and the potential energy term. We define the local interface energy as: where σ αβ is the interface energy of the boundary between α and β and η is the numerical interface width. The first term in equation (5) is the gradient energy and the second term is a double-obstacle interface potential. The latter determines a stable smooth interface with a sinus profile and we set the interface thickness to five grid points (η = 5∆x). On the interface between the fields α and β, the driving force for grain growth can be calculated by: ) .
We utilize the TDGL equation to model the time-evolution of the magnetization field: in which L is the kinetic coefficient. Note that the TDGL equation gives the same result as the Landau-Lifshitz equation under overdamping [31], and this allows us to reduce the computational effort as the Larmor precession is not important for our study.

MPF modeling for magnetic energy
On the length scale of interest, the magnetic energy can be divided into the exchange energy, the magnetocrystalline anisotropic energy, the Zeeman energy, and the demagnetizing energy. We take the common form of these energies from [20,32]. The key point of coupling the phase field equation and the TDGL equation is to build an MPF model for the magnetic energy. Therefore we define a linear magnetization model as: The first term on the right-hand side describes the spontaneous magnetization M sp only present in an FM phase, while the second term is zero in an FM material and corresponds to the magnetic susceptibility χ otherwise (M sp = 0, ξ > 0 for paramagnetic, and M sp = 0, ξ < 0 for diamagnetic materials). This model can describe the magnetization field of the whole system and no extra numerical steps are required to distinguish the phase difference. Furthermore, we can keep the M field mathematically continuous and differentiable on the grain boundaries. The local magnetization field in a unit volume is given by: where α is the index of phase fields, M sp is the magnitude of spontaneous magnetization, and m is its orientation. Within each phase field (ϕ α = 1) we assume constant values of M sp and ξ and on interfaces and junction points the local magnetization is the mixture of the involved phase fields. Specifically, if all phases are para-or dia-magnetic, equation (9) will reduce to the so-called Wiedemann's additive law of the magnetic susceptibility [33]. With this model for M local , we can determine the mean value of magnetic energy contribution from all phase fields. The gradient of the local magnetization field determines the local exchange energy density as: where A α is the exchange energy coefficient of ϕ α and M i,j stands for ∂M i /∂x j . Note that the induced magnetization does not contribute to the exchange interaction. The local Zeemann energy, i.e. the coupling between the local magnetization and the applied external field (H ext ), can be written as: ) .
Similarly, the local demagnetizing energy related to the demagnetizing field H d is defined as: In principle, our method allows us to treat arbitrary magnetic anisotropies. We restrict this study to the uniaxial anisotropy which is typical for Sm 2 Co 17 [34] and can be calculated as: in which K is the anisotropic coefficient and θ is the angle between the magnetization direction m and the easy axis (E A ).
Locally the anisotropy is given as: By summing up equations (10)- (12) and (14), and integrating them over volume, we obtain the local magnetic energy functional: (15) Inserting this equation in equation (7), we obtain the local effective field and the time-evolution of the magnetization in FM grains. On the other hand, the magnetization of nonferromagnetic phases is directly induced by the magnetic fields. We insert equation (15) in equation (6) to calculate the driving force for grain boundary motion as: If the magnetic properties do not differ in adjacent grains, there will be no driving force contributed by related terms. For FM systems, the driving force mainly depends on the relative angle between the easy axes of grains and the local magnetization. For non-FM systems (M sp = 0), equation (16) reduces to: This form is similar to the model used to discuss the magnetic driving force in bismuth, which is given by equation (1). We note that the present work focuses on the impact of the magnetic anisotropy on the grain evolution in an external field only. Depending on the material class, for example, the giant magnetostrictive materials [35], the magnetostriction may have an impact on the grain evolution and can be added in the same formalism. However, for the test case of SmCo, the magnetostriction coefficient is small [36] and the anisotropy is probably more relevant.

Simulation details
The discussed magnetic solver is implemented in the OpenPhase software project and we solve the equations using Finite Difference Method. For this purpose, we define the simulation box as a discrete space with equidistant grid spacing ∆x in all directions. The demagnetizing field in equation (12) is obtained by solving Maxwell's equations in Fourier space following [32] and using the FFTW software library [37]. To validate our implementation of the magnetic dynamics, we perform reference simulations with the widely used software Object Oriented MicroMagnetic Framework (OOMMF) [38].
During the simulation, we need to solve two equations of motion: the TDGL equation for the M field and the MPF equation for the ϕ fields. As spin dynamics are typically orders of magnitude faster than grain boundary motion, we decouple these dynamics and relax the M field to a (meta)stable state before we propagate the ϕ fields in time.
As the first test of our implementation, we model the magnetic properties of an idealized isotropic (K = 0) FM thin disk in vacuum. The diameter and thickness of the sample are 93 and 15 grid points, while the dimension of the simulation box is 121 × 121 × 41 and we apply Periodic Boundary Condition (PBC). We set the grid space to ∆x = 10 −8 m, the spontaneous magnetization strength to M sp = 8 × 10 5 A m −1 , and the exchange coefficient to A m = 1.3 × 10 −11 J m −1 .
Sm 2 Co 17 has the hexagonal Th 2 Ni 17 -type structure [24] with uniaxial magnetocrystalline anisotropy along [0001] axis [39]. We take the material parameters from [34,40] at 800 K. This temperature is typical for the annealing of SmCo materials during processing [41,42], it is sufficient to activate grain boundary motion and the material is still in its FM phase. These parameters include spontaneous magnetization: M sp = 7.96 × 10 5 A m −1 , exchange coefficient: A m = 7.8 × 10 −12 J m −1 , and uniaxial magnetocrystalline anisotropy coefficient: As a compromise between accurate sampling of the magnetic structure of the SmCo samples and the efficient modeling of micrometer-sized systems, we set the grid spacing to ∆x = 1.237 × 10 −8 m. We note that this is a rather rough approximation for the subtle details of the domain evolution in small external fields. However, for moderate field strengths, the magnetization fields in the grains are uniform and are well described by this grid (see figures 4(c) and (d)).
Note that we only include temperature by the choice of the coefficients. In addition, in our current work, we do not consider the reduction of the magnetic properties on the grain boundary, because the physical thickness of the grain boundary is on the nanometer scale and much thinner than the numerical grain boundary we apply [43]. Thus we expect a minor influence.

Model verification
We start with the validation of our implementation using a magnetic vortex state in a thin magnetic disk. For this purpose, we initialize the M field randomly. After relaxation, the magnetization vectors form the expected magnetic structure with m pointing along the z-direction in the center of the disc surrounded by a magnetic vortex (see figure 1).
For quantitative validation, the comparison between the results obtained from OpenPhase and OOMMF is shown in figure 2. We find an excellent agreement of the magnetization profiles (see panels (a)-(c)). Only the change of magnetization on the surface of the sample is more gradual in our implementation (see figure 2). This is however the expected trend as the profile of the surface in the MPF method is smooth, but sharp in OOMMF. The negligible difference in the vortex profile is caused by the PBC we apply, which allows the sample to interact with itself through the boundary and therefore enhances the demagnetizing field, as shown in figure 2(d).
In summary, the magnetic solver in our implementation perfectly reproduces the trends found by OOMMF, and the small quantitative differences can be related to the slightly different simulation setup.

Bi-crystal system
After validating our solver for the magnetization field, we now turn to the grain boundary motion. In order to disentangle the grain growth driven by the curvature of the grain boundary and the impact of the magnetic field, we first study a flat grain boundary in the rectangular sample as illustrated in figure 3.
We focus on a 60 • twist grain boundary and set the easy axis of grain 1 to E 1 A = (0, √ 3/2, −0.5) and that of grain 2 to E 2 A (0, √ 3/2, 0.5). The grain boundary is set perpendicular to the x-axis and located in the middle of the sample, as shown in figure 3. We randomly initialize the magnetization vectors which spontaneously align with the local easy axis and form antiparallel domains. This state is metastable and neither domains nor grains change further with time (see figure 4(a)).
Starting with this configuration we apply a homogeneous magnetic field (H ext ) along the easy axis of the second grain (E 2 A ). Figures 4(b)-(d) show the grain structure after 5000 steps for different field strengths. With the increase of H ext , the magnetization field becomes more uniform, and already for 800 kA m −1 , both grains are in a single domain state. Most important, grain 2 has grown at the expense of grain 1 for all field strengths. Thereby the traveling distance of the grain boundary increases strongly with the field strength for moderate field strength (1000 kA m −1 ), and it saturates beyond 5000 kA m −1 (cf figure 5(b)). These trends can be understood by the field-induced changes of the anisotropy energy density shown in figure 5(a). In the initial state (black) with H ext = 0, the magnetization vectors are either parallel or anti-parallel to their local easy axes to lower the anisotropic energy. Sharp peaks of the anisotropy energy in both grains are located at the domain walls.  In the applied field, the domain walls in grain 2 vanish, and as M is collinear to E A the energy density reaches zero for H > 200 kA m −1 . For grain 1, the magnetization gradually rotates with the field strength from its easy axis towards the external field direction, and thus the energy density increases systematically in the whole grain (cf colored lines). The energy of the system thus can be reduced if grain 2 grows at the expense of grain 1.
In the limiting case H → ∞, M is parallel to H ext . Then the maximum magnetic driving force in this sample can be estimated as: where θ H is the angle between the local easy axis and the external field. Equation (18) can also be obtained from the Stoner-Wohlfarth model.

Polycrystal film
In real polycrystalline or multi-phase systems, the interfaces and their junction points have a finite curvature which induces additional driving forces on the grain boundaries. Thus the . Both grains are separated by a flat wall normal to x and the external field is applied collinear to E 2 A . time-evolution of the grain structure results from the competition between this normal grain coarsening process and the magnetic driving forces discussed in the previous section.
Our implementation allows us to study this complex interplay for arbitrary crystal shapes. As an example, we study a polycrystal thin film of Sm 2 Co 17 surrounded by a vacuum with 151 × 151 × 15 grid points and 25 grains with identical phases but different orientations of the easy axes (see figure 6(a)). We compare three cases with different distributions of the easy axes. In all these cases, we initialize a distribution of the orientations around the y-axis and vary the width of the distribution. Without magnetic driving force, we observe the expected normal coarsening governed by interface energy and curvature. Grains with larger initial volume, e.g. grains 6 and 12, grow at the expense of their smaller neighbors, e.g. grains 1 and 17 (cf the changes from figures 6(a) and (b)).   Next, we apply a homogeneous magnetic field of 800 kA m −1 along the y-axis. Figures 6(c) and (d) illustrate the grain structure after 5000 steps in the external field for cases 1 and 3. In case 1, where deviation angles are less than 15 • resulting in an upper boundary of the magnetic driving force of 41 kJ m −3 (equation (18)) which is only moderate. Indeed also the grain structure after coarsening is similar to the nonmagnetic case (cf figures 6(b) and (c)). The main difference is the larger shrinkage of grain 2 in the applied field which has the largest deviation angle of all grains. For case 3, large changes are induced by the magnetic field. For example, grain 12 grows in normal coarsening but due to its deviation angle of 37 • it is fully eliminated by the external field. In addition, the magnetic field assists the growth of grains 4, 14, 18, and 21 with minimal deviation angles. Figure 7 compares initial and final volume fractions for case 1 (squares), case 2 (circles), and case 3 (crosses). The red line separates growing (above) and shrinking (below) grains. For comparison, the reference case (black pluses) shows the typical coarsening pattern without magnetic interactions. In the latter case, the growing grains have a larger initial volume fraction and are condensed in the right part of figure 7, while barely any grain with an initial volume fraction below 3.5% grows.
In the presence of the magnetic driving force, the probability for a grain to grow depends on volume and deviation angle. As shown by colored data points in figure 7, the deviation angles of most growing grains are <30 • , while those of the shrinking grains are usually >30 • . This trend is obvious for grains above an initial volume fraction of 3.5% for all three cases. A growing grain (grain 18 in case 3) below 3.5% points out that the growth for small grains is also possible if their magnetic energy is low enough to offset their disadvantage on interfacial energy. In summary, the external field increases the volume fraction of grains with a lower deviation angle. Therefore one may expect that the overall magnetocrystalline anisotropy energy of the sample and its magnetic performance are improved. Our implementation directly allows us to test this hypothesis. For this purpose, we randomly initialize the magnetization field for case 3 and determine the M − H hysteresis loop for a field along the y-direction. Figure 8(a) compares the hysteresis loops for the initial polycrystalline state (green crosses), after coarsening without magnetic field (red pluses), and after coarsening in the magnetic field (blue dots) to a single crystal with an easy axis along y (black squares). As expected, the single crystal shows a square-shaped hysteresis loop with the intrinsic remanence magnetization M ri = 782 kA m −1 and the intrinsic coercive field H ci = 572 kA m −1 .
For the polycrystalline material, the different orientations of the grain reduce M ri and H ci to 703 kA m −1 and 331 kA m −1 , respectively. The growth of disoriented grains may further weaken the coercive field. Thus after normal coarsening without magnetic driving force, some larger disoriented grains have grown and the sample properties have dropped to M ri = 678 kA m −1 and H ci = 319 kA m −1 . In full agreement with experimental findings in steel, the decrease of magnetic properties has also been related to the growth of grains with random orientation during coarsening [3]. In contrast to this, field-assisted coarsening indeed improves the magnetic performance of the sample to M ri = 745 kA m −1 and H ci = 385 kA m −1 .

Summary and conclusion
We implemented a micromagnetic model in the multi-phase field (MPF) open source code OpenPhase. To validate our implementation we simulated the magnetization in a soft FM disc and quantitatively reproduced results found by the micromagnetic code OOMMF. Importantly, our implementation allows us to assess the coupling between the micromagnetism and the grain structure evolution with time on the mesoscopic scale. Picking the example of Sm-based permanent magnets, we showed that uniaxial magnetic anisotropy in a single-phase polycrystal favors the growth of grains with a low deviation angle which allows for selective grain growth during the coarsening process and can improve the magnetic performance.
Note that our model can also be used for multi-phase materials, FM grains with other anisotropy classes, as well as mixtures of FM and para-or dia-magnetic grains. Furthermore, an extension of this model to grain boundary phases [44] is straightforward.

Data availability statement
The data cannot be made publicly available upon publication because they are not available in a format that is sufficiently accessible or reusable by other researchers. The data that support the findings of this study are available upon reasonable request from the authors.