Grain coarsening in two-dimensional phase-field models with an orientation field

TL;DR

This study investigates the limiting grain size distribution in two-dimensional phase-field models with an orientation field, revealing that resolution of small angle grain boundaries critically influences the observed distribution, aligning some models with experimental results.

Contribution

It demonstrates that insufficient resolution of small angle grain boundaries in orientation field models leads to lognormal distributions similar to experimental findings.

Findings

LGSD is sensitive to the resolution of small angle grain boundaries.

Insufficient resolution can cause models to produce lognormal distributions.

Differences in LGSD across studies may stem from boundary detection methods.

Abstract

In the literature, contradictory results have been published regarding the form of the limiting (long-time) grain size distribution (LGSD) that characterizes the late stage grain coarsening in two-dimensional and quasi-two-dimensional polycrystalline systems. While experiments and the phase-field crystal (PFC) model (a simple dynamical density functional theory) indicate a lognormal distribution, other works including theoretical studies based on conventional phase-field simulations that rely on coarse grained fields, like the multi-phase-field (MPF) and orientation field (OF) models, yield significantly different distributions. In a recent work, we have shown that the coarse grained phase-field models (whether MPF or OF) yield very similar limiting size distributions that seem to differ from the theoretical predictions. Herein, we revisit this problem, and demonstrate in the case of OF models [by R. Kobayashi et al., Physica D 140, 141 (2000) and H. Henry et al. Phys. Rev. B 86, 054117 (2012)] that an insufficient resolution of the small angle grain boundaries leads to a lognormal distribution close to those seen in the experiments and the molecular scale PFC simulations. Our work indicates, furthermore, that the LGSD is critically sensitive to the details of the evaluation process, and raises the possibility that the differences among the LGSD results from different sources may originate from differences in the detection of small angle grain boundaries.

Figures

Equations

$$\displaystyle p(x)=\frac{1}{x\sigma\sqrt{2\pi}}\mathrm{exp}\bigg{\{}-\frac{[ln(x)-\mu]^{2}}{2\sigma^{2}}\bigg{\}},$$

$$\displaystyle p(x)=(2e)^{2}\frac{2x}{(2-x)^{4}}\mathrm{exp}\bigg{\{}-\frac{4}{2-x}\bigg{\}}.$$

$$\displaystyle p(x)=\frac{\beta}{\alpha}\left(\frac{x}{\alpha}\right)^{(\beta-1)}\mathrm{exp}\bigg{\{}-\left(\frac{x}{\alpha}\right)^{\beta}\bigg{\}},$$

$$\displaystyle F=\int d\mathbf{r}\bigg{\{}\frac{\epsilon^{2}T}{2}\left(\nabla\phi\right)^{2}$$

$$\displaystyle+f_{\mathrm{bulk}}(\phi)+f_{\mathrm{ori}}(\phi,\nabla\theta)\bigg{\}},$$

$$f_{\mathrm{bulk}}(\varphi )$$

$$p(\varphi )$$

$$f_{\mathrm{ori}}^{\mathrm{KWC}}(\varphi ,\nabla \theta )$$

$$f_{\mathrm{ori}}^{\mathrm{HMP}}(\varphi ,\nabla \theta )$$

$$\frac{\partial \varphi }{\partial t}$$

$$\frac{\partial \theta }{\partial t}$$

$$\displaystyle F=\int d\mathbf{r}\bigg{\{}\frac{\epsilon^{2}T}{2}\left(\nabla\phi\right)^{2}+WTg(\phi)+[1-p(\phi)]\Delta f+$$

$$\displaystyle+\phi^{4}\left[H_{1}T|\nabla\theta|+H_{2}T(\nabla\theta)^{2}\right]\bigg{\}},$$

$$\displaystyle\dot{\phi}=M_{\phi}\bigg{\{}\epsilon^{2}T(\nabla^{2}\phi)-WTg^{\prime}(\phi)+p^{\prime}(\phi)\Delta f-$$

$$\displaystyle-4\phi^{3}\left(H_{1}T|\nabla\theta|+H_{2}T(\nabla\theta)^{2}\right)\bigg{\}}$$

$$\displaystyle\dot{\tilde{\phi}}=\tilde{M}_{\phi}\bigg{\{}\tilde{\nabla}^{2}\phi-\tilde{W}g^{\prime}(\phi)+\tilde{\Delta f}p^{\prime}(\phi)-$$

$$\displaystyle-4\phi^{3}\left(\tilde{H}_{1}|\tilde{\nabla}\theta|+\tilde{H}_{2}(\tilde{\nabla}\theta)^{2}\right)\bigg{\}}$$

$$\dot{\theta }=M_{\theta }\nabla \left\{{\varphi }^4\left[H_{1}T\frac{\nabla \theta }{|\nabla \theta |}+2H_{2}T(\nabla \theta )\right]\right\}$$

