A first-principles study of oxygen vacancy pinning of domain walls in
PbTiO3

Show preview

Hide preview

ar X
iv :c
on d-
m at
/0 30
52 24
v1 [
co nd
-m at.
mt rl-
sc i]
12 M
ay 20
03 A first-principles study of oxygen vacancy pinning of domain walls in PbTiO3
Lixin He and David Vanderbilt Department of Physics and Astronomy, Rutgers University, Piscataway, New Jersey 08854-8019
(May 12, 2003)
We have investigated the interaction of oxygen vacancies and 180◦ domain walls in tetragonal PbTiO3 using density-functional theory. Our calculations indicate that the vacancies do have a lower formation energy in the domain wall than in the bulk, thereby confirming the tendency of these defects to migrate to, and pin, the domain walls. The pinning energies are reported for each of the three possible orientations of the original Ti–O–Ti bonds, and attempts to model the results with simple continuum models are discussed.
PACS: 61.72.-y, 85.50.-n, 77.84.-s
I. INTRODUCTION
Ferroelectric materials are of intense interest for use in nonvolatile memory applications, in which the electric polarization of an array element is used to store a bit of information.1,2 However, the switchable polarization tends to decrease after many cycles of polarization rever- sal during device operation, a problem that is known as polarization fatigue. In the last decade, polarization fa- tigue in ferroelectrics has been under intensive study. Al- though several models have been proposed to explain this phenomenon,3 there is still no consensus, and many de- tails of the fatigue process remain unclear. Nevertheless, it is generally believed that defects, especially charged defects, play an important role. For example, a series of experiments has provided some understanding of how such defects may pin the domain walls.4,5 In particular, attention has been drawn to oxygen vacancies, which are often the most common and mobile defects in perovskite ferroelectrics. Moreover, it has been found that fatigue resistance can be greatly improved by replacing the Pt electrodes with conductive-oxide electrodes,6 which can be explained in terms of the ability of oxide electrodes to control the concentration of oxygen vacancies in the sam- ple. Furthermore, the fatigue rate in Pb(Zr,Ti)O3 films was found to be very sensitive to the oxygen partial pres- sure above the sample, suggesting that oxygen vacancies strongly affect the fatigue process.7
Mechanisms for polarization fatigue based on pinning of domain walls by oxygen vacancies have been dis- cussed phenomenologically by several authors.3,8,9 How- ever, in order to put such phenomenologically theories on a firm atomistic basis, it is important to have detailed first-principles calculations that can provide information about the structure and energetics of the ferroelectric domain walls, of the oxygen vacancies, and of the inter- actions between the two. While domain walls have been studied using Landau-
type continuum theories in earlier works,10–12 first- principles calculations are essential for an accurate mi- croscopic description of the domain walls.13–15 Most re- cently, Meyer and Vanderbilt studied 180◦ and 90◦ do- main walls in PbTiO3 using first-principles methods,
15
establishing the geometry of the domain walls at the
atomic level and calculating the creation energy of the domain walls. The domain wall was found to be ex- tremely narrow, with a width of the order of the lattice constant a; the positions of the atoms change rapidly in- side the domain wall and converge to their bulk value very quickly outside. As for oxygen vacancies, recent calculations have provided very useful information about the structure of oxygen-vacancy defects in this class of materials.16–18 However, to our knowledge, direct first- principles studies of the interactions between vacancies and ferroelectric domain walls have not yet appeared. Thus, the goal of this work is to use first-principles
calculations to examine how an oxygen vacancy would interact with a ferroelectric domain wall, and thus to shed some light on how oxygen vacancies might affect the switching process and cause polarization fatigue. In- deed, as the atomistic structure and polarization profile are very different in a domain wall than in the bulk of a ferroelectric domain, one may expect that the oxygen vacancies would also behave differently in such very dif- ferent environments. To explore these issues, we adopt PbTiO3 as a model system for this study, and calculate the formation energies for neutral oxygen vacancies in the bulk and in 180◦ domain walls of tetragonal PbTiO3. Our calculations indicate that the vacancies do have a lower formation energy in the domain wall than in the bulk, thereby confirming the tendency of these defects to migrate to, and pin, the domain walls. The order of mag- nitude of the computed pinning energy is∼ 10−1 eV, with substantial variations depending on geometrical configu- ration. The rest of the paper is organized as follows. In Sec. II,
we describe the technical details of our computational method and the supercells we used to model the domain walls. We present the results on the vacancy pinning energies of each of the three possible orientations of the original Ti–O–Ti bonds from first-principles calculations in Sec. III. A simple continuum model is introduced and used to help understand these results in Sec. IV. In Sec. V we discuss briefly the anticipated role of domain- wall pinning in fatigue, and discuss the case of charged vs. neutral defects. Finally, we summarize in Sec. VI.
1

