DFT calculations of strain effects on perovskite ferroelectric films

1. Introduction

The property of ferroelectricity was first observed in 1921 in Rochelle salt. Due to their ability to sustain a macroscopic polarization that can be switched by the action of an electric field, ferroelectric (FE) materials have been of technological importance in microelectronic devices and computer memories.

The structure of a perovskite FE contains cation A at the cube corners, a cation B at the center of the cube, and oxygen atoms at the face center of the cube forming a regular octahedron. Typically, perovskites have paraelectric phase at high temperature. As temperature decreases and goes through the phase transition, other phases with lower symmetry begin to form.

High temperature
.
Low temperature

The properties of perovskite ferroelectric thin films are strongly influenced by the magnitude of the epitaxial strain arising from the lattice parameter mismatch between the film and the substrate, the thermal strains becuase of the difference between the thermal expansion coefficients of the film and the substrate, the self-strain during phase transformations, and defects.

Therefore, mapping out temperature-versus-misfit strain phase diagram by calculation has critical significance to predict and interpret the experiment results as well as to find some potetial new phenomena. In the following section, taking BaTiO₃ as an example, the phase diagrams obtained from first principles calculation will be discussed.

2. First principles calculations

2.1. Full first-principles calculations

In the work done by Dieguez etc. [1], they mapped out the equilibrium structure of BaTiO₃ as a function of epitaxial strains by using parameter-free total-energy methods based on density functional theory (DFT).

The DFT calculations were carried out in the Kohn-Sham framework using the VASP software package. The electron-ion interaction was described by the projector augmented wave (PAW) method; semicore electrons are included for Ba semicore electrons are included in the case of Ba ($5s^{2}5p^{6}6s^{2}$) and Ti ($3s^{2}3p^{6}4s^{2}3d^{2})$). Also, they used the Ceperley–Alder form of the localdensity approximation (LDA) exchange-correlation functional, a 700 eV plane-wave cutoff, and a 6$\times$6$\times$6 Monkhorst–Pack sampling of the Brillouin zone. The optimization is performed in the five-atom unit cell of the cubic perovskite structure (space group Pm$\overline{3}$m) for six possible phases shown below:

At each value of misfit, the DFT calculation starts from a structure of fixed symmetry, then relax the atomic positions and the out-of-plane cell vector by displacing the Ti and O atoms until the Helmann-Feynman forces and the stress tensor components which corresponding to the thin-film condition are below some given thresholds. Repeating these steps for all the six possible phases, they obtained the total energy as a function of misfit strain for the six phases, as shown below. The misfit strain is defined as $(a-a_{0})/a_{0}$, where $a_{0}$ is theoretical in-plane lattice constant for free cubic BaTiO₃, and a is the chosen in-plane lattice constant for the optimization.

The DFT calculation shows that at zero temperature, the sequence of the most stable phase in the BaTiO₃ epitaxial thin film from large compressive strains to large tensile strains is: c$\rightarrow$r$\rightarrow$aa.

2.2. First-principles-based KSV calculations: a Landau-Devonshire model

2.2.1 Formalism
In Landau-Devonshire theory, the thermodynamic potential of a FE material is expressed by a polynomial in powers of polarization vectors. It is also possible to develop models with the same spirit, but constructed from first-principles calculations. According to King-Smith and Vanderbilt[2], the total energy can be expressed as a Taylor expansion around the cubic perovskite structure in terms of six independent components of the strain tensor and the three Cartesian components.

For bulk , the total energy per unit cell can be expressed as [3]:

(1)

\begin{align} \emph{E}=\emph{E}^{elas}(\{\eta_{i}\})+\emph{E}^{soft}(\{\emph{u}_{\alpha}\})+\emph{E}^{int}(\{\eta_{i}\},\{\emph{u}_{\alpha}\}) \end{align}

where E^elas, E^soft, E^int represents the energy purely from strain, purely from soft-mode amplitude, and the interaction between these two kinds of degrees of freedom. For crystals with cubic symmetry, the three energy components can be expressed as:

(2)

\begin{align} \emph{E}^{elas}(\{\eta_{i}\})=\frac{1}{2}\emph{B}_{11}(\eta_{1}^{2}+\eta_{2}^{2}+\eta_{3}^{2})+\emph{B}_{12}(\eta_{1}\eta_{2}+\eta_{2}\eta_{3}+\eta_{3}\eta_{1})+{\frac{1}{2}}\emph{B}_{44}(\eta_{4}^{2}+\eta_{5}^{2}+\eta_{6}^{2}) \end{align}