$$\tilde{\dot{\theta }}={\tilde{M}}_{\theta }\tilde{\nabla }\left\{{\varphi }^4\left[\frac{\tilde{\nabla }\theta }{|\tilde{\nabla }\theta |}+2\widetilde{H_2^{\theta }}\tilde{\nabla }\theta \right]\right\}$$

$$\displaystyle F=\int d\mathbf{r}\bigg{\{}\frac{\epsilon^{2}T}{2}\left(\nabla\phi\right)^{2}+WTg(\phi)-p(\phi)\Delta f+$$

$$\displaystyle+q(\phi)HT(\nabla\theta)^{2}\bigg{\}},$$

$$\displaystyle\dot{\phi}=M_{\phi}\big{\{}\epsilon^{2}T(\nabla^{2}\phi)-WTg^{\prime}(\phi)+p^{\prime}(\phi)\Delta f-$$

$$\displaystyle-q^{\prime}(\phi)HT(\nabla\theta)^{2}\big{\}},$$

$$\displaystyle\tilde{\dot{\phi}}=\tilde{M}_{\phi}\big{\{}\tilde{\nabla}^{2}\phi-\tilde{W}g^{\prime}(\phi)+\tilde{\Delta f}p^{\prime}(\phi)-$$

$$\displaystyle-q^{\prime}(\phi)\tilde{H}(\tilde{\nabla}\theta)^{2}\big{\}},$$

$$\dot{\theta }=M_{\theta }\nabla \left\{HTq(\varphi )(\nabla \theta )\right\},$$

$$\tilde{\dot{\theta }}={\tilde{M}}_{\theta }\tilde{\nabla }\left\{q(\varphi )(\tilde{\nabla }\theta )\right\},$$

Full text

I Introduction

The majority of the solid matter we use appears in a polycrystalline form (including technical alloys, concrete, polymers, minerals, drugs, sugar, salt, cholesterol, peptides, etc.); i.e., they are composed of a large number of small crystallites. The properties of such materials depend on the size-, composition-, and shape-distributions of the crystallites they consist of. This also appears to be the case for polycrystalline thin metal films, where the properties of the grain size distribution may, e.g., influence the quality of metallization of semiconductors in electronics industry (by influencing the current-carrying capability via the resistance to electromigration ref1 ; ref2 ). The distributions that characterize the polycrystalline structure can be influenced by the conditions of preparation (e.g., solidification, electrodeposition or sputtering) and subsequent processing, including heat treatments. As a result, the understanding of the grain coarsening process and the ability to predict the associated properties of the polycrystalline matter are of high technological importance, and have been the subject of intensive experimental and theoretical research.

Grain coarsening is mostly due to the motion of grain boundaries and triple junctions in a way that leads to the reduction of the excess free energy associated with the grain boundary network. While this process is fairly simple, the complexity of the geometry makes it difficult to predict the limiting grain size distribution (LGSD). Theory, experiments, and simulations agree that there exists a (time invariant) limiting grain size distribution, which evolves in a self-similar way, and the time dependence of the average grain size can be expressed as $\langle R\rangle=kt^{n}$, where $n$ is the growth exponent ref1 ; ref3 . Despite long-standing efforts, no convincing theoretical derivation of the LGSD has been proposed yet. Moreover, other phenomena such as grain rotation, elasticity, and anisotropy may affect the LGSD, adding to the difficulty of the task.

Herein, we concentrate on grain coarsening in two dimensional (2D) systems. The present results are expected to be relevant to 2D multi-grain structures, including thin metallic/ceramic films, colloidal aggregates, and plasma crystals. In the following paragraphs we briefly review the 2D results from experiment, theory, and simulation. This is a huge area indeed, and we can only summarize the main results here. Although the relevant process have been studied for some time, this brief review will show that the results are not without contradiction and that no clear understanding of the discrepancies among them has been reached yet.

A recent detailed experimental study that summarizes results on 27 sputter deposited Al and Cu films ref1 indicates a lognormal LGSD with $\sigma=0.5$ and $\mu=-0.12$ (Fig. 1), and $n\approx 1/2$ for the growth exponent. The lognormal distribution,

[TABLE]

where $x$ is the linear size of the grains evaluated from their area, was found to be fairly robust. It remains an accurate description under a broad range of experimental conditions, including purity of the sputtering target, type of substrate, film thickness, temperature of deposition, actual and homologous annealing temperatures, the time of annealing, the grain size, and the twin density within the grains ref1 . This result is in agreement with some previous experimental results on thin films ref4 ; ref5 , while other works indicate deviations from the lognormal distribution ref6 . It is also worth noting that significantly smaller $n$ values have also been reported in other works ref7 ; ref8 . The main cause for grain coarsening is the excess free energy of the grain boundary network, which relaxes via curvature driven migration of the grain boundaries ref9 . The velocity of the latter is often sought in the form $v=\kappa M$, where $\kappa$ is the curvature, and $M$ the grain boundary mobility. It is expected that other phenomena, such as grain rotation, elasticity, anisotropy, etc. may also contribute to the dynamics of grain coarsening ref10 .

