Giant room temperature elastocaloric effect in metal-free thin-film perovskites

Solid-state refrigeration which is environmentally benign has attracted considerable attention. Mechanocaloric (mC) materials, in which the phase transitions can be induced by mechanical stresses, represent one of the most promising types of solid-state caloric materials. Herein, we have developed a thermodynamic phenomenological model and predicted extraordinarily large elastocaloric (eC) strengths for the (111)-oriented metal-free perovskite ferroelectric [MDABCO](NH4)I3 thin-films. The predicted room temperature isothermal eC ΔSeC/Δσ (eC entropy change under unit stress change) and adiabatic eC ΔTeC/Δσ (eC temperature change under unit stress change) for [MDABCO](NH4)I3 are −60.0 J K−1 kg−1 GPa−1 and 17.9 K GPa−1, respectively, which are 20 times higher than the traditional ferroelectric oxides such as BaTiO3 thin films. We have also demonstrated that the eC performance can be improved by reducing the Young’s modulus or enhancing the thermal expansion coefficient (which could be realized through chemical doping, etc.). We expect these discoveries to spur further interest in the potential applications of metal-free organic ferroelectrics materials towards next-generation eC refrigeration devices.


INTRODUCTION
Solid-state refrigeration has attracted considerable attention due to its potential applications in small-scale cooling technology 1-3 . There are several physical effects that can be utilized for solidstate refrigeration, such as magnetocaloric (MC) 4,5 , electrocaloric (EC) [6][7][8][9] , and mechanocaloric (mC) effects 10,11 . In MC and EC effects, a temperature change is induced by applying a magnetic field and electric field, respectively. However, the MC and EC strengths for typical ferromagnetic/ferroelectric materials are relatively small. For instance, the EC strengths for prototype ferroelectrics such as PbZr 0.95 Ti 0.05 O 3 12 and fluorinated polymers 13 thin films are less than 1.0 K m MV −1 . To achieve a large temperature change, a strong magnetic/electric field is required, which is normally provided by a huge auxiliary device, limiting the practical applications of the MC/EC materials in cooling devices 9,14,15 .
The mC effect 10 , the application of a mechanical force that induces thermal cooling, can be divided into two categories, elastocaloric (eC) and barocaloric (BC) effects, for which the driving mechanical forces are uniaxial stress and hydrostatic pressure, respectively. The schematic illustration of an eC refrigeration cycle for thin films is shown in Fig. 1. The refrigeration cycle begins with applying a compressive stress on the film adiabatically (step a), the temperature rises from T 0 to T max due to the release of latent heat. Then followed by the heat transfer from thin films to the surroundings (step b), with a decrease of temperature while the compressive stress is kept constant. The temperature keeps dropping when the compressive stress is gradually reduced (step c), as a result of the inverse eC effect. Finally, at step d, the compressive stress is fully removed. The film absorbs heat from the surroundings where the temperature rises back to T 0 after a full refrigeration cycle.
To maximize the temperature change and improve the energy efficiency for the cooling cycle, materials with high eC strength are required. Previously, high eC strength with ΔT eC /Δσ of 0.7 K GPa −1 has been predicted for ferroelectric oxide heterostructures such as BaTiO 3 thin films 1 . One key factor limiting the eC performance of the oxide ferroelectrics is their low compressibility with strong ionic bonds. Even though a high adiabatic temperature change might be obtained by applying large external stress, the eC strength is quite small which gives rise to low energy efficiency for a refrigeration cycle.
Recently  3 . In addition to the smaller mass density, the molecular ferroelectric material has a much smaller elastic stiffness (~1 GPa) in comparison to typical ferroelectric oxides (~100 GPa). Therefore, it is rational to believe that large volume and entropy changes could be induced by applying external mechanical stimuli to [MDABCO](NH 4 )I 3 , which may generate a significant mC effect.
The very large BC effect was also reported in plastic crystals which are featured by the extensive disorder and giant compressibility 22,23 . For example, The BC strength ΔS/Δσ for NPG plastic crystals can reach 961.5 J K −1 kg −1 GPa −1 at 320 K 22 . The metal-free organic perovskite ferroelectric [MDABCO](NH 4 )I 3 is held together by both ionic and hydrogen-bonding interactions, which has considerably large compressibility. Large piezoelectric response was also predicted in [MDABCO](NH 4 )I 3 from firstprinciples calculations 24 . Moreover, [MDABCO](NH 4 )I 3 also has advantages of high polarization, mechanical flexibility, low weight, and low processing temperatures [25][26][27][28][29][30][31] , as compared to other organic ferroelectric systems.
Computational methods such as thermodynamic calculations have emerged as powerful tools to the understanding of eC effects, 1,32-36 which could help to explain and design materials with high eC effects. In this work, we employ thermodynamic calculations to study the polarization properties and phase transitions of (111)-oriented ferroelectric [MDABCO](NH 4 )I 3 thin films under different combinations of in-plane misfit strains and out-of-plane stresses, from which the eC properties are calculated. As compared to the bulk system, thin films are much smaller in size which is beneficial to the miniature of the device, while also providing an extra degree of freedom (the substrate strain) in the design of the cooling device. It is revealed that the isothermal ΔS eC /Δσ and adiabatic ΔT eC /Δσ for a [MDABCO](NH 4 )I 3 thin film under an out-of-plane compressive stress of 1.28 GPa and in-plane isotropic misfit strain of 0.01 at 300 K can reach as high as −60.0 J K −1 kg −1 GPa −1 and 17.9 K GPa −1 , respectively, much superior to the traditional ferroelectric oxides. The giant eC effect is found to be attributed to the low Young's modulus and high thermal expansion coefficient of [MDABCO](NH 4 )I 3 .

