Path integral Monte Carlo simulations of the geometrical effects in KDP crystals

Path integral Monte Carlo (PIMC) simulations with very simple models were used in order to unveil the physics behind the isotope effects in H-bonded ferroelectrics. First, we studied geometrical effects in the H-bonds caused by deuteration with a general three-site model based on a back-to-back double Morse potential plus a Morse potential between oxygens, fitted to explain different general features for a wide set of H-bonded compounds. Our model results show the Ubbelohde or geometrical effect (GE), i.e., the expansion of the H-bond with deuteration, in agreement to what is observed in H-bonded ferroelectrics with short H-bonds. Moreover, adjusting the potential parameters to ab initio results, we have developed a 1D model which considers the bilinear proton-proton interaction in mean-field to study nuclear quantum effects that give rise to the GE in KDP crystals. PIMC simulations reveal that protons tunnel more efficiently than deuterons along the 1D chain, giving rise to a strong attraction center that pulls the oxygens together. This mechanism, which is based on the correlation between tunneling and geometrial modifications of the H-bonds, leads to a strong GE in the ordered phase of the chain at low temperature which is in good agreement with the experimental data.


Introduction
or KDP is the prototype of a wide family of H-bonded ferroelectric compounds which has extensive applications as a key component in optoelectronic devices [1]. Besides the technological interest, KDP has also attracted much attention due to its rich, complex and intriguing phenomenology, e.g., the huge isotope effect that displays associated to its ferroelectric-paraelectric (FE-PE) phase transition. With deuteration, the critical temperature changes from ≈ 122 K to ≈ 210 K. The saturated polarization at low also shows a large isotope effect, increasing from ≈ 5.0 µC/cm 2 for KDP to ≈ 6.2 µC/cm 2 for a sample with 98% of deuteration [2].
The origin of these strong isotope effects is still controversial. The first explanation of the large increase of upon deuteration was given by the quantum tunneling model [3], which focuses purely on mass-dependent effects. However, increasing experimental evidence since the late nineteen eighties showed that the large isotope effect is mainly driven by geometrical modifications of the H bonds [4,5] (Ubbelohde effect [6]). The recent observation of tunneling in the PE phase of KDP by neutron Compton scattering experiments added even more controversy to the problem [7], although in deuterated KDP (DKDP), tunneling could not be detected [8].
Ab initio calculations have recently shown that tunneling and geometric effects are complementary aspects of the same phenomenon [9,10]. With a simple selfconsistent model based on ab initio results, it is demonstrated that the wave function solution of the nonlinear Schrödinger equation for deuteron/proton clusters evolves from a double peak to a broad single peak located at the center of the H-bonds as the cluster mass diminishes. This is explained by a strong nonlinear feedback between proton delocalization (tunneling) and the effective proton potential barrier in the H-bonds, which changes concomitantly with the H-bond geometry. It is concluded that such a large mass dependence can explain the large isotope effect found in KDP, via an amplified and selfconsistent geometric modification of the H bond in agreement with experiments. On the other hand, these results are in striking contrast with the very weak dependence obtained at fixed potential and geometry. Thus, the proton tunneling subunit and the host lattice are strongly coupled and the host-and-tunneling system is not separable.
Many models were successfully developed in the past to shed light into the general phenomenology of H-bonded ferroelectric materials [11][12][13][14][15][16][17][18]. In this paper, we address with very simple models the problem of geometrical effects in KDP crystals by performing path integral Monte Carlo (PIMC) simulations. First, we develop a three-site model for the H-bond to study local quantum geometric effects. This simple model already serves us to gain knowledge about the interplay between proton tunneling and H-bond geometric modifications such as the O-O distance variation. After this first insight, we develop a 1D chain model of concatenated H-bonds to study in the ordered phase the geometrical effects caused by deuteration. The model parameters are fitted using recent ab initio results [19]. We demonstrate that this simple linear model can account for the geometrical effects observed in real H-bonded ferroelectrics, which are at the root of the giant isotope effect in the critical temperature observed in the FE phase transitions of these materials. The paper is organized as follows: in the next section we explain the models used and describe details of the PIMC calculations. Section 3 describes and discusses the results obtained for the three-site model and for the linear chain. Finally, we elaborate a summary and our conclussions in section 4. three-site model. ≡ OO is the distance between oxygen nuclei. OH is the proton-oxygen distance. The variable = OO − 2 OH is defined as the distance between the two possible equilibrium positions of the proton. Then, = OO /2 − OH is the proton coordinate relative to the H-bond center. This parameter definition is also used in the linear chain model.