There appears to be no consensus regarding the theoretical form of LGSD in 2D: Kolmogorov ref11 has shown that starting from a single grain of known volume, repeated crushing will lead to a lognormal distribution of particle volumes. In Feltham’s early work ref4 a lognormal distribution was postulated. In a more sophisticated theoretical treatment, that follows the route of Lifshitz, Slyozov, and Wagner used in addressing Ostwald ripening ref12 ; ref13 , Hillert ref14 derived the distribution function for two dimensional isotropic growth in multigrain structures

[TABLE]

Regarding grain growth as a statistical phenomenon taking place via the growth of faces, Louat ref15 derived a Rayleigh type grain size distribution. A statistical extension of Feltham’s approach by Kurtz and Carpay, in turn, supports the lognormal form ref16 . The empirical Weibull distribution, corresponding to

[TABLE]

and its special case the Rayleigh distribution $(\beta=2)$ were also used to approximate the experimentally observed LGSDs ref17 . While the use of Weibull distribution is empirical in this context, the Rayleigh distribution was derived by Louat ref15 . Relying on the Neumann-Mullins growth law, in more recent works Pande and Cooper ref18 ; ref19 deduce a Fokker-Planck equation for the grain size distribution, which yields a self-similar asymptotic solution that can be reached from arbitrary initial state. They propose an approximate analytical solution for LGSD, that can be tuned by varying a single parameter between the Rayleigh distribution (where all curvature effects are neglected) and Hillert’s model (where the drift velocity due to curvature is the only driving force) ref19 . Seeking a flexible empirical LGSD distribution, Rickman et al. quantify the deviation from the lognormal distribution in terms of a cumulant expansion tailored to the lognormal distribution ref20 .

The models discussed above predict $n=1/2$. A comparison of the predicted LGSDs with parameters taken from fitting to recent experimental data of ref1 are presented in Fig. 1. In agreement with previous work ref1 ; ref5 ; ref19 , these distributions appear to be less satisfactory than the lognormal form.

Computer simulations have also been used extensively to investigate grain coarsening. The results appear again somewhat contradictory. For example early results on Monte Carlo simulations for the Q = 32 2D Potts model ref21 indicate a fair agreement with Hillert’s mean field model, yet a later Q = 72 Potts study reports a substantially different LGSD ref22 . The other 2D models predict rather similar limiting grain size distributions. This includes a grain boundary migration based numerical approach by Moldovan et al. ref10 , three different versions of the multi-phase-field theory (MPF1 ref22 , MPF2 ref23 , MPF3 ref24 ), and two orientation field based models: one by Kobayashi, Warren, and Carter (KWC) ref25 and another by Henry, Mellenthin, and Plapp (HMP) ref26 . LGSDs for the OF models were reported in Ref. ref27 . Another approach that led to similar results employed a numerical surface solver to relax the interface energy ref28 . The predicted distributions are fairly similar (with some scattering), although different models and methods are compared. This finding might originate from the fact that these models all tend to reduce the free energy associated with the grain boundary network. Apparently, LGSD from these 2D simulations fit reasonably well to the Weibull type probability density. Although the investigated 2D simulations are consistent with each other (Fig. 2), and follow the behavior obtained using the model of Mullins ref9 in Ref. ref1 , yet they significantly differ from the lognormal distribution representing recent experimental data ref1 . These findings are consistent with earlier results of Pande and Cooper ref18 ; ref19 . Interestingly, simulations performed using a dramatically different process for coarsening, i.e. grain rotation, yield a lognormal distribution ref10 , yet much different from experiments (see Fig. 3).

Remarkably, the only simulation results for LGSD that agree with the experiments are from the phase-field crystal (PFC) model, an approach that works on the molecular scale ref29 (see Fig. 3). The PFC model developed by Elder et al. ref30 can be regarded as a simple dynamical density functional theory that incorporates the crystal structure, anisotropies, elasticity, and dislocations automatically (for a review on PFC see Ref. ref31 ). As such, it contains a much richer physics that the previous models, including automatically elasticity, dislocation dynamics, grain rotation, molecular scale description of grain boundaries, etc. A possibility is that this richness of phenomena is responsible for the accurate LGSD it predicts. Herein, we raise a different possibility, which might explain at least partly the observed differences between LGSDs from simulations and experiments.

We use the orientation field models HMP and KWC to demonstrate that the limiting grain size distribution is critically sensitive to details of the evaluation of the number of the grains, especially to the resolution of the small angle grain boundaries. Variation of the misorientation, below which a grain boundary is not detected any more, yields a continuous transition between the lognormal distribution that can be observed for poorly resolved low angle grain boundaries (only a fraction of them are found), and the general behavior from 2D simulations shown in Fig. 2, obtained when the low angle grain boundaries are well resolved (the majority of them was found). This finding raises the possibility that improving the resolution of low angle grain boundaries in the experiments (and maybe in PFC simulations) might yield LGSDs falling closer to the results provided by the 2D simulations.

