Effects of diagonal strains and H-bond geometry in antiferroelectric squaric acid crystals

The proton ordering model of the phase transition and physical properties of antiferroelectric crystals of squaric acid is modiﬁed by taking into account the inﬂuence of diagonal lattice strains and of the local geometry of hydrogen bonds, namely of the distance δ between the H-sites on a bond. Thermal expansion, the spontaneous strain ε 1 − ε 3 , and speciﬁc heat of squaric acid are well described by the proposed model. However, a con-sistent description of hydrostatic pressure inﬂuence on the transition temperature is possible only with further modiﬁcations of the model.


Introduction
The crystals of squaric acid, H 2 C 4 O 4 (3,4-dihydroxy-3-cyclobutene-1,2-dione) are an epitome of two-dimensional antiferroelectrics. The hydrogen bonded C 4 O 4 groups form planes parallel to ac and stacked along the b-axis. Below the transition at 373 K, a spontaneous polarization arises in these planes, with the neighbouring planes polarized in the opposite directions. The crystal symmetry changes from centrosymmetric tetragonal, I4/m, to monoclinic, P2 1 /m, and spontaneous symmetry-changing strains ε 1 − ε 3 (orthorhombic) and ε 5 (monoclinic), both of B g symmetry, arise [1][2][3]. The local symmetry of the H 2 C 4 O 4 groups is C 1h below and above the antiferroelectric phase transition [4].
Elastic and thermoelastic properties of squaric acid are remarkably anisotropic. Compressibility and thermal expansion [5] are much higher in a direction perpendicular to the planes of hydrogen bonds than within the planes. The symmetry-changing strains ε 1 − ε 3 and ε 5 are confined to the ac plane. The anomalous parts of the diagonal strains ε 1 , ε 3 and of ε 2 , caused by electrostriction, have different signs.
There is also experimental evidence for non-equivalence of hydrogen bonds going along two perpendicular directions (e.g. [1,6]). The difference between degrees of proton ordering on these bonds is about 2% at T N − 13 K and T N − 21 K [1]. The O-H and H-site distances are also found to be slightly different for the perpendicular bonds.
Hydrostatic pressure rapidly decreases the antiferroelectric transition temperature with the slope of about 11 K/kbar [4,7,8] at pressures below about 25 kbar. At higher pressures, this dependence deviates from linearity, and around 28 kbar, the transition temperature rapidly falls to zero: a quantum paraelectric state is induced [4]. At further increase of pressure at a constant temperature (at least for temperatures between 100 and 300 K), another phase boundary is detected [4], at which the local symmetry of the H 2 C 4 O 4 groups changes from C 1h to C 4h .
Theoretical description of the antiferroelectric transition in squaric acid is usually based on some versions of the proton ordering model, either two-dimensional, invoking four-particle correlations between protons within the planes [9][10][11][12], or one-dimensional, where either non-interacting [13] or coupled [14] perpendicular pseudospin chains are considered. The four-particle model can be reduced to the model of interacting one-dimensional chains by the proper choice of the model parameters [13]. The four-particle

