PHYSICAL REVIEW B 72, 104104 共2005兲

Slater and Takagi defects in KH2PO4 from first principles

J. Lasave,1 S. Koval,1,* N. S. Dalal,2 and R. Migoni1

1Institutode Física Rosario, Universidad Nacional de Rosario, 27 de Febrero 210 Bis, 2000 Rosario, Argentina
2Department of Chemistry and NHMFL, Florida State University, Tallahassee, Florida 32306, USA
共Received 5 April 2005; revised manuscript received 1 June 2005; published 9 September 2005兲

We performed ab initio electronic structure calculations to study the stability of the Slater and Takagi
protonic defects in the ferroelectric phase of KH2PO4. We show that the Slater defects are stabilized only when
they are aligned forming chains of defects in the lattice, in accordance with the displacive model of the
transition. In contrast, the Takagi defects are not stable in the lattice. The ab initio energy estimation for the
Takagi defect 共⬇54 meV兲 is in fair good agreement with previous phenomenological estimates. However, the
calculated uncorrelated Slater defect energy 共⬇17 meV兲 is larger than that from previous model assignments
共5.2 meV兲, while the latter agrees well with our result for the correlated value per formula unit. This suggests
that the earlier estimates contain considerable contribution from correlations between configurations of differ-
ent H2PO4 groups. Our ab initio results show that the formation of chains of Slater defects is accompanied by
phosphate rotations, and a concomitant lattice contraction in the basal plane. Therefore, we propose that the as
yet unexplained anomalous lattice contraction observed in the neighborhood of the phase transition is due to an
increase of the Slater defect population.

DOI: 10.1103/PhysRevB.72.104104 PACS number共s兲: 77.84.Fa, 71.15.Nc, 77.80.⫺e

I. INTRODUCTION bond common to two polar states moves to the other bond
Potassium dihydrogen phosphate 共KH2PO4, or KDP兲, is side. This leads to the formation of a Takagi-pair defect with
the prototype of the hydrogen-bonded ferroelectric com- energy 2⑀T in two neighboring phosphates that finally remain
pounds. It has been the subject of intensive experimental and with one and three protons 关see Fig. 1共c兲兴.
theoretical studies since the discovery of its ferroelectric Later improvements of these order-disorder models came
transition at 123 K in 1935.1 Despite considerable progress with the addition of long-range interactions5,6 and the intro-
since then, a complete understanding of the ferroelectric duction of tunneling7 to account for the huge H ↔ D isotope
transition mechanism is still lacking. effect observed: Upon H → D substitution, the Tc increases
From the theoretical side, the first successful model of the from 123 to 224 K. Models that include proton-lattice inter-
KDP phase transition was developed by Slater in 1941.2 In actions, and hence may display displacive features of the
the Slater model, the ice rules are applied, and only two transition, were first introduced by Kobayashi.8 Actually,
protons are allowed to be close to each phosphate in six strong evidence that the transition involves a displacive
possible configurations 共Fig. 1兲. The lowest-energy configu- component with electronic instabilities was recently detected
rations are the polar ones 关see Fig. 1共a兲兴, with electric dipoles by nuclear magnetic resonance 共NMR兲 experiments for
pointing along the z direction, which characterize the low- KD1.8PO4 共Ref. 9兲 and squaric acid.10 Particularly, in the
temperature ferroelectric phase, with a net spontaneous po-
larization along this direction. The two polar configurations
are built with protons attached to the bottom or the top oxy-
gens in the phosphate, and differ in the sign of the corre-
sponding dipoles along z. The remaining four configurations
are the ones 共hereafter called “Slater”兲 where the two protons
are laterally attached, respectively, to two oxygens, one at the
top and the other at the bottom of the phosphate units 关Fig.
1共b兲兴. These configurations have a higher energy ⑀S 共Slater
energy兲 compared to that of the polar configurations in the
ferroelectric phosphates. The Slater model predicts a sharp
first-order ferroelectric transition, but because it is a static
model in its original form, it is difficult to use it for under-
standing, in particular, dynamic properties, such as electric
transport properties, and related protonic hopping in the low-
temperature ferroelectric phase.3
Takagi improved the theory by including the possibility of
configurations with one or three protons attached to the phos- FIG. 1. 共Color online兲 Different proton configurations for the
phate 共Takagi configurations兲 with energy ⑀T per phosphate phosphates in KDP: 共a兲 Polar, 共b兲 Slater, and 共c兲 Takagi-pair defect.
above that of the polar configurations.4 These configurations Left: Top view. Right: Frontal view. ⑀S and 2⑀T are the energies of
violate the ice rules and arise, e.g., when a proton from a H the formation of a Slater and a Takagi-pair defect, respectively.

1098-0121/2005/72共10兲/104104共8兲/$23.00 104104-1 ©2005 The American Physical Society

LASAVE et al. PHYSICAL REVIEW B 72, 104104 共2005兲