The present paper is structured as follows: In Section II, we briefly recapitulate the essence of the orientation field theories HMP and KWC, while Section III specifies the materials properties and other conditions used in the simulations. Section IV describes the details of the evaluation methods employed for determining the LGSD, with a detailed analysis of different factors that influence the limiting distribution, and discuss the consequences. Finally, in Section V, we summarize the main results and offer a few concluding remarks.

II Orientation field models

We employ here two phase-field models that describe polycrystalline solidification and grain coarsening on equal footing: the KWC and HMP models described in Refs. ref25 and ref26 , respectively. In these models the transition between the liquid and crystalline phases is monitored by a structural order parameter, the phase field $\phi(\mathbf{r},t)$, whereas the local crystallographic orientation is specified by a scalar orientation field $\theta(\mathbf{r},t)$, normalized so that $\theta\in[0,1]$, considering the crystal symmetry. For $k$-fold symmetry, the orientation angle $\vartheta$ can vary between 0 and $2\pi/k$. Owing to the $k$-fold symmetry, $\vartheta$ angles outside this region are equivalent to a specific orientation angle inside the $[0,2\pi/k]$ regime. (Then $\theta=\vartheta/(2\pi/k)\in[0,1]$ describes all orientations.) Being an angular variable $\theta=0$ and $\theta=1$ are equivalent, and the magnitude of the orientation field difference is limited: $|\Delta\theta|\leq\frac{1}{2}$. This should be considered, e.g., when evaluating the differential operators acting on the orientation field. For the sake of simplicity, isotropic systems will be considered, yet the grain boundary energy depends naturally on the misorientation in these models, as will be displayed in Fig. 4(a) below.

II.1 Free energy functional

The free energy of the respective system is a functional of the fields $\phi(\mathbf{r},t)$ and $\theta(\mathbf{r},t)$, which can be written in the following form for both orientation field approaches we use here:

[TABLE]

where the model parameters $\epsilon^{2}$ and $W$ are expressible in terms of the free energy and thickness of equilibrium solid-liquid interface, and $T$ is the temperature. $f_{\mathrm{bulk}}$ switches between the free energy densities of the bulk solid and liquid phases according to the interpolating function $p(\phi)$ as follows

[TABLE]

where the free energy densities of the bulk solid and liquid phases, $f_{\mathrm{s}}(T)$ and $f_{\mathrm{l}}(T)$, respectively, can be taken from thermodynamic databases.

The two orientation field based models we use differ in the form of the orientational contribution to the free energy, $f_{\mathrm{ori}}(\phi,\nabla\theta)$:

[TABLE]

where $r(\phi)=\phi^{4}$ and $q(\phi)=\frac{7\phi^{3}-6\phi^{4}}{(1-\phi)^{3}}$. These functions differ slightly form the original formulations ref32 . Parameters $H$, $H_{1}$, and $H_{2}$ tune the strength of the respective gradient terms, and scale the grain boundary energy.

II.2 Equations of motion (EOMs)

The time evolution of the system is assumed to follow standard variational dynamics of non-conserved fields:

[TABLE]

[TABLE]

where $M_{\phi}(=const.)$ and $M_{\theta}$ are the respective mobilities, $M_{\theta}=M_{\theta,S}+(M_{\theta,L}+M_{\theta,S})p(\phi)$, whereas $\zeta_{\phi}$ and $\zeta_{\theta}$ are noise terms that represent the fluctuations of the respective fields. Their correlators are $\langle\zeta_{i}(\mathbf{r},t),\zeta_{i}(\mathbf{r}^{\prime},t^{\prime})\rangle=\omega_{ij}2kTM_{ij}\delta(\mathbf{r}-\mathbf{r}^{\prime})\delta(t-t^{\prime})$, where $\omega_{ij}$ is the noise strength coefficient for the $i$th field in the $j$th phase ($i=\phi,\theta$ and $j=L,S$). For further details of EOMs and their dimensionless forms see Appendix I.

We found that in the orientation field models topological defects appear that resemble the disclinations observed in 2D atomistic models using the hexatic order parameter ref33 , which may influence the motion of grain boundaries and trijunctions. This phenomenon has already been detected by Warren et al. ref25 (c). The properties of these defects and possible ways to remove them via a complex orientation field, a 3-component orientation field, and other means will be addressed elsewhere ref34 . Herein, we use $\omega_{\phi,S,L}=0,$ $\omega_{\theta,S}=0$ and $\omega_{\theta,L}=0.1$, to remove the pinning effect of these topological defects in the orientation field.

II.3 Numerics

The dimensionless form of Eqs. (9) and (10) were solved numerically on rectangular grids of different sizes, using finite difference discretization combined with explicit forward Euler time stepping, while prescribing periodic boundary conditions. Parallel codes were developed for a CPU cluster and GPU cards. The computations were performed on GPU cards of various types.

