# Variational quantitative phase-field modeling of non-isothermal sintering
process

Timileyin David Oyedeji<EMAIL_ADDRESS>Yangyiwei Yang
<EMAIL_ADDRESS>Mechanics of Functional Materials Division,
Institute of Materials Science, Technische Universität Darmstadt, 64287
Darmstadt, Germany  Herbert Egger  Johann Radon Institute for Computational
and Applied Mathematics and Institute for Computational Mathematics, Johannes-
Kepler University Linz, 4040 Linz, Austria  Bai-Xiang Xu<EMAIL_ADDRESS>darmstadt.de Mechanics of Functional Materials Division, Institute of
Materials Science, Technische Universität Darmstadt, 64287 Darmstadt, Germany

###### Abstract

Phase-field modeling has become a powerful tool in describing the complex
pore-structure evolution and the intricate multi-physics in non-isothermal
sintering processes. However, the quantitative validity of conventional
variational phase-field models involving diffusive processes is a challenge.
Artificial interface effects, like the trapping effects, may originate at the
interface when the kinetic properties of two opposing phases are different. On
the other hand, models with prescribed antitrapping terms do not necessarily
guarantee the thermodynamics variational nature of the model. This issue has
been solved for liquid-solid interfaces via the development of the variational
quantitative solidification phase-field model. However, there is no related
work addressing the interfaces in non-isothermal sintering, where the free
surfaces between the solid phase and surrounding pore regions exhibit strong
asymmetry of mass and thermal properties. Also, additional challenges arise
due to the conserved order parameter describing the free surfaces. In this
work, we present a variational and quantitative phase-field model for non-
isothermal sintering processes. The model is derived via an extended non-
diagonal phase-field model. The model evolution equations have naturally
cross-coupling terms between the conserved kinetics (i.e., mass and thermal
transfer) and the non-conserved one (grain growth). These terms are shown via
asymptotic analysis to be instrumental in ensuring the elimination of
interface artefacts, while also examined to not modify the thermodynamic
equilibrium condition (characterized by dihedral angle). Moreover, we
demonstrate that the trapping effects and existence of surface diffusion in
conservation laws are direction-dependent. An anisotropic interpolation scheme
of the kinetic mobilities which differentiates the normal and the tangential
directions along the interface is discussed. Numerically, we demonstrate the
importance of the cross-couplings and the anisotropic interpolation via
presenting thermal-microstructural evolutions.

††preprint: APS/123-QED

##  I Introduction

Sintering is a typical densification technique in thermal processing of bulk
materials from packed powders [1, 2, 3]. In the present day, many new
techniques based on sintering have been proposed and broadly applied in the
industry, where the thermal bonding effect is introduced by treatments other
than direct heating, such as laser scan, electrical current and
electromagnetic field [4, 3, 5, 6]. These techniques are collectively termed
as “unconventional” sintering [7, 8]. Due to the distinct heating mechanisms
among unconventional sintering techniques, effects of non-isothermal factors
on the properties of products, like heating/cooling rate and temperature
inhomogeneity gain increasing attentions alongside the conventional ones such
as chemical composition as well as size of powders, atmosphere, and pressure.

Therefore, it is essential to identify and understand the physical effects and
interactions of these factors in the sense of bridging the process parameters,
microstructure, and properties of the materials to further tailor the
performance for applications of interest. Two major types of interfaces are
essential for sintering process, namely the free surface between pore and
substance and the grain boundary between adjacent crystal grains. There are
analytical models for describing the pores/grains evolution, the two-particle
coalescence model [9, 10], dodeca-/tetrakaidecahedron grain model [11], and
various models treating the pores/grains through assumed geometries, like
spheres or cylinders [12, 13]. Nevertheless, complex grain/pore geometry and
entangled multiple physics during sintering goes beyond the capacity of these
models.

For such purpose, phase-field modeling and simulation is promising. In the
conventional variational phase-field theory, order parameters (OP) are applied
to represent the spatio-temporal distribution of microstructure, i.e. pores
and grain orientations in the case of sintering. The thermodynamic potential
of the microstructure can then be formulated by an energy functional w.r.t.
the OPs, including the interface contribution through the corresponding
gradient terms of OPs. From non-equilibrium thermodynamics, the evolution
equations of the OPs can be derived on the basis of the variational theory. It
circumvents the necessity of interface tracking. There are variational phase-
field sintering models considering an isothermal scenario. For instance,
Kazaryan et al. [14] and Wang [15] proposed a line of phase-field model, which
was used later for studying two-particle necking and coalescence [16, 17, 18,
19, 20, 21, 22] and densification of porous microstructure [17, 23], and in
simulating the overall microstructure evolution of the particle aggregation
[24] or particle stack [15, 25]. Rigid-body motions were also incorporated
within the model [15, 25, 20, 21]. Furthermore, a phase-field sintering model
adopting the grand potential concept was also developed [26, 27]. To simulate
sintering process under highly heterogeneous thermal environment, the phase-
field sintering model coupled with transient heat and/or chemical diffusion
simulations are needed. The phase-field approach allows such consideration
through additive inclusion of the energy contributions by the related physical
fields, such as temperature or chemical concentration. In our previous work
[28], a variational non-isothermal phase-field sintering model was proposed,
which was applied for simulations of the selective laser sintering on a
single-layer and multi-layer [29] powder beds, and for sintering under
prescribed high temperature gradient [8].

On the other hand, one theoretical issue of the conventional variational
phase-field models involving thermal/chemical diffusive process is the
quantitative validity. Artificial interface effects may originate from
violation of conservation laws and discontinuity of the chemical/thermal
potentials at the interface (trapping effects) [30, 31]. These interface
effects scale with the interface width. Theoretically, via asymptotic
analysis, phase-field models should be reduced to their associated free-
boundary problems in order to ensure their quantitative validity. Based on
thin-interface limit analysis, Karma and Rappel [32, 33] first published a
quantitative phase-field model for the solidification of pure materials with
equal diffusivities in the solid and liquid phases. Moreover, by introducing
an antitrapping term in the diffusion flux equation in order to eliminate the
trapping effect, Karma [34] presented a quantitative model for the case of
isothermal solidification of alloys with negligible diffusivity in the solid
phase. Furthermore, for the case with arbitrarily different diffusivities in
opposing phases, corresponding antitrapping terms have been also proposed for
isothermal [35] and non-isothermal consideration [36]. Thereby a new parameter
relating the interface velocity and diffusion flux was further introduced to
ensure full elimination of all interface artifacts.

It should be noted that modifying a variationally derived evolution equation
by prescribed antitrapping terms do not necessarily guarantee the variational
nature of the model, which is, however, important for thermodynamics
soundness. Therefore there have been efforts to develop variational
formulations of quantitative phase-field models. Using phenomenological linear
relations, variational formulation of quantitative phase-field models have
been developed by considering kinetic cross-coupling between the conserved
diffusion fields and the nonconserved OPs (non-diagonal model) [37, 38, 39,
40]. Time evolution equations of the models then exhibit cross-coupling
kinetic terms that are formulated in a similar fashion due to Onsager’s
symmetry. Furthermore, the parameters of these coupling terms are explicitly
formulated in terms of the models parameters by considering relations between
the models and their sharp-interface counterparts. The cross-coupling term in
the diffusion equations which can be likened to the antitrapping term
alongside the coupling term in the phase-field evolution equations have been
noted to enable full elimination of artificial interface effects [40, 41]. The
non-diagonal model has been employed to investigate quantitative phase-field
simulations of dendritic growth [42] and to examine quantitative simulations
of eutectic and eutectoid transformations [43] in which the necessity of the
cross-coupling terms were substantiated in both instances.

By separately considering the thermodynamic quantities of two opposing phases
and then treating the interface as a mixture of the phases (two-phase
variational approach), Ohno et al. [44, 45] presented quantitative variational
phase-field models for binary alloy solidification with two-sided diffusion.
In the two-phase formulation the diffusion fields mixture laws are ensured at
the interface as constraints implemented by the Lagrange multiplier approach,
and the flux fields of each single-phase fields are formulated variationally.
Emergence of Lagrange multipliers in thermodynamic potential formulation gives
rise to cross-coupling terms in the model time evolution equations which serve
to eliminate the artificial interface effects. Additionally, the necessity of
an anisotropic interpolation of the diffusivity (different interpolations for
the normal and the tangential directions across the diffuse interface) is
demonstrated in eliminating the anomalous interface effects. Though the two-
phase variational approach is promising for the study of quantitative
validity, the variational nature of the model is only implicitly implemented
through variationally formulated single-phase fluxes. The variational behavior
of the final model after inserting the Lagrange multiplier still needs to be
examined. Moreover, due to the assumptions of negligible temperature jump or
chemical potential jump across the diffuse interface, the models in Refs. [44,
45] are applicable mostly for slow solidification processes.

Based on literature review, there is currently no variational quantitative
phase-field model for non-isothermal sintering. In comparison to the non-
isothermal solidification models with non-conserved OPs, additional challenge
can be expected due to the conserved OPs involved here. In this work, we
derive variational formulation of a quantitative phase-field model for non-
isothermal sintering processes where the free surfaces between the solid phase
and surrounding atmosphere/pore regions have strong asymmetry of both mass and
thermal properties. The model is derived via an extension of the non-diagonal
phase-field model. Different from the conventional variational non-isothermal
sintering phase-field models, the derived model contains cross-coupling terms
in the diffusion and phase-field evolution equations, which are essential to
ensure the quantitative validity of the model. Furthermore, we demonstrate
that the existence of the trapping effects and presence of surface diffusion
in conservation laws are direction-dependent. It hence highlights the need of
an anisotropic interpolation of the diffusivity tensor.

The paper is structured as follows. The formulations of the quantitative model
(denoted as “quantitative model” hereinafter) are derived in Section II where
the entropy functional and time evolution equations are explicitly given.
Sharp-interface description across solid free surfaces is briefly explained in
Section III. Afterwards, a linkage between model parameters and sharp-
interface equations using a reduction procedure is demonstrated in Section IV.
Section V shows the verification and importance of quantitative model followed
by comparing with the non-isothermal sintering model proposed in our former
work (denoted as “existing model” hereinafter). Conclusions are presented in
Section VI.

##  II Model Formulation

Underlying physical processes involved in non-isothermal sintering can be
classified as but not limited to: (a) the mass/heat transport, including
diffusion through sorts of paths (volume, surface, and grain boundaries) and
mass flows (viscous or fluid flow); (b) the structural relaxation, including
the rigid-body motions of powders and interface (mostly the grain boundaries)
migration. All underlying interactive processes collectively lead to two
significant phenomena: one is the densification (eliminating the pores), in
which the total surface energy should be reduced; the other is the grain
coarsening, in which the total grain-boundary energy should decrease [1, 2,
46]. In the following, we then present a framework for deriving non-isothermal
variational quantitative phase-field sintering model, with its quantitative
validity engendered by asymptotic analysis followed.

###  II.1 Entropy and free energy functionals

In this model, a conserved OP $\rho$ denoting the solid density fraction is
used to indicate the solid region ($\rho=1$) and the atmosphere/pores region
($\rho=0$) while a series of non-conserved OPs $\\{\eta_{i}\\}$ are used to
represent the different grain orientations of the solid grains. Considering a
non-isothermal scenario, the entropy functional $S$ for a subdomain $\Omega$
within the sintering system is defined as

$\begin{split}S(e,\rho,\left\\{\eta_{i}\right\\})=&\int_{\Omega}\bigg{[}\bigg{.}s(e,\rho,\left\\{\eta_{i}\right\\})-\frac{\kappa_{\rho}}{2}|\nabla\rho|^{2}\\\
&-\frac{\kappa_{\eta}}{2}\sum_{i}|\nabla\eta_{i}|^{2}\bigg{.}\bigg{]}\text{d}\Omega,\end{split}$
(1)

with

$s=\frac{1+h(\rho)}{2}s_{\mathrm{ss}}(e_{\mathrm{ss}})+\frac{1-h(\rho)}{2}s_{\mathrm{at}}(e_{\mathrm{at}})+s_{\text{cf}}(\rho,\left\\{\eta_{i}\right\\}),$
(2)

where $s$ is the local entropy density, $e$ is the internal energy density
while $\kappa_{\rho}$ and $\kappa_{\eta}$ are the gradient energy coefficients
associated with $\rho$ and $\\{\eta_{i}\\}$ respectively. $s_{\mathrm{ss}}$ is
the solid phase bulk entropy density and is dependent on the internal energy
density of the solid $e_{\mathrm{ss}}$. The bulk entropy density of the
atmosphere $s_{\mathrm{at}}$ is dependent on the internal energy density of
the atmosphere $e_{\mathrm{at}}$. $h(\rho)=2\rho-1$ is an interpolation
function. The configurational entropy $s_{\text{cf}}$ is related to the
spatial distribution of entropy density proportional to $\rho$ and
$\\{\eta_{i}\\}$. It is formulated in a form of Landau-type polynomial similar
to the one given by Ref. [15] as

$\displaystyle s_{\text{cf}}(\rho,\\{\eta_{i}\\})=$
$\displaystyle\underline{C}_{\text{cf}}\left[\rho^{2}(1-\rho)^{2}\right]\>+\underline{D}_{\text{cf}}\bigg{[}\bigg{.}\rho^{2}\>+6(1-\rho)\sum_{i}\eta^{2}_{i}\>$
(3)
$\displaystyle-4(2-\rho)\sum_{i}\eta^{3}_{i}\>+3\left(\sum_{i}\eta^{2}_{i}\right)^{2}\bigg{.}\bigg{]},$

where $\underline{C}_{\text{cf}}$ and $\underline{D}_{\text{cf}}$ are
constants. The multi-well potential in Eq. (3) can be seen to exhibit minimal
at various regions such as: atmosphere
($\rho=0,\>\\{\eta_{1}=0,\cdots,\eta_{n}=0$}), and solid grains at different
orientations ($\rho=1,\>\\{\eta_{1}=1,\cdots,\eta_{n}=0$}), $\cdots$ ,
($\rho=1,\>\\{\eta_{1}=0,\cdots,\eta_{n}=1$}). One advantage of this potential
form is that its constant parameters can be directly linked to material
properties [17].

Assuming $e$ can be expressed as

$e=\frac{1+h(\rho)}{2}e_{\mathrm{ss}}+\frac{1-h(\rho)}{2}e_{\mathrm{at}}+e_{\text{pt}}(\rho,\left\\{\eta_{i}\right\\}),$
(4)