latter case, the coexistence of order-disorder and displacive results are in agreement with neutron diffraction,11,12 Comp-
features appear to be present near the phase transition. ton scattering,32 and electron-nuclear double resonance
Recent high-resolution neutron diffraction experiments 共ENDOR兲 measurements performed recently.33
provided evidence that the structural modifications in the hy- However, despite the evident importance of proton con-
drogen bond geometry are strongly correlated with the order figurational defects in the dynamical behavior of these mate-
parameter.11,12 In this regard, these experiments show that, in rials, their energies and stabilities have been estimated so far
the paraelectric phase, the distance ␦ between the two stable only at a rather phenomenological level. To our best knowl-
minima of the proton in the H bond is strongly correlated edge, no ab initio calculations have been carried out to study
with Tc. It is also argued that the origin of the isotope effect the Slater or Takagi defects in H-bonded systems. In this
is the geometric effect and not the tunneling itself, and that work, we investigate the stability of different proton configu-
the proton and the host are definitively not separable
rations by means of ab initio calculations within the density
systems.13
functional theory 共DFT兲. Details of the ab initio method and
Phenomenological models were also proposed recently to
take into account the geometrical effects with some the calculations are presented in Sec. II. In Sec. III, we ana-
success.10,14–17 However, such model parametrizations are lyze the results obtained and finally, we discuss and elaborate
not always reliable regarding the physical meaning of the our conclusions in Sec. IV.
parameter values and their transferability upon structural or
chemical changes. II. AB INITIO METHOD AND COMPUTATIONAL
Earlier magnetic resonance studies of DKDP 共KD2PO4兲 DETAILS
showed the existence of deuteron jumping between H bonds,
that is consistent with the Slater–Takagi theory of defects.18 The first-principles calculation were carried out using the
The domain-wall motion near the transition in DKDP was SIESTA program,34 which is a fully self-consistent approach
explained in terms of the effective motion of Takagi groups, that employs a linear combination of confined pseudoatomic
leaving behind a string of Slater defects.19 The slow motion orbitals 共PAO兲 of the Sankey–Niklevsky type as basis func-
of domain walls is also believed to be the cause of an tions for the valence electrons.35 We have chosen a double-
anomaly in the dielectric constant of KDP as a function of zeta basis set with polarization functions 共DZP兲 for these
temperature, which shows a plateau region immediately be- orbitals. The interaction between ionic cores and valence
low Tc. The sudden fall of the plateaulike anomaly at TF electrons is represented by angular-dependent norm-
⬇ Tc − 60 K is described as a “kink” and is related to the conserving pseudopotentials of the Troullier–Martins type.36
freezing of the domain-wall motion.20 The stability and mo- Due to the large overlap between core and valence charge
tion of thin domain walls was recently described by a Slater– densities, we have also included nonlinear core corrections
Takagi model which also includes interactions between pro- 共NLCC兲 for a proper description of the K ion. Exchange-
tons in parallel bonds.21 correlation energy functionals are computed within the
The Slater and Takagi energies 共Fig. 1兲 were deduced ear- gradient-corrected Perdew–Burke–Ernserhof 共PBE兲 approxi-
lier using phenomenological models 共hydrogen-bond occu- mation.37
pation models兲 that were fitted to the experimentally deter- The grid computation of the numerical integrals was car-
mined polarization curves in DKDP.6 In the model of Ref. 6, ried out with an energy cutoff of 125 Ry. The orbital con-
the original Slater–Takagi model that takes into account the finement energy Ec, which measures the degree of approxi-
short-range interactions between protons in the same phos- mation of the confined PAO,34 is set to 50 meV, which gives
phate, was extended to include also long-range interactions total energies and geometries of sufficient accuracy. These
with protons in neighboring phosphates. A good fit to the approximations were also tested against results from stan-
polarization curves was achieved for DKDP, although the dard pseudopotential plane-wave calculations.24,25,27 The
calculated entropy was nearly twice that found in specific- electronic Brillouin-zone sampling is reduced to the ⌫-point,
heat experiments.6 This discrepancy was corrected by Blinc which proved to be a good approximation due to the large
and Svetina with the incorporation of the tunneling effects to supercells used in the calculation.25
this extended model.22 They determined the Slater–Takagi The ferroelectric 共FE兲 phase structure in KDP is described
energies of the different proton configurations around a given with the conventional face-centered-tetragonal 共fct兲 cell, con-
phosphate by integrating out all internal coordinates in a taining 8 formula units 共64 atoms兲, and axes rotated through
given PO4 group. Using a different fitting procedure, but 45° with respect to the conventional body-centered-
including also quantum effects from a small-tunneling ap- tetragonal 共bct兲 cell. A schematic view of the FE structure
proximation, Fairall and Reese tabulated the parameters of from the top 共projection on the xy plane兲 is shown in Fig.
the hydrogen-bond occupation 共order-disorder兲 model for 2共a兲. In this phase, two hydrogens are covalently attached to
various members of the KDP-type family.3,23 each phosphate in the polar configuration 关see Fig. 1共a兲兴,
Ab initio calculations were conducted only very recently which has the lowest possible energy.
in KDP.24–31 It was found that significant instabilities arise The main calculations were carried out with the cell pa-
only when protons are allowed to relax together with heavy rameters fixed at the experimental values at T = TKDP
c − 20 K
ions.25–27 The huge deuteration isotope effect is shown to for KDP, which is orthorhombic.38 After relaxing all of the
arise due to a self-consistent coupling, via modifications of atoms in this cell, we obtained the FE 共polar兲 phase as de-
the H-bridge covalency, between tunneling and the geometri- picted in Fig. 2共a兲, with internal parameters in general good
cal effects, thus dramatically enhancing the latter.25 These agreement with experiments.24,25,31

104104-2
SLATER AND TAKAGI DEFECTS IN KH2PO4 FROM… PHYSICAL REVIEW B 72, 104104 共2005兲

