# Grain growth competition and formation of grain boundaries during
solidification of hcp alloys

A.K. Boukellal1,†,∗ M. Sarebanzadeh1,2 A. Orozco-Caballero3 F. Sket1 J.
LLorca1,2 D. Tourret1,∗ 1 IMDEA Materials Institute, Getafe, 28906 Madrid,
Spain.  2 Department of Materials Science, Polytechnic University of
Madrid/Universidad Politécnica de Madrid, 28040 Madrid, Spain.  3 Department
of Mechanical Engineering, Chemistry and Industrial Design, Polytechnic
University of Madrid/Universidad Politécnica de Madrid, 28012 Madrid, Spain.

###### Abstract

Grain growth competition during directional solidification of a polycrystal
with hexagonal (hcp) symmetry (Mg-1wt$\%$Gd alloy) is studied by phase-field
modeling, exploring the effect of the temperature gradient $G$ on the
resulting grain boundary (GB) orientation selection. Results show that
selection mechanisms and scaling laws derived for cubic (fcc, bcc) crystals
also apply to hcp materials (within their basal plane), provided a re-
estimation of fitting parameters and re-scaling to account for the sixfold
symmetry. While grain growth competition remains stochastic with rare events
of unexpected elimination or side-branching along the developing GBs, we also
confirm an overall transition from a geometrical limit to a favorably oriented
grain limit behavior with an increase of thermal gradient within the dendritic
regime, and the progressive alignment of dendrites and GBs toward the
temperature gradient direction with an increase of $G$ during the dendritic-
to-cellular morphological transition. Comparisons with original thin-sample
directional solidification experiments show a qualitative agreement with PF
results, yet with notable discrepancies, which nonetheless can be explained
based on the stochastic variability of selected GB orientations, and the
statistically limited experimental sample size. Overall, our results extend
the understanding of GB formation and grain growth competition during
solidification of hcp materials, and the effect of thermal conditions,
nonetheless concluding on the challenges of extending the current studies to
three dimensions, and the need for much broader (statistically significant)
data sets of GB orientation selected under well-identified solidification
conditions.

###### keywords:

Solidification , Grain boundaries , Microstructure , Hcp crystal structure ,
Mg alloys.

##  1 Introduction

Metals and alloys represent a crucial class of structural materials in the
manufacturing sector. Their crystalline microstructures and mechanical
behavior are relatively well-identified and understood within a single crystal
[1, 2, 3, 4]. Yet, while most technological metallic materials are
polycrystalline, a fundamental understanding remains lacking upon key
mechanisms of the formation of heterogeneous polycrystalline microstructures
and upon their resulting mechanical behavior, in particular for materials with
complex hexagonal (hcp) crystal structures. Among these, Magnesium-based
alloys represent tremendous opportunities for weight reduction in automotive
applications [5, 6, 7], as well as a strong potential in aeronautics and, most
recently, in biodegradable implants [8, 9]. Mg-rare earth (Mg-RE) alloys, in
particular, have recently emerged with improved mechanical properties,
especially hardness, yield strength, tensile and creep resistance [8, 10, 11,
12]. The mechanical properties can be improved by heat treatment processes,
but also by controlling the microstructure during the solidification process,
especially when macrosegregation-related phenomena are involved [1, 2].

Microstructural characteristics, such as dendritic spacings and
polycrystalline grain texture, determine the properties, performance, and
lifetime of metal cast parts. Most often, solidification processing is the
first processing step accompanied by microstructure formation. Even though our
understanding of crystal growth and solidification as a whole has evolved
significantly over the past decades [13, 14, 3, 4], microstructure selection
mechanisms are still not completely understood. For instance, a comprehensive
picture of dendritic grain growth competition, i.e., of the selection of the
shape and orientation of grain boundaries (GBs), is still lacking, in
particular for hcp materials.

The overall objective of this work is to establish a comprehensive
understanding of microstructure selection and grain growth competition during
solidification of hcp alloys. To do so, we perform two-dimensional (2D) phase-
field (PF) simulations of directional solidification of a Mg-1wt$\%$Gd alloy,
using a thoroughly-validated quantitative PF model [15, 16] with an hcp
anisotropy vector formulation [17], accelerated via parallelization using
graphics processing units (GPUs). We perform a full mapping of the bi-
crystalline grain orientation, exploring a wide range of temperature gradient
$G$, and compare our results to current theories for the selection of primary
dendrite arm spacing (PDAS), dendritic growth orientation, and GB orientation
selection, so far established within the context of cubic (fcc, bcc) crystal
structures. We also compare our simulation results with original experimental
measurements in thin samples directionally solidified under well-identified
temperature gradient $G$ and growth velocity $V$.

##  2 State-of-the-art theories

Before presenting our methods (Sec. 3) and discussing our results (Sec. 4),
first we review the state-of-the-art, in particular theories for the selection
of PDAS (Sec. 2.1), growth orientation (Sec. 2.2), and GB orientation
selection during directional solidification of alloys (Sec. 2.3).

###  2.1 Selection of primary dendrite arm spacing

For a given set of processing parameters $(G,V)$, a wide range of dendritic
array spacings may be selected [18, 19, 20, 21, 22, 23, 24]. In the case of
dendrites growing in a direction tilted by an angle $\alpha$ with respect to
the temperature gradient $G$, we denote $\Lambda_{0}$ the primary dendrite arm
spacing (PDAS) measured in the normal direction to the temperature gradient,
while $\Lambda_{\alpha}$ stands for the PDAS measured normal to the dendrite
growth direction (later illustrated in Figure 2). Classical phenomenological
models typically predict $\Lambda_{0}$ in an array of dendrites aligned with
$G$, as a function of the processing parameters [25]. Gandin et al. [26]
proposed an orientation-dependent PDAS model

$\Lambda_{\alpha}(G,\alpha)=LG^{-c}\\{1+d[\cos(\alpha)^{-e}-1]\\}$ (1)

that accounts for the growth direction angle $\alpha$ with respect to $G$. The
prefactor $L$ depends on the alloy properties (e.g. phase diagram) and growth
velocity $V$ (here a constant) while $c$, $d$ and $e$ are used as fitting
parameters.

The resulting orientation-dependence of the PDAS was validated against
experimental measurements [26] and phase-field simulations [27] for thin-
sample and regular arrays [28] directional solidification of cubic-symmetry
systems (e.g transparent succinonitrile compounds). Yet, to the best of our
knowledge, it was not compared to either experimental or computational results
for systems with a hcp crystal structure.

###  2.2 Preferred crystal growth direction

At low velocity, just above the onset of morphological destabilization of a
planar solid-liquid interface, cellular patterns exhibit a mostly isotropic
shape and grow along the main temperature gradient direction. At higher
velocities, dendrite orientation transitions have been reported in various
systems as a function of alloy composition [29, 30, 31, 32] or cooling rates
[33]. Still, most crystalline solids have been reported to grow along well-
defined crystallographic directions. For cubic systems (fcc, bcc), the
preferred growth direction is usually $\langle 100\rangle$, while for
hexagonal symmetry (hcp) systems, both $\langle 10\bar{1}0\rangle$ (most
common in Zn [34, 35]) and $\langle 11\bar{2}0\rangle$ (most common in Mg [36,
37]) preferred growth directions have been reported.

Between cellular and dendritic regimes, the preferred growth direction
$\alpha$ shows a continuous transition from the orientation of the thermal
gradient to the orientation selected by the crystal structure. This transition
has been formulated as [38, 39]

$\frac{\alpha}{\alpha_{0}}=1-\frac{1}{1+f\,\rm{Pe}^{g}}$ (2)

expressing the ratio between the growth direction $\alpha$ and the crystal
lattice (e.g. $\langle 100\rangle$) orientation $\alpha_{0}$ as a function of
the spacing Péclet number, ${\rm Pe}=\Lambda_{0}V/D$. The latter stands as a
measure of the interactions between the neighboring dendrites via the solutal
diffusion field. Here, $\Lambda_{0}$ is the PDAS measured perpendicularly to
$G$, and $f$ and $g$ are two fitting parameters. At high Pe numbers, i.e.,
high PDAS and/or low $G$, $\alpha\approx\alpha_{0}$ (i.e. growth along
crystallographic directions in the dendritic regime), while $\alpha\approx
0^{\circ}$ at low Pe numbers, i.e., low PDAS and/or high $G$ (i.e. growth
aligned with $G$ in the cellular regime).

This progressive transition has been derived from observations in thin-sample
directional solidification experiments [38, 39] and was later validated
computationally using quantitative phase-field simulations [27, 40, 41, 28].
The latter have shown that the parameter $f$, initially proposed to be a
function of $\alpha_{0}$ [38, 39], could be considered a constant for a given
material, and that Eq. (2) remains valid for three-dimensional systems with
different spatial array ordering (see Supplementary Information of [28]).
However, these studies all focused on cubic-symmetry crystals and, to the best
of our knowledge, this growth orientation selection during the cellular to
dendritic transition has not been reported for hcp crystals so far.

###  2.3 Selection of GB orientations during solidification

The early classical theory of grain growth competition by Walton and Chalmers
[42] states that the grain best aligned with the thermal gradient direction,
which grows at a lower undercooling compared to grains with higher
misorientation with respect to $G$, cannot be eliminated by these higher-
misorientation neighbors. A common extrapolation of this argument leads to the
consideration that the GB orientation equals the growth orientation of the
best-oriented grain, i.e. the minimum $|\alpha|$ among the two grains. This is
known as the favorably oriented grain (FOG) criterion. The resulting GB
orientation $\theta$, illustrated by the $\theta(\alpha_{1},\alpha_{2})$ map
in Fig. 1a, is thus given by

$\theta_{\rm FOG}(\alpha_{1},\alpha_{2})=\begin{cases}\alpha_{1}&\quad\text{if
}|\alpha_{1}|\leq|\alpha_{2}|\\\ \alpha_{2}&\quad\text{if
}|\alpha_{2}|\leq|\alpha_{1}|\end{cases}$ (3)

where $\alpha_{1}$ and $\alpha_{2}$ are the growth directions of the competing
grains labelled 1 and 2. As illustrated in Figure 2, to be consistent with
previous publications [27, 43, 44, 45] we consider that grain 1 is above the
GB while grain 2 is below the GB if the growth direction goes from left to
right. Early experimental observations supported qualitatively this simple
description by Eq. (3), e.g. exhibiting progressive elimination of grains with
the highest misorientation with $G$ [46, 47]. However, phenomena contradicting
this simple theory, such as the unusual overgrowth of well-oriented dendrites
at converging GB, were observed experimentally [48, 49] and later reproduced
and explained with phase-field simulations [50, 51, 43].

Figure 1:  Grain boundary orientation angle $\theta$ shown as a function of
the misorientation of grain 1, $\alpha_{1}$, and grain 2, $\alpha_{2}$ (see
Fig. 2) using (a) the Favorably Oriented Grain (FOG) criterion (Eq. (3)) and
(b) the Geometrical Limit (GL) (Eq. (4)). For simplicity, both maps consider a
well-developed dendritic regime with $\alpha\approx\alpha_{0}$ (see Section
2.2).  Figure 2:  Configuration of the PF simulations showing two, red and
blue, grains growing from left to right (direction of the temperature
gradient, $G$). Periodic boundary conditions applied normally to the growth
direction onto to a bicrystalline domain (bottom) allow tracking two GBs, one
converging GB and one diverging GB. For each GB (top schematic), angles
$\alpha_{1}$ and $\alpha_{2}$ are defined as the orientation of the upper and
lower grains, respectively, and $\theta$ is the orientation of the GB. All
angles are measured counterclockwise with respect to the thermal gradient $G$
main direction (here aligned with the $x$ axis), i.e. here illustrated for
$\alpha_{1}>0$, $\alpha_{2}<0$ and $\theta>0$. In a grain, PDAS $\Lambda_{0}$
and $\Lambda_{\alpha}$ are measured perpendicular to $G$ and to the growth
direction, respectively.

Pineau et al. [44] proposed another asymptotic growth limit, termed the
geometrical limit (GL), which can be interpreted as considering continuous
side-branching, i.e. infinitely small secondary arm spacing, such that side-
branches instantaneously occupy the space between two grains. The resulting GB
orientation, illustrated by the $\theta(\alpha_{1},\alpha_{2})$ map in Fig.
1b, is thus expressed as [44]

$\theta_{\rm
GL}(\alpha_{1},\alpha_{2})=\begin{cases}&\tan^{-1}\Big{(}\frac{\cos\alpha_{1}+\sin\alpha_{1}-\sec\alpha_{2}}{\cos\alpha_{1}-\sin\alpha_{1}}\Big{)}\\\
&\quad\quad\text{if }(\alpha_{1}+\alpha_{2})\geq 0~{}\text{ and
}~{}(\alpha_{1}<\alpha_{2})\\\
&\tan^{-1}\Big{(}\frac{\cos\alpha_{2}+\sin\alpha_{2}-\sec\alpha_{1}}{\cos\alpha_{2}+\sin\alpha_{2}}\Big{)}\\\
&\quad\quad\text{if }(\alpha_{1}+\alpha_{2})<0~{}\text{ and
}~{}(\alpha_{1}<\alpha_{2})\\\
&\frac{\alpha_{1}+\alpha_{2}}{2}\quad\quad\text{if
}(\alpha_{1}\geq\alpha_{2})\end{cases}$ (4)

