# Modern theory of Polarization¶

Contents

This document contains the detailed explanation of the theory of the Modern theory of polarization.

## Summary¶

The standard definition of the Polarization density is not well-defined in periodic system, however, in the Modern theory of Polarization it is demonstrated that another definition exists that relates to changes in the Polarization density. This allows for a straight-forward way to calculate derivatives of the Polarization density, and is therefore useful for calculating the piezoelectric tensor which describes the change in the polarization density in response to an applied strain.

## Relevant recipes¶

The following recipes are all relying on the Modern theory of Polarization:

`asr.borncharges.main()`

## Alternative literature¶

To learn more about this please consider

Good resource for beginners: N.A. Spaldin, “A beginners guide to the modern theory of polarization”, Journal of Solid State Chemistry 195 (2012) 2–10.

Good review paper: R. Resta, “Macroscopic polarization in crystalline dielectrics: the geometric phase approach”, Rev. Mod. Phys. 66 (1994) 899–915.

Paper with historical significance: R.D. King-Smith, D. Vanderbilt, “Theory of polarization of crystalline solids”, Phys. Rev. B 47 (1993) R1651–R1654.

Paper with historical significance: R. Resta, “Macroscopic Electric Polarization as a Geometric Quantum Phase”, Eur. Phys. Lett. 22 (1993) 133–138.

## Explanation¶

The polarisation of a dielectric medium represents the dipole per unit volume. As such it would be natural to define a macroscopic (i.e. unit cell averaged) polarisation for an infinite periodic bulk system as

where \(n\) includes the ion point charges as well as the
delocalised electron charge (\(e > 0\)). It is, however, obvious
that this definition is not meaningful as it depends on the chosen unit
cell (see p. 365 in Grosso and Parravicini 2nd. ed.). Instead, it turns
out that only *changes* in polarisation are physically meaningful, and
in fact all experimental measurements of bulk polarisation indeed probe
the difference in polarisation between two states of the crystal. In
this and the next section we show, following two different routes, that
such quantities are indeed well defined.

Before we derive a unit-cell independent formula for the change in polarization, it should be noted that expression (1) cannot be directly applied to determine such a change. This point should be clear from Figure Fig. 1.

Rather than starting from the unit cell dependent formula (1), we consider the polarization of a finite piece of the bulk for which (1) is meaningful when \(V_{\textrm{cell}}\) is replaced by the total volume of the crystal, \(V\). The idea is now to calculate the change in \(\mathbf P\) induced by some change in the Hamiltonian, and then show that taking the thermodynamic limit (\(V\to \infty\)) of the polarisation change is mathematically well defined.

In the following we consider the change in polarisation when the potential is changed adiabatically from \(v_{\lambda=0}\) to \(v_{\lambda=1}\). We have

and from Eq. (1) we can write

Using first order perturbation theory we have

Using the commutator relation \([\mathbf r,H_{\lambda}]=i\hbar\mathbf p / m_e\), the off-diagonal matrix elements of the position operator can be rewritten

and we finally arrive at the expression

This quantity is well defined for any piece of material also for a periodic solid in the thermodynamic limit. It does not depend on the choice of unit cell (because it makes no reference to the unit cell) and it is independent of the phases chosen for the Bloch states.

### Polarization change from Kubo formula¶

In the previous section, the problem with the unit cell dependent expression (1), was circumvented by considering a finite piece of material and then taking the thermodynamic limit. In this section we present an alternative formulation which defines the polarization from the current flowing through a unit cell in response to a periodic adiabatic change in the potential.

Thus we consider the current flow produced by the adiabatic change in the potential from \(v_{\lambda=0}\) to \(v_{\lambda=1}\), where \(v_\lambda\) is assumed to be periodic for all \(\lambda\). The (microscopic) polarizability is related to the current density via

As a quantum mechanical operator we have \(\frac{\partial \mathbf P(\mathbf r) }{ \partial t} = [\mathbf P, H] / i\hbar\). Thus when considering off-diagonal matrix elements of \(\mathbf P\) on energy eigenstates we have

Since we are interested in the macroscopic polarisation we perform a unit cell average. Thanks to the Bloch form of the wave functions, \(\psi_{nk}(\mathbf r)=e^{i\mathbf{k}\cdot \mathbf{r}}u_{nk}(\mathbf r)\), we have