where $e_{\text{pt}}$ accounts for the spatial distribution of the internal
energy proportional to $\rho$ and $\\{\eta_{i}\\}$ and is also formulated
similar to $s_{\text{cf}}$ as

$\displaystyle e_{\text{pt}}(\rho,\left\\{\eta_{i}\right\\})=$
$\displaystyle\underline{C}_{\text{pt}}\left[\rho^{2}(1-\rho)^{2}\right]\>+\underline{D}_{\text{pt}}\bigg{[}\bigg{.}\rho^{2}\>+6(1-\rho)\sum_{i}\eta^{2}_{i}\>$
(5)
$\displaystyle-4(2-\rho)\sum_{i}\eta^{3}_{i}\>+3\left(\sum_{i}\eta^{2}_{i}\right)^{2}\bigg{.}\bigg{]},$

where $\underline{C}_{\text{pt}}$ and $\underline{D}_{\text{pt}}$ are
constants.

Following the Legendre transformation, we can obtain the free energy
functional $F$ as

$\displaystyle F(T,\rho,\left\\{\eta_{i}\right\\})=$
$\displaystyle\int_{\Omega}\bigg{[}\bigg{.}f(T,\rho,\left\\{\eta_{i}\right\\})+\frac{T\kappa_{\rho}}{2}|\nabla\rho|^{2}$
(6)
$\displaystyle+\frac{T\kappa_{\eta}}{2}\sum_{i}|\nabla\eta_{i}|^{2}\bigg{.}\bigg{]}\text{d}\Omega,$

with

$f(T,\rho,\left\\{\eta_{i}\right\\})=\frac{1+h(\rho)}{2}f_{\mathrm{ss}}(T)+\frac{1-h(\rho)}{2}f_{\mathrm{at}}(T)+e_{\text{pt}}-Ts_{\text{cf}},$
(7)

where $f_{\mathrm{ss}}$ and $f_{\mathrm{at}}$ are the free energy densities of
the solid phase and the atmosphere, respectively. $T$ is the temperature.
Substituting Eqs. (3) and (5) into Eq. (7), we obtain

$\displaystyle f(T,\rho,\left\\{\eta_{i}\right\\})=$
$\displaystyle\frac{1+h(\rho)}{2}f_{\mathrm{ss}}(T)+\frac{1-h(\rho)}{2}f_{\mathrm{at}}(T)$
(8)
$\displaystyle+\underline{C}\left[\rho^{2}(1-\rho)^{2}\right]+\underline{D}\bigg{[}\bigg{.}\rho^{2}+6(1-\rho)\sum_{i}\eta^{2}_{i}\>$
$\displaystyle-4(2-\rho)\sum_{i}\eta^{3}_{i}+3\left(\sum_{i}\eta^{2}_{i}\right)^{2}\bigg{.}\bigg{]},$

with

$\displaystyle\underline{C}$ $\displaystyle=$
$\displaystyle\underline{C}_{\text{pt}}-T\underline{C}_{\text{cf}},$
$\displaystyle\underline{D}$ $\displaystyle=$
$\displaystyle\underline{D}_{\text{pt}}-T\underline{D}_{\text{cf}}.$

###  II.2 Kinetic equations

Considering that $\rho$ and $e$ are conserved OPs, they satisfy mass and
energy conservation laws respectively:

$\dot{\rho}=-\nabla\cdot\mathbf{J}_{\rho},$ (9)
$\dot{e}=-\nabla\cdot\mathbf{J}_{e},$ (10)

where $\mathbf{J}_{\rho}$ is the mass diffusion flux and $\mathbf{J}_{e}$ is
the energy flux.

Following our previous work [8], the non-negative entropy production $\sigma$
in the subdomain can be formulated as

$\displaystyle\sigma=\int_{\Omega}\bigg{[}\mathbf{J}_{\rho}\cdot\nabla\frac{\delta
S}{\delta\rho}+\mathbf{J}_{e}\cdot\nabla\frac{\delta S}{\delta
e}+\sum_{i}\dot{\eta_{i}}\frac{\delta
S}{\delta\eta_{i}}\bigg{]}\text{d}\Omega,$ (11)

with

$\frac{\delta S}{\delta\rho}=-\frac{1}{T}\frac{\delta
F}{\delta\rho},\>\frac{\delta S}{\delta\eta_{i}}=-\frac{1}{T}\frac{\delta
F}{\delta\eta_{i}},\>\frac{\delta S}{\delta e}=\frac{1}{T},$

where $\nabla({\delta S}/{\delta\rho})$ is the driving force associated with
$\mathbf{J}_{\rho}$, $\nabla({\delta S}/{\delta e})$ is the driving force
associated with $\mathbf{J}_{e}$ and ${\delta S}/{\delta\eta_{i}}$ is the
driving force associated with $\dot{\eta_{i}}$.

In the view of the phenomenological linear laws of non-equilibrium
thermodynamics and also ensuring non-negative production of the entropy, we
can define the relationships between the fluxes, the non-conserved OPs time
evolution equations and their driving forces as

$\begin{bmatrix}\mathbf{J}_{\rho}\\\ \mathbf{J}_{e}\\\ \dot{\eta_{1}}\\\
\vdots\\\
\dot{\eta_{n}}\end{bmatrix}=\begin{bmatrix}\mathbf{L}_{\rho\rho}&\mathbf{L}_{\rho
e}&\mathbf{L}_{\rho\eta_{1}}&\cdots&\mathbf{L}_{\rho\eta_{n}}\\\
\mathbf{L}_{e\rho}&\mathbf{L}_{ee}&\mathbf{L}_{e\eta_{1}}&\cdots&\mathbf{L}_{e\eta_{n}}\\\
\mathbf{L}_{\eta_{1}\rho}&\mathbf{L}_{\eta_{1}e}&\mathbf{L}_{\eta_{1}\eta_{1}}&\cdots&\mathbf{L}_{\eta_{1}\eta_{n}}\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
\mathbf{L}_{\eta_{n}\rho}&\mathbf{L}_{\eta_{n}e}&\mathbf{L}_{\eta_{n}\eta_{1}}&\cdots&\mathbf{L}_{\eta_{n}\eta_{n}}\end{bmatrix}\begin{bmatrix}-\nabla(\frac{\mu}{T})\\\
\nabla(\frac{1}{T})\\\ \frac{\delta S}{\delta\eta_{1}}\\\ \vdots\\\
\frac{\delta S}{\delta\eta_{n}}\end{bmatrix},$ (12)

where $\mu={\delta F}/{\delta\rho}$ is defined as the chemical potential and
$n$ represents the total number of grain orientations.
$\mathbf{L}_{\rho\rho}$, $\mathbf{L}_{\rho e}$, $\mathbf{L}_{e\rho}$ and
$\mathbf{L}_{ee}$ are positively defined rank 2 tensors and for
$i=1,2,\cdots,n$ , $\mathbf{L}_{\rho\eta_{i}}$, $\mathbf{L}_{e\eta_{i}}$,
$\mathbf{L}_{\eta_{i}\rho}$ and $\mathbf{L}_{\eta_{i},e}$ are positively
defined rank 1 tensors while $\mathbf{L}_{\eta_{i}\eta_{i}}$ is a positively
defined rank 0 tensor.

Diagonal terms $\mathbf{L}_{\rho\rho}$ and $\mathbf{L}_{ee}$ are the
diffusional mobilities of mass and energy diffusion respectively. Mobility
term associated with the grain orientations, $\mathbf{L}_{\eta_{i}\eta_{i}}$
is simply a scalar function and is thereafter taken as $L_{\eta}$ where we
assume isotropic condition taking it to be the same regardless of $i$. The
non-diagonal terms in the Onsager matrix in Eq. (12) represent cross-couplings
between the various OPs. Based on the Onsager reciprocal relations, we have
$\mathbf{L}_{\rho e}=\mathbf{L}_{e\rho}$,
$\mathbf{L}_{\rho\eta_{i}}=\mathbf{L}_{\eta_{i}\rho}$ and
$\mathbf{L}_{e\eta_{i}}=\mathbf{L}_{\eta_{i}e}$. Note that the cross-coupling
between the different grain orientations is not considered, resulting in
similar $\dot{\eta}_{i}$ formulation for all $i$. Hence we consider only one
$\dot{\eta}_{i}$ whose formulation is representative for all $i$. The
quantities $\mathbf{L}_{\rho e}$ and $\mathbf{L}_{e\rho}$ are associated with
the mass flux due to temperature gradient (thermophoresis effect) and with the
energy flux due to chemical potential gradient (Dufour effect), respectively.
Examination of these effects has been done in our previous work [8] and is not
the main priority of this work. Therefore, the terms associated with
$\mathbf{L}_{\rho e}$ and $\mathbf{L}_{e\rho}$ in the fluxes formulations are
dropped. The time evolution equations can then be written as

$\dot{\rho}=\nabla\cdot\bigg{[}\mathbf{L}_{\rho\rho}\cdot\nabla\bigg{(}\frac{\mu}{T}\bigg{)}\bigg{]}+\nabla\cdot\bigg{[}\frac{1}{T}\sum_{i}\mathbf{L}_{\rho\eta_{i}}\frac{\delta
F}{\delta\eta_{i}}\bigg{]},$ (13a)
$\dot{e}=\nabla\cdot\bigg{[}\mathbf{L}_{ee}\cdot\frac{\nabla
T}{T^{2}}\bigg{]}+\nabla\cdot\bigg{[}\frac{1}{T}\sum_{i}\mathbf{L}_{e\eta_{i}}\frac{\delta
F}{\delta\eta_{i}}\bigg{]},$ (13b)
$\dot{\eta}_{i}=-\mathbf{L}_{\eta_{i}\rho}\cdot\nabla\bigg{(}\frac{\mu}{T}\bigg{)}-\mathbf{L}_{\eta_{i},e}\cdot\frac{\nabla
T}{T^{2}}-L_{\eta}\frac{1}{T}\frac{\delta F}{\delta\eta_{i}}.$ (13c)

Formulations expressed in Eqs. (12) and (13) present the fluxes and time
evolution equations of the associated OPs in terms of the driving forces.
However, for consistency with previous non-diagonal models [38, 40] as well as
ease of relating our model to the sharp-interface counterpart as will be
discussed later, we reformulate the phenomenological linear relations
employing the linear relations of the driving forces in terms of the fluxes
and time evolution equations such that

$-\nabla\left(\frac{\mu}{T}\right)=\mathbf{L}_{\rho\rho}^{-1}\cdot\mathbf{J}_{\rho}+\mathbf{L}_{\rho\eta_{i}}^{-1}\sum_{i}\dot{\eta_{i}},$
(14a) $-\frac{\nabla
T}{T^{2}}=\mathbf{L}_{ee}^{-1}\cdot\mathbf{J}_{e}+\mathbf{L}_{e\eta_{i}}^{-1}\sum_{i}\dot{\eta_{i}},$
(14b) $-\frac{1}{T}\frac{\delta
F}{\delta\eta_{i}}=\mathbf{L}_{\eta_{i}\rho}^{-1}\cdot\mathbf{J}_{\rho}+\mathbf{L}_{\eta_{i},e}^{-1}\cdot\mathbf{J}_{e}+L_{\eta}^{-1}\dot{\eta_{i}}.$
(14c)

Since the variation of mass density and internal energy is found across free
surfaces of the solid grains, the cross-coupling terms
$\mathbf{L}^{-1}_{\rho\eta_{i}}=\mathbf{L}^{-1}_{\eta_{i}\rho}$ and
$\mathbf{L}^{-1}_{e\eta_{i}}=\mathbf{L}^{-1}_{\eta_{i}e}$ should be defined
such that they are only evaluated at the free surfaces. Also, the non-
equilibrium effects associated with these cross terms scale with the diffuse
interface width $l$. Accordingly, following Refs. [37, 38, 40], we propose the
following formulations:

$\mathbf{L}_{\rho\eta_{i}}^{-1}=\mathbf{L}_{\eta_{i}\rho}^{-1}={M}_{1}(\rho)l\nabla\rho,$
(15a)
$\mathbf{L}_{e\eta_{i}}^{-1}=\mathbf{L}_{\eta_{i}e}^{-1}={M}_{2}(\rho)l\nabla\rho,$
(15b)

where ${M}_{1}$ and ${M}_{2}$ are scalar functions used to parametrize the
associated cross-coupling terms. $l\nabla\rho$ is a vector normal to the free
surfaces and has a magnitude of 1 at the center of the free surfaces assuming
the parameter $\alpha$ used to adjust the definition of $l$ in Ref. [47]
equals 2 [28]. Substituting Eq. (14) into (11) and taking into account the
aforementioned, we obtain the entropy production in the subdomain as

$\displaystyle\sigma=$
$\displaystyle\int_{\Omega}\bigg{[}\mathbf{L}_{\rho\rho}^{-1}\cdot\mathbf{J}_{\rho}\cdot\mathbf{J}_{\rho}+\mathbf{L}_{ee}^{-1}\cdot\mathbf{J}_{e}\cdot\mathbf{J}_{e}+L_{\eta}^{-1}\bigg{(}\sum_{i}\dot{\eta_{i}}\bigg{)}^{2}$
(16)
$\displaystyle+2l\nabla\rho\sum_{i}\dot{\eta_{i}}\cdot\left(M_{1}\mathbf{J}_{\rho}+M_{2}\mathbf{J}_{e}\right)\bigg{]}\text{d}\Omega.$

Furthermore, time evolution equations can be obtained as

$\dot{\rho}=\nabla\cdot\bigg{[}\mathbf{L}_{\rho\rho}\cdot\bigg{(}\nabla\left(\frac{\mu}{T}\right)\>+M_{1}l\nabla\rho\sum_{i}\dot{\eta_{i}}\bigg{)}\bigg{]},$
(17a) $\displaystyle c_{\text{r}}\dot{T}+\frac{\partial
e}{\partial\rho}\dot{\rho}+\sum_{i}\frac{\partial
e}{\partial\eta_{i}}\dot{\eta_{i}}=$
$\displaystyle\nabla\cdot\bigg{[}\mathbf{L}_{ee}\cdot\bigg{(}\frac{\nabla
T}{T^{2}}$ (17b)
$\displaystyle+M_{2}l\nabla\rho\sum_{i}\dot{\eta_{i}}\bigg{)}\bigg{]},$
$\displaystyle\hat{L}_{\eta}^{-1}\dot{\eta_{i}}=$
$\displaystyle\>\kappa_{\eta}\nabla^{2}\eta_{i}-\frac{1}{T}\frac{\partial
f}{\partial\eta_{i}}\>+\>l\nabla\rho\cdot\bigg{[}M_{1}\mathbf{L}_{\rho\rho}\cdot\nabla\left(\frac{\mu}{T}\right)\>$
(17c) $\displaystyle+M_{2}\mathbf{L}_{ee}\cdot\frac{\nabla T}{T^{2}}\bigg{]},$