where the first and second cases refer to converging GBs
($\alpha_{1}<\alpha_{2}$) on either side of the $\\{\alpha_{1}=-\alpha_{2}\\}$
symmetric configuration, while the third case refers mostly to diverging GBs
($\alpha_{1}>\alpha_{2}$).

Selection maps of $\theta(\alpha_{1},\alpha_{2})$ are expected to be periodic
in both $\alpha_{1}$ and $\alpha_{2}$, with a period $90^{\circ}$ for fourfold
(fcc,bcc) or $60^{\circ}$ for sixfold (hcp in the basal plane) symmetry. This
periodicity appears in both FOG (Fig. 1a) and GL (Fig. 1b) maps. Moreover,
both maps exhibit a discontinuity along the $\\{\alpha_{1}=-\alpha_{2}\\}$
diagonal. Indeed, from symmetry considerations, once averaged over a
statistically significant number of GBs, it is expected that
$\theta(\alpha_{1}=-\alpha_{2})=0^{\circ}$. The same discontinuity, though
less visible in Fig. 1, occurs at the $\\{\alpha_{1}=\alpha_{2}=\pm
30^{\circ}\\}$ corners. This latter discontinuity is even more pronounced when
differences between $\alpha$ and $\alpha_{0}$ are considered (cf. Sec. 2.2).
While FOG and GL theories both consider dendritic microstructures, actual
microstructures usually transition toward a degenerate (i.e., seaweed-like)
patterns as they approach $\alpha_{0}\to\pm 30^{\circ}$ ($\pm 45^{\circ}$ in
fourfold symmetry systems). The discontinuity at the four corners of the
$\theta(\alpha_{1},\alpha_{2})$ plot is thus usually replaced by a transition
region within a $(\alpha_{1},\alpha_{2})$ range directly linked to the
dendritic-to-degenerate transition region [27, 43].

GB orientation data across a broad range of crystalline orientations
$(\alpha_{1},\alpha_{2})$ is tedious to gather experimentally and hence
scarce, such that the selection of $\theta(\alpha_{1},\alpha_{2})$ has been
approached primarily by phase-field simulations. The most extensive mapping of
$(\alpha_{1},\alpha_{2})$ space in bi-crystalline samples to date [43] has
shown that the GB orientation selection deviates from the classical FOG
description, Eq. (3), in the following ways: (i) $|\theta|$ remains low for
diverging GBs with $\alpha_{1}\times\alpha_{2}<0$ (hence closer to the GL in
the bottom right quadrant of Fig. 1), (ii) the PF selection map exhibits a
narrow transition region along the symmetric $\\{\alpha_{1}=-\alpha_{2}\\}$
configuration, and (iii) a transition to degenerate seaweed-like patterns may
occur toward the $\\{\alpha_{i}=45^{\circ}\\}$ ($i=1,2$) symmetric
configuration (in fourfold cubic symmetry crystals, hence with singularity at
$45^{\circ}$).

Several studies have highlighted the important role of the temperature
gradient $G$ [27, 51, 45], in particular close to the cellular-to-dendritic
transition regime [27, 51], and also suggestted a transition from mostly-FOG
to mostly-GL behavior with a decrease in $G$ in the dendritic regime [45].
Recent studies have also highlighted the morphological complexity of three-
dimensional (3D) GBs [52, 28], as observed in directionally solidified samples
[53]. However, given the added dimensionality of the 3D GB bi-crystal problem,
the computational cost of 3D simulations still stands as bottleneck to perform
a comprehensive mapping of grain orientations and processing conditions
$(G,V)$. Moreover, except a handful of phase-field simulations by Eiken [54,
55] (to our knowledge, the earliest reports of PF-simulated grain growth
competition), the problem of bi-crystal grain growth competition in hcp
crystals remains for the most part unexplored.

##  3 Methods

Here, for the sake of reproducibility, we present details on the PF model
(Sec. 3.1) and simulations (Sec. 3.2), as well as solidification experiments
(Sec. 3.3). With the exception of the hcp six-fold symmetry of the interface
excess free energy (implemented using an auxiliary anisotropy vector [17]),
the model and simulation setup are similar to what was already employed in
previous studies [27, 43].

###  3.1 Phase-field model

We use the classical quantitative thin-interface phase-field model for
directional solidification of a binary alloy with anti-trapping current [16],
a preconditioned phase field $\psi=\sqrt{2}\tanh^{-1}(\varphi)$ [56, 27], and
a formulation and implementation including the anisotropy vector $\vec{A}$ to
facilitate the use of arbitrary interface energy anisotropy functions [17].
The resulting evolution equations for the preconditioned phase field $\psi$
and dimensionless solute concentration (supersaturation) field $U$ are

$\begin{split}\Bigg{[}1-&(1-k)\frac{x-\tilde{V}t}{\tilde{l_{T}}}\Bigg{]}\frac{\partial\psi}{\partial
t}=\\\
&\sqrt{2}\Bigg{[}\varphi-\lambda\Big{(}1-\varphi^{2}\Big{)}\Big{(}U+\frac{x-\tilde{V}t}{\tilde{l}_{T}}\Big{)}\Bigg{]}\\\
&+a_{s}^{2}\Big{[}\nabla^{2}\psi-\sqrt{2}\varphi|\vec{\nabla}\psi|^{2}\Big{]}\\\
&+2a_{s}\vec{\nabla}a_{s}\vec{\nabla}\psi+\frac{\sqrt{2}}{1-\varphi^{2}}\vec{\nabla}\cdot\vec{A}\end{split}$
(5)

and

$\displaystyle\begin{split}[1+&k-(1-k)\varphi]\frac{\partial U}{\partial
t}=\tilde{D}\vec{\nabla}\Big{[}(1-\varphi)\vec{\nabla}U\Big{]}\\\
&+\Big{[}1+(1-k)U\Big{]}\frac{(1-\varphi^{2})}{\sqrt{2}}\frac{\partial\psi}{\partial
t}\\\
&+\vec{\nabla}\left\\{\Big{[}1+(1-k)U\Big{]}\frac{(1-\varphi^{2})}{2}\frac{\partial\psi}{\partial
t}\frac{\vec{\nabla}\psi}{|\vec{\nabla}\psi|}\right\\}\end{split}$ (6)

where

$U=\frac{1}{1-k}\left[\frac{c/c_{l}^{0}}{(1-\varphi)/2+k(1+\varphi)/2}-1\right]$
(7)

with $c_{l}^{0}=c_{\infty}/k$ the solute concentration on the liquid side of a
flat interface at the solidus temperature $T_{0}$ for an alloy of nominal
concentration $c_{\infty}$. This model considers a linearized phase diagram
(or dilute alloy limit) with constant liquidus slope $m$ and interface solute
partition coefficient $k$, a solute diffusion coefficient $D$ in the liquid
phase and negligible diffusion in the solid phase. The temperature field
follows the common “frozen temperature approximation” with constant
temperature gradient $G$ and pulling (i.e. isotherms) velocity $V$ in the $x$
direction, resulting in the one-dimensional temperature field
$T(x)=T_{0}+G(x-Vt)$. The interface thickness $W$ and the phase field equation
relaxation time $\tau_{0}$ are taken as space and time units respectively [15,
16]. Thus, the dimensionless form of the diffusion coefficient $\tilde{D}$,
the growth velocity $\tilde{V}$ and the thermal length $\tilde{l}_{T}$ are
given respectively by

$\tilde{D}=\frac{D\tau_{0}}{W^{2}}=a_{1}a_{2}\frac{W}{d_{0}},$ (8)
$\tilde{V}=\frac{V\tau_{0}}{W}=a_{1}a_{2}\frac{Vd_{0}}{D}\Big{(}\frac{W}{d_{0}}\Big{)},$
(9)

and

$\tilde{l}_{T}=\frac{l_{T}}{W}=\frac{l_{T}}{d_{0}}\frac{1}{W/d_{0}},$ (10)

with $d_{0}=\Gamma/[mc_{\infty}(1-1/k)]$ the chemical capillary length at
$T_{0}$, $\Gamma$ the interface Gibbs-Thomson coefficient,
$l_{T}=[mc_{\infty}(1-1/k)]/G$ the thermal length, $a_{1}=5\sqrt{2}/8$ and
$a_{2}=47/75$ [15, 16]. In the bulk liquid ($\varphi=-1$), Eq. (6) reduces to
the solute diffusion equation, whereas at the interface using a coupling
coefficient $\lambda=a_{1}W/d_{0}$ between the phase field and the
supersaturation ensures a vanishing kinetic undercooling [57, 16]. In Eq. (5),
$a_{s}(\vec{n})=\gamma(\vec{n})/\gamma_{0}$ expresses the anisotropy of the
interface excess free energy projected on the basal plane, formulated as [58,
55]

$a_{s}(\vec{n})=\Big{[}1.008\,+\,\varepsilon_{66}\frac{\sqrt{6006}}{64\sqrt{\pi}}(n_{x}^{6}-15n_{x}^{4}n_{y}^{2}+15n_{x}^{2}n_{y}^{4}-n_{y}^{6}))\Big{]}$
(11)

where $\varepsilon_{66}$ is the strength of the surface energy anisotropy in
the basal plane and $(n_{x},n_{y})$ is the normal vector to the solid-liquid
interface calculated as
$\vec{n}=-\vec{\nabla}\psi/\lvert\vec{\nabla}\psi\lvert$. In Eq. (11), we
consider the two-dimensional sixfold anisotropy in the basal plane with an
out-of-plane anisotropy coefficient $\varepsilon_{20}\approx-0.026$ [58]. The
auxiliary anisotropy vector $\vec{A}$ is introduced for convenience and
versatility, as it makes it easier to switch between different forms of the
interface anisotropy [17, 59]. It appears from reformulating the anisotropy
(i.e., last) term of the phase-field equation as a divergence

$\displaystyle\vec{\nabla}\cdot\vec{A}=\sum_{i=x,y}{\partial_{i}\left[|\vec{\nabla}\varphi|^{2}a_{s}\frac{\partial
a_{s}}{\partial\varphi_{i}}\right]}$ (12)

with

$\displaystyle A_{i}$
$\displaystyle=|\vec{\nabla}\varphi|^{2}a_{s}\frac{\partial
a_{s}}{\partial\varphi_{i}}$
$\displaystyle=\frac{1-\varphi^{2}}{\sqrt{2}}|\vec{\nabla}\psi|^{2}a_{s}\frac{\partial
a_{s}}{\partial\psi_{\alpha}}$
$\displaystyle=\frac{1-\varphi^{2}}{\sqrt{2}}|\vec{\nabla}\psi|a_{s}\sum_{j=x,y}(n_{i}n_{j}-\delta_{i,j})\frac{\partial
a_{s}}{\partial n_{j}}$ (13)

with $\delta_{i,j}=1$ if $i=j$ and $0$ otherwise. Details and full expressions
of $\vec{A}$, as they appear in the code, are provided in A (for completeness,
also including 3D expressions).

In order to destabilize the initial planar interface and promote side-
branching, at each grid point ($i,j$) and each time iteration, we add a random
perturbation $\eta\beta_{i,j}\sqrt{\Delta t}$ to the value of $\psi_{i,j}$ at
the next time step :

$\psi_{i,j}(t+\Delta t)=\psi_{i,j}(t)+\Delta
t\frac{\partial\psi_{i,j}}{\partial t}+\eta\beta_{i,j}\sqrt{\Delta t}$ (14)

where $\beta_{i,j}$ is a random number picked from a flat distribution in the
range $[-0.5,+0.5]$ and $\eta$ is the noise amplitude [60, 61].

Equations (5) and (6) are solved using finite differences on a regular grid of
square elements with a spatial step $\Delta x$. They are integrated following
an Euler explicit time scheme with a time step $\Delta t=0.8\Delta
x^{2}/4\tilde{D}$. We use a moving domain following the most advanced dendrite
tip in the $x$ direction. For performance, Eq. (13) is set to zero in the bulk
phases, i.e. namely where $(1-\varphi^{2})<0.01$. The model is written in
computer uniﬁed device architecture (CUDA) programming language for massively
parallel computing on Graphic Processing Units (GPU). The time loop is
composed of three sequential kernel calls, calculating: (i) the $\vec{A}$ ﬁeld
at the current time step, (ii) the $\psi$ ﬁeld at the next time step, and
(iii) the $U$ ﬁeld at the next time step. An additional kernel shifting all
fields in the $x-$ direction is triggered occasionally when the highest
location of the interface in $x$ exceeds a given position $L_{\rm solid}$ in
the computational domain. Time stepping is performed at the end of the time
loop by swapping pointer addresses between arrays containing values of $\psi$
and $U$ at the current and next time steps.

###  3.2 Simulations

We study grain growth competition in a Mg-1wt%Gd alloy during directional
solidification under a constant temperature gradient $G$ and constant cooling
rate, or equivalent pulling (i.e. isotherms) velocity $V$. This velocity is
set to 166 µm/s to match the estimated velocity of the isotherms in the
experiment (Sec. 3.3) and we scan a broad range of temperature gradients from
0.3 to 10 K/mm.