Vb(x) Vb(y)
Vb(z) Vdw (x)
Vdw (z)
Vdw (y)
��
��
��
��
V(x)
V(x)
y
z
�
�
�
�
�
�
��
��
��
��
��
��
��
��
��
��
��
��
��
��
��
��
��
��
��
��
��
��
��
��
��
��
��
��
��
��
domain wall bulk regionbulk region z
x
(a) (b)
FIG. 1. Schematic illustration of the 8×√2×√2 supercell used in the calculations. (a) View in x-z plane. For purposes of illustration, the cell is shown divided into “bulk regions” (polarized alternately along ±zˆ, not shown) and “domain-wall regions” (comprised of the two primitive cells adjacent to the domain-wall center, which lies on a Pb–O plane), although the actual behavior is slightly more gradual. Vb and Vdw refer to possible locations of vacancies in bulk and domain-wall regions respectively. (b) View in y-z plane (with
√ 2× √ 2 supercell illustrated by the dashed lines), showing the distance between
vacancies and their periodic images (∼ a √ 2).
II. METHODOLOGY
Our calculations are based on standard density- functional theory (DFT) within the local-density ap- proximation (LDA). We use a planewave pseudopo- tential code implemented by the Vienna ab-initio Simulation Package (VASP).19,20 Vanderbilt ultrasoft pseudopotentials21 are used, with Pb 5d and Ti 3s, 3p electrons included explicitly. The 29 Ry cut-off used here was well tested and determined to be adequate for this material and these pseudopotentials.15 The structure is considered to be relaxed when the forces are less than 0.05 eV/A˚; the change of total energy at this time is typ- ically less than 1 meV. We first carried out reference calculations for domain
walls without vacancies, following closely the approach of Ref. 15. We describe the 180◦ domain wall using an 8×1×1 supercell (long direction along xˆ) with 2×4×4 k- point sampling; the tetragonal axis and polarization are along ±zˆ, and two yz-oriented domain walls divide the supercell into up and down domains of equal width. We also carried out reference calculations for vacancies
in the bulk. There are three kinds of oxygen vacancies that can be formed, depending on whether the removed oxygen atom had its Ti–O bonds along xˆ, yˆ, or zˆ, which we denote as V(x), V(y), or V(z), respectively. These bulk vacancies were studied using a 2×2×2 supercell with the tetragonal axis (and polarization) chosen along zˆ, and with 2×2×2 k-point sampling. (Thus, V(x) and V(y) are equivalent by symmetry in the bulk.) To obtain the va- cancy energy, we first calculate the total energy of the pure tetragonal 2×2×2 supercell as the reference energy. One oxygen is then removed from the supercell in such a way that the supercell remains net neutral, and the resulting structure is relaxed. We keep the volume and shape of the supercell fixed and allow only the ion po- sitions to relax. The vacancy energy is then calculated by comparing the total energy difference before and after removing the oxygen atoms from the supercell.
In order to study the interaction of the vacancies with the 180◦ domain wall, we constructed an 8×√2×√2 su- percell by doubling the 8×1×1 supercell in the y-z plane in a c2×2 or √2×√2R(45◦) arrangement, and then re- moving one oxygen atom (and its periodic images) from the interior of the domain-wall structure. This is illus- trated schematically in Fig. 1, referring here to vacan- cies labeled “Vdw” (i.e., located in the domain wall) in Fig. 1(a). The vacancies form sheets in the y-z plane, with individual vacancies separated by a distance of about a
√ 2, as indicated in Fig. 1(b). Even sparser ar-
rangements would be desirable, but the supercell already contains 80 atoms, and computational limitations make it difficult to increase the separation further. A 1×3×3 k- point mesh is used for these supercells, and once again the ion positions are allowed to relax while keeping the super- cell volume and shape fixed. (A single isolated domain wall decorated by vacancies could not expand or contract in the yˆ or zˆ directions because of the epitaxy constraint to the bulk, and our experience, consistent with Ref. 15, indicates that the lattice would expand very little in the xˆ direction if allowed to do so.) Finally, in order to calculate the energy difference for
an oxygen vacancy to be inside the domain wall, rela- tive to being in bulk, we find that it is advantageous to carry out corresponding 8×√2×√2 supercell calculations in which the oxygen vacancy has been removed from a bulk-like region of the supercell. In this way, we reduce the systematic errors associated with supercell shape, k- point sampling, etc. Our results for the binding energies of oxygen vacancies in the domain walls will normally be based on calculations of this kind, except where unavail- able as noted below.
III. RESULTS
We first report our calculations on isolated vacancies in bulk ferroelectric PbTiO3 in the 2×2×2 supercell. Using
2