Three-site model
We developed a three-site (3S) model which represents a single O-H-O cluster embedded in the Hbonded ferroelectric as it is sketched in figure 1. With the aim to model linear H-bonds, a Double Morse (or back-to-back) potential (see e.g., [20][21][22][23][24]) is usually used, which is essentially the superposition of two Morse potentials representing what the proton feels while interacting with both oxygens: where is the O-O distance, and represents the H position relative to the H-bridge center (see figure 1). If we assume that is fixed, there is a critical value = 2( −1 ln 2+ 0 ) such that for < the potential profile is a single well with a minimum at = 0. On the contrary, for > we have a symmetric doublewell potential, with a local maximum at = 0 and minima at = Notice that the energy barrier for the proton jump from one side to the other of the H-bond diminishes concomitantly with the O-O distance , vanishing for < . Actually, we are interested in the 43708-2 proton/deuteron tunneling regime, thus we would need that the equilibrium distance remains in the region where the proton barrier exists, that is > . However, simulations at low temperature with the potential described in equation 2.1, relaxing both variables and , yield to a collapse of the potential barrier and the equilibrium energy profile displays one minimum only. Therefore, it is mandatory to introduce a new interaction which preserves the system from the O-O distance collapse. This O-O potential will represent the interaction between both oxygens and the lattice. The following Morse potential between oxygens is chosen [19]: We adopted a Morse potential to describe the O-O interaction with the lattice because this kind of anharmonic potential enables the system to explore with sufficient probability O-O distances larger than 0 , in such a way that the collapse tendency to a single well is drastically diminished. This is in contrast to the case of a harmonic potential for the O-O interaction, where in this case the O-O collapse is inevitable. The complete potential for the 3S model is as follows: This correlation is precisely the important ingredient necessary for the existence of the Ubbelohde or the geometrical effect observed in compounds with strong H-bonds.

1D model of concatenated H-bonds
Going a step beyond the simple three-site model, we have developed a one dimensional chain model of concatenated H-bonds to study the GE in a more realistic way in the ordered phase. This . The supercell of dimension = 200 is subjected to periodic boundary conditions. In the simulation, is allowed to relax at zero stress, as well as each coordinate and of each unit cell . For instance, this chain represents a model approximation to the 1D H-bonded ferroelectric CsH 2 PO 4 (CDP) if the model chain oxygen is interpreted as a PO 4 unit plus an ordered hydrogen covalently bonded to the phosphate at any temperature, and the model hydrogen is the one that is disordered at high temperature in CDP [25]. Then, the global motion of hydrogens in our linear model in the ordered phase, from one minimum to the other along the H-bonds of the chain, could be related to the FE mode that accounts for the spontaneous polarization arising along the direction at low in CDP [25]. Alternatively, the chain model may also represent an approximation to the study of the GE in KH 2 PO 4 (KDP) if the model effective oxygen now represents a KDP cluster of two phosphate units including seven protons moving coordinately as a local FE mode [9,10]. In all these cases, we must adopt a convenient effective mass for the effective model hydrogen/deuteron considering that the real displacements of H(D) are accompanied with the heavier atom motions [9,10,19].
The total potential energy for the linear chain (1D) model is defined as: where 3 is the unit cell local potential defined exactly in the same way for the 3S model, as is shown in equation (2.3), and the last term is the short-range interaction energy between protons/deuterons 43708-3 stemming from the ice rules restrictions, i.e., in this 1D model, only one proton is attached to each oxygen. The last sum in equation 2.4 is restricted to nearest neighbours for each index . There is no long-range part in this model, which precludes a phase transition in one dimension. However, the last bilinear term is treated in mean-field, which enables the system to have a second order phase transition at finite temperature [26]. Therefore, the 1D model total potential, is written in the following way [27]: where ≡ 1/ is the time and lattice average of the positions for each unit cell taken at each MC step in the simulation.

