# A defect-chemistry-informed phase-field model of grain growth in oxide
electroceramics

Kai Wang  Roger A. De Souza  Xiang-Long Peng  Rotraut Merkle  Wolfgang
Rheinheimer  Karsten Albe  Bai-Xiang Xu  Mechanics of Functional Materials
Division, Institute of Materials Science, Technische Universität Darmstadt,
Darmstadt, 64287, Germany  Institute of Physical Chemistry, RWTH Aachen
University, Aachen, 52056, Germany  Max Planck Institute for Solid State
Research, Heisenbergstraße 1, Stuttgart, 70569, Germany  University of
Stuttgart, Institute for Manufacturing Technologies of Ceramic Components and
Composites (IFKB), Stuttgart, 70569, Germany  Materials Modelling Division,
Institute of Materials Science, Technische Universität Darmstadt, Darmstadt,
64287, Germany

###### Abstract

Dopants can significantly affect the properties of oxide ceramics through
their impact on the property-determined microstructure characteristics such as
grain boundary segregation, space charge layer formation in the grain boundary
vicinity, and the resultant microstructure features like bimodality due to
abnormal grain growth. To support rational oxide ceramics design, we propose a
multiphysics-based and defect-chemistry-informed phase-field grain growth
model to simulate the microstructure evolution of oxide ceramics. It fully
respects the thermodynamics of charged point defects (oxygen vacancies and
dopants) in both the grain interior and boundaries and considers the competing
kinetics of defect diffusion and grain boundary movement. The proposed phase-
field model is benchmarked against well-known simplified bicrystal models,
including the Mott-Schottky and Gouy-Chapman models. Various simulation
results are presented to reveal the impacts of defect formation energy
differences between the grain interior and the grain boundary core on the key
microstructural aspects. In particular, simulation results confirm that the
solute drag effect alone can lead to bimodal grain size distribution, without
any contribution from grain misorientation and other anisotropy.
Interestingly, abnormal grain growth simulations demonstrate that grain
boundary potentials can vary substantially: grain boundaries of larger grains
tend to have lower potentials than those of smaller grains. Such heterogeneous
grain boundary potential distribution may inspire a new material optimization
strategy through microstructure design. This study provides a comprehensive
framework for defect-chemistry-consistent investigations of microstructure
evolution in polycrystalline oxide ceramics, offering fundamental insights
into in-situ processes during critical manufacturing stages.

###### keywords:

Phase-field method , Defect-chemistry , Space charge layer , Abnormal grain
growth

##  1 Introduction

In recent decades, oxide ceramics, such as SrTiO3 (STO) and BaTiO3 (BTO), have
undergone tremendous development and are increasingly crucial in applications,
such as solid oxide fuel cells [1, 2], sensors [3, 4], capacitors [5, 6],
actuators [7, 8], and memristors [9, 10]. Oxide ceramics can be fabricated via
sintering or deposition processes [11, 12, 13, 14]. When manufactured through
sintering, oxide ceramics exhibit a polycrystalline microstructure, with their
macroscopic properties heavily influenced by the synthesis conditions [15]. To
tailor the physical and chemical properties of oxide ceramics, heterovalent
substitution (doping) is an effective method for introducing additional
charged defects through ionic compensation during synthesis. The state-of-the-
art designs are mostly empirical, since knowledge of the correlation among
doping strategies, processing and the resultant microstructure is limited. The
influential mechanisms of doping on microstructure in oxide ceramics are
manifold. First of all, the addition of doping elements in the parent material
systems triggers defect chemistry reaction and results in multiple types of
charged point defects or defect complexes, which interact with the grain
boundary (GB) chemically, electrically and elastically and form space charge
layers (SCLs) around GB. More specifically, the equilibrium concentrations of
the point defects within the bulk are thermodynamically determined through
defect chemistry [16, 17]. The scenario near inside the GB core and in the GB
vicinity are different. In undoped and acceptor doped perovskites such as Fe-
doped STO, the GB core is positively charged with excess oxygen vacancies,
notated $\text{V}_{\text{O}}^{\text{..}}$ in Kröger-Vink notation [18],
stemming from the negative formation energy difference of oxygen vacancy
between GB core and bulk, i.e. $\Delta
g_{\text{V}_{\text{O}}^{..}}=\mu^{0}_{\text{V}_{\text{O}}^{..},\text{c}}-\mu^{0}_{\text{V}_{\text{O}}^{..},\text{b}}<0$.
To compensate the positively charged GB core, a negatively charged region
forms adjacent to it, called the space charge layer (SCL).

Sharp interface models have been used to describe the SCL and the related
grain boundary potential. Among them are the Mott-Schottky (MS) model and the
Gouy-Chapman (GC) model. The MS model assumes that the dopant concentration is
constant, and only the oxygen vacancies are able to redistribute, in order to
form the SCLs. In contrast, the GC model has both defects as mobile and thus
allows both to contribute to charge compensation. Unlike oxygen vacancies,
acceptor dopant forms an accumulation region within the SCL and in the GB core
due to their negative charge [Fig. 1(e)]. The spatial inhomogeneity introduced
by SCL formation exerts a profound impact on the transport dynamics of charge
species. Evidence shows that cation diffusion parallel to the GB is expedited
[19, 20], whereas charge transport across the GB is significantly hindered
[21].

SCLs determine not only the GB local properties of the synthesized material,
but also impact the grain growth dynamics already during the synthesis, e.g.,
the well-known solute drag effect. Grain growth in oxide ceramics, driven by a
reduction in surface area and controlled by the diffusion of charged species,
plays a crucial role in determining the final grain size distribution and the
concentration distributions of various point defects. Segregation in a moving
GB core necessitates additional free energy dissipation and results in the
solute drag effect [22, 23]. This effect significantly interferes with the
kinetics of grain growth and can independently trigger abnormal grain growth
even in the absence of any texture or pinning sites [24]. Dopant segregation
to GB core has been experimentally evidenced, especially at high temperatures
[25, 26]. Experiments also show that even 0.2% amount of Fe in STO can pin
certain portion of GBs [Fig. 1(b)]. This trend increases with 2% Fe [Fig.
1(c)]. When increasing the dopant concentration to 5% [Fig. 1(d)], the grain
growth is completely suppressed [25, 27]. In addition, the grain growth
coefficient is found to follow a non-Arrhenius behavior between
$1350^{\circ}$C and $1425^{\circ}$C. Two types of GB core with different
mobilities are proposed to explain the non-Arrhenius behavior [28, 25].
Despite various evidence supporting the impact of SCLs on GB properties and
grain growth dynamics, experimental characterization of the SCL remains a
challenge [29]. The complexity arises from the inherent structural and
compositional heterogeneities present in polycrystalline oxide ceramics. One
of the primary obstacles in this characterization is the limitation of in-situ
transmission electron microscopy (TEM) investigations. Such studies are
typically unfeasible due to the high sintering temperatures required for these
ceramics, which exceed the operational limits of most TEM setups.

It is hence desirable to employ theoretical models and numerical simulations
to evaluate SCLs and abnormal grain growth. De Souza derived numerical
solutions based on defect chemistry in acceptor-doped oxide systems at
equilibrium state [30]. The impacts of different charge species, the formation
energy difference of oxygen vacancy and dopant between GB core and bulk phase,
temperature and partial pressure of oxygen on grain boundary potential are
comprehensively studied in a bicrystal with planar GB core. However, applying
these numerical solutions to polycrystalline cases is impractical. Under this
circumstance, the phase-field method with a diffuse interface emerges as a
powerful tool to solve different moving boundary problems [31, 32]. There are
a series of phase-field sintering models, e.g. to consider the solute drag
effect in metals [24] and ceramics [33] and for non-isothermal sintering [34,
35, 36, 37]. But these models do not consider the charged defects,
electrostatics or defect chemistry, which are important for oxide ceramics.
Recently, Vikrant et al. developed a multi-physics phase-field model with an
emphasis on GB misorientation and investigate the SCL formation at equilibrium
and quasi-equilibrium state [38]. Then, this model is applied to explain the
abnormal grain growth with a bimodal pattern as observed in experiments [39].
However, in this phase-field model the defect chemistry is not explicitly
included, and thus there is no simulation results dedicated to the influence
of defect energetics on SCLs and abnormal grain growth.

Figure 1: Microstructures in undoped (a) and Fe-doped STO (b-d) after heating
to 1300${}^{\circ}\text{C}$ for 10h. With increasing Fe concentration, less
grain growth occurs. A fraction of very small grains appears in the
microstructures for 0.2% Fe. For 2% Fe, (c), more such small grains are
evident, until for 5% Fe, all grains remain small. Note the different scale in
(a) [25]. (e) Fe segregation to the grain boundaries as evident by STEM-EDS
[26]. Sources: 10.1016/j.jeurceramsoc.2022.11.074,
10.1016/j.actamat.2024.119941

In the present work, we develop a phase-field grain growth model founded on
the principle of defect-chemistry to simulate SCL formation and abnormal grain
growth during sintering of oxide electroceramics. By taking into account the
formation energy differences and the equilibrium defect concentrations in bulk
and GB core, the thermodynamic descriptions of point defects are included in
the phase-field model, along with the kinetics of grain boundary motion and
defect diffusion. The spatio-temporal distributions of defect concentration
and electrostatic potential in the grain interior and at GB are then solved
self-consistently.

The outline of this paper is as follows. Section 2 presents the fundamentals
of defect chemistry and the classical sharp interface bicrystal models, i.e.,
the GC and MS models, for they serve as important benchmarks in later
sections. Section 3 describes the grain growth phase-field model in the
principle of defect chemistry and electrostatics. Then, we demonstrate the
implementation of the phase-field model via the finite element method in
Section 4. The capabilities of this model are discussed for bicrystal and
polycrystalline cases in Section 5. Thereby the solute drag induced abnormal
grain growth simulation results are systematically presented and discussed,
particularly in terms of the impact of segregation energies. It is followed by
a concluding Section 6.

##  2 Defect chemistry of oxide ceramics

The basics of defect chemistry and the sharp interface bicrystal models (the
MS and the GC models) for the equilibrium state were well described in Refs.
[30, 40, 16]. Thereby the GB is treated as sharp interface with GB core
properties. Thus the concentration profiles of oxygen vacancy and dopant are
discontinuous, with concentration jumps across the GB core. For the bicrystal,
it becomes a one-dimensional case.

Assume a coordinate system $x$, with the origin $x=0$ located at the GB core.
Moreover, point defects are assumed to be dilute and non-interacting in both
GB core and bulk phases. Therefore, the electrochemical potentials of oxygen
vacancy and dopant in the bulk and in the core region can be formulated [30]
as

$\mu_{\text{V}_{\text{O}}^{..},\text{b}}(x)=\mu_{\text{V}_{\text{O}}^{..},\text{b}}^{0}+k_{\text{B}}T\ln\left[\frac{c_{\text{V}_{\text{O}}^{..},\text{b}}(x)}{N_{\text{V}_{\text{O}}^{..},\text{b}}-c_{\text{V}_{\text{O}}^{..},\text{b}}(x)}\right]+z_{\text{V}_{\text{O}}^{..}}e\phi(x),$
(1)

$\mu_{\text{dop},\text{b}}(x)=\mu_{\text{dop},\text{b}}^{0}+k_{\text{B}}T\ln\left[\frac{c_{\text{dop},\text{b}}(x)}{N_{\text{dop},\text{b}}-c_{\text{dop},\text{b}}(x)}\right]+z_{\text{dop}}e\phi(x),$
(2)

$\mu_{\text{V}_{\text{O}}^{..},\text{c}}=\mu_{\text{V}_{\text{O}}^{..},\text{c}}^{0}+k_{\text{B}}T\ln\left[\frac{c_{\text{V}_{\text{O}}^{..},\text{c}}}{N_{\text{V}_{\text{O}}^{..},\text{c}}-c_{\text{V}_{\text{O}}^{..},\text{c}}}\right]+z_{\text{V}_{\text{O}}^{..}}e\Phi_{0},$
(3)

$\mu_{\text{dop},\text{c}}=\mu_{\text{dop},\text{c}}^{0}+k_{\text{B}}T\ln\left[\frac{c_{\text{dop},\text{c}}}{N_{\text{dop},\text{c}}-c_{\text{dop},\text{c}}}\right]+z_{\text{dop}}e\Phi_{0},$
(4)

where $\mu_{\text{V}_{\text{O}}^{..},\text{b}}^{0}$,
$\mu_{\text{V}_{\text{O}}^{..},\text{c}}^{0}$, $\mu_{\text{dop},\text{b}}^{0}$
and $\mu_{\text{dop},\text{c}}^{0}$ denote the standard formation energies of
oxygen vacancy and dopant in the bulk and in the core, respectively. The
symbols $N_{\text{V}_{\text{O}}^{..},\text{b}}^{0}$,
$N_{\text{V}_{\text{O}}^{..},\text{c}}^{0}$, $N_{\text{dop},\text{b}}^{0}$ and
$N_{\text{dop},\text{c}}^{0}$ indicate the number of available sites per unit
volume of oxygen vacancy and dopant in the bulk and core, respectively.
$z_{\text{V}_{\text{O}}^{..}}$ and $z_{\text{dop}}$ denote the charge numbers
of oxygen vacancy and acceptor dopant. $\phi(x)$ is the electrostatic
potential, and $\Phi_{0}$ is the electrostatic potential in the GB core.
Similarly, $c_{\text{dop},\text{c}}$ and
$c_{\text{V}_{\text{O}}^{..},\text{c}}$ are the concentrations of dopant and
oxygen vacancy in the GB core. Thus, we have $\Phi_{0}=\phi(0)$,
$c_{\text{dop},\text{c}}=c_{\text{dop},\text{b}}(0)$ and
$c_{\text{V}_{\text{O}}^{..},\text{c}}=c_{\text{V}_{\text{O}}^{..},\text{b}}(0)$.
Additionally, $k_{\text{B}}$, $e$ and $T$ are the Boltzmann constant,
elementary charge and temperature, respectively. The formation of a SCL
requires the segregation of oxygen vacancy in the core. The smaller value of
$\mu_{\text{V}_{\text{O}}^{..},\text{c}}^{0}$ than
$\mu_{\text{V}_{\text{O}}^{..},\text{b}}^{0}$, i.e.
$\mu_{\text{V}_{\text{O}}^{..},\text{c}}^{0}-\mu_{\text{V}_{\text{O}}^{..},\text{b}}^{0}<0$,
drives the formation of SCL. At the equilibrium state, when electrochemical
potentials should be everywhere the same (in other words, no gradient), we
have

$\mu_{\text{V}_{\text{O}}^{..},\text{c}}=\mu_{\text{V}_{\text{O}}^{..},\text{b}}(x)=\mu_{\text{V}_{\text{O}}^{..},\text{b}}(\infty),$
(5)
$\mu_{\text{dop},\text{c}}=\mu_{\text{dop},\text{b}}(x)=\mu_{\text{dop},\text{b}}(\infty),$
(6)

Hence, equilibrium concentration profiles of oxygen vacancy and dopant in the
SCL take the form

$c_{\text{V}_{\text{O}}^{..},\text{b}}(x)=c_{\text{V}_{\text{O}}^{..},\text{b}}(\infty)\exp\left[\frac{-z_{\text{V}_{\text{O}}^{..}}e\\{\phi(x)-\phi(\infty)\\}}{k_{\text{B}}T}\right],$
(7)
$c_{\text{dop},\text{b}}(x)=\frac{c_{\text{dop},\text{b}}(\infty)N_{\text{dop},\text{b}}\exp\left[\frac{-z_{\text{dop}}e\\{\phi(x)-\phi(\infty)\\}}{k_{\text{B}}T}\right]}{N_{\text{dop},\text{b}}+c_{\text{dop},\text{b}}(\infty)\left\\{\exp\left[\frac{-z_{\text{dop}}e\\{\phi(x)-\phi(\infty)\\}}{k_{\text{B}}T}\right]-1\right\\}}.$
(8)

In analogy, equilibrium concentrations of oxygen vacancy and dopant in the
core are

