PHYSICAL REVIEW B
103
, 144201 (2021)
Understanding the morphotropic phase boundary of perovskite solid solutions as a frustrated state
Ying Shi Teh
,
1
Jiangyu Li,
2
and Kaushik Bhattacharya
1
,
*
1
Division of Engineering and Applied Science, California Institute of Technology, Pasadena, California 91125, USA
2
Department of Materials Science and Engineering, Southern University of Science and Technology, Shenzhen 518055, Guangdong, China
(Received 18 December 2020; revised 7 March 2021; accepted 23 March 2021; published 5 April 2021)
This paper presents a lattice model incorporating random fields and long-range interactions where a frustrated
state emerges at a specific composition but is suppressed elsewhere. The model is motivated by perovskite solid
solutions and explains the phase diagram in such materials including the morphotropic phase boundary (MPB)
that plays a critical role in applications for its enhanced dielectric, piezoelectric, and optical properties. Further,
the model also suggests the possibility of phenomena by exploiting MPBs.
DOI:
10.1103/PhysRevB.103.144201
I. INTRODUCTION
Perovskites are materials with a chemical composition
ABO
3
, where A and B are typically transition metals that have
a crystal structure similar to that of the mineral perovskite,
CaTiO
3
(
Pm
3
m
), Fig.
1(a)
(e.g., Refs. [
1
,
2
]). Some of these
materials undergo a displacive phase transition on cooling
from the high temperature perovskite structure to a low tem-
perature low symmetry structure that is not centrosymmetric
and thus displays ferroelectricity or ferromagnetism. There-
fore these materials are widely used in capacitor, ultrasonic,
optical, sensor, and actuator applications for their dielectric
and piezoelectric properties.
The basic structure is extremely stable, and it is possible to
have solid solutions A(C
x
D
1
−
x
)O
3
where two metallic species
C and D substitutionally occupy the B site of the lattice.
The low temperature structure in these compounds depends
on composition and changes across the
morphotropic phase
boundary
(MPB) which is largely independent of temperature.
Figure
1(b)
shows the phase diagram of the widely used
piezoelectric lead zirconate titanate (PbZr
x
Ti
1
−
x
O
3
or PZT)
that has a MPB at
x
=
0
.
52 with a ferroelectric rhombohedral
(
R
3
m
) structure in the Zr-rich compositions
1
and ferroelectric
tetragonal (
P
4
mm
) structure in the Ti-rich compositions. The
dielectric and piezoelectric properties as well as the ability to
pole a ceramic increase dramatically at the MPB [
3
–
5
], and
therefore PZT is widely used close to this composition. The
search for lead-free dielectric and piezoelectric materials has
also focused on solid solutions with MPBs (e.g., Ref. [
6
]).
The nature of the MPB has been extensively studied since
the discovery of a bridging monoclinic phase (
Cm
)atthe
MPB by Noheda
et al.
[
7
] using x-ray powder diffraction.
It has been suggested that the presence of a bridging phase
enables a larger intrinsic [
4
,
8
,
9
] and extrinsic [
5
] piezoelectric
effect at the MPB. Since then, various structures have been
*
Corresponding author: bhatta@caltech.edu
1
PZT shows a second rhombohedral (
R
3
c
) phase at low tempera-
ture at high Zr compositions, but we focus on the compositions near
the MPB.
observed: the combination of two monoclinic phases (
Cm
and
Ic
[
10
]or
Cm
and
Pm
[
11
]), combination of tetragonal
(
P
4
mm
) and monoclinic (
Cm
)[
12
], and combination of rhom-
bohedral (
R
3
m
) and monoclinic (
Cm
)[
13
]. This uncertainty
has been attributed to the disorder in the composition which
further results in a disorder in the structure and the difficulty
of resolving local structures [
14
]. This role of disorder is also
supported by first principles calculations [
15
,
16
]. This, how-
ever, raises the question as to why the disorder does not affect
the structure away from the MPB. Another interesting obser-
vation concerns the domain patterns. Classical well-defined
domain patterns are observed away from the MPB, but highly
fragmented domain patterns are observed near the MPB [
17
].
It has been argued that this fragmented domain pattern also
contributes to the high piezoelectric response [
18
].
In short, critical questions remain open. Why is the effect
of compositional disorder suppressed to form an unambiguous
structure away from the MPB but suddenly revealed at the
MPB? Is there a definitive crystal structure at the MPB? Why
do domain patterns become fragmented near the MPB? Does
compositional disorder play a role in the ease of poling at the
MPB? Can the MPB be exploited to create new phenomena?
These questions have proved to be challenging. The disor-
dered nature of the solid solution requires a large ensemble
that takes it beyond the scope of first principles calculations
without the introduction of a composition-dependent hybrid
pseudopotential [
19
] that does not include randomness. On
the other hand, phase-field Landau-Ginzburg methods can
provide insight into domain patterns but require a phenomeno-
logical model of the phase transition.
In this paper, we propose a model based on the random-
field Ising model with long-range interactions that incorpo-
rates the basic elements of the underlying physics. The B
sites of a perovskite form a reference cubic lattice that is
occupied randomly by atoms of either C or D species. The
local quantum mechanical interactions create a propensity for
the unit cell to break cubic symmetry depending on the species
at the B site. The ferroelectric, ferromagnetic, and ferroelastic
polarizations lead to long-range interactions. We create an
effective Hamiltonian with these physics.
2469-9950/2021/103(14)/144201(7)
144201-1
©2021 American Physical Society
TEH, LI, AND BHATTACHARYA
PHYSICAL REVIEW B
103
, 144201 (2021)
(
)
(
)(
)(
)
MPB
A
B
O
FIG. 1. (a) Perovskite structure. (b) Phase diagram of PZT
adapted from Cross [
3
]. F
R
,F
T
,andP
C
denote ferroelectric
rhombohedral, tetragonal, and paraelectric phases, respectively.
(c) Schematic illustration of the lattice model proposed in this work.
In the first part of the paper, we show that our model
provides insights to the questions concerning the MPB of
ferroelectric solid solutions like PZT. In particular, the long-
range interactions overwhelm the local disorder in the C-rich
and D-rich compositions with the exchange of stability tak-
ing place at a specific composition where the material is
frustrated. This frustration manifests itself as the MPB. The
frustrated state also enables easy poling as observed. In the
second part, we use the model to explore the possibility of ob-
taining a strongly coupled multiferroic material system using
the insights obtained from the first part.
II. MPB IN FERROELECTRIC MATERIALS
A. Model
Consider a
d
-dimensional periodic lattice (
d
=
2 or 3) with
N
lattice points as shown in Fig.
1(c)
. Each lattice point
i
is characterized by fixed (quenched) chemical composition
(
c
i
) of either type C [
c
i
=
0 indicated by a red open circle
in Fig.
1(c)
]ortypeD[
c
i
=
1 indicated by a blue closed
circle]. Each lattice point carries a dipole state [
p
i
indicated
by the arrows in Fig.
1(c)
] that can take one of a number of
orientations
C
∪
D
determined by the Hamiltonian
W
(
{
p
i
}
;
{
c
i
}
)
=
N
i
=
1
h
i
−
1
2
<
i
,
j
>
J
e
ij
p
i
·
p
j
+
W
e
.
(1)
The first term encodes the information that lattice site of type
C (respectively, D) energetically prefers the set of dipole states
C
indicated by the red arrows (respectively,
D
indicated by the
blue arrows), though they can take states in
D
(respectively,
C
)
with an energetic cost
h
C
D
>
0 (respectively,
h
D
C
>
0):
h
i
=
⎧
⎨
⎩
h
C
D
c
i
=
0&
p
i
∈
D
,
h
D
C
c
i
=
1&
p
i
∈
C
,
0
otherwise.
(2)
The second term, where the sum is limited to nearest
neighbors, is the exchange energy. Ferroelectricity requires
noncentrosymmetric displacements of ions which arise from
a delicate balance between short-range repulsions of electron
clouds and a short-range portion of the Coulomb interac-
tion with a range of the lattice constant [
20
]. The exchange
term captures this net ferroelectric effect as it promotes
like neighbors when the exchange constant
J
e
ij
>
0. The ex-
change constant in this model is composition dependent and
takes the form of
J
e
ij
=
1
2
[
J
e
(
c
i
)
+
J
e
(
c
j
)]
,
J
e
(
c
)
=
J
C
e
c
=
0
,
J
D
e
c
=
1
.
(3)
The third term represents the electrostatic contribution that
includes the long-range dipole-dipole interaction scaled by
D
e
(which incorporates the dipole strength, lattice constant, and
electromagnetic constants) and the influence of the applied
external electric field
E
.
W
e
=
W
em
(
{
p
i
}
,
D
e
,
E
)
:
=
D
e
1
(
d
−
1)
N
i
,
j
=
1
,
R
1
r
d
ij
p
i
·
p
j
−
d
(
p
i
·
r
ij
)(
p
j
·
r
ij
)
r
2
ij
+
2
π
d
N
i
=
1
|
p
i
|
2
−
E
·
N
i
p
i
.
(4)
Given a lattice where the composition of each site is ran-
domly assigned subject to a fixed average, we use a Markov
chain Monte-Carlo (MCMC) method with cooling to obtain
the equilibrium distribution at a given temperature. The state
is initialized by randomly assigning a polarization from
C
∪
D
. Adapting the Metropolis-Hastings algorithm to our mul-
tistate setting, a site is chosen at random and its dipole state
is updated to one of the
N
s
states according to the transition
probability
P
s
=
exp(
−
β
W
(
s
)
)
N
s
r
=
1
exp(
−
β
W
(
r
)
)
,
s
=
1
, ...,
N
s
,
(5)
where
β
is the inverse temperature and
N
s
is the cardinal-
ity of
C
∪
D
. We avoid the system getting trapped in local
minima at low temperatures by starting at a high temperature
(
β
=
0) and slowly cooling (increasing
β
) to the temperature
of interest, while performing enough MCMC steps to reach
equilibrium at each temperature.
Note that the dipole-dipole term is conditionally conver-
gent. So, we adopt Ewald summation [
21
–
23
] to separate
it into a short-range contribution that converges rapidly in
real space, a long-range contribution that converges rapidly in
Fourier space, a self energy, and a surface term that depends
on the boundary condition. Further, notice that for any flip,
W
e
≈−
E
i
·
p
i
where
E
i
=−∇
p
i
W
e
is the electric field.
While
E
i
has to be recomputed after each flip, the error is
small for individual flips. Therefore we update
E
i
only every
√
N
flips where
N
is the size of the lattice. We can then per-
form
√
N
flips independently and in parallel thereby enabling
acceleration on a graphical processing unit (GPU). Further,
E
i
is readily computed using fast Fourier transforms which can
also be implemented on GPUs. Further details are provided in
Supplemental Material [
24
].
B. Results
We study an example motivated by PZT though the
results are generic. Here, the C and D lattice points rep-
resent lead zirconate (PZ) and lead titanate (PT) unit
cells containing Zr and Ti atoms, respectively. Recall
that the former prefers rhombohedral or
111
polariza-
tion states while the latter prefers tetragonal or
100
144201-2
UNDERSTANDING THE MORPHOTROPIC PHASE BOUNDARY ...
PHYSICAL REVIEW B
103
, 144201 (2021)
FIG. 2. Emergence of a morphotropic phase transition (MPB)
as a competition between short-range (compositional) disorder and
long-range (exchange and electrostatic) interactions. (a) Order pa-
rameter vs inverse temperature for various compositions. (b) Dipole
orientation in the ordered phase at various compositions. (c) Phase
diagram showing the MPB. (d) Domain patterns at various com-
position. (e) Experimentally observed domain patterns at various
compositions (reprinted with permission from Woodward, Knudsen,
and Reaney [
17
]). (f) Average domain size vs composition [average
(red) and ten realizations (gray)]. (g) Small-scale oscillations in
terms of the finest Haar wavelet horizontal (H), vertical (V), and
diagonal (D) detail coefficients vs composition.
polarization states. We consider two dimensions
d
=
2 so that
C
=
(
p
C
o
/
√
2)
{
[1
,
1]
,
[1
,
−
1]
,
[
−
1
,
1]
,
[
−
1
,
−
1]
}
while
D
=
p
D
o
{
[1
,
0]
,
[0
,
1]
,
[
−
1
,
0]
,
[0
,
−
1]
}
. Unless otherwise speci-
fied, we set
p
C
0
=
p
D
0
=
h
C
D
=
h
D
C
=
J
C
e
=
J
D
e
=
D
e
=
1,
unit lattice distance between any two nearest neighboring
sites, and conduct our simulations on a 256
2
lattice. The
lattice size is chosen to be sufficiently large upon convergence
studies [
24
] to depict randomness in composition as well as
for us to visualize domain formation. In each simulation, we
begin with an inverse temperature of
β
=
0 and repeatedly
increase its value with a small step size of
β
=
0
.
05 until we
reach
β
=
5. At each temperature value, at least 2000 Monte
Carlo (MC) sweeps (each sweep consists of
N
=
256
2
steps)
are performed with a total of approximately 2
×
10
5
sweeps.
Figure
2
shows the results of these simulations. Fig-
ure
2(a)
shows the evolution of the order parameter [
ξ
=
1
κ
max
κ
max
κ
=
1
κ
C
(
κ
), where
C
(
κ
)
=
p
i
·
p
j
is the correlation
function over any two sites
i
and
j
that satisfy
κ
−
1
<
|
r
i
−
r
j
|
κ
] in a series of simulations with varying average
composition. The material is disordered at high temperature
but becomes ordered at low temperatures. The phase transition
is somewhat diffuse due to the compositional disorder. Fig-
ure
2(b)
shows the nature of the ordered phase. Remarkably,
we find that all dipoles are in the rhombohedral (
C
) states until
a critical composition of about 50% beyond which all dipoles
are in the tetragonal (
D
) states. Indeed, at a composition of
33
.
3%, a third of the sites would prefer tetragonal dipoles.
However, the exchange and electrostatic interaction with the
neighbors overwhelm this preference and instead force it
into the rhombohedral state. The opposite happens at a com-
position of 66
.
7%. The exchange of stability between the
rhombohedral and tetragonal states takes place at a well-
defined critical composition. This observation is extremely
robust: Figure
2(b)
includes results from 10 realizations.
In short, we see the
emergence of the morphotropic phase
boundary
(MPB). Remarkably, the disorder is completely sup-
pressed at all compositions except the MPB.
The resulting phase diagram is shown in Fig.
2(c)
(where
the order-disorder transition temperature is taken to be that of
the maximum curvature of the
ξ
−
β
curve). It is structurally
consistent with experimental observations with the paraelec-
tric phase at high temperature and different ordered phase
at low temperatures depending on composition. The details
depend on the parameters as we shall show later.
The resulting domain patterns are also interesting. Fig-
ure
2(d)
shows the typical domain patterns at three different
compositions (Animations of the simulations are provided in
Ref. [
24
]). We see the 2D analogs of 71
◦
or 109
◦
domain walls
with (10) normals at composition of 40%, and we see the 2D
analog of 90
◦
domain walls with (11) normals at composition
of 60%. However, at the MPB (composition 50%), we see a
mixture of rhombohedral (
C
) and tetragonal (
D
) states with a
highly fragmented and frustrated domain pattern. Figure
2(f)
shows that the average domain size falls precipitously at the
MPB compared to that at all other compositions. Figure
2(g)
shows that the oscillations in polarization in the horizon-
tal, vertical, and diagonal directions are also magnified at
the MPB. Specifically, we take the Haar wavelet transform
of the domain pattern, and Fig.
2(g)
shows the normalized
sum of squares of the horizontal, vertical, and diagonal de-
tail coefficients obtained from level one (finest level) Haar
wavelet decomposition averaged over 10 realizations. All
these observations are consistent with experimental observa-
tions. Figure
2(e)
reproduces the experimental observations of
Woodward, Knudsen, and Reaney [
17
]: classical well-defined
domain patterns are observed away from the MPB, but highly
fragmented domain patterns are observed near the MPB as
in our simulations [Fig.
2(d)
]. We note that the fragmented
patterns in our simulations at the MPB are found to be an-
gled at around 22
.
5
◦
to accommodate divergence-free rotation
of dipole states across the boundary when changing from a
tetragonal state to a rhombohedral state or vice versa. Such
domain walls typically do not exist in three dimensions lead-
ing to further fragmentation [
25
].
We comment that both the long-range interaction and the
disorder in composition are necessary for this behavior. In the
absence of the long-range dipole-dipole interaction (i.e., when
D
e
=
0), the average number of dipoles track the composition
as shown in Fig.
3(a)
. The lower critical dimension of a two-
state random field Ising model is two, i.e., the lattice remains
disordered even at low temperatures [
26
,
27
]. This is consistent
with the observations in the first row of Fig.
3(b)
where the
material remains disordered. We see some domains, with
C
-
rich states and
D
-rich states and meandering domain walls.
However, the microstructure is still disordered with all states
present, and domains with
C
-rich and
D
-rich states coexist.
This is why the average number of dipoles track the compo-
sition. Crystallographically distinguished domain walls arise
144201-3
TEH, LI, AND BHATTACHARYA
PHYSICAL REVIEW B
103
, 144201 (2021)
FIG. 3. The effect of different dipole-dipole interaction strength
D
e
. (a) and (b) show the dipole orientation and domain patterns in the
ordered phase at various compositions and
D
e
values.
J
e
=
1in(b).
with small
D
e
and become sharper as it increases [Fig.
3(b)
]
along with a complete suppression of the disorder within the
domains. Finally, at large
D
e
, we see less fragmented domains
at the MPB because the role of disorder becomes comparably
smaller. The complex domain patterns at the MPB also do not
appear when the composition is not random [
24
].
Figure
4
shows how the phase diagram changes as we
change the relative magnitudes of
h
C
D
,
h
D
C
and
J
C
e
,
J
D
e
while
keeping their averages the same. The columns show that the
order-disorder transformation temperature becomes composi-
tion dependent when
J
C
e
=
J
D
e
with the slope of the transition
temperature dependent on the relative magnitudes of
J
C
e
and
J
D
e
while the rows show that the MPB composition depends
critically on the relative magnitude of
h
C
D
and
h
D
C
. Figure
5
shows that the MPB composition also depends sensitively on
the relative polarization magnitudes
p
C
0
and
p
D
0
.
An explicit solution of this problem remains open, but we
can understand various aspects qualitatively. We first consider
the order-disorder phase transition temperature
T
c
. Recall that
according to the Onsager solution to the (two state) Ising
FIG. 4. The effect of different energetic penalties and exchange
constants for different composition types of phase diagram.
h
C
D
and
J
C
e
are varied with both (
h
C
D
+
h
D
C
)
/
2and(
J
C
e
+
J
D
e
)
/
2kept
constant at 1.
FIG. 5. The effects of dipole magnitude (
p
C
0
,
p
D
0
).
p
C
0
is varied
with (
p
C
0
+
p
D
0
)
/
2 kept constant at 1.
model on a square lattice,
T
c
=
2
J
e
/
ln(1
+
√
2)
=
2
.
27 for
J
e
=
1. At the two pure end states, there is no disorder, but
we differ from the Ising model in two important ways. First,
we have eight allowable states and we expect this to depress
T
c
. Second, we have dipole-dipole interactions; the energy
due to this is high in the disordered state but is driven to
almost zero by the formation of domains in the ordered states.
Consequently we expect the
T
c
to increase with
D
e
.Wehave
confirmed this: For
J
e
=
1, we find
T
c
to be 1
.
3
,
2
.
0
,
2
.
3
,
2
.
9
for
D
e
=
0
,
0
.
5
,
1
,
2, respectively, for the end states. In any
case, for
D
e
fixed, we expect
T
c
to be proportional to
J
e
asweseeinFig.
4
for the end states. In fact, this is also
true in the intermediate compositions since the disorder in
chemical composition has little effect on
T
c
. The average
J
e
is
(1
−
λ
)
J
C
e
+
λ
J
D
e
for a disordered lattice with volume fraction
λ
of D sites. Consequently we expect
T
c
to be proportional to
λ
asweseeinFig.
4
.
We now turn to the MPB composition
λ
c
. We consider very
low temperature where the entropy can be neglected. Since the
dipole-dipole energy suppresses the disorder, we simply need
to compare the energy between the configuration where every
lattice point is in a
C
state with one where every lattice point is
in a
D
state. Further, the formation of domains can reduce the
dipole-dipole energy. Therefore, to leading order, the energy
in these configurations is
W
C
=
λ
h
D
C
−
2
(1
−
λ
)
J
C
e
+
λ
J
D
e
p
C
0
2
,
W
D
=
(1
−
λ
)
h
C
D
−
2
(1
−
λ
)
J
C
e
+
λ
J
D
e
p
D
0
2
.
(6)
Equating the two, we obtain
λ
c
=
h
C
D
+
2
J
C
e
p
C
0
2
−
p
D
0
2
h
C
D
+
h
D
C
+
2
J
C
e
−
J
D
e
p
C
0
2
−
p
D
0
2
.
(7)
This is consistent with the results in Fig.
4
[e.g., in the middle
column where
J
C
e
=
J
D
e
=
p
C
0
=
p
D
0
=
1,
λ
c
=
h
C
D
/
(
h
C
D
+
h
D
C
)
=
h
C
D
/
2, which is 0
.
35
,
0
.
5
,
0
.
65 in the three rows] and
Fig.
5
[
h
C
D
=
h
D
C
=
J
C
e
=
J
D
e
=
1, so
λ
c
=
1
/
2
+
2(
p
C
0
−
p
D
0
) which is 0
.
42
,
1
,
0
.
58 in the three cases].
The situation is largely similar with some difference in
detail in three dimensions (see Ref. [
24
]). We observe that
the MPB emerges with local and exchange energies even in
the absence of long-range dipole-dipole interactions (
D
e
=
0).
This is consistent with the fact that local critical dimension
of a random-field Ising model is two. In a random-field Ising
model with two states and positive exchange constant, the
lattice is disordered in two dimensions (i.e., does not undergo
the order-disorder transition) while the lattice can be ordered
144201-4