# A quantitative variational phase field framework

Arnab Mukherjee1,2<EMAIL_ADDRESS>James A. Warren2 Peter W.
Voorhees1,3,4 1Center for Hierarchical Materials Design, Northwestern
University, Evanston, Illinois 60208, USA  2Material Measurement Laboratory,
National Institute of Standards and Technology, Gaithersburg, Maryland 20899,
USA  3Department of Materials Science and Engineering, Northwestern
University, Evanston, Illinois 60208, USA  4Department of Engineering Sciences
and Applied Mathematics, Northwestern University, Evanston, Illinois 60208,
USA

###### Abstract

The finite solid-liquid interface width in phase field models results in non-
equilibrium effects, including solute trapping. Prior phase field modeling has
shown that this extra degree of freedom, when compared to sharp-interface
models, results in solute trapping that is well captured when realistic
parameters, such as interface width, are employed. However, increasing the
interface width, which is desirable for computational reasons, leads to
artificially enhanced trapping thus making it difficult to model departure
from equilibrium quantitatively. In the present work, we develop a variational
phase field model that guarantees a temporal decrease in the free energy with
independent kinetic equations for the solid and liquid phases. Separate
kinetic equations for the phase concentrations obviate the assumption of point
wise equality of diffusion potentials, as is done in previous works. Non-
equilibrium effects such as solute trapping, drag and interface kinetics can
be introduced in a controlled manner in the present model. In addition, the
model parameters can be tuned to obtain “experimentally-relevant” trapping
while using significantly larger interface widths than prior efforts. A
comparison with these other phase field models suggests that interface width
of about three to twenty-five times larger than current best-in-class models
can be employed depending upon the material system at hand leading to a speed-
up by a factor of $W^{(d+2)}$, where $W$ and $d$ denote the interface width
and spatial dimension, respectively. Finally the capacity to model non-
equilibrium phenomena is demonstrated by simulating oscillatory instability
leading to the formation of solute bands.

##  I Introduction

Rapid solidification refers to the use of large cooling rates or undercoolings
resulting in growth speeds of solid-liquid front typically of the order of a
$\text{\,}\mathrm{cm}\text{\,}{\mathrm{s}}^{-1}$ or larger. Rapid
solidification processing encompasses a wide range of industrially-relevant
manufacturing processes such as laser-induced surface melting [1], spray
forming [2], welding [3] and glass formation [4] amongst others. With the
recent surge in additive manufacturing processes, there has been a renewed
interest in experimental [5, 6, 7], theoretical development and numerical
modeling [8, 9, 10, 11, 12] of rapid solidification of metallic alloys.

Rapid solidification is characterized by the breakdown of the local
interfacial equilibrium assumption. Owing to the high growth speeds, the
solute atoms do not have enough time to diffuse through the characteristic
solid-liquid interface width $\delta$ which is of the order of inter-atomic
spacing. As a result solute is incorporated in the growing solid phase at a
concentration significantly different from that predicted by the phase
diagram. This non-equilibrium or incomplete solute partitioning at the solid-
liquid interface is referred to as “solute trapping.” Cognizance of solute
trapping holds the key in understanding the growth morphology [13, 14],
formation of secondary phases [6] and micro-segregation pattern of the
resulting microstructure [7].

Many of the insights into the phenomenon of solute trapping have followed from
analyses modeling the liquid-solid interface using either a sharp-interface
[15, 16, 17, 18, 19] or a semi-sharp interface or continuum model [20, 21].
Typically, the sharp interface models employ a solution to the diffusion
equation to describe the bulk phases which are subjected to external and
interfacial conditions. The interfacial conditions are described by response
functions, one each for the interfacial concentration and temperature [22,
16]. While most sharp-interface theories differ in the details of the
interface response functions, most of them predict a velocity-dependent
partition coefficient $k_{v}$, which increases monotonically from the
equilibrium partition coefficient $k_{e}$ to $1$ as velocity increases from
zero and tends to infinity, and an interface temperature that varies with
velocity. The continuum models go a step further to incorporate the variation
of the chemical potential and diffusivity through the interface which leads to
non-monotonic dependence of $k_{v}$ on the velocity [22]. Free energy
dissipated by solute diffusion through the interface during interface motion
leads to solute drag [21, 16]. While the sharp-interface and continuum models
provide a thermodynamically-consistent framework to study trapping, tracking
the morphologically complex interface makes implementing such models in higher
dimensions an arduous task.

Phase field models have enjoyed great success in simulating microstructure
evolution in myriad problems in solidification [23, 24] and solid-state phase
transformation [25, 26, 27]. Unlike the sharp-interface models, the efficacy
of the phase field method lies in its treatment of the interface. An indicator
variable $\phi$, known as the phase field is introduced to demarcate the
regions of different phases, say, solid and liquid by assigning constant
values (say $\phi=0,1$) in the bulk phases. The phase field varies smoothly
between the bulk values in a thin transition region known as the diffuse
interface width $W$. While the indicator function circumvents the need for
interface tracking, the diffuse interface thickness introduces an additional
length scale in the formulation. $W$ can be treated as either (i) a physical
entity $\delta$, in which case the value is typically of the order of Angstrom
or a few nanometers, or (ii) a mere mathematical construct, such that $W$ can
be enlarged by orders of magnitude than the physical solid-liquid interface
thickness. The value of $W$ is restricted by the condition to resolve the
smallest microstructural feature. Ascribing the latter notion to $W$,
microstructures of the order of few $\text{\,}\mathrm{\SIUnitSymbolMicro m}$
can be simulated within reasonable times. However, the computation results are
now dependent on $W$ and an additional step of assessing the $W$ dependence on
results needs to be performed. This can be done numerically by performing
convergence tests or analytically by the method of matched asymptotic
analysis.

The finite interface width in phase field models introduces a number of non-
equilibrium effects including solute trapping [28, 29]. Depending upon the
solidification conditions the solute trapping might be physically realistic or
an undesirable artifact. It is computationally beneficial to use as large an
interface width as possible, while preserving the required physical behavior
of the model. For low speed directional solidification, solute trapping is
negligible and can be eliminated using a judicious choice of interpolation
functions [30, 29] and incorporation of an “antitrapping current” [28, 31],
the form of which is selected based on the thin-interface analysis to recover
local equilibrium conditions. Since the antitrapping current does not
contribute to the free energy dissipation such models are termed as non-
variational models. So, while solute trapping can and does occur during rapid
solidification, phase field models with large values of $W$ chosen for
computational expedience will yield artificially enhanced trapping [32, 33].

Ahmad and co-workers [34] studied solute trapping in context of the Wheeler-
Boettinger-McFadden (WBM) phase field model [35]. In the WBM model, the
interface is can be interpreted as a mixture of solid and liquid with the same
composition. Thus, a single variable $c$ tracks the solute concentration in
both phases. Selecting $W$ of the order of 0.1 $\text{\,}\mathrm{nm}$, the
variation of the partition coefficient with velocity was found to be in good
agreement with the sharp-interface continuous growth model (CGM). In addition,
the interface temperature agrees well with CGM with partial solute drag.
Trapping has also been studied for the Kim-Kim-Suzuki (KKS) phase field model
[36]. Unlike the WBM model, the KKS phase field model treats the interface as
a mixture of solid and liquid with distinct concentrations $c_{s}$ and $c_{l}$
respectively. However, the phase concentrations are mutually dependent via the
constraint of equal diffusion potential at each spatial point. The imposition
of pointwise equality of diffusion potential, however, does not preclude
solute trapping and the far-field diffusion potentials in the solid and liquid
exhibit a jump that was shown to be consistent with the high velocity
asymptotics of the WBM model. Similar conclusions have been deduced by a
number of studies employing different versions of the above models on
different material systems [37, 38, 39, 33, 40, 41, 42]. A common feature of
the above models is that the partition coefficient was found to be sensitive
to the chosen interface width, diffusivity and free energy interpolation
functions. The interpolation functions (say $h(\phi)$), which connect bulk
quantities through the interface, are usually constructed based on the
restriction that $h(0)=0$, $h(1)=1$ and $h^{\prime}(\phi=0,1)=0$. Since the
interpolation functions are not associated with any inherent length scale, $W$
is the only free parameter in the model to tailor partitioning. It is evident
that additional degrees of freedom are needed if certain control over solute
trapping and drag is to be introduced in the phase field models. Moreover, the
choice of $W$ puts a severe limitation on the computation time, especially in
higher dimensions. The recent advances have therefore focused on addressing
the issue of introducing additional degrees of freedom such that interface
thickness can be increased, while still maintaining experimentally/physically-
relevant trapping.

A finite interface dissipation model was recently introduced by Steinbach et
al. [43]. Similar to the KKS model, distinct phase concentrations are
introduced. However, rather than assuming local equilibrium condition,
separate kinetic equations for the phase concentrations $c_{s}$ and $c_{l}$
were derived and the rate of exchange of atoms across phases is determined by
a kinetic parameter termed as the “interface permeability” which provides the
necessary degree of freedom. The solute trapping in the model is controlled
through a combination of values of the interface permeability and interface
thickness. However, the model is not analytically tractable making it
difficult to select appropriate model parameters. In addition, it is unclear
whether the interface permeability can be tuned to additionally control solute
drag or interface temperature. Pinomaa and Provatas [44] have extended the
thin-interface asymptotics to map the dilute alloy model of Echebarria-Folch-
Karma-Plapp (EFKP) [45] to the sharp-interface CGM. A desired solute
partitioning is attained by tuning the “anti-trapping current.” A recent model
by Ji et al. [46] controls the partitioning by introducing an additional free
parameter in the diffusivity interpolation function to increase the
diffusivity at the interface. For different choices of the interface width,
this additional parameter can be selected to give rise to same trapping
behaviour. Both the above models are valid for dilute alloys and one-sided
diffusivity i.e. the diffusion in the solid phase is zero.

A new class of phase field models termed as “non-diagonal” have also been
introduced independently by several authors [47, 48, 49, 50, 51, 52].
Diagonality refers to the fact that in the previous phase field models such as
WBM and KKS, the temporal derivative of concentration and the phase field
depends only on the functional derivative of the free energy with respect to
concentration (driving force for diffusion) and phase field (driving force for
phase transformation) respectively. In other words, the diffusion flux is
proportional to driving force for diffusion and rate of interface motion is
coupled to the driving force for phase transformation. The dissipative flux in
the bulk and the interface is coupled in a single diffusion flux term.
However, the “non-diagonal” models introduce a kinetic cross-coupling between
solute diffusion and phase-transformation. This cross-coupling term introduces
an additional degree of freedom such that the bulk diffusion and interfacial
solute redistribution can be controlled independently. Furthermore, the
dynamic cross-coupling restores thermodynamic consistency in phase field
models with the cross-coupling terms naturally emerging in the kinetic
equations of concentration and phase field respectively. The non-diagonal
models have opened the possibility to model close-to-equilibrium
solidification processes quantitatively by successfully eliminating all the
undesirable interface kinetic effects in systems with finite diffusion in
solid [53, 54, 52]. The full potential of non-diagonal models, however, have
not been fully explored for non-equilibrium solidification problems. A recent
phase field model has been proposed by Ohno and co-workers [52] where the
kinetic cross-coupling terms and their forms naturally arise in the model
derivation. However, the model is only applicable to a low speed
solidification regime. A local equilibrium condition is imposed to relate the
phase concentrations. In addition, there are no degrees of freedom to control
non-equilibrium effects. With this motivation in mind, we develop a
variational phase field model that guarantees a decrease in free energy with
sufficient degrees of freedom to incorporate non-equilibrium effects relevant
to rapid solidification processing in a controlled manner. A key feature of
the present model is that the phase concentrations are treated as mutually
independent rather than relating them via the local equilibrium condition.
Separate kinetic equations with interfacial source terms are derived for the
phase concentrations which can be tuned to control solute trapping and drag.
Unlike the KKS model and its derivatives, separate kinetic equations alleviate
the need for the iterative solution of non-linear system of equations to
relate the phase concentrations which can be computationally time consuming
depending upon the alloy model and complexity of the system (for instance
multi-phase multi-component systems). A number of previous phase field models
are shown to be special cases of this new formulation.

The article is organized as follows: the model framework derived from
variational principles is presented in Section II. We perform a detailed
asymptotic analysis to derive the expression for the diffusion potential jump.
The equations are then specialized for a dilute alloy in Section III. It is
shown in Section IV that a number of previous phase field models can be
derived as special cases from the present model. Numerical results are
presented in Section V providing insights on tailoring a desired level for
solute trapping and drag. The model parameters can be tuned readily from
available experimental datasets, atomistic simulations or sharp-interface
theories to capture “correct” level of partitioning. The model is applied to
Si–$9\text{\,}\%$As and Ni–$5\text{\,}\%$Cu and the efficacy of the model is
demonstrated by numerical studies on the partition coefficient and interface
temperature. A comparative study with the previously developed phase field
model suggests that interface width can be enlarged by a factor of three to
twenty-five times in the velocity regime of
$0.01\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ to
$1\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$. For a finite-difference
algorithm on a uniform grid spacing this upscaling expedites the computation
by a factor $W^{(d+2)}$, where, $d$ is the dimension. Preliminary results on
an oscillatory instability leading to the formation of solute bands are also
presented. Finally, the paper is concluded by a summary of the main findings
in Section VI.

##  II Model

We consider a binary two-phase system (say, solid and liquid) such that the
phase field variable $\phi$ assumes a value of $1$ in the solid and $0$ in the
liquid. The interface, where $\phi$ varies smoothly between the bulk values is
treated as a mixture of phases, each with a distinct composition $c_{s}$ and
$c_{l}$. The overall solute concentration is defined as

$c=c_{s}g(\phi)+c_{l}\\{1-g(\phi)\\},$ (1)

where $g(\phi)$ is a smoothed step-function satisfying the constraints
$g(\phi=0)=0$, $g(\phi=1)=1$ and $g^{\prime}(\phi=0,1)=0$. The free energy
functional is written as a contribution arising from the usual interfacial and
bulk free energy density and taking in account the constraint Eq. (1) as

$\overline{\mathcal{F}}=\mathcal{F}+\int_{V}\lambda\Big{[}c-\\{c_{s}g(\phi)+c_{l}(1-g(\phi))\\}\Big{]}\mathrm{d}V,$
(2)

where

$\mathcal{F}=\int_{V}\Big{[}Hf_{\mathrm{dw}}(\phi)+\frac{\sigma}{2}|\nabla\phi|^{2}+f_{s}(c_{s})g(\phi)+f_{l}(c_{l})\\{1-g(\phi)\\}\Big{]}\mathrm{d}V.$
(3)

$f_{\mathrm{dw}}(\phi)$ is the double well potential with a quartic form
$\phi^{2}(1-\phi)^{2}$ with minima at $\phi=0$ and $1$. $H$ is the barrier
height. $\sigma$ is the gradient energy coefficient penalizing sharp gradients
and imparting smoothness to the phase field profile. $f_{s}(c_{s})$ and
$f_{l}(c_{l})$ are the bulk free densities of the solid and liquid
respectively. $\lambda$ is the Lagrange multiplier to ensure that constraint
Eq. (1) is obeyed. To select suitable choices of the evolution equations of
the field variables, we consider the temporal change of free energy functional
as

$\frac{\mathrm{d}\overline{\mathcal{F}}}{\mathrm{d}t}=\int_{V}\Bigg{[}\frac{\delta\overline{\mathcal{F}}}{\delta
c_{s}}\frac{\partial c_{s}}{\partial
t}+\frac{\delta\overline{\mathcal{F}}}{\delta c_{l}}\frac{\partial
c_{l}}{\partial
t}+\frac{\delta\overline{\mathcal{F}}}{\delta\phi}\frac{\partial\phi}{\partial
t}+\frac{\delta\overline{\mathcal{F}}}{\delta c}\frac{\partial c}{\partial
t}+\frac{\delta\overline{\mathcal{F}}}{\delta\lambda}\frac{\partial\lambda}{\partial
t}\Bigg{]}\mathrm{d}V.$ (4)