III Materials properties and other conditions

In the simulations, we used the physical properties of pure Ni. The volumetric free energy difference between the liquid and the solid was estimated using Turnbull’s linear approximation: $\Delta f=\Delta H_{f}(T_{f}-T)/T_{f}$ ref35 , where $\Delta H_{f}=2.61\times 10^{9}$ J/m3 and $T_{f}=1728$ K are the volumetric heat of fusion and the melting point, respectively. The thickness of the equilibrium solid-liquid interface were taken as $d=2$ nm, which is of the order of magnitude of results from molecular dynamics simulations ref36 . The free energy of the solid-liquid interface, $\gamma_{SL}=0.364$ J/m2 was taken from the compilation ref37 ; whereas a molar volume $V_{m}=6.59\times 10^{-6}$ m3/mol was employed. Isothermal computations were performed at $T=974$ K. The diffusion coefficient and the characteristic length applied in making the EOMs dimensionless (see Appendix I) were $D_{L}=10^{-9}$ m2/s and $\xi=40$ nm, respectively. If not stated otherwise, the dimensionless spatial and time steps were chosen as $\Delta\tilde{x}=3.125\times 10^{-3}$ and $\Delta\tilde{t}=1.4844\times 10^{-8}$, which ensured the numerical stability of the solution. The dimensionless phase-field mobility was taken as $\tilde{M}_{\phi}^{HMP}=\tilde{M}_{\phi}^{KWC}=0.9$.

III.1 Synchronizing the HMP and KWC models

Herein, we synchronize the HMP and KWC models so that they display similar properties (grain-boundary energy and time scale for grain-rotation) under the same circumstances; i.e., they can be viewed as models describing polycrystalline solidification/grain coarsening in the same matter.

To match the properties of the solid-solid interface, we have first chosen the dimensionless parameters (see Appendix) $\alpha_{2}=1.0$, $\tilde{M}_{\theta,S}^{HMP}=28.8$, and $\tilde{M}_{\theta,L}^{HMP}=92.68$, and then varied the parameters of the KWC model until similar grain boundary energies, interface profiles, and time scales were obtained. This is realized by the following set of dimensionless parameters $\alpha_{1}=10.64$, $\alpha_{2}=30.0$, $\tilde{M}_{\theta,S}^{KWC}=7.2\times 10^{-4}$, and $\tilde{M}_{\theta,L}^{KWC}=1080$. These data ensure that in the bulk crystal ($\phi=1$) the orientation mobility is negligible in both cases.

The respective grain boundary energies and interfacial profiles are shown in Fig. 4. The two models provide indeed fairly coherent predictions. One of the important features is the depression of the phase field at the grain boundary, which indicates that the crystalline order is disturbed within the grain boundary. For high misorientations the phase field may get as low as $\phi\sim 0.3$ here, but may reach $\phi=0$ at the melting point that implies the presence of a highly disordered structure at the grain boundary, a behavior similar to the one observed in molecular dynamics simulations ref38 . The relationships between the depth of the phase-field depression and the misorientation are shown for the two models in Fig. 5.

Next, we compare grain growth in the models. With the present choice of mobilities, the growth rates of the two models are of comparable magnitude [see Fig. 6]. The simulation data were fitted by the formula $\langle R\rangle=A(t-t_{0})^{n}$, where $t_{0}$ is the hypothetical starting point of growth corresponding to $\langle R\rangle=0$. The fitting was performed so that only data beyond the end of an apparent initial transient period ($8\times 10^{5}$ and $10^{6}$ time steps for the HMP and KWC models, respectively) were considered. This procedure yielded growth exponents comparable to those from other models (HMP: $n=0.45\pm 0.01$; KWC: $n=0.51\pm 0.01$).

III.2 Evaluation of LGSD

In determining the LGSD, we used typically eight $4096^{2}$ simulations performed with different initializations of the random number generator. We used one of the following methods to generate the initial grain size distribution: (i) added noise to the equation of motion of the phase field to initiate nucleation and growth; (ii) we placed randomly oriented particles of 30 pixel radius randomly in the simulation box and let them grow; (iii) we also explored the case when the randomly oriented particles were placed on a square grid. The respective simulations yielded similar LGSDs, indicating that the long time behavior is not sensitive to the initial conditions.

Unless stated otherwise, the grains were identified by the watershed algorithm of MATLAB ref39 . The watershed algorithm finds ”catchment basins” and ”watershed ridge lines” in an image by treating it as a surface, where light pixels represent high elevations and dark pixels represent low elevations. It has been applied to evaluate the grain size distribution from the $(1-\phi)$ maps. In this case, the catchment basins correspond to the grains and the watershed ridge lines are the grain boundaries.