FIG. 3. 共a兲 Total energy as a function of the hydrogen displace-

ment ␦1 in the reference bond 共solid curve and circles兲. Dashed
lines and small left-triangle symbols correspond to the higher-
energy continuation of the Slater and polar energy curves. 共b兲 Pro-
ton arrangement factor K␦ as a function of ␦1. Different symbols
and lines correspond to the eight different phosphates in the con-
FIG. 2. 共Color online兲 Conventional fct cell of KDP from a top ventional cell. 共c兲 Variation of the angle parameter ␪ with respect to
view showing: 共a兲 The FE phase, where each phosphate is in the the polar state, averaged over the group of phosphates G1
polar configuration 共P structure兲. The arrow shows the direction of = 共8 , 5 , 6 , 4兲 共triangles兲 and G2 = 共7 , 2 , 3 , 1兲 共squares兲 关see Fig. 2共b兲
the reference H1 shift to produce a Takagi pair. The labeling of the for the labels assigned to the phosphates兴. ⌬␪ is in degrees. Lines
four H bonds surrounding one of the phosphates is indicated. ␦1 is are only a guide for the eyes.
the distance between positions, equidistant from the middle of the
reference O–H–O bond labeled by 1. 共b兲 Formation of a chain of reference bond where the hydrogen is moved by constraining
four Slater defects connected through periodic boundary conditions it at different positions along the bond in order to create the
共S structure兲. Dashed arrows indicate the direction of hydrogen mo- defect. In the case of formation of polar, Takagi, and Slater
tion to form the chain. Solid arrows show the dipole moments de- configurations 共see Fig. 1兲, K␦i will approximately assume the
veloped by each phosphate in Slater configuration, which point ap- values of 2, 1, and 0, respectively.
proximately in the xy plane. The labels assigned to each of the
phosphates in the cell are in brackets.
III. RESULTS
In order to study the formation of configurational defects A. Stability of defects
involving local distortions of the polar state, we define the
magnitude K␦i for each phosphate as follows: We first explore the possibility of formation of a Takagi-
pair defect 关see Fig. 1共c兲兴. To this purpose, the hydrogen H1
K␦i = 冏 2␦1
冏
␦ i1 + ␦ i3 − ␦ i2 − ␦ i4
, 共1兲
of a reference H bond is moved in the direction of the arrow
关see Fig. 2共a兲兴 toward the opposite side of the bond, and its
position is constrained at different values of ␦1 in the polar
where ␦ ij is the ␦ value corresponding to the jth H bond 共j phase 共P structure兲. All of the remaining atoms of the cell are
= 1–4兲 connecting the ith phosphate of the fct cell 共i = 1–8兲 allowed to relax in such a way that the computed total energy
with its neighbors 关see Fig. 2共a兲兴. ␦ ij is considered to be is minimized.
negative if the corresponding H lies on the ith phosphate side In Fig. 3共a兲, we show the energy as a function of the
of the jth bond, and positive otherwise. ␦1 is the ␦ value at a hydrogen H1 displacement 共solid curve兲. As the hydrogen

104104-3
LASAVE et al. PHYSICAL REVIEW B 72, 104104 共2005兲

