Coexisting conventional and inverse mechanocaloric effects in ferroelectrics

The mechanocaloric effect is the temperature change of a material upon application or removal of an external stress. Beyond its fundamental interest, this caloric response represents a promising and ecofriendly alternative to current cooling technologies. To obtain large mechanocaloric effects, we need materials whose elastic properties (e.g., strain, elastic compliance) are strongly temperature dependent. This is the case of ferroelectric perovskite oxides, where the development of the spontaneous electric polarization is accompanied by significant strains and lattice softening. Thus, in this work we study the mechanocaloric properties of model ferroelectric PbTiO$_{3}$, by means of predictive atomistic ("second-principles") simulations and a perturbative formalism here introduced. Our calculations reveal relatively large effects (up to $-$4 K for relatively small applied compressions of $-$0.1 GPa) and several striking features. In particular, we find that the mechanocaloric response is highly anisotropic in the ferroelectric phase, as it can be either conventional (temperature increases upon compression) or inverse (temperature decreases) depending on the direction of the applied stress. We discuss and explain these surprising results, which compare well with existing experimental information. Our analysis suggests that the coexistence of conventional and inverse mechanocaloric responses is probably common among ferroelectrics and materials displaying a negative thermal expansion.


I. INTRODUCTION
A material subject to an external stress changes its temperature and/or entropy by virtue of the socalled mechanocaloric effect (also denoted elastocaloric or barocaloric, respectively, when the applied stress is uniaxial or hydrostatic) [1][2][3][4]. Mechanocaloric effects in solid-state systems can be large, reaching values from 5 K to 40 K in shape-memory alloys [5], super-ionic conductors [6], or plastic crystals [7]. This remarkable performance makes mechanocalorics a viable alternative for applications in solid-state cooling, where it is pressing to replace current technologies (e.g., compressed gas) that are polluting, noisy, and difficult to downscale.
Mechanocaloric effects are large in ferroelastic materials, which are strongly responsive to mechanical perturbations and often display stress-driven phase transformations, particularly close to the ferroelastic transition point [8]. Interestingly, ferroelectrics (e.g., perovskite oxides like BaTiO 3 and PbTiO 3 ) are known to be very sensitive to mechanical perturbations too [9], and can in principle be expected to present a considerable mechanocaloric response. Nevertheless, while ferroelectrics have been thoroughly studied in regards to electrocaloric properties [10], their mechanocaloric performance has been scarcely investigated so far [11][12][13][14].
In this article we present a theoretical study of the intrinsic mechanocaloric response of prototype ferroelectric perovskite PbTiO 3 . Our work is based on predictive atomistic ("second-principles") simulations [15][16][17] of the temperature-dependent elastic properties of PbTiO 3 , combined with a perturbative formalism that we introduce here and should be useful in the broader context of mechanocaloric investigations, including experimental ones. Among other remarkable features, we find that the ferroelectric state of PbTiO 3 presents a acute mechanocaloric anisotropy, to the extent that the sign of the adiabatic temperature change depends on the direction along which the stress is applied. We argue that this feature is probably common among ferroelectric compounds, as well as in materials displaying negative thermal expansion.

II. METHODS
In the following we describe our perturbative formalism to compute and interpret the mechanocaloric response, as well as the second-principles simulation methods that allow us to study the specific case of PbTiO 3 .

