Negative electrocaloric effect in nonpolar phases of perovskite over wide range of temperature

TL;DR

This study reveals extensive negative electrocaloric effects in nonpolar phases of perovskite materials across a wide temperature range, driven by octahedral tilting entropy, with potential for high-efficiency refrigeration.

Contribution

It uncovers the microscopic mechanism of negative ECE in nonpolar perovskites and demonstrates a giant negative ECE at room temperature using first-principles Monte Carlo simulations.

Findings

Negative ECE observed in multiple perovskite phases.

Negative ECE driven by octahedral tilting entropy change.

Giant negative ECE of 8.6 K at room temperature.

Abstract

The electrocaloric effect (ECE) offers a promising alternative to the traditional gas compressing refrigeration due to its high efficiency and environmental friendliness. The unusual negative electrocaloric effect refers to the adiabatic temperature drops due to application of electric field, in contrast with the normal (positive) ECE, and provides ways to improve the electrocaloric efficiency in refrigeration cycles. However, negative ECE is unusual and requires a clear understanding of microscopic mechanisms. Here, we found unexpected and extensive negative ECE in nonpolar orthorhombic, tetragonal, and cubic phases of halide and oxide perovskite at wide range of temperature by means of first-principle-based large scale Monte Carlo methods. Such unexpected negative ECE originates from the octahedral tilting related entropy change rather than the polarization entropy change under the application of electric field. Furthermore, a giant negative ECE with temperature change of 8.6 K is found at room temperature. This giant and extensive negative ECE in perovskite opens up new horizon in the research of caloric effects and broadens the electrocaloric refrigeration ways with high efficiency.

Tables

Table 1. Table 1: Parameters for Landau-like model. Values are given in atomic units.

	$a_u$	$b_u$	$a_s$	$b_s$	$a_{\omega }$	$b_{\omega }$	$a_{\theta }$	$b_{\theta }$
\ceCsPbI3	-0.000486	0.00270	-0.000530	0.00270	-0.1017	0.9914	-0.1828	4.102
\cePbSc_0.5Ta_0.5O3	-0.0198	0.0519	-0.00369	0.0519	-0.1981	3.777	-0.4538	13.375
	$c_1$	$c_2$	$d$	$f_1$	$f_2$	$c_1/a_u$	$f_1/a_u$	$f_2/a_u$
\ceCsPbI3	0.0162	4.216	-0.2574	0.07393	0.07118	-33.341	-152.155	-146.496
\cePbSc_0.5Ta_0.5O3	0.3112	10.482	0.1620	0.8357	0.7556	-15.714	-42.196	-38.153

Equations

$$\begin{array}{rl}{H}_{\text{model}}& =a_u{u}^2+b_u{u}^4+a_s{s}^2+b_s{s}^4\\ & +a_{\omega }{\omega }^2+b_{\omega }{\omega }^4+a_{\theta }{\theta }^2+b_{\theta }{\theta }^4\\ & +c_1{u}^2{s}^2+c_2{\omega }^2{\theta }^2+ds\omega \theta \\ & +f_1{\theta }^2(u^2+s^2)+f_2{\omega }^2(u^2+s^2)\\ & +E_{\text{elas}}+E_{\text{polar-elas}}+E_{\text{tilt-elas}}\\ & -2Z^{\ast }Eu,\end{array}$$

$$\begin{array}{rl}{E}_{\text{polar-elas}}& =B_{1xx}{\eta }_a(u^2+s^2)+B_{1yy}({\eta }_a+{\eta }_b)(u^2+s^2)\\ & +B_{4yz}{\eta }_c(u^2+s^2),\end{array}$$

$$\begin{array}{rl}{E}_{\text{tilt-elas}}& =\frac{1}{2}{C}_{1xx}(2{\eta }_a{\theta }^2+{\eta }_b{\omega }^2)\\ & +C_{1yy}({\eta }_a{\theta }^2+{\eta }_a{\omega }^2+{\eta }_b{\theta }^2)\\ & +C_{4yz}{\eta }_c{\theta }^2.\end{array}$$

$$\begin{array}{r}\Delta S=\Delta S_{\text{polar}}+\Delta S_{\text{tilt}}+\Delta S_{\text{coup}}+\Delta S_{\text{field}},\end{array}$$

$$\begin{array}{rl}\Delta S_{\text{polar}}& =\frac{1}{T}[a_u\Delta (u^2)+b_u\Delta (u^4)\\ & +a_s\Delta (s^2)+b_s\Delta (s^4)+c_1\Delta (u^2{s}^2)\\ & +E_{\text{elas, H}}+E_{\text{polar-elas}}],\end{array}$$

$$\begin{array}{rl}\Delta S_{\text{tilt}}& =\frac{1}{T}[a_{\omega }\Delta ({\omega }^2)+b_{\omega }\Delta ({\omega }^4)\\ & +a_{\theta }\Delta ({\theta }^2)+b_{\theta }\Delta ({\theta }^4)+c_2\Delta ({\theta }^2{\omega }^2)\\ & +E_{\text{tilt-elas}}],\end{array}$$

$$\begin{array}{rl}\Delta S_{\text{coup}}& =\frac{1}{T}[d\Delta (s\theta \omega )+f_1\Delta ({\theta }^2(u^2+s^2))\\ & +f_2\Delta ({\omega }^2(u^2+s^2))],\end{array}$$