The model
There are two formula units in the low-temperature phase unit cell of squaric acid. In our model, the unit cell consists of two C 4 O 4 groups and four hydrogen atoms ( f = 1, 2, 3, 4, see figure 1) attached to one of them (the A type group). All hydrogens around the B type groups are considered to belong to the A type groups, with which the B groups are hydrogen bonded. Note that the two C 4 O 4 groups of each unit cell belong to different neighboring layers. The center of each hydrogen bond lies exactly above the 33704-2 center of the hydrogen bond in the layer below it (as seen along the b axis). The bonds around each A type group are numbered counterclockwise. As it is usual in the proton ordering models, we consider interactions between protons leading to an ordering in their system. Motion of protons in double-well potentials is described by pseudospins, whose two eigenvalues σ = ±1 are assigned to two equilibrium positions of the proton. We take into account the presence of the diagonal components of the lattice strain tensor ε 1 , ε 2 , and ε 3 that are induced via thermal expansion or by application of external hydrostatic pressure.
The system Hamiltonian in the case of squaric acid includes ferroelectric intralayer long-range interactions H intra long , ensuring ferroelectric ordering within each separate layer, antiferroelectric interlayer H inter long responsible for alternation of polarizations in the stacked layers, the short-range configurational interactions between protons H short , and the so-called "seed" energy containing elastic and thermal expansion contributions associated with uniform lattice strains; c (0) i j are the corresponding "seed" elastic constants, whereas α (0) i are the "seed" thermal expansion coefficients. T 0 i determines the reference point of the thermal expansion of the crystal, which can be chosen arbitrarily. v is the unit cell volume, and N is the number of the unit cells in the crystal. Such a form of (2) later on yields standard expressions for the strains of a stressed thermally expanding solid [see equation (25)]. The short-range Hamiltonian H short describes the four-particle confirational correlations between protons sitting around each C 4 O 4 group. Similarly to how it is done for NH 4 H 2 PO 4 , the antiferroelectrics of the KH 2 PO4 type family, it is assumed that the energy of four lateral configurations ε a is the lowest of all, where two protons are in positions close to the adjacent oxygens of the C 4 O 4 group, whereas two other protons are closer to the neighboring C 4 O 4 groups (see figure 2). The next level is two diagonal configurations with the energy ε s , where the protons are close to the opposite oxygens of the C 4 O 4 group. Then, there are eight single-ionized configurations with three protons or only one close proton, having the energy ε 1 , and two double-ionized configurations (ε 0 ) with four or no protons at all close to the given C 4 O 4 group. It is believed that ε a < ε s ε 1 ε 0 . If two protons are in the most energetically favorable lateral configurations, the C 4 O 4 groups are isosceles trapezoids (point group C 1h ), although very close to squares (C 4h ). It is believed that this local distortion is caused by two single and two double alternating covalent bonds connecting the four oxygens to the carbons, and by formation of the double C-C bond within the C 4 O 4 skeleton, as shown in figures 1 and 2. Double bonds are shorter than single ones between analogous atoms. Since the origin of the skeleton distortion is the chemical bonding with the local proton configuration, rather than macroscopic uniform lattice strains, all four lateral configurations still have the same energy of short-range interactions, no matter what their orientation relatively to the crystallographic axes is (see table 1). The same holds for the diagonal (point group C 2h ), single-ionized (point group C 1 ), and double-ionized groups (point group D 2h ). It means that no splitting of short-range energy levels by the macroscopic spontaneous strain takes place, in contrast to what was assumed in earlier theories for squaric acid [9] or for KH 2 PO 4 type crystals [23,24].