with

$\hat{L}_{\eta}^{-1}=L_{\eta}^{-1}-[M^{2}_{1}l^{2}\nabla\rho\cdot\mathbf{L}_{\rho\rho}+M_{2}^{2}l^{2}\nabla\rho\cdot\mathbf{L}_{ee}]\cdot\nabla\rho.$
(18)

Hereby
$c_{\text{r}}=\frac{1+h(\rho)}{2}c_{\mathrm{ss}}+\frac{1-h(\rho)}{2}c_{\mathrm{at}}$
is the relative specific heat, where $c_{\mathrm{ss}}={\partial
e_{\mathrm{ss}}}/{\partial T}$ and $c_{\mathrm{at}}={\partial
e_{\mathrm{at}}}/{\partial T}$ are the volumetric specific heat of solid and
atmosphere, respectively.

Comparing the heat transfer equation (Eq. (17b)) to that of conventional
quantitative phase-field model [36], the second term on the right-hand-side
(RHS) can be likened to the thermal antitrapping current related to the
elimination of thermal trapping (associated with temperature jump) at the free
surfaces. Similarly, the second term on the RHS of Eq. (17a) represents some
form of antitrapping current valued only at the free surfaces. Similar to
solutal antitrapping current [34, 48] associated with solute trapping due to
jump of chemical potential, this term is termed as the mass antitrapping
current in this work. The last two terms on the RHS of the grain orientation
time evolution equations (Eq. (17c)) represent cross-coupling terms associated
with mass and energy diffusion across the free surfaces, respectively. These
terms alongside the antitrapping terms are absent in time evolution equations
of conventional non-isothermal phase-field sintering models but are very vital
in the elimination of artificial interface effects such as the trapping
effects at the free surfaces of the solid phase.

Moreover, considering no variation of solid density and thermal properties
across the grain boundaries, Eq. (17c) has no cross-coupling terms and simply
takes a form of Allen-Cahn equation at the grain boundaries. Consequently, we
limit our subsequent analysis and derivations to the free surfaces where the
cross-coupling terms are significant.

##  III Sharp-interface description across free surfaces

Considering a simple nonisothermal system consisting a sharp free surface
between a solid grain and the atmosphere, the following set of sharp-interface
equations can be described in the bulk regions:

$\frac{\partial\rho}{\partial t}=\nabla\cdot({M}_{\text{rg}}\nabla\mu),$ (19)
$c_{\text{rg}}\frac{\partial T}{\partial t}=\nabla\cdot({k}_{\text{rg}}\nabla
T),$ (20)

where for a bulk region rg (”$\mathrm{ss}$” for solid and ”$\mathrm{at}$” for
atmosphere), ${M}_{\text{rg}}$, $c_{\text{rg}}$, and ${k}_{\text{rg}}$
represent the region’s effective mass mobility coefficient, volumetric
specific heat, and effective thermal conductivity respectively. $\rho$ here
adopts the physical meaning of normalised density of the solid. Eqs. (19) and
(20) describe mass and heat transfer in the bulk regions. For the bulk
atmosphere region in particular, ${M}_{\mathrm{at}}$ describes the effective
mobility considering mass transfer mechanisms notably evaporation and
condensation. Hence, the driving force $\nabla\mu$ for mass transfer in the
atmosphere takes into account vapor pressure differences due to local
curvature [1, 2]. Similarly, ${k}_{\text{at}}$ describes effective thermal
conductivity taking into account convection and radiation.

Furthermore, energy conservation condition at the free surface can be
described as

$\displaystyle ve_{\mathrm{ss}}+{k}_{\mathrm{ss}}\left.\nabla
T\right|_{\mathrm{ss}}\cdot\mathbf{n}_{\text{sf}}=ve_{\mathrm{at}}+{k}_{\mathrm{at}}\left.\nabla
T\right|_{\mathrm{at}}\cdot\mathbf{n}_{\text{sf}}=J_{T},$ (21)

where $v$ is the velocity of the migrating free surface, and $\left.\nabla
T\right|_{\mathrm{ss}}$ and $\left.\nabla T\right|_{\mathrm{at}}$ are the
spatial gradients of the temperature at the solid and atmosphere sides of the
free surface respectively. $\mathbf{n}_{\text{sf}}$ is the unit vector normal
to the free surface. $J_{T}$ is the normal heat flux flowing through the free
surface. Similarly, explicit formulation of mass conservation at the free
surface is given as

$\begin{split}v(\rho_{\mathrm{ss}}-\rho_{\mathrm{at}})=&-{M}_{\mathrm{ss}}\left.\nabla\mu\right|_{\mathrm{ss}}\cdot\mathbf{n}_{\text{sf}}+{M}_{\mathrm{at}}\left.\nabla\mu\right|_{\mathrm{at}}\cdot\mathbf{n}_{\text{sf}}\\\
&+{M}_{\text{sf}}\nabla^{2}_{\text{sf}}\mu,\end{split}$ (22)

where $\rho_{\mathrm{ss}}$ and $\rho_{\mathrm{at}}$ are the bulk densities in
the solid and atmosphere, and $\left.\nabla\mu\right|_{\mathrm{ss}}$ and
$\left.\nabla\mu\right|_{\mathrm{at}}$ are the spatial gradients of the
chemical potential at the solid and atmosphere sides of the free surface
respectively. $M_{\text{sf}}$ represents surface diffusion mobility.
$\nabla_{\mathrm{sf}}^{2}$ is surface Laplacian. The last term in Eq. (22)
describes surface diffusion typical to sharp-interface description of mass
transfer in sintering [18]. Moreover, $v$ can be defined as

$v=v_{s}+v_{b},$ (23)

where $v_{s}$ and $v_{b}$ are the velocities contributed by surface diffusion
and bulk/volume diffusion respectively and can be expressed in terms of their
corresponding mass fluxes;

$v_{s}=-V_{\mathrm{m}}\nabla_{\mathbf{sf}}\cdot\mathbf{J}_{\text{sf}},\;v_{b}=-V_{\mathrm{m}}\mathbf{J}_{b}\cdot\mathbf{n}_{\text{sf}},$
(24)

where $V_{\mathrm{m}}$ is the molar volume and $\nabla_{\mathbf{sf}}$ is the
surface gradient. $\mathbf{J}_{\text{sf}}$ is the mass flux along the free
surface associated with surface gradient of the free surface curvature
$\mathrm{k}_{\text{sf}}$;
$\mathbf{J}_{\text{sf}}\propto\nabla_{\mathbf{sf}}\mathrm{k}_{\text{sf}}$.
$\mathbf{J}_{b}$ is mass flux from the solid bulk to the free surface
associated with the gradient of the chemical potential in the solid bulk
grains $\mu_{\mathrm{ss}}$; $\mathbf{J}_{b}\propto\nabla\mu_{\mathrm{ss}}$,
[18, 49].

In addition, the chemical potential and temperature at the free surface obey
the following relations:

$\mu|_{\mathrm{ss}}=\mu|_{\mathrm{at}},$ (25a)
$T|_{\mathrm{ss}}-T|_{\mathrm{at}}=J_{T}R_{s},$ (25b)

where $\mu|_{\mathrm{ss}}$ and $\mu|_{\mathrm{at}}$ represent chemical
potentials at the solid and atmosphere sides of the free surface,
respectively. $T|_{\mathrm{ss}}$ and $T|_{\mathrm{at}}$ represent the
temperatures at the solid and atmosphere sides of the free surface,
respectively. $R_{s}$ represents Kapitza-type thermal resistance. In this
work, we assume negligible $R_{s}$, thereby Eq. (25a) and Eq. (25b) indicate
imposed zero chemical potential and temperature jumps at the free surface.

Moreover, we infer that jump in chemical potential $\delta\mu$ across the free
surface is conjugated to $v$ and also that the temperature jump $\delta T$
across the free surface is conjugated to $J_{T}$. The kinetic boundary
conditions can then be expressed in the framework of phenomenological linear
relations as [50, 51]

$\delta\mu=\mathcal{A}v+\mathcal{B}J_{T},$ (26) $\delta
T=\mathcal{B}v+\mathcal{C}J_{T},$ (27)

where $\mathcal{A},\mathcal{B}$ and $\mathcal{C}$ are kinetic coefficients of
the positive-definite Onsager matrix. Entropy production at the free surface
$\sigma_{s}$ can be formulated as

$\sigma_{s}=v\delta\mu\>+J_{T}\delta T.$ (28)

Substituting Eqs. (26) and (27) into (28), we obtain

$\sigma_{s}=\mathcal{A}v^{2}+\mathcal{C}J_{T}^{2}+2\mathcal{B}vJ_{T}.$ (29)

##  IV Thin interface limit: Linking model with sharp-interface description

In this section, still considering a system consisting a free surface between
a solid grain and the atmosphere, we establish the relationships between
$\mathcal{A}$, $\mathcal{B}$ and $\mathcal{C}$ and the phase-field parameters
following the reduction procedure presented in Ref. [38]. Considering a 1D
system with the free surface centered at $x=0$ (shown in Fig. 1a), we have
$\rho$ and $\eta$ vary from a semi-finite solid region ($-\infty$) to a semi-
finite atmosphere region ($+\infty$). For simplicity, notation
$(\cdot)^{\prime}$ is adopted to represent the derivative w.r.t. the spatial
coordinate $x$. It is worth noting that we consider the profile of $\rho$
between two bulk values that are slightly deviated from the ideal ones, i.e.,
$\rho_{\mathrm{ss}}$ in the substance and $\rho_{\mathrm{at}}$ in the
pore/atmosphere. The origin and the thermodynamic outcome of these deviated
bulk values of $\rho$ are explicitly examined and discussed in the Appendix.

Figure 1:  (a) Asymptotic schematic across a free surface; blue solid line
represents phase-field profile at equilibrium and red large-dashed line
represents sharp-interface profile; blue dotted line shows the profile of
$\rho$ with slightly deviated values; (b) Asymptotic schematic across a grain
boundary; blue solid and dashed lines represent phase-field profiles and red
large-dashed and dotted lines represent sharp-interface profiles.

According to the phase-field method, the entropy production (Eq. (16)) for the
system considered can be formulated as

$\displaystyle\sigma=$
$\displaystyle\int^{-l/2}_{-\infty}\bigg{[}\frac{J_{\rho}^{2}(x)}{{L}^{\mathrm{ss}}_{\rho\rho}}+\frac{J_{e}^{2}(x)}{{L}^{\mathrm{ss}}_{ee}}\bigg{]}\mathop{}\\!\mathrm{d}x+\int^{\infty}_{l/2}\bigg{[}\frac{J_{\rho}^{2}(x)}{{L}^{\mathrm{at}}_{\rho\rho}}+\frac{J_{e}^{2}(x)}{{L}^{\mathrm{at}}_{ee}}\bigg{]}\mathop{}\\!\mathrm{d}x$
$\displaystyle+\int^{l/2}_{-l/2}\bigg{[}\frac{J_{\rho}^{2}(x)}{{L}_{\rho\rho}}+\frac{J_{e}^{2}(x)}{{L}_{ee}}+L_{\eta,\text{sf}}^{-1}\dot{\eta}^{2}+2l\rho^{{}^{\prime}}(x)\dot{\eta}\left(M_{1}J_{\rho}(x)+M_{2}J_{e}(x)\right)\bigg{]}\mathop{}\\!\mathrm{d}x,$

where ${L}^{\mathrm{ss}}_{\rho\rho}={L}_{\rho\rho}(\rho=\rho_{\mathrm{ss}})$
and ${L}^{\mathrm{at}}_{\rho\rho}={L}_{\rho\rho}(\rho=\rho_{\mathrm{at}})$ are
the effective mass mobilities in the corresponding regions. Also,
${L}^{\mathrm{ss}}_{ee}={L}_{ee}(\rho=\rho_{\mathrm{ss}})$ and
${L}^{\mathrm{at}}_{ee}={L}_{ee}(\rho=\rho_{\mathrm{at}})$ are the effective
energy mobilities in the corresponding regions. $L_{\eta,\text{sf}}$ is
mobility of $\eta$ at the free surface. It can be noted that in the bulk
regions $(|x|>l/2)$, only the fluxes $J_{\rho}(x)$ and $J_{e}(x)$ contribute
to entropy production as $\dot{\eta}$ and $\rho^{{}^{\prime}}(x)$ both go to
zero. The entropy production of the system considered can be formulated within
the sharp-interface description as

$\int^{0}_{-\infty}\bigg{[}\frac{J_{\rho}^{2}(x)}{{L}^{\mathrm{ss}}_{\rho\rho}}+\frac{J_{e}^{2}(x)}{{L}^{\mathrm{ss}}_{ee}}\bigg{]}\mathop{}\\!\mathrm{d}x+\int^{\infty}_{0}\bigg{[}\frac{J_{\rho}^{2}(x)}{{L}^{\mathrm{at}}_{\rho\rho}}+\frac{J_{e}^{2}(x)}{{L}^{\mathrm{at}}_{ee}}\bigg{]}\mathop{}\\!\mathrm{d}x+\sigma_{s}.$
(31)

Comparing Eqs. (IV) and (31), we obtain entropy production at the free surface
within the phase-field model as

$\displaystyle\sigma_{s}=$
$\displaystyle\int^{l/2}_{-l/2}\bigg{[}\frac{J_{\rho}^{2}(x)}{{L}_{\rho\rho}}+\frac{J_{e}^{2}(x)}{{L}_{ee}}+L_{\eta,\text{sf}}^{-1}\dot{\eta}^{2}+2l\rho^{{}^{\prime}}(x)\dot{\eta}\bigg{(}M_{1}J_{\rho}(x)+M_{2}J_{e}(x)\bigg{)}\bigg{]}\mathop{}\\!\mathrm{d}x$
(32)
$\displaystyle-\int^{0}_{-l/2}\bigg{[}\frac{J_{\mathrm{ss}(\rho)}^{2}}{{L}^{\mathrm{ss}}_{\rho\rho}}+\frac{J_{\mathrm{ss}(e)}^{2}}{{L}^{\mathrm{ss}}_{ee}}\bigg{]}\mathop{}\\!\mathrm{d}x-\int^{l/2}_{0}\bigg{[}\frac{J_{\mathrm{at}(\rho)}^{2}}{{L}^{\mathrm{at}}_{\rho\rho}}+\frac{J_{\mathrm{at}(e)}^{2}}{{L}^{\mathrm{at}}_{ee}}\bigg{]}\mathop{}\\!\mathrm{d}x,$

