Polarization branches and optimization calculation strategy applied to ABO3 ferroelectrics

TL;DR

This paper examines Berry phase polarization calculations in ferroelectrics, highlighting ambiguities and proposing an optimized strategy to reduce computational effort and obtain clear polarization results.

Contribution

It introduces a new procedure to minimize calculations and resolve ambiguities in Berry phase polarization analysis of ABO3 ferroelectrics.

Findings

Identified multivalued polarization issues in calculations.

Proposed an optimization strategy to clarify polarization results.

Reduced computational complexity for large super-cell analyses.

Abstract

Berry phase polarization calculations have been investigated for several ferroelectric materials from the point of view of practical calculations. It was shown that interpretation of the results is particular to each case due to the multivalued aspect of polarization in the modern theory. Almost all of the studied examples show ambiguous polarization results which can be difficult to solve especially for super-cells containing large number of atoms. For this reason, a procedure has been proposed to minimize the number of calculations required to produce an unambiguous polarization result from Berry phase polarization investigations.

Figures

Tables

Table 1. Table 1: Lattice parameters and atomic coordinates of the optimized P ​ b ​ T ​ i ​ O 3 𝑃 𝑏 𝑇 𝑖 subscript 𝑂 3 PbTiO_{3} , B ​ a ​ T ​ i ​ O 3 𝐵 𝑎 𝑇 𝑖 subscript 𝑂 3 BaTiO_{3} and K ​ N ​ b ​ O 3 𝐾 𝑁 𝑏 subscript 𝑂 3 KNbO_{3} unit-cells.

a (Å)		c (Å)	$z_A$	$z_B$	$z_{O_1}$	$z_{O_{2,3}}$
$BaTi{O}_3$	3.9925	4.0365	0.00	0.4785	0.0253	0.5105
$PbTi{O}_3$	3.8775	4.2070	0.00	0.5380	0.1166	0.6211
$KNb{O}_3$	3.9970	4.0630	0.00	0.4770	0.0170	0.519

Table 2. Table 2: Spontaneous polarization obtained using the proposed optimization strategy for all the studied materials shown in Figures 3 - 6 .

			$P_{CS}$ ($C{m}^{-2}$)	$d\lambda$	$P_{d\lambda }$ ($C{m}^{-2}$)	$n$	$\Lambda$ ($C{m}^{-2}$)	$P_{lin}\left(1\right)$ ($C{m}^{-2}$)	$P_{FS}$ ($C{m}^{-2}$)	$k$	$P_s$ ($C{m}^{-2}$)
DOWN	BTO	Ba	$0.000$	$5\%$	$-0.019$	$0$	$-0.0038$	$-0.380$	$-0.350$	$0$	-0.350
		Ti	$0.000$	$5\%$	$-0.019$	$0$	$-0.0038$	$-0.380$	$1.685$	$-1$
	PTO	Pb	$0.000$	$5\%$	$-0.050$	$0$	$-0.010$	$-1.000$	$1.185$	$-1$	$-0.945$
		Ti	$0.000$	$5\%$	$-0.050$	$0$	$-0.010$	$-1.000$	$1.185$	$-1$
	KNO	K	$-0.501$	$5\%$	$-0.522$	$0$	$-0.0042$	$-0.921$	$-0.870$	$0$	$-0.369$
		Nb	$-0.501$	$5\%$	$-0.522$	$0$	$-0.0042$	$-0.921$	$1.135$	$-2$
	PZT	None	$0.000$	$10\%$	$-0.093$	$0$	$-0.0093$	$-0.930$	$0.186$	$-2$	$-0.878$
UP	BTO	Ba	$0.000$	$5\%$	$0.019$	$0$	$0.0038$	$0.380$	$0.350$	$0$	$0.350$
		Ti	$0.000$	$5\%$	$0.019$	$0$	$0.0038$	$0.380$	$-1.685$	$1$
	PTO	Pb	$0.000$	$5\%$	$0.050$	$0$	$0.010$	$1.000$	$-1.185$	$1$	$0.945$
		Ti	$0.000$	$5\%$	$0.050$	$0$	$0.010$	$1.000$	$-1.185$	$1$
	KNO	K	$-0.501$	$5\%$	$0.522$	$-1$	$0.0042$	$-0.081$	$0.870$	$-1$	$0.369$
		Nb	$-0.501$	$5\%$	$0.522$	$-1$	$0.0042$	$-0.081$	$-1.134$	$1$
	PZT	None	$0.000$	$10\%$	$0.093$	$0$	$0.0093$	$0.93$	$-0.186$	$2$	$0.878$

