ARTICLE pubs.acs.org/JPCC
The Isotope-Effect in the Phase Transition of KH2PO4: New Insights from Ab Initio Path-Integral Simulations Varadharajan Srinivasan†,§ and Daniel Sebastiani*,‡,^ †
Dept. of Chemistry, IISER, Bhopal 462023, Madhya Pradesh, India Dahlem Center for Complex Quantum Systems, Free University Berlin, Arnimalle 14, 14195 Berlin, Germany § Department of Chemistry, Princeton University, Princeton, New Jersey 08544, United States ^ Max-Planck Institut f€ur Polymerforschung, 10 Ackermannweg, D 55128 Mainz, Germany ‡
ABSTRACT: We investigate the quantum-mechanical localization of 1H and 2H isotopes in the symmetric low-barrier hydrogenbonds of potassium dihydrogen phosphate (KDP) crystals in the paraelectric phase. The spatial density distributions of these hydrogen atoms are suspected to be responsible for the surprisingly large isotope effect observed for the ferroelectric phase transition in KDP. We employ ab initio path integral molecular dynamics simulations to obtain the nuclear real-space and momentum-space densities n(R) and n(k) of 1H and 2H, of which the latter densitites are compared to experimental neutron compton scattering data. Our results suggest a qualitative difference in the nature of the paraelectric phase in KDP between the two isotopes. Whereas both paraelectric states result from quantum delocalization, the essential difference is the change from a probably coherent to incoherent tunneling behavior of the hydrogen atoms across the hydrogenbonds.
’ INTRODUCTION Potassium dihydrogen phosphate (KDP) is the prototype of a wide class of hydrogen-bonded ferroelectrics that find use in nonlinear optics applications.1 Its phosphate units are interconnected by a network of hydrogen-bonds (H-bonds). The ferroelectric to paraelectric phase transition in these ionic crystals is associated with a shift from asymmetric OH 3 3 3 O H-bonds (Figure 1) to a symmetric configuration (on averge). Earlier theories for the phase transition associated the symmetry of the paraelectric phase with a classical disordering of protons over the two donor sites in each H-bond.24 However, the observation of a very large isotope effect (Tc[1H] = 122 K to Tc[2H] = 230 K5), strongly suggested the presence of quantum effects. As a result, the observed disorder was attributed6 solely to 1H-tunneling between the two equivalent donor sites in each H-bond.7,8 The increased mass due to deuteration reduces the tunnel-splitting. As the asymmetric OH 3 3 3 O bond is favored by electrostatic effects, it would then dominate the tunneling delocalization over a larger range of temperatures and lead to a higher Tc. This tunneling model assumed that the mass effect is dominant and that the effective potential energy surface (PES) for the 1 H/2H motion remains the same. However, pressure-dependent studies on KDP9 and related systems indicate10 that the tunneling hypothesis alone is insufficient to account for the observed increase in Tc. An additional Ubbelohde9 effect is assumed to contribute significantly to the increase in Tc. Nonetheless, the r 2011 American Chemical Society
effect is dominantly of quantum nature and therefore requires a full quantum description of the nuclei.11,12 In this article, we probe the nature of the paraelectric phases of protonated and deuterated KDP using a fully consistent quantum mechanical description of both nuclei and electrons. For the first time, real space and reciprocal (momentum-) space density distributions of the hydrogen nuclei are computed at realistic temperatures and in the actual periodic crystal structure, via a recently developed modification of the ab initio path-integral molecular dynamics (PIMD) method.1315 Our calculations show that the nuclear quantum delocalization of protons/deuterons strongly influences the effective potential felt by the nuclei and thus could significantly increase the isotope effect.
’ COMPUTATIONAL DETAILS Solving the nuclear Schr€odinger equation on a high-dimensional density-functional theory (DFT) based PES is currently unfeasible. Our calculations are based on the path integral technique, which represents a convenient way to map the nuclear quantum problem to a statistical ensemble average of a set of nonquantum replica of the same system. The general pathintegral method is described in detail elsewhere.1621 It has Received: March 18, 2011 Revised: May 26, 2011 Published: May 26, 2011 12631
dx.doi.org/10.1021/jp202584p | J. Phys. Chem. C 2011, 115, 12631–12635
The Journal of Physical Chemistry C
ARTICLE
^ and ÆR|kæ = eikR, we obtain energy operator T Z nnuc ðkÞ ¼ eðβ=PÞT ðkÞ d3n R ð1Þ d3n R ð2Þ eik 3 ðR
ð1Þ
R ð2Þ Þ
FðR ð1Þ , R ð2Þ Þ
ð4Þ
0
Figure 1. Left: The paraelectric KDP unit cell along the polar c axis. The H-bonds between phosphates lie along the a and b axes, respectively. Right: two phosphate molecules with their potassium counterions that form the basic building block of a KDP unit cell. The ferroelectric polarization is caused by a symmetry breaking in the KPK distances along the c axis. Bottom panel with a single OH 3 3 3 O bond illustrates the proton-localization coordinate δ and the acceptordonor distance 2R (text).
already been applied numerous times to obtain valuable insight into fundamentally and technologically relevant materials.2226 In the common path integral formulation of a quantum system of n nuclei, the canonical partition function Z = Tr exp {βH^ } = ^[T] is written as an integral over discretized paths Tr F {R(1)R(2)...R(P); R(i) ∈R3n} in coordinate-space, 8 !9 Z < βX = P 2 2 mP Z d3nP R exp V ðR ðpÞ Þ þ 2 2 ΔR ðpÞ : P p¼1 ; 2p β ð1Þ where β1 = kBT and P ∈ N is a large integer, ΔR(p) = R(p) R(pþ1) and R(Pþ1) = R(1). This partition function is equivalent to that of n (classical) ring polymers (one per particle), each composed of P elements (beads) connected by classical springs with spring constant D = (mP2)/(p2β2). The ensemble of n particles for a given bead is called a replica. The particles within the pth replica interact through the potential V (R(p)), which in our case is the total electronic energy. Every replica is a classical representation to the quantum system under study, whereas the ensemble of the P coupled replica is isomorphic to the quantum system of n atoms. The isomorphism is exploited to compute the solution of the n-particle quantum problem via a thermodynamical phase space sampling of the n coupled classical ring polymers. This sampling can be done by means of Monte Carlo techniques or molecular dynamics (MD) simulations, our choice being the latter (using the CPMD code27). Expectation values of all quantum observables, which are local in direct space can be expressed via the statistical density of beads, in particular the probability density F ½T Þ nnuc ðRÞ ¼ TrðjRæÆRj^ ð1Þ
^ Þæring polymer ¼Æδ ^ 3n ðR R
ð2Þ
ð3Þ
Properties that are represented in the space of the momenta p = pk, however, require the momentum density nnuc(k) = Tr (|kæÆk|^ F[T]), and thus a modification to the original pathintegral formulation. Since |kæ is an eigenstate of the kinetic
with the off-diagonal density matrix element F(R,R ). The latter can be sampled within the classical path-integral isomorphism via a open polymer, where the classical spring between replica R(1) and R(2) has been eliminated.21 Using this open path integral formalism at P = 32 within a gradient-corrected DFT framework,28,29 we have computed the direct- and reciprocal space density distributions n(r) and n(k) of the hydrogen atoms in KDP and DKDP crystals by opening the paths for only these atoms.42 The simulation temperatures, that is 140 and 300 K, were chosen in the para-electric sectors of the experimental phase diagrams of KDP and DKDP, respectively. Our supercell for the paraelectric phases contains four independent phosphate units on body-centered tetragonal sites (8 hydrogens per cell). Protontunneling could be coupled to collective motion of protons with wavelengths of several unit cells;11,12,30 our supercell describes only phonons of infinite wavelength. The lattice parameters for each isotope were taken to be the experimentally measured values in the paraelectric phases close to the Tc. For KDP a = 7.4264 Å, c = 6.931 Å, and for DKDP a = 7.459 Å, c = 6.957 Å. Note that the difference in volume is only about 1.3%. The total simulated time was about 10 ps (at a 0.1 fs time-step) ensuring good statistics. We believe that this computational setup yields a reasonable statistical convergence of the path integral molecular dynamics trajectory and its derived properties. It is clear, however, that an explicit computation of more subtle features in the momentum distribution would require enhanced sampling. This particularly applies to the intensely discussed existence of a node in n(k),11,12 which could indicate a truly bimodal proton distribution.
’ RESULTS As illustrated in Figure 1, the KDP crystal contains a network of OH 3 3 3 O hydrogen bonds between the phosphate units, orthogonal to the c axis along which the electric polarization is oriented in the ferroelectric phase. We characterize the H-bonds by the oxygenoxygen distance 2R and a hydrogen-localization coordinate δ = R ROl H 3 ROl Ou which quantifies the asymmetry of the OH 3 3 3 O bond. Here, the subscripts l(u) stand for the lower (upper) oxygen, ROl H refers to the instantaneous vector from the lower oxygen to the hydrogen, and ROl Ou is the unit vector pointing from the lower to the upper oxygen. Each MD step provides a set of P values for δ, whose statistical distribution n(δ) over the PIMD run defines a hydrogenlocalization function (HLF). The HLF represents the projection of the real-space density of the hydrogens along the OH 3 3 3 O bond. A negative (positive) value of δ implies that the hydrogen is closer to the lower (upper) oxygen Ol (Ou) of the H-bond (Figure 1). Part a of Figure 2 shows the HLF for the 8 protons in KDP resulting from the PIMD simulations. The broad HLF can result from either quantum delocalization or a classical hopping across a barrier in the H-bond. However, a conventional Car Parrinello (CP) simulation at the same conditions showed no hopping at all (insert in part a of Figure 2) but, instead, resulted in a symmetry-broken ferroelectric state. Hence, in the time-scale of 12632
dx.doi.org/10.1021/jp202584p |J. Phys. Chem. C 2011, 115, 12631–12635
The Journal of Physical Chemistry C
ARTICLE
Figure 2. (a) Individual hydrogen localization functions HLF (left) and 2R-δ correlation (right, averaged over all hydrogens in the simulation cell) from 10 ps of ab initio path-integral MD simulations for 1H-KDP at T = 140 K (top) and 2H-KDP at 300 K (bottom) (paraelectric phases). The insert shows the HLF from a conventional MD simulation of the same system; (b) Average momentum distributions computed and experimental (NCS).7,8.
Figure 3. (a) Computed momentum and real-space (insert) densities of several deuterons in DKDP at T = 300 K, and their spreads in Å (insert); (b) Evolution of the n[2H](δ) of one specific 2H in subsequent 1 ps windows; and (c) the computed spreads Æk2æ1 of n(k) for 1H (left) and 2H (right) for a free particle, a simple harmonic oscillator approximation (SHOA) and KDP at different temperatures.
this simulation, a classical proton cannot overcome the OH 3 3 3 O barrier in KDP at 140 K. This implies that the broad HLFs are due to quantum delocalization of the protons in the H-bonds, leading to the symmetric paraelectric phase. The proton distribution is apparently unimodal suggesting that the symmetrization could be possibly due to zero-point motion broadening and not tunneling. The present series of calculations does not allow an unquestionable conclusion on the long-debated issue of a unimodal versus a bimodal proton density distribution.31 The accuracy of the underlying potential energy surface, the approximations involved in the implementation of the path integral isomorphism, and also the convergence of the actual path integral MD simulations should be scrutinized further. However, our simulations can provide a piece of the mosaic, which can be combined with accurate measurements of structural data of highpressure and variable temperature KDP systems.3236 The computed HLFs of the 8 deuterons in paraelectric DKDP at 300 K are shown in the bottom left of part a of Figure 2. In contrast to KDP, we see a diversification of the HLF ranging from well-localized off-centered deuterons to considerably delocalized ones. An overall average in them leads to a symmetric bimodal distribution. Thus, unlike in KDP, this symmetry results from a distribution of more localized deuterons between the two OH 3 3 3 O equilibrium sites. As complementary fingerprints of
the paraelectric KDP (DKDP), the 1H (2H) nuclear momentum distributions n(k) along the OH 3 3 3 O bond are shown in part b of Figure 2 along with the corresponding experimental data from neutron Compton scattering experiments.7,8,37,38 The qualitative agreement with experiment confirms the 1H -tunneling picture in KDP. The n(k) of the deuterons is broader, corresponding to a tighter real-space localization. In fact, the experimental n[2H](k) corresponds to an even higher degree of spatial localization. It should be pointed out that, whereas our results reproduce the experimentally observed trend in localization we do not observe any nodes in the n[1H](k) as seen in the experiment.7,8,38 This could be due to a well-known issue with DFT-GGA, which tends to underestimate proton-transfer barriers in hydrogen-bonded systems. We have further computed σ = (Æk2æ1)1/2 as a measure of the instantaneous real-space hydrogen localization (σ scale-bar in part c of Figure 3). Whereas the deuterons are more localized than the protons, the spread of either species lies between that of a free particle and a simple harmonic oscillator (SHO), which is shown in part c of Figure 3 as a simple model for a (slightly red-shifted) OH stretching mode of 2500 cm1. The small difference in n(k) between 2H and 1H indicates that neither one is a classical particle. Interestingly, deuteration at the temperature and volume conditions of the paraelectric KDP simulation resulted in only 12633
dx.doi.org/10.1021/jp202584p |J. Phys. Chem. C 2011, 115, 12631–12635
The Journal of Physical Chemistry C partial localization (DKDP at 140 K). The localization is enhanced upon increasing the temperature. The geometric effect manifests itself in experiments as an increased OH 3 3 3 O distance (2R in Figure 1) in DKDP with respect to KDP in the paraelectric phase.9 Our PIMD simulations reproduced this experimental trend with an average 2R of 2.49 Å in KDP and 2.60 Å in DKDP. Furthermore, we observed a correlation between the instantaneous 2R values and the localization δ of the hydrogens, shown in part a of Figure 2. At shorter 2R, the hydrogen prefers the center of the H-bond, whereas for larger 2R the most-probable positions of the hydrogen are off-center. In fact, for very short H-bonds (2R < 2.45 Å) the off-center peaks are absent suggesting the disappearance of a symmetry-broken phase at high pressures.39 Thus, the probability of finding the hydrogen at the center of the H-bond is effectively modulated by the acceptordonor distance 2R, which determines the tunneling barrier.40,41 The symmetric tunneling regime corresponds to a rather thin range of 2R values where the δ coordinate has a distinctively larger spread about the center. This range (occurring around 2.49 Å) is wider for KDP than DKDP due to the higher mass of the 2H. We characterize structural fluctuations in paraelectric DKDP by means of momentum distributions for the individual 2H as well as the corresponding real-space localization functions, shown in Figure 3. The real-space picture shows that there are both centered and off-centered deuterons; the momentum distributions, however, all have very similar widths. In some cases, deuterons with a broader real-space distribution (e.g., the green line in part a of Figure 3) exhibit an even broader momentumspace distribution than relatively more localized ones (e.g., the blue line in part a of Figure 3). This apparent inconsistency is resolved by following the evolution of the real-space density of the hydrogen plotted in green (part b of Figure 3). Initially, the 2H is localized at negative δ values (black, red, violet lines) but it becomes delocalized after a few picoseconds (blue line) and subsequently localizes at positive δ (yellow and brown lines). At the end of our run, further delocalization starts (orange line). This plot illustrates that the 2 H is comparably localized most of the time but the association to the oxygens may change over time. We note, however, that this evolution does not represent the physical time evolution of the 2H quantum state. It merely allows us to monitor the correlation between the spread of a 2H and the OO distance in its hydrogen-bond. As the 2H moves toward the bond-center, the corresponding 2R also decreases (correlation plot in part a of Figure 2). This not only reduces the barrier to cross the bond-center but also leads to a larger real-space width in n[2H](δ).40 Thus, the HLF is the widest as the 2H passes through the center of the H-bond. Beyond this point, 2R increases again simultaneously localizing the 2H . Hence, the oscillations in the acceptordonor distance 2R, effected by rotational/vibrational modes of the phosphate groups, are ultimately responsible for the dynamic disordering of the deuterons in paraelectric DKDP. The symmetrization in the paraelectric phase of DKDP is thus brought about by an incoherent tunneling mechanism via a coupling to the modes of the lattice.
’ CONCLUSIONS We find that the paraelectric phases in KDP and DKDP are stabilized by proton quantum delocalization in the OH 3 3 3 O
ARTICLE
hydrogen bonds. The essential difference in the nature of the two isotopic crystals is that, whereas in KDP a symmetric state arises possibly due to a coherent state, in DKDP incoherent tunneling of the deuterons effected by vibrations of the lattice are responsible for the symmetrization. The isotope effect then follows from the need for thermal population of the lattice modes to enable dynamical disordering in DKDP. In other words, the so-called geometric effect is triggered by the quantum properties of the hydrogen, which is different between 1H and 2H.
’ AUTHOR INFORMATION Corresponding Author
*E-mail:
[email protected].
’ ACKNOWLEDGMENT We acknowledge helpful discussions with George Reiter, Roberto Car, and Joe Morrone. This work was financially supported by the DFG under grants Se 1008/5-1 and Se 1008/6-1, as well as a IMPRS scholarship from the Max-Planck-Institute for Polymer Research. ’ REFERENCES (1) Lines, M. E.; Glass, A. M. Principles and Applications of Ferroelectrics and Related Materials; Clarendon Press: Oxford, 1977. (2) Slater, J. C. J. Chem. Phys. 1941, 9, 16. (3) Takagi, Y. J. Phys. Soc. Jpn. 1948, 3, 271–272. (4) Takagi, Y. J. Phys. Soc. Jpn. 1948, 3, 273–274. (5) Bantle, W. Helv. Phys. Acta 1942, 15, 373. (6) Blinc, R. J. Phys. Chem. Solids 1960, 13, 204–211. (7) Reiter, G. F.; Mayers, J.; Platzman, P. Phys. Rev. Lett. 2002, 89, 135505. (8) Reiter, G. F.; Mayers, J.; Platzman, P. Ferroelectrics 2002, 268, 283–288. (9) Nelmes, R. J. Ferroelectrics 1987, 71, 87–123. (10) McMahon, M. I.; Nelmes, R. J.; Kuhs, W. F.; Dorwarth, R.; Piltz, R. O.; Tun, Z. Nature 1990, 348, 317–319. (11) Koval, S.; Lasave, J.; Kohanoff, J.; Migoni, R. Ferroelectrics 2010, 401, 103–109. (12) Colizzi, G.; Kohanoff, J.; Lasave, J.; Migoni, R. Ferroelectrics 2010, 401, 200–206. (13) Burnham, C. J.; Reiter, G. F.; Mayers, J.; Abdul-Redah, T.; Reichert, H.; Dosch, H. Phys. Chem. Chem. Phys. 2006, 8, 3966–3977. (14) Morrone, J.; Srinivasan, V.; Sebastiani, D.; Car, R. J. Chem. Phys. 2007, 126, 234504. (15) Morrone, J. A.; Car, R. Phys. Rev. Lett. 2008, 101, 017801. (16) Feynman, R. P.; Hibbs, A. R. Quantum Mechanics and Path Integrals; McGraw-Hill: New York, 1965. (17) Tuckerman, M. E.; Berne, B. J.; Martyna, G. J.; Klein, M. L. J. Chem. Phys. 1993, 99, 2796–2808. (18) Tuckerman, M.; Marx, D.; Klein, M.; Parrinello, M. J. Chem. Phys. 1996, 104, 5579–5588. (19) Cao, J. S.; Voth, G. A. J. Chem. Phys. 1994, 100, 5093–5105. (20) Voth, G. A. Adv. Chem. Phys. 1996, 93, 135–218. (21) Ceperley, D. Rev. Mod. Phys. 1995, 67, 279–355. (22) von Rosenvinge, T.; Tuckerman, M. E.; Klein, M. L. Faraday Discuss. 1997, 106, 273–289. (23) Hayes, R. L.; Paddison, S. J.; Tuckerman, M. E. J. Phys. Chem. B 2009, 113, 16574–16589. (24) Ceriotti, M.; Miceli, G.; Pietropaolo, A.; Colognesi, D.; Nale, A.; Catti, M.; Bernasconi, M.; Parrinello, M. Phys. Rev. B 2010, 82, 174306. (25) Ludue~ na, G. A.; Wegner, M.; Bjalie, L.; Sebastiani, D. ChemPhysChem 2010, 11, 2353–2360. 12634
dx.doi.org/10.1021/jp202584p |J. Phys. Chem. C 2011, 115, 12631–12635
The Journal of Physical Chemistry C
ARTICLE
(26) Luduena, G. A.; Sebastiani, D. J. Phys. Chem. Lett. 2010, 1, 3214–3218. (27) Hutter, J. Computer code CPMD, version 3.12; Copyright IBM Corp. and MPI-FKF: Stuttgart, 1990-2008, www.cpmd.org. (28) Becke, A. D. Phys. Rev. A 1988, 38, 3098. (29) Lee, C.; Yang, W.; Parr, R. G. Phys. Rev. B 1988, 37, 785–789. (30) Koval, S.; Kohanoff, J.; Lasave, J.; Collizi, G.; Migoni, R. L. Phys. Rev. B 2005, 71, 184102. (31) Kuhs, W.; Nelmes, R.; Tibballs, J. J. Phys. C: Solid State Phys. 1983, 16, 1029–1032. (32) Nelmes, R.; Tun, Z.; Kuhs, W. Ferroelectrics 1987, 71, 125–141. (33) Nelmes, R.; McMahon, M.; Piltz, R.; Wright, N. Ferroelectrics 1991, 124, 355–360. (34) Tibballs, J.; Nelmes, R. J. Phys. C: Solid State Phys. 1982, 15, L849–L853. (35) Tun, Z.; Nelmes, R. J.; Kuhs, W. F.; Stansfield, R. F. D. J. Phys. C: Solid State Phys. 1988, 21, 245–258. (36) Nelmes, R. J.; Meyer, G. M.; Tibbals, J. E. J. Phys. C: Solid State Phys. 1982, 15, 59–75. (37) Andreani, C.; Colognesi, D.; Mayers, J.; Reiter, G. F.; Senesi, R. Adv. Phys. 2005, 54, 377–469. (38) Reiter, G. F.; Shukla, A.; Platzman, P. M.; Mayers, J. New J. Phys. 2008, 10, 013016. (39) Samara, G. A. Phys. Rev. Lett. 1971, 27, 103–106. (40) Koval, S.; Kohanoff, J.; Migoni, R. L.; Tosatti, E. Phys. Rev. Lett. 2002, 89, 187602. (41) Zhang, Q.; Chen, F.; Kioussis, N.; Demos, S. G.; Radhousky, H. B. Phys. Rev. B 2002, 65, 024108. (42) All hydrogen atoms in a cell were represented by open paths neglecting any direct HH interactions. Because the H-bonds in this system are well-separated, this is a good approximation. See ref 9 for other path-opening strategies in more complicated cases.
12635
dx.doi.org/10.1021/jp202584p |J. Phys. Chem. C 2011, 115, 12631–12635