$$\Delta S_{\text{field}}=\frac{2Z^{\ast }}{T}\left[\int u\mathrm{d}E-\Delta (Eu)\right].$$

$$\Delta S_{\text{antipolar}}=\frac{1}{T}\left[a_s\Delta (s^2)+b_s\Delta (s^4)\right].$$

$$\begin{array}{rl}\mathrm{d}{H}_{\text{model}}& =\{2a_{u}u+2(c_1{s}^2+f_1{\theta }^2+f_2{\omega }^2)u+4b_u{u}^3\\ & +2[B_{1xx}{\eta }_a+B_{1yy}({\eta }_a+{\eta }_b)+B_{4yz}{\eta }_c]u-2Z^{\ast }E\}\mathrm{d}u,\end{array}$$

$${\left(\frac{\partial P}{\partial T}\right)}_E>0.$$

$$\alpha =-\frac{T}{C_E}{\left(\frac{\partial P}{\partial T}\right)}_E,$$

$${\omega }_{K,\alpha }^{\text{clus}}=\frac{1}{N_t}\sum _j|{\omega }_{R,\alpha }(j)|,$$

$$E_{\text{tot}}=E_{\text{dipole}}(\{{\mathbf{u}}_i\},\{{\mathbf{v}}_i\},{\eta }_H)+E_{\text{AFD}}(\{{\mathbf{u}}_i\},\{{\mathbf{v}}_i\},\{{\mathbf{\omega }}_i\},{\eta }_H),$$

$$\Delta S=\frac{1}{T}\left(\Delta H+\int \Omega \mathbf{P}\cdot \mathbf{E}\right),$$

$$H=E_{\text{tot}}+\frac{21}{2}N{k}_{B}T+p\Omega -\Omega \mathbf{E}\cdot \mathbf{P},$$

$$\begin{array}{rl}\Delta S{|}_{\mathbf{E}=0}& =\int \frac{1}{T}\mathrm{d}H\\ & =\Delta \left(\frac{H}{T}\right)+\int \frac{H}{T^2}\mathrm{d}T.\end{array}$$

$$E_{\text{dpl}}=\sum _{ij\alpha \beta }{Q}_{\alpha \beta }({\mathbf{r}}_i-{\mathbf{r}}_j)u_{\alpha }({\mathbf{r}}_i)u_{\beta }({\mathbf{r}}_j)$$

$$E_{\text{mode-strain}}$$

$$E_{\text{AFD-strain}}$$

$$E_{\text{mode-AFD,1}}$$

$$E_{\text{mode-AFD,2}}$$

Full text

I Introduction

Refrigeration consumes energy intensively and more than 20% of the electricity generated in the world is used for cooling [1]. Solid state cooling technologies based on the caloric effects display an attractive alternative to the traditional vapor compression cycles due to the higher operating efficiency and zero greenhouse gas emission [2, 3, 4]. Caloric effects are the phenomena of temperature and entropy change induced by electric field (electrocaloric, EC) [5, 6, 7], magnetic field (magnetocaloric) [8], mechanical stresses (elastocaloric [9], barocaloric [10] and twistocaloric [11]). EC based cooling technology is promising in modern microelectronics where millions of electronic units integration to the chip lead to tremendous power dissipation density and high-efficient thermal manage systems is highly required [12, 13]. EC effects can be positive (where the entropy decreases and temperature increases with the application of electric field), and negative (where the entropy increases and the temperature decreases with the application of electric field), as shown in Fig. 1. Both of the two EC effects can be used in cooling technologies [2, 14]. Combining the positive and negative ECE is proposed to enlarge the EC temperature span and improve the EC efficiency in refrigeration cycles [15, 14, 16, 17]. Though the positive ECE have been intensively studied in single crystals [18], thin films [19], ceramic multilayer chips [20, 21], polymer [22, 6, 23], and ceramic/polymer composite [24, 25], the negative ECE is found only in a few materials, basically complex polar ralxor and antiferroelectric oxdies [26, 27, 15, 28, 29]. The presenting temperature and magnitude of negative EC coefficient are limited, and its microscopic mechanism is highly required to be revealed.

In addition to oxide perovskites, halide perovskites are also promising candidate materials in cost-effective, high-performance electronics and optoelectronics [30, 31, 32]. The temperature change under electric-field will give great impact on performance, degradation and reliability of the materials and devices. However, the response under electric field in halide perovskite and related materials is little studied.

Here, we investigated the EC effect in nonpolar halide (\ceCsPbI3) and oxide (\ceBaCeO3) perovskite using first-principles-based large scale Monte Carlo (MC) simulation methods. Perovskite \ceCsPbI3 and \ceBaCeO3 can adopt nonpolar structures of orthorhombic, tetragonal, and cubic phases above room temperature. All the nonpolar phases possess negative ECE at small electric-field. The giant negative ECE with temperature change -8.6 K at room temperature is found, which is larger than the typically reported positive ECE with temperature change of +5.5 K [20] or negative ECE with temperature change of -5.76 K [26]. This large and extensive negative ECE comes from the entropy change related octahedral antiferrodistortive (AFD) motion (see Fig. 2a) and the strong coupling between AFD and polar motions.