Phase transitions and phase diagram
The detailed procedure for developing a thermodynamic potential energy density function for (111)-oriented [MDABCO](NH 4 )I 3 film is described in the Theoretical Method Section and Supporting Information. The parameters for this study are all taken from the previous report 16 . We first calculate the analytical equilibrium polarization for the (111)-oriented [MDABCO](NH 4 )I 3 thin film under different in-plane misfit strains and out-of-plane stresses at room temperature. Consider a single-domain (111)-oriented thin film epitaxially grown on a cubic substrate with in-plane coherency, the schematics is given in Fig. 2a. The stress and strain components in the global coordinate system satisfy the mixed mechanical boundary condition: The calculated stress and misfit strain intervals are 0.01 GPa and 0.001, respectively. Supplementary  Fig. 2 shows the contour plots for the equilibrium polarization components P x , P y , and P z in the pseudocubic system as a function of the misfit strain and external stress. It can be observed that the P x , P y , and P z components have finite values at the low-stress region (see the R and M parts of the contour plot). In the high stress and high misfit strain regions, the cubic phase is stable where all polarization components vanish (the C-region in the phase diagram). In the O region (orthorhombic phase), the P x and P y components are in the same magnitudes while P z equals zero. Figure 2b shows the stable phases for (111)-oriented [MDABCO] (NH 4 )I 3 thin films under different in-plane misfit strains and out-ofplane stresses at room temperature. There are four stable phases in total, depending on the misfit strain and applied stress. When the in-plane misfit strains are zero, the stable phase is orthorhombic O (P x 2 = P y 2 ≠ 0, P z 2 = 0) for σ 33 < −2.0 GPa, monoclinic M (P x 2 = P y 2 > P z 2 ≠ 0) for −2.0 GPa ≤ σ 33 ≤ 0.0 GPa, and rhombohedral R (P x 2 = P y 2 = P z 2 ≠ 0) for σ 33 ≥ 0.0 GPa, respectively. In the region where both the in-plane misfit strain and the out-of-plane stress are highly negative, the cubic phase C (P x 2 = P y 2 = P z 2 = 0) is stable. In this case, the mechanical condition applied to the material is similar to a hydrostatic pressure, which could favor the paraelectric cubic phase. Notably, no stable T-phase can be observed in the entire phase diagram. This can be understood since R-phase is the ground state for [MDABCO](NH 4 )I 3 and the film is (111)-oriented. T-phase has much higher energy as compared to the O phase according to the calculations from Wang et al. 16 , which is difficult to stabilize via strain engineering in (111)-oriented [MDABCO](NH 4 ) I 3 thin films.

Elastocaloric effect
With the knowledge of the dependence of stress and epitaxial strain on the polarization, now we proceed to calculate the eC  properties. Figure 3 shows the eC entropy change ΔS eC and temperature change ΔT eC versus the out-of-plane compressive stress for various [MDABCO](NH 4 )I 3 thin films under different misfit strains (u m ). To calculate the variations of the ΔS eC and ΔT eC with the phase change, compressive stress is applied on top of the thin film. With the increase of the applied stress, phase transitions from rhombohedral (R) to monoclinic (M) or cubic (C) and finally to orthorhombic (O) occur. According to the previous report 37 , [MDABCO](NH 4 )I 3 undergoes a plastic deformation when the unidirectional compressive stress is above 1.28 GPa, which is set as the upper limit of the applied stress for the calculations in this study. Note that to sustain 1.28 GPa for the whole thin film device, one has to find either suitable oxide wafers or ductile metals as substrate, or even freestanding film without substrate. Assuming a constant heat capacity of 1000 J K −1 kg −1 for [MDABCO](NH 4 )I 3 at room temperature 16 , the maximum adiabatic temperature change ΔT eC~1 8.7 K occurs with u m = −0.01. In most cases, with the highest applied stress, the adiabatic temperature change can reach 10-20 K as shown in Fig. 3a, while the maximum adiabatic temperature change is 22.9 K with a tensile misfit strain of 0.01 as shown in Table 1. Interestingly, at the large compression side (e.g., u m = −0.05 and −0.03), ΔS eC (ΔT eC ) shows a nonlinear relationship with the external stress while for all the other cases mainly linear dependency is obeserved.  3 and other ferroelectrics that has been reported recently. For the polyvinylidene fluoride (PVDF) based ferroelectric polymers, the polar to the nonpolar phase transition is very sensitive to external mechanical stresses. High eC effects have been reported and adiabatic temperature changes of~2 K were measured directly in PVDF polymers near room temperature under uniaxial stresses of~10 MPa 38,39 . Meanwhile, some typical ferroelectric perovskite thin films also present finite eC effects at room temperature. For instance, BaTiO 3 thin films show a ΔT ec5 K under huge uniaxial stress of 6.5 GPa 1,21 . On the other hand, PbTiO 3 thin film exhibit a maximum adiabatic temperature change of~20 K 18 , which has been ascribed to the mixed ferroelastic and ferroelectric behaviors. However, such a giant eC effect can only occur at~700 K, which is too high for practical applications. In comparison, [MDABCO](NH 4 )I 3 displays considerably high eC strength, which is much higher than traditional lead-free ferroelectric perovskites, as shown in the light gray area of Supplementary Fig. 4. Polymer-based ferroelectrics have been found to exhibit higher temperature change ΔT ec , but they have a very low refrigerant capacity (RC = |ΔT*ΔS|). The calculated RC performance of [MDABCO](NH 4 )I 3 thin film is 1760 J kg −1 , much superior to other ferroelectrics as reported in the previous literature. The coefficient of performance (COP) is another commonly used metric to describe the efficiency of the cooling   Figure 4 demonstrates the eC entropy change ΔS eC , and the eC temperature change ΔT eC for [MDABCO](NH 4 )I 3 thin films under different in-plane misfit strains (u m ) and out-of-plane stresses. With the increasing of the out-of-plane stress, both the magnitudes of ΔS eC and ΔT eC increase. Figure 4c, d show the contributions of volume change and polarization change to ΔS eC and ΔT eC , respectively. We choose a compressive stress of −1.28 GPa as an example, two contributions to ΔS eC can be identified as follows. The ordering of the dipoles by the application of external stress under isothermal conditions leads to a reduction or increase in the dipole configurational entropy (ΔS P ), which is associated with the phase transitions. The other contribution is the volume change induced entropy change (ΔS V ) given by β l ρ σ 1 þ σ 2 þ σ 3 ð Þ , where the thermal expansion coefficient β l and the mass density ρ, have both been evaluated in the previous paper 16 . For [MDABCO](NH 4 )I 3 film, these two parts compete with each other when compressive stress is applied but the ΔS V dominates the total entropy due to its characteristic of low Young's modulus, which is only one percent of the Young's modulus for typical ferroelectric oxides. The low Young's modulus prompts the high thermal expansion coefficient and thus results in the big entropy reduction when applying the compressive stress (see Fig. 4a, c).

Calculations of elastocaloric effect for inorganic oxide ferroelectrics
Moreover, the impact of physical properties in materials on the eC strength is investigated, by calculating some common (111)oriented ferroelectric oxide films using the same method. Material parameters are taken from existing publications for BaTiO 3 (BTO) 40

DISCUSSION
Following the previous results, we further identify several intrinsic physical properties, namely the elastic compliance s ij , thermal expansion coefficient β l on the eC effects of the metal-free ferroelectrics, hoping to find design principles for high eC materials. The calculated eC temperature changes ΔT eC is plotted in Fig. 6 and Supplementary Fig. 6, by varying the values of the above-mentioned material properties. Figure 6a shows a twodimensional plot of the temperature change ΔT eC as a function of both elastic compliance and thermal expansion. It is indicated that generally, eC decreases with increasing of elastic compliance, while increases with increasing thermal expansion coefficient. A relatively stronger dependency on β l can be observed. For instance, ΔT eC shows a dramatic increase to~100 K with increasing of the β l to five times of [MDABCO](NH 4 )I 3 with fixing of the elastic compliance ( Supplementary Fig. 6).
Note that the elastic compliance and thermal expansion coefficient are not independent physical parameters for a certain material. We map out the elastic compliance and thermal expansion coefficients of several oxide ferroelectrics, PVDF, and [MDABCO](NH 4 )I 3 in Fig. 6b. The thermal expansion coefficient increases along with the increases of elastic compliance for almost all ferroelectric materials. Compared to oxide ferroelectric films, metal-free ferroelectric [MDABCO](NH 4 )I 3 film has higher elastic compliance that will surpress the eC effect, but the higher thermal expansion coefficient has a higher influence on the eC performance, resulting in one order of magnitude higher in eC strength. Furthermore, modifying the cation [MDABCO] 2+ and anion I − in the [MDABCO](NH 4 )I 3 crystal could change the physical properties including the thermal expansion coefficient and elastic compliance 25,45 . Thus, theoretically, by rational design of the composition of the [MDABCO](NH 4 )I 3 film, the eC performance in ferroelectrics can be further improved by enhancing the thermal expansion coefficient of the organic ferroelectrics.
In summary, we have developed a thermodynamic description for (111)-oriented organic ferroelectric [MDABCO](NH 4 )I 3 thin film and investigated its phase transition behavior and eC properties under different misfit strains and applied external stresses. We found that there are two contributions to the eC strength, namely the polarization changes induced entropy change (ΔS P ) and the volume change induced entropy change (ΔS V ), which competes with each other under compressive stress. For organic ferroelectric materials which are relatively soft, the entropy change ΔS V dominates, which is an order of magnitude higher than traditional ferroelectric oxides. It   3 , which is only one percent of Young's modulus of conventional oxides. By analyzing the relationship of eC with thermal expansion coefficient and elastic compliance, we further discovered that the thermal expansion coefficient has a dominant role in determining the eC. Based on our theoretical analysis, we predict that it is possible to further improve the eC properties by enhancing the thermal expansion coefficient through chemical modification of the [MDABCO](NH 4 )I 3 . We envision this work will attract broad interest from the community of solid eC materials.