where for a region rg (“$\mathrm{ss}$” for solid and “$\mathrm{at}$” for
atmosphere), $J_{\text{rg}(\rho)}$ and $J_{\text{rg}(e)}$ represent the
region’s bulk mass and energy fluxes respectively.

For the purpose of making direct relations between Eq. (32) and its sharp-
interface counterpart, Eq. (29), we express $J_{\rho}(x)$, $J_{e}(x)$ and
$\dot{\eta}$ in terms of $v$ and $J_{T}$. First, we make analysis considering
only fluxes that are flowing through the free surface along $x$ direction (i.e
normal to the free surface) thereby we tentatively drop the contribution of
the surface diffusion flux $J_{\text{sf}}$ to $v$ since it is tangential to
the free surface. Second, we employ a quasisteady approximation that assumes
large gradients of $\rho$, $e$ and $\eta$ across the free surface such that we
define their time derivatives as:

$\dot{\rho}\approx-v\rho^{{}^{\prime}}(x),\>\dot{e}\approx-
ve^{{}^{\prime}}(x),\>\dot{\eta}\approx-v\eta^{{}^{\prime}}(x).$ (33)

We integrate both sides of the conservation laws,
$\dot{e}=-J_{e}^{{}^{\prime}}(x)$ and $\dot{\rho}=-J_{\rho}^{{}^{\prime}}(x)$
after substituting Eq. (33);

$\int_{J_{\mathrm{ss}(\rho)}}^{J_{\mathrm{at}(\rho)}}\mathop{}\\!\mathrm{d}J_{\rho}=v\int^{\rho_{\text{eq}}^{\mathrm{ss}}}_{\rho_{\text{eq}}^{\mathrm{at}}}\mathop{}\\!\mathrm{d}\rho,\quad\int_{J_{\mathrm{ss}(e)}}^{J_{\mathrm{at}(e)}}\mathop{}\\!\mathrm{d}J_{e}=v\int_{e_{\text{eq}}^{\mathrm{ss}}}^{e_{\text{eq}}^{\mathrm{at}}}\mathop{}\\!\mathrm{d}e$
(34)

with the boundary values as

$\begin{split}J_{\mathrm{ss}(\rho)}\approx
v\rho_{\text{eq}}^{\mathrm{ss}},&\quad J_{\mathrm{at}(\rho)}\approx
v\rho_{\text{eq}}^{\mathrm{at}},\\\ J_{\mathrm{ss}(e)}\approx
ve^{\mathrm{ss}}_{\text{eq}}-J_{T},&\quad J_{\mathrm{at}(e)}\approx
ve^{\mathrm{at}}_{\text{eq}}-J_{T},\end{split}$ (35)

where $\rho^{\text{eq}}_{\text{rg}}$ and $e^{\mathrm{at}}_{\text{rg}}$
($\text{rg}=\mathrm{ss},~{}\mathrm{at}$) are the equilibrium conserved OP and
internal energies, respectively. The integrals in Eq.(34) yield

$J_{e}(x)\approx ve(x)-J_{T},\quad J_{\rho}(x)\approx v\rho(x).$ (36)

Furthermore, we adopt the sigmoid formulation for the profiles of $\rho$(x)
and $\eta(x)$ in this work as

$\displaystyle\rho(x)$
$\displaystyle=\frac{1}{2}\left[\left(\rho_{\mathrm{ss}}+\rho_{\mathrm{at}}\right)+\left(\rho_{\mathrm{ss}}-\rho_{\mathrm{at}}\right)\tanh\frac{2x}{l}\right],$
(37) $\displaystyle\eta(x)$
$\displaystyle=\frac{1}{2}\left[1+\tanh\left(\frac{2x}{l}\right)\right]$ (38)

with the diffuse interface width $l$. Taking into account all the
aforementioned, we obtain $\sigma_{s}$ to be

$\displaystyle\sigma_{s}=$
$\displaystyle\int^{l/2}_{-l/2}\bigg{[}\frac{(v\rho(x))^{2}}{{L}_{\rho\rho}}-\frac{(v\rho_{\mathrm{ss}})^{2}}{2{L}^{\mathrm{ss}}_{\rho\rho}}-\frac{(v\rho_{\mathrm{at}})^{2}}{2{L}^{\mathrm{at}}_{\rho\rho}}\bigg{]}\mathop{}\\!\mathrm{d}x+\int^{l/2}_{-l/2}\bigg{[}\frac{(ve(x)-J_{T})^{2}}{{L}_{ee}}-\frac{(ve_{\mathrm{ss}}-J_{T})^{2}}{2{L}^{\mathrm{ss}}_{ee}}-\frac{(ve_{\mathrm{at}}-J_{T})^{2}}{2{L}^{\mathrm{at}}_{ee}}\bigg{]}\mathop{}\\!\mathrm{d}x$
(39)
$\displaystyle-\int^{l/2}_{-l/2}4l\rho^{{}^{\prime}}(x)\eta^{{}^{\prime}}(x)v\bigg{[}M_{1}(v\rho(x))+M_{2}(ve(x)-J_{T})\bigg{]}\mathop{}\\!\mathrm{d}x$
$\displaystyle+\int^{l/2}_{-l/2}L_{\eta,\text{sf}}^{-1}v^{2}(\eta_{\text{eq}}^{{}^{\prime}}(x))^{2}\mathop{}\\!\mathrm{d}x.$

It should be noted that the integration range of Eq. (39) can also be taken
from $-\infty$ and $+\infty$ without $\sigma_{s}$ changing. In this regard, we
extend the integration interval from $[-l/2,+l/2]$ to $[-\infty,+\infty]$ in
the following discussion. Comparing Eqs. (29) and (39), we obtain:

$\displaystyle\mathcal{A}=$
$\displaystyle\int^{\infty}_{-\infty}\bigg{[}\frac{\rho^{2}(x)}{{L}_{\rho\rho}}-\frac{(\rho_{\mathrm{ss}})^{2}}{2{L}^{\mathrm{ss}}_{\rho\rho}}-\frac{(\rho_{\mathrm{at}})^{2}}{2{L}^{\mathrm{at}}_{\rho\rho}}\bigg{]}\mathop{}\\!\mathrm{d}x-4\int^{\infty}_{-\infty}M_{1}l\rho^{{}^{\prime}}(x)\eta^{{}^{\prime}}(x)\rho(x)\mathop{}\\!\mathrm{d}x$
(40)
$\displaystyle+\int^{\infty}_{-\infty}\bigg{[}\frac{e^{2}(x)}{{L}_{ee}}-\frac{(e_{\mathrm{ss}})^{2}}{2{L}^{\mathrm{ss}}_{ee}}-\frac{(e_{\mathrm{at}})^{2}}{2{L}^{\mathrm{at}}_{ee}}\bigg{]}\mathop{}\\!\mathrm{d}x-4\int^{\infty}_{-\infty}M_{2}l\rho^{{}^{\prime}}(x)\eta^{{}^{\prime}}(x)e(x)\mathop{}\\!\mathrm{d}x$
$\displaystyle+\int^{\infty}_{-\infty}L_{\eta,\text{sf}}^{-1}[\eta^{{}^{\prime}}(x)]^{2}\mathop{}\\!\mathrm{d}x,$

$\mathcal{B}=\int^{\infty}_{-\infty}2M_{2}l\rho^{{}^{\prime}}(x)\eta^{{}^{\prime}}(x)\mathop{}\\!\mathrm{d}x-\int^{\infty}_{-\infty}\bigg{[}\frac{e(x)}{{L}_{ee}}-\frac{e^{\mathrm{ss}}}{2{L}^{\mathrm{ss}}_{ee}}-\frac{e^{\mathrm{at}}}{2{L}^{\mathrm{at}}_{ee}}\bigg{]}\mathop{}\\!\mathrm{d}x,$
(41)

$\mathcal{C}=\int^{\infty}_{-\infty}\bigg{[}\frac{1}{{L}_{ee}}-\frac{1}{2{L}^{\mathrm{ss}}_{ee}}-\frac{1}{2{L}^{\mathrm{at}}_{ee}}\bigg{]}\mathop{}\\!\mathrm{d}x,$
(42)

The explicit formulations of $\mathcal{A}$, $\mathcal{B}$, and $\mathcal{C}$
implies that the phase-field parameters can be carefully tuned so as to obtain
$\mathcal{A}=0$, $\mathcal{B}=0$ and $\mathcal{C}=0$ which guarantees
$\delta\mu=0$ and $\delta T=0$ across a migrating free surface. However, it is
important to note that even if $\delta\mu=\delta T=0$ is guaranteed, Almgren
[30] showed that for phase-field models, conservation law at interfaces with
opposing phases having asymmetry mobility coefficients has in existence two
additional terms: interface stretching term and surface diffusion term. In the
sintering system, interface stretching represents excess mass and internal
energy along the arclength of the free surfaces of the solid phase [48] and
these excesses can both be respectively eliminated if
$\int^{\infty}_{-\infty}\mathop{}\\!\mathrm{d}x[\rho_{\text{eq}}-\rho^{\mathrm{ss}}_{\text{eq}}/2-\rho^{\mathrm{at}}_{\text{eq}}/2]=0$
and
$\int^{\infty}_{-\infty}\mathop{}\\!\mathrm{d}x[e_{\text{eq}}-e^{\mathrm{ss}}_{\text{eq}}/2-e^{\mathrm{at}}_{\text{eq}}/2]=0$
[30, 40]. Taking $\rho_{\text{eq}}$ as defined in Eq. (37) ensures that the
interface excess of $\rho_{\text{eq}}$ is eliminated. Also, the interface
excess of $e_{\text{eq}}$ is eliminated if $h(\rho)$ is taken as an odd
function. Furthermore, surface diffusion terms in the mass and energy
conservation laws at the free surfaces of the solid are respectively
parameterized by the mobilities
$L^{\mathrm{sf}}_{\rho\rho}=\int^{\infty}_{-\infty}\mathop{}\\!\mathrm{d}x[L_{\rho\rho}(\rho)-L_{\rho\rho}^{\mathrm{ss}}/2-L_{\rho\rho}^{\mathrm{at}}/2]$
and
$L^{\mathrm{sf}}_{\text{ee}}=\int^{\infty}_{-\infty}\mathop{}\\!\mathrm{d}x[L_{\text{ee}}(\rho)-L_{\text{ee}}^{\mathrm{ss}}/2-L_{\text{ee}}^{\mathrm{at}}/2]$.
[30, 40].

In order to make $\delta T=0$, we need to ensure that
$\mathcal{B}=\mathcal{C}=0$. Consequently, $L_{\text{ee}}$ should be
formulated such that it gives the bulk region energy mobilities at the
corresponding regions, ensures $\mathcal{C}=0$ and also guarantee that the
model replicates the sharp-interface energy conservation law (Eq. (21)) where
there is no surface diffusion effect (i.e $L^{\mathrm{sf}}_{\text{ee}}=0$). To
achieve this, Almgren [30] proposed a mobility interpolation function which is
a combination of odd functions with parameters adjusted relative to the bulk
mobilities. This method is contended by Ohno et al. [44] as the mobility
interpolation function produces a non-monotonic function and also contributes
to a limited ratio of the possible bulk mobilities. Nevertheless, it is vital
to note that while simultaneous elimination of $\delta T$ and surface
diffusion effect somewhat put constraints on a scalar formulation of
$L_{\text{ee}}$, the emergence of both effects is actually direction dependent
[52]. $L_{\text{ee}}$ formulation constraint due to $\delta T$ (Eq. (42))
emerges under the consideration of flux components normal to the free surfaces
as seen in the analysis done above while the integral associated with the
surface diffusion effect modification of energy conservation emanates due to
consideration of flux components in tangential direction to the free surfaces
[52]. Therefore, ensuring $\mathcal{C}=0$ and eliminating surface diffusion
term in energy conservation equation are respectively pertinent only at the
normal and tangential directions of the free surfaces. Considering all the
aforementioned and also taking into account the physical context of the energy
mobility, we propose an anisotropic $\mathbf{L}_{ee}$ for the full sintering
description and relate it to the anisotropic thermal conductivity as

$\begin{split}\mathbf{L}_{\text{ee}}&=\left[k_{\perp}\mathbf{N}_{\mathrm{sf}}+k_{\parallel}\mathbf{T}_{\text{sf}}+k_{\text{gb}}\mathbf{T}_{\text{gb}}\right]T^{2}\\\
&=L^{\perp}_{\text{ee}}\mathbf{N}_{\mathrm{sf}}+L^{\parallel}_{\text{ee}}\mathbf{T}_{\text{sf}}+L^{\text{gb}}_{\text{ee}}\mathbf{T}_{\text{gb}},\end{split}$
(43)

with

$k_{\perp}=\left[\frac{1+g(\rho)}{2{k}_{\mathrm{ss}}}+\frac{1-g(\rho)}{2{k}_{\mathrm{at}}}\right]^{-1},$
(44)
$k_{\parallel}=\frac{1+g(\rho)}{2}k_{\mathrm{ss}}+\frac{1-g(\rho)}{2}k_{\mathrm{at}},$
(45) $k_{\text{gb}}=16\sum_{i\neq j}\eta_{i}^{2}\eta_{j}^{2}k_{\text{gb}},$
(46)

and

$\displaystyle\mathbf{N}_{\mathrm{sf}}$
$\displaystyle=\mathbf{n}_{\text{sf}}\otimes\mathbf{n}_{\text{sf}},$ (47)
$\displaystyle\mathbf{T}_{\mathrm{sf}}$
$\displaystyle=\mathbf{I}-\mathbf{n}_{\text{sf}}\otimes\mathbf{n}_{\text{sf}},$
$\displaystyle\mathbf{T}_{\mathrm{gb}}$
$\displaystyle=\mathbf{I}-\mathbf{n}_{\text{gb}}\otimes\mathbf{n}_{\text{gb}}.$