II Results

II.1 Phenomenological Model

To understand the ECE in perovskites, the Landau-like model is constructed with four types of non-vanishing order parameters and the homogeneous strain. The four order parameters are (i) AFD motion at M point $\omega_{M,z}=\omega$, (ii) AFD motion at R point $\omega_{R,x}=\omega_{R,y}=\theta$, (iii) antipolar soft mode motion at X point $u_{X,x}=u_{X,y}=s$, and (iv) polar soft mode motion at $\Gamma$ point $u_{x}=u_{y}=u$ (which is fully driven by external electric field). According to the symmetry, the components of the homogeneous strain are $\eta_{1}=\eta_{2}=\eta_{a},\eta_{3}=\eta_{b},\eta_{4}=\eta_{c},\eta_{5}=\eta_{6}=\eta_{d}$. Therefore, using $\omega,\theta,s,u$ and $\eta_{a},\eta_{b},\eta_{c},\eta_{d}$ as order parameters, the Landau-like model could be written as

[TABLE]

where $\tilde{E}_{x}=\tilde{E}_{y}=E$ is the components of external electric field, and $E_{\text{elas}}$ is the energy associated with the homogeneous strain whose details can be found in Ref. [33], $E_{\text{polar-elas}}$ and $E_{\text{tilt-elas}}$ have following forms

[TABLE]

and

[TABLE]

Note that the parameters in this model of Eq. (1) can be directly derived from the parameters used in our effective Hamiltonian (see Appendix) which has been confirmed by first principles calculations [34]. For example, $b_{u}=b_{s}=4\alpha+\gamma$ and $f_{2}=2E_{xxyy}$, where $\alpha,\gamma,E_{xxyy}$ are the parameters of effective Hamiltonian calculations [34]. Table 1 shows the values of the parameters in the model of Eq. (1).

To further understand the ECE from the entropy concept quantitatively, the isothermal entropy change induced by applying electric field is calculated within the model of Eq. (1) using Eq. (15). The isothermal entropy change can be split into four parts

[TABLE]

where $\Delta S_{\text{polar}}$ is the entropy change only comes from the soft mode, antipolar motion, strain and their coupling

[TABLE]

$\Delta S_{\text{tilt}}$ is the contribution only from AFD motion and its coupling with strain

[TABLE]

$\Delta S_{\text{coup}}$ is the entropy change comes from the coupling terms between soft mode and AFD

[TABLE]

$\Delta S_{\text{field}}$ is the contribution directly related to the electric field

[TABLE]

The notation $\Delta(g)$ means the change of value $g$ in the isothermal process. Particularly, the entropy change associated with antipolar motion $\Delta S_{\text{antipolar}}$, which is also a part of $\Delta S_{\text{polar}}$, is given by

[TABLE]

II.2 Structures at finite temperature

Using first-principles-based effective Hamiltonian and hybrid Monte Carlo metheds (see Appendix A and B), we first investiagted structures and EC properties of halide perovskite structure of \ceCsPbI3 at finite temperature. Figure 2b shows average AFD motions of in-phase tilting at M point ($\omega_{M,z}$) and antiphase tilting at R point ($\omega_{R,x}$) in supercell as a function of temperature (see Fig. 2a about schematics of such motions, and the corresponding definition in Appendix). At high temperature larger than 420 K, the statistical averages of $\bm{\omega}_{M}$ and $\bm{\omega}_{R}$ are zero, characterizing a cubic $Pm\bar{3}m$ phase (C phase). As the temperature decreases to 420 K, $\omega_{M,z}$ shows a non-zero average value while $\bm{\omega}_{R}$ is still zero, characterizing a tetragonal $P4/mbm$ phase (T phase) with iodine octahedra tilting pattern $a^{0}a^{0}c^{+}$ (Glazer’s notation [35]). As for the temperature smaller than 390 K, both $\omega_{R,x}$ and $\omega_{M,z}$ have non-zero values. The $\omega_{R,x}$ and $\omega_{R,y}$ always have almost identical values below 390 K, characterizing a orthorhombic $Pnma$ phase with tilting pattern $a^{-}a^{-}c^{+}$ (O phase). The phase diagram is consistent with experimental measurement [36] and previous calculations [34]. Note that Figs. 2b,c also show the average absolute value of AFD $|\omega|$, which possesses a finite value at all investigated temperatures, even for the C phase where both $\bm{\omega}_{M}$ and $\bm{\omega}_{R}$ vanish, indicating the existence of disordered local AFD in C phase.

Figure 2c shows AFD vectors $\omega_{M,z}$ and $\omega_{R,x}$ as a function of temperature under electric field of $\tilde{E}$=0.43 MV/cm. This phase diagram is very similar to that in the absence of the electric field. The phase transition from C phase to T phase occurs at the temperature of 415 K, and the transition from T phase to O phase occurs at 385 K. Both of the transition temperatures are slightly smaller than that at zero electric field by 5 K. Interestingly, the polarization shown in Fig. 2c driven by external electric field decreases with the decreasing of temperature. This is different from prototypical ferroelectrics where the polarization increases with the decrease of temperature and becomes large at lower temperature. This unusual phenomenon of polarization driven by electric field with respect to temperature implies the abnormal temperature response under electric field. Note that due to the soft mode motion under finite electric field, the symmetry group of the three phases under the application of electric field are no longer exactly $Pm\bar{3}m$, $P4/mbm$ and $Pnma$, but adopt lower symmetry. However, the notations of C, T and O phase are still used, since the AFD motions at small fields are similar to that in the absence of electric field (see Fig. 2), and the soft mode motion is rather small. As the electric field increase above 2.0 MV/cm, no phase transition is observed, and the AFD vectors $\bm{\omega}_{R},\bm{\omega}_{M}$ remain zero in the whole investigated temperature range.