A. Perturbative formalism
We adapt to the mechanocaloric case a perturbative formalism recently introduced by some of us in the context of electrocaloric effects [18]. Our starting point is the usual thermodynamic expression for the adiabatic temperature change caused by an applied stress [4], where the integral runs from zero to the final stress σ a , η a denotes the strain conjugate to the applied stress, T the temperature and C σ the specific heat at constant stress. Here we introduce the thermal expansion coefficient α a , indicating that it depends on temperature and arXiv:2109.00988v1 [cond-mat.mtrl-sci] 2 Sep 2021 stress. The subscript a labels stress and strain components in Voigt notation [19] and we assume a Cartesian reference. Note that we write ∆T (σ a ) to indicate explicitly that the temperature change can be anisotropic, i.e., it may depend on the direction of the applied stress. Finally, we adopt the usual sign convention that a negative external stress corresponds to a compression of the material. As usually done in the literature [2], in the following we work with an approximate version of the temperature change, where the "0" superscripts indicate values calculated at zero applied stress and we assume that C σ and α a are independent from the mechanocaloric temperature change.
(This is a reasonable approximation since, in most situations of interest, |∆T (σ a )| T (0) ≈ 300 K.) Besides being simpler and easier to treat, this formula allows us to focus on how the thermal expansion α a controls the effect. Note that this is the quantity that can be expected to be most strongly T -and stress-dependent; and the only one whose sign is undetermined and free to change.
Now we write the strain in terms of the stress as a Taylor series, where η (0) a stands for the spontaneous strain at zero stress; S (0) ab is the compliance tensor at zero stress, with Taking the temperature derivative of Eq. (3), we obtain where α (0) a is the thermal expansion vector of the unperturbed material and we have which we can view as stress-induced thermal expansion coefficients. Using these simple expressions, now we evaluate Eq. (2) in several cases.

Uniaxial stress
Imagine we apply stress along a specific Cartesian axis a, so that σ b = δ ab σ a where δ ab is the Kronecker delta.
We can then substitute Eq. (5) into Eq. (2) and solve the integral to obtain where we introduce ∆T (n) (σ a ) to denote the contribution of n-th order to the total adiabatic temperature change.

Biaxial stress
Suppose now that we apply stress along two Cartesian directions simultaneously. To fix ideas, let us imagine we work with σ 1 and σ 2 , both varying from zero to a final value σ, the generalization of our arguments being trivial.
When we work with two stress components, the onedimensional integrals presented above are not directly applicable, and we need to take some extra intermediate steps. Let us discuss two equivalent ways to compute the mechanocaloric temperature change in this case.
Method I. We can rotate our Cartesian coordinates so that one of the axes coincides with the direction of the applied stress, and thus recover a one-dimensional problem. In this particular example, we can use the unitary transformation where σ = 2 −1/2 (σ 1 + σ 2 ) defines the stress we intend to apply. (The inverse transformation gives σ 1 = σ 2 = 2 −1/2 σ , assuming σ ⊥ = 0.) The strain η , conjugate to σ , is similarly given by η = 2 −1/2 (η 1 + η 2 ), and we also have α = 2 −1/2 (α 1 + α 2 ) for the corresponding thermal expansion coefficient. We have thus defined an equivalent one-dimensional problem where the applied stress σ varies from zero to √ 2σ. Now, to compute the temperature change, we need to solve the integral where Here we have used the relation α 21 , which relies on the fact that the compliance matrix S is symmetric. Also, for conciseness, we have truncated the series to the linear order in the applied stress. We can finally evaluate the integral to obtain 2 )σ Method II. Alternatively, we can imagine that the stress is applied in two steps. (The result of our integral for the temperature change does not depend on how we apply the stress, a property that relies on the fact that the underlying free energy is (assumed to be) an exact differential.) In the first step, we have σ 1 varying from zero to σ, with σ 2 = 0 throughout. The corresponding temperature change is given by the original one-dimensional integral In the second step, we have a constant σ 1 = σ, while σ 2 varies from zero to σ. The corresponding onedimensional integral is where one must note that Hence, inserting Eq. (14) into Eq. (13), we get (15) Finally, we calculate the total temperature change by adding the individual variations, to obtain which coincides exactly with the outcome of Method I described above (Eq. (11)). It is worth noting that this two-step approach is, in essence, directly applicable to any multicaloric effect, where, instead of considering two different components of the same field (stress in our case), one applies two fields of different nature (e.g., mechanical and electric).

Hydrostatic pressure
We finish with the case of an isotropic stress (hydrostatic pressure) where all σ 1 = σ 2 = σ 3 vary from zero to σ. The two methods described above for the biaxial case are readily applicable, the final result for the adiabatic temperature change being