Introducing the individual variation with respect to $c_{s}$, $c_{l}$, $\phi$,
$c$ and $\lambda$ form Eqs. (2) and (3) we obtain

$\frac{\mathrm{d}\overline{\mathcal{F}}}{\mathrm{d}t}=\int_{V}\Bigg{\\{}\Bigg{[}\frac{\delta\mathcal{F}}{\delta
c_{s}}-\lambda g(\phi)\Bigg{]}\frac{\partial c_{s}}{\partial
t}+\Bigg{[}\frac{\delta\mathcal{F}}{\delta
c_{l}}-\lambda\\{1-g(\phi)\\}\Bigg{]}\frac{\partial c_{l}}{\partial
t}+\Bigg{[}\frac{\delta\mathcal{F}}{\delta\phi}-\lambda(c_{s}-c_{l})g^{\prime}(\phi)\Bigg{]}\frac{\partial\phi}{\partial
t}+\lambda\frac{\partial c}{\partial t}\Bigg{\\}}\mathrm{d}V.$ (5)

The conservation of mass is captured through the solute concentration $c$
continuity equation as

$\frac{\partial c}{\partial t}=-\nabla\cdot J.$ (6)

The net flux $J$ is assumed to be a simple interpolation of the solid and
liquid flux as

$J=J_{s}g(\phi)+J_{l}\\{1-g(\phi)\\}.$ (7)

The temporal evolution of variables $c_{s}$ and $c_{l}$ can be postulated by
noting that these quantities need only be conserved in their respective bulk
phases when $c=c_{s}$ as $\phi=1$ and $c=c_{l}$ as $\phi=0$. In the
interfacial region, a source term is permitted. A possible choice can be
written as

$\frac{\partial c_{s}}{\partial t}=-\nabla\cdot
J_{s}+S_{s}(c_{s},c_{l},\phi)g^{\prime}(\phi)\frac{\partial\phi}{\partial t},$
(8)

and

$\frac{\partial c_{l}}{\partial t}=-\nabla\cdot
J_{l}+S_{l}(c_{s},c_{l},\phi)g^{\prime}(\phi)\frac{\partial\phi}{\partial t}.$
(9)

The source terms in Eqs. (8) and (9) are proportional to
$\frac{\partial\phi}{\partial t}$ to take into account that the solute
redistribution is dependent on the interface velocity. The functions $S_{s}$
and $S_{l}$ in the source terms introduce a degree of freedom, which in
principle can be tailored to map the phase field model onto a desired set of
sharp-interface relations. We will see later that the source functions $S_{s}$
and $S_{l}$ contribute to the free energy dissipation. However, the free
energy dissipation does not depend explicitly on the forms or the sign of the
source functions. The exact form of the interfacial source terms $S_{s}$ and
$S_{l}$ will be deduced later based on the asymptotic analysis of the model.
Substituting Eqs. (6), (7), (8) and (9) in Eq. (5), applying the divergence
theorem and neglecting the surface integrals, we have

$\displaystyle\frac{\mathrm{d}\overline{\mathcal{F}}}{\mathrm{d}t}=\int_{V}\Bigg{[}\nabla\Big{(}\frac{\delta\mathcal{F}}{\delta
c_{s}}\Big{)}-\lambda g^{\prime}(\phi)\nabla\phi\Bigg{]}\cdot
J_{s}+\Bigg{[}\nabla\Big{(}\frac{\delta\mathcal{F}}{\delta
c_{l}}\Big{)}+\lambda g^{\prime}(\phi)\nabla\phi\Bigg{]}\cdot J_{l}$
$\displaystyle+\Bigg{[}\frac{\delta\mathcal{F}}{\delta\phi}-\lambda(\tilde{c_{s}}-\tilde{c_{l}})g^{\prime}(\phi)+\Big{(}S_{s}\frac{\delta\mathcal{F}}{\delta
c_{s}}+S_{l}\frac{\delta\mathcal{F}}{\delta
c_{l}}\Big{)}g^{\prime}(\phi)\Bigg{]}\frac{\partial\phi}{\partial t},$ (10)

where
$\tilde{c_{s}}-\tilde{c_{l}}=c_{s}-c_{l}+S_{s}g(\phi)+S_{l}\\{1-g(\phi)\\}$.
We note that the terms associated with Lagrange multiplier in the square
brackets that are dotted to the fluxes $J_{s}$ and $J_{l}$ are directed normal
to the interface (because of the presence of the operator $\nabla\phi$).
Therefore, decomposing the gradient operators and fluxes in the normal and
tangential directions with respect to the interface as
$\nabla(\delta\mathcal{F}/\delta
c_{s,(l)})=\nabla_{n}(\delta\mathcal{F}/\delta
c_{s,(l)})+\nabla_{t}(\delta\mathcal{F}/\delta c_{s,(l)})$ and
$J_{s,(l)}=J_{s,(l)}^{n}+J_{s,(l)}^{t}$ where the subscripts $\mathrm{n}$ and
$\mathrm{t}$ denote the normal and tangential directions with respect to the
interface, respectively. The normal and tangential derivatives are expressed
in terms of projection operators as
$\nabla_{n}=(\mathbf{n}\otimes\mathbf{n})\nabla$ and
$\nabla_{t}=(\mathbf{1}-\mathbf{n}\otimes\mathbf{n})\nabla$ respectively.
$\mathbf{n}$ denotes the interface normal defined by
$\mathbf{n}=-\frac{\nabla\phi}{|\nabla\phi|}$ and $\mathbf{1}$ is a unit
tensor allows Eq. (II) to be re-written as

$\displaystyle\frac{\mathrm{d}\overline{\mathcal{F}}}{\mathrm{d}t}=\int_{V}$
$\displaystyle\Bigg{[}\nabla_{n}\Big{(}\frac{\delta\mathcal{F}}{\delta
c_{s}}\Big{)}-\lambda g^{\prime}(\phi)\nabla\phi\Bigg{]}\cdot
J_{s}^{n}+\nabla_{t}\Big{(}\frac{\delta\mathcal{F}}{\delta c_{s}}\Big{)}\cdot
J_{s}^{t}$ $\displaystyle+$
$\displaystyle\Bigg{[}\nabla_{n}\Big{(}\frac{\delta\mathcal{F}}{\delta
c_{l}}\Big{)}+\lambda g^{\prime}(\phi)\nabla\phi\Bigg{]}\cdot
J_{l}^{n}+\nabla_{t}\Big{(}\frac{\delta\mathcal{F}}{\delta c_{l}}\Big{)}\cdot
J_{l}^{t}$ $\displaystyle+$
$\displaystyle\Bigg{[}\frac{\delta\mathcal{F}}{\delta\phi}-\lambda(\tilde{c_{s}}-\tilde{c_{l}})g^{\prime}(\phi)+\Big{(}S_{s}\frac{\delta\mathcal{F}}{\delta
c_{s}}+S_{l}\frac{\delta\mathcal{F}}{\delta
c_{l}}\Big{)}g^{\prime}(\phi)\Bigg{]}\frac{\partial\phi}{\partial t},$ (11)

To guarantee a temporal decrease in the free energy we choose the forms of
fluxes as

$J_{s}^{n}=-{M_{s}^{n}(\phi)}\Bigg{[}\nabla_{n}\Big{(}\frac{\delta\mathcal{F}}{\delta
c_{s}}\Big{)}-\lambda g^{\prime}(\phi)\nabla\phi\Bigg{]},$ (12)
$J_{s}^{t}=-M_{s}^{t}(\phi)\nabla_{t}\Big{(}\frac{\delta\mathcal{F}}{\delta
c_{s}}\Big{)},$ (13)
$J_{l}^{n}=-{M_{l}^{n}(\phi)}\Bigg{[}\nabla_{n}\Big{(}\frac{\delta\mathcal{F}}{\delta
c_{l}}\Big{)}+\lambda g^{\prime}(\phi)\nabla\phi\Bigg{]},$ (14)
$J_{l}^{t}=-M_{l}^{t}(\phi)\nabla_{t}\Big{(}\frac{\delta\mathcal{F}}{\delta
c_{l}}\Big{)},$ (15)

and,

$\frac{\partial\phi}{\partial
t}=-M_{\phi}\Bigg{[}\frac{\delta\mathcal{F}}{\delta\phi}-\lambda(\tilde{c_{s}}-\tilde{c_{l}})g^{\prime}(\phi)+\Big{(}S_{s}\frac{\delta\mathcal{F}}{\delta
c_{s}}+S_{l}\frac{\delta\mathcal{F}}{\delta
c_{l}}\Big{)}g^{\prime}(\phi)\Bigg{]}.$ (16)

such that $\frac{\mathrm{d}\overline{\mathcal{F}}}{\mathrm{d}t}\leq 0$ as long
as $M_{s}^{n}(\phi)$, $M_{s}^{t}(\phi)$, $M_{l}^{n}(\phi)$, $M_{l}^{t}(\phi)$
and $M_{\phi}>0$. $M_{s(l)}^{n}(\phi)$ and $M_{s(l)}^{t}(\phi)$ are the atomic
mobilities of the component in the solid (liquid) phase in the normal and
tangential directions respectively. It is to be noted that the free energy
dissipation does not depend on the form or the sign of source functions
$S_{s}$ and $S_{l}$. The normal and tangential decomposition of the fluxes
provides an additional degree of freedom in choosing different interpolation
functions of the mobility in the normal and tangential directions. We will see
later that this additional degree of freedom offers independent control over
two thin-interface effects namely Kapitza jump and surface diffusion.

To evaluate the Lagrange multiplier, we employ the overall solute conservation
Eq. (6) at each point in the domain. Substituting Eq. (1), (7) and the kinetic
Eqs. (8) and (9) in Eq. (6) we obtain

$(\tilde{c_{s}}-\tilde{c_{l}})\frac{\partial\phi}{\partial
t}=-(J_{s}^{n}-J_{l}^{n})\cdot\nabla\phi.$ (17)

The above relation (which is a Stefan-like condition) relates the three
dissipative fluxes at each spatial point in the interface. Further,
substituting $J_{s}^{n}$ and $J_{l}^{n}$ from Eqs. (12) and (14) we have

$\displaystyle\lambda g^{\prime}(\phi)=\frac{1}{\overline{M}(\phi)}\Bigg{[}$
$\displaystyle-\frac{(\tilde{c_{s}}-\tilde{c_{l}})}{|\nabla\phi|^{2}}\frac{\partial\phi}{\partial
t}$
$\displaystyle+\frac{\Bigg{\\{}M_{s}^{n}(\phi)\nabla_{n}\Big{(}\frac{\delta\mathcal{F}}{\delta
c_{s}}\Big{)}-M_{l}^{n}(\phi)\nabla_{n}\Big{(}\frac{\delta\mathcal{F}}{\delta
c_{l}}\Big{)}\Bigg{\\}}}{|\nabla\phi|}\cdot\frac{\nabla\phi}{|\nabla\phi|}\Bigg{]},$
(18)

where,

$\overline{M}(\phi)=M_{s}^{n}(\phi)+M_{l}^{n}(\phi).$ (19)

Substituting Eq. (II) in Eqs. (12)-(16) we obtain

$\displaystyle J_{s}=$
$\displaystyle-\frac{M_{s}^{n}(\phi)M_{l}^{n}(\phi)}{\overline{M}(\phi)}\nabla_{n}\Bigg{(}\frac{\partial
f_{s}}{\partial c_{s}}g(\phi)+\frac{\partial f_{l}}{\partial
c_{l}}\\{1-g(\phi)\\}\Bigg{)}$ $\displaystyle-
M_{s}^{t}(\phi)g(\phi)\nabla_{t}\frac{\partial f_{s}}{\partial
c_{s}}-\frac{(\tilde{c_{s}}-\tilde{c_{l}})M_{s}^{n}(\phi)}{\overline{M}(\phi)}\frac{1}{|\nabla\phi|}\frac{\nabla\phi}{|\nabla\phi|}\frac{\partial\phi}{\partial
t}$ (20) $\displaystyle J_{l}=$
$\displaystyle-\frac{M_{s}^{n}(\phi)M_{l}^{n}(\phi)}{\overline{M}(\phi)}\nabla_{n}\Bigg{(}\frac{\partial
f_{s}}{\partial c_{s}}g(\phi)+\frac{\partial f_{l}}{\partial
c_{l}}\\{1-g(\phi)\\}\Bigg{)}$ $\displaystyle-
M_{l}^{t}(\phi)\\{1-g(\phi)\\}\nabla_{t}\frac{\partial f_{l}}{\partial
c_{l}}+\frac{(\tilde{c_{s}}-\tilde{c_{l}})M_{l}^{n}(\phi)}{\overline{M}(\phi)}\frac{1}{|\nabla\phi|}\frac{\nabla\phi}{|\nabla\phi|}\frac{\partial\phi}{\partial
t}$ (21)

where, the normal and tangential components have been added to obtain $J_{s}$
and $J_{l}$ and

$\displaystyle\Bigg{[}\frac{1}{M_{\phi}}+\frac{(\tilde{c_{s}}-\tilde{c_{l}})^{2}}{\overline{M}(\phi)}\frac{1}{|\nabla\phi|^{2}}\Bigg{]}$
$\displaystyle\frac{\partial\phi}{\partial t}=$
$\displaystyle\frac{M_{s}^{n}(\phi)}{\overline{M}(\phi)}\Bigg{[}\sigma\nabla^{2}\phi-
Hf_{\mathrm{dw}}^{\prime}(\phi)-\Big{(}f_{s}-f_{l}-(c_{s}-c_{l})\frac{\partial
f_{s}}{\partial c_{s}}$ $\displaystyle\hskip 42.67912pt-
S_{l}\\{1-g(\phi)\\}(\frac{\partial f_{s}}{\partial c_{s}}-\frac{\partial
f_{l}}{\partial c_{l}})\Big{)}g^{\prime}(\phi)\Bigg{]}+$
$\displaystyle\frac{M_{l}^{n}(\phi)}{\overline{M}(\phi)}\Bigg{[}\sigma\nabla^{2}\phi-
Hf_{\mathrm{dw}}^{\prime}(\phi)-\Big{(}f_{s}-f_{l}-(c_{s}-c_{l})\frac{\partial
f_{l}}{\partial c_{l}}$ $\displaystyle\hskip 42.67912pt-
S_{s}g(\phi)(\frac{\partial f_{l}}{\partial c_{l}}-\frac{\partial
f_{s}}{\partial c_{s}})\Big{)}g^{\prime}(\phi)\Bigg{]}+$
$\displaystyle\frac{(\tilde{c_{s}}-\tilde{c_{l}})}{\overline{M}(\phi)}\frac{1}{|\nabla\phi|}\Bigg{(}M_{s}^{n}(\phi)g(\phi)\nabla\frac{\partial
f_{s}}{\partial c_{s}}-M_{l}^{n}(\phi)\\{1-g(\phi)\\}\nabla\frac{\partial
f_{l}}{\partial c_{l}}\Bigg{)}\cdot\frac{\nabla\phi}{|\nabla\phi|}.$ (22)

Inserting Eqs. (II) and (II) in Eqs. (7) and (8) the individual kinetic
equations for the phase concentrations can be derived as

$\displaystyle\frac{\partial c_{s}}{\partial t}=\nabla\cdot\Bigg{[}$
$\displaystyle\frac{M_{s}^{n}(\phi)M_{l}^{n}(\phi)}{\overline{M}(\phi)}\nabla_{n}\Bigg{(}\frac{\partial
f_{s}}{\partial c_{s}}g(\phi)+\frac{\partial f_{l}}{\partial
c_{l}}\\{1-g(\phi)\\}\Bigg{)}$
$\displaystyle+M_{s}^{t}(\phi)g(\phi)\nabla_{t}\frac{\partial f_{s}}{\partial
c_{s}}$
$\displaystyle+\frac{(\tilde{c_{s}}-\tilde{c_{l}})}{|\nabla\phi|}\frac{M_{s}^{n}(\phi)}{\overline{M}(\phi)}\frac{\nabla\phi}{|\nabla\phi|}\frac{\partial\phi}{\partial
t}\Bigg{]}+S_{s}g^{\prime}(\phi)\frac{\partial\phi}{\partial t},$ (23)