$c_{\text{V}_{\text{O}}^{..},\text{c}}=\frac{N_{\text{V}_{\text{O}}^{..},\text{c}}c_{\text{V}_{\text{O}}^{..},\text{b}}(\infty)\exp\left[-\frac{\Delta
g_{\text{V}_{\text{O}}^{..}}+z_{\text{V}_{\text{O}}^{..}}e\Phi_{0}}{k_{\text{B}}T}\right]}{N_{\text{V}_{\text{O}}^{..},\text{b}}+c_{\text{V}_{\text{O}}^{..},\text{b}}(\infty)\exp\left[-\frac{\Delta
g_{\text{V}_{\text{O}}^{..}}+z_{\text{V}_{\text{O}}^{..}}e\Phi_{0}}{k_{\text{B}}T}\right]},$
(9)
$c_{\text{dop},\text{c}}=\frac{N_{\text{dop},\text{c}}c_{\text{dop},\text{b}}(\infty)\exp\left[-\frac{\Delta
g_{\text{dop}}+z_{\text{dop}}e\Phi_{0}}{k_{\text{B}}T}\right]}{N_{\text{dop},\text{b}}+c_{\text{dop},\text{b}}(\infty)\exp\left[-\frac{\Delta
g_{\text{dop}}+z_{\text{dop}}e\Phi_{0}}{k_{\text{B}}T}\right]},$ (10)

where $\Delta
g_{\text{V}_{\text{O}}^{..}}=\mu_{\text{V}_{\text{O}}^{..},\text{c}}^{0}-\mu_{\text{V}_{\text{O}}^{..},\text{b}}^{0}$
and $\Delta
g_{\text{dop}}=\mu_{\text{dop},\text{c}}^{0}-\mu_{\text{dop},\text{b}}^{0}$.
All parameters are tabulated in Table 1.

Table 1: Summary of the parameters in the sharp interface method Parameter | Definition | Unit
---|---|---
$\mu_{\text{V}_{\text{O}}^{..},\text{b}}^{0}$ | Formation energy of oxygen vacancy in the bulk phase | eV
$\mu_{\text{V}_{\text{O}}^{..},\text{c}}^{0}$ | Formation energy of oxygen vacancy in the core region | eV
$\mu_{\text{dop},\text{b}}^{0}$ | Formation energy of dopant in the bulk phase | eV
$\mu_{\text{dop},\text{c}}^{0}$ | Formation energy of dopant in the core region | eV
$\Delta g_{\text{V}_{\text{O}}^{..}}$ |  The difference of formation energy of oxygen vacancy between core and bulk, i.e. $\Delta g_{\text{V}_{\text{O}}^{..}}=\mu_{\text{V}_{\text{O}}^{..},\text{c}}^{0}-\mu_{\text{V}_{\text{O}}^{..},\text{b}}^{0}$ | eV
$\Delta g_{\text{dop}}$ |  The difference of formation energy of dopant between core and bulk, i.e. $\Delta g_{\text{dop}}=\mu_{\text{dop},\text{c}}^{0}-\mu_{\text{dop},\text{b}}^{0}$ | eV
$N_{\text{V}_{\text{O}}^{..},\text{b}}$ | Maximum number density of oxygen vacancy in the bulk phase | $\frac{1}{\text{m}^{3}}$
$N_{\text{V}_{\text{O}}^{..},\text{c}}$ | Maximum number density of oxygen vacancy in the core region | $\frac{1}{\text{m}^{3}}$
$N_{\text{dop},\text{b}}$ |  Maximum number density of dopant in the bulk phase. $N_{\text{dop},\text{b}}=1/a^{3}$, with $a$ being the perovskite lattice constant | $\frac{1}{\text{m}^{3}}$
$N_{\text{dop},\text{c}}$ | Maximum number density of dopant in the core region | $\frac{1}{\text{m}^{3}}$
$z_{\text{V}_{\text{O}}^{..}}$ | Valency number of oxygen vacancy, $z_{\text{V}_{\text{O}}^{..}}=2$ | -
$z_{\text{dop}}$ | Valency number of dopant, $z_{\text{dop}}=-1$ for $\text{Fe}_{\text{Ti}}^{{}^{\prime}}$ in Fe-doped STO | -
$c_{\text{dop},\text{b}}(\infty)$ | Number density of dopant in the bulk far from the core region | $\frac{1}{\text{m}^{3}}$
$c_{\text{V}_{\text{O}}^{..},\text{b}}(\infty)$ | Number density of oxygen vacancy in the bulk far from the core region | $\frac{1}{\text{m}^{3}}$
$k_{\text{B}}$ | Boltzmann constant $k_{\text{B}}=1.38\times 10^{-23}$ | $\frac{\text{J}}{\text{K}}$
$T$ | Temperature | K
$\epsilon_{0}$ | Vacuum permittivity $\epsilon_{0}=8.854\times 10^{-12}$ | $\frac{\text{C}}{\text{Vm}}$
$\epsilon_{\text{r}}$ | Relative permittivity | -
$w_{\text{c}}$ | width of the core | m

###  2.1 The Mott-Schottky model

The Mott-Schottky model assumes that the dopant is immobile and its
concentration remains fixed and constant everywhere, i.e. $\nabla
c_{\text{dop}}=0$. Note that this assumption simplifies the scenario but
breaks the equilibrium condition of the electrochemical potential of the
dopants. Based on the assumption, Poisson’s equation can be written as

$\begin{split}\epsilon_{0}\epsilon_{\text{r}}\frac{\partial^{2}\phi}{\partial
x^{2}}&=-\rho(x)\\\
&=-e[z_{\text{dop}}c_{\text{dop},\text{b}}(\infty)+z_{\text{V}_{\text{O}}^{..}}c_{\text{V}_{\text{O}}^{..},\text{b}}(x)].\end{split}$
(11)

Substituting
$c_{\text{dop},\text{b}}(\infty)=2c_{\text{V}_{\text{O}}^{..},\text{b}}(\infty)$
into Eq. $\eqref{MS equation}$, we have

$\epsilon_{0}\epsilon_{\text{r}}\frac{\partial^{2}\phi}{\partial
x^{2}}=2ec_{\text{V}_{\text{O}}^{..},\text{b}}(\infty)\left\\{1-\exp\left[\frac{-2e\\{\phi(x)-\phi(\infty)\\}}{k_{\text{B}}T}\right]\right\\}.$
(12)

With
$l_{\text{D}}=\sqrt{\frac{\epsilon_{0}\epsilon_{\text{r}}k_{\text{B}}T}{2e^{2}c_{\text{dop},\text{b}}(\infty)}}$
being the Debye length and $\alpha=\frac{e}{k_{\text{B}}T}$, Eq.$\eqref{MS2}$
becomes

$\frac{\partial^{2}\phi}{\partial x^{2}}=\frac{1}{2\alpha
l_{\text{D}}^{2}}\left\\{1-\exp\left[\frac{-2e\\{\phi(x)-\phi(\infty)\\}}{k_{\text{B}}T}\right]\right\\},$
(13)

###  2.2 The Gouy-Chapman model

In the Gouy-Chapman model both point defects (oxygen vacancies and dopant) are
mobile. The electrochemical potential of the acceptor dopant is also
considered to be constant, as is that of the oxygen vacancies. Thus, Poisson’s
equation becomes

$\begin{split}\epsilon_{0}\epsilon_{\text{r}}\frac{\partial^{2}\phi}{\partial
x^{2}}&=-\rho(x)\\\
&=-e[z_{\text{dop}}c_{\text{dop},\text{b}}(x)+z_{\text{V}_{\text{O}}^{..}}c_{\text{V}_{\text{O}}^{..},\text{b}}(x)].\end{split}$
(14)

Substituting Eq. (7) and (8) into (14), we have

$\begin{split}\epsilon_{0}\epsilon_{\text{r}}\frac{\partial^{2}\phi}{\partial
x^{2}}&=-e\left\\{z_{\text{dop}}\frac{c_{\text{dop},\text{b}}(\infty)N_{\text{dop},\text{b}}\exp\left[\frac{-z_{\text{dop}}e\\{\phi(x)-\phi(\infty)\\}}{k_{\text{B}}T}\right]}{N_{\text{dop},\text{b}}+c_{\text{dop},\text{b}}(\infty)\left\\{\exp\left[\frac{-z_{\text{dop}}e\\{\phi(x)-\phi(\infty)\\}}{k_{\text{B}}T}\right]-1\right\\}}\right.\\\
&\left.+z_{\text{V}_{\text{O}}^{..}}c_{\text{V}_{\text{O}}^{..},\text{b}}(\infty)\exp\left[\frac{-z_{\text{V}_{\text{O}}^{..}}e\\{\phi(x)-\phi(\infty)\\}}{k_{\text{B}}T}\right]\right\\}.\end{split}$
(15)

The numerical solutions of the MS and GC models have been achieved through an
iteration process [30]. In the present work, we propose an alternative
numerical method via finite element method (FEM). By considering the symmetric
configuration of a bicrysal (ignoring the grain misorientation), the semi-
infinite simulation domain of $[0,\infty)$ is sufficient. The location of the
GB core is located at the left boundary of the semi-infinite domain ($x=0$).
As for the boundary conditions, the electrostatic field is grounded at the
right side, i.e., $\phi|_{x=\infty}=0$, while at the left side, a mixed
boundary condition is given by

$\frac{\partial\phi}{\partial
x}\bigg{|}_{x=0}=-\frac{Q_{\text{c}}(\phi_{0})}{2\epsilon_{0}\epsilon_{\text{r}}},$
(16)

where
$Q_{\text{c}}(\phi_{0})=ew_{\text{c}}(z_{\text{V}_{\text{O}}^{..}}c_{\text{V}_{\text{O}}^{..},\text{c}}+z_{\text{dop}}c_{\text{dop},\text{c}})$
and $\phi_{0}=\phi|_{x=0}$. We employed the boundary element in the finite
element method to regard directly the mixed type of boundary condition defined
in the above equation. The obtained numerical solutions agree fully with these
from Ref. [30] (See Supplementary).

##  3 Phase-field grain growth model with defect chemistry and electrostatics

###  3.1 Free energy density functional and electrochemical potentials

We include here the defect chemistry into the grain growth model for oxide
ceramics, with Fe-doped STO as a prototype. Upon grain contact in Fe-doped
STO, the lower formation energy of oxygen vacancy in the core compels point
defects to segregate towards the GB core. Afterwards, the depletion or
accumulation zone appears in the vicinity of GB core. The SCL forms through a
diffusion process. The handling of SCL is similarly addressed in Ref. [30].
The segregation of point defects in the core serves to dissipate the total
free energy of the system and also induces solute drag effect when the core
moves.

Assuming an Fe-doped STO polycrystalline system consisting of $n$ grains, we
choose Kim-Kim-Suzuki (KKS) phase-field model and utilize a set of non-
conserved order parameters (OPs) ${\eta_{i}}$ to distinguish between different
grains [41]. Within grain $i$, $\eta_{i}=1$ and $\eta_{j}=0$ for $j<n$ and
$j\neq i$. Within the core region of grain $i$ and $j$, the order parameters
$\eta_{i}$ and $\eta_{j}$ spatially vary from 0 to 1. In the latter context,
$\eta$ denotes the entire set of OPs {$\eta_{i}$}. When considering
electrochemical contributions, we employ two conserved concentration fields
along with an electrostatic potential field. Other normalized field variables
are $C_{\text{V}_{\text{O}}^{..}}$ for oxygen vacancy concentration,
$C_{\text{dop}}$ for dopant concentration, and $\phi$ for the electrostatic
potential. In Fe-doped STO, $\text{Fe}^{3+}$ and $\text{Fe}^{4+}$ can coexist.
At high temperatures and low oxygen partial pressures, a higher proportion of
$\text{Fe}^{3+}$ is observed and thus the dopant
$\text{Fe}_{\text{Ti}}^{{}^{\prime}}$ has negative charge. For example, during
the sintering process at 1623 K in air, approximately 90% of the iron is
estimated to be in the $\text{Fe}^{3+}$ state, according to the bulk defect
chemical model (see Fig. A1 in Ref. [42]). Thus, in this study, we assume the
valence state of Fe to be 3+ and therefore the valance of dopant
$\text{Fe}_{\text{Ti}}^{{}^{\prime}}$ is -1, i.e. $z_{\text{dop}}=-1$. Then,
the formulation of total free energy density functional is given by

$\mathscr{F}=\int_{\Omega}(f^{\text{grad}}+f^{\text{loc}}+f^{\text{ech}})\text{d}\Omega.$
(17)

$f^{\text{grad}}$ represents the gradient term, it is given by

$f^{\text{grad}}=\frac{1}{2}\kappa\sum_{i=1}^{n}|\nabla\eta_{i}|^{2},$ (18)

where $\kappa$ is the coefficient of the gradient term. $f^{\text{loc}}$ is
the well potential. According to Ref. [43], its formulation is given by

$f^{\text{loc}}=\omega\sum_{i=1}^{n}\left(\frac{\eta_{i}^{4}}{4}-\frac{\eta_{i}^{2}}{2}+\gamma\sum_{i=1}^{n}\sum_{j>i}^{n}\eta_{i}^{2}\eta_{j}^{2}+\frac{1}{4}\right),$
(19)

with $\gamma=1.5$ for producing symmetric profile of $\eta$ and isotropic
grain growth.

The determination of parameters such as $\kappa$ and $\omega$ can be linked to
characteristics such as core width ($w_{\text{c}}$) and interface energy
($\sigma$). Their relations can be expressed as

$\sigma=\frac{\sqrt{2}}{3}\sqrt{\kappa\omega},$ (20)

and

$w_{\text{c}}=\sqrt{\frac{8\kappa}{\omega}}.$ (21)

In the present work, we only consider the isotropic interface energy, i.e. it
is not dependent on GB misorientation. The third term in free energy
functional is the electrochemical contribution. Here, we treat the
electrochemical free energy density as a mixture of bulk phase and core by an
interpolation function. Thus, $f^{\text{ech}}$ is

$f^{\text{ech}}=[1-h(\eta)]f^{\text{ech}}_{\text{b}}(C_{\text{V}_{\text{O}}^{..},\text{b}},C_{\text{dop},\text{b}})+h(\eta)f^{\text{ech}}_{\text{c}}(C_{\text{V}_{\text{O}}^{..},\text{c}},C_{\text{dop},\text{c}})-\frac{\epsilon_{0}\epsilon_{\text{r}}}{2}(\nabla\phi)^{2},$
(22)

where the subscripts b and c indicate the bulk phase and core, respectively.
$h(\eta)$ is the interpolation function and is given by

$h(\eta)=\frac{4}{3}\left[1-4\sum_{i=1}^{n}\eta_{i}^{3}+3\left(\sum_{i=1}^{n}\eta_{i}^{2}\right)^{2}\right].$
(23)

Here, $h(\eta)$ exhibits non-monotonic behavior. Specifically, for an $i/j$
GB, $h(\eta_{i})=0$ when $\eta_{i}=1$ (representing the bulk phase), whereas
$h(\eta_{i})=1$ when $\eta_{i}=0.5$ (representing the core). This
interpolation function allows us to differentiate between the bulk phase and
the core region. $\epsilon_{0}$ is the vacuum permittivity;
$\epsilon_{\text{r}}$ is the relative permittivity. For STO,
$\epsilon_{\text{r}}=9\times 10^{4}/(T-35)$ [30]. $\epsilon_{\text{r}}$ is
only dependent on temperature in the present work. Hence,
$\epsilon_{\text{r}}$ is constant across the SCL zone for specific
temperature. Furthermore,
$f^{\text{ech}}_{\text{b}}(C_{\text{V}_{\text{O}}^{..},\text{b}},C_{\text{dop},\text{b}})$
and
$f^{\text{ech}}_{\text{c}}(C_{\text{V}_{\text{O}}^{..},\text{c}},C_{\text{dop},\text{c}})$
in Eq. 22 denote the electrochemical free energy contributions from the bulk
and core, respectively, with $C_{\text{V}_{\text{O}}^{..},\text{b}}$ and
$C_{\text{dop},\text{b}}$ being the phase concentrations of oxygen vacancy and
dopant in the bulk phase, while $C_{\text{V}_{\text{O}}^{..},\text{c}}$ and
$C_{\text{dop},\text{c}}$ being the phase concentrations of oxygen vacancy and
dopant in the core. When $h(\eta)=1$, $C_{\text{dop},\text{c}}$ and
$C_{\text{V}_{\text{O}}^{..},\text{c}}$ correspond to
$c_{\text{dop},\text{c}}$ and $c_{\text{V}_{\text{O}}^{..},\text{c}}$. When
$h(\eta)=0$, $C_{\text{dop},\text{b}}$ and
$C_{\text{V}_{\text{O}}^{..},\text{b}}$ are same as
$c_{\text{dop},\text{b}}(x)$ and $c_{\text{V}_{\text{O}}^{..},\text{b}}(x)$.
Local concentration is treated as the mixture of phase concentrations in bulk
and core. The relations between local concentration and phase concentration
are expressed as

