# Berry-phase theory of proper piezoelectric response

## I Introduction

The calculation of spontaneous polarization and piezoelectric response within the framework of first-principles methods of electronic structure theory has proven to be a rather subtle problem. In a landmark paper, Martin  showed that the piezoelectric tensor is well-defined as a bulk quantity in a crystalline insulator. However, at that time it was far from clear whether the spontaneous polarization itself could be regarded as a bulk property in the same sense, and calculations of piezoelectric constants by finite differences of spontaneous polarization were therefore not possible.

The situation changed in 1993 with the development of the “Berry-phase” theory of polarization , which provided a direct and straightforward method for computing the electric polarization. (For a useful review, see Ref. .) Nevertheless, some subtleties remain regarding the computation of the piezoelectric tensor components by finite differences . First, the Berry-phase theory gives the polarization as a multivalued quantity, and the piezoelectric response that would be computed from a given one of the many branches is not invariant with respect to choice of branch. Second, a distinction is made between the “proper” and “improper” piezoelectric response , and it might not be clear which of these is to be associated with the finite-difference calculation.

The purpose of the present paper is to elucidate the physics of the spontaneous polarization, the piezoelectric response, and the relations between the two. It is clarified that the improper polarization is the one given by the naive finite-difference approach, and that while this quantity is indeed branch-dependent, the proper polarization, which should be compared with experiment, is not. As a result of this analysis, a simplified recipe for the direct finite-difference computation of the proper piezoelectric response is given.

## II Berry-phase theory of polarization

We consider a periodic insulating crystal in zero macroscopic electric field, and assume that the electronic ground state can be described by a one-electron Hamiltonian $`H`$ as in density-functional or Hartree-Fock theory. The eigenstates of $`H`$ are the Bloch functions $`\psi _{n𝐤}`$ with energies $`ϵ_{n𝐤}`$, and it is conventional to define the cell-periodic Bloch functions

$$u_{n𝐤}(𝐫)=e^{i𝐤𝐫}\psi _{n𝐤}(𝐫)$$
(1)   
having periodicity $`u_{n𝐤}(𝐫)=u_{n𝐤}(𝐑+𝐫)`$, where $`𝐑`$ is any lattice vector. The contribution of the $`n`$’th occupied band to the spontaneous electric polarization of the crystal can then be written 

$$𝐏_n=\frac{ie}{(2\pi )^3}d^3ku_{n𝐤}|_𝐤|u_{n𝐤}.$$
(2)   
We take the convention that $`n`$ runs over bands and spin, so a factor of two would need to be inserted in Eq. (2) to account for paired spins. The total spontaneous polarization is then given by

$$𝐏=\frac{e}{\mathrm{\Omega }}\underset{\tau }{}Z_\tau 𝐫_\tau +\underset{n\mathrm{occ}}{}𝐏_n,$$
(3)   
where $`Z_\tau `$ and $`𝐫_\tau `$ are the atomic number and cell position of the $`\tau `$’th nucleus in the unit cell, and $`\mathrm{\Omega }`$ is the unit cell volume.

Strictly speaking, Eq. (2) applies only to an isolated band, i.e., a band for which $`ϵ_{n𝐤}`$ does not become degenerate with any other band at any point in the Brillouin zone. This restriction is not essential; methods for extending the analysis to composite groups of occupied bands containing arbitrary degeneracies and crossings have been developed as described in Refs. . However, for simplicity of presentation, it will be assumed here that only isolated bands are present. For the same reason, spin degeneracy is suppressed throughout.

There is a certain arbitrariness inherent in Eq. (2) associated with the freedom to choose the phases of the Bloch functions $`\psi _{n𝐤}`$. For, suppose we make a different choice

$$|\stackrel{~}{\psi }_{n𝐤}=e^{i\beta (𝐤)}|\psi _{n𝐤}.$$
(4)   
We shall refer to this as a “gauge transformation” of the Bloch functions. Note that the choice of $`\beta (𝐤)`$ is restricted by the fact that $`𝐤`$ and $`𝐤+𝐆`$ label the same wavefunction (where $`𝐆`$ is a reciprocal lattice vector), so that $`\beta (𝐤+𝐆)\beta (𝐤)`$ must be an integer multiple of $`2\pi `$ for any $`𝐆`$. Thus, the most general form of $`\beta (𝐤)`$ is

$$\beta (𝐤)=\beta _{\mathrm{per}}(𝐤)+𝐤𝐑$$
(5)   
where $`\beta _{\mathrm{per}}`$ is a periodic function in $`k`$-space and $`𝐑`$ is some real-space lattice vector. Letting $`\stackrel{~}{𝐏}_n`$ be the result of inserting the $`\stackrel{~}{u}_{n𝐤}`$ in place of the $`u_{n𝐤}`$ in Eq. (2), and using Eqs. (1), (4), and (5), one finds