B. Simulations for PbTiO3
In this work we evaluate the quantities controlling the mechanocaloric temperature change (specific heat, spontaneous thermal expansion, etc.) by running Monte Carlo simulations of an atomistic second-principles model of PbTiO 3 . This model, described in detail in Ref. 15, has successfully been used in several theoretical investigations of PbTiO 3 and related compounds [20][21][22], including a recent study of the electrocaloric response [18]. Its only noteworthy deficiency is that it predicts a ferroelectric transition temperature that is too low as compared with the experimental one (510 K vs 760 K). However, as we show below, our results for the key quantities of interest here (most notably, the thermal expansion and specific heat) compare very well with experimental values upon a simple temperature shift that makes the theoretical and experimental Curie points coincide; hence, this quantitative discrepancy is not critical for the present purposes.
In the simulations we typically use a periodically-repeated supercell containing 10×10×10 instances of the 5-atom perovskite unit cell; we run 10,000 Monte Carlo sweeps for thermalization and 100,000 more to compute averages. Near the ferroelectric transition temperature we find it necessary employ a larger 12×12×12 supercell, to reduce finite-size effects; in that case we run 10,000 thermalization sweeps followed by 75,000 production sweeps. For the sake of simplicity, in the results presented below the ferroelectric polarization is always chosen to lie along the z axis.
To obtain thermodynamic properties (e.g., equilibrium strain, specific heat, elastic compliance tensor) from our Monte Carlo simulations, we proceed in the usual manner, and rely on well-known linear-response formulas (see, e.g., the description in Ref. [18]). To illustrate the linearresponse approach, let us consider the case of the elastic compliance, which we have not been able to find in previous literature. The compliance component S ab is calculated as where η a is the average equilibrium value of the strain at given conditions of temperature and applied stress. More specifically, η a is where Z is the partition function Here i runs over microscopic states of the simulated sys- is the equilibrium volume of the simulation supercell at the considered conditions of temperature and stress, and U i and η i are, respectively, the energy and strain of state i. Then, by taking the derivative of η a with respect to the applied stress, we can readily obtain which is our linear-response formula for the compliance. As usually done in simulation studies of piezoelectric properties that involve similar quantities [23], here we neglect the stress derivative of the volume; we explicitly checked this is an accurate approximation.
To finish this section, let us note that in this study we define strains by taking as a reference the lattice constant of the cubic phase of PbTiO 3 , as obtained from a symmetry-constrained first-principles simulation at 0 K. This reference value is l 0 = 3.88Å. Hence, for the computed normal strains η a (where a = 1, 2, 3 labels the pseudo-cubic perovskite axes), the corresponding lattice constants are given by Note that we can use this simple expression because shear strains (η a with a = 4, 5, 6) are zero at all temperatures in PbTiO 3 . Finally, in the following we choose the Cartesian axes to coincide with the pseudo-cubic directions of the perovskite lattice. Also, we assume that the tetragonal ferroelectric phase of PbTiO 3 is characterized by a positive polarization along the third Cartesian axis.

III. RESULTS AND DISCUSSION
We now present our results, emphasizing and discussing the most interesting aspects revealed by our calculations. We also compare with available experimental information.