In Eq (43), $L^{\perp}_{\text{ee}}$ is the energy mobility in normal direction
to the free surfaces defined to ensure $\mathcal{C}=0$,
$L^{\parallel}_{\text{ee}}$ is the energy mobility in the tangential direction
to the free surfaces formulated to ensure $L^{\mathrm{sf}}_{\text{ee}}=0$ in
the energy conservation law, and $L^{\text{gb}}_{\text{ee}}$ represents the
energy mobility in the grain boundary. Similarly, $k_{\perp}$ and
$k_{\parallel}$ represent the thermal conductivities at the normal and
tangential directions to the free surfaces respectively while $k_{\text{gb}}$
represent the thermal conductivity in the grain boundary. ${k}_{\mathrm{ss}}$
and ${k}_{\mathrm{at}}$ are respectively the effective thermal conductivities
in the solid phase and atmosphere region and $k_{\text{gb}}$ is the effective
thermal conductivity in the grain boundary. Surface and grain boundary normal
vectors are calculated from the gradient of corresponding OPs, e.g.,
$\mathbf{n}_{\text{sf}}\equiv\nabla\rho/|\nabla\rho|$. $\mathbf{I}$ is the
identity tensor and $\otimes$ represents the dyadic product. $g(\rho)=2\rho-1$
is an odd function that satisfies $g(\rho=\rho_{\mathrm{ss}})=1$ and
$g(\rho=\rho_{\mathrm{at}})=-1$.

Noting that $L^{\perp}_{\text{ee}}=k_{\perp}T^{2}$ and therefore substituting
Eq. (44) into (41), we obtain

$\mathcal{B}=2\chi M_{2}-\frac{\beta
l}{2T^{2}}\left(\frac{1}{2{k}_{\mathrm{ss}}}-\frac{1}{2{k}_{\mathrm{at}}}\right),$
(48)

with

$\chi=l\int^{\infty}_{-\infty}\rho^{{}^{\prime}}(x)\eta^{{}^{\prime}}(x)\mathop{}\\!\mathrm{d}x=2(\rho_{\mathrm{ss}}-\rho_{\mathrm{at}})/3,$
(49)
$\beta=\frac{e_{\text{ht}}}{l}\int^{\infty}_{-\infty}[h(\rho)g(\rho)-1]\mathop{}\\!\mathrm{d}x=-e_{\text{ht}}(\rho_{\mathrm{ss}}-\rho_{\mathrm{at}})^{2},$
(50)

where $e_{\text{ht}}=e_{\mathrm{ss}}-e_{\mathrm{at}}$. The functions defined
in Eq. (38) and (37) is adopted to calculate integrals in Eqs. (49) and (50).

Therefore, to obtain $\mathcal{B}=0$, we take

$M_{2}=\frac{\beta l}{4\chi
T^{2}}\left(\frac{1}{2{k}_{\mathrm{ss}}}-\frac{1}{2{k}_{\mathrm{at}}}\right).$
(51)

Following [17, 8, 53] whereby the different mass diffusion routes in sintering
process i.e bulk/volume diffusion, surface diffusion along the free surfaces
and grain boundary diffusion are taking into account, we propose an
anisotropic $\mathbf{L}_{\rho\rho}$ and relate it to the anisotropic
diffusivity as

$\begin{split}\mathbf{L}_{\rho\rho}&=\left[D_{\mathrm{v}}\mathbf{I}+D_{\mathrm{sf}}\mathbf{T}_{\text{sf}}+D_{\text{gb}}\mathbf{T}_{\text{gb}}\right]/s_{\mathrm{v}}\\\
&=L^{\mathrm{v}}_{\rho\rho}\mathbf{I}+L^{\mathrm{sf}}_{\rho\rho}\mathbf{T}_{\text{sf}}+L^{\text{gb}}_{\rho\rho}\mathbf{T}_{\text{gb}},\end{split}$
(52)

with

$D_{\mathrm{v}}=\left[\frac{1+g(\rho)}{2{D}_{\mathrm{ss}}}+\frac{1-g(\rho)}{2{D}_{\mathrm{at}}}\right]^{-1},$
(53) $D_{\mathrm{sf}}=16\rho^{2}(1-\rho)^{2}D_{\text{sf}},$ (54)
$D_{\text{gb}}=16\sum_{i\neq j}\eta_{i}^{2}\eta_{j}^{2}D_{\text{gb}},$ (55)

and the volumetric entropy as

$s_{\mathrm{v}}=\frac{1}{T}\frac{\partial\mu}{\partial\rho},$ (56)

where the linear approximation is sometimes taken as
$s_{\mathrm{v}}\approx\mathcal{R}/V_{\mathrm{m}}$ with the ideal gas constant
$\mathcal{R}$ and molar volume $V_{\mathrm{m}}$ [17, 54, 8]. In Eq. (52),
$L^{\mathrm{v}}_{\rho\rho}$ is the mass mobility in the normal direction to
the free surfaces associated with the bulk/volume diffusion in solid phase and
atmosphere region, $L^{\mathrm{sf}}_{\rho\rho}$ is the mass mobility in the
tangential direction to the free surfaces associated with mass transport via
surface diffusion. Consideration of $L^{\mathrm{sf}}_{\rho\rho}$ ensures that
the model replicates the sharp-interface mass conservation law Eq. (22) where
surface diffusion is considered. $L^{\text{gb}}_{\rho\rho}$ represent the mass
mobility in the grain boundary. Similarly, $D_{\mathrm{v}}$ represents the
volume diffusivity, which is interpolated by the effective diffusivities in
the solid phase (${D}_{\mathrm{ss}}$) and atmosphere region
(${D}_{\mathrm{at}}$). $D_{\text{sf}}$ and $D_{\text{gb}}$ are the effective
diffusivities in the free surfaces and grain boundary, respectively.

We propose $M_{1}$ to have a similar formulation as $M_{2}$ in Eq. (51);

$M_{1}=-\frac{3ls_{\mathrm{v}}}{16}\left(\frac{A_{\mathrm{ss}}}{2{D}_{\mathrm{ss}}}-\frac{A_{\mathrm{at}}}{2{D}_{\mathrm{at}}}\right),$
(57)

with

$\displaystyle A_{\mathrm{ss}}$
$\displaystyle=\rho_{\mathrm{ss}}+\rho_{\mathrm{at}},$ (58) $\displaystyle
A_{\mathrm{at}}$ $\displaystyle=3A_{\mathrm{ss}}-2.$ (59)

Substituting Eqs. $\eqref{an15}$ and $\eqref{an16}$ into $\eqref{a27}$, we
then obtain

$\mathcal{A}=\frac{\psi
L_{\eta,\text{sf}}^{-1}}{l}-\frac{l(\rho_{\mathrm{ss}}-\rho_{\mathrm{at}})^{2}}{4}\left[\frac{s_{\mathrm{v}}A_{\mathrm{ss}}}{{D}_{\mathrm{ss}}}+\frac{\zeta}{T^{2}}\left(\frac{1}{2{k}_{\mathrm{ss}}}+\frac{1}{2{k}_{\mathrm{at}}}\right)\right],$
(60)

with

$\zeta=\frac{e^{2}_{\text{ht}}}{l}\int^{\infty}_{-\infty}[1-h^{2}(\rho)]\mathop{}\\!\mathrm{d}x=e^{2}_{\text{ht}},$
(61)

and

$\psi=l\int^{\infty}_{-\infty}(\eta^{{}^{\prime}}(x))^{2}\mathop{}\\!\mathrm{d}x=2/3,$
(62)

also taking into account:

$\begin{split}\int^{\infty}_{-\infty}\left[(\rho_{\text{eq}}(x))^{2}+g(\rho_{\text{eq}})(\rho_{\text{eq}}(x))^{2}-1\right]\mathop{}\\!\mathrm{d}x\\\
=-\frac{3l}{4}A_{\mathrm{ss}}(\rho_{\mathrm{ss}}-\rho_{\mathrm{at}})^{2},\end{split}$
(63)