the theoretical values of the lattice constants (a = 3.86 A˚, c/a = 1.0466) as obtained in Ref. 15, our calcula- tions show that the oxygen vacancies of type V(z) are more stable than V(x) and V(y) by 0.3 eV, similar to what was found by Park and Chadi for somewhat dif- ferent conditions.18 As indicated in Sec. II, however, we prefer to calculate pinning energies for oxygen vacancies in domain walls by using reference bulk vacancy calcu- lations in an 8×√2×√2 supercell that can be used with and without domain walls, in order to provide maximum cancellation of systematic errors. Such calculations will be the basis for the the results given in this section. How- ever, we will return to the use of the 2×2×2 supercell in Sec. IV for some calculations relevant to the modeling of vacancies in a perturbed environment. We then confirm the structure of the relaxed 180◦ do-
main wall. Our results for the 180◦ domain wall are very close to those of the previous calculation.15 The domain wall is located on the Pb-O plane. The atomic displace- ments converge rapidly to their bulk values, with most of the displacements ocurring within one unit cell of the do- main wall center. Consequently, the polarization reverses sharply within approximately one unit cell of the domain wall, and quickly saturates to a bulk value further away. Thus, the domain wall is extremely narrow, only about two lattice constant wide, and we find that it displays a significant x-z shear in addition to the reduced polariza- tion. With this justification, we can heuristically divide the supercell into two different regions, a bulk region and a domain-wall region, as sketched in oversimplified form in Fig. 1(a). To minimize the interactions between neighboring va-
cancies, the domain-wall supercell is doubled as described in Sec. II. The shape of the supercell in the y-z plane is indicated by the dashed lines in Fig. 1(b). The to- tal energy E0 of this domain wall structure (80-atom
8×√2×√2 supercell) is calculated for use as the refer- ence energy. To calculate the oxygen vacancy energy, we remove one
oxygen from the supercell and relax the resulting struc- ture. Recall that there are three types of oxygen vacan- cies, V(x), V(y), and V(z), according to the orientation of the Ti–O–Ti bond from which the oxygen atom has been removed. For each type, we choose one vacancy as close as possible to the center of the bulk region, and another as close as possible to the central domain-wall plane. We label these as ‘b’ (bulk) and ‘dw’ (domain wall), respectively. Thus, Vdw(x) is an x-oriented va- cancy at the domain wall, etc. In all calculations, we keep the supercell charge-neutral. All unrelaxed vacancies have an My mirror symme-
try through the defect site. In addition, the unrelaxed
Vdw(x) and Vb(x) defects have extra C (y) 2 and Mx sym-
metries respectively, because they lie precisely in, or else half-way between, the domain-wall planes. When relax- ing the structure, we observe that the two Ti ions that neighbor the vacancy site relax most, as expected. In the
TABLE I. Vacancy formation and pinning energies. ∆Edw is the energy to create a vacancy in the domain-wall region of the supercell. ∆Eb is the reference energy to create a va- cancy either in a supercell without domain walls (“Ref. with- out DW”), or in the bulk-like region of a supercell with do- main walls (“Ref. with DW”), and Ep = ∆Eb −∆Edw is the corresponding pinning energy. All energies are in eV.
Ref. without DW Ref. with DW ∆Edw ∆Eb Ep ∆Eb Ep
Unrelaxed V(x) 10.726 10.777 0.051 10.764 0.038 V(y) 10.943 10.946 0.003 10.929 −0.014 V(z) 10.939 11.102 0.163 11.098 0.159
Relaxed V(x) 10.445 10.567 0.122 10.542 0.097 V(y) 10.660 10.743 0.083 10.736 0.076 V(z) 10.459 10.720 0.261 — —
bulk region, the relaxation of Vb(x) and Vb(z) does not lead to any symmetry breaking, while for Vb(y) there are Pb ion displacements that break the My symmetry and contribute about 10 meV to the relaxation energy. In the domain-wall region, Vdw(y) also has similar distor- tions, but the energy is only lowered by about 5 meV. As for Vdw(x), which starts with a relatively high unrelaxed symmetry (four-element point group C2h), the relaxation
breaks the My and C (y) 2 symmetries so that the only re-
maining symmetry is inversion (point group C1h), and the total energy is lowered by about 20 meV. The total energy of the oxygen-vacancy supercell is
then calculated for each configuration, and a vacancy formation energy is computed using an appropriate ref- erence. The results are listed in Table I. Here ∆Edw = Ev∈dw + Eoxy − Edw, where Ev∈dw is the energy of the supercell with the vacancy in the domain wall, Eoxy is the energy of a free oxygen atom, and Edw is the en- ergy of a vacancy-free supercell with domain walls; and either ∆Eb = Ev+dw + Eoxy − Edw (“Ref. with DW”) or ∆Eb = Ev + Eoxy − E0 (“Ref. without DW”), where Ev+dw is the energy of a vacancy in the bulk-like region of a supercell containing domain walls, Ev is the energy of a vacancy in a domain-wall-free supercell, and E0 is the energy of the corresponding bulk supercell containing neither domain walls nor vacancy. As shown in Table I, all the vacancies in the domain
wall have lower formation energies than their counter- parts in the bulk. This indicates that there is an attrac- tive interaction between the vacancy and the domain wall and implies a positive pinning energy Ep = ∆Eb−∆Edw. Where available, the values for Ep obtained using the
reference supercell with domain walls is to be preferred (last two columns of Table I), since one expects a more systematic cancellation of errors in this case. However, in some situations this turns out not to be possible. For the vacancy V(z), in particular, a problem arises. We find that when we attempt to compute the relaxed en-
3