$C_{\text{V}_{\text{O}}^{..}}=[1-h(\eta)]C_{\text{V}_{\text{O}}^{..},\text{b}}+h(\eta)C_{\text{V}_{\text{O}}^{..},\text{c}},$
(24)
$C_{\text{dop}}=[1-h(\eta)]C_{\text{dop},\text{b}}+h(\eta)C_{\text{dop},\text{c}}.$
(25)

In principle of the thermodynamic model extended by Ref. [30], the expressions
for the electrochemical free energies
$f^{\text{ech}}_{\text{c}}(C_{\text{V}_{\text{O}}^{..},\text{c}},C_{\text{dop},\text{c}})$
and
$f^{\text{ech}}_{\text{b}}(C_{\text{V}_{\text{O}}^{..},\text{b}},C_{\text{dop},\text{b}})$
are given by

$\begin{split}f^{\text{ech}}_{\text{c}}(C_{\text{V}_{\text{O}}^{..},\text{c}},C_{\text{dop},\text{c}})&=\frac{1}{V_{\text{m}}}\left\\{g^{0}_{\text{V}_{\text{O}}^{..},\text{c}}C_{\text{V}_{\text{O}}^{..},\text{c}}+g^{0}_{\text{dop},\text{c}}C_{\text{dop},\text{c}}\right.\\\
&+RT\left[C_{\text{V}_{\text{O}}^{..},\text{c}}\ln(C_{\text{V}_{\text{O}}^{..},\text{c}})+(\tilde{N}_{\text{V}_{\text{O}}^{..},\text{c}}-C_{\text{V}_{\text{O}}^{..},\text{c}})\ln(\tilde{N}_{\text{V}_{\text{O}}^{..},\text{c}}-C_{\text{V}_{\text{O}}^{..},\text{c}})\right.\\\
&\left.\left.+C_{\text{dop},\text{c}}\ln(C_{\text{dop},\text{c}})+(\tilde{N}_{\text{dop},\text{c}}-C_{\text{dop},\text{c}})\ln(\tilde{N}_{\text{dop},\text{c}}-C_{\text{dop},\text{c}})\right]\right\\}\\\
&+\frac{\mathcal{F}}{V_{\text{m}}}\left(z_{\text{V}_{\text{O}}^{..}}C_{\text{V}_{\text{O}}^{..},\text{c}}+z_{\text{dop}}C_{\text{dop},\text{c}}\right)\phi,\end{split}$
(26)
$\begin{split}f^{\text{ech}}_{\text{b}}(C_{\text{V}_{\text{O}}^{..},\text{b}},C_{\text{dop},\text{b}})&=\frac{1}{V_{\text{m}}}\left\\{g^{0}_{\text{V}_{\text{O}}^{..},\text{b}}C_{\text{V}_{\text{O}}^{..},\text{b}}+g^{0}_{\text{dop},\text{b}}C_{\text{dop},\text{b}}\right.\\\
&+RT\left[C_{\text{V}_{\text{O}}^{..},\text{b}}\ln(C_{\text{V}_{\text{O}}^{..},\text{b}})+(\tilde{N}_{\text{V}_{\text{O}}^{..},\text{b}}-C_{\text{V}_{\text{O}}^{..},\text{b}})\ln(\tilde{N}_{\text{V}_{\text{O}}^{..},\text{b}}-C_{\text{V}_{\text{O}}^{..},\text{b}})\right.\\\
&\left.\left.+C_{\text{dop},\text{b}}\ln(C_{\text{dop},\text{b}})+(\tilde{N}_{\text{dop},\text{b}}-C_{\text{dop},\text{b}})\ln(\tilde{N}_{\text{dop},\text{b}}-C_{\text{dop},\text{b}})\right]\right\\}\\\
&+\frac{\mathcal{F}}{V_{\text{m}}}\left(z_{\text{V}_{\text{O}}^{..}}C_{\text{V}_{\text{O}}^{..},\text{b}}+z_{\text{dop}}C_{\text{dop},\text{b}}\right)\phi,\end{split}$
(27)

where $g^{0}_{\text{V}_{\text{O}}^{..},\text{c}}$,
$g^{0}_{\text{V}_{\text{O}}^{..},\text{b}}$, $g^{0}_{\text{dop},\text{c}}$,
$g^{0}_{\text{dop},\text{b}}$ are standard formation energies of oxygen
vacancy and dopant in bulk and core.
$\tilde{N}_{\text{V}_{\text{O}}^{..},\text{c}}$,
$\tilde{N}_{\text{V}_{\text{O}}^{..},\text{b}}$,
$\tilde{N}_{\text{dop},\text{c}}$, $\tilde{N}_{\text{dop},\text{b}}$ are
dimensionless numbers of available sites per unit volume for oxygen vacancy
and dopant in bulk and core. $z_{\text{V}_{\text{O}}^{..}}$ and
$z_{\text{dop}}$ are valance numbers of oxygen vacancy and dopant.
$\mathcal{F}$, $V_{\text{m}}$, $R$, $T$ are the Faraday constant, molar
volume, gas constant and temperature, respectively. The electrochemical
potentials of oxygen vacancy and dopant in bulk and core are obtained by

$\mu_{\text{V}_{\text{O}}^{..},\text{b}}^{\text{PF}}=\frac{\delta\mathscr{F}}{\delta
C_{\text{V}_{\text{O}}^{..},\text{b}}}=\frac{g_{\text{V}_{\text{O}}^{..},\text{b}}^{0}}{V_{\text{m}}}+\frac{RT}{V_{\text{m}}}\ln{\left[\frac{C_{\text{V}_{\text{O}}^{..},\text{b}}}{\tilde{N}_{\text{V}_{\text{O}}^{..},\text{b}}-C_{\text{V}_{\text{O}}^{..},\text{b}}}\right]}+\frac{\mathcal{F}}{V_{\text{m}}}z_{\text{V}_{\text{O}}^{..}}\phi,$
(28)
$\mu_{\text{V}_{\text{O}}^{..},\text{c}}^{\text{PF}}=\frac{\delta\mathscr{F}}{\delta
C_{\text{V}_{\text{O}}^{..},\text{c}}}=\frac{g_{\text{V}_{\text{O}}^{..},\text{c}}^{0}}{V_{\text{m}}}+\frac{RT}{V_{\text{m}}}\ln{\left[\frac{C_{\text{V}_{\text{O}}^{..},\text{c}}}{\tilde{N}_{\text{V}_{\text{O}}^{..},\text{c}}-C_{\text{V}_{\text{O}}^{..},\text{c}}}\right]}+\frac{\mathcal{F}}{V_{\text{m}}}z_{\text{V}_{\text{O}}^{..}}\phi,$
(29) $\mu_{\text{dop},\text{b}}^{\text{PF}}=\frac{\delta\mathscr{F}}{\delta
C_{\text{dop},\text{b}}}=\frac{g_{\text{dop},\text{b}}^{0}}{V_{\text{m}}}+\frac{RT}{V_{\text{m}}}\ln{\left[\frac{C_{\text{dop},\text{b}}}{\tilde{N}_{\text{dop},\text{b}}-C_{\text{dop},\text{b}}}\right]}+\frac{\mathcal{F}}{V_{\text{m}}}z_{\text{dop}}\phi,$
(30) $\mu_{\text{dop},\text{c}}^{\text{PF}}=\frac{\delta\mathscr{F}}{\delta
C_{\text{dop},\text{c}}}=\frac{g_{\text{dop},\text{c}}^{0}}{V_{\text{m}}}+\frac{RT}{V_{\text{m}}}\ln{\left[\frac{C_{\text{dop},\text{c}}}{\tilde{N}_{\text{dop},\text{c}}-C_{\text{dop},\text{c}}}\right]}+\frac{\mathcal{F}}{V_{\text{m}}}z_{\text{dop}}\phi,$
(31)

which are consistent with the standard formations of electrochemical
potentials in defect chemistry, i.e. Eq. 1 to Eq. 4.

###  3.2 Governing equations

In the framework of KKS phase-field model, the evolution of phase-field OPs
obeys the Allen-Cahn equation. For grain $i$, the evolution of $\eta_{i}$ in
time is

$\frac{1}{L}\frac{\partial\eta_{i}}{\partial
t}=\frac{\delta\mathscr{F}}{\delta\eta_{i}},$ (32)

where $L$ is the GB mobility. Substituting Eq. (17) into (32), we can obtain

$\begin{split}\frac{1}{L}\frac{\partial\eta_{i}}{\partial
t}&=\nabla\cdot\kappa\nabla\eta_{i}-\omega\frac{\partial
f^{\text{loc}}(\eta)}{\partial\eta_{i}}\\\ &+\frac{\partial
h(\eta)}{\partial\eta_{i}}\left[f_{\text{b}}^{\text{ech}}-f_{\text{c}}^{\text{ech}}-\frac{\partial
f_{\text{b}}^{\text{ech}}}{\partial
C_{\text{V}_{\text{O}}^{..},\text{b}}}(C_{\text{V}_{\text{O}}^{..},\text{b}}-C_{\text{V}_{\text{O}}^{..},\text{c}})-\frac{\partial
f_{\text{b}}^{\text{ech}}}{\partial
C_{\text{dop},\text{b}}}(C_{\text{dop},\text{b}}-C_{\text{dop},\text{c}})\right].\end{split}$
(33)

At the quasi-equilibrium state, the core moves with a constant velocity,
$v_{c}$, we have $\dot{\eta}=-v_{c}\eta^{\prime}(x)$, thus Eq. (32) becomes

$\begin{split}\frac{v_{c}}{L}\frac{d\eta}{dx}&=\nabla\cdot\kappa\nabla\eta_{i}-\omega\frac{\partial
f^{\text{loc}}(\eta)}{\partial\eta_{i}}\\\ &+\frac{\partial
h(\eta)}{\partial\eta_{i}}\left[f_{\text{b}}^{\text{ech}}-f_{\text{c}}^{\text{ech}}-\frac{\partial
f_{\text{b}}^{\text{ech}}}{\partial
C_{\text{V}_{\text{O}}^{..},\text{b}}}(C_{\text{V}_{\text{O}}^{..},\text{b}}-C_{\text{V}_{\text{O}}^{..},\text{c}})-\frac{\partial
f_{\text{b}}^{\text{ech}}}{\partial
C_{\text{dop},\text{b}}}(C_{\text{dop},\text{b}}-C_{\text{dop},\text{c}})\right],\end{split}$
(34)

Conserved concentration fields $C_{\text{V}_{\text{O}}^{..}}$ and
$C_{\text{dop}}$ follow Fick’s law, which are expressed as

$\frac{\partial C_{\text{V}_{\text{O}}^{..}}}{\partial
t}=\nabla\left(M_{\text{V}_{\text{O}}^{..}}\nabla\frac{\delta\mathscr{F}}{\delta
C_{\text{V}_{\text{O}}^{..}}}\right),$ (35) $\frac{\partial
C_{\text{dop}}}{\partial
t}=\nabla\left(M_{\text{dop}}\nabla\frac{\delta\mathscr{F}}{\delta
C_{\text{dop}}}\right),$ (36)

where $M_{\text{V}_{\text{O}}^{..}}$ and $M_{\text{dop}}$ are mobilities of
oxygen vacancy and dopant, which are given by
$M_{\text{V}_{\text{O}}^{..}}=\frac{D_{\text{V}_{\text{O}}^{..}}}{\partial^{2}f^{\text{ech}}/\partial
C_{\text{V}_{\text{O}}^{..}}^{2}}$ and
$M_{\text{dop}}=\frac{D_{\text{dop}}}{\partial^{2}f^{\text{ech}}/\partial
C_{\text{dop}}^{2}}$, with $D_{\text{V}_{\text{O}}^{..}}$ and $D_{\text{dop}}$
being the diffusivities of oxygen vacancy and dopant, respectively. Moreover,
local equilibrium between bulk phases and core should be satisfied everywhere.
Therefore, we have following constrains

$\delta\mathscr{F}/\delta
C_{\text{V}_{\text{O}}^{..}}=\delta\mathscr{F}/\delta
C_{\text{V}_{\text{O}}^{..},\text{b}}=\delta\mathscr{F}/\delta
C_{\text{dop},\text{b}},$ (37) $\delta\mathscr{F}/\delta
C_{\text{dop}}=\delta\mathscr{F}/\delta
C_{\text{V}_{\text{O}}^{..},\text{c}}=\delta\mathscr{F}/\delta
C_{\text{dop},\text{c}}.$ (38)

Therefore, $\delta\mathscr{F}/\delta C_{\text{V}_{\text{O}}^{..}}$ and
$\delta\mathscr{F}/\delta C_{\text{dop}}$ in Eqs. 35 and 36 can be substituted
by $\mu_{\text{V}_{\text{O}}^{..},\text{b}}^{\text{PF}}$,
$\mu_{\text{dop},\text{b}}^{\text{PF}}$ or
$\mu_{\text{V}_{\text{O}}^{..},\text{c}}^{\text{PF}}$,
$\mu_{\text{dop},\text{c}}^{\text{PF}}$. In addition, the governing equation
of electrostatic potential field is obtained by Poisson’s equation

$\epsilon_{0}\epsilon_{\text{r}}\nabla^{2}\phi+\frac{\mathcal{F}}{V_{\text{m}}}(z_{\text{V}_{\text{O}}^{..}}C_{\text{V}_{\text{O}}^{..}}+z_{\text{dop}}C_{\text{dop}})=0,$
(39)

The parameters utilized in phase-field model are tabulated in table 2.