$\int^{\infty}_{-\infty}\left[(\rho_{\text{eq}}(x))^{2}[1-g(\rho_{\text{eq}})\right]\mathop{}\\!\mathrm{d}x=\frac{l}{4}A_{\mathrm{at}}(\rho_{\mathrm{ss}}-\rho_{\mathrm{at}})^{2}.$
(64)

The functions defined in Eq. (38) and (37) are again adopted to calculate
integrals in Eq. (62) and Eqs. (61) - (64).

Therefore, in order to ensure $\mathcal{A}=0$, we take

$L_{\eta,\text{sf}}^{-1}=\frac{l^{2}(\rho_{\mathrm{ss}}-\rho_{\mathrm{at}})^{2}}{4\psi}\left[\frac{s_{\mathrm{v}}A_{\mathrm{ss}}}{{D}_{\mathrm{ss}}}+\frac{\zeta}{T^{2}}\left(\frac{1}{2{k}_{\mathrm{ss}}}+\frac{1}{2{k}_{\mathrm{at}}}\right)\right].$
(65)

The mobility $L_{\eta,\text{gb}}^{-1}$ of $\\{\eta_{i}\\}$ can be obtained
from the physical grain boundary mobility $G_{\text{gb}}^{\text{eff}}$ and
grain boundary energy $\gamma_{\text{gb}}$ as [55, 8]

$L_{\eta,\text{gb}}^{-1}=\frac{\kappa_{\eta}}{G_{\text{gb}}^{\text{eff}}\gamma_{\text{gb}}}.$
(66)

Recalling Eqs. (17c) and (18), it is worth mentioning that this mobility is
defined under the driving force represented by entropy, which should be
distinguished from the original formulation in Ref. [55] as here
$\kappa_{\eta}$ adopts the dimension of the entropy per length. Accordingly,
for $L_{\eta}^{-1}$ as regards the full sintering description, we then take

$L_{\eta}^{-1}=16\rho^{2}(1-\rho)^{2}L_{\eta,\text{sf}}^{-1}+L_{\eta,\text{gb}}^{-1}.$
(67)

Recalling the anisotropic definitions of $\mathbf{L}_{ee}$ and
$\mathbf{L}_{\rho\rho}$ in Eqs. (43) and (52), calculation of
$\hat{L}_{\eta}^{-1}$ in Eq. (18) can be further simplified as

$\hat{L}_{\eta}^{-1}=L_{\eta}^{-1}-l^{2}|\nabla\rho|^{2}[M^{2}_{1}L_{\rho\rho}^{\mathrm{v}}+M_{2}^{2}L_{ee}^{\perp}],$
(68)

as $L_{ee}^{\perp}$ and $L_{\rho\rho}^{\mathrm{v}}$ are respectively one of
the eigen-values of $\mathbf{L}_{ee}$ and $\mathbf{L}_{\rho\rho}$,
corresponding to the eigen-direction of $\mathbf{n}_{\mathsf{sf}}$
($\mathbf{n}_{\mathsf{sf}}\equiv\nabla\rho/|\nabla\rho|$).

It is worth noting that $M_{1}$ and $M_{2}$ are derived based on the constant
postulate, i.e., $M_{1}$ and $M_{2}$ are spatio-temporal independent constants
for a sintering system with known mass diffusivities and thermal
conductivities of substance and atmosphere as well as given diffuse interface
width, since the spatio-temporal dependency of all OP-related terms (Eqs.
(49)-(50), (62), and (61)-(64)) vanish after integral. More importantly, the
quantitative phase-field model degenerates to the conventional one when the
system has no differences in mass diffusivity and thermal conductivity between
solid and atmosphere. In that sense, when
${D}_{\mathrm{ss}}={D}_{\mathrm{at}}$, $M_{1}=0$ and also when
${k}_{\mathrm{ss}}={k}_{\mathrm{at}}$, $M_{2}=0$, demonstrating that the
antitrapping terms in Eqs (17a) and (17b) and cross-coupling term in Eq. (17c)
reduce to zero. In addition, we note that variational quantitative phase-field
models such as the one presented in this work do not generally demonstrate
high numerical accuracy [33, 44, 45]. Correct mapping of the variational model
onto the associated sharp-interface equations only guarantees its quantitative
validity and not its numerical efficiency needed for realistic utilization
[44]. Therefore, a nonvariational form of the model might be best suited for
practicability. The nonvariational form can be simply developed via
modification of model parameters and functions while ensuring that the thin-
interface asymptotic remains consistent.

Table 1:  Set of dimensionless quantities and parameters employed for the simulations in this work. $L_{x}$ |  $L_{y}$ |  $g_{0}$ |  $\underline{C}$ |  $\underline{D}$ |  $e_{\text{ht}}$ |  $L_{\eta,\text{gb}}$ |  $s_{\mathrm{v}}$
---|---|---|---|---|---|---|---
60 |  50 |  0.01 |  1 |  0.062 |  1 |  1 |  1
Figure 2:  (a) Schematic of the simulation setup in 2D simulation of an
elliptical inclusion (b) Comparison of chemical potential jump $\delta\mu$
across the free surface with respect to interface width; blue triangle symbols
represent model with $M_{1}=0$ while red diamond symbols represent model with
$M_{1}\neq 0$. Plots of $\mu$ and $\rho$ across the free surface as a function
of $x$ with $l=1$, $D_{\mathrm{at}}/D_{\mathrm{ss}}=2$ for (c) $M_{1}\neq 0$
and (d) $M_{1}=0$. $\delta\mu$ is obtained using an extrapolation of $\mu$ at
the center of the free surface $\rho=0.5$. $\Delta\mu$ is the chemical
potential difference between the bulk values.

##  V Results and discussion

###  V.1 Model verification for an elliptical inclusion

In order to examine the capability of the model in ensuring $\delta\mu=\delta
T=0$ at the free surface, we perform diffusion-driven reshaping simulations of
an elliptical inclusion with major axis $A$ and minor axis $B$ morphing into a
circle. We set up a simulation domain with the lengths $L_{x}$ and $L_{y}$ in
$x$ and $y$ direction, respectively. The domain is further subjected to an
initial temperature gradient $\nabla T=g_{0}$ along the x-axis. $\rho$ is
taken to vary smoothly from one in the inclusion ($\rho_{\mathrm{ss}}=1$) to
zero outside ($\rho_{\mathrm{at}}=0$) with $l=1$. A full schematic of the
simulation setup is given in Fig. 2a. The normalized values of the employed
model parameters are given in Table 1.

First, we consider a case of asymmetric mass transport where
$D_{\mathrm{at}}/D_{\mathrm{ss}}=2$. We set
$k_{\mathrm{at}}/k_{\mathrm{ss}}=1$, hence only employing the mass
antitrapping current term associated with mass diffusion while the thermal
antitrapping current is tentatively dropped. Profile of chemical potential
$\mu(x)$ across the moving free surface is presented for the cases $M_{1}=0$
and $M_{1}\neq 0$ in Figs. 2c and 2d, respectively. An extrapolation of
$\mu(x)$ gives the chemical potential jump ($\delta\mu$) at the center of the
free surface $\rho=0.5$. Typically, $\delta\mu\neq 0$ implies an exchange of
mass between the solid and atmosphere, which can be likened to the trans-
interface diffusion phenomenon. However, no mass exchange is expected between
the solid and atmosphere regions during sintering. Therefore, $\delta\mu=0$
should be held in phase-field simulations in order to achieve realistic mass
diffusion. It is obvious from Figs. 2c and 2d that the case with $M_{1}=0$
shows a significantly larger $\delta\mu$ compared to the one with $M_{1}\neq
0$, in which the relatively small $\delta\mu$ is attributed to possible
numerical errors. The results demonstrate that the mass antitrapping current
parameterized by $M_{1}$ is necessary in order to eliminate the artificial
diffusion flux across the interface during sintering for cases of asymmetric
mass transport. Figs. 2c and 2d also show the gap in space $\delta x$ between
the center of the free surface and the point where the extrapolations of $\mu$
meet. Note that $\delta x=0$ when $\delta\mu=0$, indicating the coherence
between the numerically predicted interface by $\rho\approx 0.5$ and the
theoretical sharp interface where $\delta\mu=0$. Similar to $\delta\mu$, the
numerical results demonstrate a significantly larger $\delta x$ for the case
with $M_{1}=0$ compared to the one with $M_{1}\neq 0$, implying an apparent
deviation in the position between the predicted interface and theoretical
sharp interface.

We also note the existence of another chemical potential drop $\Delta\mu$
across the free surface, characterizing the differences between the bulk
values, as depicted in Figs. 2c and 2d. This $\Delta\mu$, which is identical
for both cases at a time point, were numerically examined to be the outcome of
the deviated bulk values of $\rho$, i.e., $\rho_{\mathrm{ss}}$ and
$\rho_{\mathrm{at}}$ that are slightly deviated from ideal (equilibrium) one
and zero respectively, as listed in Table. S1. Such chemical potential drop
generally does not appear in the conventional sharp-interface interpretation
of the sintering [1, 2]. Meanwhile, the deviated bulk values of the conserved
mass OP have been depicted in previous works [56, 57, 58, 59, 60] with
theoretical and numerical analyses given in Refs. [59] and [60], which are
further discussed in the Appendix.

Furthermore, Fig. 2b shows a comparison of $\delta\mu$ vs. diffuse interface
width $l$ between the cases with/without $M_{1}$ parameterized. It can be
observed that both cases present the convergence $\delta\mu\to 0$ as $l\to 0$,
replicating the sharp-interface condition when $l$ tends to infinitesimal.
However, as $l$ increases, $\delta\mu$ present a relatively rapid growth in
the case with $M_{1}=0$ compared to the one with $M_{1}\neq 0$, demonstrating
that the employment of mass antitrapping current parameterized by $M_{1}$ can
significantly reduce the artificial interface effect (here the growing
$\delta\mu$) along with increasing diffuse interface width. In this sense,
mass antitrapping current allows reasonable quantitative simulations,
especially at larger interface widths. Furthermore, we note that the
convergence of both models might be well investigated considering a steady
state free surface velocity. This we hope to report in our upcoming work.

Additionally, we investigate a case of asymmetric heat transport with
$k_{\mathrm{at}}/k_{\mathrm{ss}}=0.05$. Similar to previous simulation, we
examine the thermal antitrapping term associated with heat transport. The mass
antitrapping term is tentatively dropped by setting
$D_{\mathrm{at}}/D_{\mathrm{ss}}=1$. Simulations are performed for existing
model (i.e., $M_{2}=0$) and quantitative model with $M_{2}\neq 0$. Further
details of results are given in supplementary material (Fig. S2). $\delta T=0$
realized at the sharp-interface is expected to be obtained during phase-field
simulations in order to guarantee quantitative simulations. For model with
$M_{2}=0$, however, emerging $\delta T\neq 0$ demonstrates the importance of
the thermal antitrapping current. Here, it is important to note that measured
$\delta T$ has a relatively low magnitude compared to the bulk temperature at
the free surface. Importance of thermal antitrapping term $M_{2}$ in
eliminating temperature jump for asymmetric heat transport has also been
demonstrated in Ref. [41] where non-diagonal phase-field model was also used.

Figure 3:  Snapshots of grains coalescence of two spherical grains with
distinct sizes. $D_{\mathrm{at}}/D_{\mathrm{ss}}=2$, $l=2$ and
$k_{\mathrm{at}}/k_{\mathrm{ss}}=0.05$ are set. Comparison is made between
models with ATs and without ATs. Temperature isolines are also indicated.
$t^{*}=10^{3}$ unit. Figure 4:  Surface profile ($\rho=0.5$) colored by the
curvature and the mass diffusion fluxes, i.e., the fluxes before
($\mathbf{J}_{\mathrm{\rho}}$) and after
($\mathbf{J}_{\mathrm{\rho}}^{\prime}$) correction with the antitrapping
contribution ($\mathbf{J}_{\mathrm{\rho,AT}}$) , at two distinct sites. The
length of visualized arrows have been scaled according to the magnitude of the
fluxes uniformly.

###  V.2 Comparison between the quantitative and the existing models

In this section, we perform simulations for grain coalescence of two spherical
grains with distinct sizes. Comparisons of microstructure and temperature
distribution are made between quantitative model where antitrapping terms are
taken into account and existing model where these terms are not considered.
The two models are referred to as model with ATs and model without ATs in the
following discussions.

We set up a simulation domain with the lengths $L_{x}$ and $L_{y}$ in $x$ and
$y$ direction, respectively. Similar to previous setup, the domain is
subjected to $\nabla T=g_{0}$ along the x-axis. Simulations are performed for
asymmetric mass and heat transport where $D_{\mathrm{at}}/D_{\mathrm{ss}}=2$,
$k_{\mathrm{at}}/k_{\mathrm{ss}}=0.05$ with $l=2$. A full schematic of the
simulation setup is supplemented in Fig. S4a.

Transient microstructures and temperature profiles for both models are
compared and presented in Fig. 3. First, we observe that mass transport was
faster for model without ATs compared to model with ATs. At $t/t^{*}=0.633$
and $t/t^{*}=0.815$ in Figs. 3c and 3d respectively, a more coalesced grain is
obtained for model without ATs compared to model with ATs. The difference in
progress of coalescence can be further explained by the visualization of mass
diffusion fluxes at the free surface $\rho=0.5$ as presented in Fig. 4. The
free surface profile is colored by the local curvature calculated as
$-\nabla\cdot\mathbf{n}_{\mathrm{sf}}$. Furthermore, mass diffusion fluxes are
indicated at two distinct points; a concave point and a convex point.
Typically, mass flux at any point on the free surface is expected to be
correctly captured along the tangential direction to the free surface at that
point. It can be clearly observed that $\mathbf{J}_{\rho}$ which is the mass
diffusion flux without the mass antitrapping current deviates in direction
from the tangential direction (dashdot lines) to the free surface at both
concave and convex points. The mass antitrapping flux $\mathbf{J}_{\rho,AT}$
introduced in quantitative model can be seen flowing through the free surface
in the normal direction from the solid grain region to the atmosphere. The
combined mass flux
$\mathbf{J}^{{}^{\prime}}_{\rho}=\mathbf{J}_{\rho}+\mathbf{J}_{\rho,AT}$ shows
a corrected mass flux flowing along the tangential direction to the free
surface. Therefore, the deviation of $\mathbf{J}_{\rho}$ from its appropriate
direction is due to the existence of chemical potential jump at the free
surface. $\mathbf{J}_{\rho,AT}$ serves to eliminate this chemical potential
jump which consequently corrects this deviation. Accordingly, this
demonstrates the faster mass transport observed for model without ATs.
Chemical potential jump at the free surface tend to act as an extra driving
force for grain coalescence leading to faster mass diffusion. The elimination
of this jump via the antitrapping current leads to a slower mass transport for
model with ATs.

Figure 5:  Time evolution of (a) sintering neck $\lambda$ and (b) dihedral
angle $\Phi$ during the non-isothermal sintering process, as shown in Fig. 3.
The time points reaching maximum values are indicated by colored vertical
lines. The equilibrium dihedral angle $\Phi_{\mathrm{eq}}$, calculated from
the surface and grain boundary energies, is also indicated by black dotted
line in (b). The total simulation time $t^{*}=2084$ unit.

Furthermore, results in Fig. 3 also show the comparison of the temperature
profiles obtained for both models. The distribution of the temperature
isolines shows faster heat transport for model without ATs compared to model
with ATs. An example is given for isoline $T=0.71$. Even though it initially
tends to migrate towards high-$T$ side, at $t/t^{*}=0.023$, this tendency
breaks for model without ATs, where the isoline starts to move towards low-$T$
side, but continues for model with ATs. Result at $t/t^{*}=0.815$ in Fig. 3d
indicates a colder grain for model with ATs compared to model without ATs.
Similar to the mass transport fluxes explanation for both models, the
temperature jump at the free surface can be seen as an extra driving force for
heat transport in the model without ATs. This jump is eliminated for model
with ATs via the thermal antitrapping term thereby obtaining a slower heat
transport. The presented thermal-microstructure evolution once more
demonstrates the importance of the antitrapping currents for mass and heat
diffusion.

We further examine the in-process sintering neck $\lambda$ and dihedral angle
$\Phi$ of the simulation as presented in Fig. 3. $\lambda$ and $\Phi$ are
calculated by

$\displaystyle\lambda$ $\displaystyle=\int_{\Omega}\frac{16\sum_{i\neq
j}\eta_{i}^{2}\eta_{j}^{2}}{l}\mathop{}\\!\mathrm{d}\Omega,$ (69)
$\displaystyle\Phi$ $\displaystyle=\arctan\left(\frac{\partial
C_{\eta_{i}}}{\partial x}\right)_{\mathrm{neck}}-\arctan\left(\frac{\partial
C_{\eta_{j}}}{\partial x}\right)_{\mathrm{neck}},$

where $C_{\eta_{i}}$ and $C_{\eta_{j}}$ are the fitted semi-circular
tendencies by coordinates of contour $\eta_{i}=0.5$ and $\eta_{j}=0.5$,
respectively. $\frac{\partial C_{\eta_{i}}}{\partial x}$ and $\frac{\partial
C_{\eta_{j}}}{\partial x}$ then provide the slopes of $C_{\eta_{i}}$ and
$C_{\eta_{j}}$. In this sense, $\Phi$ is calculated using the difference
between these two angles of slope at the neck point, as shown in inset of Fig.
5b, adapted from Ref. [61]. Meanwhile, the equilibrium dihedral angle
$\Phi_{\mathrm{eq}}$ can be also evaluated by the surface
($\gamma_{\mathrm{sf}}$) and grain boundary ($\gamma_{\mathrm{gb}}$) energies,
i.e.,

$\Phi_{\mathrm{eq}}=2\arctan\frac{\gamma_{\mathrm{gb}}}{2\gamma_{\mathrm{sf}}}.$
(70)

It is worth noting that $\Phi$ approaches $\Phi_{\mathrm{eq}}$ when two
particles with identical size are sintered isothermally, as $\lambda$ reaches
the maximum and stays constant, i.e., the system reaches equilibrium [8]. With
varying interface width $l$, $\Phi$ deviates from theoretically-determined
$\Phi_{\mathrm{eq}}$ (Eq. (70)) as shown in Fig. S3. This deviation is reduced
in a similar fashion for both models with/without ATs as $l$ decreases. This
implies no modification to thermodynamic equilibrium condition (characterized
by $\Phi_{\mathrm{eq}}$) by applying the kinetic antitrapping terms. For two
non-identical grains, the time evolution of $\lambda$ and $\Phi$ are presented
in Fig. 5. Comparison is made for model with ATs and model without ATs. It can
be observed that for both models, $\Phi$ approaches $\Phi_{\mathrm{eq}}$ at
the points where $\lambda$ attains maximum values. However, the progress of
$\Phi$ towards $\Phi_{\mathrm{eq}}$ is faster for model without ATs compared
to model with ATs. This implies that while thermodynamic conditions are
attained for both models, the antitrapping terms tend to modify the progress
of neck growth and grain coalescence by removing the extra flux perpendicular
to the free surface, as evidently shown in Fig. 3 and Fig. 4.

Figure 6:  (a) Comparison of ${J}_{e,x}$ vs $l/R$ for different thermal
conductivity interpolations: direct interpolation (red diamonds), inverse
interpolation (blue circles), special function interpolation (green hexagons)
and anisotropic interpolation (orange triangles) for
${k}_{\mathrm{at}}/{k}_{\mathrm{ss}}=10$ where corresponding color lines are
fitted simulation data lines and black dashdot lines represent the plot of a
case where no artificial interface effects exist. (b) Temperature isolines
across the free surface for different thermal conductivity interpolations at
$l/R=0.05$ (c) Time evolution of ${J}_{\rho,x}$ using different diffusivity
interpolations: isotropic (dashdot lines), anisotropic (solid lines) for
varied values of $D_{\text{sf}}/D_{\mathrm{ss}}$. Local zooms around the neck
of the grains using ($\mathrm{d_{1}}$) isotropic interpolation of diffusivity
($\mathrm{d_{2}}$) anisotropic interpolation of diffusivity. $t^{*}=14\times
10^{3}$ unit.

###  V.3 Importance of anisotropic interpolations of the mobility tensor

Here, we particularly demonstrate the importance of the anisotropic
interpolations of the kinetic mobilities. First, we investigate a steady-state
heat transfer case. The numerical validation test proposed by Nicoli et al.
[52] is used and extended. We consider a square simulation domain defined as
$[0,1]$ and $[0,1]$ in $x$ and $y$ direction, respectively and subjected to
$\nabla T=-2$ along the x-axis. The domain consists of a stationary disk-
shaped solid grain with radius $R$ surrounded by an atmosphere region. A
schematic of the simulation setup is supplemented in Fig. S5a. For
${k}_{\mathrm{at}}/{k}_{\mathrm{ss}}=10$, four cases of thermal conductivity
interpolations are examined. We consider the form of interpolation utilized in
current phase-field models of non-isothermal sintering [28, 8]. This
interpolation form is similar to the formulation given in Eq. (45) and is
thereafter referred to as the direct interpolation. Also, we consider another
form of interpolation given in Ref. [62] to ensure heat flux conservation
across the interface. This form of interpolation, thereafter referred to as
the inverse interpolation, has its formulation as in Eq. (44). Moreover, a
form of interpolation proposed by Almgren [30] was used in Refs. [41] and [43]
for their non-diagonal phase-field models. The interpolation thereafter called
the special function (SF) interpolation is also examined and can be expressed
as

$\frac{1}{k(\rho)}=\bigg{(}\frac{1}{2{k}_{\mathrm{ss}}}+\frac{1}{2{k}_{\mathrm{at}}}\bigg{)}+p_{\text{sf}}(\rho)\bigg{(}\frac{1}{2{k}_{\mathrm{ss}}}-\frac{1}{2{k}_{\mathrm{at}}}\bigg{)},$
(71)

with

$p_{\text{sf}}(\rho)=(2\rho-1)[1+4a\rho(1-\rho)],$ (72)

where $a\approx 0.90$ for ${k}_{\mathrm{at}}/{k}_{\mathrm{ss}}=$ 10\. Here, it
is important to note that the formulation in Eq. (71) is adopted from Ref.
[43] because $\phi$ in Ref. [43] varies from 0 to 1 similar to $\rho$. Lastly,
we consider the anisotropic form of thermal conductivity proposed in this work
as expressed in Eq. (43). Artificial interface effects are quantified by
obtaining the average heat flux, ${J}_{e,x}$ in the domain at $x=1$. The plot
of ${J}_{e,x}$ against normalized interface widths $l/R$ are presented in Fig.
6a for different interpolation forms. The black line in Fig. 6a indicates a
reference case where no artificial interface effect exist, i.e
${J}_{e,x}\mid_{l=0}\>=\>{J}_{e,x}\mid_{l>0}$. Fig. 6b shows the temperature
isolines across the free surface for the different interpolations forms. As
shown in Fig. 6a, the direct and inverse interpolations show significant
deviations from the reference case implying the deficit of these interpolation
forms in eliminating interface effects. The SF interpolation also shows
considerable deviation from the reference case. This deviation, which might be
attributed to the nonmonotonic form of $p_{\text{sf}}(\rho)$ [44] reinforces
the limitation of the SF interpolation. On the other hand, the results by
using the anisotropic interpolation shows very convincing agreement. The
outstanding performance of the anisotropic form of interpolation necessitates
its consideration for subsequent non-diagonal phase-field modeling. Note that
for common sintering scenarios, (result is supplemented in Fig. S5b for
${k}_{\mathrm{at}}/{k}_{\mathrm{ss}}=0.2$), the SF interpolation might be
utilized for quantitative simulations. The anisotropic form of interpolation,
however, finds great importance in other processes where
${k}_{\mathrm{at}}/{k}_{\mathrm{ss}}$ is higher such as the case studies in
Ref. [63].

Also, we consider mass transport during grain coalescence of two identical
spheres using two mass diffusivity interpolations. A full schematic of the
simulation setup is supplemented in Fig. S4b. We make comparisons between the
anisotropic interpolation presented in this work, Eq. (52) and an isotropic
interpolation expressed as [15]:

$D=p_{\mathrm{ss}}(\rho){D}_{\mathrm{ss}}+p_{\mathrm{at}}(\rho){D}_{\mathrm{at}}+p_{\text{sf}}(\rho)D_{\text{sf}}+p_{\text{gb}}(\eta_{i})D_{\text{gb}},$
(73)

where $p_{\mathrm{ss}}(\rho)$ and $p_{\mathrm{at}}(\rho)$ are interpolation
functions valued as one only in the solid phase and atmosphere region
respectively. Average mass flux ${J}_{\rho,x}$ is obtained across a grain with
the plots of ${J}_{\rho,x}$ against normalized time $t$ shown in Fig. 6c.

The thin-interface limit analysis showed that eliminating chemical potential
jump across the free surface does not require a specific mass diffusivity
interpolation. Correspondingly, it has been derived in Refs. [64, 65] that the
Cahn-Hilliard equation recovers the sharp-interface limit equation of motion
for surface diffusion regardless of mass diffusivity form. Therefore, the
${J}_{\rho,x}$ vs $t$ curve is expected to be the same for both interpolations
of mass diffusivity, since theoretically no artificial interface effect is
related to the diffusivity interpolation. However, as shown in Fig. 6c, there
exist surprisingly ${J}_{\rho,x}$ numerical deviations. This can be explained
by close comparison of flux details at the free surface region. Fig.
6$\mathrm{d_{1}}$ and Fig. 6$\mathrm{d_{2}}$ demonstrate the calculated flux
by using the isotropic and anisotropic diffusivity form, respectively. We
observe that the anisotropic form of diffusivity delivers more reasonable
description of the directions of the fluxes. Around the free surface in Fig.
6$\mathrm{d_{1}}$ where we used isotropic diffusivity, there exist non-
tangential fluxes at the free surface where only tangential fluxes are
expected to contribute to surface diffusion. On the other hand, in Fig.
6$\mathrm{d_{2}}$ where the anisotropic diffusivity form is used, only fluxes
that are tangential to the free surface region exist to describe surface
diffusion. Accordingly, it is imperative that, while asymptotic analysis
confers no restriction on the diffusivity form as regards effecting
quantitative simulations in mass diffusion, the anisotropic diffusivity form
however makes it possible that the directions of fluxes are effectively
described analogous to the sharp-interface description.

##  VI Conclusions

In this work, we have developed a variational quantitative phase-field model
for non-isothermal sintering processes following the non-diagonal phase-field
approach introduced in Refs. [37, 39]. The model was formulated to eliminate
artificial interface effects due to the diffuse-interface description of the
free surfaces. Moreover, model formulations are derived in a variational
manner guaranteeing thermodynamic consistency. The proposed model differs from
conventional non-isothermal sintering models owing to that fact that cross-
coupling terms between conserved kinetics (mass and heat transfer) and the
non-conserved kinetics (grain growth) are taken into the account. These terms
parameterized by functions $M_{1}$ and $M_{2}$ can be likened to antitrapping
currents in quantitative phase-field modeling. The above-mentioned terms are
particularly essential for correct projection of the model to its sharp-
interface descriptions. Also, we derive formulations of $M_{1}$ and $M_{2}$ in
terms of the model parameters using an asymptotic analysis procedure presented
in Ref [38]. In addition we showed that anisotropic interpolations of kinetic
mobilities are also important to ascertain the elimination of artificial
interface effects at the free surface.

Numerical tests were done to highlight the importance of these cross-coupling
terms. The results presented showed the emergence of chemical potential jump
($\delta\mu$) and temperature jump ($\delta T$) at the free surface when
$M_{1}=0$ and $M_{2}=0$. $\delta\mu\neq 0$ and $\delta T\neq 0$ negates the
sharp-interface sintering description. However, employing $M_{1}\neq 0$ and
$M_{2}\neq 0$ as described in quantitative model eliminates these jumps. The
convergence behavior of $\delta\mu$ in respect to interface width ($l$) was
presented for model with $M_{1}=0$ and model with $M_{1}\neq 0$. For both
models, $\delta\mu\to 0$ as $l\to 0$ demonstrating their efficacy at
relatively smaller $l$. The major usefulness of quantitative model is seen as
$l\gg 0$ where $\delta\mu$ is significantly large for model with $M_{1}=0$
compared to model where $M_{1}\neq 0$. Additionally, the difference in
transient microstructure and temperature profiles were examined for model with
antitrapping currents and model without antitrapping currents. It was seen
that the antitrapping currents helps to eliminate extra driving forces brought
about by $\delta\mu\neq 0$ and $\delta T\neq 0$ at the free surface. Moreover,
it was demonstrated that the antitrapping currents only modify the sintering
kinetics and have no impact on the thermodynamic conditions.

Furthermore, we demonstrated numerically how the anisotropic interpolation of
kinetic mobilities delivers effective description of diffusion fluxes
comparable to sharp-interface description. Therefore, the proposed model can
serve as a great tool in studying quantitative simulations of non-isothermal
sintering and other related solid-state processes. A major outlook of this
work is to further investigate the convergence of interface velocity with
respect to interface width obtained using proposed model.

###### Acknowledgements.

Authors acknowledge the financial support of German Science Foundation (DFG)
in the Priority Program 2256 (SPP 2256, project number 441153493) and
Collaborative Research center Transregio 270 (CRC-TRR 270, project number
405553726, sub-projects A06). The authors also greatly appreciate their access
to the Lichtenberg High-Performance Computer and the technique supports from
the HHLR, Technische Universität Darmstadt.

## DATA AVAILABILITY

The authors declare that the data supporting the findings of this study are
available within the paper. Source codes of MOOSE-based application NIsoS and
related utilities are provided in the online repository
bitbucket.org/mfm_tuda/nisos.git.

*

##  Appendix A Deviation of conserved order parameter

The analyses of Cahn-Hilliard dynamics in Refs. [59] and [60] have showed that
usage of finite interface width combined with comparable curvature radius
induces deviation of the conserved order parameter ($\rho$ in this work) in
the bulk regions. It has been demonstrated that the equilibrium bulk values of
$\rho$ are contingent on the interface having negligible volume compared to
the bulk region, so that only the local free energy finds minimization.
Although this condition is viable for planar interfaces, it is not maintained
for curved interfaces with concentrated energy. In this sense, the total free
energy can be reduced by shrinking the area enclosed by the interface, which
subsequently shifts the bulk values of $\rho$ from the equilibrium ones due to
the finite volume precept [59]. Here, we define the deviated quantities of
$\rho$ from its equilibrium values (in this work
$\rho_{\mathrm{ss}}^{\mathrm{eq}}=1$ and $\rho_{\mathrm{at}}^{\mathrm{eq}}=0$)
as $\Delta\rho_{\mathrm{ss}}=\rho_{\mathrm{ss}}-1$ and
$\Delta\rho_{\mathrm{at}}=\rho_{\mathrm{at}}$. Both $\rho_{\mathrm{ss}}$ and
$\rho_{\mathrm{at}}$ are read from numerical results in Fig. 2a with $M_{1}=0$
when the particle is in elliptical and circular shapes. The tendencies of
$\Delta\rho_{\mathrm{ss}}$ and $\Delta\rho_{\mathrm{at}}$ vs. $l$ are
respectively shown in Fig. A1a. Similar to results obtained in Ref. [59],
$\Delta\rho_{\mathrm{ss}}$ and $\Delta\rho_{\mathrm{at}}$ increase with
increasing $l$. When in the elliptical shape (implying a non-equilibrium
condition), $\Delta\rho_{\mathrm{at}}>\Delta\rho_{\mathrm{ss}}$ holds for
almost every selected $l$, while
$\Delta\rho_{\mathrm{ss}}\approx\Delta\rho_{\mathrm{at}}$ when in the circular
shape (implying an equilibrium condition). These differences can be attributed
to the curvature dependency of the analytical profile of $\rho$ [60]. It
should be noted that $\Delta\rho_{\mathrm{ss}}$ and $\Delta\rho_{\mathrm{at}}$
exist even for symmetric mobilities with sufficiently large $l$, which is
distinctive to the known interface effects (like trap effects) that are
incited by asymmetric kinetic mobilities. Moreover, the comparison of
analytical values of $\Delta\rho_{\mathrm{ss}}$ and $\Delta\rho_{\mathrm{at}}$
and numerical values should be examined in further studies.

As one of the significant outcomes, deviated bulk values of $\rho$ incite
deviated chemical potential $\mu$ from its equilibrium ones in the bulk
regions, which may result in the unexpected chemical potential drop as an
extra driving force across the free surface. To examine this point, we define
the deviated quantities of $\mu$ in a similar fashion as
$\Delta\rho_{\mathrm{ss}}$ and $\Delta\rho_{\mathrm{at}}$, i.e.,
$\Delta\mu_{\mathrm{ss}}=\mu_{\mathrm{ss}}(\rho_{\mathrm{ss}})-\mu^{\mathrm{eq}}=\mu_{\mathrm{ss}}(\rho_{\mathrm{ss}})$
and
$\Delta\mu_{\mathrm{at}}=\mu_{\mathrm{at}}(\rho_{\mathrm{at}})-\mu^{\mathrm{eq}}_{\mathrm{at}}=\mu_{\mathrm{at}}(\rho_{\mathrm{at}})$,
noting that
$\mu^{\mathrm{eq}}_{\mathrm{ss}}(\rho^{\mathrm{eq}}_{\mathrm{ss}}=1)=\mu^{\mathrm{eq}}_{\mathrm{at}}(\rho^{\mathrm{eq}}_{\mathrm{at}}=0)=0$.
In Fig. A1b, we present a similar tendency of $\Delta\mu_{\mathrm{ss}}$ and
$\Delta\mu_{\mathrm{at}}$ vs. $l$ when the particle is in the elliptical
shape, where both $\Delta\mu_{\mathrm{ss}}$ and $\Delta\mu_{\mathrm{at}}$ grow
along with increasing $l$ and
$\Delta\mu_{\mathrm{ss}}>\Delta\mu_{\mathrm{at}}$ is depicted for every
selected $l$, implying the existing chemical potential drop
$\Delta\mu=\Delta\mu_{\mathrm{at}}-\Delta\mu_{\mathrm{ss}}>0$ across the free
surface at the semi-major axis, as shown in Fig. 2a. Notably, when the
particle is in the circular shape,
$\Delta\mu_{\mathrm{ss}}=\Delta\mu_{\mathrm{at}}$ is formed without the
dependency of $l$, indicating no chemical potential drop across the free
surface, i.e., $\Delta\mu=\Delta\mu_{\mathrm{at}}-\Delta\mu_{\mathrm{ss}}=0$.
This also demonstrates that the existing deviation in $\mu$ incited by
$\Delta\rho_{\mathrm{ss}}$ and $\Delta\rho_{\mathrm{at}}$ does not affect the
supposing equilibrium condition, as the particle stops morphing in the
circular shape.

Additionally, since the antitrapping coefficients $M_{1}$ and $M_{2}$ are
dependent on the bulk values $\rho_{\mathrm{ss}}$ and $\rho_{\mathrm{at}}$, we
examine the variations of $M_{1}$ and $M_{2}$ with $\Delta\rho_{\mathrm{ss}}$
and $\Delta\rho_{\mathrm{at}}$ up to 0.1, as seen in Fig. A2. In Fig. A2a, we
take $\Delta\rho_{\mathrm{at}}=0$ and examine the variations of $M_{1}$ and
$M_{2}$ with $\Delta\rho_{\mathrm{ss}}$. Similarly, we take
$\Delta\rho_{\mathrm{ss}}=0$ and examine the variations of $M_{1}$ and $M_{2}$
with $\Delta\rho_{\mathrm{at}}$ in Fig. A2b. Then, we present the variations
of $M_{1}$ and $M_{2}$ with
$\Delta\rho_{\mathrm{ss}}=\Delta\rho_{\mathrm{at}}$ in Fig. A2c. It
demonstrates that $M_{1}$ presents a linear tendency vs. increasing deviations
of all cases. $M_{2}$, however, decreases along with growing
$\Delta\rho_{\mathrm{ss}}$ but increases with growing
$\Delta\rho_{\mathrm{at}}$. For
$\Delta\rho_{\mathrm{ss}}=\Delta\rho_{\mathrm{at}}$, $M_{2}$ stays constant.
This can be explained via Eq. (51) where $M_{2}$ is proportional to
$(\rho_{\mathrm{ss}}-\rho_{\mathrm{at}})$, which is reduced to one when
$\Delta\rho_{\mathrm{ss}}=\Delta\rho_{\mathrm{at}}$ as
$(\rho_{\mathrm{ss}}-\rho_{\mathrm{at}})=[(\Delta\rho_{\mathrm{ss}}+1)-\Delta\rho_{\mathrm{at}}]=1$.

Figure A1:  The deviated bulk values of (a) $\rho$ and (b) $\mu$ with respect
to the diffuse interface width $l$. The deviated bulk values are read from the
numerical results presented in Fig. 2a. Figure A2:  The antitrapping
coefficients $M_{1}$ and $M_{2}$ with respect to (a) varying
$\Delta\rho_{\mathrm{ss}}$ when $\Delta\rho_{\mathrm{at}}=0$, (b) varying
$\Delta\rho_{\mathrm{at}}$ when $\Delta\rho_{\mathrm{ss}}=0$, and (c) varying
simutaneously $\Delta\rho_{\mathrm{ss}}$ and $\Delta\rho_{\mathrm{at}}$ while
holding $\Delta\rho_{\mathrm{ss}}=\Delta\rho_{\mathrm{at}}$.

## References

  * German [2014] R. German,  _Sintering: from empirical observations to scientific principles_ (Butterworth-Heinemann, 2014).
  * Kang [2004] S.-J. Kang,  _Sintering: Densification, Grain Growth and Microstructure._ (Elsevier, 2004).
  * W. and J. [2010] N. W. and P. J., 4 - computer modelling of sintering: theory and examples, in  _Sintering of Advanced Materials_ , Woodhead Publishing Series in Metals and Surface Engineering, edited by Z. Z. Fang (Woodhead Publishing, 2010) pp. 86 – 109.
  * Yu  _et al._ [2017] M. Yu, S. Grasso, R. Mckinnon, T. Saunders, and M. J. Reece, Review of flash sintering: materials, mechanisms and modelling, Advances in Applied Ceramics 116, 24 (2017).
  * Z.A.  _et al._ [2006] M. Z.A., A.-T. U., and O. M., The effect of electric field and pressure on the synthesis and consolidation of materials: A review of the spark plasma sintering method., Journal of Materials Science 41, 763–777 (2006).
  * Gu  _et al._ [2012] D. D. Gu, W. Meiners, K. Wissenbach, and R. Poprawe, Laser additive manufacturing of metallic components: materials, processes and mechanisms, Int. Mater. Rev. 57, 133 (2012).
  * Pinto  _et al._ [2021] F. Pinto, R. Silva, and F. La Porta, Overview of conventional and unconventional sintering methods, in  _Green Sustainable Process for Chemical and Environmental Engineering and Science_ (Elsevier, 2021) pp. 15–35.
  * Yang  _et al._ [2020] Y. Yang, T. D. Oyedeji, P. Kühn, and B.-X. Xu, Investigation on temperature-gradient-driven effects in unconventional sintering via non-isothermal phase-field simulation, Scripta Materialia 186, 152 (2020).
  * Frenkel [1945] J. Frenkel, Viscous flow of crystalline bodies under the action of surface tension, J. phys. (USSR) 9, 385 (1945).
  * Kuczynski [1949] G. Kuczynski, Self-diffusion in sintering of metallic particles, JOM 1, 169 (1949).
  * Smith [1948] C. Smith, Grains, phases, and interfaces: an interpretation of microstructure, Trans. TMS-AIME 175, 15 (1948).
  * Coble [1961] R. L. Coble, Sintering crystalline solids. i. intermediate and final state diffusion models, Journal of applied physics 32, 787 (1961).
  * Mackenzie and Shuttleworth [1949] J. Mackenzie and R. Shuttleworth, A phenomenological theory of sintering, Proceedings of the Physical Society. Section B 62, 833 (1949).
  * Kazaryan  _et al._ [1999] A. Kazaryan, Y. Wang, and B. R. Patton, Generalized phase field approach for computer simulation of sintering: incorporation of rigid-body motion, Scr. Mater. 41, 487 (1999).
  * Wang [2006] Y. U. Wang, Computer modeling and simulation of solid-state sintering: A phase field approach, Acta Materialia 54, 953 (2006).
  * Kumar  _et al._ [2010] V. Kumar, Z. Z. Fang, and P. C. Fife, Phase field simulations of grain growth during sintering of two unequal-sized particles, Mater. Sci. Eng. A 528, 254 (2010).
  * Ahmed  _et al._ [2013] K. Ahmed, C. A. 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).
  * Deng [2012] J. Deng, A phase field model of sintering with direction-dependent diffusion, Materials Transactions 53, 385 (2012).
  * Yang  _et al._ [2018] Y. Yang, M. Yi, B.-X. Xu, and L.-Q. Chen, Phase-field modeling of non-isothermal grain coalescence in the unconventional sintering techniques (2018), arXiv:1806.02799 [cond-mat.mtrl-sci] .
  * Biswas  _et al._ [2018] S. Biswas, D. Schwen, and V. Tomar, Implementation of a phase field model for simulating evolution of two powder particles representing microstructural changes during sintering, J. Mater. Sci. 53, 5799 (2018).
  * Bis [2018] Phase field modeling of sintering: Role of grain orientation and anisotropic properties, Comput. Mater. Sci. 148, 307 (2018).
  * Chockalingam  _et al._ [2016] K. Chockalingam, V. G. Kouznetsova, O. van der Sluis, and M. G. Geers, 2D Phase field modeling of sintering of silver nanoparticles, Comput. Methods Appl. Mech. Eng. 312, 492 (2016).
  * Ahmed  _et al._ [2014] K. Ahmed, J. Pakarinen, T. Allen, and A. El-Azab, Phase field simulation of grain growth in porous uranium dioxide, J. Nucl. Mater. 446, 90 (2014).
  * Mukherjee  _et al._ [2011] R. Mukherjee, T. Chakrabarti, E. A. Anumol, T. A. Abinandanan, and N. Ravishankar, Thermal stability of spherical nanoporous aggregates and formation of hollow structures by sintering-a phase-field study, ACS Nano 5, 2700 (2011).
  * Ter [2021] Three-dimensional phase field sintering simulations accounting for the rigid-body motion of individual grains, Comput. Mater. Sci. 186, 109963 (2021).
  * Hötzer  _et al._ [2019] J. Hötzer, M. Seiz, M. Kellner, W. Rheinheimer, and B. Nestler, Phase-field simulation of solid state sintering, Acta Mater. 164, 184 (2019).
  * Greenquist  _et al._ [2020] I. Greenquist, M. R. Tonks, L. K. Aagesen, and Y. Zhang, Development of a microstructural grand potential-based sintering model, Comput. Mater. Sci. 172, 109288 (2020).
  * Yang  _et al._ [2019] Y. Yang, O. Ragnvaldsen, Y. Bai, M. Yi, and B.-X. Xu, 3D non-isothermal phase-field simulation of microstructure evolution during selective laser sintering, npj Comput. Mater. 5, 81 (2019).
  * Zhou  _et al._ [2021] X. Zhou, Y. Yang, S. Bharech, B. Lin, J. Schröder, and B.-X. Xu, 3d-multilayer simulation of microstructure and mechanical properties of porous materials by selective sintering, GAMM-Mitteilungen 44, e202100017 (2021).
  * Almgren [1999] R. Almgren, Second-order phase field asymptotics for unequal conductivities, SIAM J. Appl. Math. 59, 2086 (1999).
  * McFadden  _et al._ [2000] G. McFadden, A. Wheeler, and D. Anderson, Thin interface asymptotics for an energy/entropy approach to phase-field models with unequal conductivities, Physica D: Nonlinear Phenomena 144, 154 (2000).
  * Karma and Rappel [1996] A. Karma and W.-J. Rappel, Phase-field method for computationally efficient modeling of solidification with arbitrary interface kinetics, Phys. Rev. E 53, R3017 (1996).
  * Karma and Rappel [1998] A. Karma and W.-J. Rappel, Quantitative phase-field modeling of dendritic growth in two and three dimensions, Phys. Rev. E 57, 4323 (1998).
  * Karma [2001] A. Karma, Phase-field formulation for quantitative modeling of alloy solidification, Phys. Rev. Lett. 87, 115701 (2001).
  * Ohno and Matsuura [2009] M. Ohno and K. Matsuura, Quantitative phase-field modeling for dilute alloy solidification involving diffusion in the solid, Phys. Rev. E 79, 031603 (2009).
  * Ohno [2012] M. Ohno, Quantitative phase-field modeling of nonisothermal solidification in dilute multicomponent alloys with arbitrary diffusivities., Physical review. E, Statistical, nonlinear, and soft matter physics 86 5 Pt 1 (2012).
  * Brener and Boussinot [2012] E. A. Brener and G. Boussinot, Kinetic cross coupling between nonconserved and conserved fields in phase field models, Phys. Rev. E 86, 060601 (2012).
  * Boussinot and Brener [2013] G. Boussinot and E. A. Brener, Interface kinetics in phase-field models: Isothermal transformations in binary alloys and step dynamics in molecular-beam epitaxy, Phys. Rev. E 88, 022406 (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, Phys. Rev. E 87, 012402 (2013).
  * Boussinot and Brener [2014] G. Boussinot and E. A. Brener, Achieving realistic interface kinetics in phase-field models with a diffusional contrast, Phys. Rev. E 89, 060402 (2014).
  * Boussinot  _et al._ [2017] G. Boussinot, E. A. Brener, C. Hüter, and R. Spatschek, Elimination of surface diffusion in the non-diagonal phase field model, Continuum Mechanics and Thermodynamics 29, 969 (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, Phys. Rev. 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, Phys. Rev. B 103, 184111 (2021).
  * 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, Phys. Rev. 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, Phys. Rev. E 96, 033311 (2017).
  * Olevsky [1998] E. A. Olevsky, Theory of sintering: from discrete to continuum, Materials Science and Engineering: R: Reports 23, 41 (1998).
  * Kim  _et al._ [1999] S. G. Kim, W. T. Kim, and T. Suzuki, Phase-field model for binary alloys, Phys. Rev. E 60, 7186 (1999).
  * Echebarria  _et al._ [2004] B. Echebarria, R. Folch, A. Karma, and M. Plapp, Quantitative phase-field model of alloy solidification, Phys. Rev. E 70, 061604 (2004).
  * Maximenko and Olevsky [2004] A. L. Maximenko and E. A. Olevsky, Effective diffusion coefficients in solid-state sintering, Acta Materialia 52, 2953 (2004).
  * Balibar  _et al._ [2005] S. Balibar, H. Alles, and A. Y. Parshin, The surface of helium crystals, Rev. Mod. Phys. 77, 317 (2005).
  * Brener and Temkin [2012] E. A. Brener and D. E. Temkin, Onsager approach to the one-dimensional solidification problem and its relation to the phase-field description, Phys. Rev. E 85, 031601 (2012).
  * Nicoli  _et al._ [2011] M. Nicoli, M. Plapp, and H. Henry, Tensorial mobilities for accurate solution of transport problems in models with diffuse interfaces, Phys. Rev. E 84, 046707 (2011).
  * Tonks  _et al._ [2015] M. R. Tonks, Y. Zhang, A. Butterfield, and X.-M. Bai, Development of a grain boundary pinning model that considers particle size distribution using the phase field method, Modelling and Simulation in Materials Science and Engineering 23, 045009 (2015).
  * Zhang  _et al._ [2012] L. Zhang, M. R. Tonks, P. C. Millett, Y. Zhang, K. Chockalingam, and B. Biner, Phase-field modeling of temperature gradient driven pore migration coupling with thermal conduction, Computational materials science 56, 161 (2012).
  * Moelans  _et al._ [2008] N. Moelans, B. Blanpain, and P. Wollants, Quantitative analysis of grain boundary properties in a generalized phase field model for grain growth in anisotropic systems, Phys. Rev. B 78, 024113 (2008).
  * Jacqmin [1999] D. Jacqmin, Calculation of two-phase navier–stokes flows using phase-field modeling, Journal of computational physics 155, 96 (1999).
  * Yue  _et al._ [2004] P. Yue, J. J. Feng, C. Liu, and J. Shen, A diffuse-interface method for simulating two-phase flows of complex fluids, Journal of Fluid Mechanics 515, 293 (2004).
  * Feng  _et al._ [2005] J. J. Feng, C. Liu, J. Shen, and P. Yue, An energetic variational formulation with phase field methods for interfacial dynamics of complex fluids: advantages and challenges, in  _Modeling of soft matter_ (Springer, 2005) pp. 1–26.
  * Yue  _et al._ [2007] P. Yue, C. Zhou, and J. J. Feng, Spontaneous shrinkage of drops and mass conservation in phase-field simulations, Journal of Computational Physics 223, 1 (2007).
  * Dadvand  _et al._ [2021] A. Dadvand, M. Bagheri, N. Samkhaniani, H. Marschall, and M. Wörner, Advected phase-field method for bounded solution of the cahn-hilliard navier-stokes equations, Physics of Fluids 33, 053311 (2021).
  * Moelans  _et al._ [2009] N. Moelans, F. Wendler, and B. Nestler, Comparative study of two phase-field models for grain growth, Computational Materials Science 46, 479 (2009).
  * Aalilija  _et al._ [2021] A. Aalilija, C.-A. Gandin, and E. Hachem, A simple and efficient numerical model for thermal contact resistance based on diffuse interface immersed boundary method, Int. J. Therm. Sci. 166, 106817 (2021).
  * Yang  _et al._ [2022] Y. Yang, M. Fathidoost, T. D. Oyedeji, P. Bondi, X. Zhou, H. Egger, and B.-X. Xu, A diffuse-interface model of anisotropic interface thermal conductivity and its application in thermal homogenization of composites, Scripta Materialia 212, 114537 (2022).
  * Gugenberger  _et al._ [2008] C. Gugenberger, R. Spatschek, and K. Kassner, Comparison of phase-field models for surface diffusion, Phys. Rev. E 78, 016703 (2008).
  * Ahmed  _et al._ [2016] K. Ahmed, T. Allen, and A. El-Azab, Phase field modeling for grain growth in porous solids, Journal of Materials Science 51, 1261–1277 (2016).