ergy Ev+dw of the supercell in which the vacancy has been placed as far from the domain walls as possible, the domain wall actually shifts its position during the relaxation in order to coincide with the vacancy, thus spontaneously converting the supercell from Ev+dw to Ev∈dw. This will be discussed further in Sec. III A. For this case, our best value for Ep of 261 meV is obtained by falling back to the use of a reference supercell without domain walls (middle columns of Table I). By looking at other cases (i.e., V(x), V(y), and unrelaxed cases), it can be seen that an uncertainty of approximately 25 meV is introduced by the use of the less preferable reference supercell. It may be noted in Table I that the formation ener-
gies ∆Eb are significantly different for Vb(x) and Vb(y) (∼200 meV), although by symmetry they should be equal in a true bulk environment. Here the differences arise from a supercell size effect connected with the arrange- ment of vacancies into sheets of fairly high density in the y–z plane. For a sufficiently large supercell we would expect these energies to become equal, because the lo- cal environments of Vb(x) and Vb(y) would be almost identical and the interactions between them would be negligible. However, in our case the vacancies are only separated by a distance of about a
√ 2, so the interactions
are not negligible. On the other hand, the vacancies in the domain walls should have similar interactions, and we can expect some cancellation of errors when arriving at the pinning energy. Thus, we have more confidence in the Ep values than in the relative formation energies of V(x), V(y) and V(z).
A. Domain-wall shift
As indicated in the previous subsection, when we at- tempt to calculate the energy Ev+dw of the Vb(z) va-
cancy in the bulk-like region of an 8×√2×√2 supercell containing domain walls, the domain wall spontaneously shifts toward the vacancy during the relaxation process. We carried out further tests using a 10×√2×√2 super- cell and observed the same phenomenon, as shown in Fig. 2. First, we put the vacancy on a TiO2 plane inside the domain wall (x = 4.5a) and allow relaxation. We observe that the Pb-centered domain wall shifts towards the vacancy and becomes a Ti-centered domain wall. We then put the vacancy on a TiO2 plane near the domain wall (x = 3.5a), and observe that the domain wall center (originally at x = 5a) again shifts to the left and becomes centered roughly at the TiO2 plane at x = 3.5a, ending up with almost same total energy as before. When we attempt to put the vacancy even farther from the domain wall at x = 2.5a, we find that a new domain wall forms at the vacancy position. Clearly, the domain wall is try- ing to shift its position in each case so as to minimize the polarization at the position of the vacancy, thereby demonstrating directly the pinning effect of oxygen va-
0 1 2 3 4 5 6 7 8 9 x [units of a]
−0.1
−0.08
−0.06
−0.04
−0.02
0
0.02
0.04
0.06
0.08
0.1
z [un
its of
c]
FIG. 2. Circles: layer-by-layer Pb-atom z coordinates for a relaxed 10×
√ 2× √ 2 supercell with domain walls centered at
x = 0 and x = 5a. Diamonds: same, but with a layer of V(z) vacancies inserted at x = 4.5a. Squares: same, but with the vacancy layer at x = 3.5a. The domain wall can be seen to shift toward the vacancies.
cancies and ferroelectric domain walls. This effect is most pronounced for the case of Vb(z)
because it has the strongest pinning energy, as can be seen in Table I. In the case of V(y), we find a similar but weaker effect. That is, if the V(y) vacancy is placed close enough to the domain wall (e.g., at x = 4.5a in Fig. 2), a similar shift can occur; but no shift occurs if the defect is placed farther away. To pursue the calculation of Ev+dw in order to obtain
∆Eb for Vb(z), it would be necessary for us to use a
supercell larger than √ 2×√2 in the y–z plane. However,
as this would be computationally prohibitive, we have instead chosen to recalculate the bulk vacancy energies using an 8×√2×√2 supercell without domain walls, as indicated in the previous subsection. This provides an alternative reference energy which, though less accurate, is available in all cases.
IV. MODEL CALCULATION
Our first-principles calculations give us a rough picture of the interactions between oxygen vacancies and domain walls. These calculations show that the domain walls can indeed be pinned by oxygen vacancies. We obtain esti- mates for the pinning energies, and find that the pinning energy for V(z) is much larger than for V(x) or V(y). We would like to understand better the physics underly- ing these results, and to appreciate which results might generalize to other situations (e.g., other ferroelectrics materials, or other domain-wall structures such as 90◦
boundaries).
4