We have found that even merging eight $4096^{2}$ simulations (that contained initially $\sim 17,200$ randomly oriented grains) the grain size histograms show some visible scattering that makes difficult the comparison of LGSDs not very far from each other. It was, however, detected that at late stages of the evolution of the grain boundary network (in which we are interested anyway), fitting an appropriate analytic formula to the histograms (lognormal or Weibull), the results are essentially indistinguishable (i.e., they match to several states of the grain boundary network with a similar accuracy). In the case of such states of the system, it makes sense to merge the reduced distributions, to reduce the scattering of the histograms. One needs to be careful, however, since this is possible only at the end of the relaxation process, where the limiting size distribution is established. We have employed this technique for the evaluation of all LGSDs shown. For the details see Appendix II.

We have also evaluated the probability distribution of the misorientations along the grain boundaries identified by the watershed algorithm. For the details see Appendix III.

IV Results and discussion

The aim of the present investigation is to clarify whether differences in the accuracy of the evaluations can significantly influence the results for the LGSD. This was motivated by the fact that seemingly small changes in our evaluation procedures for the LGSD could yield significantly different distributions starting from the same raw data. Here, we pinpoint the reasons for this behavior and show that it is generic for orientation-field models. More precisely, the shape of the LGSD critically depends on the detection of small-angle grain boundaries. Beyond the models treated here, this observation raises several questions that will be discussed below.

IV.1 Methods for distinguishing the grains

There are different possibilities to evaluate the grain size distribution from the phase- and orientation fields. Since the orientation field varies continuously across the grain boundaries and thus the same orientations may occur in a grain and at grain boundaries elsewhere, the recognition of grains is more complicated from the orientation field. Thus, in a previous work ref27 , we have opted for the identification of the grains via the depression of the phase field at the grain boundary, and used the watershed algorithm to locate them, yielding the red curves in Fig. 2. Herein we introduce a different approach that is also based on the phase-field map, yet relies on a threshold value in recognizing where the grain boundaries lie, allowing thus for a continuous tuning of the fraction of small angle grain boundaries considered in the evaluation.

IV.2 Tunable approach to distinguish grains

As demonstrated by Figs. 4 and 5, the depth and width of the phase-field depression at the grain boundary depends on the misorientation. As a result, recognition of the small angle grain boundaries becomes more difficult, as they become less visible for $\Delta\theta\rightarrow 0$. The evaluation procedure is illustrated in Fig. 7. Figs. 7(a) and 7(b) show the original orientation field and phase-field distributions. Here, we use a discretized gray-scale representation of the phase field that has 256 shades. We use the watershed algorithm of MATLAB ref39 to determine the grain boundary network. This discretization of the gray hue removes a minor ($\sim 10^{-4}$) scattering of the phase field in the vicinity of the grain boundaries that originate from the fluctuations of the orientation field in the liquid phase. (This scattering would otherwise produce very small grain sizes, when applying the watershed algorithm.) This procedure results in well discernible grain boundaries even in the case of small angle grain boundaries, where the phase-field depression is small [see Fig. 7(b)]. One expects that in the experiments recognition of small angle grain boundaries is more difficult. A similar situation occurs here, as the phase-field depression varies with misorientation. In order to understand how the LGSD depends on the recognition of the small angle grain boundaries, we have processed the gray-scale image further by converting it to black and white using a threshold value for phase field, $\phi_{th}$. If $\phi\leq\phi_{th}$ the actual pixel is considered as belonging to a grain boundary, and is painted black, the rest is colored white. Increasing $\phi_{th}$, an increasing number of pixels are recognized as belonging to grain boundaries, and an increasing fraction of the shallower depressions representing the small angle grain boundaries are detected [see Figs. 7(c)–7(f)]. Evidently, with increasing $\phi_{th}$ an increasing fraction of the simulation box becomes black. The theoretical upper limit to detect any difference using this method is $\phi_{th}>1-1/256$. Once the grain boundary network was computed, the grains that are defined as white areas, whose pixels can all be visited from its other pixels without crossing a grain boundary, using the watershed algorithm. The grains identified so are then divided into size categories, and presented in the form of histograms approximating the respective probability density distribution.

IV.2.1 Varying the threshold

As one may expect on the basis of Fig. 5, the LGSD depends on the choice of the threshold $\phi_{th}$. The LGSDs obtained using the threshold values applied in Figs. 7(c)–7(f) are compared with each other and the experimental distribution in Fig. 8. We have also evaluated the probability distribution of the misorientations at the grain boundaries (see Fig. 9). A comparison of the distributions corresponding to $\phi_{th}=0.79$ and $\phi_{th}=0.99$ clearly indicate that the major difference between the respective misorientation distribution is indeed a lack of small angle misorientations in the former case. It appears that with increasing $\phi_{th}$, we see a transition from a distribution falling close to the experimental lognormal distribution towards the generic Mullins-type LGSD the majority of the 2D simulations predict (cf. LGSD colored blue and red in Fig. 2). The distributions obtained using $\phi_{th}=0.99$ are in fact very close to those obtained from the grayscale image of the phase field map using the watershed algorithm (see Fig. 10), implying that with this threshold the majority of the low angle grain boundaries was found.