and,

$\displaystyle\frac{\partial c_{l}}{\partial t}=\nabla\cdot\Bigg{[}$
$\displaystyle\frac{M_{s}^{n}(\phi)M_{l}^{n}(\phi)}{\overline{M}(\phi)}\nabla_{n}\Bigg{(}\frac{\partial
f_{s}}{\partial c_{s}}g(\phi)+\frac{\partial f_{l}}{\partial
c_{l}}\\{1-g(\phi)\\}\Bigg{)}$
$\displaystyle+M_{l}^{t}(\phi)(1-g(\phi))\nabla_{t}\frac{\partial
f_{l}}{\partial c_{l}}$
$\displaystyle-\frac{(\tilde{c_{s}}-\tilde{c_{l}})}{|\nabla\phi|}\frac{M_{l}^{n}(\phi)}{\overline{M}(\phi)}\frac{\nabla\phi}{|\nabla\phi|}\frac{\partial\phi}{\partial
t}\Bigg{]}+S_{l}g^{\prime}(\phi)\frac{\partial\phi}{\partial t}.$ (24)

Eqs. (II), (II) and (II) constitute the equations of motions for the phase
field model. The phase concentrations $c_{s}$ and $c_{l}$ are mutually
independent. This independence is in contrast to models like KKS, where the
assumption of pointwise equality of diffusion potential
$\mu_{s}=\frac{\partial f_{s}}{\partial c_{s}}=\mu_{l}=\frac{\partial
f_{l}}{\partial c_{l}}$ is the second equation that must be solved at each
point in space and time to relate the phase concentrations. Rather, the
kinetic Eqs. (II) and (II) can be solved to evaluate the respective phase
concentrations. At the same time, the present model retains the desirable
feature of the KKS model that at equilibrium, it decouples the bulk and the
interfacial contributions of the free energy. This implies that the
equilibrium phase field profile is given by a hyperbolic tangent function with
interface width $W$ given by

$W=\sqrt{\frac{\sigma}{H}}.$ (25)

The interfacial energy can be expressed as

$\gamma=IWH,$ (26)

where $I$ is a constant which depends on the chosen double-well potential. The
details of the equilibrium properties of the model is provided in the
supplementary materials for the sake of brevity. Using Eqs. (II) and (II) in
Eq. (6), the overall solute conservation equation is

$\displaystyle\frac{\partial c}{\partial t}=\nabla\cdot\Bigg{[}$
$\displaystyle q_{n}(\phi)\nabla_{n}\Bigg{(}\frac{\partial f_{s}}{\partial
c_{s}}g(\phi)+\frac{\partial f_{l}}{\partial c_{l}}\\{1-g(\phi)\\}\Bigg{)}$
$\displaystyle+M_{s}^{t}(\phi)g(\phi)^{2}\nabla_{t}\frac{\partial
f_{s}}{\partial
c_{s}}+M_{l}^{t}(\phi)\\{1-g(\phi)\\}^{2}\nabla_{t}\frac{\partial
f_{l}}{\partial c_{l}}$
$\displaystyle+\frac{(\tilde{c_{s}}-\tilde{c_{l}})}{|\nabla\phi|}a(\phi)\frac{\nabla\phi}{|\nabla\phi|}\frac{\partial\phi}{\partial
t}\Bigg{]},$ (27)

where,

$q_{n}(\phi)=\frac{M_{s}^{n}(\phi)M_{l}^{n}(\phi)}{\overline{M}(\phi)}$ (28)

and,

$a(\phi)=\frac{M_{s}^{n}(\phi)g(\phi)-M_{l}^{n}(\phi)\\{1-g(\phi)\\}}{\overline{M}(\phi)}$
(29)

Several important features of the governing equations are noteworthy. First,
the term proportional to $\frac{\partial\phi}{\partial t}$ (last term) in Eq.
(II) is directly related to solute trapping as it is dependent on the velocity
(through $\frac{\partial\phi}{\partial t}$), the local concentration
difference between solid and liquid and directed normal to the interface
(through the operator $\frac{\nabla\phi}{|\nabla\phi|}$). Thus, this is an
interface dissipation flux, having the mathematical characteristics of an
anti-trapping term that has been employed in previous works [28, 45]. The
anti-trapping term is split into two parts in the respective phase
concentrations in Eqs. (II) and (II). The corresponding conjugate term
proportional to $(M_{s}^{n}(\phi)\nabla\frac{\partial f_{s}}{\partial
c_{s}}-M_{l}^{n}(\phi)\nabla\frac{\partial f_{l}}{\partial
c_{l}})\cdot\mathbf{n}$ naturally arises during the evolution of the phase
field Eq. (II). Secondly, the diffusion flux in Eq. (II) is decomposed into
normal and tangential direction. The normal and tangential derivatives are
associated with diffusion potential jump and surface diffusion respectively.
The two interpolation functions can be chosen to eliminate two important thin
interface effects: namely the Kapitza resistance, or jump in the diffusion
potentials at a stationary interface when there is a mass flux through the
interface, and surface diffusion as reported in Ref.[55] and as will be
discussed next.

###  II.1 Choice of interpolation functions

The normal and tangential derivatives offer independent control over the
diffusion potential jump due to solute transport across a stationary interface
and surface diffusion respectively. It was shown by Nicoli et al.[55] that
both these thin interface effects can be remedied by a choice of inverse and
linear interpolation of the mobilities in the normal and tangential direction
respectively. The rationale behind such choices also follow from the
asymptotic analysis discussed in the next Section II.2. While the mobilities
were chosen to be tensorial in the former work, the decomposition of the
derivatives arise naturally in the model formulation in the present work. We
note that tensorial mobilities have also been employed to capture surface
diffusion in phase field models [56, 57]. Following the above arguments we
choose $M_{s}^{n}(\phi)=M_{s}/h(\phi)$ and
$M_{l}^{n}(\phi)=M_{l}/\\{1-h(\phi)\\}$ with $h(\phi)$ as an anti-symmetric
function which can be chosen to be the usual cubic
$h(\phi)=\phi^{2}(3-2\phi)$, or the quintic potential
$h(\phi)=\phi^{3}(6\phi^{2}-15\phi+10)$. Substituting these forms in Eq. (28)
we obtain

$\frac{1}{q_{n}(\phi)}=\frac{h(\phi)}{M_{s}}+\frac{1-h(\phi)}{M_{l}}.$ (30)

The mobilities in the tangential direction are also chosen as
$M_{s}^{t}(\phi)=M_{s}/h(\phi)$ and $M_{l}^{t}(\phi)=M_{l}/\\{1-h(\phi)\\}$.
In addition we also select $g(\phi)=h(\phi)$ to interpolate phase
concentrations, free energies and fluxes in Eqs. (1), (3) and (7). The above
choice of interpolation functions naturally leads to a linear interpolation of
the diffusion potentials of solid and liquid phase in the tangential direction
as $M_{s}h(\phi)\nabla_{t}(\partial f_{s}/\partial
c_{s})+M_{l}\\{1-h(\phi)\\}\nabla_{t}(\partial f_{l}/\partial c_{l})$. The
form of the anti-trapping term in Eq. (29) can be re-written as

$a(\phi)=\frac{(M_{s}-M_{l})h(\phi)\\{1-h(\phi)\\}}{M_{s}\\{1-h(\phi)\\}+M_{l}h(\phi)}.$
(31)

It is worth noting since $h(\phi)$ is used to interpolate the phase
concentrations in Eq. (1), the interface stretching term that modifies mass
conservation is automatically negated. Thus, all the thin-interface effects
can be eliminated by a single $h(\phi)$ function.

###  II.2 Asymptotic analysis and determining the source terms $S_{s}$ and
$S_{l}$

Typically the source terms should be determined by performing an asymptotic
analysis of the model. However, thin-interface asymptotic analysis (i.e. up to
second order) of the model remains a mathematically challenging task. To make
the problem mathematically amenable and for the purpose of estimating the
choice of $S_{s}$ and $S_{l}$, we temporarily assume pointwise equality of the
diffusion potential i.e. $\mu_{s}=\mu_{l}$ in Eqs. (II), (II), (II) and (II).
This allows for substantial progress in the thin-interface asymptotics, and a
choice of source term can be deduced. The asymptotic analysis builds upon and
follows closely the expansion scheme of Refs. [29, 45, 52, 58]. The lowest
order solutions are identical to the previous works and differences arise only
at higher order expressions, namely, first order jump in diffusion potential,
first order correction to Gibbs Thomson condition and first order correction
to Stefan’s condition. We restrict the discussion only to these new findings
here and the details of the asymptotic analyses (up to first order for the
model $\mu_{s}\neq\mu_{l}$ and up to second order for the model
$\mu_{s}=\mu_{l}$) are presented in detail in the supplementary materials
along with numerical comparisons of the solution obtained for both the cases
(to ascertain the range of the validity of the assumption). We emphasize that
$\mu_{s}=\mu_{l}$ is not a constraint when solving the dynamic equations.

The diffusion potential jump up to first order in the asymptotic expansion can
be obtained as

$\displaystyle\mu^{-}-\mu^{+}=$ $\displaystyle\frac{WV}{D_{l}\partial
c_{l,0}/\partial\mu_{0}}(c_{l,0}-c_{s,0})(F_{1}^{-}-F_{1}^{+})+\frac{WV}{D_{l}\partial
c_{l,0}/\partial\mu_{0}}(F_{2}^{-}-F_{2}^{+})$
$\displaystyle+\frac{W}{D_{l}\partial
c_{l,0}/\partial\mu_{0}}q_{n}^{-}\frac{\partial\mu_{0}}{\partial
r}\Big{|}^{-}(G^{-}-G^{+}),$ (32)

where,

$F_{1}^{\pm}=\int_{0}^{\pm\infty}\Big{[}\overline{p}(\phi_{0})-\overline{p}(\phi_{0}|^{\pm})\Big{]}\mathrm{d}\eta,$
(33)
$\overline{p}(\phi_{0})=\frac{g(\phi_{0})-1-a(\phi_{0})}{q_{n}(\phi_{0})/M_{l,0}},$
(34)
$F_{2}^{\pm}=\int_{0}^{\pm\infty}\Big{[}S_{s}g(\phi_{0})+S_{l}\\{1-g(\phi_{0})\\}\Big{]}\frac{a(\phi_{0})}{q_{n}(\phi_{0})/M_{l,0}}$
(35)

and

$G^{\pm}=M_{l,0}\int_{0}^{\pm\infty}\Bigg{[}\frac{1}{q_{n}(\phi_{0})}-\frac{1}{q_{n}^{\pm}}\Bigg{]}\mathrm{d}\eta.$
(36)

$M_{l,0}$ is the lowest order liquid mobility defined as
$M_{l,0}=D_{l}\partial c_{l,0}/\partial\mu$.

The superscripts $+(-)$ denote that the quantities have been evaluated in the
liquid (solid) side of the interface. The subscript $0$ represent the zeroth
(lowest) order solution and $\phi_{0}$ represents the equilibrium phase field
profile which is of the form of hyperbolic tangent and $\eta$ denotes the co-
ordinate system normal to the interface.

The source term in the phase field model provides an additional degree of
freedom to recover a desired sharp-interface relation. From the analysis it is
evident from Eq. (II.2) that the choice of source terms
$S_{s}=S_{l}=-A(c_{s}-c_{l})$ is sufficient to impart a degree of control over
solute trapping in the phase field model. With the above choice, the diffusion
potential jump up to first order in asymptotics can be simplified as,

$\mu^{-}-\mu^{+}=\frac{WV}{D_{l}\partial
c_{l,0}/\partial\mu}(c_{l,0}-c_{s,0})(F^{-}-F^{+})+\frac{W}{D_{l}\partial
c_{l,0}/\partial\mu}q_{n}^{-}\frac{\partial\mu_{s}}{\partial
r}\Big{|}^{-}(G^{-}-G^{+}),$ (37)

where,

$F^{\pm}=\int_{0}^{\pm\infty}\Big{[}\tilde{p}(\phi_{0})-\tilde{p}(\phi_{0}|^{\pm})\Big{]}\mathrm{d}\eta,$
(38)
$\tilde{p}(\phi_{0})=\frac{g(\phi_{0})-1-(1-A)a(\phi_{0})}{q_{n}(\phi_{0})/M_{l,0}}$
(39)

The second term in the diffusion potential jump Eq. (37) is associated with a
stationary interface and arises due to solute flux across the interface. Since
we are unaware of any computational or experimental evidence of the importance
of this term in rapid solidification processing, we choose to eliminate the
term by satisfying the condition $G^{-}=G^{+}$. The inverse interpolation
function given by Eq. (30) naturally satisfies the constraint. A numerical
comparison of the two models (shown in the supplementary material) suggest
that the steady state diffusion potential in both the models are identical.
The agreement between the predictions of the two models is better (within
$10\%$) at higher values of $A$ for the entire velocity range considered. Even
with lower values of $A$ the agreement is good until velocities of about
$0.1\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$. A good agreement between
both the models motivates us to utilize the same interpolation functions to
negate the Kapitza jump for the case $\mu_{s}\neq\mu_{l}$.

Once the second term in Eq. (37) is eliminated, solute trapping can be
controlled by the first term. The first term is the diffusion potential jump
due to interface velocity and the origin of solute trapping in rapid
solidification. In the limit where the asymptotics has been performed,
$\mu_{s}=\mu_{l}$, the magnitude of this jump can be controlled by choosing an
appropriate value of $A$ as evident from Eqs. (38) and (39). We find this to
be the case as well when $\mu_{s}\neq\mu_{l}$ and by choosing this form of the
source term and verifying numerically that trapping can indeed be controlled
by the parameter $A$. While other choices of source terms are certainly
possible, we choose the above form for the sake of simplicity. Moreover, such
a simple expression for the sources considerably simplifies the subsequent
asymptotic analysis. Thus the terms involving $\tilde{c}_{s}-\tilde{c}_{l}$
can be replaced by $(1-A)(c_{s}-c_{l})$ in Eqs. (II), (II) and (II).

Once the choice of source terms is determined, it is easy to evaluate the
first order correction to the Gibbs Thomson condition which is dependent on
the form of the sources. The Gibbs Thomson condition up to first order in
asymptotics is obtained as

$\displaystyle\mu^{\pm}=$
$\displaystyle\frac{f_{s}(c_{s,0})-f_{l}(c_{l,0})}{c_{s,0}-c_{l,0}}+\frac{\gamma\kappa}{c_{s,0}-c_{l,0}}$
$\displaystyle+\frac{V}{(c_{s,0}-c_{l,0})}\frac{a_{1}}{M_{\phi}W}\Big{[}1-a_{2}^{\pm}\frac{W^{2}M_{\phi}}{D_{l}}\frac{(c_{s,0}-c_{l,0})^{2}}{\partial
c_{l,0}/\partial\mu}\Big{]}$ $\displaystyle+\frac{W}{D_{l}\partial
c_{l,0}/\partial\mu}q_{n}^{-}\frac{\partial\mu_{0}}{\partial
r}\Big{|}^{-}A\Delta F$ (40)

where,

$a_{1}=I=\int_{-\infty}^{+\infty}\left(\frac{\partial\phi_{0}}{\partial\eta}\right)^{2}\mathrm{d}\eta,$
(41) $a_{2}^{\pm}=\frac{F^{\pm}+M-(1-A)^{2}N+(1-A)P}{I},$ (42)

with the constants defined as