Alloy and simulation parameters are summarized in Table 1. Phase diagram
features (i.e. temperatures, concentrations, and resulting $k$ and $m$) for a
dilute Mg-Gd alloy were taken from ASM handbooks [62, 63]. The solid-liquid
interface excess free energy, $\gamma_{0}=89.9~{}$mJ/m2, and its anisotropy
coefficients, $\varepsilon_{66}\approx 0.003$ and
$\varepsilon_{20}\approx-0.026$, come from molecular dynamics calculations for
pure Mg [58]. The Gibbs-Thomson coefficient was also calculated for pure Mg as
$\Gamma=\gamma_{0}T_{M}/(\rho L_{f})$, with a density $\rho=1.59\times
10^{6}~{}$g/m3 and a latent heat of fusion $L_{f}=349~{}$J/g [63]. The
diffusion coefficient for Gd in liquid Mg was approximated as $D\approx D_{\rm
Mg}r_{\rm Gd}/r_{\rm Mg}$ [64], considering Mg self-diffusion coefficient
$D_{\rm Mg}\approx 6.3\times 10^{-9}~{}$m2/s at $T=1000~{}$K [65] and atomic
radii $r_{\rm Gd}=233$ pm and $r_{\rm Mg}=145$ pm [66].

Table 1:  Alloy, process, and simulation parameters.  Parameter | Symbol | Value | Unit
---|---|---|---
Alloy Gd concentration | $c_{\infty}$ | $1.0$ | wt% Gd
Solute partition coefficient | $k$ | $0.4$ | -
Liquidus slope | $m$ | $-1.66$ | K/wt% Gd
Pure solvent (Mg) melting temperature | $T_{M}$ | $923$ | K
Alloy liquidus temperature | $T_{L}$ | $921.34$ | K
Alloy solidus temperature | $T_{0}$ | $918.85$ | K
Anisotropy strength in the basal plane [58] | $\varepsilon_{66}$ | $0.003$ | $-$
Liquid Gd diffusion coefficient [65, 64, 66] | $D$ | $3.9\times 10^{-9}$ | m2/s
Gibbs-Thomson coefficient [63, 58] | $\Gamma$ | $1.5\times 10^{-7}$ | K.m
Solute capillary length at $T_{0}$ | $d_{0}$ | $0.06$ | µm
Pulling (isotherms) velocity | $V$ | 166 | µm/s
Thermal gradient | $G$ | $0.3-10.0$ | K/mm
Diffuse interface width | $W$ | $24.74$ | $d_{0}$
Grid step size | $\Delta x$ | $1.0$ | $W$
|  | $=1.5$ | µm
Noise amplitude [60, 61] | $\eta$ | 0.01 | $-$

The initial configuration of the system is close to that used in previously
reported phase-field studies of grain growth competition [27, 43]. Two grains
A and B are placed on the left (cold) side of a two-dimensional (2D) domain.
We use the notation A and B only to refer to two grains within a simulation
(see, e.g., bottom of Fig. 2), while we use numbers 1 and 2 to denote grains
for each individual GB (top of Fig. 2). In the latter, by convention, grain 1
is above and grain 2 is below the GB if the growth direction goes left-to-
right. Hence, in each simulation containing two GBs, grains A and B
alternatively play the role of grains 1 and 2 (or vice-versa) depending on
which GB is analyzed.

The nearly-flat solid-liquid interface is initialized at a given undercooling
($\Delta=0.15$) located at $x=L_{\rm solid}$ in the simulation domain. Grains
A and B initially share equally the solid domain, with GBs located along the
$(y=L_{y}/4)$ and $(y=3L_{y}/4)$ lines. In order to break the one-dimensional
symmetry of the system, the solid-liquid interface location has a small
sinusoidal perturbation of amplitude $\Delta x$ and 5 periods (i.e. bumps)
across the $y$ dimension. The preconditioned phase field $\psi$ is initialized
as a signed distance function from the interface location, with $\psi>0$ in
the solid ($x<L_{\rm solid}$) and $\psi<0$ in the liquid ($x>L_{\rm solid}$).
The concentration field, $c$, is initialized with equilibrium concentrations
on both sides of the interface, homogeneous $c$ in the solid, and in the
liquid an exponential profile $c\sim\exp\\{xV/D\\}$ decaying toward
$c_{\infty}$ when $x\to\infty$. The $c$ ($\langle 0001\rangle$) axis of both
grains is perpendicular to the simulation plane, and their $a$ ($\langle
10\bar{1}0\rangle$) axis makes an angle $\alpha_{A}$ or $\alpha_{B}$,
respectively, with the $G$ direction, i.e. with the $x$ axis (see Fig. 2). We
use an auxiliary integer field to track the grain identity ($-1$ in grain A,
$+1$ in grain B, and $0$ in the liquid). The grain index field is propagated
using the sign of the average over the eight neighboring points whenever a
liquid point reaches $(1-\varphi)^{2}>0.01$ [27, 43].

Simulation sizes are summarized in Table 2. The domain length $L_{x}$ is
chosen long enough (i) to track the emergence of tertiary branches on the
solid side and (ii) to allow the solute field to relax to $c_{\infty}$ on the
liquid side. The domain width $L_{y}$ is chosen long enough to host at least
ten PDAS. Hence, $L_{y}$ is higher for low $G$. No-flux and periodic boundary
conditions are imposed along $x$ and $y$ directions, respectively. The most
advanced interface location in the $x$ direction is maintained at a distance
$L_{\rm solid}$ of the leftmost ($x-$) domain boundary as listed in Table 2.

Table 2:  Physical dimensions of the system for different thermal gradients and corresponding numbers of mesh points in the numerical domain.  $G$ (K/mm) | $N_{x}\times N_{y}$ |  $L_{x}\times L_{y}$ (${\rm mm}^{2}$) |  $L_{\rm solid}$ (µm) | Time (s)
---|---|---|---|---
0.3, 0.5, 1.0, 2.0, 3.0 | $960\times 2560$ | $1.4\times 3.8$ | 1150 | 250
10.0 | $576\times 1024$ | $0.9\times 1.5$ | 600 | 80

This configuration, illustrated in Fig. 2, allows tracking the evolution of
one converging GB and one diverging GB for each simulation. At each GB, the
upper ($y+$) and lower ($y-$) grains are labelled as grain 1 and grain 2,
respectively. At the beginning of each simulation, we impose internal walls
(with a no-flux condition in the $y$-direction on both sides) separating the
grains in order to discard early transient interactions and thus focus on
grain growth competition at a steady velocity $V$ (see Supplementary Materials
of [43]). Once both grains reach a steady tip undercooling, we gradually
remove the “wall” separating the two grains and start recording the GB
location. Simulation time (Table 2) is selected long enough to track at least
ten branching (diverging GB) and ten elimination (converging GB) events or
until one grain is eliminated. As the moving domain follows the most advanced
primary tips, we store the coordinates of the two GBs between grains A and B
at the grid points leaving the domain on the $x-$ side. The average
orientation $\theta$ of each GB is thus extracted via a linear least-square
fitting of the GB location points, automated using gnuplot scripts. The
lateral ($y$) mobility of the GB, or its roughness, is thus represented by the
standard deviation of the linear fit.

We aim to scan the full range of bi-crystal configurations
$(\alpha_{1},\alpha_{2})$ for different temperature gradients. Accounting for
the crystal sixfold symmetry in the basal plane, we scan all combinations
($\alpha_{A},\alpha_{B}$) between $-30^{\circ}$ and $+30^{\circ}$ with steps
of $5^{\circ}$ for six values of $G$ from 0.3 to 10 K/mm with $V=166~{}$µm/s
(see Fig. 2 and Table 2). To account for the stochasticity of GB orientation
selection, in particular with respect to side-branching, we perform ten
simulations for each configuration using different random number seeds for
$\beta_{i,j}$ in Eq. (14). Simulations with $\alpha_{1}=\alpha_{2}$ were only
performed once, and primarily used to extract the growth orientation $\alpha$
for a given crystal orientation (in such case, $\alpha\approx\theta$ is
directly obtained via the same least-square fitting) as well as the primary
spacing. Simulation results are gathered in maps of the GB angle $\theta$ as a
function of $(\alpha_{1},\alpha_{2})$ similar to those represented for FOG and
GL maps in Fig. 1. Taking advantage of symmetries, one can construct the
entire map by scanning about one quarter of the $(\alpha_{1},\alpha_{2})$
space (also see Supplementary Material of [43]).

As summarized in Table 3, for each $G$, the resulting map corresponds to 10
runs for 42 $(\alpha_{A},\alpha_{B})$ configurations with
$\alpha_{A}\neq\alpha_{B}$, while the six configurations with
$\alpha_{A}=\alpha_{B}$ were performed only once. The full mapping for the six
explored $G$ thus amounts to $6\times 426=2556$ simulations. Each simulation
was performed on one Nvidia GPU (GeForce RTX 3090, RTX 2080Ti, or Titan Xp).
The average compute time for one simulation was about 62 minutes, with
individual simulations completed between 6.5 minutes (in the case of quick
elimination of a grain) and 3.5 hours (for low $G$ and, typically,
$\alpha_{1}\approx\alpha_{2}$). The cumulative simulation wall time of all
single-GPU-parallelized simulations reported here thus amounts to over 2640
hours (i.e., 110 days).

Table 3:  Number of simulations performed for each $(\alpha_{A},\alpha_{B})$ configuration.  | $\alpha_{B}(^{\circ})$
---|---
$\alpha_{A}(^{\circ})$ | 30 | 25 | 20 | 15 | 10 | 5 | 0
$-30$ | 10 |  |  |  |  |  |
$-25$ | 10 | 10 |  |  |  |  |
$-20$ | 10 | 10 | 10 |  |  |  |
$-15$ | 10 | 10 | 10 | 10 |  |  |
$-10$ | 10 | 10 | 10 | 10 | 10 |  |
$-5$ | 10 | 10 | 10 | 10 | 10 | 10 |
0 | 10 | 10 | 10 | 10 | 10 | 10 | 1
5 | 10 | 10 | 10 | 10 | 10 | 1 |
10 | 10 | 10 | 10 | 10 | 1 |  |
15 | 10 | 10 | 10 | 1 |  |  |
20 | 10 | 10 | 1 |  |  |  |
25 | 10 | 1 |  |  |  |  |

###  3.3 Experiments

While the article is mainly focused on simulation results, we also performed
experiments in order to discuss, specifically, the selection of GB orientation
during directional solidification of a Mg-1 wt$\%$Gd alloy. We built a custom
furnace for directional solidification of thin samples. The furnace was
designed to later allow in-situ X-ray radiography imaging, but the experiments
discussed here were completed without in-situ imaging, i.e. analyses were
performed ex-situ after directional solidification.

The furnace, illustrated in Fig. 3, is composed of pairs of identical heaters
for “cold” and “hot” zones separated by 6.8 mm from each other. Between the
heaters, we insert a graphite crucible consisting of two rectangular plate-
like halves and a third one with a central rectangular deep pocket (Fig. 3a).
A sample of thickness 1 mm is covered in boron nitride, placed into the
pocket, and sandwiched between the two crucible plates. The temperature field
is monitored by four thermocouples located at different heights and embedded
in the central layer of the crucible. The resulting temperature measurements
enable the feedback control of the independently controlled flat heaters, as
well as the evaluation of the thermal conditions during the solidification
experiment.

Figure 3:  Experimental setup for directional solidification of thin sample
under Ar gas atmosphere.

As illustrated in Fig. 3b, the system (crucible + sample + heaters +
thermocouples) is inserted into an outer chamber. Compressed air is flowing
between the two chambers, providing thermal insulation and minimizing damage
to electrical connections. Special feedthrough assemblies were designed to
provide an efficient seal on the wires passing through the inner chamber and
to restrain them from moving under pressure and vacuum.

We performed directional solidification experiments on thin samples ($15\times
10\times 1$ mm3) of Mg-1 wt$\%$Gd. Each experiment starts by setting an
initial thermal gradient $G$, heating the sample until it melts completely,
then cooling down both heaters at a constant rate $R$ while maintaining the
thermal gradient $G$ during the process. We aimed for a value of $G\approx
3~{}$K/mm and explored cooling rates ranging from 2 to 40 K/min Here, we
specifically focus on the analysis of one sample that exhibited a few GBs
involving grains with a low misorientation of their $c$ axis and the normal
direction to the sample plane (hence representative of the 2D simulation
configuration). The thermocouple readings indicate that the sample was
solidified under a temperature gradient $G\approx 3~{}$K/mm and a velocity of
isotherms $V\approx 166~{}$µm/s.

Electropolishing was performed on the solidified specimen after gently
grinding with 1200 grit size abrasive SiC paper. The surface of the sample was
electropolished at $-30^{\circ}$C with a voltage of 50 V using a solution of
15 vol$\%$ nitric acid and 85 vol$\%$ ethanol for 30 s. The electropolished
sample was then cleaned by immersion in an ultrasound bath of pure ethanol for
2 minutes. The resulting sample was smooth enough for microstructural
characterization in a scanning electron microscope (SEM).

