# Defect design in ferroelectrics – new insights on agglomeration

Sheng-Han Teng  Anna Grünebohm

###### Abstract

Functional properties of ferroelectrics and their change with time depend
crucially on the defect structure. In particular, point defects and bias
fields induced by defect dipoles modify the field hysteresis and play an
important role in fatigue and aging. However, a full understanding on how
order, agglomeration and strength of defect dipoles affect phase stability and
functional properties is still lacking. To close these gaps in knowledge, we
screen these parameters by ab initio based molecular dynamics simulations with
the effective Hamiltonian method for the prototypical ferroelectric material
(Ba,Sr)TiO3. Our findings suggest that the active surface area of the defects,
rather than the defect concentration is the decisive factor. For a fixed
defect concentration, clustering reduces the active surface area and thus the
defect-induced changes of phase stability and field hysteresis. Particularly
planar agglomerates of defects appear as promising route for the material
design as their impact on the field hysteresis can be controlled by the field
direction and their impact on the phase stability shows a cross-over with the
strength of the defect dipoles. For this agglomeration, we furthermore find
that pinched field hysteresis, which is beneficial, e.g. for energy storage,
can be achieved for a wide range of defect dipole strengths and thus is not
too sensitive to the choice of dopants.

###### keywords:

Ferroelectric, Phase diagrams , Field hysteresis , Point defects , Molecular
Dynamics Simulations

††journal: Acta Materialia

organization=Interdisciplinary Centre for Advanced Materials Simulation
(ICAMS) and
Center for Interface-Dominated High Performance Materials (ZGH), Ruhr-
University Bochum, addressline=Univeristätstr. 150, city=Bochum,
postcode=44780, country=Germany

##  1 Introduction

The exceptional dielectric, caloric, and piezoelectric functional responses of
ferroelectric perovskites in the vicinity of structural phase transitions, as
well as their adjustable field hysteresis in the ferroic phases, make these
materials important for a wide range of applications [1, 2] including energy
storage and electronic devices. [3, 4] Adjusting transition temperatures and
functional properties is commonly done by chemical substitution.[5, 6, 7]
Thereby, one has to distinguish between isovalent substitution, e.g.
substitution of Ba in BaTiO3 by Sr or Ca ($\rm Sr_{Ba}^{\times}$, $\rm
Ca_{Ba}^{\times}$), or Ti by Zr ($\rm Zr_{Ti}^{\times}$), and aliovalent
substitution, e.g. substitution of Ti by Fe or Nb ($\rm Fe_{Ti}^{\prime}$,
$\rm Nb_{Ti}^{\cdot}$). In the former case, the substituents introduce local
distortions and commonly reduce transition temperatures. [7, 8] Although
changes of the local polarization [8] and even non-polar cells may occur,[9,
10, 6] the dynamics of the local switching are not changed significantly for
small defect concentrations, as all dipoles switch by ionic displacements
happening on the ps-scale.[11, 12] In the latter case, depending on the type
and distribution of the dopant and the thermal treatment of the system,
changes of transition temperatures and internal bias fields with shifting and
pinching of hysteresis loops have been reported. [13, 14, 15] As internal
restoring forces and pinched hystereses are beneficial for energy storage [3]
as well as large reversible piezoelectric[16] and electrocaloric responses,
[17, 18] this field arose considerable attention. It is accepted that these
observations are related to charge compensation by vacancies and the formation
of defect complexes with local dipole moments, such as
(Fe${}_{Ti}^{\prime}$-V${}_{O}^{\cdot\cdot}$)⋅ and
(Mn${}_{Ti}^{\prime\prime}$-V${}_{O}^{\cdot\cdot}$)×.[16, 19] Switching of
these defect complexes requires slow ion migration which may take more than
$10^{7}$ s at room temperature.[19, 20, 21] Thus these defects align with the
surrounding polarization during aging and later act as restoring forces on
domain and phase structure. It is furthermore accepted that changes of the
functional properties with time, e.g., functional fatigue, are related to ion
migration and redistribution of defects. [22, 13] In experimental samples, the
distribution of defects depends on processing, thermal treatment, field-
loading and internal and external boundaries, e.g. electrodes, grains, or
dislocations [23, 13]. For example, the uniform random distribution after ion-
bombardment, [24] planar arrays of oxygen vacancies in reduced atmosphere,[25]
defect planes in a severely fatigued sample,[26] clustering of oxygen
vacancies,[27] and agglomeration of defects at dislocation lines [28] have
been reported.

To tap the full potential of defect-engineering and to improve the
understanding of functional fatigue, a microscopic understanding is mandatory
Since most microscopic modeling so far focuses on the alignment of single
defect dipoles, even after decades of research, this understanding is still
incomplete. In the present study, we use an ab initio derived effective
Hamiltonian to screen the impact of defect distribution on ferroelectric phase
diagrams and field hystereses in BaTiO3 solid solutions. Thereby we focus on
point defects without local polarization and defect dipoles which re-orient
orders of magnitude slower than the switching by Ti off-centering and can thus
be approximated by fixed dipoles,[19, 21] see Fig.1. We show that although the
impact of the defects on phase transition and field hysteresis depends on
defect strength and concentration, the distribution and agglomeration of
defects play an even larger role.

##  2 Methods

###  2.1 Molecular dynamics simulations

We use the effective Hamiltonian from Ref. 29, which is parameterized based on
first-principles calculations for the solid solution of (Ba,Sr)TiO3:[30]

$\displaystyle H^{eff}$
$\displaystyle=\frac{M_{dipole}^{*}}{2}\sum_{\bm{R},i}\dot{u}_{i}^{2}(\bm{R})$
$\displaystyle+V^{self}(\\{\bm{u}\\})+V^{dpl}(\\{\bm{u}\\})+V^{short}(\\{\bm{u}\\})$
$\displaystyle+V^{elas,homo}({\eta_{\alpha}})+V^{elas,inho}(\\{\bm{w}\\})$
$\displaystyle+V^{coup,homo}(\\{\bm{u}\\},{\eta_{\alpha}})+V^{coup,inho}(\\{\bm{u}\\},\\{\bm{w}\\})$
$\displaystyle+V^{mod,inho}(\\{\bm{w}\\},\\{s\\})-Z^{*}\sum_{\bm{R},i}\mathcal{E}\cdot
u_{i}(\bm{R}),$ (1)

where $\\{\bm{u}\\}$, $\\{\bm{w}\\}$, $\eta_{\alpha}(\alpha=1,...,6)$,
$\\{s\\}$, $Z^{*}$ and $\mathcal{E}$ are local soft mode vectors in Cartesian
coordinates $i(=x,y,z)$, local acoustic displacement vectors, the tensor of
homogenous strain, the relative number of Ba and Sr atoms in each unit cell
(u.c.), the Born effective charge associated with the soft mode, and an
external electrical field, respectively. The local mode gives the local
polarization as $\bm{P}=eZ^{*}\bm{u}/a^{3}$ where $e$ is the elementary charge
and $a$ is the lattice constant. Without loss of generality, we study
Ba0.7Sr0.3TiO3 with its Curie temperature close to room temperature using a
random distribution of Sr in combination with the chemical modulation
$V^{mod,inho}$ which depends on the Sr concentration in the unit cell and a
pressure correction of $-1.2$ GPa, cf. Ref. 30. The other terms in Eq. (1) are
the self-energy, $V^{self}$, the long-range dipole-dipole interaction,
$V^{dpl}$, the short-range harmonic interaction between local modes,
$V^{short}$, the elastic energy from homogeneous and inhomogeneous strains,
$V^{elas,homo}$ and $V^{elas,inho}$, the coupling terms between local soft
mode and strain, $V^{coup,inho}$ and $V^{coup,inho}$. As the strain is
optimized internally, only the dynamics of the dipole moments are explicitly
treated in molecular dynamics simulations, using $M^{*}_{dipole}$ the
effective mass of the local soft mode. This coarse-graining is illustrated in
Fig. 1.

