Electrocaloric response of KNbO3 from a first-principles effective Hamiltonian

Materials Science and Engineering B 196 (2015) 40–43

Contents lists available at ScienceDirect

Materials Science and Engineering B journal homepage: www.elsevier.com/locate/mseb

Short communication

Electrocaloric response of KNbO3 from a first-principles effective Hamiltonian J.A. Barr, S.P. Beckman ∗ Department of Materials Science and Engineering, Iowa State University, Ames, IA 50011, United States

a r t i c l e

i n f o

Article history: Received 30 July 2014 Received in revised form 25 January 2015 Accepted 10 February 2015 Available online 5 March 2015 Keywords: KNbO3 Electrocaloric effect Effective Hamiltonian Perovskite Ferroelectric

a b s t r a c t The electrocaloric response of KNbO3 is calculated using a first-principles based effective-Hamiltonian. Both indirect and direct molecular dynamics simulations are used to determine the adiabatic temperature change as a function of temperature and applied electric field. Both the magnitude of the electrocaloric response, Tmax , and the temperature where the maximum is observed, T* , are comparable for both molecular dynamics approaches. The results are well represented by a simple second-order regression. © 2015 Elsevier B.V. All rights reserved.

1. Introduction Solid-state cooling is a promising alternative to traditional, lowefficiency, refrigeration techniques [1,2]. One approach involves using the electrocaloric effect (ECE) in which an electric field is used to drive a pyroelectric material through an order–disorder transformation. As a result of this adiabatic change in entropy there is a change in temperature, T. A refrigeration cycle is achieved by cycling the applied electric field as the pyroelectric material alternates between hot and cold thermal reservoirs. There are many theoretical [3–7] and experimental [8–10] studies of the ECE in BaTiO3 (BTO). There also are engineering studies examining the feasibility of implementing the electrocaloric effect for cooling [11]. For example, a prototype electrocaloric oscillatory refrigeration device has been shown to produce a rise in temperature of 2.25 K upon the application of an 80 MV/m field and a 2.15 K decrease in temperature upon removal of the field [12]. Other materials have received significantly less attention than BaTiO3 . The perovskite KNbO3 (KNO) is known to have good optical and piezoelectric properties with potential applications as a pyroelectric detector, a display device, and an ultrasonic transducer [13–15]. Here we demonstrate the use of a first-principles based effective Hamiltonian molecular dynamics simulation to calculate the ECE of KNO.

∗ Corresponding author. Tel.: +1 510 708 4390; fax: +1 515 294 5444. http://dx.doi.org/10.1016/j.mseb.2015.02.004 0921-5107/© 2015 Elsevier B.V. All rights reserved.

Following this brief introduction we present the methods used to calculate the ECE. In Section 3 the results are presented and discussed and in Section 4 the paper is summarized. 2. Methods A course-grained model of the perovskite compound is used in which the energy is expressed in terms of a local soft-mode amplitude vector and a local acoustic displacement vector. This effective Hamiltonian is derived in Refs. [16–18] and the parameter set for KNO is found in Ref. [19]. The effective Hamiltonian is given by H eff =

∗  Mdipole

2

R,˛

u˙ 2˛ (R) +

∗  Macoustic

2

w˙ ˛2 (R)

R,˛

+V self ({u}) + V dpl ({u}) + V short ({u}) + V elas,homo (1 , . . ., 6 ) +V elas,inho ({w}) + V coup,homo ({u}, 1 , · · ·, 6 ) +V coup,inho ({u}, {w}).

(1)

Here, w(R) is the local acoustic displacement vector and u( R) is the local soft mode vector. M∗ acoustic



M∗

dipole

2



u˙ 2 (R) R,˛ ˛

and

w˙ 2 (R) represent the kinetic energies of the local 2 R,˛ ˛ acoustic and local soft mode vectors. The local mode self energy is given by Vself ({ u}), the long range dipole–dipole interaction is given by Vdpl ({ u}), and the short range interaction between local soft modes is given by Vshort ({ u}). The homogenous strain elastic energy is given by Velas, homo (1 , . . ., 6 ), while the inhomogeneous strain elastic energy is given by V elas,inho ({w}). Finally, the

J.A. Barr, S.P. Beckman / Materials Science and Engineering B 196 (2015) 40–43

41

coupling between local soft modes and homogenous strain is given by Vcoup,homo ({ u}, 1 , . . ., 6 ), while the coupling between local soft modes and inhomogeneous strains is given by V coup,inho ({u}, {w}). The framework for the molecular dynamics (MD) simulation and effective Hamiltonian is encoded in the software package feram [18,19].1 Methods for predicting the ECE are described in detail in Ref. [20]. The “indirect method” involves maintaining the system in a canonical ensemble to determine the polarization as a function of temperature and electric field. The adiabatic temperature change is determined from the Maxwell relation, T =