Electron backscatter diffraction (EBSD) was performed to obtain
crystallographic orientation maps using a ﬁeld-emission gun SEM (Helios
Nanolab 600i) equipped with an Oxford-HKL EBSD detector. The scanning step
size was selected to obtain sufficient spatial resolution to identify GBs. The
accelerating voltage was set to 20 kV and the current to 2.7 nA. An automated
mapping was performed that allowed joining several EBSD maps with an imposed
overlap (usually 10-20$\%$) to obtain a wide EBSD map. The diffraction data
was post-processed using Matlab toolbox MTEX [67]. The maps were cleaned using
a spline ﬁlter to ﬁll non-indexed points with the nearby pixels and eliminate
measurement errors. The threshold misorientation angle for grain
reconstruction was set to 2°.

The EBSD map was then smoothed and individual GBs were studied in order to
select regions of interest for further analysis. In particular, we focused on
regions containing GBs between grains whose $c$ axis deviation from the $z$
axis (noted $\Phi$) is lower than $25^{\circ}$ (filtered using Matlab MTEX
toolbox). Such GBs are the ones most comparable with our 2D simulations with
sixfold anisotropy within the basal plane. Each GB orientation was then
calculated by fitting the location of GB points to straight lines with gnuplot
using the same convention as in Fig. 2.

##  4 Results and discussion

###  4.1 Differences between cubic and hexagonal systems

Before presenting a quantitative analysis of our results, we highlight some of
the main differences expected between a cubic system, such as those analyzed
in previous articles [27, 43, 45], and the hcp alloy system studied here.
These differences are illustrated in Fig. 4, comparing one of the simulations
from Ref. [43] and one of the current simulations, both of them with
$\alpha_{A}=-10^{\circ}$ (blue grain) and $\alpha_{B}=+15^{\circ}$ (red
grain).

Figure 4:  Illustration of grain growth competition mechanisms at converging
(top) and diverging (bottom) GBs, highlighting differences between cubic
(left) and hexagonal (right) alloy systems. The left-hand-side is one of the
simulations analyzed in Ref. [43], corresponding to a succinonitrile-0.4
wt%acetone alloy at $G=30$ K/cm and $V=25$ µm/s, using simulation parameters
from Refs [27, 43]. The right-hand-side simulation is one of the simulations
analyzed here. Both simulations have $\alpha_{A}=-10^{\circ}$ (blue grain) and
$\alpha_{B}=+15^{\circ}$ (red grain).

As mentioned in Sec. 3.2, a first difference is the range of angles that needs
to be mapped, namely $[-45^{\circ};+45^{\circ}]$ for cubic systems and
$[-30^{\circ};+30^{\circ}]$ for hexagonal ones. In terms of growth
competition, the sixfold symmetry does not lead to a substantial difference at
converging GBs (top GBs in Fig. 4) controlled by dendrite impingement. In
contrast, diverging GBs (bottom GBs in Fig. 4), piloted by side-branching
competition, are expected to exhibit more prominent differences due to the
distinct angle between primary dendrites and higher-order side-branches
(illustrated by purple arrows in red grains of Fig. 4). On the one hand, in
cubic (fourfold) systems, the angle between primary and secondary branches is
$\approx 90^{\circ}$. As a result, for low misorientation angles, side-
branches face a relatively small temperature gradient in their respective
growth direction. In the extreme case of a growth direction of primary
dendrites oriented along the temperature gradient, secondary branches
essentially grow along isotherms. On the other hand, in hexagonal (sixfold)
systems, secondary branches grow at an angle $\approx 60^{\circ}$ with respect
to the primary dendrites. Consequently, even at low misorientation angle,
side-branches face a relatively higher temperature gradient in their
respective growth direction, which is expected to lead to increased side-
branching activity compared to cubic systems.

Note that the angle between the growth directions of primary and secondary
branches is most often not exactly $90^{\circ}$ or $60^{\circ}$, due to the
deviation from crystal orientation ($\alpha\neq\alpha_{0}$) described in Sec.
2.2 and discussed in Sec. 4.2.

###  4.2 Primary dendrite arm spacing

Fig. 5 shows the variation of selected PDAS $\Lambda_{\alpha}$ in PF
simulations (symbols) as a function of (a) $G$ and (b) $\alpha_{0}$ together
with a set of fitted curves (lines) using a slightly modified version of Eq.
(1):

$\Lambda_{\alpha}(G,\alpha_{0})=LG^{-c}\\{1+d[\cos({45\alpha_{0}}/{30})^{-e}-1]\\}$
(15)

that includes an empirical correction factor $45/30$ in the cosine function to
account for the sixfold symmetry of the hcp crystal. For simplicity, we
express $\Lambda_{\alpha}$ directly as a function of the crystal orientation
$\alpha_{0}$, rather than the dendrite growth direction $\alpha$, considering
that this will likely only affect the value of the fitted coefficients but not
the validity of the relation. All curves correspond to Eq. (15) with $L\approx
0.212\,$mm0.66K0.34, $c\approx 0.34$, $d\approx 0.55$ and $e\approx 3.0$,
obtained from a global fit of all data points for $0.3\leq G\leq 3.0$ K/mm.
The highest gradient, $G=10$ K/mm, was excluded from the fit on purpose, as it
leads to intermediate patterns between dendritic and cellular (see Sec. 4.4),
hence not entirely relevant to the $\Lambda_{\alpha}(G,\alpha_{0})$ law
derived specifically for dendritic arrays.

Figure 5:  Primary arm spacing ($\Lambda_{\alpha}$) as a function of (a) the
thermal gradient $G$ for different crystal orientations $\alpha_{0}$, and (b)
$\alpha_{0}$ for different $G$ values, showing PF results (symbols) against
fitted Eq. (15) (lines) with $L=0.212\,$mm0.66K0.34, $c=0.34$, $d=0.55$ and
$e=3.0$. In both plots, colors represent the crystal orientation $\alpha_{0}$
from $0^{\circ}$ (dark blue) to $25^{\circ}$ (dark red) while symbols
represent temperature gradient $G$ from 0.3 K/mm (big $\bigtriangledown$) to
3.0 K/mm (big $\bigtriangleup$).

The overall variation of $\Lambda_{\alpha}$ with $G$ and $\alpha_{0}$ is in
good agreement with previous works in cubic systems [26, 27]. Indeed, spacings
decrease with increasing $G$ with $\Lambda\sim G^{-0.34}$ (Fig. 5a) and they
increase with increasing $\alpha_{0}$ (Fig. 5b). Due to the low number of
measured spacings (and hence limited statistical significance), we did not
perform a fully quantitative analysis of spacing distributions, nor of the
full PDAS stability limits, as, e.g., in [68, 69, 70] (out of scope here).
However, our observations suggest that the distribution of spacing is wider
for lower values of $G$. This supports our earlier results on Al-Cu alloys
[70], which showed that the spacing stability range $[\Lambda_{\rm
min},\Lambda_{\rm max}]$ becomes wider as $G$ decreases (in the dendritic
regime) and that $\Lambda_{\rm max}/\Lambda_{\rm min}\to 2$ as $G$ increases
(toward the cellular regime). Overall, Fig. 5 illustrates that
$\Lambda_{\alpha}$ seems relatively weakly affected by $\alpha_{0}$ for
moderate misorientations (at most 25 $\%$ for $\alpha_{0}\leq 25^{\circ}$) in
comparison to its strong dependence upon $G$.

The effect of $G$ observed here appears weaker compared to previous results
for cubic-symmetry alloys [71, 18, 26, 70]. Indeed, at a given $\alpha_{0}$,
varying $G$ by an order of magnitude induces a variation of $\Lambda_{\alpha}$
(i.e., Pe) by a factor of $10^{0.34}\approx 2.2$, while it typically does so
by a factor of $10^{0.5}\approx 3.2$ for alloys with a cubic symmetry. This
weaker dependence upon $G$ might be attributed to the fact that secondary
arms, growing at an angle of 60∘ from the primary trunks, face a higher
temperature gradient in their growth direction than in cubic systems, for
which secondary arms might even grow along isotherms when
$\alpha_{0}=0^{\circ}$. In the case of hcp systems, the interactions between
dendrites through the solute field might then be more intense than in cubic
systems, and the emergence of tertiary branches as new primary dendrites might
also occur earlier (i.e. at lower $\Lambda$), resulting in a PDAS stability
domain narrower in hcp systems. While PDAS selection is not directly the main
focus of the present article, it might be of interest to address these
interpretations in future studies, i.e. by comparing PDAS stability domains
for model alloys with similar thermophysical parameters but different
crystalline symmetries.

###  4.3 Dendritic growth direction

Fig. 6 shows the ratio $\alpha/\alpha_{0}$ between growth direction, $\alpha$,
and crystal orientation, $\alpha_{0}$, measured in the PF simulations as a
function of the spacing Péclet number, ${\rm Pe}=\Lambda_{0}V/D$, for
$G=[0.3-10.0]$ K/mm (symbols) and a least-square fit (dashed black line) of
Eq. (2).

Figure 6:  Dendrite growth orientation (i.e. ratio of dendrite growth
direction angle $\alpha$ by crystal orientation $\alpha_{0}$) versus Péclet
number Pe$=\Lambda_{0}V/D$, compared to fitted Eq. (2) with $f=0.708$ and
$g=0.854$ (dashed line). Consistently with Fig. 5, different symbols represent
different temperature gradients, $G$, while colors stand for crystal
orientations, $\alpha_{0}$ .

Overall, our results show that the dependence of growth direction established
for cubic systems also applies for Mg-1wt%Gd (hcp system), with all data
points relatively well fitted by Eq. (2), here for $f=0.708$ and $g=0.854$.
Similarly as in cubic systems [38, 39, 61, 27, 45], $\alpha\to\alpha_{0}$ as
Pe increases, i.e. as $G$ decreases, with a microstructure pattern within the
well-developed dendritic regime at low $G$ (i.e. high Pe). Approximating
$f(\alpha_{0})$ by a constant value also seems to remain a valid assumption.
Fitted parameters, however, differ from those assessed in cubic systems (e.g.
$f\approx 0.2$ and $g\approx 1.7$ in [27]).

For a given temperature gradient, values of the ratio $\alpha/\alpha_{0}$ fall
within a narrow range, as summarized in Table 4. Notably, all temperature
gradients exhibit a ratio $\alpha/\alpha_{0}$ higher than 0.78 in all
simulations/configurations, except for the highest $G$ at 10 K/mm. This is a
clear indication of the latter belonging to a transition regime between
dendrites ($\alpha\approx\alpha_{0}$) and cells ($\alpha\approx 0$), while all
$G\leq 3~{}$K/mm cases are within the dendritic regime. As patterns transition
toward the cellular regime, microstructures actually exhibit increasingly
shallower cells and transition to a seaweed pattern as $\alpha_{0}$ increases,
thus preventing the identification and measurement of clear primary spacings
(hence only reported for $\alpha_{0}\leq 15^{\circ}$ at $G=10~{}$K/mm in Fig.
6). Moreover, the overall grouping of data points for $G\leq 3~{}$K/mm may
also be related to the weaker dependence of $\Lambda_{0}$ on $G$ than in cubic
systems, mentioned in Sec. 4.2, resulting in a narrower range of Pe numbers.

Table 4:  Ratio between growth angle $\alpha$ and crystal orientation angle
$\alpha_{0}$ for each temperature gradients $G$.
$G$ | $\alpha/\alpha_{0}$
---|---
(K/mm) | Min | Max | Mean
0.3 | 0.835 | 0.880 | 0.851
0.5 | 0.836 | 0.887 | 0.853
1.0 | 0.836 | 0.886 | 0.854
2.0 | 0.820 | 0.872 | 0.839
3.0 | 0.787 | 0.832 | 0.805
10.0 | 0.596 | 0.627 | 0.612

###  4.4 Microstructure overall morphologies

Figure 7:  Spatiotemporal reconstruction of PF-simulated microstructures
during directional solidification of Mg-1wt%Gd at $V=166~{}$µm/s for different
temperature gradients $G$ (rows a to f) and crystal orientations $(\alpha_{\rm
A},\alpha_{\rm B})$ (columns I to IV). The total simulated length is typically
three to four times longer than that illustrated here.

Fig. 7 shows reconstructed microstructures (i.e. spatiotemporal snapshots)
resulting from typical PF simulations for different temperature gradients
(rows) from $G=0.3~{}$K/mm (a) to 10 K/mm (f) and different bi-crystal
orientations (columns), with the orientation of the red grain fixed at
$\alpha_{\rm A}=+10^{\circ}$, and that of the blue grain with $\alpha_{\rm
B}=-25^{\circ}$ (I), $-15^{\circ}$ (II), $+15^{\circ}$ (III), and
$+25^{\circ}$ (IV). Fig. 7 highlights the well-known difference between
converging and diverging grain boundaries, forming via different mechanisms of
dendrite impingement or side-branching competition, respectively [27, 43].

The selected orientation $\theta$ of converging GBs ($\alpha_{1}<\alpha_{2}$)
– via primary dendrite impingement – remains, in all cases, close to the
growth direction $\alpha$ of the most favorably oriented grain (here, the red
grain). However, some relatively rare events of well-oriented dendrite
elimination appear – see, e.g., three elimination events in the illustrated
example for $(\alpha_{\rm A},\alpha_{\rm B})=(+10^{\circ},-15^{\circ})$ at
$G=2~{}$K/mm (dII). These events, often referred to as “unusual overgrowth” in
the literature [50, 72], result from a mechanism already explained via PF
simulations in cubic systems, namely a progressive local reduction of spacing
(i.e. increase of undercooling) of the well-oriented dendrite immediately
adjacent to the GB as misoriented dendrites impinge onto it, until its
elimination occurs. These events result in small deviations of the GB
orientation from the favorably-oriented grain growth direction.