Figure 1: Coarse-graining and modeling of defects: The atomic displacements in
all unit cells of (Ba,Sr)TiO3 are coarse-grained to free dipoles (dark blue
arrows) and local strain (optimized internally). If Ti is substituted by
another atom ($M$), the local dipole moment can be reduced. We consider the
limiting case of a non-polar cell (red dot). Defect complexes formed by dopant
$M$ and an oxygen vacancy ($V_{O}$) carry finite dipole moments which switch
orders of magnitude slower than the surrounding material and are modeled by
frozen-in dipoles (red arrows).

A supercell with $L_{x}\times L_{y}\times L_{z}=48\times 48\times 48$ u.c. is
combined with periodic boundary conditions and the dipoles are initialized
with a normal distribution with
${\left\langle{u_{x}}\right\rangle}={\left\langle{u_{y}}\right\rangle}={\left\langle{u_{z}}\right\rangle}=0.00$
Å and
${\left\langle{u_{i}^{2}}\right\rangle}-{\left\langle{u_{i}}\right\rangle}^{2}=0.03$
Å at 350 K or with
${\left\langle{u_{x}}\right\rangle}={\left\langle{u_{y}}\right\rangle}={\left\langle{u_{z}}\right\rangle}=0.10$
Å and
${\left\langle{u_{i}^{2}}\right\rangle}-{\left\langle{u_{i}}\right\rangle}^{2}=0.01$
Å at 50 K to perform cooling and heating simulations. For each temperature, we
use a time step of $\Delta t=$2 fs and thermalize the system for 40000 steps
in the canonical ensemble using the Nosé-Poincaré thermostat before we
calculate the average properties in additional 20000 steps. The final
configuration is the input for the next temperature step using $\Delta T=\pm
1$ K. We monitor the absolute values of the macroscopic polarization
|${\left\langle{P_{i}}\right\rangle}$| along direction $i$, and the mean of
the magnitude of the local polarization
${\left\langle{|P_{i}|}\right\rangle}$. Both quantities differ in the presence
of domains. Unlike stated otherwise the temperatures of ferroelectric to
ferroelectric transitions are identified from the maximum/minimum of strain in
the high/low-temperature phase and the paraelectric-ferroelectric transition
temperature is determined from the maximal derivative of polarization with
respect to temperature.

For chosen equilibrated configurations the field-hysteresis in
${\left\langle{100}\right\rangle}$-type fields with a magnitude of $150$ kV/cm
is recorded starting with the negative field direction at the ramping rate of
$1.5E16$ V/m/s. Note that the virgin curve for this first poling is not
discussed in the following. As shown in Figure 2, each hysteresis loop is
characterized by the values of polarization without external bias
($P_{r}^{\pm}$) and in the maximal applied field ($P_{max}^{\pm}$), the
coercive fields ($E_{c}^{\pm}$), which we determine by the maxima of $dP/dE$,
the widths of the hysteresis loops($\Delta
E_{c}=|E_{c}^{(\rightarrow)}-E_{c}^{(\leftarrow)}|$), the bias fields
($E_{bias}=\frac{1}{N}\sum_{i}^{N}E_{c}^{(i)}$), the energy density in a
unipolar field-cycle which can be reversibly stored
($W_{rec}^{\pm}=\int_{P_{r}^{\pm}}^{P^{\pm}_{max}}EdP$) and work losses if
switching from remanent to saturation polarization
($W_{loss}^{\pm}=\int_{P_{r}^{\mp}}^{P_{max}^{\pm}}EdP-W_{rec}^{\pm}$).
Furthermore, it is convenient to define the energy storage efficiency $\eta$
as $\eta=W_{rec}/(W_{rec}+W_{loss})$. Data as well as Jupyter notebooks used
for the analysis are available via Zenodo [31] and Gitlab [32].

Figure 2: Characteristic field hysteresis representative for the pristine
material (dark blue) compared to the pinched hysteresis in the presence of 2D
defect agglomerates (red) without remanent polarization ($P_{r}^{\pm}=0$), and
the shifted hysteresis with an internal bias field ($E_{bias}$) induced by
defect dipoles collinear with the applied field (gray). $E_{c}^{\pm}$, and
$P_{r}^{\pm}$ are coercive fields and remanent polarizations for both field
directions ($E^{\pm}$) and $|P|_{max}$ is the polarization at the maximal
field strength $E_{max}$. Hatched and colored regions illustrate work loss
($W_{loss}$) and recoverable energy ($W_{rec}$) in a unipolar field loop.

###  2.2 Defect dipoles

To depict the range of reasonable defect strengths and concentrations, we
first categorize reports on defects from literature. Typical oxygen vacancy
and doping concentration are in the range of 0–2 % (i.e. 1014–1021 cm-3).[13,
24, 7, 14] It has furthermore been reported that defect dipoles increase [33,
34] or decrease [34] the remanent polarization of the pristine material
($P_{0}$). To cover this range, we model defect strengths from $0$ to
$1.8P_{0}$ with $P_{0}=28.41$ $\mu$C/cm2 the polarization we find for
(Ba${}_{0.7},$Sr0.3)TiO3 at 210 K (tetragonal phase). If not stated otherwise,
we use $u=0.1$ Å, i.e. $0.9$ $P_{0}$ as example.

We take two types of point defects into account. Fig. 1 shows possible
microscopic realizations for these scenarios: First, substituents with
different radii or electronic structures, e.g., Ce4+ or Zr4+ on B-site, may
prevent the local off-centering in the stoichiometric material and suppress
the local dipole moment (illustrates as circle).[9, 6] Second, pairs of a
substituent with different valency and an oxygen vacancy may form defect
dipoles $M+V$, which align with the polarization direction in the tetragonal
phase during aging [20] and follow different dynamics than the intrinsic
polarization. We approximate these by frozen-in dipoles which are all aligned
along one direction, cf. Ref. 18.

To screen the impact of the defect distribution, we focus on defects randomly
distributed in the material and clustering in planes and cubes. Fig. 3
collects all defect distributions and concentrations. The largest coupling
between defects and surrounding local dipoles is expected along the
polarization direction $i$, as the defects induce a local polarization
gradient $\nabla_{i}P$ and thus a depolarization field. Therefore we define
the active surface area of a defect configuration as the area of the surface
with its normal parallel to a polarization components $i$, e.g. along the
tetragonal axis or along any of the Cartesian directions in the rhombohedral
phase, measured in u.c., i.e. unit cells with an area of about $15.6$ Å2.

[abs] (a)$\Leftrightarrow$(b)[abs] (c)$\Leftrightarrow$(d)$\Leftrightarrow$(e)

Distribution |  Concentration (%) | Active surface area
---|---|---
(each unit with an area of 15.9 Å2)
T | O | R
(a) Cube |  | 0.058–1.9
---
($a$=4–13 u.c.)
32–338 | 32–338 | 32–338
(b) Rand.3D | 0.5–2 | 1106–4424 | 1106–4424 | 1106–4424
(c) Stripe |  | 0.26–1
---
($w$=6–24 u.c.)
0 | * | 576–2304
(d) Plane |  | 2
---
($w=\infty$)
0 | 0 | 4608
(e) Rand.2D | 0.26–1 | * | * | 576–2304

Figure 3: Sampled distributions as well as considered concentrations of
defects and their active surface areas (table) for: (a) cubic agglomerates,
(b) random distribution in 3D, (c)–(e) agglomeration in one plane for (c)
compact stripes, (d) full planes and (e) randomly. Thereby $a$ and $w$ are the
size and width of the cube and the stripe, respectively. In tetragonal (T),
orthorhombic (O) and rhombohedral (R) phases the active surface area is
defined by the surface of the defects along the polarization direction. Each
u.c. contributes an active area of about (3.9435 Å)2. As discussed in the
text, some 2D configurations have no effective active area in T or O phase.

##  3 Results

###  3.1 Phase diagrams

[abs] 32128200338110622124424(a)* * * *576115223044608(b)

Figure 4: Temperature - defect concentration phase diagrams in the presence of
non-polar defects ($|u_{D}|=$0.0 Å) for (a) cubic agglomerates (solid lines)
and random distribution in 3D (dashed lines) and (b) stripes of defects (solid
lines) and random distribution in 2D (dashed lines). Red and blue lines show
transition temperatures for heating and cooling, respectively. Numbers give
the approximate number of active surface units in (a) T, O (b) and R phases.
Note that the defect agglomerates in (b) have no effective active surface area
for T and O phases as discussed in the text.