−1 Cv



E2

T E1

∂Pz dEz , ∂T

(2)

where ∂Pz /∂T is the pyroelectric coefficient, E = E2 − E1 is the change in electric field that induces the entropy change, and Cv is the heat capacity. The “direct method” involves switching from a constant temperature and field canonical ensemble to a microcanonical ensemble as the electric field is switched. This allows for direct determination of the temperature from the MD simulations. These two methods are used to calculate T for KNO. For both methods, a negative effective pressure of p = −0.001T GPa is applied to simulate thermal expansion. For the indirect method, a supercell of size N = Lx × Ly × Lz = 16 × 16 × 16 is employed. The system is allowed to equilibrate for 20,000 time steps and then time averaged for 100,000 times steps using a time step of 2 fs. The heat capacity is approximated using the value of Cv = 4.617 g cm−3 × 0.767 J g−1 K−1 or 3.54 J cm−3 K−1 [21]. The direct method requires a larger supercell to accommodate the thermal fluctuations that occur when the system is switched to the microcanonical ensemble. Here a system of size N = Lx × Ly × Lz = 96 × 96 × 96 is used. The system is allowed to equilibrate for 40,000 time steps in a canonical ensemble and then after switching to the microcanonical ensemble another 40,000 time steps are given to allow equilibrium to be regained. The temperature is then determined by time averaging for 60,000 time steps. A time step of 2 fs is used in the calculations. For the direct method, the specific heat is automatically included in the calculation, but is scaled by 1/5 due to the reduced degrees of freedom inherent in the course grained model [20]. The results reported here all have had this scaling applied. 3. Results and discussion For the indirect method simulations are run for electric fields ranging from 25 to 300 kV/cm, using an increment of 25 kV/cm, and temperatures ranging from 300 K to 1600 K, at increments of 5 K. The polarization in the z-direction is determined and is plotted as in Fig. 1(a). The Curie temperature is around 700 K, although the transformation spreads over a larger temperature range under an applied field. This compares very well with the measured transformation temperature of 701 K reported by Günter [22]. Our calculated spontaneous polarization is around 55 ␮C cm−2 , which is slightly larger than reported by Günter, 41␮C cm−2 in Ref. [23] and 30 ␮C cm−2 in Ref. [22]. Birol et al. [24] reported measuring a spontaneous polarization of 25 ␮C cm−2 in polycrystalline KNbO3 . Ge et al. [25] studied the impact of crystal size on the polarization in KNbO3 and found a spontaneous polarization of 40.3 ␮C cm−2 . Kleemann et al. [26] reported a spontaneous polarization of 37 ␮C cm−2 . The variation in these values is due to difference in the samples,

1 The feram software is distributed freely under the GNU general public license at http://loto.sourceforge.net/feram/.

Fig. 1. (a) The temperature dependence of polarization in the z-direction at constant external electric fields. A weighted spline is fit to the collected data and used to determine (b) the pyroelectric coefficient as a function of temperature at constant electric field.

and more importantly, the magnitude of the applied electric fields. Theoretical work, based on density functional theory methods with conditions of no applied field and T = 0 K, find spontaneous polarizations of 37 ␮C cm−2 [15] and 48 ␮C cm−2 [27]. These results give us confidence in the quality of the coarse-grained model used here. A weighted spline is fit to our results in Fig. 1a and ∂Pz /∂T is determined to give the pyroelectric coefficient shown in Fig. 1(b). The data in Fig. 1(b) is rearranged to give the pyroelectric coefficient as a function of applied electric field, E, at constant temperatures, as plotted in Fig. 2. The curves in Fig. 2 are integrated according to Eq. (2) to determine the T. The resulting T are presented in Fig. 3(a) for the case where E2 is held constant at 300 kV/cm and E1 is varied from 50 to 250 kV/cm. For comparison, T is also determined for the case where E held constant at 100 kV/cm and E1 is varied from 50 to 200 kV/cm, as shown in

Fig. 2. The pyroelectric coefficient is replotted as a function of electric field at constant temperatures.

42

J.A. Barr, S.P. Beckman / Materials Science and Engineering B 196 (2015) 40–43

