Electrocaloric effects in multiferroics
Zhijun Jiang, Bin Xu, Sergey Prosandeev, Yousra Nahas, Sergei Prokhorenko, Jorge ?ñiguez, L. Bellaiche
aa r X i v : . [ c ond - m a t . m t r l - s c i ] F e b Electrocaloric effects in multiferroics
Zhijun Jiang,
1, 2, 3
Bin Xu, Sergey Prosandeev, Yousra Nahas, Sergei Prokhorenko, Jorge ´I˜niguez, and L. Bellaiche MOE Key Laboratory for Nonequilibrium Synthesis and Modulation of Condensed Matter,School of Physics, Xi’an Jiaotong University, Xi’an 710049, China Physics Department and Institute for Nanoscience and Engineering,University of Arkansas, Fayetteville, Arkansas 72701, USA Key Laboratory of Computational Physical Sciences (Ministry of Education),State Key Laboratory of Surface Physics, and Department of Physics, Fudan University, Shanghai 200433, China Jiangsu Key Laboratory of Thin Films, School of Physical Science and Technology, Soochow University, Suzhou 215006, China Materials Research and Technology Department, Luxembourg Institute of Science and Technology,5 Avenue des Hauts-Fourneaux, L-4362, Esch/Alzette, Luxembourg Physics and Materials Science Research Unit, University of Luxembourg, 41 Rue du Brill, L-4422 Belvaux, Luxembourg
An atomistic effective Hamiltonian is used to compute electrocaloric (EC) effects in rare-earth substitutedBiFeO multiferroics. A phenomenological model is then developed to interpret these computations, with thismodel indicating that the EC coefficient is the sum of two terms, that involve electric quantities (polarization,dielectric response), the antiferromagnetic order parameter, and the coupling between polarization and antifer-romagnetic order. The first one depends on the polarization and dielectric susceptibility, has the analytical formpreviously demonstrated for ferroelectrics, and is thus enhanced at the ferroelectric Curie temperature. Thesecond one explicitly involves the dielectric response, the magnetic order parameter and a specific magneto-electric coupling, and generates a peak of the EC response at the N´eel temperature. These atomistic results andphenomenological model may be put in use to optimize EC coefficients. The electrocaloric (EC) effect is a phenomenon by whicha material exhibits a reversible temperature change under theapplication/removal of an electric field [1–5]. It is attractingattention due to its potential to be an efficient solid-state re-frigeration technology (see, e.g., Refs. [6–19] and referencestherein).Furthermore, multicaloric effects that are driven simultane-ously by more than one type of external physical handle, suchas electric and/or magnetic fields, mechanical stress and pres-sure [20–24], are also promising to enhance change in tem-perature [23, 24].Recently, multiferroics, which are materials that possesscoupled long-range-ordered electric and magnetic degrees offreedom [25–31], have also been mentioned as possible sys-tems to enhance the EC effects by taking advantage of suchcoupling [12, 20, 21, 32–34]. The pioneering work of Ref.[33] started from a phenomenological Landau-type equationfor which coefficients were determined from first principles toinvestigate how magnetoelectric coupling modifies the EC co-efficient. The main result was that EC effects are significantlyenhanced (by about ) thanks to magnetoelectric couplingin the case that the ferroelectric and magnetic critical temper-atures coincide. However, one has to be careful when usinga Landau-type approach because fluctuations, which can beimportant for responses, are not treated explicitly and maybe underestimated. That is why atomistic approaches incor-porating couplings between electric dipoles and spins can beuseful to also study EC effects in multiferroics, as the authorsof Ref. [33] indicated. More importantly, it is presently un-clear how to understand EC coefficients in multiferroics. Forinstance, can these coefficients be considered as composed oftwo terms, with one corresponding to that occurring in normalferroelectrics and the second one related to the coupling be-tween spins and electric dipoles? If yes, what are the precise quantities involved in the second term? Are they only mag-netoelectric, or rather also involve electric and/or magneticproperties? Answering such questions will help in designingsystems with large EC response.The aim of this Letter is to resolve all these issues by (1)conducting atomistic-based simulations; (2) developing a sim-ple model that can reproduce these simulations; and (3) usingsuch simulations and model to gain a deep microscopic in-sight. We demonstrate that the EC coefficient of multiferroicscan be thought as having two parts, each associated with dif-ferent physical quantities.Here, we adopt the effective Hamiltonian ( H eff ) approachdeveloped in Ref. [35] to study disordered Bi − x Nd x FeO (BNFO) alloys. H eff parameters are provided in the Sup-plemental Material (SM) [36]. This H eff successfully repro-duced the temperature- versus -compositional phase diagramof BNFO. It predicts a R c ground state for small Nd com-positions and a P nma phase for larger concentrations, withintermediate complex states in-between. Moreover, within thecompositional range for which the R c phase is the groundstate, the ferroelectric Curie temperature T C was numericallyfound to significantly decrease with the Nd composition whilethe T N N´eel temperature is mostly independent of concentra-tion, which also agrees with measurements [51–53]. The totalinternal energy of this H eff can be expressed as a sum of twomain terms: E tot = E BFO ( { u i } , { η H } , { η I } , { ω i } , { m i } )+ E alloy ( { u i } , { ω i } , { m i } , { η loc } ) , (1)where E BFO is the H eff of pure BiFeO [39–42] and E alloy characterizes the effect of substituting Bi by Nd ions. The H eff of BNFO contains four types of degrees of freedom: (i)the local soft mode { u i } centered on the A site of Bi or Ndions in the -atom unit cell i (which is proportional to the localelectric dipole moment of that cell [43, 44]); (ii) the strain ten-sor gathering homogeneous { η H } and inhomogeneous { η I } contributions [43, 44]; (iii) the pseudovectors { ω i } that repre-sent the oxygen octahedral tiltings [45]; and (iv) the magneticmoments { m i } centered on Fe ions [54].We employ this H eff within Monte Carlo (MC) simulationson × × supercells (containing atoms) with pe-riodic boundary conditions and inside which Bi and Nd ionsare randomly distributed over the A sublattice.
20 000
MCsweeps are used for equilibration and an additional
20 000
MC sweeps are employed to compute statistical averages atfinite temperature, to obtain converged results. We also aver-age our results over random Bi/Nd distributions, to mimicwell disordered BNFO solid solutions.Regarding the linear EC coefficient, α γ , it is the derivativeof the temperature with respect to electric field at constant en-tropy. It can be obtained from MC simulations by taking ad-vantage of the cumulant formula [16, 17, 55]: α γ = − Z ∗ a lat T ( h u γ E tot i − h u γ i h E tot i (cid:10) E tot2 (cid:11) − h E tot i + k B T ) N ) , (2)where Z ∗ is the Born effective charge associated with the localmode, a lat represents the five-atom lattice constant, T is thetemperature, u γ is the γ -component of the supercell averageof the local mode with γ = x , y , or z (note that the x , y , and z axis are chosen along the pseudocubic [100] , [010] and [001] directions, respectively), E tot is the total internal energy givenby the H eff , k B is the Boltzmann constant, N is the numberof sites in the supercell, and h i defines the average over theMC sweeps at a given temperature [56]. In the following,we will denote α the quantity defined by α x + α y + α z √ . Suchdefinition corresponds to the EC response for an electric fieldapplied along [111] , which is the maximal response within a R c state.Figure 1 shows the EC coefficient as a function of tem-perature for four different Nd compositions in disorderedBi − x Nd x FeO . The results of Fig. 1 are obtained by start-ing from K adopting a R c phase and then progressivelyheating up the BNFO solid solutions up to the composition-dependent Curie temperature, T C (for all investigated temper-atures displayed in Fig. 1, the disordered Bi − x Nd x FeO al-loys possess the R c phase from K and up to T C ). This R c state is characterized by a polarization lying along [111] and oxygen octahedra tilting in an antiphase fashion about thispolarization’s direction. These solid solutions also exhibit aG-type antiferromagnetic-to-paramagnetic transition at a N´eeltemperature, T N , which is mostly independent on the compo-sition and equal to ≃ K [35]. The SM [36] provides somefinite-temperature properties above T C .Let us first focus on Fig. 1(a) that corresponds to a concen-tration of Nd equal to . The calculated T C ≃ K and T N ≃ K of Bi . Nd . FeO are in rather good agree-ment with the measurements of T C ≃ K and T N ≃ K (a) (b)(c) (d) Bi Nd FeO Bi Nd FeO Bi Nd FeO Bi Nd FeO FIG. 1. Electrocaloric coefficient, α , as a function of the tempera-ture for different compositions in disordered Bi − x Nd x FeO alloys:(a) Bi . Nd . FeO ; (b) Bi . Nd . FeO ; (c) Bi . Nd . FeO ;and (d) Bi . Nd . FeO . The solid green lines represent thefit of the MC results by the second line of Eq. (11), i.e., α = T a ′ ( T ) C ph P s ε χ + T b ′ ( T ) C ph L s ∂L s ∂P s (cid:12)(cid:12)(cid:12) T ε χ , where a ′ ( T ) = A + A T ( A and A being fitting constants), and C ph and b ′ ( T ) are also fit-ting parameters. The solid blue lines display the fit of the MC re-sults by its first contribution, T a ′ ( T ) C ph P s ε χ . The solid brown linescorrespond to the fit of the MC results by its second contribution, T b ′ ( T ) C ph L s ∂L s ∂P s (cid:12)(cid:12)(cid:12) T ε χ (see text). [52, 53]. For any investigated temperature, α basically mono-tonically increases when the system is heated up to the N´eeltemperature. It then adopts a small peak around T N , whichis found to originate from the coupling between polarizationand magnetism – we verify this by running H eff simulationsin which the coupling between local models and magneticmoments is turned off. The EC coefficient then significantlystrengthens when increasing the temperature from the end ofthis ≃ T N -centered peak and up to T C . Our predicted bigvalue of α around T C is of the order of ≃ . × − K m/V.It is thus large and close to the experimental data of . × − K m/V at T ≃ K in PbZr . Ti . O films [7] (the largestobserved α is equal to × − K m/V and has been foundin a BaTiO single crystal, see Ref. [57]) [58]. Note that H eff techniques have been demonstrated in Refs. [16, 17] to accu-rately reproduce the EC coefficients of ferroelectrics and re-laxor ferroelectrics, such as those reported in BaTiO [57, 59]and Pb(Mg,Nb)O [49].Let us now concentrate on other compositions in disorderedBi − x Nd x FeO alloys. Figures 1(b)-1(d) show the depen-dence of the EC coefficient when the Nd composition is equalto x =0 . , . and . , respectively. The Curie temper-ature T C noticeably decreases when increasing the Nd com-position, as consistent with observations and computations[35, 51–53]. Consequently, the two critical temperatures coin-cide, i.e. T C = T N , for a Nd concentration of . . Figures1(b)-1(d) especially reveals that α at the N´eel temperature isenhanced when the Nd composition increases, but it becomesmore difficult to see its associated peak.To understand the results in Fig. 1, we use a Landau free-energy potential F ( P, L, E , T ) in which we substitute polar-ization P and G-type antiferromagnetic (AFM) moment L bytheir equilibrium values P s and L s found from minimizationof free energy: ∂F∂P (cid:12)(cid:12) P = P s , E ,T = 0 and ∂F∂L (cid:12)(cid:12) L = L s , E ,T = 0 .The minimized free energy F s ( E , T ) = F ( P s , L s , E , T ) hasthe form: F s ( E , T ) = 12 a ( T ) P s ( E , T ) + 14 βP s ( E , T ) − E P s ( E , T ) + 12 b ( T ) L s ( E , T )+ 14 κL s ( E , T ) + 12 cL s ( E , T ) P s ( E , T ) , (3)where E is the electric field.Such equation implies that the polarization implicitly de-pends on magnetism, because of the cL s ( E , T ) P s ( E , T ) term. This equation is similar to the one used in Ref. [33]. Theentropy described by this free energy F s ( E , T ) , composed ofdipoles and spins, can then be obtained as S F ( E , T ) = − dF s dT (cid:12)(cid:12)(cid:12)(cid:12) E = − a ′ ( T )2 P s ( E , T ) − b ′ ( T )2 L s ( E , T ) , (4)where a ′ = da/dT and b ′ = db/dT . Note that, here we tookinto account that P s and L s are found from minimization ofthe free energy.In the case of a magnetic phase transition and presence ofpolarization, we can consider two parts of the total entropy S ( E , T ) : A first one due to electric dipoles and spins (the ac-tive part treated by the Landau potential above, with entropy S F ( E , T ) ) and a second one due to the rest of the lattice (theinert part that can be considered to be a trivial collection ofharmonic phonons, with entropy S ph ( T ) ) [60, 61]. For anadiabatic process, we have: ∆ S ( E , T ) = ∆ S F ( E , T ) + ∆ S ph ( T ) = 0 . (5)Let C ph denote the heat capacity associated with the back-ground lattice modes. Then the change of lattice entropy froman initial state (0 , T ) to the final state ( E , T ) is given by: ∆ S ph = Z TT C ph T dT ∼ = C ph ln (cid:18) TT (cid:19) . (6)Consequently, combining Eqs. (5) and (6) leads to C ph ln (cid:18) TT (cid:19) = − ∆ S F = 12 a ′ ( P s − P ) + 12 b ′ ( L s − L ) . (7)Here P s = P s ( E , T ) , P = P s (0 , T ) , L s = L s ( E , T ) , L = L s (0 , T ) , where T is the initial temperature and T = T + ∆ T is the final temperature ( ∆ T represents the temperaturechange). Solving this equation with respect to T /T yields: ( T + ∆ T ) /T = e [ a ′ ( P s − P )+ b ′ ( L s − L ) ] / C ph . (8)For small ∆ T : ∆ T = T (cid:2) a ′ ( P s − P ) + b ′ ( L s − L ) (cid:3) C ph . (9)One can then derive the following expression for α [10, 16]: α = ∂ ∆ T∂ E (cid:12)(cid:12)(cid:12)(cid:12) S ≈ T a ′ ( T )2 C ph ∂P s ∂ E (cid:12)(cid:12)(cid:12)(cid:12) T + T b ′ ( T )2 C ph ∂L s ∂ E (cid:12)(cid:12)(cid:12)(cid:12) T . (10)Here we assumed that, since the adiabatic temperaturechange is small as compared to the temperature, the constant- S derivatives can be evaluated at a constant T = T . One canwrite: α = T a ′ ( T ) C ph P s ε χ + T b ′ ( T )2 C ph ∂L s ∂P s (cid:12)(cid:12)(cid:12)(cid:12) T ∂P s ∂ E (cid:12)(cid:12)(cid:12)(cid:12) T = T a ′ ( T ) C ph P s ε χ + T b ′ ( T ) C ph L s ∂L s ∂P s (cid:12)(cid:12)(cid:12)(cid:12) T ε χ , (11)where ε is the vacuum permittivity and χ is the dielectricsusceptibility. Finally, let us note that one could try to ap-proximate C ph by adding a k B contribution for each degreeof freedom belonging to the trivial – harmonic – part of thesystem. However, it is not obvious how to count the exactnumber of active and inactive variables in the framework ofa Landau theory; we thus treat C ph as an adjustable parame-ter. Note that we did not fit C ph alone but rather the ratio of a ′ ( T ) / C ph and b ′ ( T ) / C ph .As shown by the green lines of Fig. 1, the second line ofEq. (11) fits well the MC data, when (1) using the P s , χ , L s and ∂L s ∂P s [62] obtained by our Monte-Carlo simulations (thesefour quantities are shown in Fig. 2 for the case of a Ndcomposition); and (2) assuming that C ph and b ′ ( T ) are fittingconstants, while a ′ ( T ) = A + A T with A and A are fittingparameters [63]. Since its validity is confirmed by Fig. 1, thesecond line of Eq. (11) can now be used to gain an insight [36]into the results of Fig. 1, via the decomposition of α into itstwo terms – that are T a ′ ( T ) C ph P s ε χ and T b ′ ( T ) C ph L s ∂L s ∂P s (cid:12)(cid:12)(cid:12) T ε χ .The first contribution has precisely the analytical form of theEC coefficient for non-magnetic systems , see Refs. [16, 17]. Itis shown by blue lines in Fig. 1, and is the one that contributesthe most to the total α for any composition. Its increases withtemperature and is driven by the corresponding increase in di-electric susceptibility, however moderated by the concomitantdecrease in polarization [see Figs. 2(b) and 2(a)]. This firstcontribution implicitly depends on magnetism because of thecoupling between polarization and antiferromagnetism, as ev-idenced in the change of behavior of the polarization and inthe occurrence of a plateau in the dielectric response near T N (such behavior of χ has been reported in other multiferroics[39, 64]). The second contribution of Eq. (11) is depictedin brown lines in Fig. 1, and is basically independent on theinvestigated composition for any temperature. As evidencedin Fig. 1, it is the one responsible for the small peak of α found near the N´eel temperature. This small peak becomesmore difficult to be seen in the total EC coefficient (shown ingreen) when the Nd composition increases simply because thefirst contribution provides much larger values than the secondcontribution. Figures 2(c) and 2(d) also reveal that this smallpeak originates from the activation and then sharp increaseof the magnitude of ∂L s ∂P s near T N . This derivative for tem-peratures far away below T N is then basically a constant thatcharacterizes intrinsic magnetoelectric coupling – which is re-lated to the c constant of Eq. (3). The second term of Eq. (11)tells us that the EC coefficient of a multiferroic can be opti-mized even at temperatures far away T N in systems possessingstrong coupling between polarization and magnetic ordering.Ba(Sr,Ba)MnO films may thus be a system of choice to in-vestigate electrocaloric effects due to its strong magnetoelec-tric coupling [65–67].The now-elucidated effect of ∂L s ∂P s on α near T N can befurther used to address the finite-size effects in our compu-tations of the EC coefficient. It is known that such size effectbroadens the magnetic transition when decreasing the super-cell size (see the SM [36]) [47, 48], and we also checked thatthe magnitude of the second contribution of α around T N in-creases when increasing such size. It will thus be more real-istic, regarding what to expect in experiments, to rather adopta L s = A | T N − T | β power law (see Refs. [47, 50]) near theN´eel temperature, where A and β are coefficients. Conse-quently, we (1) chose to replace, around T N , the MC data for L s by the result given by such power law with β equal to . (mean-field value); (2) continue to still use the MC data for L s for temperatures far away (below) the N´eel temperature; and(3) extract A such by imposing that this power law of item (1)matches the MC data of item (2). Using the new resulting ∂L s ∂P s along with all the previous other quantities in Eq. (11) (includ-ing the temperature behavior of the polarization) provides thedata given in Fig. 3 for the second contribution but also to-tal EC coefficient in disordered Bi . Nd . FeO alloys. Theaforementioned change of L s ’s behavior, that is a more abruptchange near T N , leads to a narrower and stronger peak of α close to the N´eel temperature. The second contribution nowamounts for of the total EC coefficient near the magnetictransition. Such latter result is in-line with the phenomenolog-ical theory of Edstr¨om et al. [33] predicting that the magneticcontribution can reach approximately of the electric con-tribution at the magnetic transition, and thus enhance the ECeffect, in epitaxial multiferroic SrMnO systems under a ten-sile strain of . – for which T N = T C . Our study explainswhy it is the case thanks to Eq. (11) that not only reproducesatomistic results but also and especially provides an insightinto the microscopic origins of the EC effects in a multifer-roic. We also used a larger supercell and such power law of L s with different β , and found that our qualitative results are (a) (b)(c) (d) Bi Nd FeO FIG. 2. Temperature dependence of some properties in disorderedBi . Nd . FeO alloys, as obtained from our MC simulations: (a)the macroscopic polarization P s ; (b) the average between the threediagonal elements of the dielectric susceptibility; (c) the AFM vec-tor; and (d) the derivative dL s / dP s .FIG. 3. Same as Fig. 1 (a) but now using a different ∂L s ∂P s (see text)in the second line of Eq. (11). still valid for any reasonable choice of β (see Fig. S3 of theSM [36]). Note that the peak of Fig. 1(a) at the N´eel tem-perature is significantly less pronounced than in Ref. [33] fortwo possible reasons. The first one is that such peak dependson the size of the simulation supercell (see the SM [36]) andthe second one is that the magnetoelectric coupling is weakerin BiFeO [39] than in SrMnO [33]. Fluctuations within the H eff are also discussed in the SM [36].In summary, an atomistic effective Hamiltonian scheme hasbeen used to compute finite-temperature electrocaloric coeffi-cients in the rare-earth substituted BiFeO multiferroic. Theresults are then interpreted via the development of a modelthat reproduces these computational data. EC coefficients canbe decomposed in two main terms. The first term takes itslargest value at the Curie temperature and explicitly dependson the polarization and dielectric susceptibility, that are bothimplicit functions of magnetic ordering and strength becauseof magnetoelectric couplings. The second term adopts a peaknear the N´eel temperature and is proportional to the antifer-romagnetic vector, the polarization derivative of the antiferro-magnetic vector and the dielectric susceptibility. Such find-ings therefore suggest an original way to induce large ECcoefficients by simultaneous optimization of electric, mag-netic and magnetoelectric properties at a selected temperaturebelow the N´eel temperature: (1) the dielectric susceptibilityshould be large; (2) the antiferromagnetic vector should bestrong; and (3) the magnetoelectric coupling ∂L s ∂P s should belarge [68]. Our results and phenomenology should be validfor all magnetoelectric multiferroics, at the exception of thosefor which a magnetic Dzyaloshinskii-Moriya interaction in-volving the polarization (e.g., the spin-current model) is im-portant. We hope that the present article deepens the fields ofmultiferroics and important subtle cross-coupling propertiessuch as electrocaloric effects.This work is supported by the National Natural Sci-ence Foundation of China (Grants No. 11804138 and No.11825403), Shandong Provincial Natural Science Foun-dation (Grant No. ZR2019QA008), China PostdoctoralScience Foundation (Grants No. 2020T130120 and No.2018M641905), “Young Talent Support Plan” of Xi’an Jiao-tong University, Postdoctoral International Exchange Pro-gram of Academic Exchange Project, and Shanghai Post-doctoral Excellence Program. B. X. acknowledges finan-cial support from National Natural Science Foundation ofChina (Grant No. 12074277), the startup fund from Soo-chow University and support from Priority Academic ProgramDevelopment (PAPD) of Jiangsu Higher Education Institu-tions. S. Prosandeev is supported by ONR Grant N00014-17-1-2818. Y. N., S. Prokhorenko and L. B. thank theDARPA Grants No. HR0011727183-D18AP00010 (TEEprogramme) and No. HR0011-15-2-0038 (MATRIX pro-gram). J. ´I. acknowledges funding from the Luxembourg Na-tional Research Fund through the CORE program (Grant No.FNR/C18/MS/12705883 REFOX, J. ´I.). [1] M. E. Lines and A. M. Glass, Principles and Applications ofFerroelectrics and Related Materials (Oxford University Press,New York, 1977).[2] J. F. Scott, Science , 954 (2007).[3] J. F. Scott, Annu. Rev. Mater. Sci. , 229 (2011).[4] T. Correia and Q. Zhang, Electrocaloric Materials (Springer,Berlin, 2014).[5] Z. Kutnjak, B. Roˇziˇc, and R. Pirc,
Electrocaloric Effect : The-ory , Measurements , and Applications (Wiley Encyclopedia ofElectrical and Electronics Engineering, 2015).[6] K. Uchino, Ferroelectric Devices (Marcel Dekker, New York,2000).[7] A. S. Mischenko, Q. Zhang, J. F. Scott, R. W. Whatmore, andN. D. Mathur, Science , 1270 (2006).[8] S. Prosandeev, I. Ponomareva, and L. Bellaiche, Phys. Rev. B , 052103 (2008).[9] I. Ponomareva and S. Lisenkov, Phys. Rev. Lett. , 167604 (2012).[10] M. C. Rose and R. E. Cohen, Phys. Rev. Lett. , 187604(2012).[11] E. Defay, S. Crossley, S. KarNarayan, X. Moya, and N. D.Mathur, Adv. Mater. , 3337 (2013).[12] X. Moya, S. K.-Narayan, and N. D. Mathur, Nat. Mater. , 439(2014).[13] W. Geng, Y. Liu, X. Meng, L. Bellaiche, J. F. Scott, B. Dkhil,and A. Jiang, Adv. Mater. , 3165 (2015).[14] M. Marathe, A. Gr¨unebohm, T. Nishimatsu, P. Entel, and C.Ederer, Phys. Rev. B , 054110 (2016).[15] G. G. Guzm´an-Verri and P. B. Littlewood, APL Mater. ,064106 (2016).[16] Z. Jiang, S. Prokhorenko, S. Prosandeev, Y. Nahas, D. Wang, J.´I˜niguez, E. Defay, and L. Bellaiche, Phys. Rev. B , 014114(2017).[17] Z. Jiang, Y. Nahas, S. Prokhorenko, S. Prosandeev, D. Wang, J.´I˜niguez, and L. Bellaiche, Phys. Rev. B , 104110 (2018).[18] B. Nair, T. Usui, S. Crossley, S. Kurdi, G. G. Guzm´an-Verri, X.Moya, S. Hirose, and N. D. Mathur, Nature (London) , 468(2019).[19] J. Shi, D. Han, Z. Li, L. Yang, S.-G. Lu, Z. Zhong, J. Chen, Q.M. Zhang, and X. Qian, Joule , 1200 (2019).[20] M. M. Vopson, Solid State Commun. , 2067 (2012); J. Phys.D: Appl. Phys. , 345304 (2013).[21] E. Stern-Taulats, T. Cast´an, L. Ma˜nosa, A. Planes, N. D.Mathur, and X. Moya, MRS Bull. , 295 (2018).[22] Y. Liu, G. Zhang, Q. Li, L. Bellaiche, J. F. Scott, B. Dkhil, andQ. Wang, Phys. Rev. B , 214113 (2016).[23] I. Takeuchi and K. Sandeman, Phys. Today , 48 (2015).[24] H. Khassaf, T. Patel, and S. P. Alpay, J. Appl. Phys. , 144102(2017).[25] G. Catalan and J. F. Scott, Adv. Mater. , 2463 (2009).[26] T. Zhao, A. Scholl, F. Zavaliche, K. Lee, M. Barry, A. Doran,M. P. Cruz, Y. H. Chu, C. Ederer, N. A. Spaldin, R. R. Das, D.M. Kim, S. H. Baek, C. B. Eom, and R. Ramesh, Nat. Mater. ,823 (2006).[27] D. Lebeugle, D. Colson, A. Forget, M. Viret, A. M. Bataille,and A. Gukasov, Phys. Rev. Lett. , 227602 (2008).[28] R. J. Zeches et al ., Science , 977 (2009).[29] N. A. Spaldin, S.-W. Cheong, and R. Ramesh, Phys. Today ,38 (2010).[30] B. Xu, J. ´I˜niguez, and L. Bellaiche, Nat. Commun. , 15682(2017).[31] N. A. Spaldin and R. Ramesh, Nat. Mater. , 203 (2019).[32] C. Cazorla and J. ´I˜niguez, Phys. Rev. B , 174105 (2018).[33] A. Edstr¨om and C. Ederer, Phys. Rev. Lett. , 167201 (2020).[34] Y. Q. Zhao and H. X. Cao, J. Mater. Sci. , 5705 (2020).[35] B. Xu, D. Wang, J. ´I˜niguez, and L. Bellaiche, Adv. Funct.Mater. , 552 (2015).[36] See Supplemental Material at [URL will be inserted by pub-lisher] for more details about (i) the effective Hamiltonianmethod and parameters; (ii) a wide range of finite-temperatureproperties in disordered Bi . Nd . FeO solid solutions; (iii)finite size effects on the electrocaloric (EC) coefficient in dis-ordered Bi . Nd . FeO solid solutions; (iv) fluctuations inthe effective Hamiltonian; (v) power law with different β ; (vi)deep new insights from the phenomenological model; and (vii)another derivation to yield our phenomenological model, whichincludes Refs. [37-50].[37] B. Xu, D. Wang , H. J. Zhao, J. ´I˜niguez, X. M. Chen, and L.Bellaiche, Adv. Funct. Mater. , 3626 (2015).[38] K. Patel, S. Prosandeev, B. Xu, and L. Bellaiche, Phys. Rev. B , 214107 (2019). [39] I. A. Kornev, S. Lisenkov, R. Haumont, B. Dkhil, and L. Bel-laiche, Phys. Rev. Lett. , 227602 (2007).[40] S. Lisenkov, I. A. Kornev, and L. Bellaiche, Phys. Rev. B ,012101 (2009).[41] D. Albrecht, S. Lisenkov, W. Ren, D. Rahmedov, I. A. Kornev,and L. Bellaiche, Phys. Rev. B , 140401(R) (2010).[42] S. Prosandeev, D. Wang, W. Ren, J. ´I˜niguez, and L. Bellaiche,Adv. Funct. Mater. , 234 (2013).[43] W. Zhong, D. Vanderbilt, and K. Rabe, Phys. Rev. Lett. ,1816 (1994).[44] W. Zhong, D. Vanderbilt, and K. Rabe, Phys. Rev. B , 6301(1995).[45] I. A. Kornev, L. Bellaiche, P. E. Janolin, B. Dkhil, and E. Suard,Phys. Rev. Lett. , 157601 (2006).[46] D. Rahmedov, D. Wang, J. ´I˜niguez, and L. Bellaiche, Phys. Rev.Lett. , 037207 (2012).[47] M. P¨arnaste, M. van Kampen, R. Brucas, and B. Hj¨orvarsson,Phys. Rev. B , 104426 (2005).[48] J. H. Mokkath, Phys. Chem. Chem. Phys. , 6275 (2020).[49] B. Roˇziˇc, M. Kosec, H. Urˇsiˇc, J. Holc, B. Maliˇc, Q. M. Zhang,R. Blinc, R. Pirc, and Z. Kutnjak, J. Appl. Phys. , 064118(2011).[50] N. W. Ashcroft and N. D. Mermin, Solid State Physics (Saun-ders College, Philadelphia, 1976).[51] S. Karimi, I. M. Reaney, Y. Han, J. Pokorny, and I. Sterianou,J. Mater. Sci. , 5102 (2009).[52] I. Levin, S. Karimi, V. Provenzano, C. L. Dennis, H. Wu, T. P.Comyn, T. J. Stevenson, R. I. Smith, and I. M. Reaney, Phys.Rev. B , 020103(R) (2010).[53] I. Levin, M. G. Tucker, H. Wu, V. Provenzano, C. L. Dennis,S. Karimi, T. Comyn, T. J. Stevenson, R. I. Smith, and I. M.Reaney, Chem. Mater. , 2166 (2011).[54] Note that the local quantity η loc ( i ) is centered on the Fe-site i and is defined as η loc ( i ) = δR ionic P j σ j , where σ j character-izes the atomic distribution of Bi or Nd ion at the A site j andthe sum over j runs over the eight A nearest neighbors of theFe-site i . δR ionic describes the relative difference of ionic ra-dius between the Nd and Bi ions within the A sublattice.[55] S. Bin-Omran, I. A. Kornev, and L. Bellaiche, Phys. Rev. B ,014104 (2016).[56] Note that the denominator in this expression essentially corre-sponds to the calculated specific heat, which has several parts.The non-trivial one is associated from the potential energy ofthe subsystem described by our effective Hamiltonian (localdipoles, O rotations, inhomogeneous strains and spins), and isgiven by the width of the energy distribution as obtained fromour Monte Carlo simulations. The second, trivial part is associ- ated to the kinetic contribution of all the lattice degrees of free-dom in our material (15 per cell), as well as the contributionassociated to the potential energy of the variables not includedin our effective Hamiltonian (6 per cell); hence the term pro-portional to k B T / . Note also that Eq. (2) is not suitable todescribe situations in which the electric field induces a phasetransition; hence, we restrict its application to relatively smallelectric fields.[57] X. Moya, E. Stern-Taulats, S. Crossley, D. Gonz´alez-Alonso,S. Kar-Narayan, A. Planes, L. Ma˜nosa, and N. D. Mathur, Adv.Mater. , 1360 (2013).[58] The possible electronic contribution to the EC effect is not in-cluded here since this H eff does not incorporate electrons asdegrees of freedom. This is probably a small effect since themajor part of the polarization comes from ionic displacements.[59] A. Karchevskii, Sov. Phys. Solid State , 2249 (1962).[60] R. Pirc, Z. Kutnjak, R. Blinc, and Q. M. Zhang, J. Appl. Phys. , 074113 (2011).[61] R. Pirc, B. Roˇziˇc, J. Koruza, B. Maliˇc, and Z. Kutnjak, EPL , 17002 (2014).[62] For each temperature, the calculations provide L s and P s ,which allows to plot L s versus P s and then to obtain the deriva-tive of L s with respect to P s for any considered temperature. dL s dP s as a function of temperature is thus obtained.[63] Note that the linear temperature dependence of a ′ ( T ) is neededfor describing results far away from the Curie temperature,probably also because the oxygen octahedral tiltings are im-plicitly included in our model within some parameters of Eq.(11).[64] D. G. Tomuta, S. Ramakrishnan, G. J. Nieuwenhuys, and J. A.Mydosh, J. Phys.: Condens. Matter , 4543 (2001).[65] T. Bayaraa, Y. Yang, H. J. Zhao, J. ´I˜niguez, and L. Bellaiche,Phys. Rev. Mater. , 084404 (2018).[66] H. Sakai, J. Fujioka, T. Fukuda, D. Okuyama, D. Hashizume, F.Kagawa, H. Nakao, Y. Murakami, T. Arima, A. Q. R. Baron, Y.Taguchi, and Y. Tokura, Phys. Rev. Lett. , 137601 (2011).[67] L. Maurel, N. Marcano, E. Langenberg, R. Guzm´an, T.Prokscha, C. Mag´en, J. A. Pardo, and P. A. Algarabel, APLMater. , 041117 (2019).[68] In case of a ferromagnetic multiferroic, Eq. (11) needs to be al-tered by simply replacing the antiferromagnetic vector by themagnetization, M . All the conclusions indicated above thusstill hold but when considering the magnetic Curie temperaturerather than the N´eel one and when involving ∂M∂P s rather than ∂L s ∂P ss