Table 2: Summary of the parameters in the phase-field model DOFs | Definition | Unit
---|---|---
$\eta_{i}$ | phase-field order parameter of grain $i$, $i=1,2,3,...,n$ | -
$C_{\text{V}_{\text{O}}^{..}}$ |  Normalized local concentration of oxygen vacancy ($C_{\text{V}_{\text{O}}^{..}}=a^{3}c_{\text{V}_{\text{O}}^{..}}$, with $a$ being the lattice constant) | -
$C_{\text{dop}}$ | Normalized local concentration of dopant, $C_{\text{dop}}=a^{3}c_{\text{dop}}$ | -
$C_{\text{V}_{\text{O}}^{..},\text{b}}$ | Normalized concentration of oxygen vacancy in bulk phase, $C_{\text{V}_{\text{O}}^{..},\text{b}}=a^{3}c_{\text{V}_{\text{O}}^{..},\text{b}}$ | -
$C_{\text{V}_{\text{O}}^{..},\text{c}}$ | Normalized concentration of oxygen vacancy in core region, $C_{\text{V}_{\text{O}}^{..},\text{c}}=a^{3}c_{\text{V}_{\text{O}}^{..},\text{c}}$ | -
$C_{\text{dop},\text{b}}$ | Normalized concentration of dopant in bulk phase, $C_{\text{dop},\text{b}}=a^{3}c_{\text{dop},\text{b}}$ | -
$C_{\text{dop},\text{c}}$ | Normalized concentration of dopant in core region, $C_{\text{dop},\text{c}}=a^{3}c_{\text{dop},\text{c}}$ | -
${\phi}$ | Electrostatic potential | V
Parameters | Definition | Unit
$\kappa$ | Coefficient of gradient term | $\frac{\text{J}}{\text{m}}$
$\omega$ | Heigh of well potential | $\frac{\text{J}}{\text{m}^{3}}$
$L$ | GB mobility | $\frac{\text{m}^{3}}{\text{Js}}$
$g^{0}_{\text{V}_{\text{O}}^{..},\text{b}}$ |  Formation energy of oxygen vacancy in the bulk phase. $g^{0}_{\text{V}_{\text{O}}^{..},\text{b}}=\mu_{\text{V}_{\text{O}}^{..},\text{b}}^{0}N_{\text{A}}$, with $N_{\text{A}}$ being Avogadro constant | $\frac{\text{J}}{\text{mol}}$
$g_{\text{V}_{\text{O}}^{..},\text{c}}^{0}$ |  Formation energy of oxygen vacancy in the core region. $g^{0}_{\text{V}_{\text{O}}^{..},\text{c}}=\mu_{\text{V}_{\text{O}}^{..},\text{c}}^{0}N_{\text{A}}$ | $\frac{\text{J}}{\text{mol}}$
$g^{0}_{\text{dop},\text{b}}$ |  Formation energy of dopant in the bulk phase. $g^{0}_{\text{dop},\text{b}}=\mu_{\text{dop},\text{b}}^{0}N_{\text{A}}$ | $\frac{\text{J}}{\text{mol}}$
$g^{0}_{\text{dop},\text{c}}$ |  Formation energy of dopant in the core region. $g^{0}_{\text{dop},\text{c}}=\mu_{\text{dop},\text{c}}^{0}N_{\text{A}}$ | $\frac{\text{J}}{\text{mol}}$
$V_{\text{m}}$ | Molar volume | $\frac{\text{m}^{3}}{\text{mol}}$
$\tilde{N}_{\text{V}_{\text{O}}^{..},\text{c}}$ |  Dimensionless available sites of oxygen vacancy in core region. $\tilde{N}_{\text{V}_{\text{O}}^{..},\text{c}}=N_{\text{V}_{\text{O}}^{..},\text{c}}a^{3}$, with $a$ being the lattice constant | -
$\tilde{N}_{\text{V}_{\text{O}}^{..},\text{b}}$ |  Dimensionless available sites of oxygen vacancy in bulk phase. $\tilde{N}_{\text{V}_{\text{O}}^{..},\text{b}}=N_{\text{V}_{\text{O}}^{..},\text{b}}a^{3}$ | -
$\tilde{N}_{\text{dop},\text{c}}$ |  Dimensionless available sites of dopant in core region. $\tilde{N}_{\text{dop},\text{c}}=N_{\text{dop},\text{c}}a^{3}$ | -
$\tilde{N}_{\text{dop},\text{b}}$ |  Dimensionless available sites of dopant in bulk phase. $\tilde{N}_{\text{dop},\text{b}}=N_{\text{dop},\text{b}}a^{3}$ | -
$R$ | Gas constant | $\frac{\text{J}}{\text{K}}$
$T$ | Temperature | K
$\mathcal{F}$ | Faraday constant | $\frac{\text{C}}{\text{mol}}$
$D_{\text{V}_{\text{O}}^{..}}$ | Diffusivity of oxygen vacancy | $\frac{\text{m}^{2}}{\text{s}}$
$D_{\text{dop}}$ | Diffusivity of dopant | $\frac{\text{m}^{2}}{\text{s}}$

##  4 Finite element implementation of the phase-field model

In this section, we present the detailed finite element implementation for the
defect-chemistry informed phase-field model. According to the model presented
in Section 3 and Table 2, we treat the OPs $\eta_{i}$ and the field variables
$C_{\text{V}_{\text{O}}^{..}}$, $C_{\text{dop}}$,
$C_{\text{V}_{\text{O}}^{..},\text{b}}$,
$C_{\text{V}_{\text{O}}^{..},\text{c}}$, $C_{\text{dop},\text{b}}$,
$C_{\text{dop},\text{c}}$, $\phi$ as degrees-of-freedom (DOFs). Moreover, to
solve Eqs. (35) and (36), we use the mixed finite element formulation and
treat the electro-chemical potentials $\mu_{\text{V}_{\text{O}}^{..}}$ and
$\mu_{\text{dop}}$ as additional DOFs [44].

The strong forms of the corresponding governing equations have already been
given in Eqs. (24), (25) and (33) to (39). According to the FEM [45], we
formulate the residuals from the weak forms of these governing equations by
introducing corresponding test functions
$\psi_{C_{\text{V}_{\text{O}}^{..},\text{c}}}$,
$\psi_{C_{\text{dop},\text{c}}}$, $\psi_{\eta_{i}}$,
$\psi_{C_{\text{V}_{\text{O}}^{..}}}$, $\psi_{C_{\text{dop}}}$,
$\psi_{\mu_{\text{V}_{\text{O}}^{..}}}$, $\psi_{\mu_{\text{dop}}}$,
$\psi_{C_{\text{V}_{\text{O}}^{..},\text{b}}}$,
$\psi_{C_{\text{dop},\text{b}}}$ and $\psi_{\phi}$. Firstly, the residuals for
the weak forms of Eqs. (24), (25) and (33) are

$R_{C_{\text{V}_{\text{O}}^{..},\text{c}}}=\int_{V}\psi_{C_{\text{V}_{\text{O}}^{..},\text{c}}}\\{C_{\text{V}_{\text{O}}^{..}}-[1-h(\eta)]C_{\text{V}_{\text{O}}^{..},\text{b}}-h(\eta)C_{\text{V}_{\text{O}}^{..},\text{c}}\\}\text{d}V,$
(40)

$R_{C_{\text{dop},\text{c}}}=\int_{V}\psi_{C_{\text{dop},\text{c}}}\\{C_{\text{dop}}-[1-h(\eta)]C_{\text{dop},\text{b}}-h(\eta)C_{\text{dop},\text{c}}\\}\text{d}V,$
(41)

$\begin{split}R_{\eta_{i}}&=\int_{V}{\psi_{\eta_{i}}}\frac{\eta_{i}^{t_{n+1}}-\eta_{i}^{t_{n}}}{\Delta
t}\text{d}V+L\int_{V}\psi_{\eta_{i}}\omega\frac{\partial
f^{\text{loc}}(\eta)}{\partial\eta_{i}}\text{d}V+L\int_{V}\kappa\eta_{i,j}\psi_{\eta_{i},j}\text{d}V-L\int_{\Gamma}\kappa\eta_{i,j}\psi_{\eta_{i}}\hat{n}_{j}\text{d}\Gamma\\\
&-L\int_{V}\psi_{\eta_{i}}\frac{\partial
h(\eta)}{\partial\eta_{i}}\left[f_{\text{b}}^{\text{ech}}-f_{\text{c}}^{\text{ech}}-\frac{\partial
f_{\text{b}}^{\text{ech}}}{\partial
C_{\text{V}_{\text{O}}^{..},\text{b}}}(C_{\text{V}_{\text{O}}^{..},\text{b}}-C_{\text{V}_{\text{O}}^{..},\text{c}})-\frac{\partial
f_{\text{b}}^{\text{ech}}}{\partial
C_{\text{dop},\text{b}}}(C_{\text{dop},\text{b}}-C_{\text{dop},\text{c}})\right]\text{d}V.\end{split}$
(42)

Note that we use the backward Euler method for the time integration. In the
very last equation, $\eta_{i}^{t_{n}}$ and $\eta_{i}^{t_{n+1}}$ indicate
$\eta_{i}$ at the time step $t_{n}$ and $t_{n+1}$, respectively. Thereby,
$\Delta t=t_{n+1}-t_{n}$, and $\hat{n}_{i}$ is the normal vector to the
boundary $\Gamma$ of the subdomain. Then, through introducing the addition
coupling fields $\mu_{\text{V}_{\text{O}}^{..}}=\delta\mathscr{F}/\delta
C_{\text{V}_{\text{O}}^{..}}$ and $\mu_{\text{dop}}=\delta\mathscr{F}/\delta
C_{\text{dop}}$, we obtain the weak forms of Eqs. (35) and (36) as

$R_{C_{\text{V}_{\text{O}}^{..}}}=\int_{V}\psi_{C_{\text{V}_{\text{O}}^{..}}}\left(\mu_{\text{V}_{\text{O}}^{..}}-\frac{\delta\mathscr{F}}{\delta
C_{\text{V}_{\text{O}}^{..},\text{b}}}\right)\text{d}V.$ (43)

$R_{C_{\text{dop}}}=\int_{V}\psi_{C_{\text{dop}}}\left(\mu_{\text{dop}}-\frac{\delta\mathscr{F}}{\delta
C_{\text{dop},\text{b}}}\right)\text{d}V.$ (44)

$R_{\mu_{\text{V}_{\text{O}}^{..}}}=\int_{V}\psi_{\mu_{\text{V}_{\text{O}}^{..}}}\frac{C_{\text{V}_{\text{O}}^{..}}^{t_{n+1}}-C_{\text{V}_{\text{O}}^{..}}^{t_{n}}}{\Delta
t}\text{d}V+\int_{V}M_{\text{V}_{\text{O}}^{..}}\mu_{V_{O}^{..},i}\psi_{\mu_{\text{V}_{\text{O}}^{..}},i}\text{d}V-\int_{\Gamma}M_{\text{V}_{\text{O}}^{..}}\mu_{V_{O}^{..},i}\psi_{\mu_{\text{V}_{\text{O}}^{..}}}\hat{n}_{i}\text{d}\Gamma.$
(45)

$R_{\mu_{\text{dop}}}=\int_{V}\psi_{\mu_{\text{dop}}}\frac{C_{\text{dop}}^{t_{n+1}}-C_{\text{dop}}^{t_{n}}}{\Delta
t}\text{d}V+\int_{V}M_{\text{dop}}\mu_{dop,i}\psi_{\mu_{\text{dop}},i}\text{d}V-\int_{\Gamma}M_{\text{dop}}\mu_{dop,i}\psi_{\mu_{\text{dop}}}\hat{n}_{i}\text{d}\Gamma.$
(46)

The weak forms of the local equilibrium constrain in Eqs. (37) and (38) are

$R_{C_{\text{V}_{\text{O}}^{..},\text{b}}}=\int_{V}\psi_{C_{\text{V}_{\text{O}}^{..},\text{b}}}\left(\frac{\delta\mathscr{F}}{\delta
C_{\text{V}_{\text{O}}^{..},\text{b}}}-\frac{\delta\mathscr{F}}{\delta
C_{\text{V}_{\text{O}}^{..},\text{c}}}\right)\text{d}V.$ (47)

$R_{C_{\text{dop},\text{b}}}=\int_{V}\psi_{C_{\text{dop},\text{b}}}\left(\frac{\delta\mathscr{F}}{\delta
C_{\text{dop},\text{b}}}-\frac{\delta\mathscr{F}}{\delta
C_{\text{dop},\text{c}}}\right)\text{d}V.$ (48)

Finally, the weak form of Eq (39) is

$R_{\phi}=\int_{V}\epsilon_{0}\epsilon_{\text{r}}\phi_{,i}\psi_{\phi,i}\text{d}V-\int_{\Gamma}\epsilon_{0}\epsilon_{\text{r}}\phi_{,i}\psi_{\phi}\hat{n}_{i}\text{d}\Gamma-\int_{V}\psi_{\phi}\frac{\mathcal{F}}{V_{\text{m}}}(z_{\text{V}_{\text{O}}^{..}}C_{\text{V}_{\text{O}}^{..}}+z_{\text{dop}}C_{\text{dop}})\text{d}V.$
(49)

Following the Galerkin approach [46], the test functions are discretized as
$\psi_{C_{\text{V}_{\text{O}}^{..},\text{c}}}=N_{C_{\text{V}_{\text{O}}^{..},\text{c}}}^{I}\psi_{C_{\text{V}_{\text{O}}^{..},\text{c}}}^{I}$,
$\psi_{C_{\text{dop},\text{c}}}=N_{C_{\text{dop},\text{c}}}^{I}\psi_{C_{\text{dop},\text{c}}}^{I}$,
$\psi_{\eta_{i}}=N_{\eta_{i}}^{I}\psi_{\eta_{i}}^{I}$,
$\psi_{C_{\text{V}_{\text{O}}^{..}}}=N_{C_{\text{V}_{\text{O}}^{..}}}^{I}\psi_{C_{\text{V}_{\text{O}}^{..}}}^{I}$,
$\psi_{C_{\text{dop}}}=N_{C_{\text{dop}}}\psi_{C_{\text{dop}}}^{I}$,
$\psi_{\mu_{\text{V}_{\text{O}}^{..}}}=N_{\mu_{\text{V}_{\text{O}}^{..}}}\psi_{\mu_{\text{V}_{\text{O}}^{..}}}^{I}$,
$\psi_{\mu_{\text{dop}}}=N_{\mu_{\text{dop}}}\psi_{\mu_{\text{dop}}}^{I}$,
$\psi_{C_{\text{V}_{\text{O}}^{..},\text{b}}}=N_{C_{\text{V}_{\text{O}}^{..},\text{b}}}\psi_{C_{\text{V}_{\text{O}}^{..},\text{b}}}^{I}$,
$\psi_{C_{\text{dop},\text{b}}}=N_{C_{\text{dop},\text{b}}}\psi_{C_{\text{dop},\text{b}}}^{I}$
and $\psi_{\phi}=N_{\phi}\psi_{\phi}^{I}$, where $I$ denotes the node index
and Einstein summation convention is used.
$\psi_{C_{\text{V}_{\text{O}}^{..},\text{c}}}^{I}$,
$\psi_{C_{\text{dop},\text{c}}}^{I}$, $\psi_{\eta_{i}}^{I}$,
$\psi_{C_{\text{V}_{\text{O}}^{..}}}^{I}$, $\psi_{C_{\text{dop}}}^{I}$,
$\psi_{\mu_{\text{V}_{\text{O}}^{..}}}^{I}$, $\psi_{\mu_{\text{dop}}}^{I}$,
$\psi_{C_{\text{V}_{\text{O}}^{..},\text{b}}}^{I}$,
$\psi_{C_{\text{dop},\text{b}}}^{I}$ and $\psi_{\phi}^{I}$ are the nodal
weights, while $N_{C_{\text{V}_{\text{O}}^{..},\text{c}}}^{I}$,
$N_{C_{\text{dop},\text{c}}}^{I}$, $N_{\eta_{i}}^{I}$,
$N_{C_{\text{V}_{\text{O}}^{..}}}^{I}$, $N_{C_{\text{dop}}}^{I}$,
$N_{\mu_{\text{V}_{\text{O}}^{..}}}^{I}$, $N_{\mu_{\text{dop}}}^{I}$,
$N_{C_{\text{V}_{\text{O}}^{..},\text{b}}}^{I}$,
$N_{C_{\text{dop},\text{b}}}^{I}$ and $N_{\phi}^{I}$ are the shape functions
for the corresponding variables. Then we can write down the discretized
residuals as