33704-3
To proceed from the representation of proton configuration energies to the pseudospin representation, we use a standard procedure, originally developed for the KH 2 PO 4 type crystals [25,26], where the Hamiltonian of the short-range correlations between protons, surrounding each A type C 4 O 4 group is written as follows: whereN i (yq) is the operator of the four-particle configuration i; s f = ±1 is the sign of the eigenvalue of the σ yq f operator in this particular configuration; E i is the energy of the configuration. It is assumed that s f = +1 if a proton at the f th bond is localized at the H-site proximal to the particular A type C 4 O 4 group, and s f = −1 if the proton is localized at the other (distal) H-site of the same bond. Using equation (3), we arrive at the following expression for the Hamiltonian The short-range interaction constants 33704-4 Table 1. Proton configurations and their energies; ε a < ε s ε 1 ε 0 .
are linear functions of the Slater-Takagi type energy parameters Note that the model of non-interacting perpendicular one-dimensional chains is obtained from equation (18) at Φ = V = 0, i.e. at ε = w 1 = 2w, where the following order of the configuration energies should be assumed [13] ε a < ε 1 < ε s = ε 0 . It can be shown, as it has been done for the KH 2 PO 4 ferroelectrics, that the contributions of the correlations from the A and B type groups to the total thermodynamic potential are equal. The Hamiltonian of the short-range interactions in this case can be written as follows: where the expression for H A q y is given by equation (4). The long-range interactions in the system Hamiltonian (1) are the dipole-dipole interactions, in analogy to the case of KH 2 PO 4 [25], renormalized by the proton-lattice coupling. The MFA is successfully used for these interactions in the calculations for hydrogen-bonded ferroelectrics, if the fluctuations in the vicinity of the phase transition temperature neglected in this approximation are not the subject of special interest.
Within the MFA, we obtain the following expressions for the long-range intralayer 33704-5 and for interlayer interactions. Here N y is the total number of the layers. The internal mean fields are as follows: The following symmetry of the pseudospin mean values is assumed for the antiferroelectrically ordered two-sublattice model in the absence of external electric field Here, k 2 = (0, b 2 , 0); b 2 is the basic vector of the reciprocal lattice; the factor exp[ik 2 R y ] = ±1 denotes two sublattices of an antiferroelectric, R y is the position vector of the y-th layer, and The last relation reflects the slight non-equivalence, mentioned in introduction, of hydrogen bonds, linking C 4 O 4 groups along the a and c axes. Even though each particular interaction parameter J intra f f (qq ) and J inter f f (y y ; qq ) is obviously changed by the strains ε 1 and ε 3 , the symmetry of the long-range interaction matrices Fourier transforms over the bond indices f and f , as can be easily checked, remains unchanged in the presence of the orthorhombic strain ε 1 − ε 3 : both for J intra (0) and J inter f f (k 2 ). The symmetry (12) is obvious for strictly square-shaped C 4 O 4 groups (point group C 4h ). This is a statistically average symmetry of the paraelectric phase, where the hydrogens are placed, also statistically, in the middle of the hydrogen bonds.
Taking into account equations (10)-(12), we can write the Hamiltonians of the long-range interactions as follows: where We also took into account the fact that N y N q A = N.
For the sake of simplicity, we shall hereafter neglect the weak non-equivalence of the perpendicular chains of hydrogen bonds and use a single order parameter instead of equation (11).
The four-particle cluster approximation will be used for the short-range interactions, described by the Hamiltonian (4). With the long-range interactions taken into account in the mean field approximation, 33704-6 the thermodynamic potential of the system should be written as follows: Here The four-particle cluster Hamiltonian is where The fields ∆ yq f are the effective cluster fields that describe short-range interactions of the spin σ yq f with the particles from outside the cluster q. They are determined from the self-consistency condition stating that pseudospin mean values calculated with the four-particle (18) and with the one-particle Hamiltonians should coincide. We get Taking into account equations (13), (15), (18), (19), the thermodynamic potential per one unit cell is obtained in the following form where D = a + cosh 2z + 4b cosh z + 1, a = exp(−βε), b = exp(−βw).
In the earlier theories [18,27], the short-range Slater-Takagi energies in the KH 2 PO 4 family crystals were considered as quadratic functions of the distance δ. For the squaric acid, we shall employ the same scheme. Using the term of the relative deviation of δ from its value δ 0 at ambient pressure (we shall call it a displacement µ ) we assume that Here, the quadratic in µ terms, omitted in [18], are included into consideration. For the parameter of long-range (dipole-dipole) interactions ν, we take into account both the dependence of the dipole moments on δ and the changes in the interaction parameter due to the overall crystal deformation [18] and associated with changes in the equilibrium distances between protons (dipoles)

33704-7
It should be underlined that none of the earlier theories [18,27] described the thermal expansion of the crystals; therefore, the deformational effects there were only those caused by external pressures. On the contrary, in the present model the strains ε i are induced both by temperature changes and by external pressures. In the mean field approximation, the expansion (23) gives rise to the terms of electrostriction type in the Hamiltonian, linear in the strains and quadratic in the sublattice polarization (order parameter η).
In [18,27], the distance δ was treated as a pressure dependent and temperature independent model parameter, with the linear pressure dependence chosen either from the available experimental data or by fitting the theory to experiment for the transition temperatures. In the present work, we shall use a similar approach and assume δ to vary according to its experimentally observed above the transition linear temperature [1] and external hydrostatic pressure p [17] dependences where T N0 is the transition temperature at ambient pressure. It means that the anomalous temperature behavior of δ below the transition point and its jump at T N are neglected. Minimization of the thermodynamic potential (20) with respect to the order parameter η and strains ε i ∂g ∂η = 0, ∂g ∂ε i = 0 yields the following equations From the above it follows that in equilibrium where s (0) ki is the matrix of "seed" elastic compliances, inverse to c (0) i j . One can see that in the paraelectric phase (η = 0), the microscopic contributions to the strains vanish, whereas in the ordered phase they are proportional to η 2 , indicating the electrostriction type contributions governed by the parameters ψ i .
Finally, the molar entropy of the proton subsystem is as follows:

Calculations
In the calculations, the thermodynamic potential is minimized numerically with respect to the order parameter η. At the same time, the strains ε i are determined from equation (25).
The quantities that should be described include: the temperature curves at ambient pressure of There are four different "seed" elastic constants c (0) 11 , c (0) 22 , c (0) 12 , and c (0) 13 . We take c (0) 11 to be equal to the experimental value of c 11 [28] above the transition point. Experimental elastic constant c 22 was found to slightly decrease with an increasing temperature [28,29]. The "seed" c (0) 22 is chosen accordingly, coinciding with the data of [28]. Finally, c (0) 12 and c (0) 13 , for which no convincing experimental data are available, were chosen to provide a correct fit to the experimental [30] pressure dependence of the lattice constants a and b at 292 K.
The parameters of the short-range correlations ε 0 and w 0 govern the temperature behavior of the order parameter η (in particular, the magnitude of its jump at the transition ∆η c and steepness of its rise to saturation with lowering temperature) and the value of the transition temperature at ambient pressure T N0 . The latter is also extremely sensitive to the value of the long-range interactions parameter ν 0 . Hence, the set of ε 0 , w 0 , and ν 0 is chosen to yield T N0 = 373.5 K, ∆η c ≈ 0.57, as well as a correct temperature curve of η between the transition and saturation. Contributions of the double-ionized configurations are neglected by putting w 1 → ∞.
The "seed" thermal expansion coefficients α (0) i as well as the parameters ψ i are determined by fitting the theoretical temperature dependences of the lattice strains to experimental data [5,31]. In fact, α (0) i should be simply equal to the corresponding paraelectric experimental values, as is seen directly in equation (25). The parameters ψ i , on the other hand, are unambiguously determined by fitting the calculated anomalous parts of the strains to the experiment below the transition temperature, using the same equation (25). As ψ i are relevant for the ordered phase only, they do not have to adhere to the symmetry of the paraelectric phase; hence, we can take ψ 1 ψ 3 , as is indeed required by the just described fitting.
As we have already mentioned, the temperatures T 0 i determine the reference point of the thermal expansion of the crystal, which can be set arbitrarily. Thus, T 0 i are not fitting parameters of the model and, therefore, can also be chosen arbitrarily. In our calculations we chose them to yield zero values of the lattice strains ε i just above the transition temperature at ambient pressure. In fact, T 0 1 = T 0 2 = T N0 = 373.5 K, as seen from equation (25).
As already described, we take δ to vary according to its experimentally observed linear temperature and external hydrostatic pressure (24). The coefficients δ T and δ p are deduced from the data of [1] and [17].
The unit cell volume is v = 2 · 10 −28 m 3 just above the transition point at ambient pressure, as determined from the data of [5]. The final values of the model parameters, used in our calculations, are given in table 2.
In figure 3 we show the calculated temperature dependence of the order parameter η and the spontaneous strain ε 1 − ε 3 at ambient pressure. Experimental points for η were obtained from the 13 C NMR measurements. A clear first order phase transition is observed, with the jump of the order parameter ∆η c 0.57. The spontaneous strain ε 1 − ε 3 is negative below the transition and, as it follows from equation (25), it is proportional to the square of the order parameter η 2 .  The temperature dependences of the diagonal lattice strains ε i and the corresponding thermal expansion coefficients are shown in figure 4. The coefficients were calculated by numerical differentiation of the strains with respect to temperature.  A clear anisotropy of the thermoelastic properties of squaric acid within the ac plane and in the perpendicular direction is seen. At temperature lowering, the strain ε 2 has a downward jump at the transition and a negative anomalous part in the ordered phase. The strains ε 1 and ε 3 , on the other hand, have upward jumps and positive anomalous parts. As is shown above [see equation (25)], the anomalous contributions to the macroscopic strains are strictly proportional to the square of the order parameter η 2 . The thermal expansion coefficients α 1 = α 3 in the paraelectric phase are by one order of magnitude smaller than α 2 ; their anomalies at the transition point are of different signs.
The specific heat of the proton subsystem is calculated by numerical differentiation of the entropy (26) with respect to temperature where M = 114.06 g/mol is the molar mass of squaric acid. The corresponding temperature curve is given in figure 5. The experimental points for the anomalous part of the specific heat are obtained by subtracting the regular part, best described as a slightly non-linear curve C reg = −0.48803 + 0.00738T − 7.3 · 10 −6 T 2