II.3 Negative ECE in broad range of temperature

Figure 3a shows the calculated EC coefficient $\alpha$ as functions of temperature $T$ and the applied electric field $\tilde{E}$ along the pseudocubic [110] direction, in which the white solid line is the isoline of zero EC coefficient $\alpha$. The EC coefficient is positive above this isoline (red color), indicating positive ECE, which generally exists in normal paraelectric materials [37]. The EC coefficient is negative in a large region below such isoline (blue color). Interestingly, this negative ECE exists in all the nonpolar phases at a large range of temperature from 300 K to $\simeq 550$ K, very different from that in the previous studies, where negative ECE is typically found in antiferroelectrics phase in a relatively narrow range of temperature [26, 27]. Figure 3b shows the calculated relative entropy diagram $\tilde{S}(T,\tilde{E})$, where the isentropic lines are depicted by the white solid lines. The adiabatic system evolves in one isentropic line in a reversible process. On tracing the isentropic lines in the direction of increasing field (i.e. from the bottom to the top in Fig. 3b), the left (respectively, right) bending of isentropic lines imply negative (respectively, positive) ECE. Figure 3c shows the maximal negative adiabatic temperature change ($\Delta\tilde{T}_{S}$) as a function of starting temperature, which is the temperature change ($\Delta\tilde{T}_{S}$) when tracing the isentropic lines in Fig. 3b under electric field. The magnitude of negative adiabatic temperature change is large (up to -8.6 K) at low temperature about 300 K, and decreases with the increase of temperature. The negative temperature change is larger than 2 K over a broad range of temperature from 300 K to 415 K.

From Eq. (1), as the external electric field only couples with soft mode motion directly, the enthalpy change induced by change of soft mode under small electric field (where negative ECE is observed) has the form

[TABLE]

where the last term is directly induced by the electric field and is constant under a given magnitude of electric field, and the third term could be neglected for small $u$. As shown in Table 1, the parameter $a_{u}$ for the first term is negative, indicating the instability of soft mode distortion. For the second term, the large positive $c_{1}$, $f_{1}$ and $f_{2}$ imply the strong coupling between AFD ($\theta$ and $\omega$), antipoloar ($s$) and soft mode ($u$) motion, and the antipolar and AFD motion tend to restrain the soft mode distortion. The fourth term typically has small magnitude comparing with the second one. The values of the same parameters for \cePbSc_0.5Ta_0.5O3 (PST) [38] which possesses positive ECE are also listed in Table 1 for comparison. One can see the $f_{1}/a_{u},f_{2}/a_{u},c_{1}/a_{u}$ for PST have much smaller magnitudes than those in \ceCsPbI3, indicating the coupling between soft mode motions and AFD as well as antipolar motion are not strong enough to suppress the soft mode motions in PST, which does not possess negative ECE.

The strong competition between soft mode motion and AFD motion ($f_{1}\theta^{2}u^{2}$ and $f_{2}\omega^{2}u^{2}$) can be confirmed by AFD and polarization in Fig. 2c. At the small electric filed 0.43 MV/cm, the magnitude of $\omega_{M,z}$ and $\omega_{R,x}$ increase with the decreasing of temperature in T and O phase, while the polarization decreases with the decreasing of temperature from about 500 K. Therefore, the derivative of polarization $P$ (or the soft mode $u$) with respective to temperature below $\simeq 500$ K at small constant field is positive

[TABLE]

This contrasts to the fact in prototypical ferroelectrics where the polarization increases with the decrease of temperature. As suggested by Maxwell relations, the EC coefficient could be written as [2, 39]

[TABLE]

where $C_{E}$ is the heat capacity at constant electric field and is positive. It is thus clear that the positive value of $\left(\frac{\partial P}{\partial T}\right)_{E}$ would result in negative value of $\alpha$ (that is negative ECE). Such analysis implies that the negative ECE below $\simeq 500$ K originates from the strong coupling between AFD and soft mode.

II.4 Negative ECE in O Phase

The structure of \ceCsPbI3 is O phase below the temperature 390 K without applying electric field. The $O$ phase possesses negative ECE (blue color below the white line in Fig. 3a) at small electric field. At the electric field above 2.0 MV/cm, the O phase transforms into C phase, and the ECE becomes positive. To further investigate the negative ECE, Fig. 4a shows the adiabatic temperature change $\Delta\tilde{T}_{S}$ due to the application of electric field, from the initial temperatures of 300 K, 340 K and 380 K. For the adiabatic process starting from 300 K, the temperature decreases slightly as the electric field increases from 0 MV/cm, which is consistent with the negative EC coefficient (see Fig. 3a). The temperature drop is about -0.82 K at 1.95 MV/cm. As the electric field further increases and passes through the phase boundary between O and C phases, the temperature change drops drastically to -8.6 K. The further increasing of electric field larger than 2.0 MV/cm results in the increase of the temperature change (positive ECE). The adiabatic process starting from 340 K behaves similarly to that starting from 300 K, possessing a temperature jump at 1.8 MV/cm with an extreme value of -7.5 K. For the starting temperature of 380 K, the temperature change drop occurs at the phase boundary line between O and T phase at lower electric field about 1.5 MV/cm, and exhibits a broader valley at the vicinity of phase boundaries.