We start our analysis with the phase stability without external electrical
fields. For defect-free (Ba0.7,Sr0.3)TiO3, we find transition temperatures
from cubic to tetragonal ($T_{C}^{C-T}$) to orthorhombic ($T_{C}^{T-O}$), and
to rhombohedral ($T_{C}^{O-R}$) phases at 293 K, 210 K, and 168 K within the
spread of experimental values of about 285–308 K, 217–222 K and 154–172 K,
respectively.[35, 36, 37, 38] These transition temperatures are modified with
the concentration of non-polar defects, see Fig. 4. Thereby, changing the
distribution of defects from random to 3D agglomerates (dashed and solid lines
in subfigure (a)) or to 2D agglomerates (subfigure (b)) has a larger impact on
the transition temperatures than changing their concentration. First,
increasing the concentration of randomly distributed defects induces an
approximately linear reduction of $T_{C}^{C-T}$ of more than 35 K for 2 %
defects. Second, agglomeration in 2D mainly reduces $T_{C}^{O-R}$, e.g., by 33
K for 1 % of defects and by 92 K for a full plane of defects. Third, cubic
agglomerates with up to 2% concentration have a minor impact on all phase
transition temperatures. Large agglomerates in 3D and 2D induce a slight
increase of $T_{C}^{T-O}$, while $T_{C}^{R-O}$ increases slightly in the
presence of large 3D agglomerates or a high concentration of randomly
distributed defects. Note that the qualitative changes of transition
temperatures for randomly distributed defects are in good agreement to the
experimental results published for the substitution of Ti by Zr.[7]

Most of these observations can be understood by the change of the active
surface areas of the defects in the different phases, see Fig. 3. First of
all, the active surface area along one polarization component impacts
$T_{C}^{C-T}$. For a given defect concentration, it is maximal for a
distribution without direct defect neighbors. Indeed, random distributions
induce the largest change of $T_{C}^{C-T}$ and of macroscopic polarization and
strain, see Fig. 5 and cf. Fig. 12 (a). If the defects cluster in one cube,
the active surface area is reduced, e.g. by factor of 10 for 1 % defects. Even
for 2 % defects in such a cluster, changes of macroscopic polarization and
strain are barely visible, cf. Fig. 13. Both, in case of random distribution
and 3D clustering, the active surface area does not change if going from T to
O or R phase and thus the lower transition temperatures are not too sensitive
to the defect concentration. On the other hand, agglomerates in 2D reduce the
symmetry of the system and one has to distinguish polarization normal to the
defect plane with a large active surface area and polarization parallel to the
defect plane. For T and O phases, the defects stabilize the parallel
polarization directions ($P_{y}$ and $P_{x}=P_{y}$) without large active areas
and thus $T_{C}^{C-T}$ and $T_{C}^{T-O}$ are not modified. Only in the R
phase, one polarization component has to point along the normal, is subject to
a maximal active surface area and $T_{C}^{O-R}$ is reduced with the defect
concentration.

[abs]  (a)$\ast$

$P_{i}$

Normalized histogram$T=150K$(c)(d)

[abs]  (b)$\ast$

$P_{i}$

Normalized histogram$T=100K$(e)(f)

Figure 5: Changes of (a)-(b) polarization components
$|{\left\langle{P_{i}}\right\rangle}|$, with $i:x,y,z$, in cooling simulations
and (c),(e) of dipole histograms, induced by non-polar defects ($|u_{D}|=$0.0
Å) for (a),(c),(d) random distribution and (b),(e),(f) stripes in $z$-planes
with varying width ($w$). For (c) 2 % defects at 150 K and (e) 1 % defects
($w=24$ u.c.) at 100 K, histograms of local $P_{z}$ (right) are compared to
the pristine reference system (left). (d) and (f) show corresponding time-
averages of local dipoles close to the defects color encoded by the sign of
$P_{z}$.

While the main trends of the transition temperatures are thus related to the
active surface areas of the defects, also their dilution impacts phase
stability and dipole configuration. Figure 5 compares the change of
macroscopic polarization between defects being (a) randomly distributed and
(b) agglomerated in 2D stripes. The underlying histograms of one polarization
component in the presence of 2% defects are shown to the right in subfigures
(c) and (e), while the histograms to the left correspond to the pristine
material at the corresponding temperatures. In case of random distribution,
about $4.4$ % of the dipoles are reduced (peak around 7.2 $\mu C/cm^{2}$).
Closer inspection shows that these are mainly nearest neighbors of defects
along the polarization direction, cf. Fig. 5 (d). As the defects start to
cluster, more and more defects are reduced along this direction and the
coupling between defects and dipoles becomes more long-ranged, see subfigure
(f). Already a stripe with $w=6$ u.c. (i.e. 0.5% defects) induces gradual
rotations of dipoles in at least 3 layers next to the agglomerate. Above a
critical width of about $w=12$ u.c., negative domains with a diameter of $8$
u.c. form at the defect agglomerate, and for $w=24$ u.c. the dipoles normal to
the defect are modified in the whole system by the formation of domains with
R71 walls along [110], cf. Fig. 15. This increase of the interaction length
with the compactness of the defects is only induced along the polarization
directions. Thus, while a compact stripe of defects has a finite interface
along the $x$-direction, the induced modifications of $P_{x}$ are confined to
that plane and do not couple to the macroscopic polarization. The effective
active area marked as ($*$) in Fig. 4 is thus zero. Note that also large 3D
agglomerates, e.g. $a=13$ u.c. with an active surface unit of 338 u.c. (cf.
Fig. 15) induce rotations of the dipole in several unit cells and multi-domain
structures, in this case for all three polarization directions.

Going from a perfect defect plane to the pristine material is also possible by
the dilution of defects on the plane, see dashed lines in Fig. 4 (b). While a
multi-domain structure forms for more than 25 % coverage of compact
agglomerates, more than 50 % coverage is needed in a diluted plane. Thereby
also the change of $T_{C}^{O-R}$ for small defect concentrations is larger in
the former case. On the other hand, the cross-section of the diluted defect
plane is larger for a given defect concentration. Once domain walls form, the
whole system splits into 2 equidistant R71 stripes and $T_{C}^{O-R}$ decreases
faster for intermediate defect strengths. These changes of the interaction
range with the dilution of the defects can be related to the dilution of the
depolarization field: On the one hand, the depolarization field by separated
defects only couples to the local environment. On the other hand, the fields
of dense defects add up, which increases the interaction length and results in
the formation of domains above a critical size and density of the agglomerate.

Non-polar defects also modify the character of the transitions and their
thermal hysteresis. In the pristine material, all transitions are of first
order and we find thermal hysteresis of 5 K, 20 K and 30 K, correspondingly.
Generally, the hysteresis is systematically reduced with the concentration of
defects, for example to 2 K, 11 K and 11 K for 2 % of randomly distributed
non-polar defects. Particularly, the character of the O–R transition can be
modified by agglomerated defects, cf. continuous change of polarization for
wide agglomerates in Fig. 4. For these configurations, hysteresis closes for
less than 0.5 % of defects due to the reduction of $T_{C}^{O-R}$ during
heating. Analogously, thermal hysteresis is below our resolution for about 0.5
% and 1 % defects randomly distributed in a plane or clustered in a cube,
respectively. This can be understood by the discussed change of the
polarization component normal to the agglomerate which induces monoclinic
distortions and by domain wall formation in the R phase, as both these
modifications may reduce the energy barrier for the transition.[39, 2] Note
that a hysteresis of about 15 K remains at the T–O transition even for a full
defect plane due to the weak coupling between defects and the polarization
parallel to the plane.

In summary, the impact of non-polar defects on the ferroelectric phase diagram
is dominated by the effective active surface area of the defect configuration
and, in case of compact defect agglomerates, depolarization fields modify the
domain structure. Similar, but less pronounced changes of the phase stability
have to be expected for locally reduced, but finite local dipoles switching by
polar off-centering.