Equations

$$z_i(\lambda )=z_i^{CS}+(z_i^{FS}-z_i^{CS})\lambda ,$$

$$P_s=P_{FS}-P_{CS},$$

$${\mathbf{P}}_{state}={\mathbf{P}}_{Berry}+n{\mathbf{P}}_{\mathbf{q}},n\in \mathbb{Z},$$

$$P_{lin}(\lambda )=\Lambda \lambda +P_{CS},$$

$$\Lambda =\frac{P_{d\lambda }-P_{CS}}{d\lambda },$$

$$\Lambda (n)=\frac{P_{d\lambda }+n{P}_q-P_{CS}}{d\lambda },n\in \mathbb{Z}.$$

$$k=\left[\frac{P_{lin}(\lambda =100\%)-P_{FS}}{P_q}\right],k\in \mathbb{Z},$$

$$P_s=(P_{FS}+k{P}_q)-P_{CS}.$$

$$Z_X^{\ast }=\frac{\Omega }{e}\frac{P(dr)-P(0)}{dr},$$

$$P=\frac{e}{\Omega }\sum _i{Z}_i^{\ast }{z}_i,$$

Full text

††preprint: APS/123-QED

I Introduction

Ferroelectric materials, have been under intense research for the better part of the last three decades due to the tremendous integration potential into applications ranging from high-density non-volatile memories to solar cells and logic gates scott1989 ; hidaka1992 ; glinchuk2009 ; kumari2015 ; rong2016 ; scott2007 ; choi2009 ; garcia2009 ; huang2010 ; yuan2011 ; liu2014 ; chen2015 ; kim2016 . Most ferroelectrics belong to the perovskite $ABO_{3}$ family which is a versatile group of materials that can be obtained through a wide range of synthesis methods moure2015 . One important property for real life applications is the large polarization value obtained in some ferroelectric compounds such as $Pb(Zr_{x}Ti_{1-x})O_{3}$ (PZT) thin films (around $1\text{\,}\mathrm{C}\text{\,}{\mathrm{m}}^{-2}$) pintilie2007 . Unfortunately these high polarization values are mostly obtained in lead-containing materials which have a rather large environmental footprint. For this reason, there is a sustained experimental and theoretical effort to find lead-free perovskites that ideally retain the high polarization value but have a much lower toxicity. The theoretical approach for this search is based on high throughput calculations using automated scripts searching for particular material properties (such as a high polarization value) kundu2015 ; volonakis2016 ; filip2016 ; filip2016a . The most reliable theoretical model to date for computing the spontaneous polarization is an implementation of the Berry phase (BP) formalism berry1984 adapted by Resta resta1992 in 1992 to correctly define and calculate the bulk polarization value of a ferroelectric. It was later implemented by King-Smith and Vanderbilt kingsmith1993 in first-principles density-functional theory numerical calculations routine and together have become what is now called the modern theory of polarization. This method has been used to compute polarization values for various bulk ferroelectric materials. However, special attention is required for each studied material. The difficulty lies behind the definition of the spontaneous polarization within the modern theory of polarization as the time integral of the electric current appearing when the system is distorted adiabatically from a reference centrosymmetric state (CS) towards a final ferroelectric state (FS) resta1992 ; resta1993a ; resta1993b ; resta1994 ; rabe2007 . Therefore, in order to calculate the spontaneous polarization, one needs the polarization values for at least two different system states (CS and FS). Henceforth the spontaneous polarization is the difference between the two. Therefore such an approach is fundamentally different than previously used models where the bulk polarization was viewed as a collection of neatly arranged electrical dipoles rabe2007 . In addition, the calculated polarization is not single valued but a multivalued function of the system state; it can only be obtained up to an integer number (indexing polarization branches) multiplied by a certain geometrical system constant resta1992 ; resta1993a ; resta1994 ; rabe2007 . This intricacy complicates the interpretation of numerical results considerably, since the obtained polarization values for the required CS and FS may not have the same indeterminacy (i.e. they may belong to different polarization branches). As such the difference in polarization of the CS and FS is multivalued contradicting experimental spontaneous polarization measurements. This ambiguity is corrected by repeating calculations for a large number of intermediary system states, so as to identify the polarization branches of the CS and FS. However, the procedure is system dependent and can also lead to incorrect results if the proper corrections are not used, as pointed out by Neaton et. al neaton2005 .

