Category:Dielectric properties

Introduction

When an external electric field ${\displaystyle \mathbf {E} }$ acts on a medium, both the electronic and ionic charges will react to the perturbing field . For dielectric materials, in a very simplistic approach, one can think that the bound charges will create dipoles inside the medium, leading to an induced polarization, ${\displaystyle \mathbf {P} }$. The combined effects of both fields are expressed in the electric displacement field ${\displaystyle \mathbf {D} }$, given by

${\displaystyle \mathbf {D} =\mathbf {E} +4\pi \mathbf {P} }$.

If the external field is not strong enough to greatly change the properties of the dielectric medium, one can treat the induced polarization within what it called the linear response regime. Here the information on how the dielectric reacts to the external field is contained in the dielectric function, ${\displaystyle \epsilon _{\alpha \beta }}$

${\displaystyle \epsilon _{\alpha \beta }=\delta _{\alpha \beta }+4\pi {\frac {\partial P_{i}}{\partial E_{j}}}}$,

which leads to (assuming that the system has time-reversal)

${\displaystyle D_{\alpha }(\omega )=\epsilon _{\alpha \beta }(\omega )E_{\beta }(\omega )}$.

Depending on the nature of the external field there are different approaches for the calculation of ${\displaystyle \epsilon }$. If ${\displaystyle \mathbf {E} }$ is static, then one can rely on perturbative methods based on finite differences or Density functional perturbation Theory. However, if one wishes to reproduce results where a time-dependent ${\displaystyle \mathbf {E} }$-field was used (e.g. optical absorption, reflectance, magento-optical Kerr effect), such methods can no longer be employed as the response will depend on the frequency of the incident field. For these cases one must employ methods based on time-dependet linear response (e.g. Green-Kubo) or Many-body Perturbation Theory.

Bellow we present an overview of all possible cases where VASP employs either one of such methods for the calculation of ${\displaystyle \epsilon }$.

Methods for computing ${\displaystyle \epsilon }$

Static response: Density functional perturbation Theory (DFPT) and Finite differences based methods

LEPSILON

By setting LEPSILON=.True., VASP uses DFPT to compute the static ion-clamped dielectric matrix with or without local field effects. Derivatives are evaluated using Sternheimer equations, avoiding the explicit computation of derivatives of the periodic part of the wave function. This method does not require the inclusion of empty states via the NBANDS parameter.

At the end of the calculation the both the values of ${\displaystyle \epsilon }$ including (LRPA=.True.) or excluding (LRPA=.False.) local-field effects are printed in the OUTCAR file. Users can perform a consistency check by comparing the values with no local field to the zero frequency results for ${\displaystyle \epsilon }$ obtained from a calculation with LOPTICS=.True..

LCALCEPS

With LCALCEPS=.True., the dielectric tensor is computed from the derivative of the polarization, using

${\displaystyle \epsilon _{ij}^{\infty }=\delta _{ij}+{\frac {4\pi }{\epsilon _{0}}}{\frac {\partial P_{i}}{\partial {\mathcal {E}}_{j}}}\qquad {i,j=x,y,z}.}$

However, here the derivative is evaluated explicitly by employing finite-differences. The direction and intensity of the perturbing electric field has to be specified in the INCAR using the EFIELD_PEAD variable. As with the previous method, at the end of the calculation VASP will write the dielectric tensor in the OUTCAR file. Control over the inclusion of local-field effects is done with the variable LRPA.

Dynamical response: Green-Kubo and Many-body Perturbation theory Methods

LOPITCS

The variable LOPTICS allows for the calculation of the frequency dependent dielectric function once the ground state is computed. It uses the explicit expression to evaluate the imaginary part of ${\displaystyle \epsilon }$

${\displaystyle \epsilon _{\alpha \beta }^{(2)}\left(\omega \right)={\frac {4\pi ^{2}e^{2}}{\Omega }}\mathrm {lim} _{q\rightarrow 0}{\frac {1}{q^{2}}}\sum _{c,v,\mathbf {k} }2w_{\mathbf {k} }\delta (\epsilon _{c\mathbf {k} }-\epsilon _{v\mathbf {k} }-\omega )\times \langle u_{c\mathbf {k} +\mathbf {e} _{\alpha }q}|u_{v\mathbf {k} }\rangle \langle u_{v\mathbf {k} }|u_{c\mathbf {k} +\mathbf {e} _{\beta }q}\rangle ,}$

while the real part is evaluated using the Kramers-Kroning relation. At this level there are no effects coming from local fields.

This method requires a relatively large number of empty states, controlled by the variable NBANDS in the INCAR file and it should be checked for convergence.

Furthermore, the INCAR should also include values for CSHIFT (the broadening applied to the Lorentzian function which replaces the ${\displaystyle \delta }$-function), and NEDOS (the frequency grid for ${\displaystyle \omega }$).

ALGO = TDHF

This option performs a time-dependent Hartree-Fock or DFT calculation. It follows the Casida equation and uses a Fourier transform of the time-evolving dipoles to compute ${\displaystyle \epsilon }$.

The number of NBANDS controls how many bands are present in the time-evolution. This does not need to be as high as when LOPTICS is active, and only a few bands above the band gap need to be included.

The choice of time-dependent kernel is controlled by AEXX, HFSCREEN, and LFXC variables. For calculations using hybrid functionals, AEXX controls the fraction of exact exchange used in the exchange correlation potential, while HFSCREEN specifies the range-separation parameter. For a pure TDDFT calculation, LFXC uses the local exchange-correlation kernel in the time-evolution equations.

ALGO = TIMEEV

Uses a delta-pulse electric field to probe all transitions and calculate the dielectric function by following the evolution in time of the dipole momenta. This algorithm is able to fully reproduce the absorption spectra from standard Bethe-Salpeter calculations by setting the correct time-dependent kernel with LHARTREE=.True. and LFXC=.True.

The time step is controlled automatically by the CSHIFT and PREC variables. This means that the smaller the value of CSHIFT and the more accurate the level of precision chosen by the user, the higher the number of time steps that VASP will perform, and so the higher the cost of the calculation.

The number of valence and conduction bands involved in the time propagation are set by the NBANDSO and NBANDSV variables, respectively. Once again, users are advised to choose a small number of bands near the band gap if they wish to reproduce optical measurements.

Finally, the maximum energy used in both the Fourier transform and in calculating the frequency dependent dielectric function is set by OMEGAMAX, and the sampling of the frequency grid is controlled by NEDOS.

ALGO = CHI

Here the frequency dielectric function is computed at the independent particle level but starting from a GW calculation. VASP will compute the polarizability ${\displaystyle \chi }$ by setting ALGO=Chi in the INCAR file and then use