To figure out the entropy change in the phenomena of negative ECE in O phase, Fig. 4b shows the isothermal entropy change computed from the Monte Carlo (MC) simulations by Eq. (15) [$\Delta S$ (MC)], the model of Eq. (4) [$\Delta S$ (model)], and the different terms in Eq. (4), as a function of electric field in an isothermal process at 300 K. The $\Delta S$ (model) curve consistents very well with the $\Delta S$ (MC) curve, indicating the validity of the model of Eq. (1) and Eq. (4). The $\Delta S_{\text{antipolar}}$ (Eq. (9)) is close to zero. The $\Delta S_{\text{field}}$ and $\Delta S_{\text{polar}}$ curves show negative values in the whole investigated field range, while the $\Delta S_{\text{tilt}}$ and $\Delta S_{\text{coup}}$ display positive values. The positive total entropy change mainly comes from the large $\Delta S_{\text{tilt}}$. Furthermore, the positive $\Delta S_{\text{coup}}$ is the entropy change from coupling between AFD and polar or antipolar motion. Therefore, the positive entropy change $\Delta S$ (model) and the resulting negative ECE come from the AFD related entropy change. The conclusion that the negative ECE is induced by the AFD related entropy change contrasts to the literature [26, 2] in which antipolar entropy change is proposed to leading to the negative ECE. Note that the AFD ($\omega$ and $\theta$) does not interplay with electric field directly (see the full analytical form of the effective Hamiltonian [34]), but couples with $u$ by $\omega^{2}u^{2}$ and $\theta^{2}u^{2}$, which leads to the AFD related entropy change due to $u$ (or polarization) responsible for the electric field indirectly.

On the other hand, near 2.4 MV/cm, the $\Delta S_{\text{tilt}}$ and $\Delta S_{\text{coup}}$ terms exhibit a sudden positive jump, leading to a jump of total entropy. In Fig. 4c, the enthalpy also exhibits a sudden jump near 2.4 MV/cm. Such jumps of entropy and enthalpy are caused by the phase transition from O phase to C phase driven by the external electric field, which is endothermic in nature. Further more, such phase transition takes place from an ordered state (with non-vanishing statistical average of $\omega_{M,z}$ and $\omega_{R,x}$ in O phase) to a disordered state (with vanishing order of AFD pattern as well as statistical average of $\omega_{M,z}$ and $\omega_{R,x}$ in C phase). The disordering of AFD motion contributes greatly to the overall entropy, resulting in the large negative ECE near the transition (Fig. 4a). This endothermic phase transition induced negative ECE also occurs in \cePbZrO3 [27].

In order to gain a microscopic understanding of the AFD entropy change, the local clusters AFDR,cl (or AFDM,cl) are defined to describe the local ordering of the AFD pattern. At finite temperature under electric field, the AFDs are not fully ordered, but partially ordered. The AFDR,cl (respectively, AFDM,cl) is defined as the local clusters in which the AFD vector in each unit cell $\bm{\omega}_{R}(i)=(-1)^{n_{x}(i)+n_{y}(i)+n_{z}(i)}\bm{\omega}_{i}$ [respectively, $\bm{\omega}_{M}(i)=(-1)^{n_{x}(i)+n_{y}(i)}$] are nearly parallel to each other. The local clusters AFDK,cl (where $K$ corresponds to $R$ or $M$ here) are practically identified by comparing directions of $\bm{\omega}_{K}(i)$ with their nearest neighbors [40], two neighbor AFD vectors are considered to belong to the same cluster if the cosine of the angle between them is larger than 0.85. The local clusters AFDK,cl can be described by size $\langle s\rangle$ (which is the number of units in the cluster) and the average magnitude of AFD vector in the cluster $\omega_{K,\alpha}^{\text{clus}}$ defined by

[TABLE]

where $\alpha$ is the cartesian direction index, $j$ sums over all the unit cells that belongs to a local AFD cluster with size larger than 1, and $N_{t}$ is the number of such unit cells in the simulated supercell. Note that the cluster can propagates from one side of the supercell to its opposite side, which is the so-called AFD percolation, the cluster of AFD percolation has infinite size under periodic boundary conditions.

As shown in Figs. 4d, e and f, the size $\langle s\rangle$ of the AFD cluster decreases when the electric field increase. Almost all unit cells belong to the clusters of percolation in the absence of field, some AFD vector do not belong to the AFD cluster and percolation at electric field of 2.33 MV/cm, and there is no percolation and most AFD vector do not belong to clusters at the electric field of 2.69 MV/cm. The ordered AFD vectors in the absence of field become disordered at high electric field, and the AFD entropy increases with the increase of electric field.

II.5 Negative ECE in T phase