indirect methods the maximum T is found near T = 800 K when the field is selected to be relatively low, ranging from 50 to 150 kV/cm. The direct method predicts T that are between 25% and 30% larger than the indirect. This is possibly due to error introduced by the use of a temperature independent heat capacity used in the indirect method. However, the qualitative trends between the two are in reasonable agreement. These results are sufficiently well-behaved to be represented by a simple polynomial regression fit. The predicted maximum, Tmax , and the temperature where it occurs, T* , can be written as a function of the electric fields. For (E1 < E2 ) the resulting expressions are, Tmax (E1 , E2 ) = 8.23 − 4.58 × 10−2 E1 − 6.44 × 10−5 E21 − 8.05 × 10−4 E2 + 8.29 × 10−5 E1 E2 + 3.60 × 10−5 E22 ,

(3)

and T ∗ (E1 , E2 ) = 657.69 + 8.71 × 10−1 E1 − 5.42 × 10−4 E21 + 6.92 × 10−1 E2 − 5.69 × 10−4 E1 E2 − 1.08 × 10−3 E22 .

Fig. 3. (a) The predicted electrocaloric effect by the indirect method where E2 is held constant at 300 kV/cm and E1 is varied. (b) The predicted electrocaloric effect by the indirect method where E is held constant at 100 kV/cm, but E1 and E2 are varied.

(4)

The coefficients of determination for these are R2 = 0.999956 and R2 = 0.999998 respectively, which gives us confidence in the expressions capabilities of accurately predict the electrocaloric behavior of KNO. The coefficients of higher order terms are all at least two orders of magnitude smaller than those reported here. 4. Summary and conclusions

Fig. 3(b). It is noteworthy that for a constant E both the magnitude of T and the temperature of the maximum T vary. This demonstrates the importance of the absolute values of E1 and E2 on the final ECE, i.e., the impact of applied field on the phase transformation, as shown in Eq. (2). As expected, it is clear that for a given electric field range the temperature with the maximum T is observed in Fig. 3 for the temperatures with the greatest peaks in Fig. 2. For example, when the field ranges from 100 to 200 kV/cm, the greatest pyroelectric coefficient occurs around 825 K in Fig. 2, close to the temperature seen in Fig. 3(b) where the maximum T occurs for the same electric field range. Increasing the electric field range, E, naturally increases the magnitude of T and also shifts it toward the temperature T = 775 K, the Curie temperature. The direct method ECE calculations are shown in Fig. 4, which can be directly compared to Fig. 3(b). For both the direct and 0 50 ... 150 kV/cm 100 ... 200 kV/cm 150 ... 250 kV/cm 200 ... 300 kV/cm

-2 -4

ΔT [K]

-6 -8

-10 -12 -14 600

700

800

900

1000

1100

1200

1300

Temperature [K] Fig. 4. The ECE as calculated using the direct method. E is held constant at 100 kV/cm, but E1 and E2 are varied.

We have predicted the ECE in KNO using a first-principles based effective Hamiltonian that has been implemented in a molecular dynamics framework. The minor differences between the two molecular dynamics method are possibly due to approximations, including the use of a constant temperature, experimentally determined heat capacity. Differences between these calculations and experimental measurements may arise from the course grained model’s idealization of domain switching as well as discrepancies due to the first-principle determined energies and structures, although to the best of the authors’ knowledge such measurements have yet to be reported in the literature. Qualitatively, the predicted trends are entirely consistent with those observed for BTO [3,20,28]. The overall magnitude of the ECE is directly proportional to the magnitude of E with the largest change in electric fields giving the largest ECE. The ECE depends not only on the range, but also on the absolute values of field. For example, if the range is held constant at E =100 kV/cm, then a larger T is achieved if the absolute fields, E1 and E2 , are relatively low as shown in Fig. 3(b). This also results in a maximum T at lower temperatures. It is also seen that the maximum T occurs above the Curie temperature and is shifted to higher temperatures as the magnitude of the electric field increases, agreeing with previous experimental data for single crystal BTO [28]. This similarity is consistent with BTO and KNO having similar bonding, with both having their temperature dependent ferroelectric polymorphs determined by the hybridized covalent bonding between the B- and O-sites. For quantitative comparison, switching BTO’s electric field from 60 to 160 kV/cm results in a temperature change of approximately 6 K [20] and switching KNO induces a ECE T of 6.8 K. This can be explained by examining the polarization of the first stable ferroelectric phase below the Curie temperature, which is tetragonal for both. In the case of BTO, the tetragonal phase has a polarization of approximately 37 ␮C cm−2 and KNO is approximately 55 ␮C cm−2 . This results in BTO and KNO having similar pyroelectric coefficients, yet KNO will be slightly greater, approximately −0.9 ␮C cm−2 K−1

J.A. Barr, S.P. Beckman / Materials Science and Engineering B 196 (2015) 40–43

for E = 25 kV/cm. Thus, it is demonstrated here that materials with similar bonding properties will have qualitatively and quantitatively similar ECE behavior.

[4] [5] [6] [7]

Acknowledgements The authors gratefully acknowledge Prof. Takeshi Nishimatsu’s effort developing and distributing the feram software package for the scientific community. The “free” mode of software development provided by the GNU General Public License is essential for scientific progress in the modern era. We also appreciate the many interesting discussions we have had with Prof. Nishimatsu regarding functional properties of perovskite compounds. The authors were supported by the National Science Foundation through grant DMR-1105641. JAB acknowledges travel support provided by the National Science Foundation through a supplement to grant DMR-1037898. Computational resources were provided in part by the Center for Computational Materials Science, Institute for Materials Research (CCMS-IMR) at Tohoku University and in part by the Iowa State University, Computational Ceramic Engineering (ISU-CCE) program. References [1] R. Chukka, S. Vandrangi, S. Shannigrahi, L. Chen, Europhys. Lett. 103 (2013) 47011. [2] S.-G. Lu, Q. Zhang, Adv. Mater. 21 (2009) 1983–1987. [3] S. Beckman, L. Wan, J.A. Barr, T. Nishimatsu, Mater. Lett. 89 (2012) 254–257.

[8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28]

43

H.-X. Cao, Z.-Y. Li, J. Appl. Phys. 106 (2009). G. Akcay, S.P. Alpay, J.V. Mantese, G.A. Rossetti, Appl. Phys. Lett. 90 (2007). J.H. Qiu, Q. Jiang, J. Appl. Phys. 105 (2009). J.A. Barr, T. Nishimatsu, S.P. Beckman, Computational modeling of electrocaloric effect for solid-state refrigeration, in: Nanoscale Thermoelectric Materials: Thermal and Electrical Transport, and Applications of Solid-State Cooling and Power Generation. Symposium (Mater. Res. Soc. Symposium Proceedings), vol. 1543, 2013, pp. 39. Y. Bai, K. Ding, G.-P. Zheng, S.-Q. Shi, J.-L. Cao, L. Qiao, AIP Adv. 2 (2012). X. Moya, E. Stern-Taulats, S. Crossley, D. Gonzlez-Alonso, S. Kar-Narayan, A. Planes, L. Ma-osa, N.D. Mathur, Adv. Mater. 25 (2013) 1360–1365. Y. Bai, X. Han, X.-C. Zheng, L. Qiao, Sci. Rep. 3 (2013). Y. Jia, Y. Sungtaek Ju, Appl. Phys. Lett. 100 (2012). H. Gu, X. Qian, X. Li, B. Craven, W. Zhu, A. Cheng, S.C. Yao, Q.M. Zhang, Appl. Phys. Lett. 102 (2013). B. Zysset, I. Biaggio, P. Günter, J. Opt. Soc. Am. B 9 (1992) 380–386. K. Nakamura, T. Tokiwa, Y. Kawamura, J. Appl. Phys. 91 (2002). L. Wan, T. Nishimatsu, S.P. Beckman, J. Appl. Phys. 111 (2012). R.D. King-Smith, D. Vanderbilt, Phys. Rev. B 49 (1994) 5828–5844. W. Zhong, D. Vanderbilt, K.M. Rabe, Phys. Rev. B 52 (1995) 6301–6312. T. Nishimatsu, U.V. Waghmare, Y. Kawazoe, D. Vanderbilt, Phys. Rev. B 78 (2008). T. Nishimatsu, feram at sourceforge.net, 2007–2014. T. Nishimatsu, J.A. Barr, S.P. Beckman, J. Phys. Soc. Jpn. 82 (2013). D.N. Nikogosyan, Nonlinear Optical Crystals: A Complete Survey, Springer, 2005. P. Gunter, J. Appl. Phys. 48 (1977). P. Gunter, Jpn. J. Appl. Phys. 16 (1977). H. Birol, D. Damjanovic, N. Setter, J. Am. Cream. Soc. 88 (2005). H. Ge, Y. Huang, P. Hou, H. Xiao, M. Zhu, R. Soc. Chem. Adv. 4 (2014). W. Kleemann, F.J. Schafer, M.D. Fontana, Phys. Rev. B 30 (1984) 1148–1154. P. García-Fernández, P. Aguado-Puente, J. Junquera, Phys. Rev. B 87 (2013), 085305. X. Moya, E. Stern-Taulats, S. Crossley, D. Gonzalez-Alonso, S. Kar-Narayan, A. Planes, L. Manosa, N.D. Mathur, Phys. Rev. B 25 (2013) 1360–1365.