It is at this point where a distinction between the analytical theoretical description of the modern theory of polarization and its numerical implementation and usage, should be made. While the BP theory of polarization is general to all materials, the interpretation of results obtained from the numerical implementations are strongly system dependent. This fact implies that almost always one must resort to calculating the polarization for multiple intermediary system states which can be very time consuming especially for larger systems.

The present paper aims to illustrate the BP polarization theory from the perspective of practical density functional theory hohenberg1964 calculations and to introduce a strategy for choosing the optimum calculation points required to resolve the polarization ambiguity in a minimum number of steps. The manuscript is structured as follows: in the next section, the materials used in this study are presented, followed by the computational details in Section III. Berry phase calculations are briefly described in Section IV and the polarization calculations results for the materials introduced in Section II are individually discussed in Section V. The proposed optimization strategy is introduced in Section VI using the example of a larger system. The manuscript will be concluded with a discussion of all the results obtained in this study and the implemented strategy in Section VII.

II Materials and structures

The application of the BP polarization theory can be better understood if practical calculations are discussed. For this reason, in this study, four materials have been used to exemplify various aspects of polarization calculations: $BaTiO_{3}$ (BTO), $PbTiO_{3}$ (PTO), $KNbO_{3}$ (KNO) and $Pb(Zr_{0.25}Ti_{0.75})O_{3}$ (PZT). The choice for the first three materials has been motivated by the simplicity of their respective unit-cells which allowed for a large number of calculations to be performed for intermediary system distortion states. The symmetry of the systems is tetragonal and they each have five atoms in the unit-cell (see Figure 1). The lattice parameters and atomic coordinates for the first three materials have been summarized in Table 1.

$Pb(Zr_{0.25}Ti_{0.75})O_{3}$ has been used to exemplify BP polarization calculations for a larger system. The 40 atom super-cell was obtained by arranging 8 PTO unit-cells in a $2\times 2\times 2$ grid and replacing two $Ti$ atoms with $Zr$. Being a larger system the polarization calculations took considerably longer time to complete compared to calculations for the smaller studied systems. While the relative time it takes for a polarization calculation to finish is not relevant when individual materials are tested, it becomes important when multiple such calculations are performed for an entire range of different materials.

III Computational details

All calculation for the present study were performed using the generalised gradient approximation (GGA) density functional theory hohenberg1964 as implemented in the Quantum Espresso suite giannozzi2009 . Projected augmented-wave (PAW) Perdew-Burke-Ernzerhof pseudo-potentials optimized for solids (PBEsol) perdew2008 , from the THEOS library theos , have been used for all the numerical results obtained in this study. However, it should be noted that while small differences can be obtained for the value of the spontaneous polarization, the ambiguities generated by the multivalued aspect remain. For this reason the precision of the numerical results in this paper is less important since they were used for visualisation of the peculiarities involved with the Berry phase polarization calculations.

For the first three materials, a kinetic energy cut-off of $100\text{\,}\mathrm{R}\mathrm{y}$ was chosen and the Brillouin zone was sampled in a $5\times 5\times 5$ Monkhorst-Pack monkhorst1976 uniform k-point grid for all self-consistent calculations. The polarization calculations (non-self-consistent) were performed with a $5\times 5\times 20$ k-point grid. For the BTO and PTO materials, structural optimization was performed starting from experimental data found in the literature wang2010 ; saghiszabo1998 . The optimized lattice parameters were obtained by minimizing the total energy with respect to the volume of the unit cell and the final atomic coordinates were obtained by relaxing the internal coordinates until the Hellmann-Feynman forces are converged to $2.6\text{\times}{10}^{-9}\text{\,}\mathrm{eV}\text{\,}{\mathrm{\SIUnitSymbolAngstrom}}^{-1}$. The KNO case has been treated using the same lattice parameters and atomic coordinates as obtained by Dall’Olio et. al dallolio1997 .