With this motivation, we consider a model in which the vacancy formation energy depends on the immediate vacancy environment as characterized by local polariza- tions and strains. In particular, we consider a continuum description of the polarization and strain fields in the domain wall, as would occur in a Landau-type model, and assume that the vacancy energy can be expressed as a function of the local strain and polarization only. (A more sophisticated model might involve also a depen- dence on the local gradients of these fields, but we have not pursued this here.) Thus, in general we would write the vacancy formation energy as Ev(η,P), where η and P are the strain tensor and polarization vector describing the state of the local environment before removal of the oxygen atom. However, we specialize here to the case of interest, a 180◦ wall lying in a y–z plane separating tetragonal phases with polarizations along ±zˆ on either side. Thus, we focus on only the z component of po- larization, and for convenience we define a dimensionless reduced polarization p = Pz/Pbulk. As for the strain ten- sor, we have ηyy = ηyz = ηzz = 0 by lattice continuity and ηxy = 0 by symmetry. Moreover, we find that ηxx remains quite small in the domain-wall region. Thus, we focus only on the xz shear strain component and let ηs = (a/c) ηxz. The vacancy formation energy Ev(ηs, p) is thus considered as a function of the local ηs and p. Expanding E(ηs, p) in powers of p,
Ev(ηs, p) = A(ηs) +B(ηs) p 2 +O(p4) , (1)
where odd powers in p have been dropped by symmetry, and only terms up to O(p2) are retained henceforth. We then expand A(ηs) and B(ηs) in powers of ηs as
A(ηs) = a0 + a1ηs + a2η 2 s +O(η3s) (2)
and
B(ηs) = b0 + b1ηs + b2η 2 s +O(η3s) . (3)
Dropping terms beyond quadratic order in ηs, we take the coefficients a1, a2, a3, b1, b2 and b3, in Eqs. (1-3) to constitute the parameters of our model. To obtain these six parameters, we do a series of cal-
culations at a set of different ηs values for both p = 0 and p = 1. To do this, we construct 2×2×2 supercells at different fixed values of ηs and calculate the vacancy en- ergies as in the last section, allowing relaxation of ionic positions but not strains. A 2×2×2 k-mesh is used in all these calculations. We do calculations at p = 0 by enforcing inversion symmetry of the lattice. If there is no constraint of symmetry imposed, the ions relax freely, resulting in p = 1. The results of these calculations are shown in Fig. 3. It might be thought that the coefficients a1 and b1 in
Eqs. (2-3) should vanish by symmetry, but if the oxy- gen vacancy induces a distortion which lowers the lattice symmetry as discussed in Sec. III, this may not always be true. Consider, for example, the case of V(x) at ηs = 0
0 0.02 0.04 0.06 0.08 ηs
10
10.2
10.4
10.6
10.8
11
10.6
10.65
10.7
10.75
10.8
Va ca
nc y
En er
gy (e
V)
10.3
10.4
10.5
10.6
10.7
10.8 (a)
(b)
(c)
p=0
p=1
p=0
p=1
p=0
FIG. 3. Symbols: calculated vacancy formation energies vs. shear strain as obtained from 2×2×2 supercell calculations for (a) V(x), (b) V(y), and (c) V(z). Plus signs and crosses are for p = 1 and p = 0 respectively. Lines are fits to Eqs. (1-3).
and p = 0 (inversion symmetry imposed). While the un- relaxed defect has D2h symmetry, the relaxation lowers the symmetry to C2h (E, I, My and C
y 2 ) and reduces the
energy by about 65 meV relative to the case where no symmetry breaking is allowed. Actually, there are two equivalent defects, related to each other via the broken symmetry Mx, having degenerate energies at ηs = 0 but having a1 coefficients of opposite sign (i.e., opposite re- sponse to an applied xz strain ηs.) In Fig. (3) only the energy of the more stable of these two defects is plotted as a function of ηs > 0, and the dashed line is a fit to Eq. (2). The resulting (negative) value of a1 is given in Table II. Considering the other vacancies at p = 0, vacancy V(y)
has a similar symmetry breaking but its orientation is such that the degeneracy would be split by a ηyz strain,
5