$M=\int_{-\infty}^{+\infty}\left[\int_{0}^{\eta}\tilde{p}(\phi_{0})\mathrm{d}\eta^{\prime}\right]g^{\prime}(\phi_{0})\frac{\partial\phi_{0}}{\partial\eta}\mathrm{d}\eta,$
(43)
$N=M_{l,0}\int_{-\infty}^{+\infty}\frac{1}{M_{s}^{n}(\phi_{0})+M_{l}^{n}(\phi_{0})}\mathrm{d}\eta,$
(44) $P=\int_{-\infty}^{+\infty}a(\phi_{0})\tilde{p}(\phi_{0})\mathrm{d}\eta,$
(45)

and

$\Delta F=F^{-}-F^{+}.$ (46)

The value of the first order diffusion potential is different on either side
of the interface due to difference in the values of the constants $F^{+}$ and
$F^{-}$ in the expression of $a_{2}$. The correction of the mass conservation
equation due to finite interface width can be obtained up to first order as

$\displaystyle q_{n}^{+}\frac{\partial\mu}{\partial
r}\Big{|}^{+}-q_{n}^{-}\frac{\partial\mu}{\partial r}\Big{|}^{-}=$
$\displaystyle-V(c_{l}-c_{s})-VW\kappa\left(H^{+}-H^{-}\right)$
$\displaystyle-W\frac{\partial^{2}\mu_{0}}{\partial s^{2}}D_{l}\frac{\partial
c_{l,0}}{\partial\mu}\left(S^{+}-S^{-}\right)$ (47)

where,

$H^{\pm}=\int_{0}^{\pm\infty}\left[g(\phi_{0})-g(\phi_{0}|^{\pm})\right]\mathrm{d}\eta$
(48)

and

$S^{\pm}=\frac{1}{M_{l,0}}\int_{0}^{\pm\infty}\left[q_{t}(\phi_{0})-q_{t}(\phi_{0}|^{\pm})\right]\mathrm{d}\eta,$
(49)

with

$q_{t}(\phi_{0})=M_{s}^{t}(\phi_{0})g(\phi_{0})^{2}+M_{l}^{t}(\phi_{0})\left\\{1-g(\phi_{0})\right\\}^{2}$
(50)

The term $(c_{l}-c_{s})$ on the right hand side of the above equation denotes
jump in the concentration at the interface up to first order. The second and
the third terms in Eq. (II.2) represents the correction due to interface
stretching and surface diffusion respectively. From Eqs. (49) and (50) it is
evident that anti-symmetric forms of $g(\phi_{0})$ and $q_{t}(\phi_{0})$ will
satisy the constraints $H^{+}=H^{-}$ and $S^{+}=S^{-}$. Therefore,
$M_{s}^{t}(\phi_{0})$ and $M_{l}^{t}(\phi_{0})$ are chosen to give rise to
anti-symmetric form of $q_{t}(\phi_{0})$. The choice of the interpolation
functions $g(\phi)=h(\phi)$, $M_{s}^{t}(\phi)=M_{s}/h(\phi)$ and
$M_{l}^{t}/\left\\{1-h(\phi)\right\\}$ discussed in the Section II.1 naturally
satisfies the above constraints.

###  II.3 Alternative choice of interpolation functions

All the thin-interface effects can be controlled or eliminated within the
variational framework. However, the choice of single $h(\phi)$ function to
interpolate mobilities to remove both the Kapitza jump and surface diffusion
imposes a severe restriction on the numerical efficiency of the model. For
instance, the inverse interpolation function Eq. (30) leads to an abrupt
transition between the bulk values across the interface as shown in Fig.1. The
interface diffusivity is governed by the solid diffusivity. This restricts the
choice of the interface width $W$ in the simulations which has to be
progressively decreased with decreasing solid diffusivity to resolve the
diffusion length in the solid. Therefore, we utilize additional degrees of
freedom in the model by exploiting different interpolation functions for the
mobilities in the normal and tangential direction to control the different
thin-interface effects, while still retaining the variational character of the
equations.

The problem related to the abrupt variation of mobility in inverse
interpolation function can be ameliorated by interpolating the mobilities in
$q_{n}(\phi)$ in Eq. (28) as $M_{s}^{n}(\phi)=M_{s}/h_{n}(\phi)$ and
$M_{l}^{n}(\phi)=M_{l}/\\{1-h_{n}(\phi)\\}$ to allow for a smoother variation
of $q_{n}(\phi)$ as

$\frac{1}{q_{n}(\phi)}=\frac{h_{n}(\phi)}{M_{s}}+\frac{1-h_{n}(\phi)}{M_{l}},$
(51)

where

$h_{n}(\phi)=\frac{bp(\phi)}{bp(\phi)+1-p(\phi)},$ (52)

where, $b$ is a constant which determines the variation of the mobility across
the interface and $p(\phi)$ is given by

$p(\phi)=h(\phi)+a(\phi-\frac{1}{2})\phi^{2}(1-\phi)^{2}.$ (53)

$p(\phi)$ comprises of a summation of two anti-symmetric functions with an
adjustable parameter $a$. The constant $a$ is evaluated numerically such that
the integrals $G^{-}$ and $G^{+}$ in Eq. (36) are equal. It can be verified
that the function satisfies the constraints $p(\phi=0)=p(\phi=1)=0$ and
$p^{\prime}(\phi=0,1)=0$. Usually multiple values of $a$ satisfy the
constraint. The variation of $q_{n}(\phi)$ for different values of $a$ are
shown in Fig.2 for $b=0.001$ and $D_{s}/D_{l}=0.001$. As $a$ increases the
variation of $q_{n}(\phi)$ becomes non-monotonic. Therefore, we choose from
the smaller roots of $a$ for the study. The effect of the parameter $b$ on
variation of $q_{n}(\phi)$ is shown in Fig.2(b). Values of $b$ not close to
$D_{s}/D_{l}$ gives rise to asymmetric profiles, although a symmetric profile
is not a strict requirement. However, the choice of $b$ and $a$ should lead to
$q_{n}(\phi)\geq 0$ for the phase concentration equations to be stable.

To eliminate surface diffusion, we choose $M_{s}^{t}(\phi)=M_{s}/p(\phi)$ and
$M_{l}^{t}=M_{l}/\\{1-p(\phi)\\}$ to interpolate the tangential component of
the diffusion potentials since $p(\phi)$ is anti-symmetric. In addition we
choose $g(\phi)=p(\phi)$ to interpolate the phase concentrations (Eq. (1)),
free energies (Eq. (3)) and fluxes (Eq. (7)). The modified form of $a(\phi)$
(Eq. (29)) is written as

$a(\phi)=\frac{(M_{s}-bM_{l})p(\phi)\\{1-p(\phi)\\}}{M_{s}\\{1-p(\phi)\\}+M_{l}bp(\phi)}.$
(54)

The phase concentrations can thus be written as

$\displaystyle\frac{\partial c_{s}}{\partial t}=\nabla\cdot\Bigg{[}$
$\displaystyle q_{n}(\phi)\nabla_{n}\Bigg{(}\frac{\partial f_{s}}{\partial
c_{s}}p(\phi)+\frac{\partial f_{l}}{\partial c_{l}}\\{1-p(\phi)\\}\Bigg{)}$
$\displaystyle+M_{s}\nabla_{t}\frac{\partial f_{s}}{\partial
c_{s}}+\frac{({c_{s}}-{c_{l}})(1-A)}{|\nabla\phi|}\frac{M_{s}\\{1-p(\phi)\\}}{\overline{M}(\phi)}\frac{\nabla\phi}{|\nabla\phi|}\frac{\partial\phi}{\partial
t}\Bigg{]}$
$\displaystyle-A(c_{s}-c_{l})p^{\prime}(\phi)\frac{\partial\phi}{\partial t},$
(55)

$\displaystyle\frac{\partial c_{l}}{\partial t}=\nabla\cdot\Bigg{[}$
$\displaystyle q_{n}(\phi)\nabla_{n}\Bigg{(}\frac{\partial f_{s}}{\partial
c_{s}}p(\phi)+\frac{\partial f_{l}}{\partial c_{l}}\\{1-p(\phi)\\}\Bigg{)}$
$\displaystyle+M_{l}\nabla_{t}\frac{\partial f_{l}}{\partial
c_{l}}-\frac{({c_{s}}-{c_{l}})(1-A)}{|\nabla\phi|}\frac{bM_{l}p(\phi)}{\overline{M}(\phi)}\frac{\nabla\phi}{|\nabla\phi|}\frac{\partial\phi}{\partial
t}\Bigg{]}$
$\displaystyle-A(c_{s}-c_{l})p^{\prime}(\phi)\frac{\partial\phi}{\partial t}$
(56)

where, the terms in the square bracket in Eqs. (II.3) and (II.3) are the solid
($J_{s}$) and liquid fluxes ($J_{l}$) respectively. The overall solute
concentration is governed by the equation

$\displaystyle\frac{\partial c}{\partial t}=\nabla\cdot\Bigg{[}$
$\displaystyle q_{n}(\phi)\nabla_{n}\Bigg{(}\frac{\partial f_{s}}{\partial
c_{s}}p(\phi)+\frac{\partial f_{l}}{\partial c_{l}}\\{1-p(\phi)\\}\Bigg{)}$
$\displaystyle+M_{s}^{t}p(\phi)\nabla_{t}\frac{\partial f_{s}}{\partial
c_{s}}+M_{l}^{t}\\{1-p(\phi)\\}\nabla_{t}\frac{\partial f_{l}}{\partial
c_{l}}$
$\displaystyle+\frac{({c_{s}}-{c_{l}})}{|\nabla\phi|}(1-A)a(\phi)\frac{\nabla\phi}{|\nabla\phi|}\frac{\partial\phi}{\partial
t}\Bigg{]},$ (57)

Similarly, the phase field equation is written as

$\displaystyle\Bigg{[}\frac{1}{M_{\phi}}+$
$\displaystyle\frac{(c_{s}-c_{l})^{2}(1-A)^{2}bp(\phi)\\{1-p(\phi)\\}}{\\{\overline{M}(\phi)(bp(\phi)+1-p(\phi))\\}}\frac{1}{|\nabla\phi|^{2}}\Bigg{]}\frac{\partial\phi}{\partial
t}=\sigma\nabla^{2}\phi-Hf_{\mathrm{dw}}^{\prime}(\phi)-\Delta
G_{\mathrm{driv}}p^{\prime}(\phi)$
$\displaystyle+\frac{(c_{s}-c_{l})(1-A)p(\phi)\\{1-p(\phi)\\}}{\overline{M}(\phi)}\frac{1}{|\nabla\phi|}\Bigg{(}M_{s}\nabla\frac{\partial
f_{s}}{\partial c_{s}}-bM_{l}\nabla\frac{\partial f_{l}}{\partial
c_{l}}\Bigg{)}\cdot\frac{\nabla\phi}{|\nabla\phi|},$ (58)

where,

$\displaystyle\Delta
G_{\mathrm{driv}}=f_{s}(c_{s})-f_{l}(c_{l})-\mu_{l}(c_{s}-c_{l})-$
$\displaystyle(c_{s}-c_{l})(\mu_{s}-\mu_{l})\times$
$\displaystyle\Big{\\{}\frac{M_{s}\\{1-p(\phi)\\}(A\\{1-p(\phi)\\}-1)-M_{l}p(\phi)Ap(\phi)}{\overline{M}(\phi)}\Big{\\}}$
(59)

with $\overline{M}(\phi)=M_{s}\\{1-g(\phi)\\}+bM_{l}g(\phi)$.

It can be verified that the asymptotic expressions given by Eqs. (37), (38),
(39) and (36) remain unchanged but with new definitions of $q_{n}(\phi_{0})$
and $a(\phi_{0})$ given by Eqs. (51) and (54). Owing to the numerical
efficiency, in the rest of the article we will employ Eqs. (II.3), (II.3) and
(II.3) for numerical purposes.

The expression of the driving force Eq. (II.3) warrants a discussion to make
connection with the sharp-interface picture.The driving force in the classical
expression for free energy dissipation from the CGM [16] can be written as

$\Delta G_{\mathrm{driv}}=\Delta G_{c}-\mathcal{D}\Delta G_{D}$ (60)

where, $\Delta G_{c}$ is the free energy of crystallization or the free energy
change when solid of composition $c_{s}$ is formed from liquid of composition
$c_{l}$. Thus,

$\Delta G_{c}=f_{s}-f_{l}-\mu_{l}(c_{s}-c_{l}).$ (61)

$\Delta G_{D}$ is the solute-drag free energy and is defined as the amount of
energy dissipated in solute redistribution, given by,

$\Delta G_{D}=(c_{s}-c_{l})(\mu_{s}-\mu_{l}).$ (62)

Substituting Eqs. (61) and (62) in Eq. (60) yields

$\Delta
G_{\mathrm{driv}}=f_{s}-f_{l}-\mu_{l}(c_{s}-c_{l})-\mathcal{D}(c_{s}-c_{l})(\mu_{s}-\mu_{l}).$
(63)

Comparing Eqs. (II.3) and (63) it can be seen that the drag coefficient in
this phase field model can be written as

$\mathcal{D}(\phi)=\frac{M_{s}\\{1-p(\phi)\\}(A\\{1-p(\phi)\\}-1)-M_{l}p(\phi)Ap(\phi)}{M_{s}\\{1-p(\phi)\\}+M_{l}p(\phi)}.$
(64)

It is evident that solute drag in the present model is dependent on the value
of the model parameter $A$. Thus, $A$ not only controls solute trapping but
also solute drag.

Figure 1:  The variation of the diffusivity interpolation function
$q_{n}(\phi)$ given by Eq. (30) for different ratios of solid to liquid
diffusion coefficients. The profiles exhibit a sharp transition as $\phi$
approaches $0$ (i.e. liquid). As a result the interface diffusivity
corresponding to $\phi=0.5$ is dictated by the diffusion coefficient in the
solid.  Figure 2: (a) The variation of the diffusivity interpolation function
$q_{n}(\phi)$ given by Eq. (51) for different values of the roots of $a$ with
$b=0.001$ and $D_{s}/D_{l}=0.001$ such that the constraint $G^{+}=G^{-}$ is
satisfied. Higher values of the root leads to non-monotonic variation across
the interface. (b) The variation of the $q_{n}(\phi)$ across the interface for
different values of $b$ and $D_{s}/D_{l}=0.001$. The quintic interpolation
function $h(\phi)=\phi^{3}(6\phi^{2}-15\phi+10)$ has been employed in the
plots.

##  III Dilute alloy approximation

To perform numerical computation we have to specify the thermodynamics of the
system. We specialize the equations by considering a dilute ideal solution
alloy. The free energy of the bulk solid and liquid is written as

$f_{i}(c_{i})=f_{i}^{\mathrm{A}}(T)+\epsilon_{i}c_{i}+\frac{RT}{V_{\mathrm{m}}}(c_{i}\ln{}c_{i}-c_{i})\hskip
56.9055pti=s,l$ (65)

where $f_{i}^{\mathrm{A}}$ is the free energy of the pure A in the respective
phase and the second and the third term represent the change of the internal
energy density and the dilute form of entropy of mixing respectively. We
linearize the free energy about $(T-T_{\mathrm{m}})$ to recover straight
solidus and liquidus and retaining the first order terms to obtain

$f_{i}(c_{i})=f_{i}^{A}(T_{\mathrm{m}})-s_{i}(T-T_{\mathrm{m}})+\epsilon_{i}c_{i}+\frac{RT_{\mathrm{m}}}{V_{\mathrm{m}}}(c_{i}\ln
c_{i}-c_{i})$ (66)

where $s_{i}=-\partial f_{i}^{\mathrm{A}}/\partial T$ are the entropy
densities.

At equilibrium, from the condition of equality of diffusion potential it is
easy to show

$k_{e}=\frac{c_{s}^{\mathrm{eq}}}{c_{l}^{\mathrm{eq}}}=\mathrm{exp}\Bigg{[}-\frac{\epsilon_{s}-\epsilon_{l}}{RT_{\mathrm{m}}/V_{\mathrm{m}}}\Bigg{]}.$
(67)

In addition, the equality of grand-potential i.e.

$\omega_{s}=f_{s}-\mu_{s}c_{s}=f_{l}-\mu_{l}c_{l}=\omega_{l}$ (68)