$\displaystyle
R_{C_{\text{V}_{\text{O}}^{..},\text{c}}}^{I}=\int_{V}N_{C_{\text{V}_{\text{O}}^{..},\text{c}}}^{I}\\{C_{\text{V}_{\text{O}}^{..}}-[1-h(\eta)]C_{\text{V}_{\text{O}}^{..},\text{b}}-h(\eta)C_{\text{V}_{\text{O}}^{..},\text{c}}\\}\text{d}V$
(50) $\displaystyle
R_{C_{\text{dop},\text{c}}}^{I}=\int_{V}N_{C_{\text{dop},\text{c}}}^{I}\\{C_{\text{dop}}-[1-h(\eta)]C_{\text{dop},\text{b}}-h(\eta)C_{\text{dop},\text{c}}\\}\text{d}V$
(51) $\displaystyle\begin{aligned}
R_{\eta_{i}}^{I}&=\int_{V}\frac{N_{\eta_{i}}^{I}}{\Delta
t}\text{d}V+L\int_{V}N_{\eta_{i}}^{I}\omega\frac{\partial
f^{\text{loc}}(\eta)}{\partial\eta_{i}}\text{d}V+L\int_{V}\kappa\eta_{i,j}N_{\eta_{i},j}^{I}\text{d}V-L\int_{\Gamma}\kappa\eta_{i,j}N_{\eta_{i}}^{I}\hat{n}_{j}\text{d}\Gamma\\\
&-L\int_{V}N_{\eta_{i}}^{I}\frac{\partial
h(\eta)}{\partial\eta_{i}}\left[f_{\text{b}}^{\text{ech}}-f_{\text{c}}^{\text{ech}}-\frac{\partial
f_{\text{b}}^{\text{ech}}}{\partial
C_{\text{V}_{\text{O}}^{..},\text{b}}}(C_{\text{V}_{\text{O}}^{..},\text{b}}-C_{\text{V}_{\text{O}}^{..},\text{c}})-\frac{\partial
f_{\text{b}}^{\text{ech}}}{\partial
C_{\text{dop},\text{b}}}(C_{\text{dop},\text{b}}-C_{\text{dop},\text{c}})\right]\text{d}V\end{aligned}$
(52) $\displaystyle
R_{C_{\text{V}_{\text{O}}^{..}}}^{I}=\int_{V}N_{C_{\text{V}_{\text{O}}^{..}}}^{I}\left(\mu_{\text{V}_{\text{O}}^{..}}-\frac{\delta\mathscr{F}}{\delta
C_{\text{V}_{\text{O}}^{..},\text{b}}}\right)\text{d}V$ (53) $\displaystyle
R_{C_{\text{dop}}}^{I}=\int_{V}N_{C_{\text{dop}}}^{I}\left(\mu_{\text{dop}}-\frac{\delta\mathscr{F}}{\delta
C_{\text{dop},\text{b}}}\right)\text{d}V$ (54) $\displaystyle
R_{\mu_{\text{V}_{\text{O}}^{..}}}^{I}=\int_{V}\frac{N_{\mu_{\text{V}_{\text{O}}^{..}}}^{I}}{\Delta
t}\text{d}V+\int_{V}M_{\text{V}_{\text{O}}^{..}}\mu_{\text{V}_{\text{O}}^{..},i}N_{\mu_{\text{V}_{\text{O}}^{..}},i}^{I}\text{d}V-\int_{\Gamma}M_{\text{V}_{\text{O}}^{..}}\mu_{V_{O}^{..},i}N_{\mu_{\text{V}_{\text{O}}^{..}}}^{I}\hat{n}_{i}\text{d}\Gamma$
(55) $\displaystyle
R_{\mu_{\text{dop}}}^{I}=\int_{V}\frac{N_{\mu_{\text{dop}}}^{I}}{\Delta
t}\text{d}V+\int_{V}M_{\text{dop}}\mu_{\text{dop},i}N_{\mu_{\text{dop}},i}^{I}\text{d}V-\int_{\Gamma}M_{\text{dop}}\mu_{\text{dop},i}N_{\mu_{\text{dop}}}^{I}\hat{n}_{i}\text{d}\Gamma$
(56) $\displaystyle
R_{C_{\text{V}_{\text{O}}^{..},\text{b}}}^{I}=\int_{V}N_{C_{\text{V}_{\text{O}}^{..},\text{b}}}^{I}\left(\frac{\delta\mathscr{F}}{\delta
C_{\text{V}_{\text{O}}^{..},\text{b}}}-\frac{\delta\mathscr{F}}{\delta
C_{\text{V}_{\text{O}}^{..},\text{c}}}\right)\text{d}V$ (57) $\displaystyle
R_{C_{\text{dop},\text{b}}}^{I}=\int_{V}N_{C_{\text{dop},\text{b}}}^{I}\left(\frac{\delta\mathscr{F}}{\delta
C_{\text{dop},\text{b}}}-\frac{\delta\mathscr{F}}{\delta
C_{\text{dop},\text{c}}}\right)\text{d}V$ (58) $\displaystyle
R_{\phi}^{I}=\int_{V}\epsilon_{0}\epsilon_{\text{r}}\phi_{,i}N_{\phi,i}^{I}\text{d}V-\int_{\Gamma}\epsilon_{0}\epsilon_{\text{r}}\phi_{,i}N_{\phi}^{I}\hat{n}_{i}\text{d}\Gamma-\int_{V}N_{\phi}^{I}\frac{\mathcal{F}}{V_{\text{m}}}(z_{\text{V}_{\text{O}}^{..}}C_{\text{V}_{\text{O}}^{..}}+z_{\text{dop}}C_{\text{dop}})\text{d}V$
(59)

The linearization of the residuals can be written in the element-level linear
formulation, in which the terms in the tangent matrix are calculated by
$K_{\xi,\zeta}^{IJ}=\partial R_{\xi}^{I}/\partial\zeta^{J}$, where
${\xi,\zeta}$ indicates individual DOF. The non-zero terms in the tangent
matrix are presented in supplementary. The numerical implementation of the
proposed phase-field model is carried out in Multiphysics Object-Oriented
Simulation Environment (MOOSE) framework [47, 48].

##  5 Phase-field simulation results of Fe-doped SrTiO3

We present first in Section 5.1 one-dimensional phase-field simulation results
to demonstrate its capability of reproducing the formation of a SCL at
equilibrium state. Subsequently, we investigate the influence of solute drag
effects induced by the segregation of point defects in the core on the
formation of an asymmetric SCL under quasi-equilibrium state in Section 5.2.
Finally, we conduct two-dimensional (2D) phase-field simulations to study the
microstructure evolution process in STO polycrystalline material during
sintering and analyze the resulting grain size distribution in Section 5.3.

Note that in the paper we ignore the GB orientation dependency or other
anisotropy of the segregation energies $\Delta g_{\text{V}_{\text{O}}^{..}}$
and $\Delta g_{\text{dop}}$, for we intend to reveal the most fundamental
mechanisms of SCL formation and abnormal grain growth. Without loss of
generality, the proposed phase-field model can be extended accordingly to
include such anisotropy aspects.

###  5.1 GB space charge layers at equilibrium

To demonstrate the capability of the phase-field model at thermodynamically
equilibrium state, one-dimensional phase-field simulations are conducted and
compared with the sharp interface predictions [30]. The semi-finite simulation
domain has the size of 780 nm, with the core located at the left boundary. The
grain boundary thickness is prescribed to be 0.78 nm [30]. The domain is
discretized into a mesh consisting of 2000 elements. The simulation
temperature is set at 600 K, and the relative permittivity
$\epsilon_{\text{r}}$ of strontium titanate (STO) is specified as 160,
according to [30], which is calculated via $\epsilon_{\text{r}}=9\times
10^{4}/(T-35)$ as mentioned in section 3.1. The interface energy, molar volume
and interface mobility of STO are 0.6 $\text{J}\text{m}^{-2}$, 3.6$\times
10^{-5}$ $\text{m}^{3}\text{mol}^{-1}$ and 1$\times 10^{-7}$
$\text{m}^{3}\text{J}^{-1}\text{s}^{-1}$, respectively [38]. The relative
tolerance for convergence is set to $1\times 10^{-8}$.

####  5.1.1 Benchmark case 1: The MS Model

By fixing a uniform dopant concentration, the proposed phase-field model can
be used to simulate the MS case. The results are given in Fig. 2. The initial
concentrations of oxygen vacancies and dopant far from the core (denoted as
$c_{\text{dop},\text{b}}(\infty)$ and
$c_{\text{V}_{\text{O}}^{..},\text{b}}(\infty)$, respectively) are $7\times
10^{24}\text{ m}^{-3}$ and $3.5\times 10^{24}\text{ m}^{-3}$ [30]. The number
of available sites of oxygen vacancy in bulk phase and GB core are
$N_{\text{V}_{\text{O}}^{..},\text{b}}=5.12\times 10^{28}$ $\text{m}^{-3}$ and
$N_{\text{V}_{\text{O}}^{..},\text{c}}=3\times 10^{26}$ $\text{m}^{-3}$,
respectively[30, 49]. The grain boundary potential $\Phi_{0}$ is plotted as a
function of $\Delta g_{\text{V}_{\text{O}}^{..}}$ in Fig. 2(a). Due to the
right boundary being grounded, the grain boundary potential is directly
obtained at the interface position when $\eta=0.5$, as defined in the phase-
field method. For $\Delta g_{\text{V}_{\text{O}}^{..}}$ larger than -1.25 eV,
the phase-field results exhibit good agreement with analytical calculations.
As $\Delta g_{\text{V}_{\text{O}}^{..}}$ becomes more negative, particularly
at -1.5 eV and -2.0 eV, the available oxygen vacancy sites are completely
exhausted. Therefore, the concentrations of oxygen vacancies in the core are
close to $N_{\text{V}_{\text{O}}^{..},\text{c}}$. Although the phase-field
results demonstrate convergent behavior, a discrepancy between the phase-field
results and the analytical solution emerges. To explain the aforementioned
discrepancy, we plot the concentration of oxygen vacancies in the core as a
function of $\Delta g_{\text{V}_{\text{O}}^{..}}$ in Fig. 2(b). Similar to the
grain boundary potential, $c_{\text{V}_{\text{O}}^{..}}$ is also extracted at
$\eta=0.5$. We observe that the phase-field results coincide with the
analytical predictions, even when $\Delta g_{\text{V}_{\text{O}}^{..}}=-2.0$
eV. Next, we examine the concentration of oxygen vacancy in the bulk in terms
of $x$. Here, we plot the phase-field results for $\Delta
g_{\text{V}_{\text{O}}^{..}}$=-0.5 eV, -1.0 eV and -1.5 eV in Fig. 2 (c), (d)
and (e) in the range of 0 to 50 nm, which includes half of the interface
region. For the analytical solution, we shift it to the right by a distance of
$w_{\text{c}}/2$ to match the bulk phase in the phase-field results. When
${\Delta g_{\text{V}_{\text{O}}^{..}}}=-0.5$ eV,
$c_{\text{V}_{\text{O}}^{..},\text{b}}$ obtained by phase-field simulations
agrees well with analytical predictions as shown in Fig. 2(c1). However, in
Fig. 2(d1) and (e1), for $\Delta g_{\text{V}_{\text{O}}^{..}}=-1.5$ and $-2.0$
eV, the phase-field results produce higher
$c_{\text{V}_{\text{O}}^{..},\text{b}}$. This comparison indicates that less
oxygen vacancy is segregated in the core in the phase-field simulations when
choosing more negative $\Delta g_{\text{V}_{\text{O}}^{..}}$. Additionally, we
plot the area density of the net charge,
$Q_{\mathrm{net}}=\int_{-w_{\text{c}}/2}^{w_{\text{c}}/2}(z{{}_{\text{dop}}}c_{{}_{\text{dop},\text{c}}}+z_{{}_{\text{V}_{\text{O}}^{..}}}c_{{}_{\text{V}_{\text{O}}^{..},\text{c}}})\text{d}x$,
in the core in Fig. 2(f). Less net charge is reproduced via phase-field
simulation when $\Delta g_{\text{V}_{\text{O}}^{..}}$ is smaller. In short,
this discrepancy is attributed to the diffuse interface in the phase-field
simulations. Due to its finite interface thickness and smooth interface
characteristics, the core is less positively charged compared to analytical
solutions. Consequently, the electrostatic potential profile obtained via
phase-field simulations decays faster than analytical calculations as shown in
Fig. 2(c2), (d2) and (e2). Up to this point, we clearly observe a good
agreement between the oxygen vacancy concentrations in the GB core and bulk
phase as predicted by the phase-field results and the sharp interface
predictions. The different approaches to handling interfaces in the phase-
field model and the sharp interface method can result in discrepancies when
the available densities of point defects are fully occupied.

Figure 2: Comparison of phase-field (PF) simulation outcomes with those of the
sharp interface (SI) model predictions obtained via FEM and taken from Ref.
[30] in the MS model. Figure (a) illustrates the grain boundary potential as a
function of $\Delta g_{\text{V}_{\text{O}}^{..}}$. Figure (b) depicts the
relationship between the oxygen vacancy concentration and $\Delta
g_{\text{V}_{\text{O}}^{..}}$. Figures (c), (d), and (e) present the
concentration distribution of dopant and the electrostatic potential
distribution for different oxygen vacancy segregation energy. Figure (f)
demonstrates the net area density of charge in the core as a function of
$\Delta g_{\text{V}_{\text{O}}^{..}}$.

####  5.1.2 Benchmark case 2: The GC Model

By setting the electrochemical potentials of the acceptor dopant and of the
oxygen vacancies to be constant, the presented phase-field model is then
reduced to the GC model, in order to reproduce the equilibrium SCL formation
in the GC case. The number of available sites of dopant and oxygen vacancy in
bulk and GB core are $N_{\text{V}_{\text{O}}^{..},\text{c}}=1\times 10^{27}$
$\text{m}^{-3}$, $N_{\text{V}_{\text{O}}^{..},\text{b}}=5.12\times 10^{28}$
$\text{m}^{-3}$, $N_{\text{dop},\text{c}}=1.68\times 10^{27}$ $\text{m}^{-3}$
and $N_{\text{dop},\text{b}}=1.69\times 10^{28}$ $\text{m}^{-3}$, respectively
[30, 49]. The Debye lengths
$l_{\text{D}}=\sqrt{\frac{\epsilon_{0}\epsilon_{\text{r}}k_{\text{B}}T}{2e^{2}c_{\text{dop},\text{b}}(\infty)}}$
in both MS and GC models are 5.71 nm. In Fig. 3(a), the grain boundary
potentials obtained via the phase-field simulations are compared with the
sharp interface predictions. $\Delta g_{\text{V}_{\text{O}}^{..}}$ is fixed at
-1.0 eV here.

When $\Delta g_{\text{dop}}$ is positive, dopant cations are rejected from the
core, thereby causing a more pronounced grain boundary potential. Conversely,
as $\Delta g_{\text{dop}}$ becomes negative, dopant segregation to the core
occurs, reducing the net charge and consequently decreasing the grain boundary
potential. Moreover, as $\Delta g_{\text{dop}}$ becomes sufficiently negative,
the available number sites for dopant segregating within the core become
scarce. Therefore, the grain boundary potential profile is flat when $\Delta
g_{\text{dop}}$ is positive or sufficiently negative. The comparison indicates
the phase-field results agree with the predictions of analytical solution. The
discrepancy is smaller than 7%. The concentrations of vacancies and acceptor
dopants in the core as a function of $\Delta g_{\text{dop}}$ are plotted in
Fig. 3(b).

In Fig. 3(c), (d), and (e), we present the concentration profiles of oxygen
vacancies and dopant, as well as the distribution of electrostatic potential
in the bulk phase for $\Delta g_{\text{dop}}$ = 0 eV, -0.125 eV and -0.375 eV,
respectively. Due to symmetry, we only depict half of the SCL in one bulk
phase.

At $\Delta g_{\text{dop}}$ = 0 eV, phase-field results predict the formation
of a dopant accumulation zone near the core. The dopant concentration in the
core can be much lower than in the accumulation zone, leading to the formation
of a nonmonotonic curve as shown in Fig. 3(c1). Contrarily, in Virkant’s
findings [38], the expected accumulation zone at $\Delta g_{\text{dop}}$=0 eV
is replaced by a depletion zone. As $\Delta g_{\text{dop}}$ decreases to
-0.125 eV, dopant concentration in the core increases significantly, causing
the formation of monotonic concentration profile. Further decreasing $\Delta
g_{\text{dop}}$ to -0.375 eV, the dopant concentration in the core increases
to $N_{\text{dop},\text{b}}$. It is evident that the accumulation zone becomes
lower than 0 eV, more exactly -0.125 eV. Additionally, near the core, the
concentration profile of dopant exhibits a distinct kink, indicating a notable
change in the gradient.