These results indicate that the LGSD is critically sensitive to the resolution of the small angle grain boundaries, and raise the possibility that this sensitivity may at least partly be responsible for the deviation between the results of the (fairly coherent) 2D simulations shown in Fig. 2 and other methods that lead to different LGSDs. We note in this respect that with the exception of approaches based on the OF concept; the simulations shown in Fig. 2 are immune to such errors: e.g., the multi-phase-field models use an individual phase-field for each orientation, so in these models grains can always be distinguished.

IV.2.2 Discrete orientations

The MPF simulations are performed usually with a relatively large but finite number of discrete orientations (in early studies about 30 equidistant orientations were regarded as a satisfactory approximation in 2D ref40 , however, recently $\sim 100$ orientations are considered more appropriate ref23 ). We have tried a similar approach in the orientation field models: we started simulations with initial grains of 30 equidistant orientations of equal probability. Apart from the effect of establishing a continuous transition at the grain boundaries, the initial orientations remained dominant during the process of grain coarsening. This results in grain boundaries with misorientations that are an integer multiples of the minimal misorientation $\Delta\theta_{\mathrm{min}}=1/30$. Accordingly, when varying $\phi_{th}$, the change of the LGSD happens stepwise between $\phi_{th}=0.85$ and $0.86$ for the KWC model, and between $\phi_{th}=0.81$ and $0.82$ in the case of the HMP model. This stepwise behavior is also visible, when plotting the number of detected grains as a function of the misorientation $\Delta\theta$ that corresponds to the threshold $\phi_{th}$ used in the evaluation process (Fig. 12). For grain orientations varying continuously, the number of grains also vary continuously, whereas for equidistant discrete orientations the number of grains varies stepwise. The respective misorientation distributions reflect these [see Fig. 13(a)]. Furthermore, in the case of the lower threshold ($\phi_{th}=0.85$), the first peak and part of the second is missing from the misorientation distribution, which explains the abrupt change of the LGSD, and shows again the importance of the amount of small angle grain boundaries in shaping the LGSD [see Fig. 13(b)].

IV.2.3 A few additional remarks

It is worth mentioning that although in the OF models the grain boundary energy is misorientation dependent, still it does not depend on the inclination of the interface, whereas the dependence on inclination is evidently present in the experiments and the PFC simulations. Also, all the simulations considered herein are strictly two-dimensional. Yet, in principle, it cannot be excluded that as the real experiments are quasi-two-dimensional (the thickness is not fully negligible), they have to be modeled as thin but three-dimensional cases, when addressing the LGSD. Such studies are, however, out of the scope of the present work.

We note furthermore that the present results may have implications regarding the controversy between the lognormal LGSD from experiment and PFC simulations, and the distribution that appears to be a generic (Mullins-type) solution from other 2D simulations: (a) It might be possible that the population of the small angle grain boundaries is indeed lower in the experiments and the PFC model than in the simulations shown in 2. If so, evaluation of the misorientation distribution may be useful tool for a more complete characterization of the grain coarsening process. (b) Another theoretical possibility is that the resolution of the small angle grain boundaries is less satisfactory in the cases of experiments and PFC simulations. While the TEM experiments (whose results are shown as reference in Figs. 1, 2, and 3) are expected to resolve all grain boundaries no matter what the misorientation is ref41 , other methods may be less successful in detecting the small angle grain boundaries. Concerning the PFC simulations, our experience in analyzing polycrystalline patterns in terms of the complex hexatic order parameter seems to indicate that the evaluation of the numbers of grains may be subject to uncertainties especially when the small angle grain boundaries are defined by only a few dislocations ref42 . Work is underway to clarify these issues further.

Finally, we wish to stress the point that grain coarsening is a complex process, which depends on several factors. Our work shows that minor changes in the detection of low angle grain boundaries can lead to completely different LGSDs, which may hide the effect of other factors. Although our study is limited to the OF models, it indicates that extreme care needs to be taken to ensure the accuracy of the grain size distribution.

V Summary

Using orientation field based phase-field models we have investigated how the detection of small angle grain boundaries influence the limiting grain size distribution (LGSD), towards which the grain boundary network relaxes at long times. We make the following concluding remarks:

(i) It appears that the 2D computer simulations relying on the orientation field based phase-field models predict LGSDs that are consistent with LGSDs from the Q = 72 state Potts model, Mullins’ model, three versions of the multi-phase-field theory, and with results from a numerical surface solver. These distributions, however, differ significantly from the lognormal LGSD emerging from the experiments and the 2D phase-field crystal model.

(ii) In the orientation field models we have observed that the LGSD is critically sensitive to the detection of the small angle grain boundaries: We introduced an evaluation method in which the variation of a threshold changed the fraction of small angle grain boundaries detected. It has been shown then that considering an increasing fraction of the small angle grain boundaries, the LGSD varies from a lognormal distribution falling close to those emerging from the experiments and the phase-field-crystal model, to a Mullins-type LGSD the other 2D computer simulations predict. The respective changes in the amount of small angle grain boundaries are clearly seen in the misorientation distribution along the grain boundary network.