yields the relation

$f_{s}^{\mathrm{A}}(T)-f_{l}^{\mathrm{A}}(T)=\frac{RT_{\mathrm{m}}}{V_{\mathrm{m}}}(c_{s}^{\mathrm{eq}}-c_{l}^{\mathrm{eq}}).$
(69)

The driving forces in the phase field equation Eq. (II) are given by

$f_{s}-f_{l}-(c_{s}-c_{l})\frac{\partial f_{s}}{\partial
c_{s}}=\frac{RT_{\mathrm{m}}}{V_{\mathrm{m}}}c_{l}\ln\Big{(}\frac{1}{k_{e}}\frac{c_{s}}{c_{l}}\Big{)}+\frac{RT_{\mathrm{m}}}{V_{\mathrm{m}}}\Big{[}c_{l}^{\mathrm{eq}}(k_{e}-1)-(c_{s}-c_{l})\Big{]}$
(70)

and

$f_{s}-f_{l}-(c_{s}-c_{l})\frac{\partial f_{l}}{\partial
c_{l}}=\frac{RT_{\mathrm{m}}}{V_{\mathrm{m}}}c_{s}\ln\Big{(}\frac{1}{k_{e}}\frac{c_{s}}{c_{l}}\Big{)}+\frac{RT_{\mathrm{m}}}{V_{\mathrm{m}}}\Big{[}c_{l}^{\mathrm{eq}}(k_{e}-1)-(c_{s}-c_{l})\Big{]}$
(71)

The driving force can also be expressed in terms of undercooling. For a
linearized phase diagram the straight solidus and liquidus lines have slopes
$m/k_{e}$ and $m$ respectively. The equilibrium composition $c_{l}^{e}$ can be
eliminated in favor of the temperature $T$ as

$T=T_{\mathrm{m}}-|m|c_{l}^{\mathrm{eq}},$ (72)

where $T_{\mathrm{m}}$ is the melting temperature of pure $B$. Therefore,

$f_{s}-f_{l}-(c_{s}-c_{l})\frac{\partial f_{s}}{\partial
c_{s}}=\frac{RT_{\mathrm{m}}}{V_{\mathrm{m}}}c_{l}\Big{(}\frac{1}{k_{e}}\frac{c_{s}}{c_{l}}\Big{)}+\frac{RT_{\mathrm{m}}}{V_{\mathrm{m}}}\Big{[}\frac{(T_{\mathrm{m}}-T)(k_{e}-1)}{|m|}-(c_{s}-c_{l})\Big{]}.$
(73)

Similarly,

$f_{s}-f_{l}-(c_{s}-c_{l})\frac{\partial f_{l}}{\partial
c_{l}}=\frac{RT_{\mathrm{m}}}{V_{\mathrm{m}}}c_{s}\Big{(}\frac{1}{k_{e}}\frac{c_{s}}{c_{l}}\Big{)}+\frac{RT_{\mathrm{m}}}{V_{\mathrm{m}}}\Big{[}\frac{(T_{\mathrm{m}}-T)(k_{e}-1)}{|m|}-(c_{s}-c_{l})\Big{]}.$
(74)

###  III.1 Directional solidification

In experiments, it is difficult to maintain a fixed undercooling and far
easier to control the interface speed through an arrangement known as the
directional solidification (DS). In a DS set up, a temperature gradient is
established across the sample, which is then pulled through the temperature
field at constant velocity. The two variables that control the microstructure
features, temperature gradient $G$ and pulling speed $v_{p}$ are controlled
independently. In addition, a frozen temperature approximation is employed
which assumes that the temperature field is unperturbed by the advancing
interface such that

$T(x,t)=T_{\mathrm{ref}}+G[(x-x_{o})-v_{p}t].$ (75)

$T_{\mathrm{ref}}$ is the reference temperature at the initial position
$x_{o}$ of the solid-liquid interface. It is customary to measure the local
temperature and concentration with respect to a reference temperature and
concentration. For a sample with composition $c_{\infty}$, the solidus
temperature $T_{s}$ and the corresponding alloy composition $c_{\infty}$ is
selected to be the reference temperature and concentration respectively. The
driving forces Eqs. (73) and (74) can be re-written with respect to
($c_{\infty},T_{s}$) are

$f_{s}-f_{l}-(c_{s}-c_{l})\frac{\partial f_{s}}{\partial
c_{s}}=\frac{RT_{\mathrm{m}}}{V_{\mathrm{m}}}c_{l}\Big{(}\frac{1}{k_{e}}\frac{c_{s}}{c_{l}}\Big{)}+\frac{RT_{\mathrm{m}}}{V_{\mathrm{m}}}\Bigg{[}\frac{(T-T_{s})(1-k_{e})}{|m|}+\Big{\\{}(c_{l}-c_{s})-\frac{c_{\infty}}{k_{e}}(1-k_{e})\Big{\\}}\Bigg{]}$
(76)

and

$f_{s}-f_{l}-(c_{s}-c_{l})\frac{\partial f_{l}}{\partial
c_{l}}=\frac{RT_{\mathrm{m}}}{V_{\mathrm{m}}}c_{s}\Big{(}\frac{1}{k_{e}}\frac{c_{s}}{c_{l}}\Big{)}+\frac{RT_{\mathrm{m}}}{V_{\mathrm{m}}}\Bigg{[}\frac{(T-T_{s})(1-k_{e})}{|m|}+\Big{\\{}(c_{l}-c_{s})-\frac{c_{\infty}}{k_{e}}(1-k_{e})\Big{\\}}\Bigg{]}$
(77)

##  IV Relation to other phase field models

In this section we show that under certain conditions a number of previously
developed phase field models can be derived as special cases of the present
model.

###  IV.1 Wheeler-Boettinger-McFadden (WBM) model

The WBM model [35] was the first phase field model developed to study
isothermal binary alloy solidification. In the WBM model, the interface is
defined as a mixture of solid and liquid of same composition i.e.
$c_{s}=c_{l}=c$. The free energy is thus defined as

$f(c,\phi)=f_{s}(c)g(\phi)+f_{l}(c)\\{1-g(\phi)\\}.$ (78)

The constraint Eq. (2) and decomposition to kinetic Eqs. (9) and (II) is not
required. The governing equations for $c$ and $\phi$ can be derived to
guarantee monotonic decrease of the free energy functional Eq. (5) as

$\frac{\partial c}{\partial
t}=\nabla\cdot\Big{[}q_{n}(\phi)\nabla\frac{\delta\mathcal{F}}{\delta
c}\Big{]},$ (79)

and

$-\frac{1}{M_{\phi}}\frac{\partial\phi}{\partial
t}=\frac{\delta\mathcal{F}}{\delta\phi}$ (80)

The diffusion potential $\mu=\frac{\delta\mathcal{F}}{\delta c}$ is defined as

$\mu=\frac{\partial f_{s}(c)}{\partial c}g(\phi)+\frac{\partial
f_{l}(c)}{\partial c}\\{1-g(\phi)\\}$ (81)

Alternatively, we can set $c_{s}=c_{l}=c$ in Eq. (II) to recover the phase
field kinetic equation of the WBM model as

$\frac{1}{M_{\phi}}\frac{\partial\phi}{\partial t}=\sigma\nabla^{2}\phi-
Hf_{\mathrm{dw}}^{\prime}(\phi)-\Big{\\{}f_{s}(c)-f_{l}(c)\Big{\\}}g^{\prime}(\phi).$
(82)

A direct reduction of the concentration equation to the WBM model is obtained
by substituting $c_{s}=c_{l}=c$,
$M_{s}^{n}(\phi)=M_{s}^{t}(\phi)=M_{s}^{n}/h_{n}(\phi)$,
$M_{l}^{n}(\phi)=M_{l}^{t}(\phi)=M_{s}/\\{1-h_{n}(\phi)\\}$, with
$h_{n}(\phi)=M_{s}g(\phi)/\\{M_{s}g(\phi)+M_{l}(1-g(\phi))\\}$ in Eq. (II) to
yield

$\frac{\partial c}{\partial
t}=\nabla\cdot\Big{[}q(\phi)\nabla\Big{\\{}\frac{\partial f_{s}(c)}{\partial
c}g(\phi)+\frac{\partial f_{l}(c)}{\partial
c}\\{1-g(\phi)\\}\Big{\\}}\Big{]},$ (83)

where $q(\phi)=M_{s}g(\phi)+M_{l}\\{1-g(\phi)\\}$. Thus with the choice of
$g(\phi)=h(\phi)$ the WBM model can be recovered.

###  IV.2 Kim-Kim-Suzuki (KKS) model

The KKS model [36] like the present model assumes the interface to be a
mixture of solid and liquid with distinct compositions $c_{s}$ and $c_{l}$.
However the phase concentrations are related via the equality of diffusion
potentials i.e. $\frac{\partial f_{s}}{\partial c_{s}}=\frac{\partial
f_{l}}{\partial c_{l}}$ at each spatial position in the interface. The kinetic
equation for the KKS model can be recovered by substituting $\frac{\partial
f_{s}}{\partial c_{s}}=\frac{\partial f_{l}}{\partial c_{l}}$ and setting the
source term $A=1$ in Eq. (II) to obtain

$\frac{1}{M_{\phi}}\frac{\partial\phi}{\partial t}=\sigma\nabla^{2}\phi-
Hf_{\mathrm{dw}}(\phi)-\Big{\\{}f_{s}(c_{s})-f_{l}(c_{l})-\mu_{l}(c_{s}-c_{l})\Big{\\}}g^{\prime}(\phi).$
(84)

The concentration equation for the KKS model can be recovered by additionally
assuming $M_{s}^{n}(\phi)=M_{s}^{t}(\phi)=M_{s}^{n}/h_{n}(\phi)$,
$M_{l}^{n}(\phi)=M_{l}^{t}(\phi)=M_{s}/\\{1-h_{n}(\phi)\\}$, with
$h_{n}(\phi)=M_{s}g(\phi)/\\{M_{s}g(\phi)+M_{l}(1-g(\phi))\\}$ in Eq. (II) and
letting $g(\phi)=h(\phi)$ to obtain

$\frac{\partial c}{\partial t}=\nabla\cdot\Big{[}q(\phi)\nabla\mu_{l}\Big{]}$
(85)

where $q(\phi)=M_{s}h(\phi)+M_{l}\\{1-h(\phi)\\}$. It is to be noted that the
equality of diffusion potential is solved at each spatial point and time
iteration to relate the phase concentrations. This additional step in KKS can
become a computation bottle-neck depending on the alloy model in multi-phase,
multi-component systems. In contrast in the present model, phase
concentrations are governed by separate kinetic equations and can handle
arbitrary phase diagrams with equal ease. One can thus expect to gain some
computational expedience especially for complex free energies. However, a
comparative study on the numerical performance of both the models is outside
the scope of the present work.

###  IV.3 Echebarria-Folch-Karma-Plapp (EFKP) model

EFKP [45] is a specialized form of KKS with dilute alloy approximation as
shown in Ref.[31]. To make connection with the EFKP model we define the solute
field in terms of local supersaturation
$U=(c_{l}-c_{l}^{\mathrm{eq}})/(c_{l}^{\mathrm{eq}}-c_{s}^{\mathrm{eq}})$.
Furthermore, the equality of diffusion potential for the dilute alloy
approximation relates the phase concentrations as $c_{s}/c_{l}=k_{e}$, where,
$k_{e}$ is the equilibrium partition coefficient. Assuming the source terms to
be $A=1$ Eqs. (84) and (85) can be re-written employing Eq. (76) as

$\frac{1}{M_{\phi}}\frac{\partial\phi}{\partial t}=\sigma\nabla^{2}\phi-
Hf_{\mathrm{dw}}^{\prime}(\phi)-\tilde{\lambda}\Big{[}U+\frac{T-T_{s}}{|m|c_{l}^{\mathrm{eq}}(1-k_{e})}\Big{]}h^{\prime}(\phi)$
(86)

where
$\tilde{\lambda}=\frac{RT_{\mathrm{m}}}{V_{\mathrm{m}}}c_{l}^{\mathrm{eq}}(1-k_{e})^{2}$
and

$\Big{[}1-h(\phi)(1-k_{e})\Big{]}\frac{\partial U}{\partial
t}=\nabla\cdot\Big{[}q(\phi)\nabla
U\Big{]}+[1+(1-k_{e})U]h^{\prime}(\phi)\frac{\partial\phi}{\partial t}.$ (87)

which are the variational form of the governing equations employed in the EFKP
model. The non-variational form of the model employs an anti-trapping current
that is added to Eq. (87) to maintain local equilibrium at the interface such
that the above Eq. is modified as

$\Big{[}1-h(\phi)(1-k_{e})\Big{]}\frac{\partial U}{\partial
t}=\nabla\cdot\Big{[}q(\phi)\nabla
U+\\{1+(1-k_{e})U\\}a(\phi)\frac{\nabla\phi}{|\nabla\phi|}\frac{\partial\phi}{\partial
t}\Big{]}+[1+(1-k_{e})U]h^{\prime}(\phi)\frac{\partial\phi}{\partial t}.$ (88)

The form of $a(\phi)$ is deduced based on the thin-interface analysis of the
model. The model by Pinomaa and Provatas [44] tune the form of $a(\phi)$ to
control trapping.

###  IV.4 Non-diagonal model of Ohno et al.

One can derive the non-diagonal model of Ohno et al. [52] which was employed
to study equilibrium solidification processes by setting $\mu_{s}=\mu_{l}$ and
choosing the source terms to be zero i.e. $S_{s}=S_{l}=0$ (i.e. $A=0$) with
the choice of $M_{s}^{n}(\phi)=M_{s}^{t}(\phi)=M_{s}/h(\phi)$,
$M_{l}^{n}(\phi)=M_{l}^{t}(\phi)=M_{l}/\\{1-h(\phi)\\}$ and $g(\phi)=h(\phi)$
in Eqs. (II) and (II) to obtain

$\displaystyle\Bigg{[}\frac{1}{M_{\phi}}+$
$\displaystyle\frac{(c_{s}-c_{l})^{2}h(\phi)\\{1-h(\phi)\\}}{M_{s}\\{1-h(\phi)\\}+M_{l}h(\phi)}\frac{1}{|\nabla\phi|^{2}}\Bigg{]}\frac{\partial\phi}{\partial
t}=\sigma\nabla^{2}\phi-
Hf_{\mathrm{dw}}^{\prime}(\phi)-\Big{[}f_{s}(c_{s})-f_{l}(c_{l})$
$\displaystyle-\mu_{l}(c_{s}-c_{l})\Big{]}h^{\prime}(\phi)+\frac{(c_{s}-c_{l})}{|\nabla\phi|}a(\phi)\nabla\mu_{l}\cdot\frac{\nabla\phi}{|\nabla\phi|},$
(89)

and,

$\displaystyle\frac{\partial c}{\partial t}=$
$\displaystyle\nabla\cdot\Bigg{[}q_{n}(\phi)\nabla_{n}\mu_{l}+q_{t}(\phi)\nabla_{t}\mu_{l}+\frac{({c_{s}}-{c_{l}})}{|\nabla\phi|}a(\phi)\frac{\nabla\phi}{|\nabla\phi|}\frac{\partial\phi}{\partial
t}\Bigg{]},$ (90)

where,

$q_{n}(\phi)=\frac{M_{s}M_{l}}{M_{s}\\{1-h(\phi)\\}+M_{l}h(\phi)},$ (91)
$q_{t}(\phi)=M_{s}h(\phi)+M_{l}\\{1-h(\phi)\\}$ (92)

and

$a(\phi)=\frac{(M_{s}-M_{l})h(\phi)\\{1-h(\phi)\\}}{M_{s}\\{1-h(\phi)\\}+M_{l}h(\phi)}$
(93)

Trapping is present in the model at high velocities but cannot be controlled.
The present model tunes solute trapping by the judicious choice of the
interfacial source terms.

###  IV.5 Finite Interface Dissipation (FID) Model