The oxygen vacancy concentration profiles of are plotted in Fig. 3(c2), (d2)
and (e2). Depletion appears for all $\Delta g_{\text{dop}}$. Notably, as
$\Delta g_{\text{dop}}$ decreases, there is an increase in the concentrations
within these depletion zones. In Fig. 3(c3), (d3) and (e3), we present a
comparison of the electrostatic potential across the SCL with analytical
solutions. The phase-field simulation results agree well with these analytical
predictions. Decreasing $\Delta g_{\text{dop}}$ causes lower electrostatic
potential distribution near the core. Moreover, we also present the net charge
area density $Q_{\mathrm{net}}$ in Fig. 3(f). The trend of the curve obtained
via phase-field results agrees with the analytical predictions. Analogous to
the MS case, phase-field results are smaller than analytical predictions.

Figure 3: Comparison of phase-field (PF) simulation results with the sharp
interface (SI) results via FEM method for the GC case. Figures (a) and (b)
showcase the variations in grain boundary potential, dopant concentration, and
oxygen vacancy concentration in the core as functions of $\Delta
g_{\text{dop}}$. Figures (c), (d), and (e) detail the spatial distributions of
dopant and oxygen vacancy concentrations in the bulk phase and the
electrostatic potential. Figure (f) depicts the area density of net charge in
the core, also as a function of $\Delta g_{\text{dop}}$. $\Delta
g_{\text{V}_{\text{O}}^{..}}$ was kept at -1.0 eV for all calculations.

###  5.2 GB space charge layers at quasi-equilibrium

In this section, we investigate systematically the solute drag effect caused
by an asymmetric SCL at a quasi-statically moving GB in a STO bicrystal. The
solute drag phenomenon arises when point defects, segregated within a moving
core in a quasi-equilibrium state, lead to additional free energy dissipation
to maintain motion. At low core velocities, these defects can keep pace with
the core, significantly increasing the necessary driving force. However,
beyond a critical velocity, the defects lag behind, remaining in the bulk
phases and causing the breakdown of symmetric concentration and electrostatic
potential profiles across the SCL.

To investigate the formation of asymmetric SCL due to solute drag effects
under quasi-equilibrium conditions, a constant core velocity
(${v}_{\text{c}}$) along the simulation domain is assigned as a parameter in
the simulations. The velocity of the core ranges from 0.01
$\text{nm}\,\text{s}^{-1}$ to 100 $\text{nm}\,\text{s}^{-1}$ to study its
influence on the SCL. The defect chemistry parameters
$N_{\text{V}_{\text{O}}^{..},\text{c}}$,
$N_{\text{V}_{\text{O}}^{..},\text{b}}$, $N_{\text{dop},\text{c}}$ and
$N_{\text{dop},\text{b}}$ in this case are $1\times 10^{27}$ $\text{m}^{-3}$,
$5.12\times 10^{28}$ $\text{m}^{-3}$, $1.68\times 10^{27}$ $\text{m}^{-3}$ and
$1.69\times 10^{28}$ $\text{m}^{-3}$, respectively. Simulations are conducted
at 1623 K, with a relative permittivity ($\epsilon_{\text{r}}$) of STO set to
56, based on [38]. The diffusivities of oxygen vacancies and dopant in the
bulk phase are specified as $1.35\times 10^{-9}$ $\text{m}^{2}\text{s}^{-1}$
and $1.3\times 10^{-11}$ $\text{m}^{2}\text{s}^{-1}$, respectively [38]. The
simulation domain and boundary conditions remain consistent with the previous
case.

In alignment with Vikrant [38], the total driving force $F_{\text{T}}$ for the
core motion is given by

$F_{\text{T}}=\int_{-\infty}^{\infty}\frac{v_{c}}{L}\left(\frac{\partial\eta}{\partial
x}\right)^{2}\text{d}x+\int_{-\infty}^{\infty}{\mu_{\text{V}_{\text{O}}^{..}}}\frac{\partial
C_{\text{V}_{\text{O}}^{..}}}{\partial
x}\text{d}x+\int_{-\infty}^{\infty}{\mu_{\text{dop}}}\frac{\partial
C_{\text{dop}}}{\partial x}\text{d}x.$ (60)

The second and third terms in Eq. 60 denote as the electrochemical drag force,
which accounts for the electrostatic and chemical contributions. The
integration in Eq. 60 are solved numerically.

In Fig. 4, we illustrate the concentration profiles of oxygen vacancies and
dopant, as well as the charge density and electrostatic potential as functions
of $x$. From Fig. 4(a) to Fig. 4(e), we progressively increase the velocity of
the core from 0 $\text{nm}\,\text{s}^{-1}$ to 100 $\text{nm}\,\text{s}^{-1}$.
The direction of the core motion is from left to right along $x$-axis. In the
following context, the left side of the core is named as the growing grain,
while the right side is called the shrinking grain.

In Fig. 4(a), the symmetric SCL formation is presented at $v_{c}=0$
$\text{nm}\,\text{s}^{-1}$. However, the symmetric pattern begins to
deteriorate even at very low velocity, e.g. 0.01 $\text{nm}\,\text{s}^{-1}$,
as shown in Fig. 4(b1), (b2) and (b3). Compared with the concentration
profiles at an equilibrium state, the accumulation zone of dopant in the
growing grain is reduced, while that of the shrinking grain increases. This
asymmetry arises, since during the core motion, dopant cations in the
shrinking grain are captured by core due to their low diffusivity. Driving the
motion of extra dopant cations within the core demands extra free energy
dissipation. Therefore, the dopant cations flowing from the core region into
the growing grain lead to the enhancement of the accumulation zone in the
growing grain. However, the concentration of dopant in the core decreases from
$7.8\times 10^{25}$ $\text{m}^{-3}$ ($v_{c}=0$ in Fig. 4(a1)) to $6.8\times
10^{28}$ $\text{m}^{-3}$ at $v_{c}=0.01$ $\text{nm}\,\text{s}^{-1}$. In the
growing grain far from the core, dopant concentration reaches $7.2\times
10^{25}$ $\text{m}^{-3}$, which is higher than equilibrium concentration, i.e.
$c_{\text{dop},\text{b}}(\infty)$, as shown in Fig. 4(a1). In Fig. 4(b2), we
plot the calculated charge density
$\rho=e(z_{\text{dop}}c_{\text{dop}}+z_{\text{V}_{\text{O}}^{..}}c_{\text{V}_{\text{O}}^{..}})$.
The charge density in the core decreases. Compared with the equilibrium case,
we also observe the formation of less (more) negative charged zone in the
shrinking (growing) grain, respectively. In Fig. 4(b3), the grain boundary
potential increases from 0.248 V ($v_{c}=0$ $\text{nm}\,\text{s}^{-1}$) to
0.265 V at $v_{c}=0.01$ $\text{nm}\,\text{s}^{-1}$. Moreover, an electrostatic
potential different of 0.011 V develops between the shrinking and growing
grains.

When we increase the velocity of the core to 0.5 $\text{nm}\,\text{s}^{-1}$,
dopant cations are not capable of catching the movement of core, as depicted
in Fig. 4(c1). Consequently, the dopant concentration within the core
decreases substantially, resulting in a smaller depletion zone ahead of the
core and a reduced accumulation zone behind it. Similarly, the concentration
of oxygen vacancies within the core also decreases significantly.
Correspondingly, as shown in Fig. 4(c2), the charge density in the core drops
to $8.6\times 10^{-7}$ $\text{C}\text{m}^{-3}$, which contrasts with the
increase in charge density in the core observed in Virkant’s work [38] at
intermediate velocities. We attribute the difference to the defect chemistry
incorporated in the present phase field model.

When the core velocity increases to 2 $\text{nm}\,\text{s}^{-1}$, the
segregation of dopant in the core becomes negligible, as shown in Fig. 4(d1).
The concentration of dopant on both the left and right sides of the bulk phase
returns to the equilibrium value, i.e.
$c_{\text{V}_{\text{O}}^{..},\text{b}}(\infty)$. In Fig. 4(d2), due to the
continuous decrease in charge density within the core, the negatively charged
zone is diminished. Furthermore, the grain boundary potential increases to
0.446 V and the electrostatic potential difference between the left and right
sides reaches 0.122 V.

At a high core velocity, such as 100 $\text{nm}\,\text{s}^{-1}$, the dopant
concentration profile becomes completely flat, while the oxygen vacancy
concentration profile exhibits strong asymmetry, as depicted in Fig. 4(e1).
The concentration of oxygen vacancies within the core decreases significantly
to $2.46\times 10^{26}$ $\text{m}^{-3}$. In the shrinking grain, the depletion
zone of oxygen vacancies remains sharp, about 2.5 nm, but it extends to
approximately 20 nm in the growing grain. Consequently, as shown in Fig.
4(e2), the negatively charged zone in the shrinking grain is much smaller, and
the negative charge region in the growing grain is negligible. In Fig. 4(e3),
the electrostatic potential is illustrated. The grain boundary potential
increases to 1.02 V and decays very slowly in the growing grain. In addition,
the electrostatic potential difference becomes more pronounced, reaching 0.995
V.

In summary, at small core velocity, the asymmetric formation of SCL is mainly
caused by dopant coming apart from the core. The charge density in the
negatively charged SCL and in the positively charged core both decrease. As a
result, the electrostatic potential in the shrinking grain increases to a
positive value. At high core velocity, dopant cations do not segregate to the
core. The oxygen vacancy concentration profile is strongly asymmetric. This
asymmetry causes an extended diffusion layer in the growing grain and elevates
the electrostatic potentials in both the core and the shrinking grain.

Figure 4: Asymmetric SCL formation in a STO bicrystal system with
$c_{\text{dop},\text{b}}(\infty)=7.0\times 10^{25}$ $\text{m}^{-3}$, $\Delta
g_{\text{V}_{\text{O}}^{..}}=-1.5$ eV and $\Delta g_{\text{dop}}=-0.5$ eV at
quasi-equilibrium state under different constant core velocity (a)
$v_{\text{c}}=0$ $\text{nm}\,\text{s}^{-1}$, (b) $v_{\text{c}}=0.01$
$\text{nm}\,\text{s}^{-1}$, (c) $v_{\text{c}}=0.5$ $\text{nm}\,\text{s}^{-1}$,
(d) $v_{\text{c}}=2$ $\text{nm}\,\text{s}^{-1}$, (e) $v_{\text{c}}=100$
$\text{nm}\,\text{s}^{-1}$. Each row has three sub-figures, depicting the
distribution of dopant and oxygen vacancy concentrations, total charge density
and electrostatic potential, respectively.

In order to further investigate solute drag effects and identify the critical
core velocity at which a velocity jump occurs during sintering process, we
present the core velocity as a function of the total driving force calculated
via Eq. 60 in Fig. 5. The segregation energy difference, $\Delta
g_{\text{dop}}$ and $\Delta g_{\text{V}_{\text{O}}^{..}}$, allow direct tuning
of the concentrations of oxygen vacancy and dopant in the core, thereby
impacting the critical total driving force required for a specific core
velocity. In Fig. 5(a), we hold $\Delta g_{\text{dop}}$ and
$c_{\text{dop},\text{b}}(\infty)$ constant at -0.5 eV and $7\times 10^{25}$
$\text{m}^{-3}$ and progressively decrease $\Delta
g_{\text{V}_{\text{O}}^{..}}$ from -0.5 eV to -1.75 eV. During the grain
shrinking process, the driving force resulting from the GB curvature becomes
increasingly significant. Consequently, the velocity of the GB core evolves
from a low value to a high value. At low core velocities, a more negative
$\Delta g_{\text{V}_{\text{O}}^{..}}$ induces higher concentration of oxygen
vacancy in the core consequently demanding a larger driving force to maintain
the same velocity. As the core velocity increases to 1
$\text{nm}\,\text{s}^{-1}$, the dopant cations are unable to keep pace with
the core movement, resulting in a nearly flat dopant concentration within the
core. Hence, this leads to the same linear relationship between the total
driving force and the core velocity for different $\Delta
g_{\text{V}_{\text{O}}^{..}}$ at high velocity. On the other hand, for fixed
$\Delta g_{\text{V}_{\text{O}}^{..}}$, such as $\Delta
g_{\text{V}_{\text{O}}^{..}}$=-1.75 eV, we clearly observe the velocity jumps
from 0.065 $\text{nm}\,\text{s}^{-1}$ to approximate 1.5
$\text{nm}\,\text{s}^{-1}$ at $1.45\times 10^{7}$ $\text{N}\text{m}^{-2}$. As
we increase $\Delta g_{\text{V}_{\text{O}}^{..}}$ to -1.0 eV, we hardly
observe the velocity jump. When $\Delta g_{\text{V}_{\text{O}}^{..}}$ = -0.5
eV, both the concentrations of dopant and oxygen vacancy in the core decrease
significantly. The solute drag effect therefore becomes negligible. Fig. 5(b)
shows the influence of $\Delta g_{\text{dop}}$ on solute drag effects when
$\Delta g_{\text{V}_{\text{O}}^{..}}$ and $c_{\text{dop},\text{b}}$ are held
constant at -1.5 eV and $7\times 10^{25}$ $\text{m}^{-3}$. In contrast to Fig.
5(a), we can evidently observe velocity jumps for all $\Delta g_{\text{dop}}$.
As we decrease $\Delta g_{\text{dop}}$ to less negative value, such as from
-0.2 eV to -0.3 eV, the enhancement of the solute drag effect is not
significant. In Fig. 5(c), we explore the impact of the dopant concentration
in the bulk far from the core, denoted as $c_{\text{dop},\text{b}}(\infty)$.
We fix $\Delta g_{\text{V}_{\text{O}}^{..}}$ and $\Delta g_{\text{dop}}$ at
-1.5 eV and -0.5 eV. As $c_{\text{dop},\text{b}}(\infty)$ increases, a larger
driving force is necessary to counterbalance the solute drag effect to achieve
the same core velocity. Additionally, a lower critical velocity is observed
when the velocity jump occurs at smaller values of
$c_{\text{dop},\text{b}}(\infty)$. Specifically, the critical velocity
increases from 0.065 $\text{nm}\,\text{s}^{-1}$ to 0.095
$\text{nm}\,\text{s}^{-1}$, as $c_{\text{dop},\text{b}}(\infty)$ is raised
from $4\times 10^{25}$ $\text{m}^{-3}$ to $2\times 10^{26}$ $\text{m}^{-3}$.

In the sintering process, the driving forces are initially large and decrease
as the grain size increases. In other words, the GB core velocity starts on
the upper right side of the S-shaped curves in Fig. 5 at the initial stage of
the grain growth. As the grain size increases, the GB core velocity decreases
to the critical velocity. At this point, a transition from high velocity to
low velocity occurs at the critical value. In Fig. 5(a), we observe GB core
velocity jumps when $\Delta g_{\text{V}_{\text{O}}^{..}}=$ -1.5 eV and -1.75
eV. The critical velocities for different $\Delta
g_{\text{V}_{\text{O}}^{..}}$ are approximately 0.35
$\text{nm}\,\text{s}^{-1}$. Additionally, the GB core velocity jumps to a
lower value when $\Delta g_{\text{V}_{\text{O}}^{..}}$ is more negative. In
Fig. 5(b), the critical GB core velocity increases from 0.15
$\text{nm}\,\text{s}^{-1}$ to 0.65 $\text{nm}\,\text{s}^{-1}$ as $\Delta
g_{\text{dop}}$ decreases from -0.3 eV to -0.7 eV. In Fig. 5(c), a larger
$c_{\text{dop},\text{b}}(\infty)$ leads to a higher critical velocity.