$$\stackrel{~}{𝐏}_n=𝐏_n\frac{e𝐑}{\mathrm{\Omega }}.$$
(6)   
Thus, while the contribution of this band to the electronic polarization is not absolutely gauge-invariant, it is gauge-invariant modulo $`e/\mathrm{\Omega }`$ times a real-space lattice vector. Actually, this is precisely the type of qualified invariance we should have expected. After all, the choice of the location $`𝐫_\tau `$ of the atom representing sublattice $`\tau `$ in the unit cell has a similar ambiguity; we could just as well choose $`\stackrel{~}{𝐫}_\tau =𝐫_\tau +𝐑^{}`$, where $`𝐑^{}`$ is another lattice vector, leading to precisely the same kind of “modulo $`e𝐑/\mathrm{\Omega }`$” ambiguity in the expression for $`𝐏`$ in Eq. (3).

Perhaps the most natural way to incorporate this kind of ambiguity in the definition of the polarization is to regard $`𝐏`$ as a multivalued quantity; that is, it simultaneously takes on a lattice of values given by some $`𝐏^{(b)}`$ (here ‘$`b`$’ is a “branch” label) and all its periodic images $`𝐏^{(b)}+e𝐑/\mathrm{\Omega }`$ (with $`R`$ running over all lattice vectors of the crystal). To interpret this intuitively, we can say that from the point of view of its dipolar properties, the real insulator behaves like a fictitious crystal composed of two sublattices of point $`\pm e`$ charges, with the sublattice of $`e`$ charges displaced relative to the $`+e`$ sublattice by $`\mathrm{\Omega }𝐏/e`$. That is, choosing one of the $`+e`$ charges as the origin, $`\mathrm{\Omega }𝐏/e`$ takes on a lattice of values that is precisely the lattice of positions of the $`e`$ charges.

This is illustrated in Fig. 1 for an imaginary tetragonal crystal (dimensions $`a\times a\times c`$) with one monovalent ion located at the cell corners, and a single (spinless) electron band giving rise to the distributed electron charge indicated schematically by the contours in Fig. 1(a). (We assume that $`M_z`$ mirror symmetry is broken in some way.) Eq. (2) then gives the location

$$𝐫_n=\mathrm{\Omega }𝐏_n/e$$
(7)   
of the effective unit point charge $`e`$ illustrated in Fig. 1(b). As discussed in Refs. , this location is just the charge center of the Wannier function associated with the electron band. The polarization will then take on a lattice of values having $`x`$, $`y`$, and $`z`$ components of $`m_1e/ac`$, $`m_2e/ac`$, and $`(\gamma +m_3)e/a^2`$, respectively, where the $`m_i`$ are integers. More generally, when several occupied bands are present, one can rewrite Eq. (3) as

$$𝐏=\frac{e}{\mathrm{\Omega }}\underset{\tau }{}Z_\tau 𝐫_\tau \frac{e}{\mathrm{\Omega }}\underset{n\mathrm{occ}}{}𝐫_n.$$
(8)   
In practice, one proceeds by computing the component of $`𝐏_𝐧`$ along a particular crystallographic direction $`\alpha `$ via the quantity

$$\varphi _{n,\alpha }=\frac{\mathrm{\Omega }}{e}𝐆_\alpha 𝐏_n,$$
(9)   
where $`𝐆_\alpha `$ is the primitive reciprocal lattice vector in direction $`\alpha `$. In cases of simple symmetry (e.g., tetragonal or rhombohedral ferroelectric phases), a single $`\varphi _n`$ suffices to determine $`𝐏_n`$, but in general $`𝐏_n`$ can be reconstructed from the three $`\varphi `$’s via

$$𝐏_n=\frac{1}{2\pi }\frac{e}{\mathrm{\Omega }}\underset{\alpha }{}\varphi _{n,\alpha }𝐑_\alpha ,$$
(10)   
where $`𝐑_\alpha `$ is the real-space primitive lattice vector corresponding to $`𝐆_\alpha `$. The $`\varphi _{n,\alpha }`$ are angle variables (“Berry phases”) that are well-defined modulo $`2\pi `$, given by

$$\varphi _{n,\alpha }=\mathrm{\Omega }_{\mathrm{BZ}}^1_{\mathrm{BZ}}d^3ku_{n𝐤}|i𝐆_\alpha _𝐤|u_{n𝐤},$$
(11)   
where $`\mathrm{\Omega }_{\mathrm{BZ}}=(2\pi )^3/\mathrm{\Omega }`$ is the Brillouin zone (BZ) volume.