Suppose the system is in the ground state of \(H(\lambda)\). We now consider the change in \(\mathbf P\) when the Hamiltonian is changed adiabatically to \(H(\lambda + d\lambda)\). This change can be obtained from the Kubo formula using \(\mathbf P\) as the observable and \(dH(\lambda)=\frac{\partial v_\lambda}{ \partial \lambda} d\lambda\) as the time-independent perturbation. The finite imaginary frequency \(i\eta\) in the Kubo formula ensures that the perturbation is switched on adiabatically so that the system stays in the ground state. With this we obtain

which coincide with Eq. (3).

### Polarisation change as a Berry phase on the occupied manifold¶

Eq. (3) uniquely specifies the macroscopic polarisation change due to an adiabatic change of the crystal potential. It has the drawback that it involves a sum over unoccupied states making it costly to evaluate in practice. As shown below, it is possible to obtain an expression involving only the occupied subspace. Furthermore, it is shown that the polarization change, \(\Delta \mathbf P\), can be calculated from knowing only its value at the end points of the adiabatic path \(\lambda=0..1\). This comes, however, at the price of an introduced ambiguity, namely that the polarisation change can be determined only up to an integer number of polarisation quanta, \(e L / V_\mathrm{cell}\), where \(L\) is the unit cell length. In practice, however, this is not a problem because \(|\Delta \mathbf P|\ll e L / V_\mathrm{cell}\).

We use the relations

where the cell periodic Hamiltonian is given by

It should be noted that for the above relations to hold it is essential that the cell-periodic functions, \(u^\lambda_{nk}\), are analytic with respect to \(\mathbf k\) and \(\lambda\). Substituting into Eq. (3) we obtain (after some manipulations)

XXX (show this!). It can be shown (see e.g. Grosso and Paravicini) that the above expression can be rephrased as

where

Considering the polarisation along a particular direction, say the \(z\)-axis, the derivative only connects Bloch states along \(\mathbf k_z\). In this case the BZ integral can be discretised in the directions perpendicular to \(z\), and the contribution for each \(\mathbf k_{\perp}\) becomes

where \(A\) is the area of the unit cell in \(xy\) plane. We can write this as

where

is nothing but the Berry phase picked up along the 1D BZ. As always the expression is invariant under a change in the phases of the wave functions, \(e^{i\theta(\mathbf k)}\), as long as \(\theta\) is differentiable on the BZ torus (i.e. with periodic boundary conditions). We notice, however, that in contrast to the normal Berry phase, the Hamiltonian \(H(\mathbf k,\lambda)\), from which the cell-periodic functions derive, is not cyclic over the 1D BZ because \(H(\mathbf k,\lambda)=H(\mathbf k+\mathbf G,\lambda)\) only modulus a gauge transformation, i.e. a unitary transformation of the form \(\exp(i\chi(\mathbf r))\). This means that

(which is not just a phase factor). We refer to this relation as the periodic gauge.

Now, we show that Eqs. ((4) - (6)) only determine \(\Delta P\) up to an integer number of polarisation quanta. To this end consider the special case where the Hamiltonians at \(\lambda=0\) and 1 are identical, e.g. an atom is moved along a closed loop. In this case \(u_{n\mathbf k}^{(0)}\) and \(u_{n\mathbf k}^{(1)}\) can at most differ by a phase,

Inserting this in Eq. (6) yields

Because of Eq. (7) we must have \(e^{i\theta_{n\mathbf k}}=e^{i\theta_{n,\mathbf k+\mathbf G}}\) meaning that

where \(\beta\) is BZ-periodic (and differentiable) in \(\mathbf k\). We thus conclude that for \(H(\lambda=0)=H(\lambda=1)\) we have

where \(V_\mathrm{cell} = Ac\). This shows that the polarisation change in direction \(\alpha\) is only determined up to the polarisation quantum \((e/V_{\textrm{cell}})L_{\alpha}\).

Eqs. ((4) - (5)) invites the interpretation in terms of an absolute polarisation. However, as previously discussed such a concept is not well defined. Thus \(\mathbf P(\lambda)\) only makes sense as a device to compute the change in polarisation (which when evaluated in terms of the Berry phase is defined only modulus the polarisation quantum).

### Practical calculations of the polarization¶

We now describe how the Berry phase theory can be used to calculate real world quantities in practice. Eq. (5) is slightly rewritten to make apparent its use of a trace

where it is understood that the inside of the trace is a matrix in band-indices \(n,m\) and that trace is taken over the occupied manifold of bands. \(\mathbf{n}\) is a unit-vector along the direction the polarization is calculated. The derivative of the Bloch-functions is expanded to first order in \(\mathbf{k}\)