Path integral Monte Carlo simulations
In the PIMC simulations [28], the effective short-time propagator for two adjacent points in the discretized imaginary-time path describing each quantum particle was evaluated to fourth-order accuracy with the Takahashi-Imada approximation [28][29][30]. The effective action in this case allows us to significantly reduce the Trotter number required for convergence. In all the simulations performed we have used = 128 beads for the quantum polymer associated with each atom in the O-H...O bonds, which yielded well-converged results [19,25,28]. Additionally, a normal-mode representation of the quantum polymers was used in order to ensure ergodicity in the MC sampling [28,30]. The PIMC simulations were performed at low = 50 K such that the quantum nuclear effects were predominant compared to entropic contributions in the 3S model and also with the aim to obtain GE in the ordered phase for the 1D model (the classical version of this model has a transition to a disordered paraelectric phase at ≈ 350 K). The simulations for the 3S model consisted of 1 × 10 6 MC steps preceded by 5 × 10 5 steps of thermalization.
In the 1D chain model simulations, we took 3 × 10 4 steps of thermalization plus 1 × 10 5 MC steps for computing averages. In this case, each calculation performed was an average of 20 runs with different random number generator seeds.
To characterize the degree of particle delocalization in the PIMC simulations, we studied the centroid and radius of gyration (RG) distributions for the quantum polymers [31]. The centroid is defined as the center of mass (CM) of the polymer and represents the average position of the quantum particle. The radius of gyration represents the variance of the quantum path and is a quantitative measure of how far away are the beads or monomers from the polymer center, and therefore, provides a measure of the quantum delocalization of the particle [31].

