A unified perturbative approach to electrocaloric effects

The electrocaloric effect, that is, the temperature change experienced by an insulator upon application of an electric field, offers promising ecofriendly alternatives to refrigeration. However, the theoretical treatments of this response are mostly case specific and lack a unified picture revealing the similarities and differences among the various known effects. Here, we show that the electrocaloric effect lends itself to a straightforward interpretation when expressed as a Taylor series in the external field. Our formalism explains in a unified and simple way the most notable small-field effects reported in the literature, namely the so-called normal and inverse electrocaloric responses, corresponding to an increase or decrease of temperature under applied field, as usually found in ferroelectrics or antiferroelectrics, respectively. This helps us to clarify their physical interpretation. We then discuss in detail atomistic simulations for the prototype ferroelectric PbTiO3, explicitly evaluating subtle predictions of the theory, such as the occurrence of competing contributions to the electrocaloric response. Theoretical treatments of electrocaloric effects, interesting for their potential use in refrigeration, are mostly case specific. Here, a perturbative approach provides a unified physical understanding of the normal and inverse electrocaloric response occurring in ferroelectrics and antiferroelectrics.

A n electric field can be used to change the temperature of an insulator, by virtue of the so-called electrocaloric effect. Electrocaloric effects hold the promise of an ecofriendly alternative to refrigeration, one of the most energy-consuming activities today and in the foreseeable future. Electrocaloric effects are strong in ferroelectric and antiferroelectric materials because of their anomalously large dielectric responses and the fielddriven phase transitions that can be easily induced near the Curie point 1 .
The thermodynamic description of the electrocaloric cooling cycle is well-established 1 . The key step is the adiabatic temperature change, which can be expressed as where the integral runs from zero to a final electric field E α , applied along the α Cartesian direction. S, T, and C E are, respectively, the entropy, temperature, and constant-field heat capacity. P α is the polarization component conjugate to the applied field E α , and its T-derivative is the pyroelectric coefficient π α . Note that we have used the Maxwell relation between the field derivative of the entropy and the pyroelectric vector, which is expected to hold for ergodic materials. (Thus, the formalism introduced here is, in principle, not applicable to relaxor ferroelectrics.) Note also that the repeated α index does not involve an implicit sum; we write this α, on both sides of the equation, to emphasize the fact that the electrocaloric response can be anisotropic. Also important is the isothermal entropy change which quantifies the amount of heat the electrocaloric material will exchange with the object to be cooled. The product ΔT(E α )ΔS (E α ) is usually taken as the figure of merit for electrocaloric cooling performance. These expressions are deceivingly simple, as all the quantities in the integrands depend on both field and temperature (itself field-dependent, by virtue of the electrocaloric effect). Further, they may present complex behaviors (in particular, discontinuities) across the fieldand temperature-driven phase transitions.
There are well-known strategies to solve these integrals selfconsistently 1 ; this will not be our focus here. Rather, we want to examine the ingredient that has the greatest influence on the basic features (magnitude and sign) of the electrocaloric temperature and entropy changes, namely, the pyroelectric vector π. Indeed, π fully determines ΔS(E α ); further, of the quantities contributing to ΔT(E α ), π is likely to be the most sensitive to an electric field and, more importantly, it is the only one whose sign is not defined. (Both T and C E are always positive.) Interestingly, we can rewrite Eq. (1) in differential form as This expression is the basis for the physical interpretation of electrocaloric effects 1 . For example, it is experimentally observed that, in ferroelectrics, the application of an electric field usually results in a positive temperature change, which is attributed to the fact that electric fields create order (reduce the entropy) in such systems 1 . In contrast, negative values of ΔT have been observed in antiferroelectrics 2-5 , or when a field is applied perpendicular to the polarization of a ferroelectric phase 6 , which is compatible with the intuitive notion that, in such cases, the electric field will destabilize the equilibrium state and cause disorder (increase the entropy). These diverse observations have been reproduced by atomistic simulations [7][8][9][10][11][12][13] and explained by the corresponding Landau phenomenological theories 1,2,12,14 , suggesting that we understand electrocaloric effects quite well. However, we think the situation is not fully satisfactory, for one main reason: The theoretical treatments in the literature are casespecific, and we lack a unified and simple picture revealing the similarities and differences among the various known effects. Also, we find that, while intuitively appealing, some frequent assumptions (i.e., that electric fields cause disorder in antiferroelectrics) are somewhat vague, and we miss a formalism that allows us to interpret the observed behaviors (and eventually think about new ones) in a more rigorous manner. Here, we provide such a formalism, and show how it helps us to rationalize all the (small-field) electrocaloric effects observed in the literature. Further, we show the results of atomistic simulations that allow us to explicitly evaluate intriguing predictions of the theory, such as the existence of competing contributions to the electrocaloric response.