The $`\varphi _{n,\alpha }`$ can be regarded as giving the position of the Wannier center for band $`n`$. For the toy crystal of Fig. 1, for example, and with the origin chosen at the cell corner, one would have $`\varphi _x=\varphi _y=0`$ and $`\varphi _z=2\pi \gamma `$. The practical calculation of the $`\varphi _{n,\alpha }`$ proceeds on a discrete mesh in reciprocal space, arranged as a two-dimensional grid of $`𝐆_\alpha `$-oriented strings of k-points, as described in Refs. .

## III Piezoelectric response

The piezoelectric tensor of a crystal reflects the first-order change in spontaneous electric polarization in response to a first-order deformation of the crystal. The “improper” piezoelectric tensor is defined as 

$$c_{ijk}=\frac{dP_i}{dϵ_{jk}}$$
(12)   
in terms of the deformation

$$dr_j=\underset{k}{}dϵ_{jk}r_k,$$
(13)   
where the symmetric and antisymmetric parts of $`dϵ`$ represent infinitesimal strains and rotations, respectively. On the other hand, the “proper” piezoelectric tensor can be defined as

$$\stackrel{~}{c}_{ijk}=\frac{dJ_i}{d\dot{ϵ}_{jk}},$$
(14)   
where $`𝐉`$ is the current density that flows through the bulk of the sample in adiabatic response to a slow deformation $`\dot{ϵ}=dϵ/dt`$. According to the standard references , the relation between the improper and proper piezoelectric tensors is

$$\stackrel{~}{c}_{ijk}=c_{ijk}+\delta _{jk}P_i\delta _{ij}P_k.$$
(15)   
Writing out explicit tensor components, this last equation becomes 

$`\stackrel{~}{c}_{zzz}`$ $`=`$ $`c_{zzz},`$ (16)   
$`\stackrel{~}{c}_{zxx}`$ $`=`$ $`c_{zxx}+P_z,`$ (17)   
$`\stackrel{~}{c}_{zxy}`$ $`=`$ $`c_{zxy},`$ (18)   
$`\stackrel{~}{c}_{zxz}`$ $`=`$ $`c_{zxz},`$ (19)   
$`\stackrel{~}{c}_{zzx}`$ $`=`$ $`c_{zzx}P_x,`$ (20)   
and similarly for permutations of the cartesian labels (but not for permutations of their position in the index triplet). It might seem strange at first sight that the expressions for $`\stackrel{~}{c}_{zxz}`$ and $`\stackrel{~}{c}_{zzx}`$ have a different form, but this just reflects the fact that the deformation tensor $`ϵ`$ has been allowed to contain an antisymmetric part.

Now in the Berry-phase theory, the polarization is a multivalued quantity, so that any particular value $`𝐏^{(b)}`$ has to be identified by its branch label ‘$`b`$’, and the corresponding improper piezoelectric tensor is

$$c_{ijk}^{(b)}=\frac{dP_i^{(b)}}{dϵ_{jk}}.$$
(21)   
Since $`𝐏`$ is well-defined modulo $`e𝐑/\mathrm{\Omega }`$, and both $`𝐑`$ and $`\mathrm{\Omega }`$ vary with the deformation $`ϵ`$, Eq. (21) will clearly give different results for different choices of branch. This branch-dependence is problematic; the piezoelectric tensor is measurable, and a suitable theory ought to give a unique value for it.

Before proceeding, the reader is reminded that the piezoelectric response contains, in general, a “clamped-ion” part and an “internal-strain” part . That is, one decomposes the actual deformation into a sum of two parts: a homogeneous strain in which the nuclear coordinates follow Eq. (13) exactly (clamped-ion part), plus an internal distortion of the nuclear coordinates at fixed strain (internal-strain part). Since the latter occurs at fixed strain, all the subtleties about the branch-dependence and the proper-vs.-improper distinction disappear for this case. While the computation of the internal-strain part of the piezoelectric response may be tedious (requiring an iterative set of force calculations to determine the needed internal relaxations), it is straightforward in principle. Consequently, for the remainder of this paper, the discussion refers to the clamped-ion response only unless explicitly stated otherwise.

### A Branch-invariance of proper piezoelectric response