Geometrical effect study using the three-site model
The six potential parameters of equation (   different H-bonded compounds. There is a strong correlation between the OH and OO distances for the family of H-bonded compounds. The equilibrium distance OH diminishes systematically with increasing for > [33,34], reaching a saturated value around ∞ OH ≈ 0.95 Å for very large . Therefore, we took the parameter value 0 = 0.93 Å so that the values that minimize OH ( , ) in equation (2.1) for different values of give a curve min OH = OO /2 − min as a function of that is a lower bound for the set of experimental points spread in the OH-OO correlation [20,21,33,34]. With this choice, when the nuclear quantum effects are included in the PIMC calculations, we observe a very good agreement with the experimental correlation curve using the model of equation (2.1) with the OO distance fixed [35].
On the other hand, the parameter values for the OO interaction OO ( ) [see equation (2.2)], were initially taken from reference [23]. They were further adjusted, especially the value of OO , due to the important correlation between OH and OO , such that the classic potential profile has the minimum at  We have verified that the 3S-model PIMC simulations performed at = 50 K with = 128 beads for the quantum polymer representing each atom yielded probability distributions for the H-bond parameters ( and ) and energies well converged. The low temperature of 50 K for the simulation was chosen because we are interested in the nuclear quantum effects for the H-bonds and the geometrical changes with deuteration without most of the influence of entropic contributions in the particle dynamics. The 3S model results for the probability density contours to find the system in a given ( , ) configuration are shown in figure 3 for the proton and deuteron cases. The curves are qualitatively different but both cases are found to have symmetric distributions around = 0 in the coordinate with two prominent peaks with maximum probability, which are clearly shifted in the deuterated case. The OO distance for the peak positions are in each case:   [36,37]. Thus, our simple 3S model satisfactorily reproduces the isotopic geometrical effects for these systems.
It is worth to notice that if the OO distance is not allowed to relax, then the GE is smaller. For instance, we have fixed the value OO = 2.527 Å, which corresponds to the peak in the probability distribution for the protonic system (see figure 3), and the simulations gave a change with deuteration in the OH bond of only Δ = 0.021 Å. Comparing this result with that considering the oxygen dynamics (Δ = OH − OD = 0.032 Å), we observe an increment of ≈ 50% in the isotopic geometrical effect in the case where the oxygens are allowed to relax. This can be understood in the following way: first, when the oxygens are fixed, protons, being more delocalized than deuterons, have more probability to stay closer to the middle of the O-O bond. Second, when the oxygen dynamics is included, the protons act as a strong attraction center that pulls the two bridge oxygens together, more effectively than deuterons which are more localized near the oxygen. This proton-mediated O-O contraction lowers the potential barrier, which delocalizes even more the proton, and so on, giving rise to a nonlinear selfconsistent mechanism [9,10]. For the deuteron, being less delocalized than the proton, the selfconsistent effect is weaker. This mechanism leads to an isotopic geometrical effect which is stronger than that generated by the proton/deuteron quantum delocalization at fixed potential (fixed oxygens) [9,10].
To further illustrate the microscopic mechanism that rules the GE, we have analyzed the behavior of the quantum polymers for the proton/deuteron in the simulation via an analysis of the center of mass of the quantum polymer or centroid position and the radius of gyration representing a measure of the quantum delocalization of the particle (i.e., the extension of the quantum polymer) [31]. We plot in figure 4 the instantaneous values of as a function of the proton/deuteron centroids , taken every 100 MC steps in the PIMC simulation. As can be seen in the figure, the density of points reveals that the deuteron prefers to be localized at both sides and far from the bond middle with small values of , indicating a more classical behavior in these cases. When the deuteron centroid takes the values of closer to 0 (the bond middle), it is observed an increase of indicating that the quantum polymer is delocalized and is spread through both sides of the potential barrier, signaling the presence of tunneling in this case. Notice that the largest values of are found at ≈ 0 where delocalization is maximum. On the other hand, in the proton case, tunneling is much more frequent because the region with larger density of points appears near ≈ 0 with large values of , as shown in figure 4. This is precisely 43708-6 an important ingredient for the GE: the proton spends much more time delocalized with the quantum polymer center of mass near the middle of the O-O bond, which finally produces a strong contraction of the O-O distance. On the contrary, the deuteron is much more localized at both sides and far from the bond middle which leads to a weakening of the O-O bond and to an increase of the O-O distance. This yields the isotopic geometrical effect, which is observed in the calculated probability distribution of figure 3.

Isotope effects obtained with the 1D model simulations
The previous analysis of the 3S model results, which has clearly shown the isotopic GE, was carried out based on the parametrization of the model which reproduces the universal OH-OO correlation observed for a family of diverse H-bonded compounds. In this sense, this model is quite simple and general, accounting for the geometrical effects with deuteration of a set of H-bonded ferroelectrics with strong H-bonds. Now, we focus on the development of a 1D chain model, described in section 2.2 [see equation (2.5)], which was specifically designed to explain the isotope effects in the phase transition of KDP and was fitted to ab initio results [19]. This more realistic 1D model has, in the classical nuclei version, a ferroelectric-paraelectric transition at ≈ 350 K [35]. In this paper, we have used it in the ordered phase of KDP at = 50 K to analyze the isotopic GE which is at the root of the microscopic mechanism that leads to the giant isotope effect in the critical temperature.
We start from equation 2.5 for the 1D model, which has seven parameters to be adjusted for the KDP case. The six model parameters of the local proton potential 3 for each unit cell in the chain, which is just the same that was used in the 3S model (see equation 2.3), have been adjusted to reproduce six magnitudes obtained from ab initio calculations for KDP. These magnitudes were the global energy barrier between the PE and FE states, the O-O and distances in the FE phase, the O-O distance in the PE phase, the ab initio vibrational frequency of the PO 4 rotation mode, which is equivalent to the stretching mode in the 3S model, and the energy barrier between the energy minimum and the transition state in the FE phase keeping the O-O distance fixed (see reference [19]). We adopted the model fit to the ab initio calculations that includes dispersion corrections at the vdW-DF level, which exhibit, compared to other methods, the best agreement with the experimental geometry for both KDP and deuterated KDP (DKDP) [19].
Finally, we have fitted the remaining parameter that corresponds to the proton-proton interaction term in equations (2.4) and (2.5). To this end, was adjusted to 0.55 eV/Å 2 so that the critical temperature for the FE-PE transition obtained by the 1D model simulation with classical nuclei reaches the value of ≈ 350 K, similar to the value obtained by ab initio molecular dynamics calculations with dispersion corrections at the vdW-DF level for DKDP [38].
The final values for the parameters used in the 1D model are listed in table 2. The motion of the proton/deuteron is strongly correlated with that of the heavy ions, and its mass is dressed accordingly as discussed in reference [10]. Therefore, instead of using the bare proton (deuteron) masses (2 ), we have used in the PIMC simulations the effective masses for H and D: = 2.3 and = 3 , respectively, with the proton mass [9,10,19]. We plot in figure 5 the probability distribution contours for the PIMC simulation with the 1D model, obtained in the ordered phase at = 50 K. Due to the ordered phase, only one peak is observed in the proton and deuteron distributions, which is in contrast to the symmetrical double peaks around = 0 found in the 3S model distribution results (see figure 3). The calculated distribution for the chain of protons in figure 5 is asymmetric around the peak position due to the potential anharmonicity and quantum delocalization, which is in qualitative agreement with the experimental diffraction pattern measured near 43708-7 in the FE phase of KDP [39]. The asymmetry around the peak is less pronounced in the deuterated case as shown in figure 5, because the deuteron is less delocalized than the proton.
The prominent single peak found in the distribution results for the 1D simulation is clearly shifted in the deuterated case towards larger and , revealing the existence of the isotopic geometrical effect, i.e., the expansion of the H-bonds in the chain with deuteration.  table 3 and compared with the available experimental data for KDP and DKDP [40]. We observe a good agreement between theory and experiment, although the GE is a little bit underestimated, with difference values under deuteration ≈ 25% lower than the experimental data. Table 3. Nuclear quantum calculations of the H-bond geometries for KDP and DKDP using the 1D linear model. The results, which correspond to the peak positions of figure 5, are contrasted with the experimental data of reference [40]. Distances are in Å. To get a deeper insight into the microscopic mechanism of the geometrical effect in the linear chain model, we plot in figure 6 the distribution of the instantaneous radius of gyration as a function of the centroid positions for all H-bonds in the chain, where the points are taken every 100 MC steps along the PIMC simulation. The region with largest density of points in figure 6 coincides with the position of the peaks in both proton and deuteron cases (see figure 5). We again observe an asymmetric distribution centered in one of the sides of the H-bond consistent with the ( , ) distribution pattern of figure 5. The asymmetry observed in figure 6 is more pronounced in the proton case, indicating that protons jump more often than deuterons to the other side of the O-H-O bond. The mechanism to pass through the potential barrier is to increase the radius of gyration near ≈ 0 which means that the particle tunnels through the barrier. This is helped by a strong contraction of the distance, which diminishes concomitantly with 43708-8 the potential barrier, to a lower bound of min ≈ 2.3 Å near = 0 as shown in figure 5. Thus, we conclude that tunneling is assisted by the distance modulation. However, in this ordered phase at = 50 K, the proton spends more time in one of the sides of the O-H-O bond where the behavior is more classic (low value of ). On the other hand, in the deuteron case, the particle remains localized practically all the time, with a general classical behavior with low values of . In other words, the tunneling for the deuteron is very scarce. These results are consistent with the general assumption in the tunneling model: protons are capable of tunelling while deuterons are not [3]. However, there is an essential difference: protons tunnel being assisted by the strong correlation with the O-O distance, which is the behavior that originates the geometrical effect [9,10]. Therefore, the proton has a larger probability than the deuteron to spend more time tunneling through the barrier near the middle of the O-H-O bond, and this generates a strong attraction center that pulls the two oxygens together, much more efficiently than deuterons. This "tunneling -geometrical effect" interrelation gives rise to the final geometrical effect observed in KDP crystals, that is, the H-bond expansion with deuteration, which is crucial for the isotope effects in the FE-PE phase transitions [9,35].

Summary and conclusions
We have carried out PIMC simulations with simple models to account for the geometrical effects (GE) with deuteration in H-bonded ferroelectrics such as KDP crystals. Firstly, we have developed a general three-site (3S) model consisting in a back-to-back double Morse potential for the O-H interaction and a Morse potential which represents the interaction between the oxygens and the lattice. The model was fitted to reproduce general features for a large set of different H-bonded compounds. The computed probability distribution contours in the ( , ) configuration space, with the O-O distance and the proton/deuteron distance to the middle of the O-O bond, reveal a symmetric distribution around = 0 with two peaks on either side, for both proton and deuteron cases. The results show an increase with deuteration of and for the observed peaks, i.e., a GE, which is in agreement with that observed in H-bonded compounds with strong H-bonds. Moreover, if the oxygens are not allowed to relax during the simulation, the GE in the coordinate is much smaller, which means that there is a strong correlation between and that is important for the GE. During the PIMC simulations we have also plotted the instantaneous radius of gyration vs. the centroid position of the quantum polymer representing the proton/deuteron. The results show that the proton tunnels more frequently than the deuteron (that is, it spends more time with the center of mass near = 0 with large values of ), while the deuteron is 43708-9 more localized in both sides and far from the O-H-O bond center, with small values of (i.e., a more classsical behavior). These features yield a more effective contraction of the O-O bond in the proton case, explaining the GE observed.
Secondly, we have developed a more realistic 1D model, with the same local potential for the H-bonds as that used in the 3S model, but adding also a bilinear proton-proton interaction treated in mean-field. The parameters of the 1D model were fitted to ab initio results for KDP. The bilinear interaction parameter of the model was adjusted such that the classical nuclei version of the model has a second order FE-PE phase transition at = 350 K in agreement with ab initio molecular dynamics simulations for DKDP. In this paper, by means of PIMC simulations of the 1D model, we have studied the GE caused by deuteration in the ordered phase at = 50 K. The calculated probability distribution contours show only one peak in the ( , ) configuration space for both proton/deuteron cases. The distribution is more asymmetric in the proton case due to the anharmonicity of the potential and the quantum delocalization. The distribution pattern is in qualitative agreement with the experimental distribution determined by highresolution neutron diffraction studies [39]. The probability distribution contours show a peak which shifts substantially with deuteration. The changes in H-bond geomentry caused by the GE observed in the 1D model simulations are in good agreement with the corresponding experimental data. The distribution of the radius of gyration vs. the quantum path centroids shows that the protons tunnel through the potential barrier frequently while the deuterons are much more localized in one of the sides of the O-H-O bond and practically do not tunnel, in agreement with the well-known tunneling model [3], and also with recent neutron Compton scattering experiments [7,8]. We have shown that proton tunneling is assisted by a strong contraction of the O-O distance in the 1D model. Thus, there is a strong correlation between instantaneous tunneling and geometrical effects of the H-bond that is much more efficient in the proton case than in the deuterated system, which gives in average a strong GE for the whole simulation. This mechanism is expected to be at the root of the huge isotope effect observed in H-bonded ferroelectrics of the KDP type [9,10]. Метод iнтегралiв за траєкторiями у моделюваннi Монте-Карло (IТМК) для дуже простих моделей застосовано для з'ясування фiзичних механiзмiв, що лежать в основi iзотопiчного ефекту в сегнетоелектриках з водневими зв'язками. Зумовленi дейтеруванням геометричнi ефекти у водневих зв'язках було дослiджено за допомогою загальної тривузлової моделi, в якiй використовуються подвiйний потенцiал Морзе та потенцiал Морзе мiж киснями; параметри моделi вибрано так, щоб пояснити рiзноманiтнi загальнi властивостi низки сполук з водневими зв'язками. З розрахункiв у рамках цiєї моделi випливає виникнення геометричного ефекту (ефекту Уббелоде): видовження водневого зв'язка при дейтеруваннi, i це узгоджується з тим, що спостерiгається в сегнетоелектриках з короткими водневими зв'язками. Використовуючи для параметрiв потенцiалiв результати першопринципних розрахункiв, розвинено одновимiрну модель, в якiй бiлiнiйнi протон-протоннi взаємодiї розглядаються в наближеннi середнього поля. Ця модель використовується для дослiдження квантових ефектiв у ядрах, якi призводять до виникнення геометричного ефекту в кристалах KDP. Пiдхiд IТМК дає змогу виявити, що протони тунелюють бiльш ефективно вздовж одновимiрного ланцюжка, нiж дейтрони; це спричиняє появу сильного притягувального центра, який зменшує вiдстань мiж атомами киснiв. Цей механiзм, який ґрунтується на кореляцiї мiж тунелюванням i геометричними змiнами водневих зв'язкiв, призводить до виникнення сильного геометричного ефекту в ланцюжку у впорядкованiй фазi при низьких температурах, що добре узгоджується з експериментальними даними.