Results
Formal considerations. Let us focus on how π α controls ΔT(E α ). For the sake of simplicity, we work with an approximate version of the adiabatic temperature change where the "(0)" superscript marks values evaluated at zero applied field and at the zero-field temperature T (0) . Thus, this approximate expression captures the electrocaloric effect as obtained when we neglect the field dependence of T and C E , as well as the variation of π α and C E caused by electrocaloric heating or cooling; these are common approximations in the electrocaloric literature. (Unless we work at low temperatures, we can safely assume |ΔT(E α )| ≪ T (0) , as the largest measured electrocaloric effects are typically below 20 K 15 . Also, we can expect relatively small variations in C E except in the close vicinity of phase transitions; more on this below.) Let us examine Eq. (4) analytically. We start by writing P as a Taylor series in E, where P (0) is the spontaneous polarization, χ (n) is the nth-order dielectric susceptibility tensor, and ϵ 0 is the vacuum permittivity. For convenience, we work with a Cartesian coordinate system with one axis parallel to α, the direction of the applied field. Hence, we have E β = δ βα E α , where δ βα is the Kronecker delta; we thus obtain the α component of the polarization as where only the α-diagonal tensor elements appear. This expression for P α is general and can be used to describe any phase of an insulating material, be it ferroelectric, antiferroelectric, paraelectric, or simply dielectric. Also, note that the dielectric tensors in the above equations describe the properties of a single crystal with a well-defined orientation with respect to our Cartesian setting. By taking the T-derivative of Eq. (6) at constant field, we obtain where the π (0) vector captures the T-dependence of the spontaneous polarization P (0) . For n ≥ 1, the π (n) tensors account for the field-induced pyroelectric effect, with where the pyroelectric tensor of nth-order in the field series depends on the susceptibility of (n − 1)th-order. Using Eq. (7), we resolve Eq. (4) and obtain ααα E 3 α þ ::: where we define ΔT (n) (E α ) as the nth-order contribution to the adiabatic temperature change. Also, here it is implicitly assumed that the pyroelectric coefficients are evaluated at T (0) . Let us now see how this expression allows us to understand all the known electrocaloric effects in ferroelectric and antiferroelectric compounds for small applied electric fields. Our conclusions are summarized in Table 1.
For the sake of concreteness, here we discuss ferroelectric and antiferroelectric materials with an isotropic (e.g., cubic) hightemperature paraelectric phase, as is the case of perovskite oxides (e.g., PbTiO 3 or PbZrO 3 ), noting that our arguments can be generalized.
As it is well-known, both ferroelectric and antiferroelectric phase transitions are characterized by a dielectric anomaly at zero field, i.e., a maximum of χ (0) at T C . To fix ideas, we can imagine that all the diagonal components of χ (0) follow a Curie-Weiss law approximately. (Our simulation results for PbTiO 3 -see Fig. 1are a representative case.) This maximum controls the Tdependence of χ (0) in a wide range around T C , yielding π ð1Þ αα > 0 for T < T C and π ð1Þ αα < 0 for T > T C . Hence, in particular, the paraelectric phase of all ferroelectrics and antiferroelectrics is characterized by π ð1Þ αα < 0; according to Eq. (9), this should result in ΔT(E α ) > 0 for all α directions, as it is indeed observed.
Interestingly, the situation is rather similar for an antiferroelectric state: we still have P (0) = π (0) = 0 and ΔTðE α Þ % ΔT ð2Þ ðE α Þ / π ð1Þ αα . However, now we have T < T C and, as mentioned above, π ð1Þ αα > 0. Hence, we expect ΔT(E α ) < 0, in agreement with the experimental observations. Note that in all the above cases the electrocaloric temperature change does not depend on the sign of the applied field, a feature that is expected from the symmetry of paraelectric and antiferroelectric states, and which we readily obtain from our formalism.
Suppose now that we are in a ferroelectric phase. Without loss of generality, we assume that P (0) is parallel to the z Cartesian direction, with P ð0Þ z > 0 and P ð0Þ x ¼ P ð0Þ y ¼ 0. In this case, the π (0) vector has one nonzero component (π ð0Þ z < 0) and two null ones (π ð0Þ x ¼ π ð0Þ y ¼ 0). Further, we have π ð1Þ αα > 0 for all α (from the Curie-Weiss law and the fact that T < T C ).
Let us imagine we apply a field E x along the x Cartesian direction, thus perpendicular to P (0) . This case is exactly analogous to the antiferroelectric state discussed above: according to Eq. (9), the response is controlled by π ð1Þ xx > 0. Hence, exactly as in the antiferroelectric state, we expect an inverse electrocaloric response with ΔT(E x ) < 0, an experimentally observed behavior that might seem surprising 6,16 , but is readily obtained and explained within our formalism.
Finally, suppose that we apply a field along z, the direction of the spontaneous polarization. Here, for the first time in this discussion, we have a nonzero linear contribution to the adiabatic temperature change, and we can write Concerning ΔT (2) (E z ), the situation is identical to the above cases for T < T C : it is dominated by π ð1Þ zz > 0, which yields ΔT (2) (E z ) < 0 regardless of the sign of the applied field E z .
In contrast, the linear contribution does depend on whether E z is parallel or antiparallel to spontaneous polarization. Further, ΔT (1) z , which is negative when P ð0Þ z > 0. Hence, we have: ΔT (1) (E z ) > 0 for a parallel field (E z > 0), and ΔT (1) (E z ) < 0 when the field goes against P z .
Thus, we have two qualitatively different cases. If the applied field goes against the polarization, ΔT (1) (E z ) and ΔT (2) (E z ) are both negative, and we have every reason to expect ΔT(E z ) < 0. However, for fields parallel to polarization, we have a competition between the linear and quadratic contributions to ΔT(E z ), and the net result is in principle undetermined. Interestingly, experimental studies of ferroelectric phases showwithout exception, as far as we know -that the temperature change is positive for fields parallel to the spontaneous polarization 1 , and negative for fields tending to reverse it 6,17,18 . This is in agreement with the expectations from our formalism, suggesting that, for the case of parallel fields, the linear effect (ΔT ð1Þ ðE z Þ / π ð0Þ z E z > 0) dominates over the quadratic one (ΔT ð2Þ ðE z Þ / π ð1Þ zz E 2 z < 0). Hence, as summarized in Table 1, we find that, in all the cases considered, the leading nonzero contribution to ΔT(E α ) agrees in sign with the adiabatic temperature change observed experimentally for relatively small applied fields. We should note that the formalism just introduced bears obvious similarities with previously proposed theories to discuss electrocaloric effects, e.g., in the context of antiferroelectrics 2,9 . The novelty here relies on the fact that our equations are general and can be applied to any material and phase (ferroelectric, antiferroelectric, paraelectric, or simply dielectric), revealing the way they are connected and evidencing the (somewhat trivial) origin of the so-called inverse effects. Our formalism also emphasizes that the basic electrocaloric response (i.e., the sign of the T-change) can be understood from simple universal arguments, not relying on specific atomistic or phenomenological models.
Numerical results. To gain further insight, and to evaluate the accuracy of low-order approximations to ΔT(E α ), we now compute explicitly the ΔT (n) (E α ) terms in Eq. (9) for prototype compound PbTiO 3 (PTO).
We simulate PTO using the so-called second-principles methods introduced in refs. 19 and 20 , i.e., atomistic model potentials with all parameters fitted from first principles and whose accuracy can be systematically improved (see "Methods"). More precisely, we use the model potential first introduced in ref. 19 , which has proven its accuracy in reproducing the basic ferroelectric behavior of the material as well as many subtle structural features, as e.g., related to its domain walls 21,22 . Let us stress that, in these simulations, all the degrees of freedom for the lattice (i.e., all atomic positions, all strains) are treated on equal footing; hence, our calculations include all contributions to the electrocaloric response, and there is in fact no easy way to differentiate them in the manner is often done in phenomenological Landau approaches (where it is natural to distinguish the "dipole" subsystem from the "phonon" bath 1 ). Finally, let us mention that the only noteworthy deficiency of our secondprinciples model pertains to the predicted T C , which is lower than the experimental one (510 K vs 760 K), but this is not critical for the present purposes.
We solve our PTO model as a function of temperature and applied electric field 19 by running Monte Carlo simulations, using periodically repeated supercells composed of 10 × 10 × 10 or 12 × 12 × 12 elemental perovskite units. (Larger cells are considered in the proximity of T C .) At a given T, we run 10,000 Monte Carlo sweeps for thermalization, followed by 75,000 to 100,000 sweeps to compute averages. We checked that these calculation conditions yield sufficiently accurate results.
The key quantities we monitor are the equilibrium polarization P, the dielectric susceptibility χ, the pyroelectric vector π, and the specific heat C E , which we compute using standard linearresponse formulas 23 . (The formulas used here are given in the Methods section; they closely resemble the ones employed in other studies of electrocaloric effects by Monte Carlo simulations, e.g., in refs. 12,14 .) We can thus compute the Taylor series for π (Eq. (7)) and for ΔT(E α ) (Eq. (9)). (In "Methods" and Supplementary Figs. 1, 2, and 3 we give extra details on the calculation of the ΔT (n) (E α ) contributions.) Our results are summarized in Fig. 1.
We obtain a ferroelectric phase transition at T C = 510 K (Fig. 1a), marked by a near divergence of the susceptibility (Fig. 1b); this is characteristic of a weakly first-order transformation and is in qualitative agreement with the experimental result 24 . Without loss of generality, we choose P ð0Þ ¼ ð0; 0; P ð0Þ z Þ, with P ð0Þ z > 0 below T C . Of note is the qualitative change of the phase transition as we apply an electric field: for the small applied field, the susceptibility peak reaches higher values, reflecting the fact that the transition becomes more continuous (for a discussion of this kind of effect, see, e.g., refs. 25 and 26 ); for larger fields, the transition becomes diffuse and the related anomalies in the susceptibility (Fig. 1b) and heat capacity (Fig. 1c) tend to disappear.
For clarity, Fig. 1b only shows the temperature and field dependence of one susceptibility component, χ zz . However, note that, by symmetry, the χ tensor only has diagonal components; further, for T > T C they are all equal (cubic phase), while for T < T C we have χ xx = χ yy > χ zz (tetragonal phase). (In the ferroelectric phase, the polar direction is electrically stiffer; this is a wellknown feature of ferroelectric perovskites, related to what is usually called "easy polarization rotation" 27 .) As can be seen in Supplementary Fig. 4, the three components have an approximate Curie-Weiss behavior, with a maximum at T C and a monotonic T-dependence on both sides of the transition point. Figure 1d shows our results for the low-order pyroelectric coefficients. π ð0Þ z is null above the transition point and negative below it, as expected. As for π (1) , we find the expected behavior as well: the tensor has nonzero diagonal components at all temperatures, nearly diverging at T C , and changing signs as the system goes through the transition. Note that our numerical  3 . We show the temperature and field dependence of P (a), χ (b), and specific heat C E (c). d shows the T-dependence of the π ð0Þ z , π ð1Þ xx , and π ð1Þ zz components of the pyroelectric tensors; they correspond to the case in which we have a positive spontaneous polarization along z, as discussed in the text. Also shown is ΔT(E z ) obtained for different electric fields (e), as well as the result for E z = 10 MV m −1 decomposed in the different ΔT (n) (E z ) contributions (f). a-c Different colors and symbols are used to represent different values of the applied electric field (black circles for zero field, etc.). d Different colors and symbols correspond to different tensor components (black squares for π ð0Þ z , etc.). e Different colors and symbols correspond to different electric fields (gray circles for 0.5 MV m −1 , etc.). f Different colors and symbols correspond to different orders in the applied field (black circles for the total temperature change, blue squares for the term that depends linearly on the field, etc.). e, f No results are shown for T ≳ T C whenever the applied E induces a transition. results for π ð0Þ z and π ð1Þ zz ratify the competition that occurs when a field is applied parallel to P ð0Þ z , as discussed above. Also, our results for π (1) reveal an essentially isotropic tensor at all temperatures, even in the tetragonal phase. (The tetragonal symmetry implies π ð1Þ xx ¼ π ð1Þ yy ≠ π ð1Þ zz ; yet, in Fig. 1d, a sizeable difference between the xx and zz components is found only for T ≲ T C .) Figure 1e shows the computed ΔT(E z ) for various positive fields, with E z up to 10 MV m −1 , as a function of temperature. The results are obtained by evaluating the Taylor series in Eq. (9) up to 5th order, which is enough to get a converged T-change except in the immediate vicinity of the phase transition. More specifically, Fig. 1f shows that, even for a large field of 10 MV m −1 , the electrocaloric effect at T < T C is dominated by the leading contribution ΔT (1) (E z ), higher-orders becoming significant only very close to T C (see the result at 490 K). In fact, as shown in Supplementary Fig. 2, we find that the leading contribution to the electrocaloric response-i.e., ΔT (1) (E z ) below T C and ΔT (2) (E z ) above T C , respectively-is usually a very good approximation of the total effect.
Let us stress that in Fig. 1e, f, we restrict ourselves to fields that are small enough so that no first-order phase transition is induced. This is why we do not show any data very close to T C , and why the fields considered for T > T C are relatively tiny. (Close to T C , the paraelectric phase is easily transformed into the polar one. See representative results in Supplementary Fig. 3.) This is consistent with our perturbative approach to compute ΔT(E α ), which is designed to describe the properties of the continuously deformed zero-field state. To treat a first-order transition, one could split the integral for ΔT(E α ) in low-and high-field parts, using different perturbative P(E) expansions for each of them. In addition, one should account for the latent heat associated with the discontinuous transformation [28][29][30] . The information to tackle such situations is in principle available from our Monte Carlo calculations. (We can compute latent heat from the thermalaveraged internal energies 31 ). However, we should note that, for treating first-order transitions, direct non-perturbative computational approaches based on microcanonical molecular dynamics 8 or constrained Monte Carlo simulations 32 are better suited.
The values obtained for ΔT(E z ) (e.g., a maximum of 0.25 K when T C is approached from below for a field of 0.5 MV m −1 ) are comparable with electrocaloric effects measured for PTO; for example, ref. 33 reports a temperature change of 0.1 K for T ≲ T C and a field of 0.15 MV m −1 , and a maximum T-change of 1.9 K at T C . (We also obtain ΔS(E z ) ≈ − 2500 J K −1 m −3 for a field of 0.5 MV m −1 at T ≲ T C , while the ref. 33 reports a maximum entropy change of about − 16,500 J K −1 m −3 for a field of 0.15 MV m −1 applied exactly at T C .) Our results are also consistent with other theoretical estimates of the electrocaloric effect for PTO 34 . Finally, in Supplementary Fig. 5, we evaluate the approximations made in Eq. (9)-i.e., the use of the zero-field values for T and C E , so they can be taken out of the integral-and find that their impact is negligible (even for large fields) except very close to T C .