The ECE in T phase is then investigated. As shown in Fig. 3a, the T phase possesses negative ECE at low electric field. The EC coefficient exhibits large negative value at the phases boundary of O, T, and C phases near the temperature of 390 K under electric field of 1.4 MV/cm. This large negative EC coefficient (dark blue color in Fig. 3a) exists in a wide range of electric field from 1.3 MV/cm to 2.0 MV/cm. Figure 5a shows the adiabatic temperature change as a function of electric field starting from the temperatures of 400 K, 410 K and 420 K, where the structure of \ceCsPbI3 is T phase at zero electric field. At low field, the larger the starting temperature is, the faster adiabatic temperature change ($\Delta\tilde{T}_{S}$) drops, and the smaller the maximal negative adiabatic temperature change is. For starting temperature of 400, 410 and 420 K, the maximal negative adiabatic temperature changes are -3.6 K at about 1.7 MV/cm, -3.1 K at about 1.6 MV/cm, and -1.6 K at about 1.5 MV/cm, respectively. Note that the maximal adiabatic temperature change terrace $\Delta\tilde{T}$ presents at larger range of temperature than that of O phase (see Fig. 4a), which is consistent with Fig. 3a.

To understand the unexpected negative ECE in T phase, the entropy change is analyzed with the model proposed in Eq. (1) and Eq. (4). In T phase the antipolar motion $s$ is zero and the antipolar entropy change $\Delta S_{\text{antipolar}}$ is zero. Note that the AFD at R point ($\theta$) is in principle zero in the whole supercell of T phase. However, considering the the AFD may be partially ordered and form local AFD cluster [41], here the $\theta$ is numerically defined as the average value of $\omega_{R}$ in AFDR,cl clusters, i.e. $\omega_{R,x}^{\text{clus}}$. Figure 5c shows the snapshot taken from MC simulation at 410 K at the absent of electric field, confirming the existence of AFDR,cl clusters in macroscopic T phase. It is numerically found such treatment of $\theta$ is essential for our model to reproduce the entropy change from MC calculations correctly. Figure 5b shows the isothermal entropy change computed from MC simulation [$\Delta S$ (MC)], the model proposed in Eq. (4) [$\Delta S$ (model)] and its contribution from each term, as a function of electric field at 400 K. The $\Delta S$ (model) is consistent very well with the $\Delta S$ (MC) at low field and also close to each other at high field. Similar to the case of O phase, the $\Delta S_{\text{field}}$ and $\Delta S_{\text{polar}}$ terms possess large negative values. The positive contribution mainly comes from $\Delta S_{\text{tilt}}$. It is thus clear that the fundamental origin of NECE in T phase is also the tilting entropy. Note that the contribution from $\Delta S_{\text{coup}}$ is rather small in T phase, in contrast with that in O phase, where the $\Delta S_{\text{coup}}$ has significantly positive value (see Fig. 4b). It is numerically found that such difference stems from the trilinear coupling term $d\Delta(s\theta\omega)$ of Eq. (7) in O phase. In fact, the trilinear term is the main contribution to $\Delta S_{\text{coup}}$ in O phase. However, such term vanishes in T phase, since the antipolar motion $s$ is zero in T phase, leading to the small total value of $\Delta S_{\text{coup}}$.

II.6 Electrocaloric Response in C Phase

C phase of \ceCsPbI3 also presents negative ECE. As shown in Fig. 6a, the adiabatic temperature change starting from 450 K and 500 K can be negative at low electric field. The negative adiabatic temperature starting from 450 K reaches the maximum of -0.36 K at the electric field of 1.37 MV/cm. This maximum of negative ECE occurs at the critical electric field where $\alpha=0$ shown in in Fig. 3a. The negative ECE become weaker when the starting temperature increases, characterized by the decreasing absolute values of temperature change under electric field. The negative ECE vanishes completely at the starting temperature of 550 K, where the adiabatic temperature change is positive for all investigated magnitudes of electric field, consistent with the fact that the $\alpha$ is almost always positive at temperature larger than 550 K. The ECE becoming positive at the temperature larger than 550 K is also consistent with the fact that $\left(\frac{\partial P}{\partial T}\right)_{E}$ is negative at low electric field (see Fig. 2c) which leads to the positive $\alpha$ from Eq. (12).

Similar to the local cluster AFDR,cl in T phase, there are partially ordered AFD of local clusters AFDR,cl and AFDM,cl in C phase. The average values of $\omega_{R}$ (respectively, $\omega_{M}$) inside the AFDR,cl (respectively, AFDM,cl) clusters are used as $\theta$ (respectively, $\omega$) in the models of Eq. (1) and Eq. (4) to analyze the entropy change. Figure 6b shows the isothermal entropy change computed from MC simulation [$\Delta S$ (MC)], the model in Eq. (4) [$\Delta S$ (model)] and its contribution from different terms, as a function of electric field at 450 K. Similar to T phase, the $\Delta S_{\text{field}}$ and $\Delta S_{\text{polar}}$ are negative, $\Delta S_{\text{coup}}$ is almost zero, $\Delta S_{\text{tilt}}$ is positive and very large. Clearly, it is $\Delta S_{\text{tilt}}$ that mainly results in the negative ECE in C phase, similar to that in T phase where the large positive tilting entropy change comes from the strong coupling between AFD and soft mode motions in local clusters under electric field.