[abs] 32128200338110622124424(a)576115223044608(b)

Figure 6: Temperature - defect concentration phase diagrams for polar defects
($|u_{D}|=$0.1 Å$\hat{z}$). Symbols mark transition temperatures for cooling
(blue) and heating (red). Note that no C-T transition can be detected for
large defect concentrations and magnitudes beyond the critical point. Random
ordering (solid lines) and agglomeration (dashed lines) are compared for (a)
defects in 3D and (b) defects in a $z$-plane. Numbers give the approximate
number of active surface units in all phases for (a) cube and random
distribution (italic) and (b) for stripes.

Even larger changes of the phase stability are possible in the presence of
defect dipoles related to M-VO complexes with different switching kinetics, as
shown exemplary for defects with $|u_{D}|=0.1$ (Å) $\hat{z}$, i.e.
$u_{x}=u_{y}=0$, in Fig. 6. Again, small cubic agglomerates have a minor
impact on all phases due to their small active surface area. Besides that, the
trends differ from the non-polar case. As the polarization of the T phase
aligns with the direction of the defect dipoles, all other orderings have
similar active surface area and $T_{C}^{C-T}$ increases until a critical point
is induced by about 1% defects. Additionally, the transition temperatures
under heating for $T_{C}^{T-O}$ and $T_{C}^{O-R}$ decrease for small defect
concentrations. Already 0.5 % of defects randomly distributed in 2D or 3D thus
reduce the thermal hysteresis of the T-O transition to 4 K. Analogous trends
have been reported for external fields along $\langle 100\rangle$.[40, 41, 42]
In the tetragonal phase, the defect dipoles thus act as an internal electrical
field along the direction of the defect dipoles whose strength can be adjusted
with the concentration of the defects. However, in contrast to a homogenous
external field also active surface area of defects and defect dilution
contribute to the coupling between the internal field and the free dipoles.
Thus both the change of $T_{C}^{C-T}$ and the critical defect concentration
vary slightly between stripes and random distributions. In case of 2D
agglomeration, the effective active surface area along $x$ and $y$ is again
zero and the transition temperatures to O or R phases under cooling are not
sensitive to the defect concentration. On the other hand, randomly distributed
defects have an active surface area at their non-polar interfaces and thus
induce a reduction of $T_{C}^{T-O}$ of more than $55$ K for 2 % defects and
also $T_{C}^{O-R}$ goes down with the defect concentration, see Fig. 6 (a).
For 2D agglomerates with, e.g., $|u_{D}|=0.1$ (Å) $\hat{x}$, the defects act
similar to non-polar defects, cf. Fig. 11, with polarization along $x$ and
$xy$ in the T and O phases and small changes of the transition temperatures,
while $T_{C}^{O-R}$ decreases by $36$ K ($1$ % of defects) and by $98$ K for a
full defect plane for the R phase with polarization along $z$. Comparison
between these defects and the non-polar ones shows that the former induces a
slight increase of $T_{C}^{C-T}$ and a reduction of the thermal hysteresis at
the T-O transition, while only in the latter case $T_{C}^{T-O}$ increases
slightly for large defect concentrations. All these observations can be
related to the small internal fields along $x$.

[abs]  (a)$\ast$(b)[abs] (c)(d)

Figure 7: Impact of the strength of 1 % defects ($|u_{D}|$ $\hat{z}$) on the
ferroelectric phase diagram. Colors highlight transition temperatures for
cooling (blue) and heating (red). Agglomeration (solid lines) and random
ordering (dashed lines) are compared for (a) defects in a cube or 3D and (b)
defects in one $z$-plane. In (b) horizontal ($=$) and vertical ($||$) hatches
mark T and O phases with $P_{z}=0$ and $|P_{z}|>0$, respectively. (c) and (d)
show corresponding time-averaged local dipoles at 200 K for weak defect
dipoles ($|u_{D}|=0.03$ Å $\hat{z}$) and strong defect dipoles ($|u_{D}|=0.1$
Å $\hat{z}$), respectively. ∗Note that beyond $0.04$ Å the sample with 2D
random defects is beyond the critical point of the C-T transition.

The second obvious route to modify the internal field induced by the defect
dipoles is tuning of the defect strength. For random defect distribution and
cubic agglomerates, the increase of the internal field with defect strength
induces the same changes of the phase diagram as the change of the defect
concentration, see Fig. 7 (a). However, for 2D agglomerates with polarization
along $z$, where is a cross-over from weak to strong defect dipoles, see Fig.
7 (b). Below $u_{z}=0.03$ Å, e.g. about 27% of the saturation polarization of
the pristine material, the defect dipoles act similar to non-polar defects
with polarization in the T and O phases pointing along $P_{y}$ and
$P_{y}=P_{x}$, no effective active surface area and no change of $T_{C}^{C-T}$
and $T_{C}^{T-O}$. Only in the R phase one polarization component points
parallel to the defects and $T_{C}^{O-R}$ increases with their strength.
Thereby the O–R transition becomes continuous, cf. Fig.14. For strong defects
($u>0.03$ Å), the polarization in the T phase point along $z$ and
$T_{C}^{C-T}$ increases up to the critical point, while $T_{C}^{O-R}$ and
$T_{C}^{T-O}$ decrease with the internal defect field. At the cross-over
point, one has to expect that states with $P_{y}$ and $P_{z}$ are degenerated
which may induce multi-critical points with coexisting C, T and O as well as
with coexisting R, O and T phases. The exact cross-over points are beyond our
scope, for the chosen resolution of 0.01 Å, the transition temperatures of the
C, T and O as well as of T, O and R phases indeed approach each other.

In summary, in case of 2D agglomeration, non-polar and weak defect dipoles
have the same impact on the phase stability while strong defects act as an
internal bias field, while defects randomly distributed in 3D induce an
internal bias field for all field strengths.

###  3.2 P-E Hystereses

[abs]

$P_{i}$ ($\mu C/cm^{2}$)

$E_{z}$ ($kV/cm$)(a)[abs]  $E_{x}$ ($kV/cm$)(b)[abs]  $E_{z}$ ($kV/cm$)(c)

Figure 8: Changes of the field hysteresis at $T=250$ K if going from the
pristine material (red) to systems with different distribution of 1 % defects
(green: cube; blue and cyan: stripes and randomly in $\hat{z}$-planes; purple:
randomly) for (a) non-polar defects and (b)–(c) polar defect dipoles
($|u_{D}|=$0.1 Å$\hat{z}$) for the electric field (b) perpendicular ($E_{x}$)
and (c) collinear ($E_{z}$) to the defect dipoles. Only the polarization
component $P_{i}$ along the field direction is shown.

[abs]  (a)(c)(e) (b)(d)(f)

Figure 9: Characteristics of the field hysteresis for 1 % (a,c,e) non-polar
and (b,d,f) polar defects ($|u_{D}|=$0.1 Å$\hat{z}$) at $T=250$ K for $E=150$
kV/cm sorted by efficiency. In (a)–(b) $W_{rec}$ (green) and $W_{loss}$ (gray
hatched) and the efficiency $\eta$ (black dots, right axes), in (c)–(d) the
remanent polarization $|P_{r}|$ and in (e)–(f) the width of the hysteresis
($\Delta E_{c})$ and the bias field $E_{bias}$ for different defect
distributions are compared to the pristine material. Red and blue correspond
to positive and negative field directions which differ for the collinear
alignment of defects and the external field ($E_{z}$).