Figure 5: The total driving force on a GB with a constant velocity with (a)
$\Delta g_{\text{dop}}=-0.5$ eV, $c_{\text{dop},\text{b}}(\infty)=7\times
10^{25}$ $\text{m}^{-3}$ and different $\Delta g_{\text{V}_{\text{O}}^{..}}$,
(b) $\Delta g_{\text{V}_{\text{O}}^{..}}=-1.5$ eV,
$c_{\text{dop},\text{b}}(\infty)=7\times 10^{25}$ $\text{m}^{-3}$ and
different $\Delta g_{\text{dop}}$, (c) $\Delta
g_{\text{V}_{\text{O}}^{..}}=-1.5$ eV, $\Delta g_{\text{dop}}=-0.5$ eV and
different $c_{\text{dop},\text{b}}(\infty)$. The pentagram symbols in
subfigure b correspond to the different GB core velocities studied in Fig. 4.
The green, blue and saffron pentagram symbols indicate $v_{\text{c}}$=2
$\text{nm}\,\text{s}^{-1}$, 0.5 $\text{nm}\,\text{s}^{-1}$ and 0.01
$\text{nm}\,\text{s}^{-1}$, respectively. Due to the range limit of y-axis,
$v_{\text{c}}$=0 $\text{nm}\,\text{s}^{-1}$ (Fig. 4 a) and 100
$\text{nm}\,\text{s}^{-1}$ (Fig. 4 e) are not marked.

###  5.3 Abnormal grain growth

After systematically benchmarks and knowledge gained from quasi-static single
GB cases, we investigate different abnormal grain growth patterns of materials
with different defect segregation energies $\Delta g_{\text{dop}}$, $\Delta
g_{\text{V}_{\text{O}}^{..}}$ and bulk defect concentrations such as
$c_{\text{dop},\text{b}}$. The segregation of point defects in the core plays
a significant role in influencing grain growth kinetics. To investigate the
grain size distribution of STO polycrystalline material after sintering, we
conduct 2D phase-field simulations. Note that three-dimensional phase-field
simulations of the proposed model are achievable but require significantly
more computational resources, which will be considered in the near future. In
fact, to reveal the fundamental mechanisms, 2D simulations should be
sufficient and allow for a comprehensive parameter study given the limited
computational resources.

The simulation box size is set to 390 nm $\times$ 390 nm with a grid
resolution of 1000 $\times$ 1000 points. The polycrystalline microstructure is
initialized using a Voronoi tessellation method, where each grain was assigned
a unique grain index. We started with an initial number of 1000 grains, and
the random seed used for generating the polycrystalline microstructure is set
to 42. To optimize memory usage during simulations, we also utilize grain
tracker algorithm [34]. In the simulations, we set
$N_{\text{V}_{\text{O}}^{..},\text{c}}$,
$N_{\text{V}_{\text{O}}^{..},\text{b}}$, $N_{\text{dop},\text{c}}$ and
$N_{\text{dop},\text{b}}$ as $1\times 10^{27}$ $\text{m}^{-3}$, $5.12\times
10^{28}\text{m}^{-3}$, $1.68\times 10^{27}$ $\text{m}^{-3}$ and $1.69\times
10^{28}$ $\text{m}^{-3}$. The simulations are conducted at a temperature of
1623 K, where the non-Arrhenius behavior of the grain growth is found [28].
Periodic boundary conditions are applied to all edges of the simulation domain
for the phase-fields and concentration fields. For the electrostatic potential
field, the bottom-left corner is grounded, while the other edges are set with
the periodic boundary conditions. Furthermore, all the parameters used in the
simulation are normalized. The lattice constant $a$ = $3.9\times 10^{-10}$ m
is used for the length normalization. Other normalization factors include e.g.
$f_{0}=\kappa/a^{2}$=$2.3\times 10^{9}$ $\text{J}\text{m}^{-3}$ for the energy
density, $t_{0}=1/(Lf_{0})=0.0043$ s for the interface mobility, and
$\phi_{0}=f_{0}V_{\text{m}}/\mathcal{F}$= 0.858 V for the electrostatic
potential.

We present first the microstructure evolution of the grain growth process
under specific conditions $\Delta g_{\text{V}_{\text{O}}^{..}}$ = -1.5 eV,
$\Delta g_{\text{dop}}$=-0.4 eV and $c_{\text{dop},\text{b}}=7\times 10^{25}$
$\text{m}^{-3}$ when the sintering time is
$t/\text{h}=\\{0.067,0.33,0.67,1.0\\}$. The average area of grain at time
$t=0$ is <S>, i.e., 152.1 $\text{nm}^{2}$. When $t$=0.067 h, 0.33 h and 0.67
h, the grain size increases while the area fraction of small grain decreases.
The grain size distribution starts to transition from a mono-modal to a
bimodal pattern. When $t$=1 h, the largest grain area is about 32 times larger
than <S>. The bimodal grain size distribution is clearly visible.

Figure 6: The microstructure evolution of grain growth process at (a) 0.067 h,
(b) 0.33 h, (c) 0.67 h and (d) 1 h, with defect chemistry parameters set at
$\Delta g_{\text{V}_{\text{O}}^{..}}$ = -1.5 eV, $\Delta g_{\text{dop}}$=-0.4
eV and $c_{\text{dop},\text{b}}=7\times 10^{25}$ $\text{m}^{-3}$. These
parameters are consistent with those used in Fig. 7(e). By 1h, a bimodal grain
size distribution becomes apparent.

####  5.3.1 Influence of segregation energies and bulk defect concentration

Subsequently, we investigate the influence of defect chemistry parameters on
the microstructure evolution during sintering processes. The microstructure of
grains is presented in Fig. 7(a)-(g) after 1 hour of sintering at 1623 K. In
Fig. 7(a), (b) and (c), we increase $\Delta g_{\text{V}_{\text{O}}^{..}}$ from
-1.5 eV to -0.5 eV when $\Delta g_{\text{dop}}$=-0.5 eV and
$c_{\text{dop},\text{b}}=7\times 10^{25}$ $\text{m}^{-3}$. Bimodal grain size
distribution is clearly observed in Fig. 7(a2), (b2) and (c2). As $\Delta
g_{\text{V}_{\text{O}}^{..}}$ becomes more negative, the peak in the range of
$0<\frac{S}{<S>}<4$ induced by smaller grains becomes more pronounced. As
predicted in Fig. 5(a), solute drag effect is almost negligible for $\Delta
g_{\text{V}_{\text{O}}^{..}}$=-0.5 eV. Nonetheless, the formation of solute
cloud imprinted by vanished GBs is still detectable in Fig. 7(c3) even though
the solute drag effect is very weak. Concentration of oxygen vacancy is
detailed in Fig. 7(a4), (b4) and (c4). Segregation is more significant for
$\Delta g_{\text{V}_{\text{O}}^{..}}=-1.5$ eV. Compared to the dopant, solute
cloud of oxygen vacancies are less distinct. In view of the higher diffusivity
of oxygen vacancy, the contribution of oxygen vacancy to solute drag force is
also weak. Solute clouds of oxygen vacancies are only observed near the
junctions where grains shrink rapidly with a small radius. The electrostatic
potential distributions are depicted in Fig. 7(a5), (b5) and (c5). The
negative potential bands can be also observed near the GB cores. As all the
considered cases show, the electrostatic potentials in small shrinking grains
are significantly higher than in the surrounding growing grains, leading to a
stronger electrostatic contribution to the solute drag effects. Consequently,
the shrinking speed of small grains is reduced, which induces pinning effects
on the surrounding growing grains. Note that the heterogeneous electrostatic
potential distribution arises from the velocity difference of individual GB,
in accordance to the observation in the quasi-static cases in Fig. 4, rather
than any anisotropy which is absent in the current simulations.

In Fig. 7(a), (d) and (e), we keep $\Delta g_{\text{V}_{\text{O}}^{..}}$ and
$c_{\text{dop},\text{b}}$ at -1.5 eV and $7\times 10^{25}$ $\text{m}^{-3}$ and
test different $\Delta g_{\text{dop}}$ values. Whereas, in Fig. 7(a), (f) and
(g), $\Delta g_{\text{V}_{\text{O}}^{..}}$ and $\Delta g_{\text{dop}}$ are
constant while $c_{\text{dop},\text{b}}$ increases from $2\times 10^{25}$
$\text{m}^{-3}$ to $2\times 10^{26}$ $\text{m}^{-3}$. For all simulations, the
solute drag effects from dopant are pronounced. The dopant clouds imprinted by
vanished GBs are clearly observed behind GB cores. In Fig. 7(e2) and (g2), the
bimodal grain size distributions are presented. However, when $\Delta
g_{\text{dop}}=-0.6$ eV in Fig. 7(d2) and $c_{\text{dop},\text{b}}=2\times
10^{26}$ $\text{m}^{-3}$ in Fig. 7(f2), the unimodal grain size distribution
exhibits, due to the more pronounced probability of small grain. The
simulation results in Fig. 7(f) demonstrate the high grain number, where the
large grains are surrounded by small grains.

Figure 7: Phase-field simulation results of abnormal grain growth in Fe-doped
STO after 1 h sintering under the influence of different defect chemistry
parameters.

####  5.3.2 Grain boundary potential distribution

In Fig. 8, we present the statistics of the average grain boundary potential
in the GB core, denoted as $\bar{\phi}_{\text{core}}$. The value of
$\bar{\phi}_{\text{core}}$ is calculated as an average of the electrostatic
potential along the GB core indicated by $\eta_{i}=0.5$ over individal GB
length. Fig. 8(a) illustrates the distribution of $\bar{\phi}_{\text{core}}$
for $\Delta g_{\text{V}_{\text{O}}^{..}}=-0.5$ eV, $\Delta
g_{\text{dop}}=-0.5$ eV, $c_{\text{dop}}(\infty)=7\times 10^{25}$
$\text{m}^{-3}$ (refer to Fig. 7(c) for the associated phase-field simulation
result). When the solute drag effect is negligible, $\bar{\phi}_{\text{core}}$
is smaller than 0.05 V and approximately follows the normal distribution. In
Fig. 8(b), the distribution of $\bar{\phi}_{\text{core}}$ is presented for
$\Delta g_{\text{V}_{\text{O}}^{..}}=-1.5$ eV, $\Delta g_{\text{dop}}=-0.5$
eV, $c_{\text{dop}}(\infty)=7\times 10^{25}$ $\text{m}^{-3}$ (refer to Fig.
7(a)). Here, $\bar{\phi}_{\text{core}}$ significantly increases. In addition,
the distribution of $\bar{\phi}_{\text{core}}$ changes from mono-modal to
bimodal, with one peak at 0.31 V and another at 0.35 V. In Fig. 8(c), the
parameters $\Delta g_{\text{V}_{\text{O}}^{..}}=-1.5$ eV, $\Delta
g_{\text{dop}}=-0.5$ eV, $c_{\text{dop}}(\infty)=2\times 10^{26}$
$\text{m}^{-3}$ are chosen (refer to Fig. 7(f)). In this case, the solute drag
effect is the most pronounced in this study. Not only
$\bar{\phi}_{\text{core}}$ increases but also the bimodal distribution of
$\bar{\phi}_{\text{core}}$ becomes more distinct.

Figure 8: Probability of the average grain boundary potential in the GB cores
with different defect chemistry parameters (a)$\Delta
g_{\text{V}_{\text{O}}^{..}}=-0.5$ eV, $\Delta g_{\text{dop}}=-0.5$ eV,
$c_{\text{dop}}(\infty)=7\times 10^{25}$ $\text{m}^{-3}$ (refer to Fig. 7(c)),
(b) $\Delta g_{\text{V}_{\text{O}}^{..}}=-1.5$ eV, $\Delta
g_{\text{dop}}=-0.5$ eV, $c_{\text{dop}}(\infty)=7\times 10^{25}$
$\text{m}^{-3}$ (refer to Fig. 7(a)), (c) $\Delta
g_{\text{V}_{\text{O}}^{..}}=-1.5$ eV, $\Delta g_{\text{dop}}=-0.5$ eV,
$c_{\text{dop}}(\infty)=2\times 10^{26}$ $\text{m}^{-3}$ (refer to Fig. 7(f)).

In addition, two distinct types of GBs can be observed Fig. 9. The first type,
characterized by a lower moving velocity, exhibits nearly symmetric
concentration profiles across the GB core with a significantly higher
segregation concentration of dopant. However, within a polycrystalline system,
the accumulation zone of dopant at the space charge region transitions into a
depletion zone. Conversely, the second type of GB core, which moves at a
higher velocity, shows lower dopant concentrations due to the break of the
symmetric SCL. Similarly, the oxygen vacancy distribution in this second GB
core is also asymmetric, as illustrated in Fig. 7(e). Fig. 7(f) presents the
electrostatic potential distributions across these two types of GB cores,
where the second GB core displays a larger grain boundary potential and a more
pronounced electrostatic potential difference.

Figure 9: Two types of GBs are observed in the rectangular region highlighted
in Fig. 7(f). The distributions of dopant concentration, oxygen vacancy
concentration, and electrostatic potential in this region are shown in (a),
(b), and (c), respectively. Particularly, the dopant concentration profile
within the grain interior in (a) exhibits strong imprints of vanished GBs
during sintering. Two types of GB cores are distinguished by dashed yellow
lines in (a). The first type, characterized by a low moving velocity, has a
higher segregation concentration of dopants. In contrast, the second type,
with a higher core velocity, shows a lower and asymmetric dopant concentration
across the GB core. Furthermore, the dopant concentration, oxygen vacancy
concentration, and electrostatic potential along the two types of GB cores are
plotted in (d), (e), and (f), respectively.

In order to further evidence the two types of GB core, 14 different GBs are
selected from the rectangular region in Fig. 7 (f). The relations between
$v_{\text{c}}$ and $F_{\text{T}}$ as well as $\bar{\phi}_{\text{core}}$ and
$F_{\text{T}}$ for each GB core are represented in Fig. 10(b). Here,
$F_{\text{T}}$ is the total driving force, numerically calculated via Eq. 60.
The solid lines indicate the phase-field simulation results in the bicrystal
case (refer to the purple line in Fig. 5(c)), while the symbols represent
different GBs in the polycrystalline case. In Fig. 10(b), the relationships of
$F_{\text{T}}$ vs. $v_{\text{c}}$ and $F_{\text{T}}$ vs.
$\bar{\phi}_{\text{core}}$ in the polycrystalline case almost follow those in
the bicrystal case. There is also a velocity jump and an electrostatic
potential jump when the total driving force is approximately $2\times 10^{7}$
$\text{N}/\text{m}^{2}$. In the high GB velocity region, GB6, GB7 and GB8,
with less dopant segregation, have much higher electrostatic potentials than
the others, aligning well with the predictions in the bicrystal case. In the
low GB velocity region, GB3, GB13 and GB14 also agree well with the results
from the bicrystal quasi-static case. Whereas, GBs with higher velocities,
such as GB2, GB4, GB11 and GB12, do not fit well into the S curve of a
bicrystal case. This discrepancy can be attributed to multiple reasons. For
instance, these GBs have higher GB curvatures and can interact with other GBs,
which go beyond the bicystal model and can only be resolved by the grain
growth simulations presented here.

Figure 10: GB core velocity $v_{\text{c}}$ and average grain boundary
potential $\bar{\phi}_{\text{core}}$ as functions of the total driving force
$F_{\text{T}}$. The 14 different GBs, shown in Fig. 9(a), are chosen from the
rectangular region in Fig. 7(f) in polycrystalline simulations. The relations
between $v_{\text{c}}$ and $F_{\text{T}}$ as well as
$\bar{\phi}_{\text{core}}$ and $F_{\text{T}}$ are obtained from phase-field
simulations in bicrystal and polycrystalline Fe-doped STO. Solid lines
indicate the simulation results for the bicrystal case, while dashed lines are
the extrapolations of the simulation results in bicystal case at very low GB
velocity regime. Other symbols represent 14 different GBs chosen in the
polycrystalline case.

##  6 Conclusion

In this work, we develop a phase-field model to investigate SCL formation
during grain growth in oxide electroceramics based firmly on principles of the
defect chemistry. The main results of the paper include:

  * 1.

The proposed phase-field model reproduces a series of benchmark results at
equilibrium of Fe-doped STO bicrystal. Simulation results, including the
symmetric SCLs, grain boundary potentials, concentrations of oxygen vacancies
and acceptor dopants in the GB core and the bulk phase as well as the
electrostatic distributions, show good agreement with the analytical
predictions in the MS and GC models.

  * 2.

