Different Methods to Calculate Dielectric Constants

Introduction

Dielectric materials are used to control and store charges and electric energies.They play a key role in modern electronics and electric power system.As the requirements grow for compact,low-cost electronic and electrical power system,as well as for very high energy and power capacitive storage systems,the development of high power and energy density dielectric materials becomes a major promising field. Dielectric constant is an important property for dielectric materials.In literature,there are several methods to calculate dielectric constant,the most popular one at present is density functional perturbation theory (DFPT).In this term paper,two methods will be presented,one is DFPT,the other one is beyond DFT,molecular dynamics.

Theory

In electrostatics,

(1)
\begin{align} D_i=E_i+4 \pi P_i=\epsilon_{ij}E_j \end{align}

This is the fundamental equation from which dielectric constant is calculated.The total polarization may usually be separated into three parts: electronic, ionic and dipolar, which results in the fact that dielectric constant also has three contributions.The electronic contribution arises from the electron cloud shifts from the nucleus,and the ionic contribution is caused by relative displacements between positive and negative ions. Permanent dipoles in polar systems result in the dipolar contribution to dielectric constant.

Density Functional Perturbation Theory (DFPT)
The Hohenberg and Kohn theorem states that all the physical properties of a system of interacting electrons are uniquely determined by its ground-state charge density distribution.This property holds independently of the precise form of the electron-electron interaction.This face was used by Kohn and Sham to map the problem of a system of interacting electrons onto an equivalent noninteracting problem.Since the electrons do not interact with each other, a one-electron Schrodinger equation can be written.

