Atomistic Description for Temperature–Driven Phase Transitions in BaTiO{}_{3}

# Atomistic Description for Temperature–Driven Phase Transitions in BaTiO3

## Abstract

Barium titanate (BaTiO) is a prototypical ferroelectric perovskite that undergoes the rhombohedral–orthorhombic–tetragonal–cubic phase transitions as the temperature increases. In this work, we develop a classical interatomic potential for BaTiO within the framework of the bond–valence theory. The force field is parameterized from first–principles results, enabling accurate large–scale molecular dynamics (MD) simulations at finite temperatures. Our model potential for BaTiO reproduces the temperature–driven phase transitions in isobaric–isothermal ensemble () MD simulations. This potential allows the analysis of BaTiO structures with atomic resolution. By analyzing the local displacements of Ti atoms, we demonstrate that the phase transitions of BaTiO exhibit a mix of order–disorder and displacive characters. Besides, from detailed observation of structural dynamics during phase transition, we discover that the global phase transition is associated with changes in the equilibrium value and fluctuations of each polarization component, including the ones already averaging to zero, Contrary to the conventional understanding that temperature increase generally causes bond–softening transition, the polarization component exhibits a bond–hardening character during the orthorhombic to tetragonal transition. These results provide further insights about the temperature–driven phase transitions in BaTiO.

###### pacs:
Valid PACS appear here
1

## I Introduction

BaTiO is a ferroelectric perovskite with promising applications in electronic devices, such as non–volatile memory, high– dielectrics, and piezoelectric sensors Yun et al. (2002); Chen et al. (2007); Buessem et al. (1966); Karaki et al. (2007). Therefore, it is of great significance to investigate and understand the structural and electronic properties of BaTiO for designed material optimization and device engineering. First–principles density functional theory (DFT) has served as a powerful method to understand the electronic structures of ferroelectric materials Cohen (1992); Gonze and Lee (1997); Saha et al. (2000); Kolpak et al. (2008); Martirez et al. (2012); Morales et al. (2014). Due to the expensive computational cost, the application of DFT methods is currently limited to system of fairly small size at zero Kelvin. Many important dynamical properties, such as domain wall motions and temperature–driven phase transitions, are beyond the capability of conventional first–principles methods. An effective Hamiltonian method was developed to study finite–temperature properties of BaTiO Zhong et al. (1994, 1995); Nishimatsu et al. (2008); Fu and Bellaiche (2003). To apply this method, the subset of dynamical modes that determine a specific property should be known a priori. Molecular dynamics (MD) simulations with an atomistic potential accounting for all the modes offer distinct advantages, especially in providing detailed information about atomic positions, velocities and modifications of chemical bonds due to a chemical reaction or thermal excitation. The shell model for BaTiO has been developed Tinte et al. (1999); Tinte and Stachiotti (2001); Tinte et al. (2000); Zhang et al. (2010, 2009). However, due to the low mass assigned to the shell, a small time step in MD simulations is required to achieve accurate results, which limits the time and length scales of the simulations.

Recently, we developed a bond–valence (BV) model potential for oxides based on the bond valence theory Grinberg et al. (2002); Shin et al. (2008, 2005); Liu et al. (2013a, b). The force fields for many technologically important ferroelectric materials, PbTiO, PbZrO and BiFeO Liu et al. (2013a, b); Shin et al. (2005); Grinberg et al. (2002); Chen et al. (2015), have been parameterized based on results from DFT calculations. A typical force field requires no more than 15 parameters and can be efficiently implemented, which enables simulations of systems with thousands of atoms under periodic boundary conditions Liu et al. (2013c); Xu et al. (2015). The development of an accurate classical potential for BaTiO has proven to be difficult, mainly due to the small energy differences among the four phases (rhombohedral, orthorhombic, tetragonal, and cubic) Kwei et al. (1993); Diéguez et al. (2004); Von Hippel (1950). In this paper, we apply the bond–valence model to BaTiO and parameterize the all–atom interatomic potential to first–principles data. Our model potential for BaTiO is able to reproduce the rhombohedral–orthorhombic–tetragonal–cubic (R-O-T-C) phase transition sequence in isobaric–isothermal ensemble () MD simulations. The phase transition temperatures agree reasonably well with previous theoretical results Tinte et al. (1999). We further examine the temperature dependence of the local displacements of Ti atoms and discover several features of the phase transitions of BaTiO: the phase transitions of BaTiO involve both order–disorder and displacive characters; at the moment that the phase transition of the crystal occurs, all the polarization components undergo phase transitions, even for the nonpolar ones; and temperature increase can also cause bond–hardening for a certain component.