Like in cubic systems, the case of diverging GBs ($\alpha_{1}>\alpha_{2}$),
driven by the growth competition of side-branches, exhibits a greater
complexity and variability. When both grains grow in the same lateral
direction ($\alpha_{1}\times\alpha_{2}>0$), as illustrated in columns III-IV
of Fig. 7, the GB orientation remains relatively close to that of the
favorably-oriented grain. However, some deviation from it may result from rare
events of side-branching from the favorably oriented grain, in particular when
$\alpha_{1}$ and $\alpha_{2}$ are close at low $G$ (see, e.g., Fig. 7aIII).
When grains grow in opposite lateral directions
($\alpha_{1}\times\alpha_{2}<0$), emerging GBs exhibit a high lateral
mobility. As a result, the GBs may actually go in any direction, and their
shape may deviate significantly from straight lines (see, e.g., Fig. 7dII).
Because they grow at a $60^{\circ}$ angle from the primary dendrites, side-
branches always face a significant temperature gradient, compared to fourfold
(cubic) systems, where grains at $\alpha_{0}\approx 0^{\circ}$ generate side-
branches growing essentially along isotherms, in particular at low $G$ [27].
At low $G$, the resulting high side-branching activity contributes to the
increase in primary spacing, and hence to the available space between grains
for the side-branching competition to occur, thus also promoting the GB
lateral mobility and the variability of $\theta$ (see Section 4.5). At the
highest $G=10~{}$K/mm, microstructures clearly transition toward an
intermediate regime between dendritic and cellular, and a crystal orientation
$\alpha_{0}=\pm 25^{\circ}$ even leads to seaweed microstructures (Fig. 7fI
and 6fIV), not showing any clear primary or secondary branches.

The so-called rare events of elimination (at converging GBs) and of side-
branching (at diverging GBs) of well-oriented dendrites are fundamentally
different in nature. On the one hand, the “unusual overgrowth” at converging
GBs is directly related to a progressive increase of undercooling of the well-
oriented primary dendrite adjacent to the GB, as successive misoriented
dendrites impinge upon it at a relatively regular frequency. As such, it is
expected to have a frequency of occurrence that is directly related to the
primary spacings and primary growth directions in both grains, with a
relatively narrow statistical variability. Published results indicate that
this frequency tends to decrease when $\Lambda_{\alpha}$ increase or when
$|\alpha_{1}+\alpha_{2}|$ increase [43, 72]. On the other hand, unusual events
of side-branching from the well-oriented grain are expected to exhibit a much
broader level of statistical variability. This is due to the fact that
dendritic side-branching is, by essence, induced by selective amplification of
noise, and hence stochastic [73, 74, 75, 76]. However, in spite of the high
number of simulations presented here, the total number of unusual elimination
and branching events remains too limited for a significant analysis of their
frequency and statistical distribution. Further investigation remains needed
to conclusively establish whether these rare events could be rigorously
classified statistically.

At the scale of an overall grain structure (i.e. averaged over thousands,
millions, or more grain GBs), the effect of such small deviations on emerging
material properties is also likely small, in particular for converging grain
boundaries, which, for the most part, remain relatively close to straight
lines (Fig. 7). However, at the scale of individual GBs, given the intricate
relationship between GB structure and morphology with mechanical properties,
one may expect that deviations from a straight morphology may influence the
properties (e.g. increasing strength by disrupting the transmission of
dislocation or twins, or the extent of GB sliding).

###  4.5 GB orientation selection maps

Figure 8:  PF-predicted selected GB orientation $\theta$ as a function of the
crystal orientation angles ($\alpha_{1},\alpha_{2}$) for different values of
the temperature gradient $G$ for $V=166$ µm/s. Note that the the color bar
range $(-25^{\circ},+25^{\circ})$ is slightly narrower than in Fig. 1
$(-30^{\circ},+30^{\circ})$ due to the deviation of growth orientation from
crystal orientation (see Section 4.3) while the theoretical FOG and GL maps
(Fig. 1) consider $\alpha=\alpha_{0}$.  Figure 9:  Standard deviation from the
least-square fit of the PF-predicted GB trajectories ($y_{\rm GB}$) to a
straight line trajectory ($y_{\rm fit}$) normalized by the the domain width
($L_{y}$).  Figure 10:  Absolute deviation, averaged over all
$(\alpha_{1},\alpha_{2})$ orientations, of PF-predicted GB orientation,
$\theta_{\rm PF}$, compared to the theoretical FOG criterion,
$\widetilde{\theta}_{\rm FOG}$, and to the GL, $\widetilde{\theta}_{\rm GL}$,
both corrected to account for $\alpha/\alpha_{0}\neq 1$ (see text) and average
GB orientation of the overall PF map, $\overline{\theta_{\rm{PF}}}$, as a
function of the temperature gradient $G$. The color background represents the
transition between dendritic (orange) and intermediate dendritic-cellular
(purple) regimes, with illustrative microstructures shown at the top for
$G=0.3$ and 10 K/mm.

Fig. 8 shows the GB orientation $\theta$ as a function of
($\alpha_{1},\alpha_{2}$) obtained by phase-field simulations for
$G=[0.3-10]~{}$K/mm. Each data point is an average of ten simulations with
different fluctuations, except simulations for $\alpha_{1}=\alpha_{2}$
performed only once (see Table 3). Importantly, Fig. 8 data points along the
$(\alpha_{1}=-\alpha_{2})$ line, as well as ($\alpha_{1}=\alpha_{2}=\pm
30^{\circ}$) cases are set to $\theta=0^{\circ}$, due to symmetry
considerations. Indeed, for these simulations, a vertical mirror symmetry
leads to similar values of $(\alpha_{1},\alpha_{2})$ – yet interchanging the
roles of grains A and B as grain 1 or 2 – but an opposite value of $\theta$.
As such, considering $\theta\neq 0^{\circ}$ for these cases would actually
represent an arbitrary choice between these two opposite yet equivalent cases.

Fig. 9 shows the corresponding maps of the average standard deviation of the
linear fit of the $y-$position of the GB, $y_{\rm fit}$, fitted to the PF-
predicted position, $y_{\rm GB}$, normalized by the domain total height
$L_{y}$. These maps qualitatively show the nonlinearity (or roughness) of the
GB, which also illustrates the lateral ($y$) mobility of the GB as it is
forming.

Fig. 10 shows the deviation of the PF-calculated $\theta$ map over the entire
$(\alpha_{1},\alpha_{2})$ space compared to theoretical predictions using
$\alpha/\alpha_{0}$-corrected FOG (Fig. 1a) and GL (Fig. 1b) models, graphed
as a function of the temperature gradient $G$. Therein, FOG and GL grain
boundary orientations, respectively $\theta_{\rm FOG}$ and $\theta_{\rm GL}$,
were corrected by a multiplicative factor $\alpha/\alpha_{0}$ in order to
account for the deviation of dendrite growth direction ($\alpha$) from the
crystal orientation ($\alpha_{0}$). Since $\alpha/\alpha_{0}$ is significantly
dependent on $G$ only, we used a constant value for each $G$ (i.e. the mean
value listed in Table 4). The resulting corrected maps are denoted
$\widetilde{\theta}_{\mu}$, with $\mu\in\\{{\rm FOG},{\rm GL}\\}$, and the
deviation with phase-field results is quantified as $||\theta_{\rm
PF}-\widetilde{\theta}_{\mu}||$ defined as the square root of the average
value of $(\theta_{\rm PF}-\widetilde{\theta}_{\mu})^{2}$ over the
$(\alpha_{1},\alpha_{2})$ space. Fig. 10 also shows the PF-predicted value of
GB orientation averaged over the entire $(\alpha_{1},\alpha_{2})$ map
($\overline{\theta_{\rm PF}}$).

Both the GB orientation selection maps (Fig. 8) and the $||\theta_{\rm
PF}-\widetilde{\theta}_{\mu}||$ deviation plots (Fig. 10) exhibit two distinct
trends at low $G$ and at high $G$, which we attribute to a transition from a
well-developed dendritic regime ($G\leq 1$ K/mm) to an intermediate regime
between dendritic and cellular ($G\geq 2$ K/mm). This transition is
illustrated with the color gradient background of Fig. 10, where illustrated
microstructures for the two extreme values of $G$ clearly show the difference
between growth morphologies (also see Fig. 7).

Within the dendritic regime ($G\leq 1$ K/mm), the $\theta_{\rm PF}$ selection
map gets closer to the GL and further from the FOG criterion as $G$ decreases,
as the pattern becomes more pronouncedly dendritic, while the overall average
GB orientation does not really exhibits any significant variation
($9^{\circ}<\overline{\theta_{\rm PF}}<10^{\circ}$). In Fig. 8, most of the
$\theta(\alpha_{1},\alpha_{2})$ map exhibits little variation within this $G$
range, but the most striking change appears for diverging GBs with
$\alpha_{1}\times\alpha_{2}<0$, i.e. the bottom-right quadrant of maps. This
transition from FOG to GL as $G$ decreases in the dendritic regime is
consistent with prior observation from PF simulation of fourfold symmetry
systems [45]. We also confirmed it further via simulations at lower $G$ for
selected $(\alpha_{1},\alpha_{2})$ configurations.

Figure 11:  EBSD inverse pole figure maps obtained from two different regions
of a Mg-1wt%Gd sample directionally solidified at a cooling rate
$R=30~{}$K/min and $G=3$ K/mm (i.e. $V\approx 166~{}$µm/s). Grains with
relatively low $c-$axis misorientation are labelled A through G and shown next
to arrows illustrating the corresponding crystal orientation $\alpha_{0}$, and
resulting grain boundaries are labelled I through V (see Table 5).  Table 5:
Comparison of experimental GB orientations measured in Fig. 11 ($\theta_{\rm
Exp}$) with PF-predicted ones ($\theta_{\rm PF}$) interpolated from
$\theta(\alpha_{1},\alpha_{2})$ of Fig. 8e, for selected grain boundaries with
relatively low $c-$axis misorientation ($\Phi$) with respect to the normal to
the sample/imaging plane.
GB | Grains 1/2 |  $\alpha_{1}$(∘) |  $\alpha_{2}$(∘) |  $\theta_{\rm Exp}$(∘) |  $\theta_{\rm PF}$(∘) |  $\Phi_{1}$(∘) |  $\Phi_{2}$(∘) | Type | $\alpha_{1}\times\alpha_{2}$
---|---|---|---|---|---|---|---|---|---
I | A/B | 8.1 | 5.6 | 2 | 4.5 | 15.4 | 21.8 | Diverging | $>0$
II | A/C | 8.1 | $-9.5$ | $-4$ | 0.9 | 15.4 | 5.8 | Diverging | $<0$
III | C/D | $-9.5$ | $-13.6$ | $-12$ | $-8.1$ | 5.8 | 14.5 | Diverging | $>0$
IV | E/F | 1.5 | $-29.8$ | $-1.6$ | 1.1 | 20.6 | 12.3 | Diverging | $<0$
V | F/G | $-29.8$ | 7.2 | 16 | 6.2 | 12.3 | 15.5 | Converging | $<0$

For high temperature gradients ($G\geq 2~{}$K/mm), the microstructure
experiences a transition between dendritic (low $G$) and cellular (high $G$)
regimes. The deviation against FOG and GL (Fig. 10) shows a reverse trend
compared to the dendritic regime, but rather than strictly approaching the GL
at high $G$, this change is principally caused by the significant decrease in
$\overline{\theta_{\rm PF}}$ with $G$. Indeed, as $G$ increases and the
microstructure approaches the cellular regime, not only does the growth
directions $\alpha_{1}$ and $\alpha_{2}$ tend toward the temperature gradient
direction, but so does the GB orientation. This is apparent in Fig. 8 by the
significant reduction of $|\theta|$ throughout the entire map at high $G$. We
did not explore temperature gradients much higher than 10 K/mm because, as the
cells get shallower and the interface approaches a nearly-planar interface,
the dynamics of triple points, not properly accounted for in the current
single-phase-field model, starts taking a prominent role in the selection of
the GB orientation [77].

In terms of GB orientation variability and lateral mobility (Fig. 9), the
highest deviations appear systematically for diverging GBs (bottom-right
half), in particular for $(\alpha_{1}\times\alpha_{2}<0)$ configurations
(bottom-right quadrant), as well as at (or close to) symmetric configurations,
i.e. for $(\alpha_{1}\approx-\alpha_{2})$. High variability also appears close
to the $\alpha_{i}=\pm 30^{\circ}$ corners ($i\in\\{1,2\\}$), in relation to
the microstructure transition toward a seaweed pattern. As mentioned in Sec.
4.4 (Fig. 7), the side-branching activity along diverging GBs is promoted by a
low $G$. This results in a broader region in the $(\alpha_{1},\alpha_{2})$
space that exhibits a significant deviation, compared to higher $G$ where the
high variability is concentrated mostly on symmetric
$(\alpha_{1}=-\alpha_{2})$ configurations. Interestingly, at both low and high
$G$, a region of high variability also emerges along the
$(\alpha_{1}=-\alpha_{2})$ line for converging GBs (top-left half).