leaves the equilibrium position for the P structure 共␦1 smaller energy, by about 2 meV, than that shown in Fig. 3共a兲.
⬇ −0.24 Å ⬅ ␦ P兲, the energy increases. However, after sur- Hence, our conclusions do not change significantly if we
passing the middle of the H bond 共␦1 ⬇ 0.07 Å ⬅ ␦c兲, the en- relax the constraint imposed, which is however necessary to
ergy starts to decrease. In fact, on the other side of the equi- obtain the whole curves of Fig. 3, and not just the results at
librium position along the bond, there appears another the minima.
minimum of the total energy at ␦1 = + 0.23 Å ⬅ ␦S ⬇ −␦ P. For In addition to the minimum corresponding to Slater de-
values of ␦ ⬎ ␦S, as the H moves further away from the latter fects for ␦1 ⬎ ␦c 共Region II兲, we found around ␦S other lo-
minimum and approaches the opposite oxygen, the energy cally stable configurations with higher energy. These con-
increases again. figurations resemble those found in region I corresponding to
We show in Fig. 3共b兲 the parameter K␦i as a function of ␦1 the Takagi-pair formation, and lie on the continuation of the
for each of the eight phosphates in the conventional fct cell minimum-energy curve as is shown by small triangle sym-
used. As the hydrogen H1 is shifted toward the middle of the bols and dotted lines in Region II of Fig. 3共a兲. A small acti-
bond, we observe the tendency of the system to form a vation barrier connects these higher-energy states with the
Takagi pair. Particularly, at ␦1 = ␦c, K␦ assumes a value of ⬇1 minimum-energy S structures shown by the solid curve in
for the pair of phosphates connected by the hydrogen H1. Region II of Fig. 3共a兲. Conversely, for ␦1 ⬍ ␦c 共Region I兲, we
Nearest-neighbor hydrogens are also strongly perturbed at found higher-energy configurations resembling those with
this stage, and their positions are shifted with respect to those Slater defects of Region II, which are also shown as dotted
in the polar phase as is seen in Fig. 3共b兲, which is a conse- lines and small triangle symbols 关Fig. 3共a兲兴 and are a con-
quence of the strong correlations among ionic displace- tinuation of the minimum-energy curve of Region II into
ments.25 We will see below that this has its origin in phos- Region I.
phate rotations coupled to the H motion. Summarizing, we clearly observe two different regions in
With increasing H1 displacement, we observe at ␦1 = ␦c a Fig. 3共a兲: One for ␦1 ⬍ ␦c where the Takagi-pair formation
sudden change of configuration for one-half of the phos- takes place, and the other for ␦1 ⬎ ␦c characterized by the
phates in the cell, indicated by their values of K␦ ⬇ 0 for ␦1 formation of Slater defects chains 共S structures兲.
⬎ ␦c 关see Fig. 3共b兲兴. Moreover, the energy minimum seen in
Fig. 3共a兲 at ␦1 = ␦S, corresponds to these four phosphates in B. Energy of Slater and Takagi defects
Slater configuration 共K␦ ⬇ 0兲, while the other phosphates re-
main in polar configurations 共K␦ ⬇ 2兲. The four phosphates in We now discuss our ab initio study of the Slater–Takagi
the Slater configuration are linked and build a chain, which is energy parameters, ⑀S and ⑀T. As we already mentioned, a
extended to infinity through periodic boundary conditions. Takagi-pair defect alone is not stable in the lattice 共see Fig.
We call this configuration, depicted in Fig. 2共b兲, the S struc- 3兲. However, we may estimate its energy, 2⑀T, by constrain-
ture. In the polar P phase, all of the dipole moments point ing the reference hydrogen H1 to lie on the opposite side
along the ferroelectric z axis 关see Fig. 2共a兲兴. The formation of 共␦1 = −␦ P兲 in the FE phase 关Fig. 2共a兲兴. The oxygens and the
Slater defects in the S structure leads to a change of the remaining hydrogen atoms are fixed in their positions in this
dipole-moment directions, which turn perpendicular to the z phase, and only the P and K atoms are allowed to relax. In
axis in the affected phosphates 关see solid arrows in Fig. this way, we are not including the effect of correlations be-
2共b兲兴, while Takagi-pair defects do not appear in this struc- tween H1 and the other hydrogens in the lattice, which are
ture. In fact, the positioning of the H1 at the opposite side of part of the long-range interactions.22
the reference bond produces a concomitant displacement of For the calculation of the Slater energy, ⑀S, we moved H1
other three hydrogens in the directions of the dotted arrows, and H2 关see Fig. 2共a兲兴, fixed them at the opposite side of
as shown in Fig. 2共b兲. These H motions are strongly corre- the corresponding H bonds in the polar phase, and relaxed
lated and the energy is minimized by the formation of a the P and K atomic positions. Thus, we created two Takagi
chain of Slater defects. defects and a Slater defect at the middle of three neighboring
It is possible that S structures with finite length are stable, phosphates. We computed its energy ES2T, which together
the correlation length being larger than the size of the fct cell with the previously obtained value of ⑀T allows us to
used with periodic boundary conditions. In such a case, the S determine ⑀S:
structure should end with Takagi defects. Unfortunately, our ⑀S = ES2T − 2⑀T . 共2兲
computer capability does not allow us to consider larger
cells. For further verification of the correlations among Slater We have also determined for comparison the correlated
defects, we have also studied the stability of the defects Slater energy, which includes part of the long-range correla-
chain while constraining any one of these hydrogens to re- S = Etot共␦S兲 / 4. Here, Etot共␦S兲 is the total ab initio
tions, as ⑀corr
main fixed at its original position in the polar phase. In this energy at ␦S, where the chain of Slater defects becomes
case, the correlated displacement is aborted due to the fact stable 关see the minimum in Region II of Fig. 3共a兲兴. Interest-
that it would lead to the formation of a Takagi pair within the ingly, we have that ␦S ⬇ −␦ P. Similarly, a correlated Takagi
chain, which demands much more energy than the formation energy ⑀Tcorr is defined as the local-minimum energy at ␦S,
of the infinite chain. We want to stress here, that if we further represented by the dashed line in Region II of Fig. 3共a兲. In
relax the reference hydrogen from its position at ␦S, together this case, we mean “correlated” in the sense that all atoms of
with all atoms in the S structure, we again obtain approxi- the cell, except the H1, are allowed to relax to their equilib-
mately the same chain of Slater defects with a slightly rium positions.

104104-4
SLATER AND TAKAGI DEFECTS IN KH2PO4 FROM… PHYSICAL REVIEW B 72, 104104 共2005兲

TABLE I. Ab initio energy per phosphate in meV obtained for tively. The phosphates of a given group are on alternating xz
the different proton configurations studied in this work. Also shown planes along the y axis, i.e., G1 = 8,5,6,4 and G2 = 7,2,3,1 关see
for comparison are previous model results fitted to polarization numbers in brackets identifying the phosphates in Fig. 2共b兲兴.
curves in KDP.a For the sake of simplicity, let us consider the average
rotation-angle variation for each phosphate group G 共G1 or
Energy per phosphate 共meV兲 G2兲:

冉兺 冊
Parameters Ab initio Modela 1
⌬ ␪ G共 ␦ 1兲 = ⌬ ␪ i共 ␦ 1兲 . 共4兲
⑀T 53.5 60.9 4 i苸G
⑀Tcorr 25.0 In Fig. 3共c兲, we plot ⌬␪ ⬅ ⌬␪ G共␦1兲 as a function of ␦1 for
⑀S 16.9 5.2 each group of phosphates, G1 共triangles兲 and G2 共squares兲.
⑀corr
S 5.0 As ␦1 increases and H1 approaches the middle of the j
a
See Ref. 23. = 1 H bond, ⌬␪ becomes increasingly negative for G2 while
it remains almost constant for G1. This holds as long as ␦1
⬍ ␦c, as observed in Fig. 3共c兲. Thus, while a Takagi-pair
The calculated parameters are shown in Table I, together centered at H1 is formed, the phosphates of G2 linked to the
with results from a model derived to fit the polarization H-donor phosphate 关i = 4, in Fig. 2共b兲兴 rotate counterclock-
curve in KDP by Reese et al.23 We find a Takagi energy wise. This produces a shrinking of the O–O distances in the
similar to that derived by model calculations. However, we H bonds connected to the donor phosphate, which accounts
obtain a Slater energy substantially larger than that deter- for the strong correlations observed during the Takagi-pair
mined by the phenomenological calculations in Ref. 23. On formation 共Region I in Fig. 3兲. Moreover, as ␦ → ␦c, these
the other hand, the correlated parameters ⑀Tcorr and ⑀corr
S are, rotations and the related O–O bond shrinking would favor
respectively, ⬇1 / 2 and ⬇1 / 3 those of the corresponding un- the correlated proton hopping necessary to produce the chain
correlated parameters, showing the existence of strong cor- of Slater defects observed in Region II. In fact, soon after
relations in the system.25 We notice that the correlated pa- ␦ ⬎ ␦c, a drastic change in the behavior of ⌬␪ for both groups
rameter ⑀corr
S approaches the value obtained by the Reese is observed. Both groups rotate clockwise more than one
calculations 共see Table I兲. This fact leads us to speculate that degree at ␦c, as shown in Fig. 3共c兲. These last rotations are
the Slater energy parameter involved in the earlier concomitant with the correlated jumping of the four H’s 关see
models5,6,23 effectively contains a significant contribution of Fig. 2共b兲兴, which leads to the S structure at ␦S, as observed in
correlations. Figs. 3共b兲 and 3共c兲. Although both groups rotate in the same
sense, the O–O contractions observed at ␦S are explained by
C. Correlations among Slater defects and phosphate rotations the fact that ⌬␪ G1 ⬎ ⌬␪ G2.
In the process of forming a chain of Slater defects, we
observed other important structural changes. Besides the
phosphate distortions which give rise to the dipole moments D. Thermal expansion anomaly
in the xy plane 关see Fig. 2共b兲兴, we have also found that the The contractions in all of the O–O distances found in the
interphosphate O–O distances at ␦S shrink about 0.3% with S structure, would in turn affect the crystal cohesion. We
respect to those of the polar phase at ␦ P. If 具dOO典i is the show in this section that this feature can be related to an
average of the O–O distance between the oxygens in the ith anomalous lattice contraction observed as T increases ap-
phosphate and the oxygens linked to them in the neighboring proaching Tc in KDP and other related compounds.39,40 To
phosphates, we have verified that 具dOO典i共␦S兲 ⬍ 具dOO典i共␦ P兲 for investigate this feature in detail, we relaxed the cell param-
each phosphate in the cell. Therefore, since all interphos- eters and all atomic positions in both configurations studied
phate O–O distances decrease with respect to those of the in this work: 共i兲 The P structure 共FE phase兲, and 共ii兲 the S
polar phase approximately by the same amount at ␦S, we structure 关Figs. 2共a兲 and 2共b兲, respectively兴. It is worth men-
conclude that the appearance of an S structure is associated tioning that the full relaxations do not alter the described
with a phosphate rotation around the z direction. To investi- characteristics of the P and S structures. Both relaxed cells
gate this feature in more detail, we determined, for each are orthorhombic with cell parameters 共a P , b P , c P兲 and
phosphate, the smallest angle ␪ in clockwise sense from the x 共aS , bS , cS兲 for the P and S structures, respectively. To better
axis up to the projection onto the xy plane of a P–O bond compare the effect of the different structures in both cells, we
共see Fig. 2兲. This angle defines the orientation of the phos- define the pseudomonoclinic cell with parameters 共a⬘ , ␥ , c兲,
phates around the z axis as in Ref. 38. Then, for each phos- where a⬘ = 21 冑a2 + b2 and ␥ = 2 arcsin共b / 冑a2 + b2兲. Results
phate i, we computed as a function of ␦1 the difference ⌬␪ i from the relaxations are shown in Table II, together with
in the rotation angle respect to the polar phase at ␦ P: experimental data from Ref. 39 at different T as they tend to
Tc ⯝ 122 K from the FE phase.
⌬␪ i共␦1兲 = ␪ i共␦1兲 − ␪ i共␦ P兲. 共3兲
The relaxed FE cell is obtained in good agreement with
With ␦1 increasing from ␦ P, we found that there are basi- experiment. The calculated S cell is contracted in the xy
cally two phosphate groups G1 and G2, each having similar plane and elongated in the z direction approximately 0.1%
variations of phosphate rotation angles ⌬␪1 and ⌬␪2, respec- relative to the relaxed cell of the FE phase. We also note that

104104-5
LASAVE et al. PHYSICAL REVIEW B 72, 104104 共2005兲