${\displaystyle \epsilon _{\mathbf {G} \mathbf {G} '}^{-1}(\mathbf {q} ,\omega )=\delta _{\mathbf {G} \mathbf {G} '}-v(\mathbf {q} +\mathbf {G} )\chi _{\mathbf {G} \mathbf {G} '}(\mathbf {q} ,\omega )}$

to compute the dielectric function. Here ${\displaystyle v}$ is the Coulomb potential. Effects coming from local fields in the Random-Phase Approximation (RPA) can be activated by setting LRPA=.True..

Two methods for computing the polarizability are available: if the user sets LSPECTRAL=.True., VASP will avoid direct computation of ${\displaystyle \chi }$ and use a fast matrix-vector product. However, this can introduce spurious peaks at low frequencies for some values of CSHIFT and NOMEGA. The second method computes ${\displaystyle \chi }$ directly and is activated by setting LSPECTRAL=.False.. However it is much slower than the former method.

ALGO = BSE

Setting ALGO=BSE computes the macroscopic dielectric function ${\displaystyle \epsilon _{M}}$ by solving the Bethe-Salpeter equations. Here the electron-hole pairs are treated a new quasi-particle, an exciton, and the dielectric function is built using the new eigenvectors (${\displaystyle X_{\lambda }^{cv\mathbf {k} }}$) and eigenvalues (${\displaystyle \omega _{\lambda }}$)

${\displaystyle \epsilon _{M}(\mathbf {q} ,\omega )=1+v(\mathbf {q} )\sum _{\lambda \lambda '}\sum _{c,v,\mathbf {k} }\sum _{c',v',\mathbf {k} '}\langle c\mathbf {k} |e^{i\mathbf {qr} }|v\mathbf {k} \rangle X_{\lambda }^{cv\mathbf {k} }\langle c'\mathbf {k'} |e^{-i\mathbf {qr} }|v'\mathbf {k'} \rangle X_{\lambda '}^{c'v'\mathbf {k} ',*}\times S_{\lambda ,\lambda '}^{-1}\left({\frac {1}{\omega _{\lambda }-\omega -i\delta }}+{\frac {1}{\omega _{\lambda }+\omega +i\delta }}\right)~.}$

where ${\displaystyle S_{\lambda \lambda '}}$ is the overlap between exciton states of indices ${\displaystyle \lambda }$ and ${\displaystyle \lambda '}$ (in general the BSE Hamiltonian is not hermitian, so eigenstates associated to different eigenvalues are not necessarily orthogonal).

The number of valence and conduction states which are included in the BSE Hamiltonian is controlled by the variables NBANDSO and NBANDSV, respectively. Note that normally only a few bands above and below the band gap are needed to converge the optical spectrum, so users should be careful in setting up these two variables. Otherwise the calculation might run out of memory.

For comparison with optical experiments (e.g. absorption, MOKE, reflectance), ${\displaystyle q}$ is the photon momentum and usually it is taken in the ${\displaystyle \mathbf {q} \to 0}$ limit. Furthermore, the coupling between the resonant and anti-resonant terms can be switched off, in what is called the Tamm-Dancoff approximation. This approximation can be activated with the variable ANTIRES set to 0. Setting this variable to 1 or 2 will include the coupling, but increase the computational cost.

Level of approximation

Micro-macro connection

It is important to distinguish between macroscopic quantities, measured over several repetitions of the unit cell, and microscopic quantities, which include fields that change rapidly in all regions of the unit cell.

When measuring a property, it is its macroscopic version that will be represented in the data. This means that in order to compare the microscopic quantities like the dielectric function, one must average it over several repetitions of the unit cell. It is possible to show that the macroscopic dielectric function, ${\displaystyle \epsilon _{M}(\mathbf {q} ,\omega )}$ is related to the microscopic one via

${\displaystyle \epsilon _{M}(\mathbf {q} ,\omega )={\frac {1}{\epsilon _{\mathbf {G} =0,\mathbf {G} '=0}^{-1}(\mathbf {q} ,\omega )}}}$

where ${\displaystyle \epsilon _{\mathbf {G} =0,\mathbf {G} '=0}^{-1}(\mathbf {q} ,\omega )}$ is the inverse dielectric function at ${\displaystyle \mathbf {G} =\mathbf {G} '=0}$.

Note that this does not mean that ${\displaystyle \epsilon _{M}(\mathbf {q} ,\omega )=\epsilon _{\mathbf {G} =0,\mathbf {G} '=0}(\mathbf {q} ,\omega )}$! The full matrix ${\displaystyle \epsilon _{\mathbf {G} ,\mathbf {G} '}(\mathbf {q} ,\omega )}$ has to be inverted and it is the component at ${\displaystyle \mathbf {G} =\mathbf {G} '=0}$ that is used to calculate ${\displaystyle \epsilon _{M}(\mathbf {q} ,\omega )}$.

Finite momentum dielectric function

In the optical limit the momentum of the incoming photon is take as almost zero, since the wavelength of the electric field is several times larger than the dimensions of the unit cell. In this case the evaluating the dielectric function as to be taken with the limit of ${\displaystyle \mathbf {q} \to 0}$, since the Coulomb potential diverges at very small momentum. Be it for the independent particle regime of full BSE, one can write that

${\displaystyle \lim _{\mathbf {q} \to 0}{\frac {\langle c\mathbf {k} +\mathbf {q} |e^{\mathrm {i} \mathbf {q} \cdot \mathbf {r} }|v\mathbf {k} \rangle }{q}}\approx \lim _{\mathbf {q} \to 0}{\frac {\langle c\mathbf {k} +\mathbf {q} |1+\mathrm {i} \mathbf {q} \cdot \mathbf {r} |v\mathbf {k} \rangle }{q}}={\hat {\mathbf {q} }}\cdot \langle c\mathbf {k} +\mathbf {q} |\mathbf {r} |v\mathbf {k} \rangle }$

However, in some cases it might be important to analyse the effects of finite momentum excitons (e.g. see the case of the optical absorption of bulk hexagonal BN). Users can calculate the absorption spectrum at finite momentum by using the variable KPOINT_BSE in the INCAR file. Users should provide the index of the q-point of interest, which can be found either in the OUTCAR file.

Local fields in the Hamiltonian

Local fields, i.e. terms with finite ${\displaystyle \mathbf {G} }$ can be turned on of off in the Coulomb potential when evaluating the polarizability. VASP will distinguish the results from both in the OUTCAR file with:

MACROSCOPIC STATIC DIELECTRIC TENSOR (including local field effects) 

BORN EFFECTIVE CHARGES (including local field effects) 

PIEZOELECTRIC TENSOR (including local field effects) 

and

MACROSCOPIC STATIC DIELECTRIC TENSOR (excluding local field effects) 

BORN EFFECTIVE CHARGES (excluding local field effects) 

PIEZOELECTRIC TENSOR (excluding local field effects)

for both cases.

Another approximation can also be taken, where the contributions from the exchange-correlation kernel are neglected when evaluating the polarizability. This is equivalent to the so called Random-Phase Approximation (RPA) and can be activated with by setting LRPA=.True. in the INCAR.

Ion-clamped vs relaxed/dressed dielectric function

The dielectric function computed either in the static or dynamical response regimes does not consider the effects coming from changes in atomic coordinates due to the incoming electric field. This can be corrected, however, by adding computing the ion-relaxed (or "dressed") dielectric function ${\displaystyle {\bar {\epsilon }}}$

${\displaystyle {\bar {\epsilon }}_{\alpha \beta }=\epsilon _{\alpha \beta }+\Omega _{0}^{-1}Z_{m\alpha }(\Phi ^{-1})_{mn}Z_{n\beta },}$

where ${\displaystyle \Omega _{0}}$ is the volume of the unit cell, ${\displaystyle Z_{n\alpha }}$ is the Born effective charge, and ${\displaystyle \Phi _{mn}}$ is the force constants matrix.

Density-density versus current-current response functions

The inclusion of an electromagnetic field in the Hamiltonian is subject to a gauge choice. For instance, a classical electric field can be described by either a scalar potential ${\displaystyle \phi }$ or a longitudinal vector potential ${\displaystyle \mathbf {A} }$ in the incomplete Weyl gauge (${\displaystyle \phi }$ = 0). In the first case this means that the in perturbing potential the potential will couple to the electronic density, while in the second the vector potential will couple to a current. The fundamental consequence is that there will be two different response functions defined in both cases: a density-density response function for the first, ${\displaystyle \chi _{\rho \rho }}$; and a current-current response function, ${\displaystyle \chi _{jj}}$.

More generally, perturbations associated to longitudinal fields will be described by the density-density polarisability function ${\displaystyle \chi _{\rho \rho }}$ (e.g. laser fields taken in the classical limit), while transverse fields will be described by the current-current polarizability ${\displaystyle \chi _{jj}}$ (which is in fact a 3x3 tensor). The reason behind this comes from the fact that the time-dependent density is associated only to the longitudinal part of the current, via the continuity equation. This also links both response functions via

${\displaystyle q^{2}\chi _{jj}(\mathbf {q} ,\omega )=\omega ^{2}\chi _{\rho \rho }(\mathbf {q} ,\omega ),}$

and also guarantees that at finite momentum and frequency the dielectric functions obtained from either approach much match, ${\displaystyle \epsilon [\chi _{\rho \rho }]=\epsilon [\chi _{jj}]}$.

The current-current based approach is exact at both ${\displaystyle \mathbf {q} \to 0}$ and at ${\displaystyle \mathbf {q} =0}$. Specially for metals, it will reproduce the proper behaviour of the Drude tail at ${\displaystyle \omega =0}$. However, ${\displaystyle \chi _{jj}}$ is more prone to numerical instabilities.

Relation to observables

Polarizability

Optical conductivity

From Maxwell's equations and the microscopic form of Ohm's law it is possible to arrive at the following relation between the dielectric function and the optical conductivity ${\displaystyle \sigma (\omega )}$ tensors

${\displaystyle \epsilon _{\alpha \beta }(\omega )=\delta _{\alpha \beta }+\mathrm {i} {\frac {4\pi }{\omega }}\sigma _{\alpha \beta }(\omega ).}$

Optical absorption

For an electromagnetic wave travelling through a medium, one can express the electric field as ${\displaystyle \mathbf {E} (\mathbf {r} ,t)=\mathbf {E} _{0}e^{-\mathrm {i} (\omega t-\mathbf {q} \cdot \mathbf {r} )}}$, and the effects of the medium in the wave propagation are contained inside the dispersion relation ${\displaystyle \omega =\omega (\mathbf {q} )}$. Using Maxwell's equations on can arrive at

${\displaystyle q^{2}={\frac {\omega ^{2}}{c^{2}}}\epsilon (\omega ).}$

If the magnetic permeability of the material is assumed to be equal to that of vacuum, the equation above means that the refractive index can be written as ${\displaystyle n={\sqrt {\epsilon (\omega )}}={\tilde {n}}+\mathrm {i} k}$. Since ${\displaystyle n}$ is now complex, this means that the exponential factor in ${\displaystyle \mathbf {E} (\mathbf {r} ,t)}$ will have a dampening factor, ${\displaystyle e^{-{\frac {\omega }{c}}k{\hat {q}}\cdot \mathbf {r} }}$, which accounts for the absorption of electromagnetic energy by the medium. With this relation one can define the absorption coefficient, ${\displaystyle \alpha (\omega )}$ as

${\displaystyle \alpha (\omega )={\frac {2\omega }{c}}k(\omega ).}$

Reflectance

From the previous subsection, one can also define the reflectivity coefficient at normal incidence as

${\displaystyle R={\frac {(1-{\tilde {n}})^{2}+k^{2}}{(1+{\tilde {n}})^{2}+k^{2}}}.}$

This equation can be generalized for any angle of incidence ${\displaystyle \theta }$, resulting in the general form of Fresnel equations.

MOKE

A particular case of interaction between the incident electronmagnetic wave and the medium, when the former induces a finite magnetization in the latter. Here the reflected wave will gain an extra complex phase with respect to the incident ${\displaystyle \mathbf {E} }$-field. For a two-dimensional material (e.g. hexagonal BN, MoS${\displaystyle _{2}}$) this phase can be computed from the off-diagonal components of the dielectric tensor using

${\displaystyle \theta _{\mathrm {K} }(\omega )=-\mathrm {Re} \left[{\frac {\epsilon _{xy}(\omega )}{(\epsilon _{xx}(\omega )-1){\sqrt {\epsilon _{xx}}}}}\right].}$

Combination with other perturbations

Atomic displacements

Strain

The dielectric tensor can also be included in the evaluation of the elastic tensor, ${\displaystyle C_{jk}}$, (see [Static_linear_response:_theory|here] for more information on derived quantities from static linear response). While this quantity is normally evaluated at fixed ${\displaystyle \mathbf {E} }$-field, in cases where a thin film is placed between layers of insulating materials it is more convenient to evaluate the elastic tensor for fixed displacement field ${\displaystyle \mathbf {D} }$, since the boundary conditions fix the components of this vector in the direction normal to the surface.

If ${\displaystyle C_{jk}^{E}}$ is the elastic tensor defined at fixed ${\displaystyle \mathbf {E} }$-field and ${\displaystyle C_{jk}^{D}}$ the elastic tensor defined at fixed ${\displaystyle \mathbf {D} }$-field, then they are related by

${\displaystyle C_{jk}^{D}=C_{jk}^{E}+e_{\alpha j}(\epsilon )_{\alpha \beta }^{-1}e_{\beta k},}$

where ${\displaystyle e_{\alpha j}}$ is the ion-relaxed piezoelectric tensor.