(iii) Further work is needed to clarify whether in the experiments and phase-field-crystal simulation, the population of the small angle grain boundaries is indeed smaller than in other computer simulations, or perhaps other effects are responsible for the observed deviations in the respective LGSDs.

Acknowledgements.

This work has been supported by NKFIH, Hungary under project No. K-115959, by the Hungarian-French Bilateral Scientific and Technological Innovation Fund under Grant No. TÉT$\_12\_$FR-2-2014-0034; by the ESA MAP/PECS projects “MAGNEPHAS III” (Contract No 40000110756/11/NL/KML) and “GRADECET” (Contract No 40000110759/11/NL/KML), and by the EU FP7 project “EXOMET” (Contract No. NMP-LA-2012-280421, co-funded by ESA).

Appendix I: Dimensionless form of the orientation field models

V.1 Kobayashi-Warren-Carter-model:

The free energy functional:

[TABLE]

where

[TABLE]

EOM for the phase-field:

[TABLE]

measuring the length and time in units $\xi$ and $\tau=\xi^{2}/D_{L}$, respectively, one obtains the following dimensionless equation of motion:

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

EOM for the orientation field:

[TABLE]

or in dimensionless form:

[TABLE]

where $\tilde{M}_{\theta}^{KWC}=M_{\theta}\xi H_{1}T/D_{L}$ and $\tilde{H}_{2}^{\theta}=\frac{\tilde{H}_{2}}{\tilde{H_{1}}}=\frac{\alpha_{2}}{\alpha_{1}}\frac{2\delta}{\xi}$.

V.2 The Henry-Mellenthin-Plapp model:

Free energy functional:

[TABLE]

where $H=\alpha_{2}\frac{8\gamma_{SL}\delta}{T_{m}}$, whereas $\epsilon^{2}$, $W$, and $\Delta f$ are the same as for the KWC model.

EOM for the phase-field:

[TABLE]

yielding the following dimensionless equation of motion:

[TABLE]

where $\tilde{H}=\frac{H}{\epsilon^{2}}=\alpha_{2}\frac{8\gamma_{SL}\delta}{T_{m}}\frac{T_{m}}{6\sqrt{2}\gamma_{SL}\delta}=\alpha_{2}\frac{2\sqrt{2}}{3}$.

EOM for the orientation field:

[TABLE]

and in dimensionless form:

[TABLE]

while $\tilde{M}_{\theta}^{HMP}=M_{\theta}HT/D_{L}$.

Appendix II: Reducing the scattering of LGSD

Under the conditions used herein, practically steady state distribution was achieved at dimensionless time of about $\tilde{t_{1}}=2\times 10^{6}\Delta\tilde{t}$ for both the HMP and KWC models. The data representing the state of the system were saved after every $10^{5}$ time steps until reaching dimensionless time of $\tilde{t_{2}}=4\times 10^{6}\Delta\tilde{t}$. The raw distributions corresponding to $\tilde{t_{1}}$ and $\tilde{t_{2}}$ are shown in Fig. 14. In this period the average size varied from $\langle R\rangle/\Delta x=61.0$ to $100.7$ for HMP and from $\langle R\rangle/\Delta x=57.1$ to $96.0$ for KWC. The reduced distributions and the respective fitted (lognormal) distributions are presented in Fig. 15(a) and 15(b) for the HMP and KWC models. The difference of the distributions fitted at the two limiting cases is characterized by the total variational difference $\delta_{T}=\frac{1}{2}\int_{0}^{\infty}|p_{t_{1}}(x)-p_{t_{2}}(x)|dx$, yielding $\delta_{T}=0.020$ for HMP and 0.011 for KWC, respectively. The LGSDs obtained by merging the respective 21 reduced distributions are reasonably smooth (see Fig. 16), and enables the detection of small differences between LGSDs that would be hardly perceptible otherwise due to statistical scattering. Here $\phi_{th}=0.79$ was used. Somewhat larger differences were observed in the case of Weibull fits to the data from the watershed algorithm ($\delta_{T}=0.020$ for HMP and 0.031 for KWC), as there the applied Weibull functions approximate the distributions less accurately.

Appendix III: Evaluation of the misorientation distribution

In order to characterize the grain boundary network, we have evaluated the distribution of the misorientations weighted with the number of pixels occurring in the watershed map (Fig. 17). Pixels assigned to the ”catchment basins” found by the watershed method were associated with the orientations of the respective areas. The number of the pixels in the ”watershed ridge lines” were used to represent the frequency of the local misorientation. A histogram was made of the latter, an approximation of the probability density distribution of misorientations along the grain boundary network. The misorientation distribution defined so converges towards a limiting distribution at about the same time as the LGSD.