While it is true that the improper piezoelectric response depends, in general, on choice of branch, it is instead the proper piezoelectric tensor that should be compared with experiment. Figure 2 shows a sketch of one possible experimental setup, in which a block of piezoelectric material is sandwiched between shorted conducting electrodes, and the current $`I`$ that flows in response to a deformation $`ϵ`$ is measured. As suggested by Eq. (14), the proper piezoelectric response is related to the current that flows through the sample in response to the deformation, and is thus the experimentally measured quantity. Moreover, the induced current density $`𝐣(𝐫)`$ is periodic with the lattice, so that its unit cell average $`𝐉`$ in Eq. (14) is perfectly well-defined, and consequently the proper piezoelectric tensor $`\stackrel{~}{c}`$ cannot suffer from any dependence upon choice of branch.

It is straightforward to check this branch-independence of $`\stackrel{~}{c}`$ explicitly. Since the polarizations for two different branch choices are related by

$$𝐏^{(b^{})}=𝐏^{(b)}+\frac{e}{\mathrm{\Omega }}𝐑,$$
(22)   
one finds

$`dP_i^{(b^{})}`$ $`=`$ $`dP_i^{(b)}{\displaystyle \frac{e}{\mathrm{\Omega }^2}}d\mathrm{\Omega }R_i+{\displaystyle \frac{e}{\mathrm{\Omega }}}dR_i`$ (23)   
$`=`$ $`dP_i^{(b)}+{\displaystyle \frac{e}{\mathrm{\Omega }}}{\displaystyle \underset{l}{}}(dϵ_{ll}R_i+dϵ_{il}R_l).`$ (24)   
so that

$$c_{ijk}^{(b^{})}=c_{ijk}^{(b)}\frac{e}{\mathrm{\Omega }}\delta _{jk}R_i+\frac{e}{\mathrm{\Omega }}\delta _{ij}R_k,$$
(25)   
or, using Eq. (22),

$$c_{ijk}^{(b^{})}+\delta _{jk}P_i^{(b^{})}\delta _{ij}P_k^{(b^{})}=c_{ijk}^{(b)}+\delta _{jk}P_i^{(b)}+\delta _{ij}P_k^{(b)}.$$
(26)   
It is thus evident that $`\stackrel{~}{c}_{ijk}`$ as defined in Eq. (15) is indeed independent of choice of branch.

It is instructive to note that a similar argument applies to the part of the proper piezoelectric tensor arising from the ionic contribution $`𝐏_{\mathrm{ion}}=(e/\mathrm{\Omega })_\tau Z_\tau 𝐫_\tau `$ in Eq. (3). Recalling that we are working in the clamped-ion approximation, so that $`d𝐫_\tau `$ follows the form of Eq. (13), one finds immediately that $`\stackrel{~}{c}_{\mathrm{ion}}=0`$ by the same logic as for the previous paragraph.

Indeed, the same logic would apply to Eq. (8) if the Wannier centers $`𝐫_n`$ would undergo a homogeneous deformation of the type (13). In other words, the proper piezoelectric response is identically zero for a homogeneous deformation of both the ionic positions and the Wannier centers, in which case there is no charge flow through the interior of the crystal.

### B Simplified finite-difference formula

Of course, there is no reason to expect the Wannier centers $`𝐫_n`$ to follow a homogeneous deformation, so $`\stackrel{~}{c}`$ is not generally zero. But from this point of view, it becomes evident that the the proper piezoelectric response is precisely a measure of the degree to which the Wannier centers fail to follow a homogeneous deformation. Or equivalently, returning to Eq. (10), we see that the proper piezoelectric response measures just the variation of the Berry phases $`\varphi _{n,\alpha }`$ with the strain deformation. More precisely, starting from Eqs.(10), (12), and (15), one finds

$$\stackrel{~}{c}_{ijk}=\frac{1}{2\pi }\frac{e}{\mathrm{\Omega }}\underset{n,\alpha }{}\frac{d\varphi _{n,\alpha }}{dϵ_{jk}}R_{\alpha i}$$
(27)   
We have been working in the clamped-ion approximation, but in general, if there are internal relaxations accompanying the deformation, one can define a total Berry phase in direction $`\alpha `$,

$$\varphi _\alpha =\underset{\tau }{}Z_\tau 𝐆_\alpha 𝐫_\tau \underset{n}{}\varphi _{n,\alpha },$$
(28)   
so that

$$\stackrel{~}{c}_{ijk}=\frac{1}{2\pi }\frac{e}{\mathrm{\Omega }}\underset{\alpha }{}\frac{d\varphi _\alpha }{dϵ_{jk}}R_{\alpha i}.$$
(29)   
Naturally, the ionic contributions to $`d\varphi _\alpha /dϵ_{jk}`$ vanish in the clamped-ion approximation.