TABLE II. Coefficients obtained by fitting to the model calculations. (For V(x), b1 and b2 are not needed for the pinning energy and thus are not reported.)
a0 a1 a2 b0 b1 b2
V(x) 10.705 −1.434 −33.117 0.032 V(y) 10.705 0 −8.737 0.032 0 12.315 V(z) 10.389 0 −32.008 0.462 −0.994 18.307
which is absent here. For V(z) we observe no symmetry breaking at p = 0. Thus a1 vanishes for these cases. Turning now to p = 1, we find a reversed situation:
V(x) and V(y) show no breaking of their C2v symmetry at ηs = 0, while V(z) breaks from C4v to C1h after re- laxation. Thus, b1 = 0 for V(y) but not for V(z). For non-zero ηs in Fig. (3), no further symmetry breaking is observed beyond what was already present at ηs = 0. The parameters resulting from all the fits of Eqs. (1-3)
are listed in Table II. From the fact that a1 and a2 are negative, we see that the vacancies will prefer an envi- ronment of high shear strain. Similarly, since b0 is posi- tive, they will prefer an environment of low polarization. Thus, the parameters are suggestive of a tendency for the vacancies to pin the domain walls. In order to model quantitatively the vacancy formation
energy in the domain wall, we now have to estimate the values of p and ηs that occur in a vacancy-free domain wall at the location where the vacancy would occur. Since Vdw(x) lies in the Pb-O symmetry plane, p = 0 there. From the first-principles calculations of Ref. 15, we can estimate that p ≃ 0.8 at the neighboring TiO2 plane where Vdw(y) and Vdw(z) are located. The estimation of ηs is more subtle. The problem is
that the shear strain is not well defined in the domain walls, since it depends strongly on which of the sublat- tice we follow. For example, if we define the shear strain to be ηs = δz/c, where δz, is the displacement of the ad- jacent atoms of same type in the zˆ direction, we estimate that ηs(Pb) = 0.068, ηs(Ti) = 0.054, ηs(O(x)) = −0.029, ηs(O(y)) = −0.047, and ηs(O(z)) = −0.027 in the center of the domain wall. This variation reflects the reversal of the polarization-related displacements along zˆ as one passes through the domain wall along xˆ. We could de- fine a mean shear as ηˆs = (1/5)
∑ i ηs(i) ≈ 0.005, which
as expected is about half of the “geometrical offset” de- fined in Ref. 15 (the offset occurs over approximately two unit cells). An alternative choice is the root-mean- square shear strain η¯s = [(1/5)
∑ i η
2 s(i)]
1/2 = 0.048. We expect that the most reasonable choice of an effective ηeffs should lie somewhere between these two limits. For Vdw(x), the pinning energy we get from first-principles calculation is 97 meV. Comparing with Eqs. (1-3), we find that ηeffs = 0.03 gives a reasonable agreement with the first-principles result for this case, and we thus adopt this value. The shear strain should be slightly smaller at the location of Vdw(y) or Vdw(z), half a lattice constant away from the domain-wall center, but for simplicity we
TABLE III. Environmental model parameters appearing in Eq. 4, and resulting pinning energy Ep of the model compared with the best estimate from the direct DFT calculations in Table I.
ηs p E model p (eV) E
DFT p (eV)
V(x) 0.03 0.0 0.105 0.097 V(y) 0.03 0.8 0.012 0.076 V(z) 0.03 0.8 0.187 0.261
retain the same value of ηs = 0.03 for all three defects. Using the parameters from Table II and the values of
ηs and p discussed above, we may estimate the pinning energy Ep = Ev(0, 1)− Ev(ηs, p) via Ep = (1− p2)b0 − (a1 + b1p2)ηs − (a2 + b2p2)η2s . (4)
The resulting pinning energies are reported, and com- pared with the direct first-principles calculations, in Ta- ble III. Recall that the good agreement for V(x) occurs by
construction. The agreement for V(z) is fair and that for V(y) is somewhat poor, although at least we have the correct sign of the pinning energy even in the worst case of V(y). Thus, we find that the present system is sufficiently complex that our simple model description is only partially successful. There are several reasons why this may be. As dis-
cussed in Sec. III A, the domain walls may shift if the den- sity of oxygen vacancies is high, and the shift increases the pinning energy. This may help explain the under- estimation of the pinning energies for V(y) and V(z).22
Also, the use of a 2×2×2 supercell for the calculations of Fig. 3, from which the values of Table II were obtained, means that the vacancies were much closer to the dilute limit than was the case for the vacancy-in-domain-wall calculations of Table I. This is most likely the dominant source of the discrepancy for the case of Vdw(y), which is less sensitive to ηs. Indeed, the fact that a higher density of oxygen vacancies in the domain wall leads to a larger pinning energy is consistent with a picture in which oxy- gen vacancies would tend to make planar clusters in the domain wall, thus acting to increase the pinning energy.9
Future tests on larger supercells in the y–z directions might help clarify these issues, although these still re- main intractable for the time being. Another obvious source of the discrepancies may be
the limitations of the model. It is unsatisfying that the choice of ηs is so ambiguous, and it is unclear whether two variables (ηs and p) should suffice to describe the local environment. After all, the structural distortions change rapidly as one passes through the domain wall, so that it is not clear whether a Landau-type continuum model should be expected to capture the details of the energet- ics. It might be interesting to see whether a model more like the effective-Hamiltonian description of Bellaiche et al.,23–25 which includes compositional disorder in order to treat alloys, could be successfully used here.
6