###  4.6 Comparison with experiments

Fig. 11 shows EBSD maps obtained from two different regions of a Mg-1wt$\%$Gd
sample after a directional solidification at a cooling rate $R=30~{}$K/min and
$G=3~{}$K/mm (i.e. $V=166$ µm/s). These regions were selected because they
contain GBs between grains whose $c$ axis deviation from the $z$ axis
($\Phi_{i}$ in Table 5) is lower than $25^{\circ}$. These grains are labelled
A through G, and the five resulting GBs are labelled I through V. Table 5
summarizes the grain orientations $(\alpha_{1},\alpha_{2})$ within the $(x,y)$
plane, the tilt angle $\Phi$ between the $c$-axis of both grains and the $z$
coordinate axis, as well as measured GB orientation, $\theta_{\rm Exp}$,
compared to the PF-predicted orientation, $\theta_{\rm PF}$. The latter is
calculated by bilinear interpolation from the $\theta(\alpha_{1},\alpha_{2})$
map of Fig. 8e.

The comparison between simulation and experiments, summarized in Table 5, show
some qualitative agreement, but also highlights some notable discrepancies.
While three GB (I, III, V) exhibit orientations of the same sign, two actually
have opposite signs (II, IV). The latter correspond to diverging GBs with
$\alpha_{1}\times\alpha_{2}<0$, which we have shown to be the most prone to
high deviations. These GB angles also have at most $4^{\circ}$ in absolute
value, making the change of sign within the expected margin of stochastic
variability. Among the other three GBs, the two orientations closest to PF
predictions exhibit small deviations of $2.5^{\circ}$ (I) and $3.9^{\circ}$
(III), with both being diverging GBs with $\alpha_{1}\times\alpha_{2}>0$,
which are expected to show relatively low variability (Sec. 4.5). The last GB
(V), showing a discrepancy of almost $10^{\circ}$, is a converging GB with
$\alpha_{1}\times\alpha_{2}<0$ and, importantly, it involves a grain (F) with
an absolute misorientation close to $\pm 30^{\circ}$, putting it within the
region of highest variability among converging GBs (see Fig. 9e). Another
factor that could contribute to the discrepancy between simulation and
experiments is the relatively small observable length of the GB, in particular
for GBs II, IV, and V (Fig. 11).

In addition to the limited statistical sample size, another very important
source of possible discrepancy between simulation and experiments is the
inherent comparison of three-dimensional experiments with two-dimensional
simulations. Even though we filtered out the grains with a basal plane most
titled from the observation plane to focus on GBs with both
$\Phi_{i}<25^{\circ}$, the tilt angles of analyzed grains ($\Phi_{i}$ in Table
5) are not negligible (from 5.8 to $21.8^{\circ}$). Three-dimensional grain
growth competition, in addition to resulting in additional degrees of freedom
to describe the 3D bi-crystalline orientations, have been shown to result in
unexpectedly complex mechanisms [52, 28]. These not only include GB
orientation selection, but also the morphological destabilization of the GB
that cannot be represented any longer by a straight line (or plane in 3D), but
instead can also adopt highly convoluted shapes [28]. The sample thickness (1
mm) may also allow several grains to develop within the thickness, further
contributing to the 3D growth competition. This is prominently illustrated in
Fig. 11 by grain A growing in two separate locations of the sample.

In conclusion, even though the analysis of GB types allows a reasonable
qualitative understanding of the discrepancy between simulation and
experiments, we reckon that further work is required to obtain a fully
conclusive comparison between simulation and experiments. Efforts in this
direction are suggested in the conclusion section, most of which serve as the
basis for ongoing work.

##  5 Summary and perspectives

We presented an extensive numerical study of grain growth competition during
directional solidification of hcp crystals, considering a Mg-1wt$\%$Gd alloy,
and studied the effect of the temperature gradient over a wide range of $G$.
We also compared our results to classical and state-of-the-art theories for
primary spacing selection and preferred growth direction. We have shown that
these theories, and the overall understanding of GB orientation selection in
cubic (fcc,bcc) crystals may essentially be extended to the case of hcp
crystals in the case of (nearly) two-dimensional configurations (i.e. within
the basal plane), sometimes requiring scaling factors to account for the
sixfold (rather than fourfold) symmetry.

Regarding grain growth competition, we have shown that, like in cubic
crystals, converging grain boundaries and diverging grain boundaries with
$\alpha_{1}\times\alpha_{2}>0$ (see Fig. 2) tend to follow the orientation of
the most favorably oriented grain, however with some amount of deviation due
to relatively rare events of unusual overgrowth (at converging GBs) or side-
branching (at diverging GBs) of favorably oriented dendrites. The latter may
actually be promoted by the fact that secondary branches, growing at a
$60^{\circ}$ from the primary branch, actually face a higher temperature
gradient than in fourfold (cubic) systems. GB configurations exhibiting the
highest stochastic variability are diverging GBs with
$\alpha_{1}\times\alpha_{2}<0$, converging GBs close to a symmetric
configuration with $\alpha_{1}=-\alpha_{2}$, and growth competition involving
at least one grain with $\alpha_{0}\approx\pm 30^{\circ}$, resulting in a
dendrite-to-seaweed pattern transition.

Beyond these trends that could essentially be extrapolated from current
knowledge on cubic crystals, we have also studied the effect of the
temperature gradient $G$ on grain growth competition. We found that, within
the dendritic regime, a decrease of $G$ leads to the outcome of the growth
competition approaching the theoretical geometrical limit (GL, Fig 1b),
whereas when increasing $G$ the selection map tends toward the favorably-
oriented grain criterion (FOG, Fig 1a), as suggested by previous studies on a
subset of $(\alpha_{1},\alpha_{2})$ configuration in a cubic system [45]. When
increasing $G$ further (here for $G\geq 2~{}$K/mm), the microstructure
exhibits a transition from a dendritic to a cellular morphology. Then, as $G$
increases, the cell/dendrite growth direction tends to the temperature
gradient direction ($\alpha\to 0^{\circ}$) and, consequently, the selected GB
orientation also significantly decreases to align close to that of $G$.

Finally, we also compared our PF prediction with experimental measurements on
a directionally solidified thin samples of Mg-1wt$\%$Gd alloy. Within the few
GB analyzed in details, we found a limited qualitative agreement, with notable
discrepancies correlated with GB types expected to show the most stochastic
variability according to our PF results. Still, while we provide explanations
for the discrepancy between experimental and computed GB angles (in particular
three-dimensionality, stochasticity with small statistical sample size, and GB
types), a fully conclusive and quantitative experimental comparison requires
further research in the following directions.

Ongoing research includes gathering a broader statistic of GB orientation
measurements, in thinner ($\approx 200~{}$µm) samples, which will also be
imaged using X-ray in situ radiography and 3D reconstructed tomography of the
solidified samples. Another important perspective could involve the use of
fully 3D simulations. However, while studies of individually selected
orientation configurations are achievable [52, 28], a full mapping of the 3D
bicrystal orientation space, including multiple runs to account for GB
stochasticity, seems out of reach even using the most high-performance PF
implementations to date [52, 78]. Emerging methods to explore the widely
multidimensional space of alloys, processing parameters, and bi-crystal
orientations, could include data-based tools, such as machine learning.
Regardless of whether further understanding will come from PF exploration
similar to the present one or from data-based methods, a wider statistical
dataset of experimental GB orientation under well-defined (e.g. $G$, $V$)
conditions must be gathered. While this is not a conceptually challenging
task, it could be tedious to obtain, and probably for this reason has not been
undertaken, to our knowledge, by any group, for any crystal symmetry.

In summary, we have shown that several of the behaviors observed in grain
growth competition of cubic crystals can be extrapolated to hcp crystals
within their basal plane, provided the appropriate rescaling and re-evaluation
of some fitting parameters. We have also confirmed the transition from a GB
selection from FOG to GL behavior with a decrease of temperature gradient,
previously suggested for a limited set orientation configurations [45], as
well as the overall reduction of GB misorientation as $G$ increases beyond the
onset of the dendritic-to-cellular morphological transition regime.
Experimental comparison, while encouraging, have highlighted the need for a
broad statistically-relevant data set of GB orientation measurements from
well-identified solidification conditions, as well as three-dimensional
simulations, potentially supported by data-based approaches, in order to
establish a more comprehensive understanding of grain growth competition
during alloy solidification.

## Acknowledgements

This research was supported by the Spanish Ministry of Science within the
scope of the HexaGB project (grant number RTI2018-098245) and the project
MAD2D-CM of the Comunidad de Madrid, Recovery, Transformation and Resilience
Plan, and NextGenerationEU from the European Union. DT also gratefully
acknowledges support from the Spanish Ministry of Science through a Ramón y
Cajal Fellowship (RYC2019-028233-I).

##  Appendix A Anisotropy vector formalism for HCP symmetries

In the present work, as in [61, 79, 17], the evolution equation of the phase-
field is written in terms of an auxiliary anisotropy vector $\vec{A}$, where

$\vec{\nabla}\cdot\vec{A}=\sum_{\alpha}\partial_{\alpha}\Big{[}|\vec{\nabla}\varphi|^{2}a_{s}\Big{(}\frac{\partial
a_{s}}{\partial\varphi_{\alpha}}\Big{)}\Big{]},$ (16)

where $\alpha\in\\{x,y,z\\}$, and the components of vector
$\vec{A}=(A_{x},A_{y},A_{z})$ directly depend on the anisotropy function of
the surface energy, $a_{s}$, and can be computed for any arbitrary function
$a_{s}(n_{x},n_{y},n_{z})$.
Using the preconditioned phase-field

$\psi=\sqrt{2}\tanh^{-1}({\varphi})$ (17)

and its first derivatives

$\partial_{\alpha}\varphi=\frac{1}{\sqrt{2}}(1-\varphi^{2})\partial_{\alpha}\psi,$
(18)

one obtains

$A_{\alpha}=\frac{1-\varphi^{2}}{\sqrt{2}}|\vec{\nabla}\psi|^{2}a_{s}\frac{\partial
a_{s}}{\partial\psi_{\alpha}}.$ (19)

Introducing the components

$n_{\alpha}=-\frac{\partial_{\alpha}\psi}{|\vec{\nabla}\psi|}$ (20)

of the normal vector $\vec{n}$ to the solid-liquid interface, the last term in
Eq. (19) becomes

$\frac{\partial a_{s}}{\partial\psi_{\alpha}}=\sum_{\beta}\frac{\partial
a_{s}}{\partial n_{\beta}}\ \frac{\partial
n_{\beta}}{\partial\psi_{\alpha}}=\frac{1}{|\vec{\nabla}\psi|}\sum_{\beta}(n_{\alpha}n_{\beta}-\delta_{\alpha,\beta})\frac{\partial
a_{s}}{\partial n_{\beta}}$ (21)

with $\beta\in\\{x,y,z\\}$, and $\delta_{\alpha,\beta}=1$ if $\alpha=\beta$ or
$0$ otherwise. This leads to

$A_{\alpha}=\frac{(1-\varphi^{2})}{\sqrt{2}}|\vec{\nabla}\psi|a_{s}\sum_{\beta}(n_{\alpha}n_{\beta}-\delta_{\alpha,\beta})\frac{\partial
a_{s}}{\partial n_{\beta}}.$ (22)

### Case of a HCP symmetry

For HCP crystals, the anisotropy function [54] is

$\displaystyle\gamma(\vec{n})=~{}\gamma_{0}a_{s}(\vec{n})=~{}$
$\displaystyle\gamma_{0}\Big{[}1+\varepsilon_{20}y_{20}(\vec{n})+\varepsilon_{40}y_{40}(\vec{n})$
$\displaystyle+\varepsilon_{60}y_{60}(\vec{n})+\varepsilon_{66}y_{66}(\vec{n})\Big{]},$
(23)

where

$\displaystyle y_{20}$ $\displaystyle=\sqrt{\frac{5}{16\pi}}(3n_{z}^{2}-1)$
$\displaystyle=a_{20}(3n_{z}^{2}-1);$ (24) $\displaystyle y_{40}$
$\displaystyle=\frac{3}{16\sqrt{\pi}}(35n_{z}^{4}-30n_{z}^{2}+3)$
$\displaystyle=a_{40}(35n_{z}^{4}-30n_{z}^{2}+3);$ (25) $\displaystyle y_{60}$
$\displaystyle=\frac{\sqrt{13}}{32\sqrt{\pi}}(231nz^{6}-315n_{z}^{4}+105n_{z}^{2}-5)$
$\displaystyle=a_{60}(231n_{z}^{6}-315n_{z}^{4}+105n_{z}^{2}-5);$
$\displaystyle y_{66}$
$\displaystyle=\frac{\sqrt{6006}}{64\sqrt{\pi}}(n_{x}^{6}-15n_{x}^{4}n_{y}^{2}+15n_{x}^{2}n_{y}^{4}-n_{y}^{6})$
$\displaystyle=a_{66}(n_{x}^{6}-15n_{x}^{4}n_{y}^{2}+15n_{x}^{2}n_{y}^{4}-n_{y}^{6}).$

Then