## Ii Methods

The bond–valence model potential is developed based on the conservation principles of bond valence and bond–valence vector. The bond valence, , reflects the bonding strength and can be calculated based on the bond length, , with Brown (2009); Grinberg et al. (2002); Brown and Shannon (1973); Brown and Wu (1976); Liu et al. (2013a, b); Shin et al. (2008, 2005)

 Vij=(r0,ijrij)Cij (1)

where and are the labels for atoms; and are Brown’s empirical parameters. The bond–valence vector is defined as a vector lying along the bond, , where is the unit vector pointing from atom to atom . The total energy () consists of the Coulombic energy (), the short–range repulsive energy (), the bond–valence energy (), the bond–valence vector energy (), and the angle potential (Liu et al. (2013a, b); Shin et al. (2008, 2005):

 E=Ec+Er+EBV+EBVV+Ea (2)
 Ec=∑i
 Er=∑i
 EBV=∑iSi(Vi−V0,i)2 (5)
 EBVV=∑iDi(W2i−W20,i)2 (6)
 Ea=kNoxygen∑i(θi−180∘)2 (7)

where is the bond–valence sum (BVS), is the bond–valence vector sum (BVVS, shown in FIG. 1 (a), (b)), is the ionic charge, is the short–range repulsion parameter, and are scaling parameters with the unit of energy, is the spring constant and is the O–O–O angle along the common axis of two adjacent oxygen octahedra (FIG. 1 (c)). The bond-valence energy captures the energy penalty for both overbonded and underbonded atoms. The bond-valence vector energy is a measure of the breaking of local symmetry, which is important for correctly describing the ferroelectricity. and are preferred or target values of BVS and BVVS for atom in the ground–state structure, which can be calculated from DFT directly. It is noted that the and can be related to the moments of the local density of states in the framework of a tight binding model, providing a quantum mechanical justification for these two energy terms Liu et al. (2013a, b); Finnis and Sinclair (1984); Brown (2009); Harvey et al. (2006). The angle potential is used to account for the energy cost associated with the rotations of oxygen octahedra.

We followed the optimization protocol developed in previous studies Liu et al. (2013a, b). The optimal values of force-field parameters , , and , are acquired by minimizing the difference between the DFT energies/forces and the model–potential energies/forces for a database of BaTiO structures. All DFT calculations are carried out with the plane–wave DFT package Quantum–espresso Giannozzi et al. (2009) using the Perdew–Burke–Ernzerhof functional modified for solids (PBEsol) Perdew et al. (2008) and optimized norm–conserving pseudopotentials generated by the Opium package Opi (). A plane–wave cutoff energy of 50 Ry and 44 Monkhorst–Pack –point mesh Monkhorst and Pack (1976) are used for energy and force calcualtions. The database consists of 40–atom 222 supercells with different lattice constants and local ion displacements. The final average difference between DFT energy and model–potential energy is 1.35 meV/atom.

## Iii Performance of the classical potential

The optimized parameters are listed in TABLE 1. The performance of the obtained force field is examined by investigating the temperature dependence of lattice constants (, and ), component–resolved local displacements of Ti atoms (, , and ), and the three components of the total polarization (, , and ). We carry out MD simulations using a 101010 supercell (5000 atoms) with the temperature controlled via the Nosé–Hoover thermostat and the pressure maintained at 1 atm via the Parrinello–Rahman barostat Parrinello and Rahman (1980). As shown in FIG. 2, the simulations clearly reveal four distinct phases under different temperature ranges and three first–order phase transitions. Below 100 K, the displacements of Ti atoms and the overall polarization of the supercell are along [111] direction (), characteristic of the rhombohedral phase. At 100 K, the component of the total polarization, , becomes approximately 0, indicating a phase transition from rhombohedral to orthorhombic (, ). As the temperature increases further to 110 K, the total polarization aligns preferentially along direction (, ) and the lattice constants have . The supercell stays tetragonal until 160 K at which point the ferroelectric–paraelectric phase transition occurs. The phase transition temperatures match well with those predicted by the shell model Tinte et al. (1999) (TABLE 2). The underestimation of Curie temperature in MD simulations has been observed previously and is likely due to the systematic error of density functional used for force field optimization Liu et al. (2013a); Sepliarsky et al. (2005). We extract the averaged lattice constants at finite temperatures from MD simulations and find that they are in good agreement (error less than 1) with the PBEsol values (TABLE 3).

Domain walls are interfaces separating domains with different polarities. They are important topological defects and can be moved by applying external stimulus Liu et al. (2013c); Xu et al. (2015). The domain wall energy for a 180 wall obtained from our MD simulations is 6.63 mJ/m, which is comparable to PBEsol value, 7.84 mJ/m. This indicates that our atomistic potential can be used for studying the dynamics of ferroelectric domain walls in BaTiO. All these results demonstrate the robustness of this developed classical potential.

## Iv Atomistic features of different phases

To provide an atomistic description of the different phases of BaTiO, we analyze the distribution of local displacements of Ti atoms in each phase. Ti displacement is defined as the distance between the Ti atom and the center of the oxygen octahedral cage of a unit cell, which scales with the magnitude of polarization.

In FIG. 3 (a), we plot the distributions of Ti displacements (). It can be seen that in all four phases, the distribution is approximately a Gaussian curve whose peak shifts toward lower values as the temperature increases. This suggests that the temperature–driven phase transition has a displacive character. It is noted that the distribution of magnitudes is peaked at non–zero value even in the paraelectric cubic phase, suggesting that most Ti atoms are still locally displaced at high temperature, and that the overall net zero polarization is the result of an isotropic distribution of local dipoles along different directions. This confirms the order–disorder character for BaTiO at high temperature.

We can categorize the phase of each unit cell based on the local displacement of Ti atom. The categorization criteria are
(1) If Å, the unit cell is considered to be paraelectric cubic;
(2) For a ferroelectric unit cell, the –th component is considered to be ferroelectric if . The rhombohedral, orthorhombic, and tetragonal unit cells have three, two, and one ferroelectric component(s), respectively.
The results are shown in FIG. 3 (b). At 30 K, the supercell is made only from rhombohedral unit cells, showing that the rhombohedral phase is the ground–state structure. As the temperature increases, the supercell becomes a mixture of the four phases. It should be noted that the cubic unit cell with nearly–zero local Ti displacement seldom appears, because a cubic unit cell is energetically less favorable. The relative energies of the four phases of BaTiO from PBEsol DFT calculations are listed in TABLE 4. It can be seen that the energy differences between the tetragonal, orthorhombic and rhombohedral unit cells are small (within several meV per unit cell) Cohen and Krakauer (1990); Cohen (1992). Due to the thermal fluctuations, the populations of higher–energy ferroelectric phases (tetragonal and orthorhombic) increase as temperature increases. Above the ferroelectric–paraelectric transition temperature, locally ferroelectric unit cells are still favored over paraelectric due to the relatively high energy of cubic, the high–symmetry structure.

In FIG. 4, the distributions of Ti displacements along the three axes are plotted. At 100 K, BaTiO is at the rhombohedral phase and the distributions of Ti displacements are Gaussian–like. As the temperature increases, the phase changes to orthorhombic. The average of the polarization component shifts to zero, indicating a displacive phase transition. Besides, the standard deviation increases and the center of the distribution curve becomes flatter. For the cubic phase, the center of the Ti displacement distribution curve is also flat. As shown in FIG. 5, the center–flat curve is a summation of a Gaussian curve centering at zero, and a double–peak curve. The latter is characteristic of order–disorder transition Liu et al. (2013c). These results further demonstrate that phase transitions of BaTiO have a mix of order–disorder and displacive characters Gaudon (2015); Ehses et al. (1981); Jun et al. (1988); Comes et al. (1968); Itoh et al. (1985); Stern (2004); Kwei et al. (1993).

## V Features of the phase transitions

To investigate the structural dynamics during phase transitions in more detail, we conducted MD simulations with varying temperatures. In three different sets of simulations, the temperatures were increased from 100 K to 110 K (R to O), 110 K to 120 K (O to T) and 155 K to 165 K (T to C) respectively. The temperature was controlled by the Nosé–Hoover thermostat with a thermal inertia parameter =10 and the 10 K temperature change was accomplished in 60 ps. We analyze the temperature dependence of Ti displacement distributions along three axes. The dynamics of Ti displacement distributions during the phase transitions are plotted in FIG. 6. The time evolution of the average and standard deviation of the Ti displacement distributions are shown in FIG. 7.

Phase transition occurs when one component undergoes polar–nonpolar transition. The first column (from 100 K to 110 K) shows the changes of Ti displacement distributions during the rhombohedral to orthorhombic phase transition. In the and direction, the averages of the distribution shift up, which is a characteristic of displacive transition. Meanwhile, in the direction, the average becomes zero and the variance becomes significantly larger, indicating that the transition is a mix of displacive and bond–softening characters Gaffney and Chapman (2007). For the orthorhombic to tetragonal phase transition (second column), the transition of the component, which is a polar–nonpolar transition, includes both displacive and bond–softening features. For the component, the transition involves both an increase of the average and a decrease of the standard deviation. For the direction, even though the Ti displacement distribution is centered at zero above and below the transition, the Ti displacements are located closer to zero, indicating an increase in bond hardness. From 155 K to 165 K, there is also a bond–hardness–changing transition for the components ( and ) with zero averages. We collectively refer to ‘bond–softening’ and ‘bond–hardening’ as ‘bond–hardness–changing’.

Based on the features of the Ti displacement distributions at different phases, the schematic representation of the thermal excitation between different energy surfaces is presented in FIG. 8. From our results, the characteristics of BaTiO phase transition can be summarized as: (1) For BaTiO, the mechanisms of phase transitions include both bond–hardness–changing and displacive transition. The sudden shifts of the average and standard deviation correspond to displacive with some order–disorder contribution and bond–hardness–changing transitions respectively; (2) Unlike the conventional understanding that thermal excitation usually causes bond–softening, increasing temperature can also cause bond hardening. The component of polarization during the orthorhombic to tetragonal transition is an example of this case. (3) When the phase transition occurs, each component of polarization undergoes a change, even for the component(s) which is(are) non–polar before and after the transition. The transition(s) that each component undergoes are listed in TABLE 5.

## Vi Conclusion

In this work, we develop a classical atomistic potential for BaTiO based on the bond valence model. Molecular dynamics simulation with this optimized potential can not only reproduce the temperature–driven phase transitions, but can also be a powerful tool in studying the phase transition process with high temporal and spatial resolutions. The detailed analysis of the local displacements of Ti atoms reveals that in each phase (including the paraelectric phase), the majority of Ti atoms are locally displaced, and the phase transitions in BaTiO exhibit a mixture of order–disorder and displacive character. The distribution of Ti displacement is a Gaussian curve or a curve involving a Gaussian and a double peak one. Both the average and standard deviation are features of each specific phase, and they can be the order parameters for describing the Gibbs free energy. By analyzing the dynamics of Ti displacement distributions during phase transition, we discover several rules of BaTiO phase transitions: the global phase transition is associated with significant changes in each component, even for the components which are nonpolar, and the orthorhombic to tetragonal transition exhibits a bond–hardening character in the component, which is opposite to the conventional understanding that temperature increase generally causes bond–softening transition.

## Acknowledgments

Y.Q. was supported by the U.S. National Science Foundation, under Grant No. CMMI1334241. S.L. was supported by the U.S. National Science Foundation, under Grant No. CBET1159736 and Carnegie Institution for Science. I.G. was supported by the Office of Naval Research, under Grant No. N00014–12–1–1033. A.M.R. was supported by the Department of Energy, under Grant No. DE–FG02–07ER46431. Computational support was provided by the High–Performance Computing Modernization Office of the Department of Defense and the National Energy Research Scientific Computing Center of the Department of Energy.

### Footnotes

1. preprint: APS/123-QED

### References

1. W. S. Yun, J. J. Urban, Q. Gu,  and H. Park, Nano Lett. 2, 447 (2002).
2. K.-H. Chen, Y.-C. Chen, Z.-S. Chen, C.-F. Yang,  and T.-C. Chang, Appl. Phys. A Mater. Sci. 89, 533 (2007).
3. W. Buessem, L. Cross,  and A. Goswami, J. Am. Ceram. Soc. 49, 33 (1966).
4. T. Karaki, K. Yan,  and M. Adachi, Jpn. J. of Appl. Phys. 46, 7035 (2007).
5. R. E. Cohen, Ferroelectrics 136, 65 (1992).
6. X. Gonze and C. Lee, Phys. Rev. B 55, 10355 (1997).
7. S. Saha, T. P. Sinha,  and A. Mookerjee, Phys. Rev. B 62, 8828 (2000).
8. A. M. Kolpak, D. Li, R. Shao, A. M. Rappe,  and D. A. Bonnell, Phys. Rev. Lett. 101, 036102 (2008).
9. J. M. P. Martirez, E. H. Morales, W. A. Al-Saidi, D. A. Bonnell,  and A. M. Rappe, Phys. Rev. Lett. 109, 256802 1 (2012).
10. E. H. Morales, J. M. P. Martirez, W. A. Saidi, A. M. Rappe,  and D. A. Bonnell, ACS nano 8, 4465 (2014).
11. W. Zhong, D. Vanderbilt,  and K. M. Rabe, Phys. Rev. Lett. 73, 1861 (1994).
12. W. Zhong, D. Vanderbilt,  and K. M. Rabe, Phys. Rev. B 52, 6301 (1995).
13. T. Nishimatsu, U. V. Waghmare, Y. Kawazoe,  and D. Vanderbilt, Phys. Rev. B 78, 104104 (2008).
14. H. Fu and L. Bellaiche, Phys. Rev. Lett. 91, 257601 (2003).
15. S. Tinte, M. G. Stachiotti, M. Sepliarsky, R. L. Migoni,  and C. O. Rodriquez, J. Phys.: Condens. Matter 11, 9679 (1999).
16. S. Tinte and M. G. Stachiotti, Phys. Rev. B 64, 235403 (2001).
17. S. Tinte, M. Stachiotti, M. Sepliarsky, R. Migoni,  and C. Rodriguez, Ferroelectrics 237, 41 (2000).
18. Y. Zhang, J. Hong, B. Liu,  and D. Fang, Nanotechnology 21, 015701 (2010).
19. Y. Zhang, J. Hong, B. Liu,  and D. Fang, Nanotechnology 20, 405703 (2009).
20. I. Grinberg, V. R. Cooper,  and A. M. Rappe, Nature 419, 909 (2002).
21. Y.-H. Shin, J.-Y. Son, B.-J. Lee, I. Grinberg,  and A. M. Rappe, J. Phys.: Condens. Matter 20, 015224 (2008).
22. Y.-H. Shin, V. R. Cooper, I. Grinberg,  and A. M. Rappe, Phys. Rev. B 71, 054104 (2005).
23. S. Liu, I. Grinberg,  and A. M. Rappe, J. Physics.: Condens. Matter 25, 102202 (2013a).
24. S. Liu, I. Grinberg, H. Takenaka,  and A. M. Rappe, Phys. Rev. B 88, 104102 (2013b).
25. F. Chen, J. Goodfellow, S. Liu, I. Grinberg, M. C. Hoffmann, A. R. Damodaran, Y. Zhu, P. Zalden, X. Zhang, I. Takeuchi, A. M. Rappe, L. W. Martin, H. Wen,  and A. M. Lindenberg, Adv. Mater. 27, 6371 (2015).
26. S. Liu, I. Grinberg,  and A. M. Rappe, Appl. Phys. Lett. 103, 232907 (2013c).
27. R. Xu, S. Liu, I. Grinberg, J. Karthik, A. R. Damodaran, A. M. Rappe,  and L. W. Martin, Nat. Mater. 14, 79 (2015).
28. G. Kwei, A. Lawson, S. Billinge,  and S. Cheong, J. Phys. Chem. 97, 2368 (1993).
29. O. Diéguez, S. Tinte, A. Antons, C. Bungaro, J. B. Neaton, K. M. Rabe,  and D. Vanderbilt, Phys. Rev. B 69, 212101 1 (2004).
30. A. Von Hippel, Rev. Mod. Phys. 22, 221 (1950).
31. I. D. Brown, Chem. Rev. 109, 6858 (2009).
32. I. Brown and R. Shannon, Acta Crystallogr. A 29, 266 (1973).
33. I. Brown and K. K. Wu, Acta Crystallogr. B 32, 1957 (1976).
34. M. Finnis and J. Sinclair, Philos. Mag. A 50, 45 (1984).
35. M. A. Harvey, S. Baggio,  and R. Baggio, Acta Crystallogr. B 62, 1038 (2006).
36. P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, et al., J. Phys.: Condens. Matter 21, 395502 (2009).
37. J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou,  and K. Burke, Phys. Rev. Lett. 100, 136406 (2008).
38. http://opium.sourceforge.net.
39. H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976).
40. M. Parrinello and A. Rahman, Phys. Rev. Lett. 45, 1196 (1980).
41. M. Sepliarsky, M. G. Stachiotti,  and R. L. Migoni, Phys. Rev. B 72, 014110 (2005).
42. R. E. Cohen and H. Krakauer, Phys. Rev. B 42, 6416 (1990).
43. M. Gaudon, Polyhedron 88, 6 (2015).
44. K. H. Ehses, H. Bock,  and K. Fischer, Ferroelectrics 37, 507 (1981).
45. C. Jun, F. Chan-Gao, L. Qi,  and F. Duan, J. Phys. C: Solid State Phys. 21, 2255 (1988).
46. R. Comes, M. Lambert,  and A. Guinier, Solid State Commun. 6, 715 (1968).
47. K. Itoh, L. Zeng, E. Nakamura,  and N. Mishima, Ferroelectrics 63, 29 (1985).
48. E. A. Stern, Phys. Rev. Lett. 93, 037601 (2004).
49. K. J. Gaffney and H. N. Chapman, Science 316, 1444 (2007).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters