calculations as implemented in VASP [-@Kresse1993; -@Kresse1996]. Projector augmented-wave method was used [-@Blochl1994]. A plane-wave basis set with the cutoff energy 520 eV was employed, and PBEsol (Perdew-Burke-Ernzerhof revised for solids) exchange-correlation functional was adoped [-@Perdew2008]. We used the lattice constant optimized in bulk calculations of BaHfO$_3$ and the internal atomic positions were fully relaxed until the force became less than 0.01 eV/Å. Details of the analytic tight-binding approximation and effective Hamiltonian description were presented in Appendix A. We employed maximally localized Wannier functions [-@Marzari1997; -@Souza2001; -@Mostofi2008] to further analyze the results of the first-principles calculations. The Wannier functions were constructed for $d$ orbitals of Ta in one set, and $p$ orbitals of three neighboring O as well as $d$ orbitals of Ta in the other set.

**Acknowledgments**

We thank Jung Hoon Han, Changyoung Kim, Choong-Hyun Kim, Hyeong-Do Kim, Minu Kim, Hyun-Woo Lee, Hosub Jin, Seung Ryong Park, Cai-Zhuang Wang, Hong Yao, and Jaejun Yu for fruitful discussions and comments. This work was supported by the NRF of the MSIP of the Korean government Grant No. 2006-0093853 (M.K., J.I.) and IBS-R009-Y1 (S.B.C.). Research at Ames laboratory (M.K.) was supported by the U.S. Department of Energy, Office of Basic Energy Sciences, 