The $Pb(Zr_{0.25}Ti_{0.75})O_{3}$ material was also investigated as an example of a larger system described by a super-cell formed with 8 PTO unit-cells disposed in a $2\times 2\times 2$ grid where 2 $Ti$ atoms have been replaced by $Zr$. The lattice parameters of the PZT super-cell are twice the ones obtained for PTO (see Table 1). The atomic coordinates were further relaxed to account for the presence of the two $Zr$ atoms. Since the larger system will take considerably longer time for each polarization calculation run to complete, it was used to highlight the consequences of randomly choosing intermediary system distortion states between CS and FS in order to solve the polarization ambiguity. Lastly, as a language convention, a polarization calculation run is a set of two consecutive calculations for a given system state: first a self-consistent calculation on a uniform k-point grid, followed by a non-self-consistent calculation with an increased number of k points along the axis where the polarization is calculated.

IV Berry phase primer and usage

The modern theory of polarization defines ferroelectric polarization as the time integral of the current that appears through the sample when the studied system is adiabatically distorted from a reference state to a final ferroelectric state. This is valid as long as the system remains insulating in all the intermediate states along the distortion path rabe2007 . A visual representation of two such processes are given in Figure 1 where the tetragonal $ABO_{3}$ system is distorted from the CS to the FS along the z axis with respect to the A and B sublattices, respectively. For all materials studied in this paper, the polarization direction is parallel to the z axis of the unit-cell and thus the systems will only be distorted along this axis. For the cases shown in Figure 1 we can write the system distortion as a linear function of a dimensionless parameter $\lambda$:

[TABLE]

where the index $i$ spans all the atoms in the unit cell and $\lambda\in[0,1]$. The coordinates $z^{CS}_{i}$ and $z^{FS}_{i}$ represent the $z$ coordinates of atom $i$ for the centrosymmetric and the ferroelectric state, respectively, while $z_{i}(\lambda)$ is the coordinate of the same atom $i$ for an intermediary distortion $\lambda$. This means that all the atoms in the unit-cell will be moved together by a fraction $\lambda$ of their corresponding final displacements $\Delta z_{i}=z^{FS}_{i}-z^{CS}_{i}$. Following the definition of polarization in the modern theory, the difference between the values obtained for the FS and CS should represent the spontaneous polarization of the studied system.

[TABLE]

where, $P_{s}$ is the measurable spontaneous polarization, $P_{FS}$ the computed polarization in the ferroelectric state and $P_{CS}$ the computed polarization in the reference state (in our case, the centrosymmetric state). This is a measurable quantity and the modern theory mimics the generic experimental measurement method, where the system state is switched between the two stable ferroelectric states using an external electric field along a hysteresis cycle.

However, since the polarization is a multivalued function, Eq. (2) does not strictly represent the spontaneous polarization rabe2007 ; resta1993a . Indeed it can be shown that polarization is only well-defined modulo a polarization quantum given by: $P_{q}=\frac{e\,\bm{R}}{\Omega}$, where $\Omega$ is the unit-cell volume, $e$ is the unit charge and $\bf{R}$ is any lattice vector. This means that the calculated polarization for any system state $P_{state}$ is actually an entire family of values separated by integer multiples of $P_{q}$ (also called branches) and given by resta1993a ; rabe2007 :

[TABLE]

where, $P_{Berry}$ is the calculated polarization and $n$ is an integer indexing the polarization branch resta1993a ; rabe2007 . It is now clear that, if the polarization values for the CS and FS do not belong to the same branch, then the spontaneous polarization will continue to have the form given in Eq. (3), contradicting experimental findings.