METHODS Thermodynamic calculations of ferroelectric thin films
Based on the Landau-Devonshire theory, there are mainly four contributions in the total thermodynamic free energy density: (1) Polarization energy G polarization ¼ α i P i þ α ij P i P j þ α ijk P i P j P k þ α ijkl P i P j P k P l þα ijklm P i P j P k P l P m þ α ijklmn P i P j P k P l P m P n ; (2) thermal energy G thermal ¼ Àc Where α i , α ij , α ijk , α ijkl , α ijklm , and α ijklmn are the Landau coefficients, and c p , C ijkl , K ij , ε 0 , and E i represent the heat capacity, elastic modulus tensor, background dielectric constants, permittivity of vacuum, and electric fields. P i , T, and T 0 are the polarization components, temperature, and romm temperature, respectively. The eigen strain is defined as ε ij ¼ Q ijkl P k P l þ R ijklm P k P l P m þ β l ðT À T 0 Þ Here Q ijkl and R ijklm are the electromechanical strictive tensors. β l is the linear thermal expansion coefficient. Consider a paraelectric-ferroelectric phase transition from high symmetry cubic phase (with space group P432) to rhombohedral phase (space group R3), expanding the polynomials to eighth order, neglecting all the symmetry forbidden terms, one can rewrite: Where the α-terms are the simplified Landau coefficients, s 11 , s 12 , and s 44 are the elastic compliance, σ i ði ¼ 1; 6Þ are the stresses, λ eff is the high-order electrostrictive term. E i are the applied electric fields, K ij are the relative dielectric constants. We then transfer the polarization and stresses from the local coordinate system x to the global coordinate system x′ by using the transformation matrix t ij : For the epitaxial thin film with the in-plane misfit strains of ε 0 1 and ε 0 2 (the strain in global coordinate), the free energy ΔG′ can be derived using the Legendre transformation In Eq. 3, the misfit strains ε 0 1 and ε 0 2 are the same for a cubic substrate. For the epitaxial film under typical in-plane biaxial misfit strains and outplane stress in a global coordinate system, it is subjected to a mixed mechanical boundary condition with in which the ε 0 1 and ε 0 2 are the in-plane strains, the ε 0 6 is the shear strain, the σ 0 3 is the out-plane stress along [111] direction, and the σ 0 4 and σ 0 5 are the in-plane stresses.
Assuming the film is pre-polarized in the electric field, the polarization ( Pi ) under given misfit strains and out-plane applied stress at room temperature can be obtained by minimizing the thermodynamic free energy density ΔG′ with respect to polarization Pi .

Establishment of phase diagrams
The phase evolutions under applied stress at different misfit strains for [MDABCO](NH 4 )I 3 thin films are thoroughly investigated and exhibited in Fig. 2 after we obtained the magnitude of the polarization components.
Here we also employ the thermodynamic method to explore the phase evolutions for common ferroelectrics including BaTiO 3 , PbTiO 3 , PbZr 0.52 Ti 0.48 O 3 , PbZr 0.8 Ti 0.2 O 3 , K 0.5 Na 0.5 NbO 3 , and KNbO 3 . The stress versus misfit strain phase diagrams for these ferroelectric films at room temperature are displayed in Supplementary Fig. 3.

Calculation of elastocaloric effect
The total isothermal eC entropy change per mass under an application of stress σ app at room temperature can be calculated by ΔS eC ¼ À 1 ρ ∂ΔG 0 ∂T À Á σij and it is derived from two parts: (1) the polarization configuration change and (2) the stress-induced volume change Under the reversible adiabatic condition, the total entropy is conserved upon a stress change, thus we can obtain i j σ¼σapp À P 6 i j σ¼0 À Tβ l cpρ σ 1 þ σ 2 þ σ 3 ð Þ : All the physical parameters and all coefficients can be found in Supplementary Tables 1 and 2 in the Supporting Information.
The COP can be evaluated by the following equation: in which Q represents the absorbed heat and W is the specific mechanical work in the eC cooling cycle. The absorbed heat Q is derived from ∫TdS, in which the T and S are the temperature and entropy, respectively. The specific mechanical work W is calculated by 1 ρ R σdε, in which the ρ, σ, and ε are the mass density, stress, and strain, respectively.

DATA AVAILABILITY
The raw data for the thermodynamic calculations in this paper and its supplemental information files are available from the corresponding author (hongzijian100@zju. edu.cn) upon reasonable request.

CODE AVAILABILITY
The code for the thermodynamic calculations in this paper and its supplemental information files are available from the corresponding author (hongzijian100@zju. edu.cn) upon reasonable request.