TABLE II. Ab initio parameters of the pseudomonoclinic cell clusters of Slater defects in the lattice. At T = 0 K, the system
obtained by relaxation of the structures P and S. Also shown are the is in its FE ground state, but as T increases, chains of Slater
experimental parameters obtained by Kobayashi et al. at different defects within the S structure, which are locally stable,
temperatures approaching Tc.a Parameters a⬘ and c are given in Å, would have more statistical weight to effectively appear in
and ␥ in degrees. the dynamics of the system. In fact, the barrier between the P
and S configurations is of the order of ⑀corr S ⬅ Etot共␦S兲 / 4 per
Ab initio structures Experimental dataa phosphate 关see Fig. 3共a兲兴, and hence, at temperatures around
Cell 122.76 K TS ⬅ ⑀corr
S / KB ⯝ 60 K, finite-size clusters with the S structure
parameters P S 103.86 K 113.56 K 共Tc兲 would have a nonvanishing probability to coexist with P
structure clusters in the lattice. Of course, quantum delocal-
a⬘ 7.580 7.575 7.4290共7兲 7.4290共7兲 7.427共7兲 ization effects and thermal fluctuations would lead to an even
␥ 90.28 90.07 90.433 90.414 90.262 smaller effective barrier between such configurations.25,31
c 6.945 6.952 6.9264共7兲 6.9272共7兲 6.926 Then, the coexistence of Slater and polar configurations as
T → Tc 共note that TS ⯝ 1 / 2Tc兲 would in turn produce an in-
a
See Ref. 39. crease in the probability of lattice contractions in the xy
the parameter ␥S is about 0.3% smaller than ␥ P. Moreover, plane, due to the fact that the S structure exhibits smaller a⬘
␥S is almost 90°, i.e., the S cell may be considered as nearly and ␥ parameters than the P structure. This fact qualitatively
tetragonal, in contrast to the FE cell which is clearly ortho- explains the lattice contractions in the basal plane and the
rhombic 共see Table II兲. The cell contraction in the xy plane is tendency to a smaller orthorhombicity, observed experimen-
a consequence of the phosphate rotations and the concomi- tally as T → Tc 共see Table II兲.
tant O–O bond shrinking discussed in the previous subsec- The appearance of Slater defects would also elongate the
tion. parameter c, which is in accordance with the experimental
As was pointed out in Ref. 40, an anomalous contraction trend for this magnitude between 100 K and 115 K 共Table
of the cell parameters a and b is observed experimentally II兲. On the other hand, our analysis cannot account for the
when T increases tending to Tc in KDP. This anomaly was final contraction observed in the parameter c very near to Tc,
also observed in squaric acid and PbHPO4. The experimental which might be related with an ultimate PO4 reorientation
contraction in the xy plane is shown in Table II, in terms of produced dynamically near the phase transition.40,43
the parameters a⬘ and ␥ of the pseudomonoclinic cell. The
contraction in these parameters is about 0.05–0.2% as T in- IV. DISCUSSIONS AND CONCLUSIONS
creases from Tc − 20 K to Tc. It is also observed experimen-
tally as T → Tc that the orthorhombicity of the cell tends to The existence of strings of Slater configurations created
decrease, i.e., ␥ diminishes and approaches the value of 90° by the migration of Takagi defects was postulated several
共see Table II兲. decades ago and shown to be consistent with deuteron NMR
The behavior of the unit cell parameter c in the region experiments.19 These strings are visualized in phenomeno-
兵Tc − 20 K , Tc其 is somewhat different than that for the basal logical simulations, as part of domain walls responsible of
parameters. Actually, as shown in Table II, in the first part of the dielectric anomaly observed at the domain freezing tem-
this region, the parameter c grows and it finally experiments perature TF.21 Our ab initio calculations enable the formation
a contraction very near to Tc. As a result of this, the total cell of such nonpolar layers in KDP, giving further support to the
volume experiences a final contraction as Tc is approached results obtained in Ref. 21. According to our results, these
from the FE phase. layers would arise macroscopically in the lattice at tempera-
The existence of Slater defects near the transition in the tures around TS ⬇ 60 K, in remarkable coincidence with TF.
paraelectric phase was inferred from the 75As nuclear quad- The problem of how the O-H ¯ O proton exchange takes
rupole resonance 共NQR兲 experiments in KH2AsO4 and a place is a long pending issue.41 Some experiments, e.g., Ra-
family of related H-bonded arsenates.41 Molecular dynamics man studies,44 seem to indicate the importance of the H2PO4
simulations with a modified strong dipole-proton coupling order-disorder model, in which only the two up-down H2PO4
model 共MS-DPC兲 including Slater contributions16 explained configurations of C2 symmetry are involved. However, the
75
satisfactorily the NMR data of the 31P nucleus in As NQR data seem to indicate that, besides the polar states,
KD1.8H0.2PO4 共Ref. 9兲 and the 75As NQR data.42 These cal- the Slater configurations 共C1 symmetry兲 also take part in the
culations indicated also that the S-shaped anomaly in the T proton exchange.42 Although the NQR experiments were not
dependence of the local order parameter in KDP is related to carried out on KDP itself, but on other compounds
the mentioned NQR data and that it originates from the 共KH2AsO4 and related materials兲, it is claimed that the C1
change of population between polar and Slater configurations configurations play also an important role in the KDP phase
on approaching Tc. For instance, the sudden decrease of the transition as well. Moreover, it is argued that the smaller
local-order parameter found as T increases and approaches statistical weight of the C1 configurations and signal-to-noise
Tc, was correlated with the appearance of locally induced problems in the scattering experiments are the origin of the
Slater configurations.16 absence of indications of these configurations in Raman
The last discussion enables us to propose that the ob- measurements in KDP.41
served lattice contraction, in the xy plane as T grows and Regarding this controversy we can say on the one hand,
reaches Tc, can be explained in terms of the appearance of that strong correlations allow the existence of Slater defects