Eq. (3) is not only an important consequence of the modern approach to ferroelectric polarization but also the key to the correct application of the numerical implementation of this theory for practical cases. From the theoretical point of view, one can identify if the polarization values for the CS and FS belong to the same branch by looking at their difference and comparing it to the polarization quantum resta1994 ; kingsmith1993 ; rabe2007 . If the difference is much smaller than the polarization quantum then no ambiguity has appeared and Eq. (2) is valid. For the case when the difference is comparable to the polarization quantum, intermediary system distortions should be considered in order to clarify what branch the two polarization values belong to and which corrections should be made. However, there is no clear indication what it means for the polarization difference to be “much smaller” and it will be shown that this comparison cannot be used effectively to solve the ambiguity. Figure 2 summarizes the calculation steps in order to obtain the spontaneous polarization using the Berry phase polarization method.

V Numerical calculations and polarization branches

$BaTiO_{3}$

The first choice of distortion path shown in Figure 3a) returns the ideal result: All the calculated polarization values belong to a single polarization branch. It is then clear that by applying Eq. (2) with the FS and CS values one obtains the correct polarization value of $0.351\text{\,}\mathrm{C}\text{\,}{\mathrm{m}}^{-2}$. Using the checkpoints at the end of Section IV, the values for the “UP” directions are as follows: $P_{CS}=$0\text{,}\mathrm{C}\text{,}{\mathrm{m}}^{-2}, $P_{FS}=$0.351\text{\,}\mathrm{C}\text{\,}{\mathrm{m}}^{-2} and $P_{q}=$2.036\text{,}\mathrm{C}\text{,}{\mathrm{m}}^{-2}. The difference between the FS and CS values is almost 6 times smaller than the polarization quantum and this could have been considered sufficient to conclude that an unambiguous result can be obtained without any other intermediary system distortions. However, the situation changes drastically if the unit-cell of the system is now distorted with respect to the B atom site (see Figure [3](#S5.F3)b)). From the CS to around $80\%$ distortion, the calculated polarizations are neatly arranged on the same branch just like in the previous example. Nevertheless for the rest of the $20\%$ left to the final FS the values jump suddenly on a different branch. Analysing the end values of interest, one obtains: $P_{CS}=$0\text{\,}\mathrm{C}\text{\,}{\mathrm{m}}^{-2} and $P_{FS}=$-1.685\text{,}\mathrm{C}\text{,}{\mathrm{m}}^{-2}, while the polarization quantum remains the same, $P_{q}=$2.036\text{\,}\mathrm{C}\text{\,}{\mathrm{m}}^{-2}.

According to the calculations steps outlined in Figure 2 the difference is now $-1.685\text{\,}\mathrm{C}\text{\,}{\mathrm{m}}^{-2}$ which is almost the same as $P_{q}$ in absolute value. This means that an ambiguity has arisen and a correction is needed. In this case, with that many calculated values, it is easy not only to identify where the calculated polarization values reside on their corresponding branch but also what correction should be used. For this case, by adding one polarization quanta to $P_{FS}$ will bring its value on the same branch as $P_{CS}$ (for the “UP” direction). The spontaneous polarization in Eq. (2) becomes: $P_{s}=P_{FS}+P_{q}-P_{CS}=$0.351\text{,}\mathrm{C}\text{,}{\mathrm{m}}^{-2}$$.

$PbTiO_{3}$

For the first distortion path in Figure 4a) more polarization jumps appear not only towards the end of the distortion interval. This only reflects the “random” nature of these polarization jumps that are driven by the numerical implementation of the BP polarization theory. In order to verify that these jumps are not related to any numerical errors, the total energy is plotted (dashed green line) as a function of the system distortion. If the irregularities in the dependence of polarization on the system distortion had a different source, a corresponding discontinuity in the energy dependence should be observed.

This is clearly not the case in any of the studied examples. Apart from the appearance of irregular polarization jumps, the important FS and CS polarization are located on different branches similar to the case in Figure 3b): $P_{CS}=$0\text{,}\mathrm{C}\text{,}{\mathrm{m}}^{-2}, $P_{FS}=$-1.185\,210\,2\text{\,}\mathrm{C}\text{\,}{\mathrm{m}}^{-2} and $P_{q}=$2.13\text{,}\mathrm{C}\text{,}{\mathrm{m}}^{-2}. The spontaneous polarization can be obtained in a similar fashion, by adding one polarization quanta (for the “UP” direction): $P_{s}=P_{FS}+P_{q}-P_{CS}=$0.944\text{\,}\mathrm{C}\text{\,}{\mathrm{m}}^{-2}. It should be noted that the value of the obtained spontaneous polarization in this case is not “much smaller” than the polarization quanta and this can add more confusion when analysing calculated data. For the second distortion path in Figure 4b), the irregular jumps from the previous case have disappeared but the polarization ambiguity is solved just like in the previous case. In fact, if it was not for the one instance in Figure 3a), one could be inclined to conclude at this point that, by always adding one polarization quanta to $P_{FS}$ for the “UP” direction, Eq. (2) would lead to an unambiguous result. Unfortunately this is not true, as the next calculated example will show.