(3)

\begin{align} \emph{E}^{soft}(\{\emph{u}_{\alpha}\})=\kappa\emph{u}^{2}+\alpha\emph{u}^{4}+\gamma(\emph{u}_{x}^{2}\emph{u}_{y}^{2}+\emph{u}_{y}^{2}\emph{u}_{z}^{2}+\emph{u}_{z}^{2}\emph{u}_{x}^{2}) \end{align}

(4)

\begin{align} \emph{E}^{int}(\{\eta_{i}\},\{\emph{u}_{\alpha}\})=\frac{1}{2}\emph{B}_{1xx}(\eta_{1}\emph{u}_{x}^{2}+\eta_{2}\emph{u}_{y}^{2}+\eta_{3}\emph{u}_{z}^{2})+\emph{B}_{1yy}[\eta_{1}(\emph{u}_{y}^{2} +\emph{u}_{z}^{2})+\eta_{2}(\emph{u}_{z}^{2}+\emph{u}_{x}^{2})+\eta_{3}(\emph{u}_{x}^{2}+\emph{u}_{y}^{2})]+\emph{B}_{4yz}(\eta_{4}\emph{u}_{y}\emph{u}_{z}+\eta_{5}\emph{u}_{z}\emph{u}_{x}+\eta_{6}\emph{u}_{x}\emph{u}_{y}) \end{align}

In the case of epitaxial thin film, the epitaxial strain constraint imposed by the substrate is [4]:

(5)

\begin{align} \eta_{1}=\eta_{2}=\overline{\eta} \end{align}

(6)

\begin{align} \eta_{6}=0 \end{align}

Considering no shear stress, the energy expression for epitaxial thin film is given by:

(7)

\begin{align} \widetilde{\emph{G}}=\emph{E}-\sigma_{3}\eta_{3} \end{align}

By minimizing the enthalpy with respect to $\eta_{3}$, $\eta_{4}$, and $\eta_{5}$ to obtain the expression of the three variables as a function of $\overline{\eta}$, $\emph{u}_{x}$, $\emph{u}_{y}$, and $\emph{u}_{z}$, substituting the expressions into Equation above, and then make arrangements, $\widetilde{\emph{G}}$ could be expressed as:

(8)

\begin{align} \widetilde{\emph{G}}=(\emph{A}_{\overline{\eta}\overline{\eta}}\overline{\eta}^{2}+\emph{A}_{\overline{\eta}\sigma}\overline{\eta}\sigma+\emph{A}_{\overline{\sigma}\overline{\sigma}}\overline{\sigma}^{2})+(\emph{B}_{\overline{\eta}}\overline{\eta}+\emph{B}_{\sigma}\sigma+\emph{B})\emph{u}_{xy}^{2}+(\emph{C}_{\overline{\eta}\overline{\eta}}\overline{\eta}+\emph{C}_{\overline{\sigma}}\sigma+\emph{C})\emph{u}_{z}^{2} +\emph{D}\emph{u}_{xy}^{4}+\emph{E}\emph{u}_{z}^{4}+\emph{F}\emph{u}_{xy}^{2}\emph{u}_{z}^{2}+\emph{H}\emph{u}_{xy}^{4}\sin^{2}\theta\cos^{2}\theta \end{align}

The coefficients appearing in the energy expression for thin film are all related to the coefficients appearing in that for bulk. And all these coefficients can be obtained from first-principles calculations.

2.2.2 First-principles calculations of the coefficients

In their paper[2], King-Smith and Vanderbilt showed that the minimum total energy as a function of soft-mode displacement can be determined correct to fourth order in $\emph{u}_{\alpha}$, which is the soft-mode magnitude, in terms of nine independent interaction parameters $\emph{B}_{11}$, $\emph{B}_{12}$, $\emph{B}_{44}$, $\emph{B}_{1xx}$, $\emph{B}_{1yy}$, $\emph{B}_{4yz}$, $\kappa$, $\alpha$ and $\gamma$. Based on the Vanderbilt ultrasoft-pseudopotential scheme, these parameters can be calculated by using a conjugate-gradient technique for minimizing the Kohn-Sham energy functional. More details are shown in their paper[2]. Deiguez et al. repeated the calculations, and pushed the boundaries of the numerical approximations that control the accuracy of the first-principles calculations [4]. Their calculated results are shown below.