III Discussion

All the nonpolar phases of O phase ($a^{-}a^{-}c^{+}$), T phase ($a^{0}a^{0}c^{+}$) and C ($a^{0}a^{0}a^{0}$) phase in \ceCsPbI3 possess negative ECE, originating from the AFD (highly ordered or partly ordered) entropy change. The result of AFD entropy changes resulting in negative ECE is different from the hypothesis in which antipolar entropy change gives rise to the negative ECE [26, 2]. In \ceCsPbI3, the antipolar entropy change is almost zero in O phase, and the there is no macroscopic antipolar motion in T and C phases. Our results of negative ECE in the paraelectric phase of T and C phases are in contrast with the calculations of the unified perturbation model [37] in which the paraelectric compounds possess positive ECE. This is because of the parameter of the AFD cluster and its coupling with soft mode, which leads to the the dielectric susceptibility $\chi$ decreases with the decreasing of temperature (Fig. 2a), in contrast with the normal Curie-Weiss law [37].

Such negative ECE comes from the AFD entropy change in nonpolar phases can also exist in other perovskites, including both halide and oxide peroskites. We take \ceBaCeO3 as an example to investigate the ECE in oxide perovskites. As suggested in previous measurements [42, 43] and first-principles-based calculations [44], \ceBaCeO3 possesses paraelectric $Pm\bar{3}m$ phase, nonpolar tetrahedral $I4/mcm$ phase (with AFD pattern $a^{-}a^{0}a^{0}$) and antipolar orthorhombic $Pbnm$ phase (with AFD pattern $a^{-}a^{-}c^{+}$) as the temperature decreases from 1200 K to room temperature in a monodomain form. The adiabatic temperature change are calculated for different phases of \ceBaCeO3 [44] as shown in Fig. 7. The structure phase is O phase ($a^{-}a^{-}c^{+}$) at the starting temperature of 400 K, and T phase ($a^{-}a^{0}a^{0}$) at the strating temperature of 800 K, C phase ($a^{0}a^{0}a^{0}$) at the starting temperature of 1200 K. All the phases exhibit negative ECE at certain small electric field range.

The novel phenomena of negative ECE is found in a very broad range of temperature in nonpolar phases of O, T and C phases for halide and oxide perovskites. These negative ECE is induced by the octahedra AFD related entropy change under the application of electric field and by the strong coupling between AFD and polarization in the nonpolar phases. We thus expect that the presently determined negative ECE would be measured as the important thermal control in the perovskite application in electronic, photoelectric and photovoltaic devices.

Appendix A Effective Hamiltonian

We use the effective Hamiltonian methods developed in Ref. [34], along with the parameters computed from first-principle calculations. The degrees of freedom of effective Hamltonian are local soft mode vectors $\{\bm{u}_{i}\}$ of each 5-atom unit cell, local displacement vectors $\{\bm{v}_{i}\}$ related to inhomogeneous strain, pseudo vectors $\{\bm{\omega}_{i}\}$ representing the rotation of iodine octahedra [also known as antiferrodistortive (AFD) motions], and homogeneous strain tensor $\eta_{H}$. The total energy of the effective Hamiltonian has two main terms

[TABLE]

where $E_{\text{dipole}}$ is the energy of local soft modes, strains and their coupling, including five main terms [33], namely, the quartic soft mode self energy, the quadratic soft mode long range dipolar energy, the quadratic soft mode short range interaction energy (up to the third nearest neighbor), the elastic energy (including both homogeneous and inhomogeneous contributions), and the interaction between soft mode and strain (quadratic in soft mode and linear in strain); the $E_{\text{AFD}}$ is the energy of AFD motions and their coupling with strains and soft modes, including the quartic AFD onsite energy, the quadratic and quartic AFD short range interaction energy (up to the first neighbor), the interaction between AFD and strain (quadratic in AFD and linear in strain), and the biquadratic and trilinear interaction between soft mode and AFD motion. The interaction coefficient are greatly simplified due to the symmetry. The complete analytical form of the effective Hamiltonian is proposed in Ref. [34]. An additional term $-Z^{*}\sum_{i}\bm{u}_{i}\cdot\bm{E}$ is added in $E_{\text{tot}}$ to consider the effect of external electric field [45], where the summation of $i$ runs over all the unit cells in the simulated supercell, $Z^{*}$ is the Born effective charge of the soft mode, and $\bm{E}$ is the external electric field. Note that $\tilde{\bm{E}}$ is the rescaled field according the MC calculations and the measurements [39, 46, 38].