Nevertheless, we believe that our model captures some of the essential physics of the pinning mechanism of oxy- gen vacancies in the 180◦ domain walls. It gives the cor- rect sign and overall order of magnitude for the pinning energy, and correctly reflects that the pinning is strongest for Vdw(z), intermediate for Vdw(x), and weakest for Vdw(y). It also helps clarify the relative roles of strain and polarization effects in the pinning mechanism. We thus expect that it may be of some use for understanding other ferroelectric materials as well.
V. DISCUSSION
We now briefly discuss how oxygen vacancies may af- fect the ferroelectric switching process. As is well known, switching in ferroelectrics occurs not through a homoge- neous concerted reversal of the polarization in the bulk, but through the motion of domain walls separating re- gions of different polarization. Thus, insofar as these domain walls become pinned, the switching will be sup- pressed. In a pristine ferroelectric material that has a robust
hysteresis loop and a large remanant polarization, the ferroelectric domain walls are presumably only weakly pinned by some pre-existing defects. Our calculations indicate that oxygen vacancies will tend to migrate into these domain walls over time, since they experience a binding energy to the domain wall of between 100 and 250 meV. In fact, once inside the domain wall, we would ex- pect the vacancy to hop into the V(z) configuration, since this is lower in energy than the V(x) or V(y) configura- tions. Thus, the effective binding energy is ∼250 meV, the value associated with the V(z) vacancy. If a significant number of vacancies accumulate in the
domain wall, they in turn can act to pin the domain wall. When the density of such vacancies is low, they will not pin the domain walls strongly, and switching will still be able to occur. As the areal density of vacancies increases, however, an increasing fraction of domain walls (or in- creasingly large portions of individual domain walls) may become immobile, resulting in the decay of the switchable polarization. Of course, there are many limitations of our theoreti-
cal analysis, and the real experimental situation could be much more complicated. We have studied neutral oxy- gen vacancies, which correspond to vacancies of charge +2e neutralized by electrons residing in nearby states of mainly Ti 3d character. In the absence of neutraliz- ing electrons, the vacancies may pin the domain walls even more strongly. (In this case, the situation becomes more complicated, since we can expect that the domain walls may acquire a tilt in order to compensate the va- cancy charges. That is, if the 180◦ domain wall is not ex- actly parallel to the polarization, there is a bound charge ∆P · nˆ that can help neutralize the vacancy charges, an effect which may contribute to the strength of the pinning
effect.5) Alternatively, the oxygen vacancies may tend to aggregate into clusters9 or to form defect complexes of various kinds. Finally, the polarization switching process is rather a complicated dynamical process that we have not attempted to model in detail. Nevertheless, we believe that our first-principles results
serve as a first step towards understanding the possible role of oxygen vacancies in the pinning of ferroelectric domain walls. While we have investigated only one class of defects that may be involved in pinning, at least our calculations provide a lower limit for the strength of the pinning effect, which may be stronger if other defects or defect complexes play the dominant role. Our re- sults may also serve as input for more complex model- ing and simulation. For example, it could be used to extend a model such as that of Ahluwalia and Cao, who have done simulations of domain-wall formation in a 2D model simulation,26 by the inclusion of vacancies into the simulation. Our results might also be useful in the formu- lation of an effective-Hamiltonian approach27 that could be used to carry out finite-temperature simulations of the domain-wall behavior. Such studies might help quantify the effective strength of the pinning effect under more realistic conditions.
VI. SUMMARY
In summary, we have used first-principles density- functional calculations to investigate the interaction of oxygen vacancies and 180◦ domain walls in tetragonal PbTiO3. Our calculations indicate that the vacancies do have a lower formation energy in the domain wall than in the bulk, thereby confirming the tendency of these defects to migrate to, and pin, the domain walls. The pinning energies are calculated for each of the three pos- sible orientations of the original Ti–O–Ti bonds, and are found to be 97 meV, 76 meV and 261 meV for V(x), V(y) and V(z) respectively. We also introduce a sim- ple continuum model with only two parameters (p, ηs) to model the results. This simple model gives pinning energies that agree qualitatively with the first-principles calculations, and we expect that it may prove useful for other ferroelectric systems as well. This work was supported by ONR grant N0014-97-1-
0048.
1 J. F. Scott, Phys. World 8, 46 (1995). 2 O. Auciello et al., Phys. Today 51, 22 (1998). 3 A. K. Tagantsev et al., J. App. Phys. 90, 1387 (2001). 4 W. L. Warren, D. Dimos, B. A. Tuttle et al., Appl. Phys. Lett. 65, 1018 (1994)
7