leading to the approximate expression for the polarization

(Here we have removed \(\mathbf{n}\cdot\) as it should be clear that the polarisation along a specific direction is obtained by dotting with \(\mathbf{n}\)). In principle, this expression can be straightforwardly implemented numerically. However, it appears that the result depend on the (arbitrary) phases of the Bloch states. Eq. (8) requires that the \(u_{n\mathbf k}\) are differentiable with respect to \(\mathbf k\). But the wave functions obtained from practical DFT codes come with arbitrary phases. To show that the result is in fact independent of the phases, we use that the logarithm of a matrix, \(S\), which is close to the unit matrix, to first order is

which allows us to write

Now we can use the fact that the trace of a logarithm of a matrix is equal to the logarithm of the determinant

(which can be confirmed by inserting the eigen-representation of \(S\)) yielding

Finally we can pull the sum into the logarithm by converting it to a product

This expression shows that the polarization is in fact independent of the arbitrary phases of the wave functions. It is implemented in the GPAW code, and one example of its use will be illustrated in the next section.

### The Born charge¶

Now we consider the induced polarization when displacing an atom in a crystal from its equilibrium position. If the atom is ionized and thus have donated or accepted a finite number of electrons (like in the NaCl crystal), the induced polarisation can be simply given by the charge of the ion multiplied by the displacement \(\delta \mathbf{R}\)

Here \(Z_\mathrm{ion}\) is a number describing the net-charge associated with the ion. If the electrons are strongly bound to the ion they will follow the displacement of the ion and \(Z\) will be expected to be an integer, however, in the general case where electrons do not strictly follow the displacement of the ion, \(Z\) will be a fractional number known as the Born charge. The Born charge of a given atom, \(a\), in a crystal is a tensor defined as

where \(i,j=x,y,x\) denote the direction. In this equation it is understood that atom \(a\) in all unit cells are displaced such that the assumption of a periodic perturbation behind Eq. (9) is satisfied. At this point it is instructive to recall the definition of the electronic dielectric tensor and susceptibilites that we have studied so far in the course:

and

In writing these relations we have suppressed the \(q\)- and \(\omega\)-dependence of the response functions. The important point is the “el” superscript, which indicates that the induced polarization is created by the electrons moving in the frozen crystal, i.e without allowing the atoms to move. To obtain total dielectric tensor and susceptibilities we must add the ionic part describing the additional polarization due to the vibrating lattice. The calculation of the ionic contribution to the dielectric function requires the vibrational frequencies of the lattice, i.e. the phonons, and the Born charges, as input. If you would like to see how this goes, consult page 417-419 and 423-424 in GP.

In practice, formula (10) is evaluated as a finite difference

Finally, we need to use Eq. (9) to calculate the polarisation at a finite displacement of the atoms. However, it is important to remember that the complex logarithm has a branch cut which typically lies from \([-\infty, 0]\), which can lead to discontinuous jumps of the integrand in Eq. (9) yielding unphysical results (the integrand should be continuous). An example is shown in Fig. Fig. 2 for the two-dimensional material MoS\(_2\) where the integral (over \(\mathbf{k}_\perp\)) is one-dimensional and therefore can be easily plotted. Here it is clear that the branch cut of the logarithm is being crossed leading to discontinuous jumps in the integrand (blue line). This can be fixed by comparing neighbouring k-points in the integrand and adding or subtracting a multiple of \(2\pi\) to ensure that the Berry phases change slowly as a function of \(\mathrm{k}_\perp\) (orange lines). Using this scheme we find that two-dimensional MoS\(_2\) in the H-phase has the following Born charges: \(Z^\mathrm{Mo}_{[xx, yy, zz]} = [-1.07, -1.07, -0.13]\) and \(Z^\mathrm{S}_{[xx, yy, zz]} = [0.53, 0.53, 0.07]\) (all off-diagonal elements are zero). Now it can be seen that \(Z^\mathrm{S} \approx -Z^\mathrm{Mo} / 2\) which is actually a variant of a deeper principle known as the acoustic sum rule which says that \(\sum^A Z^A_{ij} = 0\) (when the net-charge of the total cell is zero), where \(A\) runs over all atoms in the unit cell. It is interesting to note that the Born charges of S are positive and those of Mo are negative while the opposite is found for the net charge of the atoms in the equilibrium structure (S takes electron density from Mo). This shows that the concept of Born charges on covalently bonded structures like MoS\(_2\) is highly non-trivial.