Simulations are carried out comprehensively for a GB core moving with constant
velocities under quasi-equilibrium state. At low core velocity, the
accumulation zone in the shrinking grain decreases, while it increases in the
growing grain, leading to the asymmetric SCL formation. With increasing GB
core velocity, the dopant segregation diminishes. However, the movement of
oxygen vacancy within the GB core only breaks at extremely high velocity due
to the high diffusivity of oxygen vacancy. In addition, the electrostatic
potential in the shrinking grain becomes higher than that in the growing
grain.

  * 3.

Abnormal grain growth is observed even without distinguishing any grain
orientation or considering any anisotropic GB mobility in the phase-field
simulations at 1623K. It indicates that the solute drag effect alone can lead
to bimodal grain size distribution. As the formation energy differences of
oxygen vacancies and acceptor dopants ($\Delta g_{\text{V}_{\text{O}}^{..}}$
and $\Delta g_{\text{dop}}$) become more negative, the grain growth velocity
decreases. This trend is also observed with increasing dopant concentration
($c_{\text{dop},\text{b}}$). Notably, the dopant concentration profile within
the grain interior distinctly reflects the significant imprint of the grain
boundaries that disappeared during sintering.

  * 4.

Simulation results reveal two distinct types of GB cores. The first type, with
nearly symmetric concentration profiles at the GB, exhibits stronger solute
drag effects, resulting in a much lower moving velocity and a tendency to
remain pinned in position. The second type of GB core shows asymmetric
concentration profiles and a higher moving velocity, contributing
predominantly to grain coarsening.

  * 5.

Abnormal grain growth resulting from solute drag effects leads to a large
variation of grain boundary potentials. Typically, smaller grains exhibit
higher grain boundary potentials at the later stages of grain growth and form
more blocking GBs in term of ionic conductivity. When considering the total
conductivity of an electrolyte material, the electrical current needs to find
a ”detour” around these blocking GBs by passing through the surrounding larger
grains, with less blocking GBs. It implies that to increase the conductivity
one may also optimize the microstructure so that blocking GBs can be
circumvented more easily. In contrast, if blocking GBs are desired, such as in
capacitors, the grain size should be kept small to prevent the development of
detour paths.

The proposed phase-field model offers naturally a series of promising
extensions. We can e.g. incorporate additional species, such as neutral
dopants, electrons, and electron holes, into the current phase-field model to
quantitatively investigate the formation of the SCL. Furthermore, GB
orientation plays a crucial role in SCL formation, and both the interface
energy and the defect chemistry parameters $\Delta
g_{\text{V}_{\text{O}}^{..}}$ and
$\omega_{c}N_{\text{V}_{\text{O}}^{..},\text{c}}$ depend on GB orientation
[49]. These complexities can be included in the phase-field model in the near
future. Moreover, we assume in this paper only the $\text{Fe}^{3+}$ and ignore
$\text{Fe}^{4+}$, which is reasonable for high temperature scenarios. But, the
situation regarding solute drag effect changes at low temperature when a
significant amount of $\text{Fe}^{4+}$ is present. The $\text{Fe}^{3+}$ ions
in the GB core can not easily follow the movement of the GB core due to their
low diffusivity. In contrast, the rapid movement of electron holes can adjust
the local $\text{Fe}^{3+}$/$\text{Fe}^{4+}$ ratio and establish an equilibrium
charge distribution. The solute drag effect may differ when a significant
amount of $\text{Fe}^{4+}$ is considered.

## Acknowledgement

The financial support of German Science Foundation (DFG) in the framework of
the Collaborative Research Centre 1548 (CRC 1548, project number 463184206)
and the Project 471260201 are acknowledged. The authors K. Wang and B.-X. Xu
greatly appreciate the access to the Lichtenberg II High-Performance Computer
(HPC) and the technique supports from the HHLR, Technische Universität
Darmstadt. The computing time on the HPC is granted by the NHR4CES Resource
Allocation Board under the project “special00007”.

## References

  * [1] T. Ishihara, H. Matsuda, Y. Takita, Doped $\text{LaGaO}_{3}$ perovskite type oxide as a new oxide ionic conductor, J. Am. Chem. Soc. 116 (9) (1994) 3801–3803.
  * [2] K.-D. Kreuer, Proton-conducting oxides, Annu. Rev. Mater. Res. 33 (1) (2003) 333–359.
  * [3] Z. Wang, Z. Liu, G. Zhao, Z. Zhang, X. Zhao, X. Wan, Y. Zhang, Z. L. Wang, L. Li, Stretchable unsymmetrical piezoelectric $\text{BaTiO}_{3}$ composite hydrogel for triboelectric nanogenerators and multimodal sensors, ACS nano 16 (1) (2022) 1661–1670.
  * [4] N. Y. Chan, M. Zhao, J. Huang, K. Au, M. H. Wong, H. M. Yao, W. Lu, Y. Chen, C. W. Ong, H. L. W. Chan, et al., Highly sensitive gas sensor by the $\text{LaAlO}_{3}/\text{SrTiO}_{3}$ heterostructure with Pd nanoparticle surface modulation, Adv. Mater. 26 (34) (2014) 5962–5968.
  * [5] S. W. Lee, J. H. Han, S. Han, W. Lee, J. H. Jang, M. Seo, S. K. Kim, C. Dussarrat, J. Gatineau, Y.-S. Min, et al., Atomic layer deposition of $\text{SrTiO}_{3}$ thin films with highly enhanced growth rate for ultrahigh density capacitors, Chem. Mater. 23 (8) (2011) 2227–2236.
  * [6] C. Hou, W. Huang, W. Zhao, D. Zhang, Y. Yin, X. Li, Ultrahigh energy density in $\text{SrTiO}_{3}$ film capacitors, ACS Appl. Mater. Interfaces 9 (24) (2017) 20484–20490.
  * [7] A. Kursumovic, E. Defay, O. J. Lee, C.-F. Tsai, Z. Bi, H. Wang, J. L. MacManus-Driscoll, A new material for high-temperature lead-free actuators, Adv. Funct. Mater. 23 (47) (2013) 5881–5886.
  * [8] S.-H. Baek, C.-B. Eom, Epitaxial integration of perovskite-based multifunctional oxides on silicon, Acta Mater. 61 (8) (2013) 2734–2750.
  * [9] R. Muenstermann, T. Menke, R. Dittmann, R. Waser, Coexistence of filamentary and homogeneous resistive switching in Fe-doped $\text{SrTiO}_{3}$ thin-film memristive devices, Adv. Mater 22 (43) (2010) 4819–4822.
  * [10] Z. Hu, Q. Li, M. Li, Q. Wang, Y. Zhu, X. Liu, X. Zhao, Y. Liu, S. Dong, Ferroelectric memristor based on $\text{Pt/BiFeO}_{3}/\text{Nb-doped}$ $\text{SrTiO}_{3}$ heterostructure, Appl. Phys. Lett. 102 (10) (2013).
  * [11] T. P. Mishra, S. Wang, C. Lenser, D. Jennings, M. Kindelmann, W. Rheinheimer, C. Broeckmann, M. Bram, O. Guillon, Ultra-fast high-temperature sintering of strontium titanate, Acta Mater. 231 (2022) 117918.
  * [12] M. J.-H. Reavley, H. Guo, J. Yuan, A. Y. R. Ng, T. Y. K. Ho, H. T. Tan, Z. Du, C. L. Gan, Ultrafast high-temperature sintering of barium titanate ceramics with colossal dielectric constants, J. Eur. Ceram. Soc. 42 (12) (2022) 4934–4943.
  * [13] M. Sayer, K. Sreenivas, Ceramic thin films: fabrication and applications, Science 247 (4946) (1990) 1056–1060.
  * [14] H. Li, X. Hu, Y. Wei, Z. Yu, X. Zhang, R. Droopad, A. Demkov, J. Edwards Jr, K. Moore, W. Ooms, et al., Two-dimensional growth of high-quality strontium titanate thin films on Si, J. Appl. Phys. 93 (8) (2003) 4521–4525.
  * [15] R. Chaim, M. Levin, A. Shlayer, C. Estournès, Sintering and densification of nanocrystalline ceramic oxide powders: a review, Adv. Appl. Ceram. 107 (3) (2008) 159–169.
  * [16] J. Maier, Physical chemistry of ionic materials: ions and electrons in solids, John Wiley & Sons, 2023.
  * [17] D. M. Smyth, The defect chemistry of metal oxides, 2000.
  * [18] F. Kröger, H. Vink, Relations between the concentrations of imperfections in crystalline solids, in: Solid State Phys., Vol. 3, Elsevier, 1956, pp. 307–435.
  * [19] J. P. Parras, R. A. De Souza, Grain-boundary diffusion of cations in fluorite-type oxides is faster but not always easier, Acta Mater. 195 (2020) 383–391.
  * [20] R. Sažinas, I. Sakaguchi, M.-A. Einarsrud, T. Grande, 134Ba diffusion in polycrystalline $\text{BaMO}_{3}$ (M= Ti, Zr, Ce), AIP Adv. 7 (11) (2017).
  * [21] R. Waser, R. Hagenbeck, Grain boundaries in dielectric and mixed-conducting ceramics, Acta Mater. 48 (4) (2000) 797–825.
  * [22] J. W. Cahn, The impurity-drag effect in grain boundary motion, Acta Metall. 10 (9) (1962) 789–798.
  * [23] M. Hillert, B. Sundman, A treatment of the solute drag on moving grain boundaries and phase interfaces in binary alloys, Acta Metall. 24 (8) (1976) 731–743.
  * [24] S. G. Kim, Y. B. Park, Grain boundary segregation, solute drag and abnormal grain growth, Acta Mater. 56 (15) (2008) 3739–3753.
  * [25] M. P. Zahler, S. M. Kraschewski, H. Stoermer, D. Gerthsen, M. Baeurer, W. Rheinheimer, Grain growth and segregation in Fe-doped $\text{SrTiO}_{3}$: Experimental evidence for solute drag, J. Eur. Ceram. Soc. 43 (4) (2023) 1613–1624.
  * [26] D. Jennings, M. P. Zahler, D. Wang, Q. Ma, W. Deibert, M. Kindelmann, C. Kübel, S. Baumann, O. Guillon, J. Mayer, et al., Grain boundary segregation in iron doped strontium titanate: From dilute to concentrated solid solutions, Acta Mater. 273 (2024) 119941.
  * [27] W. Rheinheimer, M. J. Hoffmann, Grain growth in perovskites: What is the impact of boundary transitions?, Curr. Opin. Solid State Mater. Sci. 20 (5) (2016) 286–298.
  * [28] W. Rheinheimer, M. J. Hoffmann, Non-arrhenius behavior of grain growth in strontium titanate: New evidence for a structural transition of grain boundaries, Scr. Mater. 101 (2015) 68–71.
  * [29] J. Fleig, The grain boundary impedance of random microstructures: numerical simulations and implications for the analysis of experimental data, Solid State Ion. 150 (1-2) (2002) 181–193.
  * [30] R. A. De Souza, The formation of equilibrium space-charge zones at grain boundaries in the perovskite oxide $\text{SrTiO}_{3}$, Phys. Chem. Chem. Phys. 11 (43) (2009) 9939–9969.
  * [31] K. Wang, G. Boussinot, C. Hüter, E. A. Brener, R. Spatschek, Modeling of dendritic growth using a quantitative nondiagonal phase field model, Phys. Rev. Mater. 4 (3) (2020) 033802.
  * [32] K. Wang, G. Boussinot, E. A. Brener, R. Spatschek, Quantitative nondiagonal phase field modeling of eutectic and eutectoid transformations, Phys. Rev. B 103 (18) (2021) 184111.
  * [33] V. Rehn, J. Hötzer, W. Rheinheimer, M. Seiz, C. Serr, B. Nestler, Phase-field study of grain growth in porous polycrystals, Acta Mater. 174 (2019) 439–449.
  * [34] Y. Yang, O. Ragnvaldsen, Y. Bai, M. Yi, B.-X. Xu, 3D non-isothermal phase-field simulation of microstructure evolution during selective laser sintering, npj Comput. Mater. 5 (1) (2019) 81.
  * [35] Y. Yang, T. D. Oyedeji, P. Kühn, B.-X. Xu, Investigation on temperature-gradient-driven effects in unconventional sintering via non-isothermal phase-field simulation, Scr. Mater. 186 (2020) 152–157.
  * [36] Y. Yang, M. Fathidoost, T. D. Oyedeji, P. Bondi, X. Zhou, H. Egger, B.-X. Xu, A diffuse-interface model of anisotropic interface thermal conductivity and its application in thermal homogenization of composites, Scr. Mater. 212 (2022) 114537\.
  * [37] T. D. Oyedeji, Y. Yang, H. Egger, B.-X. Xu, Variational quantitative phase-field modeling of nonisothermal sintering process, Phys. Rev. E 108 (2) (2023) 025301.
  * [38] K. Vikrant, W. Rheinheimer, R. E. García, Electrochemical drag effect on grain boundary motion in ionic ceramics, npj Comput. Mater. 6 (1) (2020) 165.
  * [39] K. Vikrant, W. Rheinheimer, H. Sternlicht, M. Bäurer, R. E. García, Electrochemically-driven abnormal grain growth in ionic ceramics, Acta Mater. 200 (2020) 727–734.
  * [40] X. Tong, D. S. Mebane, R. A. De Souza, Analyzing the grain-boundary resistance of oxide-ion conducting electrolytes: Poisson-Cahn vs Poisson-Boltzmann theories, J. Am. Ceram. Soc. 103 (1) (2020) 5–22.
  * [41] S. G. Kim, W. T. Kim, T. Suzuki, Phase-field model for binary alloys, Phys. Rev. E 60 (6) (1999) 7186.
  * [42] R. Merkle, W. Sitte, J. Maier, Water incorporation into materials with three mobile carriers: Two-fold relaxation of the electromotive force in $\text{Fe-doped SrTiO}_{3}$ and importance of hole trapping, Solid State Ion. 347 (2020) 115174.
  * [43] N. Moelans, B. Blanpain, P. Wollants, Quantitative analysis of grain boundary properties in a generalized phase field model for grain growth in anisotropic systems, Phys. Rev. B 78 (2) (2008) 024113.
  * [44] Y. Bai, Y. Zhao, W. Liu, B.-X. Xu, Two-level modeling of lithium-ion batteries, J. Power Sources 422 (2019) 92–103.
  * [45] T. J. Hughes, The finite element method: linear static and dynamic finite element analysis, Courier Corporation, 2012.
  * [46] S. C. Brenner, The mathematical theory of finite element methods, Springer, 2008\.
  * [47] D. Schwen, L. Aagesen, J. Peterson, M. Tonks, Rapid multiphase-field model development using a modular free energy based approach with automatic differentiation in moose/marmot, Comput. Mater. Sci. 132 (2017) 36–45.  doi:https://doi.org/10.1016/j.commatsci.2017.02.017.
URL https://www.sciencedirect.com/science/article/pii/S0927025617300885

  * [48] G. Giudicelli, A. Lindsay, L. Harbour, C. Icenhour, M. Li, J. E. Hansel, P. German, P. Behne, O. Marin, R. H. Stogner, J. M. Miller, D. Schwen, Y. Wang, L. Munday, S. Schunert, B. W. Spencer, D. Yushu, A. Recuero, Z. M. Prince, M. Nezdyur, T. Hu, Y. Miao, Y. S. Jung, C. Matthews, A. Novak, B. Langley, T. Truster, N. Nobre, B. Alger, D. Andrš, F. Kong, R. Carlsen, A. E. Slaughter, J. W. Peterson, D. Gaston, C. Permann, 3.0 \- MOOSE: Enabling massively parallel multiphysics simulations, SoftwareX 26 (2024) 101690.  doi:https://doi.org/10.1016/j.softx.2024.101690.
URL https://www.sciencedirect.com/science/article/pii/S235271102400061X

  * [49] R. De Souza, E. Dickey, The effect of space-charge formation on the grain-boundary energy of an ionic solid, Phil. Trans. R. Soc. A 377 (2152) (2019) 20180430.