2.2.3 Testing this model: the case of zero perpendicular stress

Figures below show the energy curves as a function of strain for six phases calculated by the first-principles-based KSV theory (right) compared with the full DFT results (left).

The agreement between the two figures is very good, except small differences that may come from the different DFT calculation methods and the intrinsic errors of the Taylor expansion of the energy in Landau-Devonshire model.

2.2.4 Finite-temperature calculations

The first-principles-based calculations can be also applied to finite temperature calculations. But compared with the zero-temperature condition, the total energy of the system at finite temperature would include more terms coming from the coupling between short-range and long-range local modes.

The calculations include three steps [5]:
(1) Construct an effective Hamiltonian to describe the important degrees of freedom of the system.

The central quantity for studying the equilibrium properties of a system at finite temperature is its partition function [5]. To construct an appropriate parameterized Hamiltonian, two fundamental approximations had been used:

(a) The first approximation is to use an energy surface represented by a low-order Taylor expansion. In the calculation of ref. [5], the expansion is up to fourth-order terms.
(b) The second approximation is to express the energy surface only as a function of the soft-mode amplitudes and strain. This approximation helps reducing the number of degrees of freedom per cell from 15 to 6.

With these approximations, the Hamiltonian at finite temperature consists of five parts: a local-mode self-energy, a long-range dipole-dipole interaction, a short-range interaction between soft modes, and elastic energy, and an interaction between the local modes and local strain.

(2) Determine all the parameters of the effective Hamiltonian from high-accuracy ab initio LDA calculations

The parameters in the effective Hamiltonian have already determined in the zero-temperature calculation. Therefore, the total-energy functional of the system is fully specified by a set of parameters which can be obtained from first-principles calculations. By using the Vanderbilt ultrasoft pseudopotential and a preconditioned conjugate gradient method which was also mentioned above, a generalized Kohn-Sham functional is directly minimized.

(3) Carry out Monte Carlo (MC) simulation to map out the finite-temperature phase diagram.

The constraint of fixed in-plane strain was imposed by fixing three of the six elements of the homogeneous strain tensor during the MC simulations. For each value of in-plane strain, MC thermal averages are obtained for the unconstrained components of the homogeneous strain and the average polarization, and phase transitions are identified by monitoring the symmetries of these quantities. Finally the phase diagram can be obtained, and showed below.

Compared with analytical approaches, MC simulation requires only the ability to compute changes in total-energy as the configuration is changed. With suitable analysis of statistical errors and finite-size effects, the results of MC simulation can be made arbitrarily accurate.

2.3 Comparison of these two calculation methods [3, 4]

The full first-principles calculations is straightforward to get accurate results for arbitrary misfit strains at T=0. The challenge for this approach is in predictions of temperature-driven behavior. Also, the quantitative accuracy of predicted transition temperatures is relatively low when compared with experimental phase diagram.

In comparison, the first-principles-based KSV calculations, the information needed for the model fit can readily be calculated consistently and accurately within the framework of the chosen first-principles approach. Once the fitting has been completed, such models can be used to compute detailed properties much faster than the full first-principles calculations. Thus , they can be applied to larger systems where the full DFT calculations would be impractical.

3. Future work

For real films, some aspects are important to be introduced in the DFT calculations[6]:

(1) The interaction of epitaxial strain with other factors influencing the films, like free surfaces or interfaces with the substrate to produce the observed behavior of real films.
(2) Extension of first-principles approaches to multi-domain films.
(3) Exploration of the epitaxial strain phase diagrams in more materials with potential desirable properties or new phenomena.

4. References

[1] Oswaldo Dieguez et al. Phys. Rev. B 69 (2004) 212101
[2] R. D. King-Smith and David Vanderbilt Phys. Rev. B 49 (1994) 5828
[3] Oswaldo Dieguez and David Vanderbilt Phase Transitions 81 (2008) 607
[4] Oswaldo Dieguez et al. Phys. Rev. B 72 (2005) 144101
[5] W. Zhong et al. Phys. Rev. B 52 (1995) 6301
[6] Karin M. Rabe Curr. Opin. Soid State Mater. Sci. 9 (2005) 122