How do the different defect configurations affect the field response of the
material? As reference, red lines in Fig. 8 show the field hysteresis of the
pristine material in the T phase at 250 K. The $P-E$ curve is symmetric,
centered at the origin, and has large coercive fields and remanent
polarization. Depending on their ordering and strength, defects considerably
modify these characteristics, see other colors in Fig. 8 and the corresponding
strain–field ($\sigma$-E) curves in Fig. 17. Non-polar or weak defect dipoles,
reduce remanent polarization ($|P_{r}^{\pm}|$) and particularly the width of
the hysteresis loops ($\Delta E_{c}$) for all orderings while the hysteresis
stays symmetric as shown in subfigures (a). Analogous to the phase diagrams,
the differences between the configurations are related to the active surface
area of the defects. The active surface area and the impact on the field
hysteresis is small for cubic agglomerates and 2D agglomerates with the field
applied in the defect plane. In these cases the changes of $\Delta E_{c}$ (and
$P_{r}$) are below $18$ % ($3$ %) and $20$ % ($3$ %), respectively, cf. Fig.
16. Larger modifications are induced by randomly ordered defects that have a
larger active surface area. The largest changes are induced by 2D agglomerates
if the field is applied along $z$. This defect distribution induces
depolarization fields and acts as a reversible restoring force on the
polarization direction. Thus, the polarization switches between the T phase
with $P_{x}$ and a monoclinic state with $\pm P_{z}>P_{x}$, cf. Fig. 18. The
response $P_{z}(E_{z})$ is thus pinched to a double-loop hysteresis. Thereby,
for the chosen defect concentration of 1%, the randomly distributed defects in
one plane impact the surrounding dipoles more than stripes, cf. discussion on
Fig. 5, and stabilize a reversible linear change of $P_{z}$ with the field for
more than $\pm 96$ kV/cm.

As shown in subfigure (b), strong defect dipoles along $z$ furthermore induce
full pinching of the hysteresis for a perpendicular field along $x$. In this
case, internal and external fields compete with each other for all orderings.
For both 2D agglomerates and a random distribution of defects, the remanent
polarization along the field direction is zero and the strain response is
maximized, cf. Fig. 17 and Fig. 18. There is no hysteresis for randomly
distributed defects and a double hysteresis loop is induced in field ranges of
less than 10 kV/cm by compact stripes. Strong defect dipoles collinear to the
field, shown in subfigure (c), reduce the width of the field hysteresis, e.g.
by about 28 % for the cubic agglomerates and by about 45 % for random ordering
in 3D. Furthermore, these defects act as internal bias fields and shift the
center of the hysteresis to negative $E_{z}$. Thereby larger changes are
induced in $E_{c}^{+}$, i.e. parallel to the defect direction, and $E_{c}^{+}$
is shifted to negative field directions for random defect distributions. For
the cube, with its small active surface area only this positive coercive field
is reduced. Increasing the active surface area, not only $E_{c}^{+}$ shifts
further, but also $P_{r}^{-}$ is reduced. These trends can be understood by
pinned dipoles. Already for the cube, its neighbors along the polarization
direction do not switch for the chosen field strength and thus act as
nucleation centers for the macroscopic switching towards the positive
polarization direction lowering the coercive field. Going from cube to random
in 3D and finally to agglomeration in 2D, the interaction length between
defects and surrounding dipoles and thus the pinned volume increase, cf. Fig.
19. Additionally, the increase of $|E_{c}^{-}|$ induced by the internal bias
field, increases from stripe to random in plane to random in 3D. This again
shows that the active surface area governs the main trends of polarization and
phase stability, while the subtle details also depend on the compactness or
dilution of defects.

Figure 9 compares the characteristics of the hysteresis loops for all defect
configurations. For the pristine material, a work loss of about 3 J/cm-3 is
dissipated if the polarization switches from $P_{r}^{-}$ to the saturation
polarization $P_{max}^{x}$ and vice versa for negative field directions, while
the recoverable work as well as reversible changes of polarization and strain
without switching are negligible, in qualitative agreement with previous
reports.[43] On the one hand, non-polar or weak defects in most configurations
induce a small reduction of the loss only. The efficiency $\eta$ is thus
approximately constant and increases slightly for random defects. On the other
hand, large efficiencies are possible for the double loop hysteresis found for
2D ordering. Thereby, the smallest hysteresis in case of defects randomly
distributed in one plane allows for an efficiency of about 94 %. Analogously,
pinched hysteresis with reduced losses and maximal $\Delta P$ and $\eta$ are
stabilized by strong defects along $z$ in a field along $x$ as long as the
active surface area is not reduced by clustering in 3D. Maximal recoverable
energy and efficiency of 1.07 J/cm-3 and 96 % are possible for defects
randomly distributed in the sample. Finally, for strong defects and fields
along $z$, one has to distinguish between positive and negative field
directions. In the former case shown in red, no polarization switching can be
induced during unipolar field loops for randomly distributed defects and thus
the used definition of loss and efficiency would predict optimal
configurations. However, this is not relevant as the possible recoverable
energy as well as changes of $P$ and strain are negligible. In case of the
negative field shown in blue, maximal values of $W_{rec}$ and a moderate
efficiency are possible due to the shifted and reduced hysteresis for random
defect distributions.

In summary, the field characteristics depend strongly on the defect
arrangement and the relative direction of the field. Weak defects induce the
largest modification of the hysteresis if they agglomerate in 2D, while the
pinching of the hysteresis for strong defect dipoles perpendicular to the
field is maximal for randomly oriented dipoles. Finally, strong defects
parallel to the field induce the largest shift of hysteresis for random
distribution, while 2D agglomerates induce the largest asymmetry.

##  4 Conclusion

While it is commonly accepted that point defects and defect complexes in
ferroelectric perovskites modify ferroelectric phase diagrams and field
hysteresis, the impact of defect distribution and agglomeration is so far
underrepresented. In this study, we have used $ab\ initio$ based simulation to
close this knowledge gap.

We have analyzed a range of defect concentrations and strengths of defect
dipoles. Thereby we found a cross-over between non-polar or weak defect
dipoles, which mainly induce local depolarization fields, and strong defect
dipoles larger than about 30% of the saturation polarization of the tetragonal
phase, which induce internal electrical fields. Both types of defects reduce
thermal hysteresis and coercive field and affect the phase diagram by less
than 1% defect. We have shown that, both for phase stability and field
hysteresis, the spatial distribution of defects and their agglomeration are
thereby more important than their exact concentration or strength.
Specifically, depolarization and internal fields are associated with the
active surface area of the defect configuration. Thus, randomly distributed
defects have e.g. the largest impact on the paraelectric to ferroelectric
phase transition. Agglomeration or clustering reduces this active surface area
and thus the impact of the defects on the material. Thus, large changes of
functional properties are to be expected for the formation or decomposition of
defect clusters with time. Most interesting are defect agglomerates in two
dimensions, which reduce the symmetry of the system. In this case, independent
of the defect type, i.e. the strength of the defect dipole, particular
polarization directions are favored by the defects. Thus phase stability and
field hysteresis depend on the relative direction of defect agglomerate and
external field. Pinched hysteresis loops and large recoverable stored energy
are possible for one direction, while the defects are rather invisible in the
other direction. The control of defect dipole strengths (e.g., by using
different dopants), agglomeration (e.g., by ion bombardment), as well as the
prevention of ion migration (e.g., by heat treatment and interfaces) are thus
promising routes to tap the full potential of defect engineering in
ferroelectric materials. One promising route to tailor material properties by
defects may thus be the design of 2-dimensional defect-rich areas.
Furthermore, it was interesting to further explore if there are atomistic
realizations of defect dipoles close to the cross-over strength, which could
be utilized to introduce multi-critical points in the phase diagram.

## Acknowledgments

Funding: This work was supported by the German research foundation (DFG) GR
4792/2. We thank Prof. Dr. Markus Stricker for his helpful comments and
fruitful discussions.

## Declaration of generative AI in scientific writing

During the preparation of this work, the authors used DeepL and Grammarly for
wording suggestions and checking of language. After using these tools, we
reviewed and edited the content as needed and take full responsibility for the
content of the publication.

## References

  * [1] R. W. Whatmore, Y.-M. You, R.-G. Xiong, C.-B. Eom, 100 years of ferroelectricity—A celebration, APL Mater. 9 (7) (2021) 070401.  doi:10.1063/5.0059208.
URL https://aip.scitation.org/doi/full/10.1063/5.0059208

  * [2] A. Grünebohm, M. Marathe, R. Khachaturyan, R. Schiedung, D. C. Lupascu, V. V. Shvartsman, Interplay of domain structure and phase transitions: Theory, experiment and functionality, J. Phys.: Condens. Matter 34 (7) (2021) 073002.  doi:10.1088/1361-648X/ac3607.
URL https://doi.org/10.1088/1361-648x/ac3607

  * [3] V. Veerapandiyan, F. Benes, T. Gindel, M. Deluca, Strategies to Improve the Energy Storage Properties of Perovskite Lead-Free Relaxor Ferroelectrics: A Review, Materials 13 (24) (2020) 5742.  doi:10.3390/ma13245742.
URL https://www.mdpi.com/1996-1944/13/24/5742

  * [4] R. Khosla, S. K. Sharma, Integration of Ferroelectric Materials: An Ultimate Solution for Next-Generation Computing and Storage Devices, ACS Appl. Electron. Mater. 3 (7) (2021) 2862–2897.  doi:10.1021/acsaelm.0c00851.
URL https://doi.org/10.1021/acsaelm.0c00851

  * [5] V. F. Michel, T. Esswein, N. A. Spaldin, Interplay between ferroelectricity and metallicity in BaTiO3, J. Mater. Chem. C 9 (27) (2021) 8640–8649.  doi:10.1039/D1TC01868J.
URL https://pubs.rsc.org/en/content/articlelanding/2021/tc/d1tc01868j

  * [6] V. Veerapandiyan, M. N. Popov, F. Mayer, J. Spitaler, S. Svirskas, V. Kalendra, J. Lins, G. Canu, M. T. Buscaglia, M. Pasciak, J. Banys, P. B. Groszewicz, V. Buscaglia, J. Hlinka, M. Deluca, Origin of Relaxor Behavior in Barium-Titanate-Based Lead-Free Perovskites, Adv. Electron. Mater. 8 (2) (2022) 2100812.  doi:10.1002/aelm.202100812.
URL https://onlinelibrary.wiley.com/doi/abs/10.1002/aelm.202100812

  * [7] M. Acosta, N. Novak, V. Rojas, S. Patel, R. Vaish, J. Koruza, G. A. Rossetti, J. Rödel, BaTiO3-based piezoelectrics: Fundamentals, current status, and perspectives, Appl. Phys. Rev. 4 (4) (2017) 041305.  doi:10.1063/1.4990046.
URL https://aip.scitation.org/doi/10.1063/1.4990046

  * [8] A. Dimou, A. Biswas, A. Gruenebohm, Ab initio based study on atomic ordering in (Ba, Sr)TiO3 (Aug. 2023).  arXiv:2308.13263, doi:10.48550/arXiv.2308.13263.
URL http://arxiv.org/abs/2308.13263

  * [9] N. A. Hill, Why Are There so Few Magnetic Ferroelectrics?, J. Phys. Chem. B 104 (29) (2000) 6694–6709.  doi:10.1021/jp000114x.
URL https://doi.org/10.1021/jp000114x

  * [10] Y.-B. Ma, C. Molin, V. V. Shvartsman, S. Gebhardt, D. C. Lupascu, K. Albe, B.-X. Xu, State transition and electrocaloric effect of BaZrxTi1-xO3: Simulation and experiment, Journal of Applied Physics 121 (2) (2017) 024103\.  doi:10.1063/1.4973574.
URL https://doi.org/10.1063/1.4973574

  * [11] Y. W. Li, J. F. Scott, D. N. Fang, F. X. Li, 90-degree polarization switching in BaTiO3 crystals without domain wall motion, Appl. Phys. Lett. 103 (23) (2013) 232901.  doi:10.1063/1.4832784.
URL https://aip.scitation.org/doi/10.1063/1.4832784

  * [12] H.-Y. Huang, M. Wu, L.-J. Qiao, The behaviour of 180∘ polarization switching in BaTiO3 from first principles calculations, Comp. Mater. Sci. 82 (2014) 1–4.  doi:10.1016/j.commatsci.2013.09.013.
URL https://www.sciencedirect.com/science/article/pii/S0927025613005375

  * [13] Y. A. Genenko, J. Glaum, M. J. Hoffmann, K. Albe, Mechanisms of aging and fatigue in ferroelectrics, Mater. Sci. Eng. B 192 (2015) 52–82.  doi:10.1016/j.mseb.2014.10.003.
URL https://www.sciencedirect.com/science/article/pii/S0921510714002189

  * [14] M. Krupska-Klimczak, P. Gwizd, I. Jankowska-Sumara, D. Sitko, P. Jeleń, An Overview of Some Nonpiezoelectric Properties of BaTiO3 Ceramics Doped by Eu Ions, Materials 15 (15) (2022) 5363.  doi:10.3390/ma15155363.
URL https://www.mdpi.com/1996-1944/15/15/5363

  * [15] Z. Li, J. Yu, S. Hao, P.-E. Janolin, Enhancing properties of lead-free ferroelectric BaTiO3 through doping, J. Eur. Ceram. Soc. 42 (12) (2022) 4693–4701.  doi:10.1016/j.jeurceramsoc.2022.05.023.
URL https://www.sciencedirect.com/science/article/pii/S0955221922003673

  * [16] X. Ren, Large electric-field-induced strain in ferroelectric crystals by point-defect-mediated reversible domain switching, Nat. Mater. 3 (2) (2004) 91–94.  doi:10.1038/nmat1051.
URL https://www.nature.com/articles/nmat1051

  * [17] J. Li, J. Lv, D. Zhang, L. Zhang, X. Hao, M. Wu, B.-X. Xu, M. Otonicar, T. Lookman, B. Dkhil, X. Lou, Doping-induced Polar Defects Improve the Electrocaloric Performance of Ba0.9Sr0.1Hf0.1Ti0.9O3, Phys. Rev. Appl. 16 (1) (2021) 014033.  doi:10.1103/PhysRevApplied.16.014033.
URL https://link.aps.org/doi/10.1103/PhysRevApplied.16.014033

  * [18] A. Grünebohm, T. Nishimatsu, Influence of defects on ferroelectric and electrocaloric properties of BaTiO3, Phys. Rev. B 93 (13) (2016) 134101.  doi:10.1103/PhysRevB.93.134101.
URL https://link.aps.org/doi/10.1103/PhysRevB.93.134101

  * [19] L. Zhang, E. Erdem, X. Ren, R.-A. Eichel, Reorientation of (MnTi${{}^{\prime\prime}}-$VO$\bullet\bullet$)× defect dipoles in acceptor-modified BaTiO3 single crystals: An electron paramagnetic resonance study, Appl. Phys. Lett. 93 (20) (2008) 202901.  doi:10.1063/1.3006327.
URL https://aip.scitation.org/doi/10.1063/1.3006327

  * [20] P. Erhart, P. Träskelin, K. Albe, Formation and switching of defect dipoles in acceptor-doped lead titanate: A kinetic model based on first-principles calculations, Phys. Rev. B 88 (2) (2013) 024107\.  doi:10.1103/PhysRevB.88.024107.
URL https://link.aps.org/doi/10.1103/PhysRevB.88.024107

  * [21] Z. Zhang, Z. Wen, T. Li, Z. Wang, Z. Liu, X. Liao, S. Ke, L. Shu, Flexoelectric aging effect in ferroelectric materials, J. Appl. Phys. 133 (5) (2023) 054102\.  doi:10.1063/5.0134531.
URL https://aip.scitation.org/doi/full/10.1063/5.0134531

  * [22] C. Verdier, D. C. Lupascu, H. von Seggern, J. Rödel, Effect of thermal annealing on switching dynamics of fatigued bulk lead zirconate titanate, Appl. Phys. Lett. 85 (15) (2004) 3211–3213.  doi:10.1063/1.1801679.
URL https://doi.org/10.1063/1.1801679

  * [23] D. C. Lupascu, Agglomeration and Microstructural Effects, in: D. C. Lupascu (Ed.), Fatigue in Ferroelectric Ceramics and Related Issues, Springer Series in Materials Science, Springer, Berlin, Heidelberg, 2004, pp. 63–79.  doi:10.1007/978-3-662-07189-2_3.
URL https://doi.org/10.1007/978-3-662-07189-2_3

  * [24] S. Saremi, R. Xu, F. I. Allen, J. Maher, J. C. Agar, R. Gao, P. Hosemann, L. W. Martin, Local control of defects and switching properties in ferroelectric thin films, Phys. Rev. Mater. 2 (8) (2018) 084414.  doi:10.1103/PhysRevMaterials.2.084414.
URL https://link.aps.org/doi/10.1103/PhysRevMaterials.2.084414

  * [25] T. Suzuki, M. Ueno, Y. Nishi, M. Fujimoto, Dislocation Loop Formation in Nonstoichiometric (Ba,Ca)TiO3 and BaTiO3 Ceramics, J. Am. Ceram. Soc. 84 (1) (2001) 200–206.  doi:10.1111/j.1151-2916.2001.tb00631.x.
URL
https://ceramics.onlinelibrary.wiley.com/doi/abs/10.1111/j.1151-2916.2001.tb00631.x

  * [26] S. M. Yang, T. H. Kim, J.-G. Yoon, T. W. Noh, Nanoscale Observation of Time-Dependent Domain Wall Pinning as the Origin of Polarization Fatigue, Adv. Funct. Mater. 22 (11) (2012) 2310–2317.  doi:10.1002/adfm.201102685.
URL https://onlinelibrary.wiley.com/doi/abs/10.1002/adfm.201102685

  * [27] D. Akbarian, D. E. Yilmaz, Y. Cao, P. Ganesh, I. Dabo, J. Munro, R. V. Ginhoven, A. C. T. van Duin, Understanding the influence of defects and surface chemistry on ferroelectric switching: A ReaxFF investigation of BaTiO3, Phys. Chem. Chem. Phys. 21 (33) (2019) 18240–18249.  doi:10.1039/C9CP02955A.
URL https://pubs.rsc.org/en/content/articlelanding/2019/cp/c9cp02955a

  * [28] K. K. Adepalli, J. Yang, J. Maier, H. L. Tuller, B. Yildiz, Tunable Oxygen Diffusion and Electronic Conduction in SrTiO3 by Dislocation-Induced Space Charge Fields, Advanced Functional Materials 27 (22) (2017) 1700243.  doi:10.1002/adfm.201700243.
URL https://onlinelibrary.wiley.com/doi/abs/10.1002/adfm.201700243

  * [29] W. Zhong, D. Vanderbilt, K. M. Rabe, First-principles theory of ferroelectric phase transitions for perovskites: The case of BaTiO3, Phys. Rev. B 52 (9) (1995) 6301–6312.  doi:10.1103/PhysRevB.52.6301.
URL https://link.aps.org/doi/10.1103/PhysRevB.52.6301

  * [30] T. Nishimatsu, A. Grünebohm, U. V. Waghmare, M. Kubo, Molecular Dynamics Simulations of Chemically Disordered Ferroelectric (Ba,Sr)TiO3 with a Semi-Empirical Effective Hamiltonian, J. Phys. Soc. Jpn. 85 (11) (2016) 114714.  doi:10.7566/JPSJ.85.114714.
URL https://journals.jps.jp/doi/10.7566/JPSJ.85.114714

  * [31] S.-H. Teng, Impact of defects on the phase diagram and field response of ferroelectric (Ba,Sr)TiO3 (Mar. 2024).  doi:10.5281/zenodo.8418579.
URL https://zenodo.org/records/8418579

  * [32] S.-H. Teng, Sheng-Han Teng / p1-defect_study $\cdot$ GitLab (Mar. 2024).
URL https://gitlab.ruhr-uni-bochum.de/tengssh/p1-defect_study

  * [33] S. Liu, R. E. Cohen, Multiscale simulations of defect dipole–enhanced electromechanical coupling at dilute defect concentrations, Appl. Phys. Lett. 111 (8) (2017) 082903.  doi:10.1063/1.4989670.
URL https://aip.scitation.org/doi/10.1063/1.4989670

  * [34] R. Machado, A. Di Loreto, A. Frattini, M. Sepliarsky, M. G. Stachiotti, Site occupancy effects of Mg impurities in BaTiO3, J. Alloy. Compd. 809 (2019) 151847.  doi:10.1016/j.jallcom.2019.151847.
URL https://www.sciencedirect.com/science/article/pii/S0925838819330865

  * [35] Y. Zhou, D. Yang, G. Jiang, Effect of Yttrium Doping on Microstructure and Temperature Resistance Characteristics of Barium Strontium Titanate Ferroelectric Ceramics, in: 2022 IEEE International Conference on High Voltage Engineering and Applications (ICHVE), 2022, pp. 1–4.  doi:10.1109/ICHVE53725.2022.9961754.
  * [36] B. Vigneshwaran, P. Kuppusami, S. Ajithkumar, H. Sreemoolanadhan, Study of low temperature-dependent structural, dielectric, and ferroelectric properties of BaxSr(1-x)TiO3 (x = 0.5, 0.6, 0.7) ceramics, J Mater Sci: Mater Electron 31 (13) (2020) 10446–10459.  doi:10.1007/s10854-020-03593-3.
URL https://doi.org/10.1007/s10854-020-03593-3

  * [37] J. W. Liou, B. S. Chiou, Dielectric characteristics of doped Ba1-xSrxTiO3 at the paraelectric state, Mater. Chem. Phys. 51 (1) (1997) 59–63.  doi:10.1016/S0254-0584(97)80267-2.
URL https://www.sciencedirect.com/science/article/pii/S0254058497802672

  * [38] V. V. Lemanov, E. P. Smirnova, P. P. Syrnikov, E. A. Tarakanov, Phase transitions and glasslike behavior in Sr1-xBaxTiO3, Phys. Rev. B 54 (5) (1996) 3151–3157.  doi:10.1103/PhysRevB.54.3151.
URL https://link.aps.org/doi/10.1103/PhysRevB.54.3151

  * [39] B. Noheda, D. E. Cox, G. Shirane, S.-E. Park, L. E. Cross, Z. Zhong, Polarization Rotation via a Monoclinic Phase in the Piezoelectric 92% PbZn1/3Nb2/3O3-8% PbTiO3$, Phys. Rev. Lett. 86 (17) (2001) 3891–3894.  doi:10.1103/PhysRevLett.86.3891.
URL https://link.aps.org/doi/10.1103/PhysRevLett.86.3891

  * [40] M. Marathe, A. Grünebohm, T. Nishimatsu, P. Entel, C. Ederer, First-principles-based calculation of the electrocaloric effect in BaTiO3: A comparison of direct and indirect methods, Phys. Rev. B 93 (5) (2016) 054110\.  doi:10.1103/PhysRevB.93.054110.
URL https://link.aps.org/doi/10.1103/PhysRevB.93.054110

  * [41] M. Marathe, D. Renggli, M. Sanlialp, M. O. Karabasov, V. V. Shvartsman, D. C. Lupascu, A. Grünebohm, C. Ederer, Electrocaloric effect in BaTiO3 at all three ferroelectric transitions: Anisotropy and inverse caloric effects, Phys. Rev. B 96 (1) (2017) 014102.  doi:10.1103/PhysRevB.96.014102.
URL https://link.aps.org/doi/10.1103/PhysRevB.96.014102

  * [42] N. Novak, R. Pirc, Z. Kutnjak, Impact of critical point on piezoelectric and electrocaloric response in barium titanate, Phys. Rev. B 87 (10) (2013) 104102.  doi:10.1103/PhysRevB.87.104102.
URL https://link.aps.org/doi/10.1103/PhysRevB.87.104102

  * [43] B. Luo, X. Wang, E. Tian, L. Wu, L. Li, First-principles effective Hamiltonian simulation of ABO3-type perovskite ferroelectrics for energy storage application, J. Appl. Phys. 120 (7) (2016) 074106.  doi:10.1063/1.4961204.
URL https://aip.scitation.org/doi/10.1063/1.4961204

##  Appendix A Pristine material

The phase diagram of Ba0.7Sr0.3TiO3 as found for the pristine case (without
defect dipoles) is shown in Fig. 10. Thereby, subfigures (a) and (b) show the
change of the macroscopic polarization components and of the lattice
constants, respectively.

[abs] (a)[abs] (b)