$\displaystyle\frac{\partial a_{s}}{\partial n_{x}}$
$\displaystyle=6\varepsilon_{66}a_{66}n_{x}(n_{x}^{4}-10n_{x}^{2}n_{y}^{2}+5ny^{4});$
(28) $\displaystyle\frac{\partial a_{s}}{\partial n_{y}}$
$\displaystyle=-6\varepsilon_{66}a_{66}n_{y}(5n_{x}^{4}-10n_{x}^{2}n_{y}^{2}+n_{y}^{4});$
(29) $\displaystyle\frac{\partial a_{s}}{\partial n_{z}}$
$\displaystyle=n_{z}\Big{[}6\varepsilon_{20}a_{20}+20\varepsilon_{40}a_{40}(7n_{z}^{2}-3)$
$\displaystyle+6\varepsilon_{60}a_{60}(231n_{z}^{4}-210n_{z}^{2}+35)\Big{]}.$
(30)

This leads to the components of $\vec{A}$

$\displaystyle A_{x}=~{}$
$\displaystyle\frac{1-\varphi^{2}}{\sqrt{2}}|\vec{\nabla}\psi|a_{s}~{}\times$
$\displaystyle\Big{\\{}6\varepsilon_{66}a_{66}n_{x}(n_{x}^{4}-10n_{x}^{2}n_{y}^{2}+5n_{y}^{4})(n_{x}^{2}-1)$
$\displaystyle-6\varepsilon_{66}a_{66}n_{y}^{2}n_{x}(5n_{x}^{4}+10n_{x}^{2}n_{y}^{2}-n_{y}^{4})$
$\displaystyle+n_{z}^{2}n_{x}\Big{[}6\varepsilon_{20}a_{20}+20\varepsilon_{40}a_{40}(7n_{z}^{2}-3)$
$\displaystyle+2\varepsilon_{60}a_{60}(693n_{z}^{4}-630n_{z}^{2}+105n_{z})\Big{]}\Big{\\}};$
(31) $\displaystyle A_{y}=~{}$
$\displaystyle\frac{1-\varphi^{2}}{\sqrt{2}}|\vec{\nabla}\psi|a_{s}~{}\times$
$\displaystyle\Big{\\{}6\varepsilon_{66}a_{66}n_{y}\Big{[}n_{x}^{2}(n_{x}^{4}-10n_{x}^{2}n_{y}^{2}+5n_{y}^{4})$
$\displaystyle-(n_{y}^{2}-1)(5n_{x}^{4}-10n_{x}^{2}n_{y}^{2}+n_{y}^{4})\Big{]}$
$\displaystyle+2n_{z}^{2}n_{y}\Big{[}3\varepsilon_{20}a_{20}+10\varepsilon_{40}a_{40}(7n_{z}^{2}-3)$
$\displaystyle+\varepsilon_{60}a_{60}(693n_{z}^{4}-630n_{z}^{2}+105)\Big{]}\Big{\\}};$
(32)

$\displaystyle A_{z}=~{}$
$\displaystyle\frac{1-\varphi^{2}}{\sqrt{2}}|\vec{\nabla}\psi|a_{s}~{}\times$
$\displaystyle\Big{\\{}6\varepsilon_{66}a_{66}n_{z}\Big{[}n_{x}^{2}(n_{x}^{4}-10n_{x}^{2}n_{y}^{2}+5n_{y}^{4})$
$\displaystyle-n_{y}^{2}(5n_{x}^{4}-10n_{x}^{2}n_{y}^{2}+n_{y}^{4})\Big{]}$
$\displaystyle+2n_{z}(n_{z}^{2}-1)\Big{[}3\varepsilon_{20}a_{20}+10\varepsilon_{40}a_{40}(7n_{z}^{2}-3)$
$\displaystyle+\varepsilon_{60}a_{60}(693n_{z}^{4}-630n_{z}^{2}+105)\Big{]}\Big{\\}}.$
(33)

For convenience, these components can be reformulated as:

$\displaystyle A_{x}=~{}$
$\displaystyle\frac{1-\varphi^{2}}{\sqrt{2}}|\vec{\nabla}\psi|a_{s}~{}\times$
$\displaystyle
n_{x}\Big{\\{}-a_{66}\varepsilon_{66}n_{y}^{2}(30n_{x}^{4}-60n_{x}^{2}n_{y}^{2}+6n_{y}^{4})$
$\displaystyle+a_{66}\varepsilon_{66}(n_{x}^{2}-1)(6n_{x}^{4}-60n_{x}^{2}n_{y}^{2}+30n_{y}^{4})$
$\displaystyle+n_{z}^{2}[6a_{20}\varepsilon_{20}+a_{40}\varepsilon_{40}(140n_{z}^{2}-60)$
$\displaystyle+a_{60}\varepsilon_{60}(1386n_{z}^{4}-1260n_{z}^{2}+210)]\Big{\\}}$
(34) $\displaystyle A_{y}=~{}$
$\displaystyle\frac{1-\varphi^{2}}{\sqrt{2}}|\vec{\nabla}\psi|a_{s}~{}\times$
$\displaystyle
n_{y}\Big{\\{}a_{66}\varepsilon_{66}n_{x}^{2}(6n_{x}^{4}-60n_{x}^{2}n_{y}^{2}+30n_{y}^{4})$
$\displaystyle-
a_{66}\varepsilon_{66}(n_{y}^{2}-1)(30n_{x}^{4}-60n_{x}^{2}n_{y}^{2}+6n_{y}^{4})$
$\displaystyle+n_{z}^{2}[6a_{20}\varepsilon_{20}+a_{40}\varepsilon_{40}(140n_{z}^{2}-60)$
$\displaystyle+a_{60}\varepsilon_{60}(1386n_{z}^{4}-1260n_{z}^{2}+210)]\Big{\\}}$
(35) $\displaystyle A_{z}=~{}$
$\displaystyle\frac{1-\varphi^{2}}{\sqrt{2}}|\vec{\nabla}\psi|a_{s}~{}\times$
$\displaystyle
n_{z}\Big{\\{}a_{66}\varepsilon_{66}n_{x}^{2}(6n_{x}^{4}-60n_{x}^{2}n_{y}^{2}+30n_{y}^{4})$
$\displaystyle-
a_{66}\varepsilon_{66}n_{y}^{2}(30n_{x}^{4}-60n_{x}^{2}n_{y}^{2}+6n_{y}^{4})$
$\displaystyle+(n_{z}^{2}-1)[6a_{20}\varepsilon_{20}+a_{40}\varepsilon_{40}(140n_{z}^{2}-60)$
$\displaystyle+a_{60}\varepsilon_{60}(1386n_{z}^{4}-1260n_{z}^{2}+210)]\Big{\\}}$
(36)

According to [58],

$\varepsilon_{20}$ | $\varepsilon_{40}$ | $\varepsilon_{60}$ | $\varepsilon_{66}$
---|---|---|---
$-0.026$ | 0 | 0 | 0.003

In two dimensions, the projection of Eq. (23) on the xy-plane ($n_{z}=0$)
leads to

$\displaystyle\gamma(\vec{n})=\gamma_{0}a_{s}(\vec{n})=\gamma_{0}\Big{[}1.008+\varepsilon_{66}y_{66}(\vec{n})\Big{]}$
(37)

and the components of $\vec{A}$, as they appear in the code, are

$\displaystyle A_{x}=~{}$
$\displaystyle\frac{1-\varphi^{2}}{\sqrt{2}}|\vec{\nabla}\psi|a_{s}a_{66}\varepsilon_{66}~{}\times$
$\displaystyle
n_{x}\Big{\\{}-n_{y}^{2}(30n_{x}^{4}-60n_{x}^{2}n_{y}^{2}+6n_{y}^{4})$
$\displaystyle+(n_{x}^{2}-1)(6n_{x}^{4}-60n_{x}^{2}n_{y}^{2}+30n_{y}^{4})\Big{\\}}$
(38) $\displaystyle A_{y}=~{}$
$\displaystyle\frac{1-\varphi^{2}}{\sqrt{2}}|\vec{\nabla}\psi|a_{s}a_{66}\varepsilon_{66}~{}\times$
$\displaystyle
n_{y}\Big{\\{}n_{x}^{2}(6n_{x}^{4}-60n_{x}^{2}n_{y}^{2}+30n_{y}^{4})$
$\displaystyle-(n_{y}^{2}-1)(30n_{x}^{4}-60n_{x}^{2}n_{y}^{2}+6n_{y}^{4})\Big{\\}}$
(39)