The finite interface dissipation model by Steinbach et al. [43] similar to the
present work treat the phase concentrations $c_{s}$ and $c_{l}$ as independent
variables. While the FID model cannot be derived as a special case from the
present model, we seek to highlight the difference and similarities in the
structure of both the models. The starting point of both the present and the
FID models are Eqs. (1)-(5). The subsequent derivation of FID follows under
three assumptions: (i) the spatial domain is divided into distinct reference
volumes (RV), (ii) each RV is insulated from outside and (iii) processes
within a RV and exchange of atoms between adjacent RVs are independent and can
be superimposed. The evolution of the phase concentrations inside an RV are
chosen according to Eq. (5) as

$\frac{\partial c_{s}}{\partial t}=-P_{s}(\phi)\left[\frac{\delta F}{\delta
c_{s}}-\lambda g(\phi)\right]$ (94) $\frac{\partial c_{l}}{\partial
t}=-P_{l}(\phi)\left[\frac{\delta F}{\delta
c_{l}}-\lambda\\{1-g(\phi)\\}\right],$ (95)

where, $P_{s}(\phi)$ and $P_{l}(\phi)$ are the kinetic coefficients termed as
interface permeability for redistribution fluxes and chosen to be of the form
$P_{s}(\phi)=P/g(\phi)$ and $P_{l}(\phi)=P/\\{1-g(\phi)\\}$. Based on
assumption (ii), $\lambda$ is evaluated from the condition

$\frac{\partial c}{\partial t}=0.$ (96)

Substituting Eqs. (1), (94) and (95) in Eq. (96) we obtain

$\lambda=\frac{\partial f_{s}}{\partial c_{s}}g(\phi)+\frac{\partial
f_{l}}{\partial
c_{l}}\\{1-g(\phi)\\}-\frac{(c_{s}-c_{l})}{P}g^{\prime}(\phi)\frac{\partial\phi}{\partial
t}$ (97)

Substituting Eq. (97) in Eqs. (94) and (95) and adding the diffusion fluxes
between adjacent RVs the kinetic equations for phase concentrations can be
obtained as

$\displaystyle g(\phi)\frac{\partial c_{s}}{\partial t}=$
$\displaystyle\nabla\cdot\left(g(\phi)M_{s}\nabla\frac{\partial
f_{s}}{\partial c_{s}}\right)+Pg(\phi)\\{1-g(\phi)\\}\left(\frac{\partial
f_{l}}{\partial c_{l}}-\frac{\partial f_{s}}{\partial c_{s}}\right)$
$\displaystyle+(c_{l}-c_{s})g(\phi)g^{\prime}(\phi)\frac{\partial\phi}{\partial
t}$ (98) $\displaystyle\\{1-g(\phi)\\}\frac{\partial c_{l}}{\partial t}=$
$\displaystyle\nabla\cdot\left(\\{1-g(\phi)\\}M_{l}\nabla\frac{\partial
f_{l}}{\partial c_{l}}\right)+Pg(\phi)\\{1-g(\phi)\\}\left(\frac{\partial
f_{s}}{\partial c_{s}}-\frac{\partial f_{l}}{\partial c_{l}}\right)$
$\displaystyle+(c_{l}-c_{s})\\{1-g(\phi)\\}g^{\prime}(\phi)\frac{\partial\phi}{\partial
t}.$ (99)

It is worth noting that diffusion fluxes from only the respective phases
(first term in the above equations) are added during superposition. Adding
Eqs. (IV.5) and (IV.5) the overall solute conservation writes as

$\frac{\partial c}{\partial
t}=\nabla\cdot\left(g(\phi)M_{s}\nabla\frac{\partial f_{s}}{\partial
c_{s}}\right)+\nabla\cdot\left(\\{1-g(\phi)\\}M_{l}\nabla\frac{\partial
f_{l}}{\partial c_{l}}\right).$ (100)

The evolution equation for the phase-field follows from Eq. (5) as

$\frac{\partial\phi}{\partial t}=-M_{\phi}\left[\frac{\delta
F}{\delta\phi}-\lambda(c_{s}-c_{l})g^{\prime}(\phi)\right].$ (101)

Substituting Eq. (97) in the above equation yields the equation of motion as

$\displaystyle\left[\frac{1}{M_{\phi}}+\frac{(c_{s}-c_{l})^{2}}{P}g^{\prime}(\phi)^{2}\right]\frac{\partial\phi}{\partial
t}=\frac{\sigma}{2}\nabla^{2}\phi-Hf_{\mathrm{dw}}^{\prime}(\phi)$
$\displaystyle-\left[(f_{s}-f_{l})-\left\\{\frac{\partial f_{s}}{\partial
c_{s}}g(\phi)+\frac{\partial f_{l}}{\partial
c_{l}}\left(1-g(\phi)\right)\right\\}(c_{s}-c_{l})\right]g^{\prime}(\phi).$
(102)

Let us first compare the structure of the phase concentration equations. In
the FID model diffusion flux in the kinetic equations of the phase
concentrations are governed by gradient of the diffusion potential of the
respective phases (first term in Eqs. (IV.5) and (IV.5)). The source
expressions consist of terms driven by the difference between the diffusion
potentials of the two phases (second term) and the solute redistribution term
due to the moving interface or phase transformation (third term). This result
is borne out of the assumptions (i) and (iii). On the other hand, Eqs. (II.3)
and (II.3) in the present model suggest that the phase concentrations are
driven by the gradient of the diffusion potential of both the phases. This
difference arises because of the specific decomposition of the flux in the
interfacial region according to Eq. (7) and the subsequent derivation of the
form of the fluxes that guarantee free energy decrease. The sources in the
present work only account for solute redistribution due to phase change and
controls partitioning. The magnitude of the diffusion potential jump in FID
model is controlled by the combination of parameters $P$ and $W$. While the
spirit of gaining flexibility to tune solute trapping in both the models is
similar i.e. decomposition of independent kinetic equations for phase
concentrations leading to additional parameters namely source terms and
interface permeability, the connection of the FID model to the sharp interface
description is not as apparent as the present work. The authors presented an
expression to relate the diffusion potential jump to $P$ and $W$ but limited
the treatment to a stationary interface. This is also evident by careful
consideration of Eqs. (IV.5) and (IV.5) revealing that the magnitude of
diffusion potential jump at a stationary interface
($\frac{\partial\phi}{\partial t}=0$) is governed by $P$. It is similar to the
second term in Eq. (37) of the present analysis and arises due to mass flux
across the interface. Solute trapping during rapid solidification is a result
of suppressed solute partitioning due to high interface velocities. The source
terms in the present work have a direct connection to this phenomenon as
evident from the first term in Eq. (37). In the FID model, although solute
partitioning during rapid solidification is controlled by $P$ for enlarged
interface widths, it is difficult to separately control the two contributions
of the diffusion potential jump. A rigorous mathematical analysis of the
governing equations may help in this regard but is outside the scope of the
present work.

We next compare the equations of motion for the phase-field in both the
models. The derivation of both the models lead to a new mobility correction
term as evident on the left hand sides of Eqs. (II.3) and (IV.5). This is not
entirely fortuitous as the starting form of the phase-field equation arises
from the free energy functional Eq. (5). The differences in the rest of the
equations are due to the different assumptions made while evaluating
$\lambda$. The driving force in FID model (the terms in the square bracket on
the right hand side of Eq. (IV.5)) can be rearranged to obtain

$\Delta G_{\mathrm{driv}}=f_{s}-f_{l}-\frac{\partial f_{l}}{\partial
c_{l}}(c_{s}-c_{l})-\left(\frac{\partial f_{s}}{\partial c_{s}}-\frac{\partial
f_{l}}{\partial c_{l}}\right)(c_{s}-c_{l})g(\phi)$ (103)

which has a similar structure to the $\Delta G_{\mathrm{driv}}$ of the present
model as seen from Eq. (II.3). A comparison with the sharp interface
expression Eq. (63) suggests that the drag coefficient in the FID model is
given by $\mathcal{D}(\phi)=g(\phi)$. While the magnitude of the jump
$\left(\frac{\partial f_{s}}{\partial c_{s}}-\frac{\partial f_{l}}{\partial
c_{l}}\right)$ is controlled by $P$, its effect on solute drag is not apparent
from the above expression. In contrast in the present work, it is evident from
Eq. (64) that $A$ controls solute drag in addition to solute trapping.

Figure 3:  Effect of solute trapping parameter $A$ on the concentration
profiles with $v_{p}=0.32$ $\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ and
$W=13.5$ nm for (a) $D_{s}/D_{l}=0.1$ and $b=1.0$ and (b) $D_{s}/D_{l}=0.001$
and $b=0.00001$. In the parameter space where $M_{s}<bM_{l}$ is satisfied,
increasing the value of $A$ increases solute segregation, while the trend is
reversed when $M_{s}>bM_{l}$, due to the sign change in the expression for
$a(\phi)$ given by Eq. (54).

##  V Results and Discussion

###  V.1 Numerical technique

Table 1:  Thermophysical properties of Si–$9\text{\,}\mathrm{\char 37\relax}$As and Ni–$5\text{\,}\mathrm{\char 37\relax}$Cu and directional solidification parameters used in the present study. Parameter |  Si–$9\text{\,}\mathrm{\char 37\relax}$As [44] |  Ni–$5\text{\,}\mathrm{\char 37\relax}$Cu [59]
---|---|---
Equilibrium partition coefficient, $k_{e}$ | 0.3 | 0.5
Melting temperature, $T_{\mathrm{m}}$[$\text{\,}\mathrm{K}$] | 1685 | 1728
Equilibrium liquidus slope $m$[$\text{\,}\mathrm{K}\text{\,}{\mathrm{\char 37\relax}}^{-1}$] | -4.0 | -2.0
Alloy concentration, $c_{\mathrm{\infty}}$[$\text{\,}\mathrm{\char 37\relax}$] | 9 | 5
Gibbs Thomson coefficient, $\Gamma$[$\text{\,}\mathrm{K}\text{\,}\mathrm{m}$] | 3.4 $\times$ $10^{-7}$ | 2.6 $\times$ $10^{-7}$
Capillary length, $d_{o}$[$\text{\,}\mathrm{m}$] | 13.5 $\times$ $10^{-9}$ | 26 $\times$ $10^{-9}$
Liquid diffusion coefficient, $D_{l}$[$\text{\,}{\mathrm{m}}^{2}\text{\,}{\mathrm{s}}^{-1}$] | 15 $\times$ $10^{-9}$ | 3 $\times$ $10^{-9}$
Kinetic coefficient, $\beta$[$\text{\,}\mathrm{s}\text{\,}{\mathrm{m}}^{-1}$] | 0.595 | 0.2
Solid diffusion coefficient, $D_{s}$[$\text{\,}{\mathrm{m}}^{2}\text{\,}{\mathrm{s}}^{-1}$] |  $10^{-3}{D_{l}}$ to $0.1{D_{l}}$ |
Thermal gradient, $G$[$\text{\,}\mathrm{K}\text{\,}{\mathrm{m}}^{-1}$] | 1.4 $\times$ $10^{7}$ |
Pulling speed, $v_{p}$[$\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$] | 0.02 to 1 |

To perform numerical tests we choose the binary alloys
Si–$9\text{\,}\mathrm{\char 37\relax}$As and Ni–$5\text{\,}\mathrm{\char
37\relax}$Cu as model systems with compositions expressed as atomic
percentages. The thermodynamic and kinetic data along with the solidification
parameters are presented in Table 1. The governing Eqs. (II.3), (II.3) and
(II.3) are written in non-dimensional forms for numerical implementation. The
non-dimensionalization procedure is presented in the supplementary material.
The governing equations have been discretized using an explicit finite-
difference technique on a staggered-grid. All the spatial derivatives have
been resolved to second-order accuracy. The gradients and the divergence
operators in the diffusion equation are discretized using a combination of
forward and backward differences, with the mobility and the anti-trapping
function calculated at the the center of the two grid points. For temporal
discretization, we have employed the order Euler scheme. The grid spacing in
most simulations is selected to be $\Delta x=0.5W$, however, occasionally we
have also employed a coarser mesh with $\Delta x=0.8W$ which is the lowest
resolution for which the results are well converged. The simulation is
initiated with an equilibrium phase field profile and a homogeneous
equilibrium phase concentrations of $c_{l}=c_{\infty}$ and
$c_{s}=k_{e}c_{\infty}$ at each point. A no-flux boundary condition is
employed at the cold end (however, see the next paragraph, as this is mooted
by the shifting algorithm), while the concentration and phase field are kept
fixed to their respective values at the hot end. The temperature profile is
pulled with a velocity $v_{p}$ resulting in solute partitioning.

In addition, we have employed a shifting algorithm to reduce computation time.
This is implemented by shifting the solute and phase field profiles by an
amount equal to $\Delta x$, whenever the product of pulling speed and time
($v_{p}t$) corresponds to a multiple of this distance. This essentially
implies cutting off a layer from the cold end and adding a fresh layer at the
hot end with a value of $\phi=0$, and $c_{l}=c_{\infty}$.

Figure 4: (a) Comparison of the partition coefficient as a function of
velocity obtained from phase field simulations for different values of
parameter $A$. (b) Comparison of the interface temperature as function of
velocity computed for different values of $A$. The solid lines represent the
prediction from the CGM which are provided for the sake of comparison. $W$ is
fixed in the simulations and chosen to be $13.5$ nm and the value of $b$ is
selected as $1$ . In addition to trapping, $A$ controls the solute drag in the
system.

###  V.2 Effect of parameter $A$ on solute profiles

The effect of solute trapping parameter $A$ on the solute profiles
corresponding to $D_{s}/D_{l}=0.1$, $b=1.0$ and $v_{p}=0.32$
$\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ for the Si-9 % As system is
shown in Fig.3(a). The interface width has been fixed for all the simulations
and is chosen as $W=13.5$ nm. The maximum of the solute profiles decreases
with increase in the value of $A$ indicating a reduction of interfacial solute
segregation. The effect of $A$ on solute segregation also depends on the sign
of the expression $a(\phi)$ given by Eq. (54). If the condition $M_{s}<bM_{l}$
is satisfied (which is true for the above set of parameters), increasing the
value of $A$ increases the level of trapping. For the solute profiles shown in
Fig. 3(b), the parameters correspond to $D_{s}/D_{l}=0.001$ and $b=0.00001$
for which the condition $M_{s}>bM_{l}$ is satisfied. In such a case increasing
the value of $A$ decreases the level of solute segregation. Thus, the
parameter $A$ (or the source terms), in addition to the interfacial width $W$,
provides a degree of freedom to tailor solute partitioning.

Figure 5:  Comparison of the partition coefficient from the phase field
asymptotics obtained from the numerical solution of Eq. (106) with $A$, $b$
and $W$ as a fitting parameter and (a) the experimental results from Kittl et
al. for Si–$9\text{\,}\%$As [60] and (b) Molecular Dynamics data for
Ni–$5\text{\,}\%$Cu [61]. Figure 6:  Comparison of the partition coefficient
from the phase field asymptotics obtained from the numerical solution of Eq.
(106) and numerical solution of the governing equations for different
combinations of $A$, $b$ and $W$ for (a) Si–$9\text{\,}\%$As and (b)
Ni–$5\text{\,}\%$Cu. Same level of trapping can be maintained for larger
interface widths by appropriately tuning $A$ and $b$.

###  V.3 Effect of parameter $A$ on velocity-dependent partition coefficient
and interface temperature

The effect of model parameter $A$ on the velocity-dependent partition
coefficient $k_{v}$ is shown in Fig.4(a). The simulations have been performed
with $W=13.5$ nm, $D_{s}/D_{l}=0.1$ and $b=1.0$. The partition coefficient is
defined from the overall solute profile according to Ahmad et al. as [34]

$k_{v}=\frac{c_{s}}{c_{l}}=\frac{c_{\infty}}{c_{\mathrm{max}}}.$ (104)