T(K)
∆C P (J/g K) Figure 5. Temperature dependence of the specific heat of squaric acid. Line: the theory; symbols are experimental points derived from the data of [33] as described in text.
(J/g K), from the total specific heat as it was measured in [33]. One can see that a satisfactory agreement with experiment is obtained, even though the specific heat was not directly involved in the above described fitting procedure. The calculated hydrostatic pressure dependences of the paraelectric lattice constants are shown in figure 6. The lattice constants were determined as a = a 0 [1 + ε (1) , where a 0 = 6.137 Å, b 0 = 5.327 Å are the values of the lattice constants just above the transition point at ambient pressure [5]. A good agreement with experiment is obtained.
In figure 7 we plot the hydrostatic pressure dependence of the phase transition temperature in squaric acid. As expected, the calculated transition temperature decreases with pressure (the dashed line). Quantitatively, however, completely non-satisfactory results are obtained. With the pressure variation of µ as observed experimentally [17] and the parameters ψ i determined by fitting to the lattice strains below transition at ambient pressure, the rate, with which the calculated transition temperature decreases with hydrostatic pressure, ∂T c /∂p = −19.5 K/kbar, is nearly twice as large as the experimental one. The observed disagreement means, foremost, that equation (23) yields a too fast decrease of the long-range interaction parameter ν, to which the theoretical values of the transition temperature are most sensitive. The pressure variation of the Slater-Takagi energies is less important here.
Below we discuss a possible origin of the model inconsistency and ways to solve this problem. To this end, let us look closely at the obtained pressure dependence of ν.
It is expected that the term i ψ i ε i in presence of high hydrostatic pressure would be positive, thereby leading to an increase of the long-range interaction parameter ν due to the reduction of the average interparticle distances in the compressed crystal. However, when the values of the parameters ψ i are chosen to fit to the experimental data [5] for the anomalous spontaneous temperature behavior of the strains ε i below the transition at ambient pressure, the sum i ψ i ε i in presence of high pressure becomes negative, not slowing, as expected, but enhancing the decrease of ν caused by the decrease of the H-site distance δ. Of course, one cannot exclude a possibility that further measurements may yield the values of the strains ε 1 (T), ε 3 (T) different from those obtained in [5] and unconfirmed yet by other measurements. The reasoning below, however, is based on the assumption that the data of [5] are correct.
An obvious workaround but rather clumsy way to obtain the necessary pressure dependence of ν is to assume that there are some other high-pressure factors, not included into (22) and (23), and to include them empirically via the terms k p p, namely and We can speculate, for instance, that external pressure causes a redistribution of electron density, thereby changing the effective charges of the ions and, as a result, chanding the interactions between them. Introduction of two extra fitting parameters k p1 and k p2 indeed allows us to describe the pressure variation of the transition temperature (see figure 6, the solid line). At k p1 = k p2 = 0.0151 kbar −1 we obtain ∂T c /∂p = −10.7 K/kbar, in total agreement with experiment. Other combinations of k p1 and k p2 values can be found, also yielding the desired fit for the T c vs p dependence. Accepting this, the already obtained good agreement with experiment for the system behavior at ambient pressure is not affected.

Conclusions
We present a modification to the proton ordering model, aimed at describing the effects associated with diagonal lattice strains in H-bonded antiferroelectric crystals of squaric acid. These effects include thermal expansion of the crystals, the appearance of spontaneous strain ε 1 − ε 3 below the phase transition, and the shift of the transition temperature with hydrostatic pressure. Here, both the macroscopic lattice strains and the changes in the local geometry of hydrogen bonds are found to be essential. As usually, the quadratic dependence of the parameters of short-range and long-range interactions between protons on the H-site distance δ is assumed.