The perovskite structures are typically simulated by $12\times 12\times 12$ supercells (corresponding to 8640 atoms) or large $32\times 32\times 32$ supercells (corresponding to 163840 atoms) with Monte Carlo (MC) simulations. The following quantities are computed: (i) The AFD at R point characterizing anti-phase tilting, defined as $\bm{\omega}_{R}=\frac{1}{N}\sum_{i}\bm{\omega}_{i}(-1)^{n_{x}(i)+n_{y}(i)+n_{z}(i)}$, where the summation of $i$ runs over all of the $N$ five-atom perovskite unit cells in the simulated supercell, $\bm{\omega}_{i}$ is the AFD vector of unit cell located at $a_{0}(n_{x}(i)\hat{\bm{x}}+n_{y}(i)\hat{\bm{y}}+n_{z}(i)\hat{\bm{z}})$, $a_{0}$ is the lattice constant of the five-atom perovskite unit cell, and $\hat{\bm{x}},\hat{\bm{y}},\hat{\bm{z}}$ are the unit vectors along the pseudocubic $[100]$, $[010]$ and $[001]$ directions, respectively; (ii) the AFD at M point characterizing the in-phase tilting, defined as $\bm{\omega}_{M}=\frac{1}{N}\sum_{i}\bm{\omega}_{i}(-1)^{n_{x}(i)+n_{y}(i)}$; (iii) the average absolute value of local AFD vectors $|\omega|=\frac{1}{N}\sum_{i}|\bm{\omega}_{i}|$; (iv) the average soft mode motion $\bm{u}=\frac{1}{N}\sum_{i}\bm{u}_{i}$; (v) the polarization $\bm{P}=\frac{1}{\Omega}Z^{*}\bm{u}$, where $\Omega$ is the volume of the simulated supercell; and (vi) the antipolar soft mode motion at X point $\bm{u}_{X}=\frac{1}{N}\sum_{i}\bm{u}_{i}(-1)^{n_{z}(i)}$.

The dielectric susceptibility tensor $\chi$ shown in Fig. 2b and the EC coefficient $\alpha$ shown in Fig. 3a are computed from the MC simulations with the cumulant formula [47, 39, 48].

The isothermal entropy change is calculated using the numerical method developed in Ref. [38] as

[TABLE]

where $H$ is the enthalpy calculated by

[TABLE]

where $N$ is the number of unit cells in the supercell, $k_{B}$ is the Boltzmann constant, $T$ is the temperature, $E_{\text{tot}}$ is the energy of effective Hamiltonian provided in Eq. (14), $p$ is the external pressure ($p=0$ in this work), $\bm{E}$ is the external electric field, and $\bm{P}$ is the polarization of the supercell. The entropy change as a function of temperature with absent of electric field is computed by

[TABLE]

The relative entropy diagram $\tilde{S}(T,\tilde{E})$ is thus determined from Eqs. (15) and (17), as shown in Fig. 3b. The adiabatic temperature change $\Delta\tilde{T}_{S}$ is calculated from the relative entropy diagram. More precisely, the adiabatic temperature change $\Delta T_{S}(T_{\text{start}},\tilde{E})$ from certain temperature $T_{\text{start}}$ induced by applying certain electric field $\tilde{E}$ is determined by finding the temperature $T_{1}$ at the given field $\tilde{E}$ that satisfies $\tilde{S}(T_{1},\tilde{E})=\tilde{S}(T_{\text{start}},0)$, and the temperature change is then calculated by $\Delta\tilde{T}_{S}(T_{\text{start}},\tilde{E})=T_{1}-T_{\text{start}}$.

Appendix B Hybrid Monte Carlo Simulation

Hybrid Monte Carlo (HMC) algorithm [49, 50, 51] is implemented for large scale simulations. In each Monte Carlo sweep (MCS) of the HMC simulations, a new trial configuration of $\{\bm{u}_{i}\},\{\bm{v}_{i}\}$ and $\{\bm{\omega}_{i}\}$ is generated by performing microcanonical molecular dynamics (MD) simulation for a short period with random initial momenta, while the homogeneous strain $\eta_{H}$ is practically updated using the standard Metropolis algorithm as in Ref. [33].

Technically, the long range dipole energy

[TABLE]

and forces associated with it are computed in the reciprocal space with the help of fast Fourier transformation algorithm, as in Refs. [51, 52, 53]. Such treatment reduces the overall computation complexity of one MCS from $O(N^{2})$ in standard Metropolis algorithm to $O(N\log N)$. Moreover, in the MD simulations, the degrees of freedom of all unit cells are updated simultaneously, in contrast with Metropolis MC simulation, in which the degrees of freedom are updated in sequence. Such simultaneous update makes it easy to run the HMC simulations in parallel, especially on shared memory architechtures [51]. The reduced complexity and ability to run in parallel makes it possible to perform simulations on large supercells efficiently.

Moreover, computational efficiency tests were performed for the $E_{\text{elas,I}}$ in Ref. [33] and the following terms (in notations of Refs. [33, 34])

[TABLE]

where the summation of $i$ runs over all the unit cells in the supercell, the summation of $j$ runs over several neighbor cells around $i$, $\{\eta_{I}\}$ is the inhomogeneous strain which is computed from $\{\bm{v}_{i}\}$ [33], $\alpha,\beta,\gamma,\delta$ are cartesian directions, and $l=1,2,\cdots,6$ is the Voigt notation index. Such tests show that it is more efficient to compute these terms (as well as their associated forces) in the reciprocal space, although they have in principle complexity of $O(N)$ calculating in the real space.

Acknowledgements.

X.M., Y.Y. and D.W. thank the National Key R$\&$D Programs of China (grant NOs. 2020YFA0711504, 2022YFB3807601), the National Science Foundation of China (grant NOs. 12274201, 51725203, 51721001, 52003117 and U1932115) and the Natural Science Foundation of Jiangsu Province(grant NO. BK20200262). We are grateful to the HPCC resources of Nanjing University for the calculations.