## References

  * [1] W. Kurz, D. J. Fisher, Fundamentals of solidification, CRC Press, 1984.
  * [2] J. A. Dantzig, M. Rappaz, Solidification: -Revised & Expanded, EPFL press, 2016\.
  * [3] W. Kurz, D. J. Fisher, R. Trivedi, Progress in modelling solidification microstructures in metals and alloys: dendrites and cells from 1700 to 2000, International Materials Reviews 64 (6) (2019) 311–354.
  * [4] W. Kurz, M. Rappaz, R. Trivedi, Progress in modelling solidification microstructures in metals and alloys. part ii: dendrites from 2001 to 2018, International Materials Reviews 66 (1) (2021) 30–76.
  * [5] B. Mordike, T. Ebert, Magnesium: properties—applications—potential, Materials Science and Engineering: A 302 (1) (2001) 37–45.
  * [6] A. I. Taub, P. E. Krajewski, A. A. Luo, J. N. Owens, The evolution of technology for materials processing over the last 50 years: the automotive example, Jom 59 (2007) 48–57.
  * [7] A. A. Luo, Magnesium casting technology for structural applications, journal of Magnesium and Alloys 1 (1) (2013) 2–22.
  * [8] N. Hort, Y.-d. Huang, D. Fechner, M. Störmer, C. Blawert, F. Witte, C. Vogt, H. Drücker, R. Willumeit, K. Kainer, et al., Magnesium alloys as implant materials–principles of property design for mg–re alloys, Acta biomaterialia 6 (5) (2010) 1714–1725.
  * [9] Y. Chen, Z. Xu, C. Smith, J. Sankar, Recent advances on the development of magnesium alloys for biodegradable implants, Acta biomaterialia 10 (11) (2014) 4561–4573.
  * [10] M. Yang, F. Pan, Effects of y addition on as-cast microstructure and mechanical properties of mg–3sn–2ca (wt.%) magnesium alloy, Materials Science and Engineering: A 525 (1) (2009) 112–120.  doi:10.1016/j.msea.2009.06.040.
  * [11] Q. Peng, Y. Huang, K. U. Kainer, N. Hort, High ductile as-cast mg–re based alloys at room temperature, Materials Letters 83 (2012) 209–212.  doi:10.1016/j.matlet.2011.08.009.
  * [12] Y. Xu, F. Gensch, Z. Ren, K. U. Kainer, N. Hort, Effects of gd solutes on hardness and yield strength of mg alloys, Progress in Natural Science: Materials International 28 (6) (2018) 724–730.  doi:10.1016/j.pnsc.2018.10.002.
  * [13] W. J. Boettinger, S. R. Coriell, A. Greer, A. Karma, W. Kurz, M. Rappaz, R. Trivedi, Solidification microstructures: recent developments, future directions, Acta materialia 48 (1) (2000) 43–70.
  * [14] M. Asta, C. Beckermann, A. Karma, W. Kurz, R. Napolitano, M. Plapp, G. Purdy, M. Rappaz, R. Trivedi, Solidification microstructures and solid-state parallels: Recent developments, future directions, Acta Materialia 57 (4) (2009) 941–971.
  * [15] A. Karma, Phase-field formulation for quantitative modeling of alloy solidification, Physical review letters 87 (11) (2001) 115701.  doi:10.1103/PhysRevLett.87.115701.
  * [16] B. Echebarria, R. Folch, A. Karma, M. Plapp, Quantitative phase-field model of alloy solidification, Physical review E 70 (6) (2004) 061604.  doi:10.1103/PhysRevLett.87.115701.
  * [17] A. K. Boukellal, M. Rouby, J.-M. Debierre, Tip dynamics for equiaxed al-cu dendrites in thin samples: Phase-field study of thermodynamic effects, Computational Materials Science 186 (2021) 110051.  doi:10.1016/j.commatsci.2020.110051.
  * [18] W. Kurz, D. Fisher, Dendrite growth at the limit of stability: tip radius and spacing, Acta Metallurgica 29 (1) (1981) 11–20.  doi:10.1016/0001-6160(81)90082-1.
  * [19] K. Somboonsuk, R. Trivedi, Dynamical studies of dendritic growth, Acta Metallurgica 33 (6) (1985) 1051–1060.  doi:10.1016/0001-6160(85)90198-1.
  * [20] R. Trivedi, K. Somboonsuk, Pattern formation during the directional solidification of binary systems, Acta metallurgica 33 (6) (1985) 1061–1068.  doi:10.1016/0001-6160(85)90199-3.
  * [21] S.-Z. Lu, J. Hunt, A numerical analysis of dendritic and cellular array growth: the spacing adjustment mechanisms, Journal of Crystal Growth 123 (1-2) (1992) 17–34.  doi:10.1016/0022-0248(92)90006-5.
  * [22] J. A. Warren, J. Langer, Prediction of dendritic spacings in a directional-solidification experiment, Physical Review E 47 (4) (1993) 2702.  doi:10.1103/PhysRevE.47.2702.
  * [23] S. H. Han, R. Trivedi, Primary spacing selection in directionally solidified alloys, Acta Metallurgica Et Materialia 42 (1) (1994) 25–41.  doi:10.1016/0956-7151(94)90045-0.
  * [24] J. Hunt, S.-Z. Lu, Numerical modeling of cellular/dendritic array growth: spacing and structure predictions, Metallurgical and Materials Transactions A 27 (3) (1996) 611–623.  doi:10.1007/BF02648950.
  * [25] R. Trivedi, W. Kurz, Solidification microstructures: A conceptual approach, Acta Metallurgica et Materialia 42 (1) (1994) 15–23.  doi:10.1016/0956-7151(94)90044-2.
  * [26] C.-A. Gandin, M. Eshelman, R. Trivedi, Orientation dependence of primary dendrite spacing, Metallurgical and Materials Transactions A 27 (9) (1996) 2727\.  doi:10.1007/BF02652367.
  * [27] D. Tourret, A. Karma, Growth competition of columnar dendritic grains: A phase-field study, Acta Materialia 82 (2015) 64–83.  doi:10.1016/j.actamat.2014.08.049.
  * [28] Y. Song, F. L. Mota, D. Tourret, K. Ji, B. Billia, R. Trivedi, N. Bergeon, A. Karma, Cell invasion during competitive growth of polycrystalline solidification patterns, Nature Communications 14 (1) (2023) 2244.
  * [29] T. Haxhimali, A. Karma, F. Gonzales, M. Rappaz, Orientation selection in dendritic evolution, Nature materials 5 (8) (2006) 660–664.
  * [30] J. Dantzig, P. Di Napoli, J. Friedli, M. Rappaz, Dendritic growth morphologies in al-zn alloys—part ii: Phase-field computations, Metallurgical and Materials Transactions A 44 (2013) 5532–5543.
  * [31] M. Becker, J. Dantzig, M. Kolbe, S. T. Wiese, F. Kargl, Dendrite orientation transition in alge alloys, Acta Materialia 165 (2019) 666–677.
  * [32] L. Wang, J. J. Hoyt, N. Wang, N. Provatas, C. W. Sinclair, Controlling solid-liquid interfacial energy anisotropy through the isotropic liquid, Nature communications 11 (1) (2020) 724.
  * [33] M. Bedel, G. Reinhart, A.-A. Bogno, C.-A. Gandin, S. Jacomet, E. Boller, H. Nguyen-Thi, H. Henein, Characterization of dendrite morphologies in rapidly solidified al–4.5 wt.% cu droplets, Acta Materialia 89 (2015) 234–246.
  * [34] J. Strutzenberger, J. Faderl, Solidification and spangle formation of hot-dip-galvanized zinc coatings, Metallurgical and Materials Transactions A 29 (1998) 631–646.
  * [35] A. Sémoroz, L. Strezov, M. Rappaz, Orientation domains and texture in hot-dipped galvanized coatings, Metallurgical and Materials Transactions A 33 (2002) 2695–2701.
  * [36] K. Pettersen, N. Ryum, Crystallography of directionally solidified magnesium alloy az91, Metallurgical Transactions A 20 (5) (1989) 847–852.
  * [37] K. Pettersen, O. Lohne, N. Ryum, Dendritic solidification of magnesium alloy az91, Metallurgical Transactions A 21 (1990) 221–230.
  * [38] S. Akamatsu, T. Ihle, Similarity law for the tilt angle of dendrites in directional solidification of non-axially-oriented crystals, Phys. Rev. E 56 (4) (1997) 4479.  doi:10.1103/PhysRevE.56.4479.
  * [39] J. Deschamps, M. Georgelin, A. Pocheau, Growth directions of microstructures in directional solidification of crystalline materials, Physical Review E 78 (1) (2008) 011605.  doi:10.1103/PhysRevE.78.011605.
  * [40] Y. Song, S. Akamatsu, S. Bottin-Rousseau, A. Karma, Propagative selection of tilted array patterns in directional solidification, Physical Review Materials 2 (5) (2018) 053403.
  * [41] F. Mota, J. Pereda, K. Ji, Y. Song, R. Trivedi, A. Karma, N. Bergeon, Effect of sub-boundaries on primary spacing dynamics during 3d directional solidification conducted on declic-dsi, Acta Materialia 204 (2021) 116500.
  * [42] D. Walton, B. Chalmers, The origin of the preferred orientation in the columnar zone of ingots, Transactions of the American Institute of Mining and Metallurgical Engineers 215 (3) (1959) 447–457.
  * [43] D. Tourret, Y. Song, A. J. Clarke, A. Karma, Grain growth competition during thin-sample directional solidification of dendritic microstructures: A phase-field study, Acta Materialia 122 (2017) 220–235.  doi:10.1016/j.actamat.2016.09.055.
  * [44] A. Pineau, G. Guillemot, D. Tourret, A. Karma, C.-A. Gandin, Growth competition between columnar dendritic grains–cellular automaton versus phase field modeling, Acta Materialia 155 (2018) 286–301.  doi:10.1016/j.actamat.2018.05.032.
  * [45] E. Dorari, K. Ji, G. Guillemot, C.-A. Gandin, A. Karma, Growth competition between columnar dendritic grains–the role of microstructural length scales, Acta Materialia 223 (2022) 117395.  doi:10.1016/j.actamat.2021.117395.
  * [46] H. Esaka, J. Stramke, W. Kurz, Columnar dendrite growth in succinonitrile-acetone alloys, Film, Ecole Polytechnique Fédérale de Lausanne, LSMX, Lausanne, Switzerland (1985).
  * [47] H. Esaka, Dendrite growth and spacing in succinonitrile-acetone alloys, Tech. rep., EPFL (1986).  doi:10.5075/epfl-thesis-615.
  * [48] A. Wagner, B. Shollock, M. McLean, Grain structure development in directional solidification of nickel-base superalloys, Materials Science and Engineering: A 374 (1-2) (2004) 270–279.  doi:10.1016/j.msea.2004.03.017.
  * [49] Y. Zhou, A. Volek, N. Green, Mechanism of competitive grain growth in directional solidification of a nickel-base superalloy, Acta Materialia 56 (11) (2008) 2631–2637.  doi:10.1016/j.actamat.2008.02.022.
  * [50] J. Li, Z. Wang, Y. Wang, J. Wang, Phase-field study of competitive dendritic growth of converging grains during directional solidification, Acta Materialia 60 (4) (2012) 1478–1493.  doi:10.1016/j.actamat.2011.11.037.
  * [51] T. Takaki, M. Ohno, Y. Shibuta, S. Sakane, T. Shimokawabe, T. Aoki, Two-dimensional phase-field study of competitive grain growth during directional solidification of polycrystalline binary alloy, Journal of Crystal Growth 442 (2016) 14–24.  doi:10.1016/j.jcrysgro.2016.01.036.
  * [52] T. Takaki, S. Sakane, M. Ohno, Y. Shibuta, T. Aoki, C.-A. Gandin, Competitive grain growth during directional solidification of a polycrystalline binary alloy: Three-dimensional large-scale phase-field study, Materialia 1 (2018) 104–113.  doi:10.1016/j.mtla.2018.05.002.
  * [53] S. Hu, L. Liu, W. Yang, J. Zhang, T. Huang, Y. Wang, X. Zhou, Competitive converging dendrites growth depended on dendrite spacing distribution of ni-based bi-crystal superalloys, Journal of Alloys and Compounds 735 (2018) 1878–1884.
  * [54] J. Eiken, Phase-field simulation of microstructure formation in technical magnesium alloys, International journal of materials research 101 (4) (2010) 503–509.  doi:10.3139/146.110302.
  * [55] J. Eiken, Phase-field simulations of dendritic orientation selection in mg-alloys with hexagonal anisotropy, in: Materials Science Forum, Vol. 649, Trans Tech Publ, 2010, pp. 199–204.
  * [56] K. Glasner, Nonlinear preconditioning for diffuse interfaces, Journal of Computational Physics 174 (2) (2001) 695–711.  doi:10.1006/jcph.2001.6933.
  * [57] A. Karma, W.-J. Rappel, Quantitative phase-field modeling of dendritic growth in two and three dimensions, Physical review E 57 (4) (1998) 4323.  doi:10.1103/PhysRevE.57.4323.
  * [58] D. Sun, M. Mendelev, C. Becker, K. Kudin, T. Haxhimali, M. Asta, J. Hoyt, A. Karma, D. J. Srolovitz, Crystal-melt interfacial free energies in hcp metals: A molecular dynamics study of mg, Physical Review B 73 (2) (2006) 024116\.  doi:10.1103/PhysRevB.73.024116.
  * [59] A. K. Boukellal, A. K. S. Elvalli, J.-M. Debierre, Equilibrium and growth facetted shapes in isothermal solidification of silicon: 3d phase-field simulations, Journal of Crystal Growth 522 (2019) 37–44.
  * [60] A. Karma, W.-J. Rappel, Phase-field model of dendritic sidebranching with thermal noise, Physical review E 60 (4) (1999) 3614.  doi:10.1103/PhysRevE.60.3614.
  * [61] J. Ghmadh, J.-M. Debierre, J. Deschamps, M. Georgelin, R. Guérin, A. Pocheau, Directional solidification of inclined structures in thin samples, Acta materialia 74 (2014) 255–267.  doi:10.1016/j.actamat.2014.04.023.
  * [62] H. Okamoto, T. Massalski, et al., Binary alloy phase diagrams, ASM International, Materials Park, OH, USA 12 (1990).
  * [63] J. Valencia, P. Quested, Asm handbook, volume 15: Casting, thermophysical properties (2008).
  * [64] A. Roy, R. Chhabra, Prediction of solute diffusion coefficients in liquid metals, Metallurgical Transactions A 19 (2) (1988) 273–279.  doi:10.1007/BF02652536.
  * [65] S. Wang, D. Liu, Y. Du, L. Zhang, Q. Chen, A. Engström, Development of an atomic mobility database for liquid phase in multicomponent al alloys: focusing on binary systems, International journal of materials research 104 (8) (2013) 721–735.  doi:10.3139/146.110923.
  * [66] https://periodictable.com.
  * [67] F. Bachmann, R. Hielscher, H. Schaeben, Grain detection from 2d and 3d ebsd data—specification of the mtex algorithm, Ultramicroscopy 111 (12) (2011) 1720–1733.  doi:10.1016/j.ultramic.2011.08.002.
  * [68] A. Clarke, D. Tourret, Y. Song, S. Imhoff, P. Gibbs, J. Gibbs, K. Fezzaa, A. Karma, Microstructure selection in thin-sample directional solidification of an al-cu alloy: In situ x-ray imaging and phase-field simulations, Acta Materialia 129 (2017) 203–216.
  * [69] Y. Song, D. Tourret, F. Mota, J. Pereda, B. Billia, N. Bergeon, R. Trivedi, A. Karma, Thermal-field effects on interface dynamics and microstructure selection during alloy directional solidification, Acta Materialia 150 (2018) 139–152.
  * [70] B. Bellón, A. Boukellal, T. Isensee, O. Wellborn, K. Trumble, M. Krane, M. Titus, D. Tourret, J. Llorca, Multiscale prediction of microstructure length scales in metallic alloy casting, Acta Materialia 207 (2021) 116686.  doi:10.1016/j.actamat.2021.116686.
  * [71] J. Hunt, Solidification and casting of metals, The Metal Society, London 3 (1979).
  * [72] T. Takaki, M. Ohno, T. Shimokawabe, T. Aoki, Two-dimensional phase-field simulations of dendrite competitive growth during the directional solidification of a binary alloy bicrystal, Acta Materialia 81 (2014) 272–283.  doi:10.1016/j.actamat.2014.08.035.
  * [73] R. Pieters, J. Langer, Noise-driven sidebranching in the boundary-layer model of dendritic solidification, Physical review letters 56 (18) (1986) 1948.
  * [74] M. N. Barber, A. Barbieri, J. Langer, Dynamics of dendritic sidebranching in the two-dimensional symmetric model of solidification, Physical Review A 36 (7) (1987) 3340.
  * [75] A. Dougherty, P. Kaplan, J. Gollub, Development of side branching in dendritic crystal growth, Physical review letters 58 (16) (1987) 1652.
  * [76] E. Brener, D. Temkin, Noise-induced sidebranching in the three-dimensional nonaxisymmetric dendritic growth, Physical Review E 51 (1) (1995) 351.
  * [77] S. Ghosh, A. Karma, M. Plapp, S. Akamatsu, S. Bottin-Rousseau, G. Faivre, Influence of morphological instability on grain boundary trajectory during directional solidification, Acta Materialia 175 (2019) 214–221.
  * [78] S. M. Elahi, R. Tavakoli, A. K. Boukellal, T. Isensee, I. Romero, D. Tourret, Multiscale simulation of powder-bed fusion processing of metallic alloys, Computational Materials Science 209 (2022) 111383.  doi:https://doi.org/10.1016/j.commatsci.2022.111383.
  * [79] A. K. Boukellal, J.-M. Debierre, H. Nguyen-Thi, G. Reinhart, Scaling laws governing the growth and interaction of equiaxed al-cu dendrites: a study combining experiments with phase-field simulations, Materialia 1 ((2018)) 62.