A. Basic Monte Carlo results
In Fig. 1 we present the temperature dependence of the polarization as obtained from our Monte Carlo simulations (panel (a)) as well as our results for the specific heat (panel (b)). We obtain a weakly first-order ferroelectric phase transition at T C = 510 K, where one polarization component becomes different from zero in a discontinuous fashion. This transition, from the hightemperature cubic (P m3m) phase to a tetragonal polar (P 4mm) state, agrees with the well-known experimental behavior of PbTiO 3 [9] (albeit the underestimated transition temperature mentioned above) and with previous simulations using the same second-principles potential [15]. The specific heat presents an anomaly (a maximum) at the transition temperature, also as expected.
In Fig. 2 we present how the strain-related quantities evolve across this phase transition. Panel (a) shows the temperature evolution of the normal strains (η 1 , η 2 and η 3 ) as well as their average (η ave , which accounts -to first order -for variations in volume). (We do not show the shear strains, which are zero at all temperatures.) The symmetry breaking associated to the ferroelectric transition is clearly visible in the figure. Indeed, this transformation has an improper ferroelastic character, whereby we observe strain changes following the onset of the spontaneous polarization.
More precisely, Fig. 2(a) shows two differentiated behaviors. In the paraelectric phase, for T > T C , the three strain components are equal (η 1 = η 2 = η 3 ) and increase gradually with temperature. Hence, we have an isotropic positive thermal expansion. In contrast, when we heat up the material from low temperatures in the ferroelectric phase (T < T C ), the strain component parallel to the polarization (η 3 ) decreases with increasing temperature, while the components perpendicular to it (η 1 = η 2 ) grow. This reflects the presence of a tetragonal distortion (usually quantified by the c/a aspect ratio of the tetragonal unit cell), which can be shown to be proportional to the square of the spontaneous polarization (see, e.g., the discussion in Ref. 24). Interestingly, the combination of these normal strains results in an average deformation that decreases with increasing temperature, yielding a negative thermal expansion for T < T C . This effect has been discussed in the literature, e.g., in Refs. 25 and 26. S 12 = S 13 = S 23 ) is obvious here as well. We also find that, in the ferroelectric phase, the direction of the polarization is elastically softer than the perpendicular plane (S 33 > S 11 = S 22 ). Further, all compliance components reach maximum values around the transition point, presenting a simple monotonic variation at all other temperatures.
Finally, in panels (d) to (f) of Fig 2 we show the temperature derivatives of the previous quantities, which we obtain numerically from the data in the first three panels. These are the α (0) a and α (0) ab tensor components introduced above, the key ingredients to compute mechanocaloric temperature changes within our perturbative formalism.

B. Mechanocaloric temperature change
The results in Figs. 1(b) and 2 allow us to compute the adiabatic mechanocaloric temperature change using the formalism introduced above. Representative results are given in Fig. 3, all corresponding to the application of a compressive stress of −0.1 GPa. We consider a uniaxial stress along the polarization direction in panel (a), a biaxial stress in the plane perpendicular to P in panel (b), and a hydrostatic pressure in panel (c). In all cases we show the total ∆T as well as the individual contributions with a linear (∆T (1) ) and quadratic (∆T (2) ) dependence on the stress. Two aspects are common to all the results in Fig. 3. First, the largest temperature changes are obtained in the vicinity of the phase transition. This is consistent with the results in Fig. 2, which shows that all the relevant quantities (T -derivatives) present their largest absolute values around T C . Second, for the considered −0.1 GPa, the linear effect dominates the caloric response, even in the vicinity of the phase transition. As a matter of fact, the linear approximation to ∆T is a good one in a wide range of applied stress; indeed, as shown in Fig. 4, we need to reach compressions of about −1.5 GPa for the quadratic effect to dominate at temperatures close to T C (see results for 523 K in panel (a)), and even higher when we move away from the transition temperature (−10.7 GPa at T = 413 K, as shown in panel (b)).
We also find that, in all cases, the obtained ∆T is positive and relatively small above T C , as consistent with the computed small isotropic thermal expansion at high temperatures (see Fig. 2(a)). A marked anisotropy of the caloric responses appears when we move below T C , reflecting the lower symmetry of the tetragonal ferroelectric state.
Indeed, most importantly, we find that the ferroelectric state may present a conventional (∆T > 0) or inverse (∆T < 0) effect depending on how the compressive stress is applied. When we compress in the plane perpendicular to the polarization (Fig. 3(b)), we obtain ∆T > 0 reflecting the positive thermal expansion of the a and b lattice constants. In contrast, a compression along the polar axis ( Fig. 3(a)) yields a negative temperature change, as a result of the reduction of the c lattice constant upon heating.
Hence we obtain a remarkable result: because of the development of the ferroelectric polarization and its impact in the lattice strains, PbTiO 3 turns out to be an extremely anisotropic mechanocaloric, to the point that the sign of the temperature change depends on the direction of the applied stress. As far as we know, this is the first example of mechanocaloric material presenting such a sign anisotropy. (Examples of materials that can switch between conventional and inverse behaviors are known, though [27][28][29].) Note also that the relationship just described between polarization and strain is not exclusive of PbTiO 3 , but typical of ferroelectric perovskites (see, e.g., the behavior of the tetragonal phase of BaTiO 3 [30,31]) and other ferroelectric families (see, e.g., the case of LiTaO 3 [32]). Hence, the coexistence of conventional and inverse mechanocaloric responses may be a common feature among ferroelectrics. In addition, materials with a negative thermal expansion (where, often, some lattice constants decrease and some increase upon heating [26,33,34]) are also good candidates to present such a behavior.
It is also worth to mention that our computed temper- , we get temperature changes exceeding 3 K close to the phase transition (510 K), and still notable as we move into the ferroelectric state (e.g., above 1 K for all temperatures between 450 K and 510 K). Also remarkably, the gap between the positive and negative effects is as large as 7 K at T C . In contrast, the response to a hydrostatic pressure (Fig. 3(c)) is relatively small within the ferroelectric phase, as it suffers from the partial cancellation of the conventional and inverse effects just discussed. We have a maximum barocaloric ∆T ≈ −1 K at T C , the inverse effect being dominant. This is a direct consequence of the negative volumetric thermal expansion obtained for T < T C and shown in Fig. 2(a). Naturally, as we cross T C and the isotropic thermal expansion changes from negative to positive, so does the sign of the barocaloric response. The same applies to the elastocaloric response to an uniaxial stress along the polar axis.