where $c_{\mathrm{max}}$ denotes the maximum value in the solute concentration
profile. Different values of $A$ yield distinct $k_{v}$ curves. At low
velocities ($v\sim 0.02$ $\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$), $A$
has little effect on $k_{v}$ and all the curves tend to the equilibrium
partition coefficient ($k_{v}=k_{e}=0.3$). As the velocities increase, the
$k_{v}$ curve corresponding to higher values of $A$, in general, exhibits
higher levels of trapping. Thus, the parameter $A$ can be tuned to achieve a
desired level of velocity-dependent partitioning.

We next focus on the variation of interface temperature with the velocity. The
interface temperature $T_{I}$ corresponding to different values of $A$ are
shown in Fig.4(b). The interface temperature in phase field simulations is
defined as the temperature corresponding to $\phi=0.5$. The interface
temperature is velocity-dependent and different choice of $A$ leads to a
distinct interface temperature-velocity relationship. The interface
temperature corresponds to the solidus temperature at low velocities. With
increase in velocity, the interface temperature initially rises and then drops
at larger velocities when the effect of interface kinetics become significant.
The result is consistent with most sharp-interface models, for instance CGM
which predict a similar velocity-dependent interface temperature. The CGM
model predicts the interfacial temperature of the form [17],

$T_{I}=T_{\mathrm{m}}+\frac{mc_{\infty}}{k_{v}}\Big{[}\frac{1-k_{v}+\\{k_{v}+\mathcal{D}(1-k_{v})\\}\ln(\frac{k_{v}}{k_{e}})}{1-k_{e}}\Big{]}+\frac{V}{V_{o}}\frac{m}{1-k_{e}}.$
(105)

The quantity in the square bracket represents the velocity- dependent liquidus
slope and the third term denotes the correction due to interface kinetics.
Another prediction of the CGM is that the interface temperature depends on
solute drag through the drag parameter $\mathcal{D}$ in Eq. (105). The value
of $\mathcal{D}$ lies between $0$ and $1$ and the case $\mathcal{D}=0$ and
$\mathcal{D}=1$ corresponds to zero and full drag respectively. From Fig.4(b),
it is evident that the model parameter $A$ controls the solute drag as well.
This can also be realized from Eq. (64). Thus, a combination of value of $A$
and $W$ uniquely determines the velocity-dependent partition coefficient and
interface temperature. We contrast the prediction from the phase field model
with CGM which predicts a drag independent partition coefficient. The result
from the phase field model is consistent with the recent sharp-interface non-
equilibrium model by Hareland et al. [62] which shows that the solute drag
affects the interfacial composition of the phases. A comparison of the
interface temperature between the phase field model and the CGM indicates that
as the value of $A$ increases, the amount of solute drag decreases in the
system.

Figure 7:  (a) Comparison of the present with previously developed phase field
models and the allowable interface widths that emulates the experimental data
from Kittl et al. [60] in Si-9 % As and (b) Comparison of the interface
temperature from different phase field models.

###  V.4 Tailoring the partition coefficient: Selection of $A$

By construction, the model parameter $A$ controls the level of solute
segregation. The next question that naturally arises is how to choose a value
of $A$ that captures the “correct” amount of trapping. While no closed form
expression for the partition coefficient can be deduced for the model
$\mu_{s}\neq\mu_{l}$, we again proceed by employing the expressions derived
for the model with the assumption $\mu_{s}=\mu_{l}$. The diffusion potential
jump expression Eq. (37) can be re-written for a dilute alloy in dimensional
form as,

$\ln\left(\frac{k_{v}}{k_{e}}\right)=\frac{VW}{D_{l}}(1-k_{v})\Delta F$ (106)

where $\Delta F=F^{-}-F^{+}$, giving a transcendental equation which can be
solved numerically to give a velocity-dependent $k_{v}$ curve for a given
combination of $A$, $b$ and $W$. This implies that one can increase $W$ for
computational convenience and still obtain an “experimentally relevant
trapping” by appropriately choosing the constants $A$ and $b$. This
observation is also true for the general model where $\mu_{s}\neq\mu_{l}$.
There are several ways in which an optimal combination of the model parameters
can be chosen. The first obvious choice is to fit the $k_{v}$ from Eq. (106)
to an experimental dataset for a range of velocities. For the sake of
comparison we choose the alloy system Si–$9\text{\,}\%$As, for which the
experimental data is readily available in the literature [60]. With $A$, $b$
and $W$ as a fitting parameter, the numerical solution of Eq. (106) is fit to
the experimental data within the velocity range of
$0.1\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ to
$1\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$. A value of $A=2.0$,
$b=0.001$ and $W=13.5$ nm with $D_{s}/D_{l}=10^{-3}$ fits the data with
reasonable agreement as presented in Fig.5(a). Other combinations of the model
parameters can also be selected to obtain the same level of trapping. We
emphasize that the above procedure gives an approximate or a first estimate of
the model parameters, since Eq. (106) is derived for the model
$\mu_{s}=\mu_{l}$. In our experience, we have observed that the $k_{v}$ curve
for the model $\mu_{s}=\mu_{l}$ closely resembles the prediction of the
general model upto velocities of about $0.3$
$\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ in most cases for the
materials parameters given in Table 1.

In absence of experimental data for an alloy, one can rely on atomistic
simulation such as molecular dynamics (MD) or sharp-interface theories to
choose an appropriate value of the model parameters. Fig.5(b) represents the
fit from solution of Eq. (106) to the MD simulation data on
Ni–$5\text{\,}\%$Cu from Yang et al. [61] in the velocity range of
$v=$0\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$$ to
$2\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$. In Fig.5(b) $A=1.8$,
$b=0.00001$ and $W=7.8$ nm has been chosen to best fit the data from the MD
simulation.

Once an estimate of the model parameters are derived, numerical simulations
are performed to check the validity of the asymptotic expressions. A
comparison of the prediction of Eq. (106) and numerical results from the phase
field model for different combinations of $A$, $b$ and $W$ are shown in
Fig.6(a) for the system Si–$9\text{\,}\%$As and for the Ni–$5\text{\,}\%$Cu
system in Fig.6(b). The results for the phase field model are within $15\%$ of
the prediction of the analytical expression. At higher velocities, the choice
of $W$ is restricted by the requirement to resolve the diffusion length which
is given by the ratio of interface diffusivity and pulling speed. It is
evident from Fig.6 that the interface width $W$ can be increased without
changing the partition coefficient by simultaneously tuning the values of
parameters $A$ and $b$. It is to be noted that the choice of $A$, $b$ and $W$
should also be chosen to capture the appropriate (measured experimentally or
through atomistic simulations) velocity-dependent interface temperature.

Figure 8:  Solute profiles exhibiting solute bands corresponding to (a)
$t=1.92\times 10^{-4}$ sec, (b) $t=2.56\times 10^{-4}$ sec, (c) $t=3.84\times
10^{-4}$ sec and (d) $t=5.12\times 10^{-4}$ sec. Time periodic oscillations of
the (e) interface velocity and (f) interface temperature. (h) Variation of
oscillation wavelength with the temperature gradient.

###  V.5 Comparison with other phase field models

We next compare the present model to a number of previously developed phase
field models. The comparison is restricted to models which arise as a special
case of the present model. A direct numerical comparison with these models
allows us to quantify the progress made due to the addition of interfacial
source terms. The intention is to compare the order of interface widths $W$
that can be employed in different phase field models. We first examine the
Si–As system in Fig.7(a). The curve denoted by Ahmad et al. is obtained from
the high velocity asymptotics of the WBM model [34]. $k_{v}$ has a similar
form to the Aziz trapping function [16] given by

$k_{v}=\frac{k_{v}+V/V_{D}}{1+V/V_{D}},$ (107)

where $V_{D}$ is given by

$V_{D}=\frac{3}{8}\frac{D_{I}}{W}\frac{\ln(1/k_{e})}{1-k_{e}}$ (108)

Since $W$ is the only free parameter in the model, to emulate
$V_{D}=$0.68\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$$, one can use a
value of $W$ no larger than $3.5$ nm which was used in the numerical study of
Danilov and Nestler [37] employing the WBM model. For the KKS model, which is
a special case of the present model (with $\mu_{s}=\mu_{l}$ and $A=1$), the
asymptotic analysis reveals that the maximum $W$ that can be chosen is about
$13.5$ nm. Since the present model has an additional degrees of freedom (in
$A$ and $b$), one can employ roughly four and half times the interface
thickness ($60.75$ nm) and flush out the extra solute from the interface by
appropriately choosing $A$ and $b$. While the asymptotic analysis of the
Pinomaa model [44] predicts the choice of $W=20$ nm, numerical studies
revealed a good convergence only until velocity of about
$0.45\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ [44]. For the present
model convergence up to $0.75\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$
can be seen in Fig.6.

The choice of $W$ which is treated as a fitting parameter in phase field
models strictly depends on the solute trapping behaviour of the material
system of interest. The solute trapping tendency of the system can be
characterized by $V_{D}$. A lower value of $V_{D}$ implies a higher trapping
tendency of the system or a steeper increase of the partition coefficient with
velocity and vice versa. Indeed, if a similar analysis is repeated for the
Ni–Cu system which has
$V_{D}=$1.5\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$$, Eq. (108) predicts
$W=0.5$ nm which brings it closer to the values that were employed in
numerical studies in Ref.[34]. An interface width of about $1$ nm can be
employed for the KKS model as evident from the asymptotic analysis (with $A=1$
and assuming $\mu_{s}=\mu_{l}$) and numerical studies for the same system in
Ref.[36]. The present model allows one to choose $W$ in the range of $7.85$ nm
to $13$ nm as evident from the convergence studies presented in Fig.6(b).
Since the smallest grid spacing needed for numeric computation scales with
interface thickness, increasing the interface thickness permits use of coarser
grid spacing and also the allowable timestep width. For an explicit finite
difference algorithm, with uniform grid spacing, the number of floating point
operations increases proportionally with the number of grid points which, in
turn, scales as $1/W^{d}$, where, $d$ denotes the number of dimensions. The
number of floating point operations also depends on the number of simulation
timesteps which scales inversely with the timestep width $\Delta t$ which in
turn scales as $W^{2}$. Therefore, the number of floating point operations
scales with the interface width as $1/W^{d+2}$. Therefore, increasing the
interface width by a factor of two speeds up the computation time by a factor
of $16$ in two dimensions and $32$ in three dimensions. The above analysis
implies a tremendous gain in computation time when the model is applied to
simulate microstructures in two and three dimensions.

A comparison of the interface temperature for the Si–As system is presented in
Fig.7(b). The interface thickness for the different models are listed in
Fig.7(a). Not all phase field models have the same degree of freedom in tuning
the interface temperature. The interface temperature can be controlled in the
present model by the parameters $A$, $b$ and $W$. Different choice of the
model parameters gives rise to distinct curves as shown in Fig.7(b). The high
velocity asymptotics of the WBM model predicts the interface temperature of
the form given by Eq. (105) with partial drag ($\mathcal{D}=24/35$) embedded
in the model. The interface temperature was only presented until velocity of
$0.2\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ in the study of Pinomaa et
al. [44] and hence not included here. Different interface temperatures for
different models indicate that the solute drag is sensitive to the specifics
of the model construction. We remark that the interface temperature can be
made to fall on the high-velocity asymptotic solution by appropriately
choosing the model parameters $A=1.4$, $b=0.00001$ and $W$ to $60.75$ nm
without appreciably changing the $k_{v}$ curve as shown earlier in Fig.(6).
However, interface temperature alone cannot be changed without simultaneously
changing $k_{v}$ if $W$ is altered in the WBM and KKS model. Thus, while the
model by Ahmad et al. is restricted by the choice of $W=3.5$ nm to reproduce
the experimental $k_{v}$ curve, the present model allows to increase the
interface width by $60.75$ nm to emulate the same $k_{v}$ and interface
temperature curve.

Figure 9:  Time periodic oscillations of the (a) interface velocity and (b)
interface temperature. (c) Variation of oscillation wavelength with the
temperature gradient.

###  V.6 Oscillations and solute bands

Hitherto, we have focused on the steady state interface motion. In rapid
solidification experiments on metallic alloys such as Ag-Cu, Al-Cu and Al-Fe
systems [63, 64], found microstructures that are periodic along the growth
direction. Such structures have been termed as ”banded structure” . The bands
consists of alternate dark bands of either cellular-dendritic or eutectic
structure, depending on the alloy composition and light bands free from any
microsegregation. The origin of banding has been rationalized in terms of an
oscillatory instability of the planar solid-liquid interface with zero
wavenumber which is characterized by large variations of the interface
velocity. Incorporating non-equilibrium effects into the Mullins-Sekerka
theory [65] i.e. velocity dependent partition coefficient and interface
temperature, Coriell and Sekerka [66], Merchant and Davis [14] and Karma and
Sarkissian [67] have independently identified an oscillatory instability with
a zero wavenumber along the solid-liquid front.

Our phase field model, which naturally incorporates the non-equilibrium
effects, is a perfect tool to investigate such an instability. We present
solute profiles exhibiting solute bands in Fig.8 corresponding to $G=1.4\times
10^{6}$ $\text{\,}\mathrm{K}\text{\,}{\mathrm{m}}^{-1}$ and $v_{p}=$
$0.3\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$. The velocity is selected
such that it corresponds to $\frac{\mathrm{d}T}{\mathrm{d}V}>0$ regime of the
interface temperature-velocity curve, where such instabilities are expected.
As the solid-liquid front advances, it leaves behind a trail of solute rich
and solute poor zones. The region between the bands is free of solute
segregation. The amplitude of the peaks and troughs of the bands are almost
constant owing to the negligible solute diffusion in solid
($D_{s}/D_{l}=0.001$). Furthermore, the solid-liquid front undergoes time-
periodic oscillation as characterized by the interface velocity and
temperature shown in Fig.9(a) and (b). The liquid composition of the interface
characterized by the maximum of the front composition also undergoes time-
periodic variation as evident from the solute profiles(Fig.8(a)-(d)).The bands
are evenly spaced with a wavelength of about
$4\text{\,}\mathrm{\SIUnitSymbolMicro m}$. While we do not aim to make a
quantitative comparison with the aforementioned analytical theories here, the
qualitative dependence of the oscillation frequency on $G$ can be observed in
Fig.9(c). As predicted by Merchant and Davis, the oscillation wavelength
decreases with increase in the value of $G$.

##  VI Conclusions

We have presented a variational phase field model in which non-equilibrium
effects such as solute trapping, drag and interface kinetics can be introduced
in a controlled manner. The key feature of the model is the derivation of
separate kinetic equations for the phase concentrations with interfacial
source terms which can be selected to control solute trapping and drag.
Separate governing equations for the phase concentrations allows for a
straightforward incorporation of arbitrary alloy model without having to solve
the condition of pointwise equality of diffusion potentials which is often the
computation bottle-neck in the previous works. To perform simulations with the
present model following steps needs to be performed:

  * •

For a given alloy system, the free energies of the solid and the liquid phases
have to be specified first. While the free energies in the present work have
been selected as an ideal dilute alloy with straight solidus and liquidus, the
model is capable to take input from thermodynamic databases such as CALPHAD to
accommodate more realistic free energies. Since $\mu_{s}$ is not constrained
to be equal to $\mu_{l}$, there is no need to solve a coupled set of nonlinear
equations to determine effective phase concentrations as in the KKS model.

  * •

Next, we select the mobility interpolation functions in the normal and
tangential directions to counter Kapitza jump and surface diffusion
respectively based on the discussion in section II.3. Depending on the ratio
of the chosen diffusivity of the phases $D_{s}/D_{l}$, the parameters $b$ and
$a$ in $q_{n}(\phi)$ (Eq. (51)) are evaluated such that $G^{+}=G^{-}$ in Eq.
(36) is satisfied. $b$ is usually chosen such that the diffusivity at the
interfacial region is high which is true if $b<D_{s}/D_{l}$. The choice of $b$
and $a$ should also satisfy $q_{n}(\phi)\geq 0$ for the phase concentration
equations to be stable.

  * •