Figure 10: Temperature dependence of (a) polarization components
($|{\left\langle{P_{i}}\right\rangle}|$, where $i:x,y,z$), and (b) lattice
constants ($a,b,c$ correspond to the length along $x,y,z$) of Ba0.7Sr0.3TiO3
without defects. Colors encode cooling (blue) and heating (red) simulations.

##  Appendix B Phase diagrams

Figure 11 compares phase diagrams in the presence of (a) polar defect dipoles
in the plane ($|u_{D}|=0.1$ Å$\hat{x}$) and (b) non-polar defects ($|u_{D}|=0$
Å) for 2D agglomerates (solid lines) and random distribution in 2D (dashed
lines).

[abs] * * * *576115223044608(a)* * * *576115223044608(b)

Figure 11: Temperature - defect concentration phase diagrams in the presence
of (a) polar defect dipoles ($|u_{D}|=$0.1 Å$\hat{x}$) and (b) non-polar
defects ($|u_{D}|=$0.0 Å) for stripes of defects (solid lines) and random
distribution in 2D (dashed lines). Red and blue lines show transition
temperatures for heating and cooling, respectively. Numbers give the
approximate number of active surface units in R phases. Note that the defect
agglomerates have no effective active surface area for T and O phases as
discussed in the text.

##  Appendix C Macroscopic changes of polarization and strain by non-polar
defects

Figure 12 (a) and (b) show the change of macroscopic strain components
($\sigma_{ii}$ with $i=x,y,z$) for random distribution in 3D and 2D
agglomerates with varying concentration or width. Figure 13 (a) and (b) show
the change in macroscopic polarization and strain for cubic agglomerates.
Figure 14 (a) and (b) show the impact of defect dipole strength on the change
of macroscopic polarization for 1 % defects of 2D random distribution and 2D
agglomerate.

[abs] (a)[abs] (b)

Figure 12: Changes of strain components $\sigma_{ii}$, with $i:x,y,z$, in
cooling simulations induced by non-polar defects ($|u_{D}|=$0.0 Å) for (a)
random distribution and (b) stripes in $z$-planes with varying width ($w$).

[abs] (a)[abs] (b)

Figure 13: Changes of (a) polarization components
$|{\left\langle{P_{i}}\right\rangle}|$ and (b) strain components
$\sigma_{ii}$, with $i:x,y,z$, in cooling simulations induced by non-polar
defects ($|u_{D}|=$0.0 Å) with varying length ($a$) of cube. Note that for the
cube with $a=13$ u.c., depolarization fields induce multidomains, resulting in
a continuous O–R phase transition.

[abs] (a)[abs] (b)

Figure 14: Impact of the strength of 1 % defects ($|u_{D}|$ $\hat{z}$) on the
changes of polarization components $|{\left\langle{P_{i}}\right\rangle}|$,
with $i:x,y,z$, in cooling simulations for (a) random distribution and (b)
stripes. Colors encode the strength of the defect dipoles from 0.00 (red),
0.02 (blue), 0.03 (green), 0.04 (light blue), 0.06 (light green) to 0.10
(purple).

##  Appendix D Microscopic changes of polarization by defects

Figure 15 (a) and (b) show the microscopic polarization for a 2D agglomerate
of $w=12$ u.c. and a cubic agglomerate of $a=13$ u.c., showing the formation
of R71 domain walls above/below the agglomerates.

\begin{overpic}[abs,width=195.12767pt]{figures/phase_diagram/agg2d/dipos_ag_12_c0_100K_y40.png}
\put(0.0,123.0){\text{(a)}} \end{overpic}

\begin{overpic}[abs,width=195.12767pt]{figures/phase_diagram/cube/dipos_cc13_c0_100K_z24_in.png}
\put(0.0,123.0){\text{(b)}} \end{overpic}

Figure 15: Cross section of time-averaged local dipoles in the presence of (a)
a 2D agglomerate with 0.5 % ($w=12$ u.c.) non-polar defects at 100 K (R phase)
and (b) a cube with 2 % ($a=13$ u.c.) non-polar defects at 100 K, showing R71
domain walls. Colors indicate the sign of $P_{z}$s.

##  Appendix E Field hysteresis and snapshots

Figure 16 (a) and (b) show the change in macroscopic polarization and strain
for the pristine to system with non-polar defects in cubic, 2D agglomerate and
2D random distributions with the bipolar field $E_{x}$ applied. Figure 17 (a),
(b) and (c) show the strain–field ($\sigma$-E) curves for non-polar defects
with the electric field $E_{z}$ and polar defects with the electric field
$E_{x}$ or $E_{z}$ applied. Figure 18 (a) and (b) show all 3 polarization
components for non-polar defects with bipolar field $E_{z}$ applied, while
Fig. 18 (a) and (b) are for polar defects ($|u_{D}|=0.1$ Å$\hat{z}$) with the
bipolar electric field $E_{x}$ perpendicular to the defect dipoles. Figure 19
shows snapshots of systems with 1 % defects of cubic agglomerate, 3D random
distribution, 2D agglomerate and 2D random distribution at $E_{z}=-150$ kV/cm,
showing pinned domains near defects that could act as nucleation centers for
the macroscopic switching towards the positive polarization direction.

[abs]

$P_{i}$ ($\mu C/cm^{2}$)

$E_{x}$ ($kV/cm$)(a)[abs]

$\sigma_{ii}$ (%)

$E_{x}$ ($kV/cm$)(b)

Figure 16: Changes of the field hysteresis at $T=250$ K if going from the
pristine material (red) to systems with different distribution of 1 % defects
(green: cube; blue and cyan: stripes and randomly in $\hat{z}$-planes) for
non-polar defects for the electric field perpendicular ($E_{x}$) to the defect
dipoles. Only the polarization component $P_{i}$ and strain component
$\sigma_{ii}$ along the field direction are shown.

[abs]

$\sigma_{ii}$ (%)

$E_{z}$ ($kV/cm$)(a)[abs]  $E_{x}$ ($kV/cm$)(b)[abs]  $E_{z}$ ($kV/cm$)(c)

Figure 17: Changes of the field hysteresis at $T=250$ K if going from the
pristine material (red) to systems with different distribution of 1 % defects
(green: cube; blue and cyan: stripes and randomly in $\hat{z}$-planes; purple:
randomly) for (a) non-polar defects and (b)–(c) polar defect dipoles
($|u_{D}|=$0.1 Å$\hat{z}$) for the electric field (b) perpendicular ($E_{x}$)
and (c) collinear ($E_{z}$) to the defect dipoles. Only the diagonal strain
component $\sigma_{ii}$ along the field direction is shown.

[abs]

$P_{i}$ ($\mu C/cm^{2}$)

$E_{z}$ ($kV/cm$). $P_{x}$$P_{y}$$P_{z}$(a)[abs]

$P_{i}$ ($\mu C/cm^{2}$)

$E_{z}$ ($kV/cm$). $P_{x}$$P_{y}$$P_{z}$(b)

[abs]

$P_{i}$ ($\mu C/cm^{2}$)

$E_{x}$ ($kV/cm$). $P_{z}$$P_{y}$$P_{x}$(c)[abs]

$P_{i}$ ($\mu C/cm^{2}$)

$E_{x}$ ($kV/cm$). $P_{z}$$P_{y}$$P_{x}$(d)

Figure 18: Field dependence of polarization component $P_{i}$ with $i=x,y,z$
at $T=250$ K with 1 % defects for (a,c) stripes and (b,d) randomly in one
plane for (a,b) non-polar defects and (c,d) polar defects with ($|u_{D}|=$0.1
Å$\hat{z}$) for the electric field (a,b) along the $z$ or (c,d) perpendicular
($E_{x}$) to the defect dipoles.

[abs] (a)[abs] (b)

[abs] (c)[abs] (d)

Figure 19: Snapshots of dipoles at 250 K for 1 % defects ($|u_{D}|=$0.1
Å$\hat{z}$) in (a) cubic agglomerate of $a=10$ u.c., (b) 3D random, (c) 2D
agglomerate of $w=24$ u.c. and (d) 2D random in a plane with the applied field
$-150$ kV/cm $\hat{z}$. Colors indicate the sign of $P_{z}$s. Note that all
these snapshots are not time-averaged and are thus noisy.