(2)
\begin{align} H_{SCF}\psi_n(\mathbf{r})=\epsilon_n\psi_n(\mathbf{r}) \end{align}
(3)
\begin{align} H_{SCF}=\frac{-\hbar^2}{2m}\frac{\partial^2 }{\partial \mathbf{r}^2}+V_{SCF}(\mathbf{r}) \end{align}
(4)
\begin{align} V_{SCF}(\mathbf{r})=V(\mathbf{r})+\nu_{xc}(\mathbf{r})+e^2\int\frac{n(\mathbf{r^'})d\mathbf{r^'}}{|\mathbf{r}-\mathbf{r^'}|} \end{align}
(5)
\begin{align} n(\mathbf{r})=2\sum_{n=1}^{N/2}{|\psi_n(\mathbf{r})|}^2 \end{align}

The above equations are the foundation of density funcstional theory (DFT),and once an explicit form for the exchange-correlation energy is available,equation (1) can be solved in a self-consistent way.
In quantum mechanics,perturbation theory is a set of approximation schemes directly related to mathematical perturbation for describing a complicated quantum system in terms of a simpler one.Equation (1) describes the unperturbed system.Now introduce a perturbation to the unperturbed Hamiltonian H0.Let V be a Hamiltonian representing a weak physical disturbance, such as a potential energy produced by an external field. Let λ be a dimensionless parameter that can take on values ranging continuously from 0 (no perturbation) to 1 (the full perturbation). The perturbed Hamiltonian is

(6)
\begin{align} H = H_0 + \lambda V \end{align}

Then the Schrodinger equation becomes:

(7)
\begin{align} \left(H_0 + \lambda V \right) |{\psi}\rangle = E_n |{\psi}\rangle \end{align}

If the perturbation is sufficiently weak, En and $|{\psi}\rangle$ can be written as power series in λ:

(8)
\begin{align} E_n = E_n^{(0)} + \lambda E_n^{(1)} + \lambda^2 E_n^{(2)} + \cdots \end{align}
(9)
\begin{align} |\psi\rangle = |\psi^{(0)}\rangle + \lambda |\psi^{(1)}\rangle + \lambda^2 |\psi^{(2)}\rangle + \cdots \end{align}

where

(10)
\begin{align} E_n^{(k)} = \frac{1}{k!} \frac{d^k E_n}{d \lambda^k} \end{align}

and

(11)
\begin{align} |\psi^{(k)}\rangle = \frac{1}{k!}\frac{d^k |\psi\rangle }{d \lambda^k} \end{align}

When ''λ = 0'', these reduce to the unperturbed values, which are the first term in each series.Substitute equation (7) and (8) into equation (6),thus

(12)
\begin{matrix} \left(H_0 + \lambda V \right) \left(|\psi^{(0)}\rangle + \lambda |\psi^{(1)}\rangle + \cdots \right) \qquad\qquad\qquad\qquad\\ \qquad\qquad\qquad= \left(E_n^{(0)} + \lambda E_n^{(1)} + \lambda^2 E_n^{(2)} + \cdots \right) \left(|\psi^{(0)}\rangle + \lambda |\psi^{(1)}\rangle + \cdots \right) \end{matrix}

The first order equation is

(13)
\begin{align} H_0 |\psi^{(1)}\rangle + V |\psi^{(0)}\rangle = E_n^{(0)} |\psi^{(1)}\rangle + E_n^{(1)} |\psi^{(0)}\rangle \end{align}

Apply the above first-order perturbation theory to the unperturbed Kohn-Sham equation,and set $|\psi^{(1)}\rangle=|\bigtriangleup\psi_n\rangle$,$V=\bigtriangleup V_{SCF}$,$E_n^{(1)}=\bigtriangleup\epsilon_n$

(14)
\begin{align} (H_{SCF}-\epsilon_n)|\bigtriangleup\psi_n\rangle=-(\bigtriangleup V_{SCF}-\bigtriangleup\epsilon_n)|\psi_n}\rangle \end{align}
(15)
\begin{align} \bigtriangleup n(\mathbf{r})=4\sum_{n=1}^{N/2}{\psi_n}^*(\mathbf{r})\bigtriangleup \psi_n(\mathbf{r}) \end{align}
(16)
\begin{align} |\bigtriangleup\psi_n\rangle=\sum_{m\not=n}|\psi_m\rangle\frac{\langle\psi_m|\bigtriangleup V_{SCF}|\psi_n\rangle}{\epsilon_n-\epsilon_m} \end{align}

In equation (13),$\bigtriangleup V_{SCF}$ is a linear functional of $\bigtriangleup n(\mathbf{r})$,then equation (13) and (14) form a set of self-consistent equations for the perturbed system completely analogous to the Kohn-Sham equations in the unperturbed case,and this becomes the foundation of the density functional perturbation theory.
From equation (15),the wave function response to a given perturbation depends only on the off-diagonal matrix elements of the perturbing potential between eigenfunctions of the unperturbed Hamiltonian.When the perturbation is a macroscopic electric field, then such matrix elements are indeed well defined,as can be seen by rewriting them in terms of the commutator between r and the unperturbed Hamiltonian,which is a lattice periodic operator:

(17)
\begin{align} \langle\psi_m|\mathbf{r}|\psi_n\rangle=\frac{\langle\psi_m|H_{SCF},\mathbf{r}|\psi_n\rangle}{\epsilon_m-\epsilon_n},\forall m\not=n \end{align}

When calculating the response of a crystal to an applied electric field E0,one must consider that the screened field acting on the electrons is

(18)
\begin{align} \mathbf{E}={\mathbf{E}}_0-4\pi\mathbf{P} \end{align}

Where P is the electronic polarization linearly induced by the screened field E:

(19)
\begin{align} \mathbf{P}=-\frac{e}{V}\int_V\mathbf{r}{\bigtriangleup^E}n(\mathbf{r})d\mathbf{r} \end{align}

Combine equation (14) and (18),

(20)
\begin{align} P_\alpha=-\frac{4e}{V}\sum_{n=1}^{N/2}\langle\psi_n|r_\alpha|{\bigtriangleup^E}\psi_n\rangle=-\frac{4e}{V}\sum_{n=1}^{N/2}\sum_{m=N/2+1}^{\infty}\frac{\langle\psi_n|[H_{SCF},r_\alpha]|\psi_m\rangle}{\epsilon_n-\epsilon_m}\times\langle\psi_m|{\bigtriangleup^E}\psi_n\rangle \end{align}

Where the subscript $\alpha$ indicates Cartesian components.Then introduce the wave function ${\overline{\psi}_n}^\alpha(\mathbf{r})$ defined as

(21)
\begin{align} {\overline{\psi}_n}^\alpha(\mathbf{r})=\sum_{m\not=n}\psi_m(\mathbf{r})\frac{\langle\psi_m|[H_{SCF},r_\alpha]|\psi_n\rangle}{\epsilon_m-\epsilon_n} \end{align}

Thus,

(22)
\begin{align} P_\alpha=-\frac{4e}{V}\sum_{n=1}^{N/2}\langle\overline{\psi}_n}^\alpha|\bigtriangleup^E\psi_n\rangle \end{align}

The electronic contribution to the dielectric tensor,${\epsilon_\infty}^{\alpha\beta}$,can be derived from simple electrostatics

(23)
\begin{align} E_{0\alpha}=(E_\alpha+4\pi P_\alpha)=\sum_{\beta}{\epsilon_\infty}^{\alpha\beta}E_\beta \end{align}

Using equation (21) to calculated the polarization induced in the $\alpha$ direction when a field is applied in the $\beta$ direction,thus

(24)
\begin{align} {\epsilon_\infty}^{\alpha\beta}=\delta_{\alpha\beta}-\frac{16\pi e}{VE_\beta}\sum_{n=1}^{N/2}\langle\overline{\psi}_n}^\alpha|\bigtriangleup^E\psi_n\rangle \end{align}

The total dielectric constant has two parts,one is the electronic contribution which has been derived,and the other one is ionic contribution.It can be defined as:

(25)
\begin{align} \epsilon={\epsilon_\infty}^{\alpha\beta}+\frac{4\pi e^2}{V}\sum_{\lambda}\frac{S_m}{{\omega_\lambda}^2} \end{align}

Where Sm is the mode oscillator strength, it is proportional to the born effective charge and the displacement of each atom.And $\omega_\lambda$ is the frequency of the λth IR-active phonon normal mode.In DFPT,the ionic contribution to dielectric constant can also be calculated,the detailed procedure is presented here,as it is similar to that of electronic contribution.

Molecular Dynamics
In the Frohlich's theory,the dielectric constant is related to the total dipole moment fluctuation.To derive this relation,there are usually two steps.First,with the help of statistical mechanics,establish a relation between the mean total dipole moment in the presence of a homogeneous external field and the equilibrium fluctuations in the absence of any external perturbations.Second,classical electrostatics is used to calculated the macroscopic polarization induced by electric field.After eliminating electric field from both equations one obtains a relation between the fluctuation and the dielectric constant.
For a system of N particles with generalized coordinates $q^N=(q_1,\ldots,q_N$ and interaction energy $U(q^N)$,enclosed in a volume V at temperature T, the mean value of the total dipole moment in the direction of a homogeneous external field E0,

(26)
\begin{align} {\langle M \rangle}_{E_0}=\frac{\int dq^N exp{(-[U(q^N)-M(q^N)E_0]/k_BT)} M(q^N)}{\int dq^N exp{(-[U(q^N)-M(q^N)E_0]/k_BT)}} \end{align}

to first order in E0 is given by

(27)
\begin{align} {\langle M \rangle}_{E_0}=E_0\frac{\langle \bigtriangleup M^2 \rangle}{3k_BT} \end{align}

In classical electrostatics,

(28)
\begin{align} \epsilon-1=4\pi{\langle P \rangle}}_{E_0}/E_0 \end{align}

The polarization is the dipole moment per unit volume,so,

(29)
\begin{align} \epsilon-1=\frac{4\pi\langle \bigtriangleup M^2 \rangle}{3Vk_BT} \end{align}
(30)
\begin{align} \langle \bigtriangleup M^2 \rangle=\langle M^2 \rangle-{\langle M \rangle}^2 \end{align}

According to this definition of dielectric constant in molecular dynamics,it may provide the ionic and dipolar contributions to dielectric constant.

Results and Discussion

To compare this two methods,the dielectric constants of crystal polyethylene and ice are calculated.
Polyethylene has –CH2-CH2- as its repeating unit. It has an orthorhombic unit cell in which there are two chains and each chain has one –CH2-CH2- unit.16 Figure 1 shows the crystal structure of Polyethylene.For ice,here ice XI is considered due to its simple structure.Ice XI is an orthorhombic, low-temperature equilibrium form of hexagonal ice. It is ferroelectric. It is considered the most stable configuration of ice Ih.Study indicated that the temperature below which ice XI forms is −36 °C (240K).The crystal structure of ice XI is shown in Figure 2.

The calculated dielectric constants are shown in the table below.The molecular dynamics calculation was done at 65 K.

The experiment value of polyethylene dielectric constant is around 2, so DFPT result shows good agreement with experiment value.From the structure of polyethylene,it does not have spontaneous polarization,and thus the dipolar contribution is 0,which is in agreement with the result of MD for polyethylene.For ice XI,it has permanent dipoles along the z direction.This can be verified by the dielectric constant calculated from MD.In x and y direction,there are no spontaneous polarization and the dielectric constant are due to ionic contribution.Compare the result of ionic contribution from DFPT with that of MD,they are very close to each other.While for z direction,MD has larger value than DFPT,which results from the dipolar contribution.

Summary

Based on the theory and results of DFPT and MD,a brief comparison of these two methods can be made:

• DFPT can give relatively accurate dielectric constants for both polar and nonpolar systems, while MD is suitable for polar system
• MD can calculate the dielectric constants of amorphous and fluid systems, but for DFPT it is difficult.
• DFPT gives the dielectric constant at 0 K, while MD can calculate dielectric constant at different temperatures.
• They may be compensation to each other. For polar system: $\epsilon_{tol}$= Electronic (DFPT) + Ionic (DFPT) + Dipolar (MD)

References
(1) C. Kittel, Introduction to Solid State Physics
(2) A. Messiah, Quantum Mechanics
(3) S. Baroni, Stefano de Gironcoli, and Andrea Dal Corso. Rev. Mod. Phys. 73 (2001)
(4) X. Gonze and C. Lee. Phys.Rev.B 55 10355-10368 (1996)
(5) X. Gonze. Phys.Rev.B 55 10337-10354 (1996)
(6) M. Neumann, Molecular Physics 50 841-858 (1983)
(7) D. J. Price and C. L. Brooks III. J.Chem.Phys 121 (2004)
(8) R. Howe and R. W. Whitworth. J.Chem.Phys 90 (1989)