104104-6
SLATER AND TAKAGI DEFECTS IN KH2PO4 FROM… PHYSICAL REVIEW B 72, 104104 共2005兲

only as long line-defects. This is in favor of a displacive-like mode and the proton soft mode was observed recently in
picture of the transition, as recently indicated by NMR another ab initio calculation.46
measurements.9 On the other hand, if we move a line of four Our calculations show that as H approaches the middle of
hydrogens in one step 关dashed arrows in Fig. 2共b兲兴, and in the O–O bond, this bond contracts and the phosphates rotate.
each successive step we displace a neighboring H line in the In a dynamical sense, it implies that this concerted motion
xz plane in alternating directions along x, then after complet- would increase the hopping probability of the hydrogens be-
ing four steps in our cell, we change the polarization sense in tween the two sides of the bond. We have shown in Fig. 3共c兲
the polar phase. We found that the largest barrier for each of that the hopping actually takes place at ␦c, and the phos-
the four steps is of the order of 4⑀corr
S , which is much smaller
phates rotate again to the other direction. We want to stress
than the barrier for the global polarization flip along z. This here that this phenomenon, although shown in a static simu-
picture is in agreement with a model according to which the lation, should be related to a mechanism of assisted proton
Slater configurations appear in the proton exchange,41 as hopping. That is, protons may jump collectively to the other
verified by the NQR data.42 We conclude that in a scenario of side of the bond, helped by phosphate rotations that lower
larger simulation cells, with more than eight KH2PO4 units, the effective barrier for the H transfer.
the polarization flip would then take place in steplike form, In summary, we have shown by means of ab initio elec-
with the formation of chains of Slater defects having a pre- tronic structure calculations, that Slater defects are stabilized
ponderant role. Moreover, this process would also be in- by their nucleation in chains in the FE phase of KDP, a fact
volved in the motion of domain walls.21 that gives further support to the displacive picture of the
Raman spectra in the x共yx兲y configuration show a well- transition. The uncorrelated Takagi energy parameter
defined inelastic side peak 共S peak兲 at about 145 cm−1, in 共⬇54 meV兲 is in good agreement with previous phenomeno-
addition to the broad band in the central part, which appears logical estimates. In contrast, the corresponding Slater pa-
when T is lowered below Tc.45 The S peak has been tradi- rameter 共⬇17 meV兲 is larger in our calculation, while the
tionally assigned to the proton tunneling mode in KDP at T correlated value per phosphate agrees well with the phenom-
⬍ Tc, although Raman measurements have discarded this enological estimates, suggesting that the latter include con-
assignment.44 Actually, the S peak shows up in the z共xy兲z̄ siderable effective correlations in their evaluation. We have
spectrum of DKDP at a higher frequency than in KDP, and also shown that the chain of Slater defects induces phosphate
hence this mode was finally assigned to the libration mode of rotations, which in turn produce a lattice contraction in the
PO4 tetrahedra, with A1 symmetry.44 However, it has been xy plane. These results establish, for the first time, a connec-
subsequently pointed out in Ref. 45 that the S peak cannot tion between both phenomena, i.e., the appearance of Slater
stem from a librational PO4 mode, because this type of mode defects and the basal lattice contraction, which were ob-
should not be affected by deuteration. We have shown in our served independently in different experiments near Tc.42,39
calculations that the appearance of Slater defect chains is
ACKNOWLEDGMENTS
accompanied by phosphate rotations. The symmetry of these
rotations is similar to that of the librational A1 mode. The This work was supported in part by the Agencia Nacional
fact that Slater defects are induced near Tc,42,16 and that they de Promoción Científica y Tecnológica and CONICET under
are coupled to phosphate rotations, as demonstrated in this Grant No. PICT99 Nr. 03-07248. R.M. and S.K. acknowl-
work, may then explain why the librational mode is deutera- edge additional support from CONICET. S.K. acknowledges
tion dependent. Moreover, this argument is further supported support from Florida State University 共FSU兲, Tallahassee,
by the fact that the phosphate rotations induced by Slater USA, and also the kind hospitality received during the stay at
defects finally lead to the lattice contraction observed experi- the Department of Chemistry and Biochemistry of FSU.
mentally near Tc,39,40 as was discussed in the previous sec- S.K. also acknowledges support from Fundación Antorchas,
tion. Interestingly, a strong coupling between the librational Argentina.

10 N.
*Corresponding author. E-mail address: [email protected] Dalal, A. Klymachyov, and A. Bussmann-Holder, Phys. Rev.
1 G. Busch and P. Scherrer, Naturwiss. 23, 737 共1935兲. Lett. 81, 5924 共1998兲.
2 J. C. Slater, J. Chem. Phys. 9, 16 共1941兲. 11 M. I. McMahon, R. J. Nelmes, W. F. Kuhs, R. Dorwarth, R. O.
3
V. H. Schmidt, Ferroelectrics 72, 157 共1987兲. Piltz, and Z. Tun, Nature 共London兲 348, 317 共1990兲.
4 Y. Takagi, J. Phys. Soc. Jpn. 3, 273 共1948兲. 12 R. J. Nelmes, M. I. McMahon, R. O. Piltz, and N. G. Wright,
5 M. E. Senko, Phys. Rev. 121, 1599 共1961兲. Ferroelectrics 124, 355 共1991兲.
6 H. B. Silsbee, E. A. Uehling, and V. H. Schmidt, Phys. Rev. 133, 13 J. A. Krumhansl, Nature 共London兲 348, 285 共1990兲.