Equation (27), or its generalization (29), is the central result of this paper, and provides a simple and practical recipe for calculating the desired proper piezoelectric response. One simply computes the needed $`d\varphi /dϵ`$ by finite differences, as $`(\varphi ^{}\varphi )/(ϵ^{}ϵ)`$ for nearby strain configurations $`ϵ`$ and $`ϵ^{}`$. Then these $`d\varphi /dϵ`$ are inserted into Eq. (27) or (29) to obtain the elements of the proper piezoelectric tensor.

### C Relation to surface charges

At the end of Sec. III A, it was pointed out that a homogeneous deformation of the lattice of positive ionic and negative Wannier-center point charges would give rise to no internal current, and hence no proper piezoelectric response. This result can be made more intuitive by considering the connection between bulk polarization and surface charges .

Consider, for example, a crystallite composed of $`N\times N\times N`$ replicas of the unit cell shown in Fig. 1(b). In general there may be an arbitrariness in the choice of surface termination, as illustrated for the top surface of this crystallite in Fig. 3. For any given termination, the macroscopic surface charge density $`\sigma `$ is uniquely defined as $`𝑑z\overline{\rho }(z)`$, where $`\overline{\rho }(z)`$ is the average charge contained in a unit cell centered at vertical coordinate $`z`$ (so that $`\overline{\rho }`$ vanishes either deep in the crystal or deep in the vacuum and its integral is convergent). For the crystal of Fig. 1, one finds $`\sigma =\gamma e/a^2`$ and $`\sigma =(\gamma 1)e/a^2`$ for the terminations of type I and II of Figs. 3(a) and 3(b), respectively. Referring back to Sec. II, where it was found that $`P_z=(\gamma +m_3)e/a^2`$, one confirms that the relation 

$$\sigma =𝐏\widehat{𝐧}$$
(30)   
is satisfied for both terminations, the ambiguity of termination corresponding to the choice of branch of $`𝐏`$.

For definiteness, assume that the surface terminations are such that the top and bottom surfaces of the crystallite have charge densities $`+\gamma e/a^2`$ and $`\gamma e/a^2`$ on the top and bottom surfaces, respectively, and zero on the sides. Then the magnitude of the total charge on the top or bottom surface is just $`N^2a^2\sigma =N^2\gamma e`$, which is clearly independent of any homogeneous ($`\gamma `$-preserving) deformation of the crystal. Thus, if this crystallite were inserted between grounded capacitor plates as in Fig. 2, no current would flow through the wire as a result of the homogeneous deformation. This is consistent with the vanishing of the proper piezoelectric response associated with such a homogeneous deformation, as already illustrated via Eq. (27).

However, for the same situation, the improper piezoelectric tensor would have nonzero elements. For the chosen surface termination, the crystallite has a total dipole moment $`𝐝=N^3\gamma ec\widehat{𝐳}`$, and a polarization $`𝐏=𝐝/N^3a^2c=(\gamma e/a^2)\widehat{𝐳}`$ as expected. Clearly this $`𝐏`$ is invariant with respect to an elongation of the crystallite along the $`\widehat{𝐳}`$ axis (strain component $`ϵ_{zz}`$), but not to an elongation along the $`\widehat{𝐱}`$ or $`\widehat{𝐲}`$ axes ($`ϵ_{yy}`$ or $`ϵ_{zz}`$), thus explaining why there is a correction to $`c_{zxx}`$ but not to $`c_{zzz}`$ in Eq. (20). Similar considerations applied to shear strains and rotations explain the remaining entries in Eq. (20).

## IV Discussion

As is evident from Eqs.(15) and (20), the distinction between the proper and improper piezoelectric tensor is only present if a spontaneous polarization is present. If the spontaneous polarization is small, as for wurtzite semiconductors , it may be a good approximation to neglect the corrections to the improper tensor components. Alternatively, linear-response methods can be used to compute the proper piezoelectric response directly . However, for a finite-difference calculation of the proper piezoelectric response of a ferroelectric material, it is essential to take the corrections to the improper response explicitly into account, as was done in Ref. .

## V Summary

In this work, a simple and straightforward method for computing the proper piezoelectric response has been proposed. Instead of first computing the improper response and then the needed corrections, the proper response is computed directly from Eq. (27) or (29). It is thus clarified that the central quantities needed to determine the proper piezoelectric response are just the variations of the Berry phases with deformation.

###### Acknowledgements.

This work was supported by ONR Grant N0001497-1-0048. I wish to thank R. Cohen for calling my attention to the problem of the proper piezoelectric response and its relation to the Berry-phase theory of polarization.