To capture a specific partition coefficient and interface temperature as a
function of velocity, the source term $A$, the diffusivity shape parameter $b$
and interface width $W$ can be tuned by comparing the numerical solution of
asymptotics Eq. (106) for partition coefficient to an experimental dataset,
atomistic simulations or sharp-interface theories. While this method provides
various combinations of the model parameters, the choice of $A$, $b$ and $W$
should additionally capture the desired interface temperature (or solute
drag).

Numerical performance of the model was demonstrated through study of the
partition coefficient and interface temperature at different interface widths.
A comparison with the previous phase field models reveal that “experimentally-
relevant” trapping can be achieved with interface widths of about $25$ times
larger than the WBM model [34], about four times that can be employed in the
KKS [36] and about $3$ times the Pinomaa [44] model depending upon the
material system at hand. Our phase field model was able to capture the
oscillatory instability during rapid solidification leading to solute bands
which are in qualitative agreement with the analytical and experimental
predictions. Future two-dimensional studies is currently underway to
quantitatively identify the regimes of cellular, dendrite and bands depending
on the rapid solidification parameters (such as temperature gradient and
pulling speeds).

###### Acknowledgements.

This work was performed under financial assistance award 70NANB14H012 from the
U.S. Department of Commerce, National Institute of Standards and Technology as
part of the Center for Hierarchical Materials Design (CHiMaD).

## References

  * Cui  _et al._ [2021] X. Cui, S. Zhang, C. Zhang, J. Chen, J. Zhang, and S. Dong, Additive manufacturing of 24CrNiMo low alloy steel by selective laser melting: Influence of volumetric energy density on densification, microstructure and hardness, Materials Science and Engineering: A 809, 140957 (2021).
  * Wolf  _et al._ [2020] W. Wolf, L. P. M. E. Silva, G. Zepon, C. S. Kiminami, C. Bolfarini, and W. J. Botta, Single step fabrication by spray forming of large volume Al-based composites reinforced with quasicrystals, Scripta Materialia 181, 86 (2020).
  * Oliveira  _et al._ [2020] J. Oliveira, T. Santos, and R. Miranda, Revisiting fundamental welding concepts to improve additive manufacturing: from theory to practice, Progress in Materials Science 107, 100590 (2020).
  * Kobold  _et al._ [2017] R. Kobold, W. Kuang, H. Wang, W. Hornfeck, M. Kolbe, and D. M. Herlach, Dendrite growth velocity in the undercooled melt of glass forming Ni50Zr50 compound, Philosophical Magazine Letters 97, 249 (2017).
  * Ghosh  _et al._ [2018] S. Ghosh, L. Ma, L. E. Levine, R. E. Ricker, M. R. Stoudt, J. C. Heigel, and J. E. Guyer, Single-track melt-pool measurements and microstructures in Inconel 625, JOM 70, 1011 (2018).
  * Chen  _et al._ [2021] M. Chen, Q. Du, H. Fu, X. He, Z. Huang, R. Shi, and J. Xie, The formation of secondary phase in sub-rapid solidification process of Al–Mg–Si alloys, Materialia 15, 101022 (2021).
  * Enrique  _et al._ [2021] P. D. Enrique, A. Keshavarzkermani, and N. Y. Zhou, Effect of microsegregation on high-temperature microstructure evolution in rapid solidification processed Nb-rich Ni superalloys, Advanced Engineering Materials 23, 2001396 (2021).
  * Keller  _et al._ [2017] T. Keller, G. Lindwall, S. Ghosh, L. Ma, B. M. Lane, F. Zhang, U. R. Kattner, E. A. Lass, J. C. Heigel, Y. Idell, _et al._ , Application of finite element, phase-field, and CALPHAD-based methods to additive manufacturing of Ni-based superalloys, Acta Materialia 139, 244 (2017).
  * Ghosh  _et al._ [2017] S. Ghosh, L. Ma, N. Ofori-Opoku, and J. E. Guyer, On the primary spacing and microsegregation of cellular dendrites in laser deposited Ni–Nb alloys, Modelling and Simulation in Materials Science and Engineering 25, 065002 (2017).
  * Karayagiz  _et al._ [2020] K. Karayagiz, L. Johnson, R. Seede, V. Attari, B. Zhang, X. Huang, S. Ghosh, T. Duong, I. Karaman, A. Elwany, _et al._ , Finite interface dissipation phase field modeling of Ni–Nb under additive manufacturing conditions, Acta Materialia 185, 320 (2020).
  * Chadwick and Voorhees [2021] A. F. Chadwick and P. W. Voorhees, The development of grain structure during additive manufacturing, Acta Materialia 211, 116862 (2021).
  * Wang  _et al._ [2019] X. Wang, P. Liu, Y. Ji, Y. Liu, M. Horstemeyer, and L. Chen, Investigation on microsegregation of IN718 alloy during additive manufacturing via integrated phase-field and finite-element modeling, Journal of Materials Engineering and Performance 28, 657 (2019).
  * Kurz and Trivedi [1994] W. Kurz and R. Trivedi, Rapid solidification processing and microstructure formation, Materials Science and Engineering: A 179, 46 (1994).
  * Merchant and Davis [1990] G. Merchant and S. Davis, Morphological instability in rapid directional solidification, Acta Metallurgica et Materialia 38, 2683 (1990).
  * Aziz [1982] M. J. Aziz, Model for solute redistribution during rapid solidification, Journal of Applied Physics 53, 1158 (1982).
  * Aziz and Kaplan [1988] M. J. Aziz and T. Kaplan, Continuous growth model for interface motion during alloy solidification, Acta Metallurgica 36, 2335 (1988).
  * Aziz and Boettinger [1994] M. J. Aziz and W. Boettinger, On the transition from short-range diffusion-limited to collision-limited growth in alloy solidification, Acta Metallurgica et Materialia 42, 527 (1994).
  * Ågren [1989] J. Ågren, A simplified treatment of the transition from diffusion controlled to diffusion-less growth, Acta Metallurgica 37, 181 (1989).
  * Jackson  _et al._ [2004] K. A. Jackson, K. M. Beatty, and K. A. Gudgel, An analytical model for non-equilibrium segregation during crystallization, Journal of Crystal Growth 271, 481 (2004).
  * Baker and Cahn [1971] J. Baker and J. Cahn, Solidification, American Society for Metals, Metals Park, OH 23 (1971).
  * Hillert and Sundman [1977] M. Hillert and B. Sundman, A solute-drag treatment of the transition from diffusion-controlled to diffusionless solidification, Acta Metallurgica 25, 11 (1977).
  * Baker [1970] J. C. Baker,  _Interfacial partitioning during solidification_ , Ph.D. thesis, Massachusetts Institute of Technology (1970).
  * Boettinger  _et al._ [2002] W. J. Boettinger, J. A. Warren, C. Beckermann, and A. Karma, Phase-field simulation of solidification, Annual Review of Materials Research 32, 163 (2002).
  * Tourret  _et al._ [2022] D. Tourret, H. Liu, and J. Llorca, Phase-field modeling of microstructure evolution: Recent applications, perspectives and challenges, Progress in Materials Science 123, 100810 (2022).
  * Chen [2002] L.-Q. Chen, Phase-field models for microstructure evolution, Annual Review of Materials Research 32, 113 (2002).
  * Moelans  _et al._ [2008] N. Moelans, B. Blanpain, and P. Wollants, An introduction to phase-field modeling of microstructure evolution, Calphad 32, 268 (2008).
  * Steinbach [2013] I. Steinbach, Phase-field model for microstructure evolution at the mesoscopic scale, Annual Review of Materials Research 43, 89 (2013).
  * Karma [2001] A. Karma, Phase-field formulation for quantitative modeling of alloy solidification, Physical Review Letters 87, 115701 (2001).
  * Almgren [1999] R. F. Almgren, Second-order phase field asymptotics for unequal conductivities, SIAM Journal on Applied Mathematics 59, 2086 (1999).
  * Karma and Rappel [1998] A. Karma and W.-J. Rappel, Quantitative phase-field modeling of dendritic growth in two and three dimensions, Physical Review E 57, 4323 (1998).
  * Ohno and Matsuura [2009] M. Ohno and K. Matsuura, Quantitative phase-field modeling for dilute alloy solidification involving diffusion in the solid, Physical Review E 79, 031603 (2009).
  * Zhang  _et al._ [2013] L. Zhang, E. V. Danilova, I. Steinbach, D. Medvedev, and P. K. Galenko, Diffuse-interface modeling of solute trapping in rapid solidification: Predictions of the hyperbolic phase-field model and parabolic model with finite interface dissipation, Acta Materialia 61, 4155 (2013).
  * Mullis  _et al._ [2010] A. Mullis, J. Rosam, and P. Jimack, Solute trapping and the effects of anti-trapping currents on phase-field models of coupled thermo-solutal solidification, Journal of Crystal Growth 312, 1891 (2010).
  * Ahmad  _et al._ [1998] N. Ahmad, A. Wheeler, W. J. Boettinger, and G. B. McFadden, Solute trapping and solute drag in a phase-field model of rapid solidification, Physical Review E 58, 3436 (1998).
  * Wheeler  _et al._ [1992] A. A. Wheeler, W. J. Boettinger, and G. B. McFadden, Phase-field model for isothermal phase transitions in binary alloys, Physical Review A 45, 7424 (1992).
  * Kim  _et al._ [1999] S. G. Kim, W. T. Kim, and T. Suzuki, Phase-field model for binary alloys, Physical Review E 60, 7186 (1999).
  * Danilov and Nestler [2006] D. Danilov and B. Nestler, Phase-field modelling of solute trapping during rapid solidification of a Si–As alloy, Acta Materialia 54, 4659 (2006).
  * Conti [1997] M. Conti, Solute trapping in directional solidification at high speed: A one-dimensional study with the phase-field model, Physical Review E 56, 3717 (1997).
  * Kim and Kim [2001] S. G. Kim and W. T. Kim, Phase-field modeling of rapid solidification, Materials Science and Engineering: A 304, 281 (2001).
  * Glasner [2001] K. Glasner, Solute trapping and the non-equilibrium phase diagram for solidification of binary alloys, Physica D: Nonlinear Phenomena 151, 253 (2001).
  * Wheeler  _et al._ [1993] A. A. Wheeler, W. J. Boettinger, and G. B. McFadden, Phase-field model of solute trapping during solidification, Physical Review E 47, 1893 (1993).
  * Boettinger and Warren [1999] W. J. Boettinger and J. A. Warren, Simulation of the cell to plane front transition during directional solidification at high velocity, Journal of Crystal Growth 200, 583 (1999).
  * Steinbach  _et al._ [2012] I. Steinbach, L. Zhang, and M. Plapp, Phase-field model with finite interface dissipation, Acta Materialia 60, 2689 (2012).
  * Pinomaa and Provatas [2019] T. Pinomaa and N. Provatas, Quantitative phase field modeling of solute trapping and continuous growth kinetics in quasi-rapid solidification, Acta Materialia 168, 167 (2019).
  * Echebarria  _et al._ [2004] B. Echebarria, R. Folch, A. Karma, and M. Plapp, Quantitative phase-field model of alloy solidification, Physical Review E 70, 061604 (2004).
  * Ji  _et al._ [2022] K. Ji, E. Dorari, A. J. Clarke, and A. Karma, Microstructural pattern formation during far-from-equilibrium alloy solidification, arXiv preprint arXiv:2209.11352 (2022).
  * Brener and Boussinot [2012] E. A. Brener and G. Boussinot, Kinetic cross coupling between nonconserved and conserved fields in phase field models, Physical Review E 86, 060601 (2012).
  * Boussinot and Brener [2014] G. Boussinot and E. A. Brener, Achieving realistic interface kinetics in phase-field models with a diffusional contrast, Physical Review E 89, 060402 (2014).
  * Wang  _et al._ [2013] H. Wang, F. Liu, G. Ehlen, and D. Herlach, Application of the maximal entropy production principle to rapid solidification: a multi-phase-field model, Acta Materialia 61, 2617 (2013).
  * Fang and Mi [2013] A. Fang and Y. Mi, Recovering thermodynamic consistency of the antitrapping model: A variational phase-field formulation for alloy solidification, Physical Review E 87, 012402 (2013).
  * Ohno  _et al._ [2016] M. Ohno, T. Takaki, and Y. Shibuta, Variational formulation and numerical accuracy of a quantitative phase-field model for binary alloy solidification with two-sided diffusion, Physical Review E 93, 012802 (2016).
  * Ohno  _et al._ [2017] M. Ohno, T. Takaki, and Y. Shibuta, Variational formulation of a quantitative phase-field model for nonisothermal solidification in a multicomponent alloy, Physical Review E 96, 033311 (2017).
  * Wang  _et al._ [2020] K. Wang, G. Boussinot, C. Hüter, E. A. Brener, and R. Spatschek, Modeling of dendritic growth using a quantitative nondiagonal phase field model, Physical Review Materials 4, 033802 (2020).
  * Wang  _et al._ [2021] K. Wang, G. Boussinot, E. A. Brener, and R. Spatschek, Quantitative nondiagonal phase field modeling of eutectic and eutectoid transformations, Physical Review B 103, 184111 (2021).
  * Nicoli  _et al._ [2011] M. Nicoli, M. Plapp, and H. Henry, Tensorial mobilities for accurate solution of transport problems in models with diffuse interfaces, Physical Review E 84, 046707 (2011).
  * Gugenberger  _et al._ [2008] C. Gugenberger, R. Spatschek, and K. Kassner, Comparison of phase-field models for surface diffusion, Physical Review E 78, 016703 (2008).
  * Ahmed  _et al._ [2013] K. Ahmed, C. Yablinsky, A. Schulte, T. Allen, and A. El-Azab, Phase field modeling of the effect of porosity on grain growth kinetics in polycrystalline ceramics, Modelling and Simulation in Materials Science and Engineering 21, 065005 (2013).
  * Provatas and Elder [2011] N. Provatas and K. Elder,  _Phase-field methods in materials science and engineering_ (John Wiley & Sons, 2011).
  * Algoso  _et al._ [2003] P. Algoso, W. Hofmeister, and R. Bayuzick, Solidification velocity of undercooled Ni–Cu alloys, Acta Materialia 51, 4307 (2003).
  * Kittl  _et al._ [2000] J. Kittl, P. Sanders, M. Aziz, D. Brunco, and M. Thompson, Complete experimental test of kinetic models for rapid alloy solidification, Acta Materialia 48, 4797 (2000).
  * Yang  _et al._ [2011] Y. Yang, H. Humadi, D. Buta, B. B. Laird, D. Sun, J. J. Hoyt, and M. Asta, Atomistic simulations of nonequilibrium crystal-growth kinetics from alloy melts, Physical Review Letters 107, 025505 (2011).
  * Hareland  _et al._ [2022] C. A. Hareland, G. Guillemot, C.-A. Gandin, and P. W. Voorhees, The thermodynamics of non-equilibrium interfaces during phase transformations in concentrated multicomponent alloys, Acta Materialia , 118407 (2022).
  * Boettinger  _et al._ [1984] W. Boettinger, D. Shechtman, R. Schaefer, and F. Biancaniello, The effect of rapid solidification velocity on the microstructure of Ag-Cu alloys, Metallurgical Transactions A 15, 55 (1984).
  * Kurz and Trivedi [1996] W. Kurz and R. Trivedi, Banded solidification microstructures, Metallurgical and Materials Transactions A 27, 625 (1996).
  * Mullins and Sekerka [1964] W. W. Mullins and R. Sekerka, Stability of a planar interface during solidification of a dilute binary alloy, Journal of Applied Physics 35, 444 (1964).
  * Coriell and Sekerka [1983] S. Coriell and R. Sekerka, Oscillatory morphological instabilities due to non-equilibrium segregation, Journal of Crystal Growth 61, 499 (1983).
  * Karma and Sarkissian [1993] A. Karma and A. Sarkissian, Interface dynamics and banding in rapid solidification, Physical Review E 47, 513 (1993).