A165 共1964兲. 14 H. Sugimoto and S. Ikeda, Phys. Rev. Lett. 67, 1306 共1991兲.
7
R. Blinc, J. Phys. Chem. Solids 13, 204 共1960兲. 15
A. Bussmann-Holder and K. H. Michel, Phys. Rev. Lett. 80, 2173
8 K. Kobayashi, J. Phys. Soc. Jpn. 24, 497 共1968兲.
共1998兲.
9 16
A. Bussmann-Holder, N. Dalal, R. Fu, and R. Migoni, J. Phys.: D. Merunka and B. Rakvin, Solid State Commun. 129, 375
Condens. Matter 13, L231 共2001兲. 共2004兲.

104104-7
LASAVE et al. PHYSICAL REVIEW B 72, 104104 共2005兲

17 D. Merunka and B. Rakvin, Chem. Phys. Lett. 393, 558 共2004兲. Phys. Rev. B 71, 184102 共2005兲.
18
V. H. Schmidt and E. A. Uehling, Phys. Rev. 126, 447 共1962兲. 32
G. F. Reiter, J. Mayers, and P. Platzman, Phys. Rev. Lett. 89,
19
V. H. Schmidt, Phys. Rev. 164, 749 共1967兲. 135505 共2002兲.
20 Y. N. Huang, X. Li, Y. Ding, Y. N. Wang, H. M. Shen, Z. F. 33 B. Rakvin and N. S. Dalal, Phys. Rev. B 44, R892 共1991兲.

Zhang, C. S. Fang, S. H. Zhuo, and P. C. W. Fung, Phys. Rev. B 34 P. Ordejón, E. Artacho, and J. M. Soler, Phys. Rev. B 53, R10441

55, 16159 共1997兲. 共1996兲; D. Sánchez-Portal, P. Ordejón, E. Artacho, and J. M.

21 V. H. Schmidt, G. Bohannan, D. Arbogast, and G. Tuthill, J. Phys. Soler, Int. J. Quantum Chem. 65, 453 共1997兲.
Chem. Solids 61, 283 共2000兲. 35
O. F. Sankey and D. J. Niklewski, Phys. Rev. B 40, 3979 共1989兲.
22
R. Blinc and S. Svetina, Phys. Rev. 147, 430 共1966兲. 36
N. Troullier and J. L. Martins, Phys. Rev. B 43, 1993 共1991兲.
23 W. Reese, Phys. Rev. 181, 905 共1969兲; C. W. Fairall and W. 37 J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77,

Reese, Phys. Rev. B 11, 2066 共1975兲. 3865 共1996兲.

24 S. Koval, J. Kohanoff, R. L. Migoni, and A. Bussmann-Holder, 38 R. J. Nelmes, Z. Tun, and W. F. Kuhs, Ferroelectrics 71, 125

Comput. Mater. Sci. 22, 87 共2001兲. 共1987兲.

25 39
S. Koval, J. Kohanoff, R. L. Migoni, and E. Tosatti, Phys. Rev. J. Kobayashi, Y. Mesu, I. Mizutani, and Y. Enomoto, Phys. Status
Lett. 89, 187602 共2002兲. Solidi A 3, 63 共1970兲.
26 Q. Zhang, F. Chen, N. Kioussis, S. G. Demos, and H. B. Ra- 40 A. Katrusiak, Phys. Rev. B 48, 2992 共1993兲.

dousky, Phys. Rev. B 65, 024108 共2001兲. 41 R. Blinc and B. Zeks, Ferroelectrics 72, 193 共1987兲.
27 S. Koval, J. Kohanoff, and R. L. Migoni, Ferroelectrics 268, 239 42 R. Blinc and J. L. Bjorkstam, Phys. Rev. Lett. 23, 788 共1969兲.

共2002兲. 43 M. Ichikawa, D. Amasaki, T. Gustafsson, and I. Olovsson, Phys.

28
C. S. Liu, N. Kioussis, S. G. Demos, and H. B. Radousky, Phys. Rev. B 64, R100101 共2001兲.
Rev. Lett. 91, 015505 共2003兲. 44 Y. Tominaga, M. Kashara, H. Urabe, and I. Tatsuzaki, Solid State
29
C. S. Liu, Q. Zhang, N. Kioussis, S. G. Demos, and H. B. Ra- Commun. 47, 835 共1983兲; Y. Tominaga, H. Urabe, and M.
dousky, Phys. Rev. B 68, 224107 共2003兲. Tokunaga, ibid. 48, 265 共1983兲.
30 G. Colizzi, J. Kohanoff, J. Lasave, S. Koval, and R. L. Migoni, 45 Y. Takagi, Ferroelectrics 72, 67 共1987兲.

Ferroelectrics 301, 61 共2004兲. 46 G. Colizzi, Ph.D. thesis, Queen’s University of Belfast, UK,
31 S. Koval, J. Kohanoff, J. Lasave, G. Colizzi, and R. L. Migoni,
2004.

104104-8