$KNbO_{3}$

Figure 5 shows the same type of polarization calculations as the previous two cases, however the results are entirely different, starting with the CS polarization value. Contrary to the results in Figures 3 and 4 the CS polarization in Figure 5 does not vanish for the centrosymmetric state. This result appears to be in contradiction with the expectation that for a centrosymmetric state of a system the polarization must vanish, yet this is only true if the polarization was single valued! rabe2007 ; spaldin2012 . Once again it is shown that using the current approach it is impossible to successfully apply the same strategy to different materials. Another difference in the results obtained for KNO is the polarization jump that appears in the near vicinity of the CS. This only appears for one polarization direction but not the other. This is a fortuitous result because it reinforces the fact that in order for Eq. (2) to be applied, both the CS and the FS polarizations must belong to the same branch. For this case the CS and FS values for the “UP” direction in Figure 5a) were: $P_{CS}=$-0.501\text{,}\mathrm{C}\text{,}{\mathrm{m}}^{-2}, $P_{FS}=$0.870\text{\,}\mathrm{C}\text{\,}{\mathrm{m}}^{-2} and $P_{q}=$1.002\text{,}\mathrm{C}\text{,}{\mathrm{m}}^{-2}. By correcting either value, the spontaneous polarization can be obtained to be $P_{s}=$0.369\text{\,}\mathrm{C}\text{\,}{\mathrm{m}}^{-2}.

VI Optimization strategy for Berry phase polarization calculations

The examples shown in Figures 3-5 have been presented using a large number of calculations in order to help visualise polarization branches. In a regular scenario, some cases can be clarified using a smaller number of calculations. On average, depending on the choice of the distortion parameter value $\lambda$, four calculations should be enough to correctly identify the polarization branches and obtain the spontaneous polarization for the first three cases presented in Section V. Yet the only common feature that remains is that each material requires an individual treatment to obtain the spontaneous polarization. This process is far more tedious for larger system sizes. In order to illustrate our optimized method we show the case of PZT.

The polarization calculations for PZT are shown in Figure 6a)-d) for an ever increasing number of intermediary system distortion states. Due to the site and the large number of atoms in the super-cell it becomes difficult to fix a sublattice and perform the distortion with respect to it. As an alternative, in this case, all atoms in the unit cell have been moved as described in Eq. (1) from the CS to the corresponding FS. The calculated polarization values for the CS and FS are: $P_{CS}=$0.0\text{,}\mathrm{C}\text{,}{\mathrm{m}}^{-2}, $P_{FS}=$-0.186\text{\,}\mathrm{C}\text{\,}{\mathrm{m}}^{-2} and the polarization quanta $P_{q}=$0.532\text{,}\mathrm{C}\text{,}{\mathrm{m}}^{-2}$$. This seems to be an ambiguous situation and generally for an unknown material it is difficult to know for sure just from the difference between the FS and CS polarization values. For this reason, an extra point has been introduced in Figure 6a) corresponding to a $50\%$ system distortion. The purpose for adding more intermediary states is to obtain a clear image of the polarization branches in order to identify what correction is needed for the values of interest. In this case, the extra point does not help to clarify the situation completely since no polarization branches are revealed and we make another choice for the next system distortion. This time, a $25\%$ distortion is added to the plot in Figure 6b) and the result does not bring the conclusion any closer. Moving on, a $75\%$ distortion is added next in Figure 6c) with the same result. It is only when smaller distortion steps are added in Figure 6d) that a real image of the PZT polarization branches emerges and a correction can be made by adding two polarization quanta to the FS value in order to bring it on the same branch as the CS case.