Discussion
In view of the above, let us comment on the usual interpretation of electrocaloric effects in terms of field-induced order or disorder. For illustrative purposes, it is convenient to pay attention to the isothermal entropy change, which, from Eqs. (2) and (7), can be written as ðE α Þ þ ::: The last line suggests that we can think of ΔS(E z ) as depending linearly on E z , the corresponding proportionality constant having spontaneous ($ π ð0Þ α ) and field-induced ($ π ð1Þ αα E α ) pyroelectric contributions.
Let us consider a ferroelectric state with P ð0Þ z > 0, and imagine we apply an electric field E z > 0. For simplicity (and without loss of generality), let us also assume that the dielectric response is dominated by the linear effect χ ð0Þ zz , higher-order terms being negligible to a good approximation. In such a situation, it is physically sound to assume that the applied field creates order, as it contributes to further align the local electric dipoles in the ferroelectric state and results in a larger order parameter P z . This field-induced ordering is captured by χ ð0Þ zz > 0, and we know that the effect gets stronger as we approach T ≲ T C . As a result of this ordering, we expect the entropy to decrease in an isothermal process (ΔS(E z ) < 0), and the temperature to rise (ΔT(E z ) > 0) if the process is adiabatic. These are clear expectations that one would hardly question.
Let us inspect how these expected variations of S and T come about in our formalism. Following the simplification mentioned above (linear dielectric response), we can write ΔS(E z ) ≈ ΔS (1) (E z ) + ΔS (2) (E z ). (The discussion for ΔT(E z ) is analogous.) The first term (ΔS ð1Þ ðE z Þ ¼ π ð0Þ z E z ) results in a reduction of the entropy, as we have π ð0Þ z <0 for the spontaneous pyroelectric effect. This seems consistent with the ordering argument given above. However, it is important to stress that this term does not contain any information about the dielectric response to the applied field, or about the order the field may create. (There is no generally expected thermodynamic relationship between ∂P z /∂E z and ∂P z /∂T.) Instead, the entropy change is fully determined by the T-dependence of the spontaneous polarization: P ð0Þ z > 0 decreases as we approach T C , a behavior that is normal.
The response to an applied field does control the quadratic contribution ΔS (2) (E z ), where the spontaneous pyroelectric effect in ΔS (0) (E z ) is replaced by the field-induced one ($ π ð1Þ zz E z ). As mentioned above, it is clear that the applied field strengthens the dipole order of the ferroelectric state, as quantified by its order parameter: we have P ð0Þ z > 0 and a field-induced polarization change ϵ 0 χ ð0Þ zz E z > 0. However, the corresponding entropy change is positive, since π ð1Þ zz / ∂χ ð0Þ zz =∂T > 0 for T < T C . Thus, the same physical mechanism (captured by χ ð0Þ zz ) results in field-induced order and a positive contribution to the entropy. Indeed, in what concerns the isothermal entropy change, what matters is not the ordering of the dipoles at a given T. Instead, we have to pay attention to how the field-induced effect changes with temperature. In a ferroelectric state, the ordering caused by the field grows as we heat up towards T C , that is, we have a positive field-induced pyroelectric effect. This somewhat anomalous behavior (the higher the temperature, the greater theinducedorder) is opposite to the spontaneous pyroelectric effect, and is the one causing ΔS (2) (E z ) > 0 (and ΔT (2) Hence, the quadratic contribution to the electrocaloric effect in a ferroelectric state, dominated by the linear susceptibility χ ð0Þ zz , seems paradoxical: the field creates order (as quantified by the order parameter) but the entropy increases. The key to the puzzle is that polarization (order) created at a given T is not what controls the entropy; its T-derivative is. This point is hardly new, as it is quite clear from the well-known equations governing electrocaloric effects (Eqs. (1) and (2)); yet, it becomes particularly apparent in our perturbative approach.
Keeping the above in mind, it is easy to interpret the electrocaloric effects in paraelectric and antiferroelectric states, which are dominated by the quadratic contribution. In the paraelectric case, the field-induced order decreases as T increases: at T > T C we have π ð2Þ It is tempting to interpret this result by focusing on the response at a given T (the field creates order, hence the entropy should get reduced), but that would not be correct. Instead, the focus should be on the field-induced pyroelectric effect which, in this case, behaves in the normal way, i.e., we have a weakerinducedorder at a higher temperature.
Finally, the inverse electrocaloric behavior of antiferroelectric states is also readily explained within this picture. In that case, we are at T < T C , and have ΔSðE z Þ % ΔS ð2Þ ðE z Þ / ðπ ð2Þ zz E z ÞE z > 0 (and ΔT(E z ) ≈ ΔT (2) (E z ) < 0). According to the above discussion, this inverse behavior stems from the fact that the applied field creates polarization more efficiently as T grows toward T C . Consequently, this disproves the frequent interpretation that the inverse electrocaloric response of antiferroelectric states is caused by fieldinduced disorder at a given T. As a matter of fact, the above analysis shows that arguments focusing on the field-induced (dis) order at a given Tand overlooking the T-derivativeare incorrect.
To conclude this part, let us comment on a ferroelectric type absent in our discussion: relaxors. As mentioned above (Eq. (1)), ergodicity breaks down in these materials, and the use of the Maxwell relation our formalism relies on is questionable. Nevertheless, some experimental 35 and theoretical 14 reports suggest that, in fact, Eq. (1) approximately captures the electrocaloric response of these systems at a quantitative level; further, it seems to yield qualitatively correct results. Hence, we think that, in the future, it may be worth considering whether our formalism, and the kind of rules summarized in Table 1, might be useful in investigations of relaxors.
In summary, we have introduced a perturbative approach to the electrocaloric effect. This formalism can be applied in quantitative simulation studies, in a straightforward way as long as field-induced (first-order) phase transitions are not present. More specifically, our simulations for ferroelectric PbTiO 3 show that, except in the vicinity of the phase transition, a low-order approximation captures the electrocaloric response with great quantitative accuracy. Most importantly, our formalism unifies and clarifies the physical interpretation of all the small-field electrocaloric effects (normal and inverse) observed in ferroelectric and antiferroelectric materials.

Methods
Second-principles simulations. For the second-principles simulations, we used the software called SCALE-UP, which implements the effective-potential approach described in detail in refs. 19,20 . The procedure to obtain this package is described by the developers in 36 . The reader may be interested to know that the same second-principles scheme to simulate lattice-dynamical properties, introduced by Wojdeł et al. in ref. 19 , has also been implemented in the package multibinit 37 .
As for the visualization and analysis of the simulation data (e.g., for taking numerical derivatives), we used standard free open-source tools that can be readily obtained to run under linux; in particular, much of the work was done using the QtiPlot package 38 .
Physical properties from the Monte Carlo simulations. Here, we explain how to compute all the quantities that are relevant in this study-i.e., χ, π, and C E -from Monte Carlo simulations, at given conditions of temperature and applied electric field. Let us begin with the electric susceptibility, where α and α 0 label spatial directions in Cartesian coordinates, and 〈. . . 〉 indicates the thermal average (equilibrium value) as obtained from our Monte Carlo simulations. This thermal average can be written as where Z is the partition function given by The sums above run over all the microscopic states of the material, labeled by i; U i , and P i are the energy and polarization of state i, respectively; β ¼ ðk B TÞ À1 , where k B is Boltzmann's constant and T is the temperature; finally, E is the applied electric field and V is the volume of the simulated system (as usually done 23 , we take the equilibrium volume and assume that volume fluctuations are negligible). Computing analytically the field derivative, we obtain Analogously, we can derive the pyroelectric tensor, π. We have is the equilibrium energy. Finally, for the specific heat C E we have where N is the number of unit cells in the simulation supercell. Here, the first term is related to the potential energy, while the second accounts for the (trivial) kinetic contribution (note we have 3 × 5 degress of freedom in the five-atom perovskite cell). Calculating the derivative, we finally get Computing the adiabatic temperature change. To compute ΔT (n) (E α ) to arbitrary order (Eq. (9)), we proceed as follows.
ΔT (1) (E α ) is simply proportional to the spontaneous pyroelectric efffect (π ð0Þ =C ð0Þ E , Eq. (9)), which we obtain directly from zero field Monte Carlo simulations, by evaluating the thermal averages as described above.
The higher-order contributions are controlled by field-induced pyroelectricity and can be deduced from the field-dependent susceptibility (see Eqs. (8) and (9)). More precisely, for T < T C , we fit χ zz = χ zz (E z ) up to third order in the field E z , which allows us to compute terms up to ΔT (5) (E z ) in Eq. (9). For T > T C , we fit χ zz = χ zz (E z ) up to fourth order in the field E z , which allows us to compute terms up to ΔT (6) (E z ). Supplementary Fig. 1 summarizes the outcomes of this fitting exercise.
Let us note that the slightly different fitting choices, above and below T C , were adopted in view of the obtained results; in both cases, we tried to push the fit to order as high as possible, while avoiding overfitting. It is also important to keep in mind that, above T C , the χ ðnÞ z:::z terms with odd n are zero by symmetry (i.e., χ ð1Þ zzz ¼ χ ð3Þ zzzzz ¼ ::: ¼ 0); in turn, this implies that the odd-order terms in the ΔT(E z ) expansion are also null above T C .
Finally, we compute numerically (by finite differences) the temperature derivatives of the susceptibility coefficients, and thus obtain the pyroelectric contants π ðnÞ z:::z corresponding to Eq. (7). The ΔT (n) (E z ) are then trivially calculated.
Representative results for the adiabatic temperature change are shown in Supplementary Fig. 2. It should be stressed that, because our theory is a perturbative one, behaviors associated to discontinuous phase transitions are left out of the fit. A clear example of this occurs at temperatures immediately above T C , where small fields are able to produce a transformation, as shown in Supplementary Fig. 3. Hence, the computed dielectric constants at 510 K are not fitted at all; and, above that temperature, the fitting range is restricted to small electric fields (see Supplementary Fig. 1).

Data availability
All the relevant data are available from the authors upon reasonable request.

Code availability
The second-principles calculations are carried out using the package SCALE-UP which is proprietary software. Some of the visualizations were done using the free package QtiPlot 38 .