5 W. L. Warren, D. Dimos, B. A. Tuttle et al., J. Appl. Phys. 77, 6695 (1995).
6 H. N. Al-Shareef, B. A. Tuttle, W. L. Warren et al., J. Appl. Phys. 79, 1013 (1996).
7 M. Brazier, S. Mansour and M. McElfresh, Appl. Phys. Lett. 74, 4032 (1999).
8 M. Dawber and J. F. Scott, Appl. Phys. Lett. 76, 1060 (2000).
9 J. F. Scott and M. Dawber, Appl. Phys. Lett. 76, 3801 (2000).
10 V. A. Zhirnov, Sov. Phys. JETP 35, 822 (1952). 11 L. N. Bulaevskii, Sov. Phys. Solid State 5, 2329 (1960). 12 W. Cao, G. R. Barsch, and J. A. Krumhansl, Phys. Rev. B 42 6396 (1990); W. Cao and L. E. Cross, Phys. Rev. B 44, 5 (1991).
13 J. Padilla, W. Zhong and D. Vanderbilt, Phys. Rev. B 53, R5969 (1996).
14 S. Po¨ykko¨ and D. J. Chadi, Appl. Phys. Lett. 75, 2830 (1999); J. Phys. and Chem. of Solids 61, 291 (2000).
15 B. Meyer and D. Vanderbilt, Phys. Rev. B 65, 104111 (2002).
16 R. I. Eglitis, N. E. Christensen, E. A. Kotomin, A. V. Post- nikov and G. Borstel, Phys. Rev. B 56 8599 (1997).
17 E. A. Kotomin, R. I. Eglitis, A. V. Postnikov G. Borstel and N. E. Christensen, Phys. Rev. B 60 1 (1999).
18 C. H. Park and D. J. Chadi, Phys. Rev. B 57, 13961 (1998). 19 G. Kresse and J. Hafner, Phys. Rev. B 47,R558 (1993). 20 G. Kresse and J. Furthmu¨ller, Phys. Rev. B 54, 11169 (1996).
21 D. Vanderbilt, Phys. Rev. B 41, 7892 (1990). 22 If we insert p=0 into Eq. (4), reflecting the idea that the domain wall may shift to become centered on the vacancy layer, we find Ep=0.040 and 0.490 eV for V(y) and V(z) respectively. This improves the agreement for V(y), but overshoots for V(z), in comparison with the direct DFT results.
23 L. Bellaiche, A. Garc´ıa and D. Vanderbilt, Phys. Rev. Lett. 84, 5427 (2000).
24 R. Hemphill, L. Bellaiche, A. Garc´ıa and D. Vanderbilt, Appl. Phys. Lett. 77, 3642 (2000).
25 A. M. George, J. I´n˜iguez and L. Bellaiche, Nature 413, 55 (2001).
26 R. Ahluwalia and W. Cao, J. Appl. Phys. 89, 8105 (2001). 27 J. Padilla, W. Zhong, and D. Vanderbilt, Phys. Rev. B 53, R5969 (1996).
8