Figure 6d) provides the best clue to introducing an optimization procedure. The image that can be formed from all the results presented in this study is that the polarization values belonging to the same branch will have a monotonous increase in absolute value starting from the reference state towards the ferroelectric one. This result has been discussed at length by Resta et. al resta1993a ; rabe2007 and it is known as the linear response theory, which states that the variation of polarization on the internal distortion should be linear. This has been shown to be true for KNO resta1993a and it can be observed to hold for the rest of the examples in this study, however it does not apply to all materials in general. Similar polarization studies on $BiFeO_{3}$, reveal a deviation from the linear dependence beyond $\approx 20\%$ distortion neaton2005 ; rabe2007 ; spaldin2012 . Nevertheless, in most cases one could safely use the following approximation for the polarization dependence on the unit-cell distortion parameter:

[TABLE]

where, $P_{CS}$ is the calculated polarization of the reference state of the studied system, which in the examples of this study was the centrosymmetric state and $\Lambda$ is the slope of the linear response. In order to be able to use this approximation the slope parameter $\Lambda$ must be computed through a “finite difference” route:

[TABLE]

where $P_{d\lambda}$ is the calculated polarization for an “infinitesimal” system distortion $d\lambda$ in the vicinity of the reference CS. However, from the definition of polarization in the BP theory, both polarization values in Eq. (5) are determined up to an integer value of polarization quanta.

Therefore the difference is still multivalued and, by extension, the slope parameter has a similar behaviour. The problem can be solved if the entire family of values for $P_{d\lambda}$ is introduced in Eq. (5), while the polarization for the reference state is kept at the calculated value $P_{CS}$ and its respective branch:

[TABLE]

It is now possible to find the integer number $n_{min}$ such that $\left|\Lambda(n_{min})\right|$ in Eq. (6) is minimum. This last operation essentially ensures that the calculated polarization for the “infinitesimal” distortion $d\lambda$ is brought to the same branch as the polarization for the reference CS. This can be easily verified by calculating several values for $\Lambda(n)$ and then using Eq. (4) to plot the corresponding dependencies. Figure 7 shows the linear polarization dependence for different values of the branch index $n$ corresponding to the calculated $P_{d\lambda}$ polarization. The continuous green line is obtained for the minimum $\left|\Lambda(n_{min})\right|$, where $n_{min}=0$. Using $n_{min}$, Eq. (6) returns the slope of the linear dependence of polarization independent of branch! With this information, Eq. (4) can be used to extrapolate the polarization value at $100\%$ distortion in the FS and subsequently obtain the integer number $k$ of quanta to correct the branch ambiguity:

[TABLE]

where, $P_{FS}$ is the calculated polarization. Finally, the spontaneous polarization can be obtained using Eq. (2) and the corrected $P_{FS}$ value:

[TABLE]

All the steps of the proposed optimization procedure have been summarized in a simple work flow in Figure 8. For materials such as the ones presented in this study, where the polarization direction is parallel to one of the principal axes, the total number of calculations is reduced to only three system states (steps 1, 3 and 7). The rest of the steps can be performed with simple arithmetic operations that can be easily included in a high throughput calculation automated script. For a more general case, the ferroelectric polarization can take any direction in the crystal and the operations proposed by this strategy must be repeated for each of the three cartesian axes.

VII Results

The most important result obtained in this study is the optimization strategy summarized in Figure 8. By comparison with the original method for Berry phase calculations in Figure 2 the proposed approach does not contain any interrogations steps and algorithm branches. Also, the method is a series of nine consecutive steps and only three calculation runs. The original approach however, may need considerably more calculations if multiple polarization jumps occur for the chosen distortion path of the system. Analysing Eq. (5) more closely, one can recognise a similitude with the method for calculating Born effective charges spaldin2012 ; rabe2007 ; resta1993a . By definition, Born effective charges represent the change in polarization divided by how much an ion is displaced spaldin2012 ; rabe2007 ; resta1993a which translates to the first derivative of polarization with respect to atomic displacement while keeping the rest of the ions and external macroscopic fields fixed resta1993a . In a practical calculation, this definition is transformed into a finite difference formula. Considering one of the materials studied in this paper, these charges are calculated by starting with the CS and then displacing atom $X$ by a small amount (typically around $0.005\text{\,}\mathrm{\SIUnitSymbolAngstrom}$) along the polarization direction, the Berry phase method is used to calculate the polarization. The Born effective charge associated to atom $X$ for a displacement along polarization direction is obtained as:

[TABLE]

where $dr$ is the infinitesimal displacement of atom $X$, $\Omega$ is the volume of the unit-cell and $e$ is the electron charge. After computing all the Born effective charges for the atoms in the unit-cell then the following linear approximation for the spontaneous polarization can be written:

[TABLE]

where the sum is over all the atoms in the unit cell and $z_{i}$ is the displacement of each atom in the FS (for this case, the displacements are only along the $z$ axis). It can be observed that the polarization approximation in Eq. (4) and Eq. (5) follows a similar approach, but with the infinitesimal simultaneous displacement of all atoms.

Using the steps summarized in Figure 8 one can obtain the spontaneous polarization for all the cases presented in this study in Figures 3-6 using only three of the calculated points for each distortion path. The results can be analysed in Table 2 where the values obtained for each step are shown separated between the two polarization directions (formally denoted as “UP” and “DOWN”). The third column of Table 2 indicates the sublattice with respect to which the systems were distorted. As described in Section III, for the PZT case all atoms were moved from their corresponding centrosymmetric positions in the super-cell toward the FS positions following the same rule in Equation (1).

The application of the proposed optimization strategy is straightforward yet there are some aspects that should be pointed out. First of all, the optimization strategy is based on the linear response theory of polarization as a function of the system distortion as discussed by Resta et. al resta1993a . This was shown to be true for the KNO case resta1993a and it was shown in this work that similar results are obtained for PTO and BTO. However, in a study by Neaton et. al on $BiFeO_{3}$ it is shown that the polarization dependence on the system distortion is only linear to about $20\%$ distortion. For situations where the materials have strong non-linear polarization dependence on the internal distortion, the proposed optimization strategy may need more calculated points in order to obtain a non-ambiguous result. Another important point to be considered, is one of the crucial conditions for the application of the BP theory: the system must be insulating in any state on the distortion path. This aspect is difficult to control,and such materials will require a more detailed investigation.

VIII Conclusions

In summary, the results presented in the current study provide a detailed view of the numerical implementation of the modern theory of polarization from the point of view of actual calculations on various ferroelectric materials. Most of the studies have been purposely performed with an exaggerated fine distortion mesh in order to provide a better visual image of the multivalued aspect of polarization and the emergence of polarization branches. We illustrate this using $BaTiO_{3}$, $PbTiO_{3}$ and $KNbO_{3}$ as test materials. It has been shown that using the approach proposed in the usual Berry phase polarization calculations, each material must be treated individually which makes it difficult to integrate such a study in an automated work-flow. For this reason, a unified optimization procedure has been proposed that can provide a starting point for polarization investigation for ferroelectric materials. The procedure uses a minimal number of calculations in order to obtain the spontaneous polarization value, thus reducing the computational effort. We hope the present study will enhance the current efforts in the theoretical investigation of known ferroelectrics and accelerate the design of new ferroelectric materials. The current study and the proposed calculation strategy provides both a visual representation of the multivalued aspect of polarization in the modern theory and a practical approach for such calculations which should complement the theoretical descriptions found in the specialized literature.

Acknowledgments

The authors would like to acknowledge the useful discussions with Dr. Claude Ederer from ETH, Zurich that helped clarify the theoretical and practical aspects of the Berry phase polarization theory.

LDF would like to acknowledge the financial support of the Romanian Ministry of Education Executive Unit for Funding High Education, Research, Development and Innovation (MEN-UEFISCDI) through the Young Research Team Grant PNII-RU-TE-2012-3-0320 (Contract No. 11).

The authors would also like to acknowledge the financial support from the NIMP Core Program nr. PN16-480102.