C. Comparison with previous works
We can compare our theoretical predictions with experiment in different ways.
First, based on experimental results for the temperature dependence of the lattice constants of PbTiO 3 (we used those reported in Ref. [35]), we can numerically obtain the spontaneous thermal expansion vector α (0) . Then, we shift down the experimental temperatures to make the transition point coincide with the calculated T C = 510 K. With this data, combined with our calculated results for C σ , we can trivially estimate ∆T (1) , i.e., the part of the mechanocaloric response that is linear in the applied stress and which, according to our calculations, dominates the effect for moderate compressions. Figure 5 shows the results for ∆T (1) thus obtained (symbols labeled "exp. (I)"), together with our theoretical predictions. We find that the agreement between theory and the experiment-based estimate is essentially perfect, in- dicating that our simulations capture very well (save the error in T C ) the thermal evolution of PbTiO 3 's lattice. In particular, the change of sign of the barocaloric effect across T C is readily recovered when using the experimental data for α (0) , which displays the same reversal across the Curie point (from negative to positive thermal expansion) as computed.
Second, we deduce ∆T (1) exclusively from experimental information, by combining the structural results of Ref. 35 with the calorimetric data of Ref. 11 (shifted too, so that the experimental and theoretical T C 's coincide). We thus obtain the symbols labeled "exp. (II)" in Fig. 5, which are in excellent agreement with our theoretical values and our first experimental estimate ("exp. (I)"). This indicates that the (temperature-shifted) specific heat is in excellent agreement with our predictions.
Indeed, for example, Refs. 11 and 36 report values between 2.6 MJ m −3 K −1 and 2.9 MJ m −3 K −1 for C σ near room temperature, while our computed value at 300 K is about 3.4 MJ m −3 K −1 . Further, after the mentioned temperature shift, we obtain C σ = 3.4 MJ m −3 K −1 at 300 K from the data in Ref. 11, in perfect agreement with our calculation. This suggests that our computational approach predicts PbTiO 3 's thermal properties with remarkable accuracy.
Finally, the barocaloric effect in PbTiO 3 was reported in Ref. 11, which in principle allows us to make a direct comparison. Below T C , this article describes an inverse behavior that is similar to the one we predict, with quantitative results (e.g., ∆T ≈ −2 K at the transition temperature for a pressure of −0.03 GPa) in reasonable agreement with our computed values. However, Ref. 11 also reports that the negative temperature change persists above T C for all the measured pressures (e.g., ∆T ≈ −0.1 K is measured at about 25 K degrees above T C for −0.03 GPa), which would suggest that PbTiO 3 presents a negative thermal expansion in its paraelectric phase. This is in direct contradiction with abundant ex-perimental [37][38][39] and theoretical [15,40,41] studies of the structural evolution of this compound. Hence, we are not sure about the status of the barocaloric results of Ref. 11 or to what extent we should expect agreement with our predictions.
Additionally, we should stress we have not found any experimental investigation addressing the (strong) anisotropy of the mechanocaloric response in the polar phase of PbTiO 3 . This question, which constitutes one of our most interesting predictions, remains to be explicitly verified experimentally.
Finally, let us note that our results are consistent with other theoretical studies of PbTiO 3 . For instance, using a first-principles-based effective Hamiltonian, Lisenkov et al. [12] predict ∆T = +6 K under the application of a tensile stress of +0.2 GPa along the polar direction and near the transition temperature; this is in acceptable agreement with our computed ∆T ≈ −4 K upon an uniaxial compression of −0.1 GPa close to T C . Our results are also consistent the effective-Hamiltonian study of Barr et al. [13], who report a temperature change of about −5 K upon release of a tensile load of +0.2 GPa close to T C . Similarly, a phenomenological study of PbTiO 3 and related ferroelectrics [14] reports temperature increases of up to +2 K for a tensile uniaxial pressure of +0.1 GPa along the polarization direction and close to T C . Hence, all these literature results are quantitatively and qualitatively consistent with the dominant linear effect revealed in the present work, whereby the sign of ∆T ≈ ∆T (1) depends on the nature (tensile or compressive) of the applied stress.

IV. SUMMARY AND CONCLUSIONS
In this work we have used predictive atomistic simulation ("second-principles") methods to investigate the intrinsic mechanocaloric response of prototype ferroelec-tric PbTiO 3 . Notably, we find that the effects can be quite large in the vicinity of the ferroelectric Curie point (up to −4 K for compressions of only −0.1 GPa), even though no phase transition is induced in the material (hence, our obtained values do not have any latent-heat contribution).
Remarkably, we reveal that the mechanocaloric response is strongly anisotropic in the ferroelectric phase of the compound. More precisely, we find that the effect is conventional (temperature increases under compression) if a stress is applied in the plane perpendicular to the spontaneous polarization, and inverse (temperature decreases) if we compress along the polarization direction. As far as we know, such a coexistence of conventional and inverse responses had never been observed or predicted before in any material, and remains to be explicitly confirmed experimentally. Nevertheless, we find that our results are compatible with available experimental information on the elastic properties of PbTiO 3 , which suggests that the predicted coexistence is real.
Our theoretical calculations rely on a perturbative formalism that we introduce here and which should be useful in the broader context of mechanocaloric studies, including experimental ones. (We illustrate this explicitly when checking our predictions against experimental information.) We should emphasize that this perturbative theory applies whenever the external stress is not as large as to induce a discontinuous phase change in the material. Hence, the caloric response it captures is eminently reversible and does not include any latent-heat contribution. While, admittedly, the best mechanocaloric materials largely base their performance on the latent heat released (absorbed) at stress-driven first-order phase transitions [1][2][3][4], the present scheme allows us to inspect in detail the behavior within the range of continuous deformations, and is the key to the most interesting conclusions of this work. For example, thanks to this perturbative formalism, we can determine the dominance of the lowest-order temperature change (linear in the applied stress) for the moderate compressions used here (typically, −0.1 GPa), and thus predict (and explain) the coexistence of conventional and inverse mechanocaloric responses in the ferroelectric state.
Along the same lines, our perturbative theory plainly shows that the key to the predicted coexistence of conventional and inverse effects lies in the differentiated temperature dependence of the strains in the polar phase of PbTiO 3 : the strain parallel to the spontaneous polarization decreases upon heating, while the perpendicular strains increase. Interestingly, this feature is shared by many ferroelectrics [30][31][32]; and a similar anisotropic behavior of the strains is typical of compounds exhibiting a negative thermal expansion [26,33,34]. Thus, our analysis suggests that such compounds are likely to present a coexistence of conventional and inverse mechanocaloric responses.
We hope this work will encourage further investigations of the (potentially exotic and large) mechanocaloric properties of materials seldom considered to that end, such as ferroelectrics. We also hope that the simple perturbative formulas introduced here will be useful in future mechanocaloric studies, both theoretical and exper-