DFT and TDDFT Calculations of Ground and Excited States of Photoelectron Emission

Abstract

The Density-Functional Theory (DFT) is a reformulation of the quantum study of a correlated N-body system into a simpler system with independent equations being solved iteratively. The DFT considers only ground states of the systems. The extension to the time-dependent case of this theory is the Time-Dependent Density-Functional Theory (TDDFT) that also takes into account the excited states of the system. These calculations are very interesting in photonics areas. In fact, the interaction between electrons and light in the vicinity of solid surfaces and nanostructures is important as pathway to integrate photonics and electronics. The capability to couple light and electrons in purposefully designed device depends on the capability of creating such devices and the understanding of the underlying science.

1. Introduction

In the fields of photonics, researchers use the dielectric function and the electronic density in their calculations and investigations [1, 2, 3, 4, 5, 6, 7, 8] to determine the physical and optical characteristics of materials such as noble metals (Au and Ag). In fact, the determination of the electronic structure of a material can give us all the physical and optical information about it. DFT and TDDFT are used to determine the electronic structure of materials such as noble metals in their ground and excited states. The purpose of these calculations is to determine the electron density $\rho \left(\overrightarrow{r}\right)$. Indeed, it informs us about all the physical properties of the studied system. We simplify this very complex calculation by using several approximations and theorems. We propose at the end a flowchart to carry out a computer simulation allowing to highlight the values of the electron density $\rho \left(\overrightarrow{r}\right)$.

2. DFT calculation steps

The DFT is based on wave Schrodinger equation

$${\widehat{\mathrm{{\rm H}}}}_T{\psi }_{\left(\left\{{\overrightarrow{R}}_I\right\},\left\{\overrightarrow{r_i}\right\},t\right)}=i\frac{\partial }{\partial t}{\psi }_{\left(\left\{{\overrightarrow{R}}_I\right\},\left\{\overrightarrow{r_i}\right\},t\right)}$$E1

We consider that the system contains N cores and M electrons.

${\widehat{\mathrm{{\rm H}}}}_T$: is the total Hamiltonian representing N cores and M electrons.

${\psi }_{\left(\left\{{\overrightarrow{R}}_I\right\},\left\{\overrightarrow{r_i}\right\},t\right)}$: is the wave function representing N cores and M electrons.

$\left\{{\overrightarrow{R}}_I\right\}$and $\left\{\overrightarrow{r_i}\right\}$: represent respectively the set of nuclear and electronic coordinates.

In the case of stationary processes, the Schrodinger becomes:

$${\widehat{\mathrm{{\rm H}}}}_T{\psi }_{\left(\left\{{\overrightarrow{R}}_I\right\},\left\{\overrightarrow{r_i}\right\}\right)}=i\mathrm{E}{\psi }_{\left(\left\{{\overrightarrow{R}}_I\right\},\left\{\overrightarrow{r_i}\right\}\right)}$$E2

where E represent the energy of the system described by ${\psi }_{\left(\left\{{\overrightarrow{R}}_I\right\},\left\{\overrightarrow{r_i}\right\}\right)}$

$${\widehat{\mathrm{{\rm H}}}}_T={\widehat{\mathrm{T}}}_T+{\widehat{V}}_T$$E3

where ${\widehat{\mathrm{T}}}_T$ is the total kinetic energy of the system operator, and ${\widehat{V}}_T$ is the operator describing all Colombian interaction.

${\widehat{\mathrm{{\rm H}}}}_T$ can be written as:

$${\widehat{\mathrm{{\rm H}}}}_T={\widehat{\mathrm{T}}}_n+{\widehat{\mathrm{T}}}_e+{\widehat{V}}_{n-e}+{\widehat{V}}_{e-e}+{\widehat{V}}_{n-n}$$E4

where

${\widehat{\mathrm{T}}}_n=\frac{-{\mathrm{\hslash }}^2}{2}{\sum }_I\frac{{\nabla }^2{\overrightarrow{R}}_I}{M_n}$: is the kinetic energy of N cores with mass $M_n$.

${\widehat{\mathrm{T}}}_n=\frac{-{\mathrm{\hslash }}^2}{2}{\sum }_i\frac{{\nabla }^2{\overrightarrow{r}}_i}{m_e}$: is the kinetic energy of M electrons with mass $m_e$.

${\widehat{V}}_{n-e}=\frac{1}{4{\mathit{\pi \epsilon }}_0}{\sum }_{i,j}\frac{e^2{Z}_i}{\left|{\overrightarrow{R}}_i-{\overrightarrow{r}}_j\right|}$: is the core-electron attractive Colombian interaction.

${\widehat{V}}_{e-e}=\frac{1}{8{\mathit{\pi \epsilon }}_0}{\sum }_{i\ne j}\frac{e^2}{\left|{\overrightarrow{r}}_i-{\overrightarrow{r}}_j\right|}$: is the electron-electron repulsive Colombian interaction.

${\widehat{V}}_{n-n}=\frac{1}{8{\mathit{\pi \epsilon }}_0}{\sum }_{i\ne j}\frac{e^2{Z}_i{Z}_j}{\left|{\overrightarrow{R}}_i-{\overrightarrow{R}}_j\right|}$: is the core-core repulsive Colombian interaction.

Eq. (4) takes into account N cores and M electrons. Some simplifications must be made:

• First approximation level: Born-Oppenheimer approximation. (see paragraph 2-1)

• Second approximation level: Hartree-Fock approximation. (see paragraph 2-2)

• Third approximation level: Inherent approximation to the resolution of equation. (see paragraph 2-3)

2.1 Adiabatic born-Oppenheimer approximation

It offers the opportunity to treat separately the cores and the electrons.

So we have

$${\mathrm{\Psi }}_R\left(\overrightarrow{R},\overrightarrow{r}\right)={\mathrm{\Phi }}_R\left(\overrightarrow{R}\right){\mathrm{\Psi }}_R\left(\overrightarrow{r}\right)$$E5

where

${\mathrm{\Phi }}_R\left(\overrightarrow{R}\right)$: is the wave function describing the cores.

${\mathrm{\Psi }}_R\left(\overrightarrow{r}\right)$: is the wave function describing the electrons.

In this approximation, the interaction electron-phonon is neglected.

Using Eqs. (3) and (4), the new Hamiltonian becomes

$${\widehat{\mathrm{{\rm H}}}}_T={\widehat{\mathrm{T}}}_e+{\widehat{V}}_{n-e}+{\widehat{V}}_{e-e}+\left({\widehat{V}}_{n-n}={\widehat{V}}_{\mathit{ext}}=\mathit{Cst}\right)$$E6

so

$$\begin{array}{ll}{\widehat{\mathrm{{\rm H}}}}_T=& \frac{-{\mathrm{\hslash }}^2}{2}\sum _i\frac{{\nabla }^2{\overrightarrow{r}}_i}{m_e}+\frac{1}{4{\mathit{\pi \epsilon }}_0}\sum _{i,j}\frac{e^2{Z}_i}{\left|{\overrightarrow{R}}_i-{\overrightarrow{r}}_j\right|}+\frac{1}{8{\mathit{\pi \epsilon }}_0}\sum _{i\ne j}\frac{e^2}{\left|{\overrightarrow{r}}_i-{\overrightarrow{r}}_j\right|}\\ & +\left(\frac{1}{8{\mathit{\pi \epsilon }}_0}\sum _{i\ne j}\frac{e^2{Z}_i{Z}_j}{\left|{\overrightarrow{R}}_i-{\overrightarrow{R}}_j\right|}=\mathit{Cst}\right)\end{array}$$E7

The nuclear kinetic energy term independent from electrons is canceled (equal to zero). The attractive potential energy (electron-core) term is removed and the repulsive potential energy (core-core) becomes a constant simply evaluated for determined geometry. Using the Hohenberg and Kohn Theorems (HK theorems) [9], the formulation of the Schrödinger equation can be now based on the electron density $\rho \left(\overrightarrow{r}\right)$. This is due to the two HK theorems. In fact, the first KH theorem indicated that the total energy of the system in the ground state is a single universal function of the electronic density.

$$E=E\left[\rho \left(\overrightarrow{r}\right)\right]$$E8

HK gives

$$E\left[\rho \left(\overrightarrow{r}\right)\right]=F_{\mathit{HK}}\left[\rho \left(\overrightarrow{r}\right)\right]+\int {\widehat{V}}_{\mathit{ext}}\left(\overrightarrow{r}\right)\rho \left(\overrightarrow{r}\right)d\overrightarrow{r},$$E9

where $F_{\mathit{HK}}\left[\rho \left(\overrightarrow{r}\right)\right]$ represent the HK universal functional and ${\overrightarrow{V}}_{\mathit{ext}}\left(\overrightarrow{r}\right)$ represent the external potential. The second HK theorem says that it is an analogue variational principle to the originally proposed in the approach of Hartree-Fock for the functional of the wave function ($\frac{\partial E\left[\psi \right]}{\partial \psi }=0$), but this time applied to the electronic density functional.

$${\left.\frac{\partial E\left[\rho \left(\overrightarrow{r}\right)\right]}{\partial \rho \left(\overrightarrow{r}\right)}\right|}_{{\rho }_0\left(\overrightarrow{r}\right)}=0$$E10

where ${\rho }_0\left(\overrightarrow{r}\right)$ is the exact electronic density of the system in the ground state. This theorem can be reformulated like: For a potential ${\widehat{V}}_{\mathit{ext}}\left(\overrightarrow{r}\right)$ and with a number M of electrons, the total energy of the system reaches its minimum value when $\rho \left(\overrightarrow{r}\right)$ is corresponding to the exact density ${\rho }_0\left(\overrightarrow{r}\right)$ of the ground state.

Considering the two HK theorems, the resolution of Schrödinger equation consists in seeking the minimization of $E\left[\rho \left(\overrightarrow{r}\right)\right]$: $\frac{\partial E\left[\rho \left(\overrightarrow{r}\right)\right]}{\partial \rho \left(\overrightarrow{r}\right)}=0$ by applying the constraint of the conservation of the total number of particles: $\int \rho \left(\overrightarrow{r}\right)d\overrightarrow{r}=M$.

This problem can be resolved using Lagrange multipliers

$$G\left[\rho \left(\overrightarrow{r}\right)\right]=\int \rho \left(\overrightarrow{r}\right)d\overrightarrow{r}-M$$E11

The conservation of the total number of particles constraint becomes $G\left[\rho \left(\overrightarrow{r}\right)\right]=0$, and by introducing auxiliary function $A\left[\rho \left(\overrightarrow{r}\right)\right]$ with:

$$A\left[\rho \left(\overrightarrow{r}\right)\right]=E\left[\rho \left(\overrightarrow{r}\right)\right]-\mathit{\mu G}\left[\rho \left(\overrightarrow{r}\right)\right]$$E12

where $\mu$ is the Lagrange multiplier, the problem to resolve becomes

$$\partial A\left[\rho \left(\overrightarrow{r}\right)\right]=\int \frac{\partial A\left[\rho \left(\overrightarrow{r}\right)\right]}{\partial \rho \left(\overrightarrow{r}\right)}\partial \mathit{\rho dr}=0$$E13

$$\partial \left\{E\left[\rho \left(\overrightarrow{r}\right)\right]-\mu \left[\int \rho \left(\overrightarrow{r}\right)d\overrightarrow{r}-M\right]\right\}=0$$

We must calculate the functional derivative of $A\left[\rho \left(\overrightarrow{r}\right)\right]$:

$$\frac{\partial A\left[\rho \left(\overrightarrow{r}\right)\right]}{\partial \rho \left(\overrightarrow{r}\right)}=\frac{\partial }{\partial \rho \left(\overrightarrow{r}\right)}\left\{E\left[\rho \left(\overrightarrow{r}\right)\right]-\mu \left[\int \rho \left(\overrightarrow{r}\right)d\overrightarrow{r}-M\right]\right\}=\frac{\partial E\left[\rho \left(\overrightarrow{r}\right)\right]}{\partial \rho \left(\overrightarrow{r}\right)}-\mu \frac{\partial }{\partial \rho \left(\overrightarrow{r}\right)}\left[\int \rho \left(\overrightarrow{r}\right)d\overrightarrow{r}\right]$$

so we have

$$\frac{\partial A\left[\rho \left(\overrightarrow{r}\right)\right]}{\partial \rho \left(\overrightarrow{r}\right)}=\frac{\partial E\left[\rho \left(\overrightarrow{r}\right)\right]}{\partial \rho \left(\overrightarrow{r}\right)}-\mu$$E14

By replacing this expression in the equation of $\partial A\left[\rho \left(\overrightarrow{r}\right)\right]$, we obtained

$$\partial A\left[\rho \left(\overrightarrow{r}\right)\right]=\int \left[\frac{\partial E\left[\rho \left(\overrightarrow{r}\right)\right]}{\partial \rho \left(\overrightarrow{r}\right)}-\mu \right]\partial \mathit{\rho dr}=0$$

$$\int \frac{\partial E\left[\rho \left(\overrightarrow{r}\right)\right]}{\partial \rho \left(\overrightarrow{r}\right)}\partial \mathit{\rho dr}=\int \mu \partial \mathit{\rho dr}$$E15

$$\int \frac{\partial E\left[\rho \left(\overrightarrow{r}\right)\right]}{\partial \rho \left(\overrightarrow{r}\right)}\partial \mathit{\rho dr}=\mu$$

Using Eq. (9) and calculating the functional derivative of $E\left[\rho \left(\overrightarrow{r}\right)\right]$ we obtained

$$\frac{\partial E\left[\rho \left(\overrightarrow{r}\right)\right]}{\partial \rho \left(\overrightarrow{r}\right)}={\widehat{V}}_{\mathit{ext}}\left(\overrightarrow{r}\right)+\frac{\partial F_{\mathit{HK}}\left[\rho \left(\overrightarrow{r}\right)\right]}{\partial \rho \left(\overrightarrow{r}\right)}$$E16

By replacing Eq. (16) in Eq. (15) we obtained

$$\mu =\frac{\partial E\left[\rho \left(\overrightarrow{r}\right)\right]}{\partial \rho \left(\overrightarrow{r}\right)}={\widehat{V}}_{\mathit{ext}}\left(\overrightarrow{r}\right)+\frac{\partial F_{\mathit{HK}}\left[\rho \left(\overrightarrow{r}\right)\right]}{\partial \rho \left(\overrightarrow{r}\right)}$$E17

This equation constitutes the fundamental DFT formalism.

It only remains to determine the expression of the unknown function $F_{\mathit{HK}}\left[\rho \left(\overrightarrow{r}\right)\right]$ in Eq. (17), that’s why we are brought to use the Kohn-Sham approximations.

2.2 Kohn-Sham (KS) approximations (equations)

This step [10] consists of two approximations to transform the theorems of Hohenberg-Kohn into a practical workable from. However, the real system studied is redefined as a fictitious fermions system without interaction and with the same electronic density $\rho \left(\overrightarrow{r}\right)$ characterizing the real system to being up the terms of inter-electronic as corrections of the other terms. Also, single particle orbitals are induced to treat the kinetic energy term of electrons more precisely than under the Thomas-Fermi theory.

2.2.1 First approximation of KS

The transformation of the real interactive system as a fictitious non-interactive system. Considering the first HK theorem, the functional $E\left[\rho \left(\overrightarrow{r}\right)\right]$ can be written as indicated in Eq. (9):

$$E\left[\rho \left(\overrightarrow{r}\right)\right]=F_{\mathit{HK}}\left[\rho \left(\overrightarrow{r}\right)\right]+\int {\widehat{V}}_{\mathit{ext}}\left(\overrightarrow{r}\right)\rho \left(\overrightarrow{r}\right)d\overrightarrow{r}$$E18

The functional $F_{\mathit{HK}}\left[\rho \left(\overrightarrow{r}\right)\right]$ is independent from ${\widehat{V}}_{\mathit{ext}}\left(\overrightarrow{r}\right)$ (external potential), and it is valid and applies regardless of the system studied. It contains a component of kinetic energy of electron ${\widehat{T}}_e\left[\rho \left(\overrightarrow{r}\right)\right]$, and another component corresponding to the mutual coulomb interaction between electrons ${\widehat{V}}_{e-e}\left[\rho \left(\overrightarrow{r}\right)\right]$. The minimization of this functional with the constraint of preservation of the number of particles M:

$\int \rho \left(\overrightarrow{r}\right)d\overrightarrow{r}=M$gives us directly the total energy of the system and the charge density of ground state from which all other physical properties can be extracted:

$$F_{\mathit{HK}}\left[\rho \left(\overrightarrow{r}\right)\right]=T_S\left[\rho \left(\overrightarrow{r}\right)\right]+E_H\left[\rho \left(\overrightarrow{r}\right)\right]+E_{\mathit{XC}}\left[\rho \left(\overrightarrow{r}\right)\right]+V_{\mathit{ext}}\left[\rho \left(\overrightarrow{r}\right)\right]$$E19

Where

$T_S\left[\rho \left(\overrightarrow{r}\right)\right]$: is the kinetic energy of non-interactive electron gas.

$E_H\left[\rho \left(\overrightarrow{r}\right)\right]=\frac{1}{2}\iint \frac{\rho \left(\overrightarrow{r}\right)\rho \left({\overrightarrow{r}}^{\text{'}}\right)}{\left|\overrightarrow{r}-{\overrightarrow{r}}^{\text{'}}\right|}d\overrightarrow{r}d{\overrightarrow{r}}^{\text{'}}$: is the chemical coulomb interaction between electrons (Hartree term).

$E_{\mathit{XC}}\left[\rho \left(\overrightarrow{r}\right)\right]$: the additional functional describing interactions inter-electronic not obtained from non-interactive system.

$V_{\mathit{ext}}\left[\rho \left(\overrightarrow{r}\right)\right]=\int {\widehat{V}}_{\mathit{ext}}\left(\overrightarrow{r}\right)\rho \left(\overrightarrow{r}\right)d\overrightarrow{r}$: is the external potential.

The functional $E_{\mathit{XC}}\left[\rho \left(\overrightarrow{r}\right)\right]$ is called “exchange-correlation energy,” this term contain all differences between fictitious non-interactive system and the real interactive system.

$$E_{\mathit{XC}}\left[\rho \left(\overrightarrow{r}\right)\right]=\left(T\left[\rho \left(\overrightarrow{r}\right)\right]-T_S\left[\rho \left(\overrightarrow{r}\right)\right]\right)+\left(V_{e-e}\left[\rho \left(\overrightarrow{r}\right)\right]-V_H\left[\rho \left(\overrightarrow{r}\right)\right]\right)$$E20

where:

$T\left[\rho \left(\overrightarrow{r}\right)\right]$: is the real kinetic energy (of real interactive system).

$T_S\left[\rho \left(\overrightarrow{r}\right)\right]$: is the non-interactive fermions system energy of KS.

$T\left[\rho \left(\overrightarrow{r}\right)\right]-T_S\left[\rho \left(\overrightarrow{r}\right)\right]$: this difference is low and is neglected.

$E_{\mathit{XC}}\left[\rho \left(\overrightarrow{r}\right)\right]$ traduced only the difference between Colombian energy of the real system $V_{e-e}\left[\rho \left(\overrightarrow{r}\right)\right]$ and the Colombian energy of non-interactive fermions system $V_H\left[\rho \left(\overrightarrow{r}\right)\right]$ of KS.

$$E_{\mathit{XC}}\left[\rho \left(\overrightarrow{r}\right)\right]\cong V_{e-e}\left[\rho \left(\overrightarrow{r}\right)\right]-V_H\left[\rho \left(\overrightarrow{r}\right)\right]$$E21

2.2.2 Second approximation of KS

This second approximation is based on the formulation of kinetic energy using an orbital approach. The exact formulation of kinetic energy T for the ground state systems is given by:

$$T=\sum _i^M{n}_{i\left.⟨{\phi }_i^{\ast }\right|}\frac{-{\mathrm{\hslash }}^2{\nabla }^2\overrightarrow{r}}{2m_e}\left|{\phi }_i\right.⟩$$E22

where ${\phi }_i$ are natural spin orbital’s and $n_i$ their respective occupation number according to the Pauli principle; $0\le n_i\le 1$.

KS have taken advantage of the fictitious non-interactive fermions system to describe the kinetic energy which is also according to the first theorem of KH as functional of the density:

$$T_s\left[\rho \left(\overrightarrow{r}\right)\right]=\sum _i^M{n}_{i\left.⟨{\phi }_i^{\ast }\right|}\frac{-{\mathrm{\hslash }}^2{\nabla }^2\overrightarrow{r}}{2m_e}\left|{\phi }_i\right.⟩$$E23

T_s is the non-interactive fermions system energy of KS as function as electron density $\rho \left(\overrightarrow{r}\right)$.

3. Kohn-Sham equations

The fundamental equation of DFT is expressed by applying variational KH principal [9] as:

$$\mu ={\widehat{V}}_{\mathit{eff}}\left(\overrightarrow{r}\right)+\frac{\partial T_s\left[\rho \left(\overrightarrow{r}\right)\right]}{\partial \rho \left(\overrightarrow{r}\right)}$$E24

where ${\widehat{V}}_{\mathit{eff}}\left(\overrightarrow{r}\right)$ is formulated as a functional of electronic density.

$${\widehat{V}}_{\mathit{eff}}\left(\overrightarrow{r}\right)={\widehat{V}}_{\mathit{eff}}\left[\rho \left(\overrightarrow{r}\right)\right]={\widehat{V}}_{\mathit{ext}}\left(\overrightarrow{r}\right)+\frac{\partial E_H\left[\rho \left(\overrightarrow{r}\right)\right]}{\partial \rho \left(\overrightarrow{r}\right)}+\frac{\partial E_{\mathit{XC}}\left[\rho \left(\overrightarrow{r}\right)\right]}{\partial \rho \left(\overrightarrow{r}\right)}$$

$${\widehat{V}}_{\mathit{eff}}\left(\overrightarrow{r}\right)={\widehat{V}}_{\mathit{ext}}\left(\overrightarrow{r}\right)+\int \frac{\rho \left({\overrightarrow{r}}^{\text{'}}\right)}{\left|\overrightarrow{r}-{\overrightarrow{r}}^{\text{'}}\right|}d{\overrightarrow{r}}^{\text{'}}+{\widehat{V}}_{\mathit{xc}}\left(\overrightarrow{r}\right)$$E25

where ${\widehat{V}}_{\mathit{xc}}\left(\overrightarrow{r}\right)$ is the exchange-correlation potential functional derivative of $E_{\mathit{XC}}\left[\rho \left(\overrightarrow{r}\right)\right]$ by $\rho \left(\overrightarrow{r}\right)$:

$${\widehat{V}}_{\mathit{xc}}\left(\overrightarrow{r}\right)=\frac{\partial E_{\mathit{XC}}\left[\rho \left(\overrightarrow{r}\right)\right]}{\partial \rho \left(\overrightarrow{r}\right)}\text{and}\rho \left(\overrightarrow{r}\right)=\sum _{i=1}^M{\left|{\phi }_i\left(\overrightarrow{r}\right)\right|}^2$$E26

The Schrödinger equation to resolve in the KS approach is:

$$\left[\frac{-{\mathrm{\hslash }}^2{\nabla }^2\overrightarrow{r}}{2m_e}+{\widehat{V}}_{\mathit{eff}}\left(\overrightarrow{r}\right)\right]\left|{\phi }_i\left(\overrightarrow{r}\right)\right.⟩={\epsilon }_i\left|{\phi }_i\left(\overrightarrow{r}\right)\right.⟩$$E27

The KS equations (Eqs. (25) and (27)) must be resolved with self-coherent method by starting with an initial electronic density. ${\widehat{V}}_{\mathit{eff}}\left(\overrightarrow{r}\right)$ is obtained for which the Schrodinger equation of Kohn-Sham (Eq. 27) is resolved, and a new electronic density is then calculated. From this new electronic density, a new ${\widehat{V}}_{\mathit{eff}}\left(\overrightarrow{r}\right)$ is determined. This process is repeated with self-coherent method up to that convergence is reached. The new electronic density obtained must be very close to the previous one corresponding to the criterion of convergence fixed before.

At this step, we only need to find the expression of $E_{\mathit{XC}}\left[\rho \left(\overrightarrow{r}\right)\right]$.Some exchange-correlation functional can be considered using different functional approximation families like: Local Density Approximation (LDA) [11], Generalized Gradient Approximation (CGA) [12], Meta Generalized Gradient Approximation (Meta CGA) [13], and Hybrid Functional [14].

Functional Family	Dependence
LDA	$\rho \left(\overrightarrow{r}\right)$
GGA	$\left|\nabla \rho \left(\overrightarrow{r}\right)\right|$, $\rho \left(\overrightarrow{r}\right)$
Meta-GGA	$\left|\nabla \rho \left(\overrightarrow{r}\right)\right|$, ${\nabla }^2\rho \left(\overrightarrow{r}\right)$, $\rho \left(\overrightarrow{r}\right)$
Hybrid	Exact Exchange, $\left|\nabla \rho \left(\overrightarrow{r}\right)\right|$, $\rho \left(\overrightarrow{r}\right)$

The order of accuracy is increasing from top to bottom. In the formalism of exchange-correlation functional, $E_{\mathit{XC}}\left[\rho \left(\overrightarrow{r}\right)\right]$ is presented like an interaction between the electronic density $\rho \left(\overrightarrow{r}\right)$ and density energy depending on $\rho \left(\overrightarrow{r}\right)$: ${\epsilon }_{\mathit{XC}}\left[\rho \left(\overrightarrow{r}\right)\right]$ with

$$E_{\mathit{XC}}\left[\rho \left(\overrightarrow{r}\right)\right]=\int {\epsilon }_{\mathit{XC}}\left[\rho \left(\overrightarrow{r}\right)\right]\rho \left(\overrightarrow{r}\right)d\overrightarrow{r}$$E28

${\epsilon }_{\mathit{XC}}\left[\rho \left(\overrightarrow{r}\right)\right]$ is considered as a summation of the contribution of exchange and correlation, (with) ${\epsilon }_{\mathit{XC}}\left[\rho \left(\overrightarrow{r}\right)\right]={\epsilon }_X\left[\rho \left(\overrightarrow{r}\right)\right]+{\epsilon }_C\left[\rho \left(\overrightarrow{r}\right)\right]$and

$$E_{\mathit{XC}}\left[\rho \left(\overrightarrow{r}\right)\right]=E_X\left[\rho \left(\overrightarrow{r}\right)\right]+E_C\left[\rho \left(\overrightarrow{r}\right)\right]=\int {\epsilon }_X\left[\rho \left(\overrightarrow{r}\right)\right]\rho \left(\overrightarrow{r}\right)d\overrightarrow{r}+\int {\epsilon }_C\left[\rho \left(\overrightarrow{r}\right)\right]\rho \left(\overrightarrow{r}\right)d\overrightarrow{r}$$E29

3.1 Local Density Approximation (LDA)

In this work we will only use the LDA approximation in our calculations. In this approximation, the electronic density can be treated locally in the form of uniform electrons of gas.

In other words, this approach is to perform the following two hypotheses:

• The exchange-correlation is dominated by the density at the point $\overrightarrow{r}$.

• The density $\rho \left(\overrightarrow{r}\right)$ is a function slowly varying with $\overrightarrow{r}$.

The fundamental hypothesis contained in the LDA formalism is to consider that the contribution of $E_{\mathit{XC}}\left[\rho \left(\overrightarrow{r}\right)\right]$ to the total energy of the system can be added by cumulated method of each portion of the non-uniform gas as it was locally uniform.

$$E_{\mathit{XC}}^{\mathit{LDA}}\left[\rho \left(\overrightarrow{r}\right)\right]=\int {\epsilon }_{\mathit{XC}}^{\mathit{LDA}}\left[\rho \left(\overrightarrow{r}\right)\right]\rho \left(\overrightarrow{r}\right)d\overrightarrow{r}$$E30

Where ${\epsilon }_{\mathit{XC}}^{\mathit{LDA}}\left[\rho \left(\overrightarrow{r}\right)\right]$ represent the exchange-correlation energy by electron in a system with electron in mutual correlation with uniform density $\rho \left(\overrightarrow{r}\right)$.

Using ${\epsilon }_{\mathit{XC}}^{\mathit{LDA}}\left[\rho \left(\overrightarrow{r}\right)\right]$, the exchange-correlation potential $V_{\mathit{XC}}^{\mathit{LDA}}\left(\overrightarrow{r}\right)$ can be obtained like:

$$V_{\mathit{XC}}^{\mathit{LDA}}\left(\overrightarrow{r}\right)=\frac{\partial \left(\rho \left(\overrightarrow{r}\right){\epsilon }_{\mathit{XC}}^{\mathit{LDA}}\left[\rho \left(\overrightarrow{r}\right)\right]\right)}{\partial \rho \left(\overrightarrow{r}\right)}$$E31

There are several forms for the term of exchange and correlation of a homogeneous electron gas, among others those of Kohn and Sham [10], Wigner [15], Perdew and Wang [16], and Hedin and Lundqvis [17], Ceperly and Alder [18]. The last one is the most used.

3.2 Approximation of Ceperley and Alder

Ceperley and Alder [18] used the exchange energy of uniform electrons gas given by Dirac formula as:

$$E_X\left[\rho \left(\overrightarrow{r}\right)\right]=-C_X\int {\rho }^{\frac{4}{3}}\left(\overrightarrow{r}\right)d\overrightarrow{r}$$E32

$${\epsilon }_X^{\mathit{LDA}}\left[\rho \left(\overrightarrow{r}\right)\right]=-C_X\int {\rho }^{\frac{1}{3}}\left(\overrightarrow{r}\right)d\overrightarrow{r}$$E33

with $C_X=\frac{3}{4}{\left(\frac{3}{\pi }\right)}^{\frac{1}{3}}$, ${\epsilon }_X^{\mathit{LDA}}\left[\rho \left(\overrightarrow{r}\right)\right]$ can be expressed like ${\epsilon }_X^{\mathit{LDA}}\left[\rho \left(\overrightarrow{r}\right)\right]=-0.738694{\rho }^{\frac{1}{3}}\left(\overrightarrow{r}\right)$.

The correlation energy ${\epsilon }_C^{\mathit{LDA}}\left[\rho \left(\overrightarrow{r}\right)\right]$ is derived on the second order of Moller-Plesset perturbation theory [19, 20, 21]:

$${\epsilon }_C^{\mathit{LDA}}\left[\rho \left(\overrightarrow{r}\right)\right]=\mathit{aLn}\left(1+\frac{b}{r_s}+\frac{b}{r_s^2}\right)$$E34

With $a=\frac{\mathit{Ln}2-1}{2{\pi }^2}=-0.01556111$ and $b=20.4562557$.

$r_s$is the Wigner-Seitz density parameter, in atomic unit we have $r_s={\left(\frac{4\mathit{\pi \rho }\left(\overrightarrow{r}\right)}{3}\right)}^{-1/3}$.

Considering the correlation term ${\epsilon }_C^{\mathit{LDA}}\left[\rho \left(\overrightarrow{r}\right)\right]$, no explicit analytical expression is known. Several different parameterizations have been proposed since the early 1970’s, and the most accurate results are based on quantum Monte Carlo simulations of Ceperley and Alder [18]. The most approximation commonly used today is that of Zunger [22] and who parameterized the Ceperley and Alder correlation functional for the spin-polarized electron gas and non-spin-polarized electron gas by the equation:

$${\epsilon }_C^{\mathit{PZ}}\left[\rho \left(\overrightarrow{r}\right)\right]=\left\{\begin{array}{l}A\times {\mathit{Lnr}}_s+B+C\times r_s\times {\mathit{Lnr}}_s+D\times r_s,r_s\le 1\\ \gamma /\left(1+{\beta }_1\times r_s+{\beta }_2\times r_s\right),r_s<1\end{array}\right.$$E35

For $r_s\le 1$, $A^U=0.0311$, $B^U=-0.048$, $C^U=0.002$ and $D^U=-0.0116$.

For $r_s<1$, $\gamma =-0.1423$, ${\beta }_1=1.0529$, and ${\beta }_2=0.3334$.

The improvements of the LDA approach must consider the gas of electrons in its real form (non-uniform and non-local), GGA, meta-GGA, and hybrid functional allow gradually approach to make in consideration their two effects.

4. Resolution of the KS Schrödinger equation

Combining between Eqs. (25) and (27) and using the expression $\rho \left(\overrightarrow{r}\right)=\sum _{i=1}^M{\left|{\phi }_i\left(\overrightarrow{r}\right)\right|}^2$, where M is the number of electrons, Eq. (27) becomes:

$$\left[\frac{-{\mathrm{\hslash }}^2{\nabla }^2\overrightarrow{r}}{2m_e}+{\widehat{V}}_{\mathit{ext}}\left(\overrightarrow{r}\right)+\int \frac{\rho \left({\overrightarrow{r}}^{\text{'}}\right)}{\left|\overrightarrow{r}-{\overrightarrow{r}}^{\text{'}}\right|}d{\overrightarrow{r}}^{\text{'}}+{\widehat{V}}_{\mathit{xc}}\left(\overrightarrow{r}\right)\right]\left|{\phi }_i\left(\overrightarrow{r}\right)\right.⟩={\epsilon }_i\left|{\phi }_i\left(\overrightarrow{r}\right)\right.⟩$$E36

The only unknown value in this expression is ${\widehat{V}}_{\mathit{xc}}\left(\overrightarrow{r}\right)$ (exchange-correlation potential).

Using LDA approximation of ${\widehat{V}}_{\mathit{xc}}\left(\overrightarrow{r}\right)$ (Eq. (30)),Where ${\epsilon }_{\mathit{xc}}^{\mathit{LDA}}\left[\rho \left(\overrightarrow{r}\right)\right]={\epsilon }_x^{\mathit{LDA}}\left[\rho \left(\overrightarrow{r}\right)\right]+{\epsilon }_c^{\mathit{LDA}}\left[\rho \left(\overrightarrow{r}\right)\right]$ with

${\epsilon }_{\mathit{xc}}^{\mathit{LDA}}\left[\rho \left(\overrightarrow{r}\right)\right]$: The exchange-correlation energy of one electron in a system of electron on mutual interaction with electronic density $\rho \left(\overrightarrow{r}\right)$.

${\epsilon }_x^{\mathit{LDA}}\left[\rho \left(\overrightarrow{r}\right)\right]$: The exchange energy of one electron.

${\epsilon }_c^{\mathit{LDA}}\left[\rho \left(\overrightarrow{r}\right)\right]$: The correlation energy of one electron.

Eq. (36) becomes:

$$\left[\frac{-{\mathrm{\hslash }}^2{\nabla }^2\overrightarrow{r}}{2m_e}+{\widehat{V}}_{\mathit{ext}}\left(\overrightarrow{r}\right)+\int \frac{\rho \left({\overrightarrow{r}}^{\text{'}}\right)}{\left|\overrightarrow{r}-{\overrightarrow{r}}^{\text{'}}\right|}d{\overrightarrow{r}}^{\text{'}}+\frac{\partial \left(\rho \left(\overrightarrow{r}\right){\epsilon }_{\mathit{xc}}^{\mathit{LDA}}\left[\rho \left(\overrightarrow{r}\right)\right]\right)}{\partial \rho \left(\overrightarrow{r}\right)}\right]\left|{\phi }_i\left(\overrightarrow{r}\right)\right.⟩={\epsilon }_i\left|{\phi }_i\left(\overrightarrow{r}\right)\right.⟩$$E37

$$\begin{array}{l}\left[\frac{-{\mathrm{\hslash }}^2{\nabla }^2\overrightarrow{r}}{2m_e}+{\widehat{V}}_{\mathit{ext}}\left(\overrightarrow{r}\right)+\int \frac{\rho \left({\overrightarrow{r}}^{\text{'}}\right)}{\left|\overrightarrow{r}-{\overrightarrow{r}}^{\text{'}}\right|}d{\overrightarrow{r}}^{\text{'}}+\frac{\partial \left(\rho \left(\overrightarrow{r}\right)\left({\epsilon }_x^{\mathit{LDA}}\left[\rho \left(\overrightarrow{r}\right)\right]+{\epsilon }_c^{\mathit{LDA}}\left[\rho \left(\overrightarrow{r}\right)\right]\right)\right)}{\partial \rho \left(\overrightarrow{r}\right)}\right]\left|{\phi }_i\left(\overrightarrow{r}\right)\right.⟩\\ ={\epsilon }_i\left|{\phi }_i\left(\overrightarrow{r}\right)\right.⟩\end{array}$$E38

$$\begin{array}{l}[\frac{-{\mathrm{\hslash }}^2{\nabla }^2\overrightarrow{r}}{2m_e}+{\widehat{V}}_{\mathit{ext}}\left(\overrightarrow{r}\right)+\int \frac{\rho \left({\overrightarrow{r}}^{\text{'}}\right)}{\left|\overrightarrow{r}-{\overrightarrow{r}}^{\text{'}}\right|}d{\overrightarrow{r}}^{\text{'}}+\frac{\partial \left(\rho \left(\overrightarrow{r}\right)\left({\epsilon }_x^{\mathit{LDA}}\left[\rho \left(\overrightarrow{r}\right)\right]\right)\right)}{\partial \rho \left(\overrightarrow{r}\right)}\\ +\frac{\partial \left(\rho \left(\overrightarrow{r}\right)\left({\epsilon }_c^{\mathit{LDA}}\left[\rho \left(\overrightarrow{r}\right)\right]\right)\right)}{\partial \rho \left(\overrightarrow{r}\right)}]\left|{\phi }_i\left(\overrightarrow{r}\right)\right.⟩\\ ={\epsilon }_i\left|{\phi }_i\left(\overrightarrow{r}\right)\right.⟩\end{array}$$E39

where

$$\begin{array}{ll}{\widehat{H}}_{\mathit{KS}}=& -\frac{1}{2}{\nabla }^2\overrightarrow{r}+{\widehat{V}}_{\mathit{ext}}\left(\overrightarrow{r}\right)+\int \frac{\rho \left({\overrightarrow{r}}^{\text{'}}\right)}{\left|\overrightarrow{r}-{\overrightarrow{r}}^{\text{'}}\right|}d{\overrightarrow{r}}^{\text{'}}+\frac{\partial \left(\rho \left(\overrightarrow{r}\right)\left({\epsilon }_x^{\mathit{LDA}}\left[\rho \left(\overrightarrow{r}\right)\right]\right)\right)}{\partial \rho \left(\overrightarrow{r}\right)}\\ & +\frac{\partial \left(\rho \left(\overrightarrow{r}\right)\left({\epsilon }_c^{\mathit{LDA}}\left[\rho \left(\overrightarrow{r}\right)\right]\right)\right)}{\partial \rho \left(\overrightarrow{r}\right)}\end{array}$$E40

$$\begin{array}{ll}{\widehat{V}}_{\mathit{eff}}\left(\overrightarrow{r}\right)=& {\widehat{V}}_{\mathit{ext}}\left(\overrightarrow{r}\right)+\int \frac{\rho \left({\overrightarrow{r}}^{\text{'}}\right)}{\left|\overrightarrow{r}-{\overrightarrow{r}}^{\text{'}}\right|}d{\overrightarrow{r}}^{\text{'}}+\frac{\partial \left(\rho \left(\overrightarrow{r}\right)\left({\epsilon }_x^{\mathit{LDA}}\left[\rho \left(\overrightarrow{r}\right)\right]\right)\right)}{\partial \rho \left(\overrightarrow{r}\right)}\\ & +\frac{\partial \left(\rho \left(\overrightarrow{r}\right)\left({\epsilon }_c^{\mathit{LDA}}\left[\rho \left(\overrightarrow{r}\right)\right]\right)\right)}{\partial \rho \left(\overrightarrow{r}\right)}\end{array}$$E41

In the LDA approximation, Eq. (38) must be the departure of our programming calculations. All terms in Eq. (38) are known and remain to define ${\widehat{V}}_{\mathit{ext}}\left(\overrightarrow{r}\right)$ and basis of functions ${\phi }_i\left(\overrightarrow{r}\right)$.

${\widehat{V}}_{\mathit{ext}}\left(\overrightarrow{r}\right)$ and ${\phi }_i\left(\overrightarrow{r}\right)$ are selected as the case we want to treat. We will give further their expressions.

Eq. (37) must be solved self-consistently by starting at a certain density. An effective potential ${\widehat{V}}_{\mathit{eff}}\left(\overrightarrow{r}\right)$ is obtained for which the Schrodinger equation of Kohn-Sham is resolved, and a new electronic density is determined. This process is repeated self-consistently until convergence is reached, and the new electronic density obtained must be very close to the previous one.

4.1 Self-coherent technical to resolve KS Schrödinger equation

The idea is that do not directly resolve Eq. (39), but to write previously the ${\varphi }_m\left(\overrightarrow{r}\right)$ in a finite basis function ${\varphi }_p^b\left(\overrightarrow{r}\right)$ as:

$${\varphi }_m\left(\overrightarrow{r}\right)=\sum _{p=1}^p{C}_p^m{\varphi }_p^b\left(\overrightarrow{r}\right)$$E42

when m = n.$\overrightarrow{k}$, $\overrightarrow{k}$: is the wave vector belonging the first Brillouin zone in the case of crystal lattice.

The resolution of Eq. (38) consists to determine the coefficients $C_p^m$ necessary to express ${\varphi }_m\left(\overrightarrow{r}\right)$ in a given basis ${\varphi }_p^b\left(\overrightarrow{r}\right)$.

We need to search a basis to get as close as possible to ${\varphi }_m\left(\overrightarrow{r}\right)$ with $p$ having a finite value.

Eq. (37) becomes:

$$\left[\begin{array}{ccc}\dots & & \dots \\ \dots & \left.⟨{\varphi }_i^b\right|{\widehat{H}}_{\mathit{KS}}\left|{\varphi }_j^b\right.⟩-{\epsilon }_m⟨{\varphi }_i^b\left|{\varphi }_j^b\right.⟩& \dots \\ \dots & & \dots \end{array}\right]\left[\begin{array}{c}{C}_1^m\\ \dots \\ C_p^m\end{array}\right]=\left[\begin{array}{c}0\\ \dots \\ 0\end{array}\right],$$E43

in which one can identify the matrix Elements of Hamiltonian of single particle and the elements of the recovery matrix, $\left(H_{\mathit{ij}}-{\epsilon }_m{S}_{\mathit{ij}}\right)C_p^m=0$, where ${\widehat{H}}_{\mathit{ij}}=\left.⟨{\varphi }_i^b\right|{\widehat{H}}_{\mathit{KS}}\left|{\varphi }_j^b\right.⟩$ and $S_{\mathit{ij}}=⟨{\varphi }_i^b\left|{\varphi }_j^b\right.⟩$ respectively represent the Hamiltonian matrix and the recovery matrix.

For a solid, these equations need to be resolved for each vector k in the Brillouin zone. This secular equation system is linear with the energy. This system constitutes a problem of determination of the proper values ${\epsilon }_m$ and proper functions ${\varphi }_i^k\left(\overrightarrow{r}\right)$ that much we know within the Hartree theory and is commonly solved from standard numerical methods. The diagonalization of the Hamiltonian matrix provides p proper values clean and p sets values of coefficients that express each p proper functions in a given basis.

More p is big, more the proper functions are precise, but the matrix diagonalization time is also particularly high.

4.2 Self-consistent cycle

Eq. (39) must be resolved in an iterative way in self-consistent cycle procedure. The procedure starts by the definition of a density of departure corresponding to a determined geometry core. Generally, the initial density is constituted from a superposition of atomic densities: ${\rho }^{\mathit{\text{in}}}={\rho }^{\mathit{\text{crystal}}}={\sum }_{\mathit{at}}{\rho }^{\mathit{at}}$.

When the elements of the Hamiltonian matrix and recovery matrix were calculated, proper vectors and proper functions are determined from the diagonalization of the matrix $\left(H_{\mathit{ij}}-{\epsilon }_m{S}_{\mathit{ij}}\right)C_p^m=0$. Following the principle of Aufbau, the orbitals are filled, a new density is determined:

$$\boxed{{\rho }^{\mathit{\text{out}}}\left(\overrightarrow{r}\right)=\sum _{\mathit{\text{out}}}{\left|{\varphi }_i\left(\overrightarrow{r}\right)\right|}^2},$$

This step concludes the first cycle. At this stage of the process, acceleration of convergence is generally used to generate a new density realized from a mixture between this output density ${\rho }^{\mathit{\text{out}}}\left(r\right)$ and density of entry of this cycle and there ${\rho }^{\mathit{\text{in}}}\left(r\right)$ . One of the simplest procedures concerning this mixture can be formulated as: ${\rho }_{i+1}^{\mathit{\text{in}}}\left(r\right)=\left(1-\alpha \right){\rho }_i^{\mathit{\text{in}}}\left(r\right)+{\mathit{\alpha \rho }}_i^{\mathit{\text{out}}}\left(r\right)$ with $0\prec \alpha \prec 1$. $\alpha$ is the mixture parameter an $i$ corresponds to iterative cycle number. Density of the new entry ${\rho }_{i+1}^{\mathit{\text{in}}}\left(r\right)$ is then introduced into a second self-consistent cycle. This process is repeated for an iterative manner until it has the convergence criterion (difference between ${\rho }^{\mathit{\text{out}}}\left(r\right)$ and ${\rho }^{\mathit{\text{in}}}\left(r\right)$) initially fixed is reached. When convergence is reached, the energy of the ground state of the system is reached (Figure 1).

Figure 1.

Schematic representation of the self-consistent cycle.

Once we obtain ${\rho }^{\mathit{\text{out}}}\left(\overrightarrow{r}\right)$ for the first cycle, we inject it like a new value of ${\rho }^{\mathit{\text{in}}}\left(\overrightarrow{r}\right)$ and we repeat this iterative operation until it is that the fixed precision at the beginning is reached.

Precision = ${\rho }^{\mathit{\text{out}}}\left(\overrightarrow{r}\right)-{\rho }^{\mathit{\text{in}}}\left(\overrightarrow{r}\right)$

$${\rho }^{\mathit{\text{in}}}\left(\overrightarrow{r}\right)={\rho }^{\mathit{\text{orbital}}}\left(\overrightarrow{r}\right)=\sum _{\mathit{at}}{\rho }^{\mathit{at}}\left(\overrightarrow{r}\right),$$

We can use ${\varphi }_k^n\left(\overrightarrow{r}\right)=e^{i\overrightarrow{k}.\overrightarrow{r}}{u}_k^n\left(\overrightarrow{r}\right)$, where $e^{i\overrightarrow{k}.\overrightarrow{r}}$ is a plane wave and $u_k^n\left(\overrightarrow{r}\right)$ a function processing Lattice periodicity (Bloch theorem $u_k^n\left(\overrightarrow{r}+\overrightarrow{k}\right)=u_k^n\left(\overrightarrow{r}\right)$).

5. Time dependent density functional theory (TDDFT) calculations

The Eigen values of Kohn-Sham does not match corresponding to the required energy for excited electrons in a correlated electron system with N electrons, because the Kohn-Sham [10] approach replaces the interacting electron problem by an independent electron problem. There is a way to study the potential and the wave function to determine the excitation energies. It is the use of the TDDFT theory which takes into account the time parameter. The excitation energy of a system can then be obtained from the response function to extreme electron density perturbations.

5.1 Rung-Gross theorems

TDDFT is based on the Rung-Gross theorem [23, 24] that is the analog of the time-dependent of Hohenberg-Kohn theorems [9]. The theorems are cited below:

5.1.1 First theorem

The time depends on the external potential ${\widehat{V}}_{\mathit{ext}}\left(\overrightarrow{r},t\right)$ is determined by the time-dependent electron density $\rho \left(\overrightarrow{r},t\right)$ to a nearly additive function for an initial state $\psi$ a state$\psi \left(t=0\right)$.

5.1.2 Second theorem

By difference with the ground state, the variational principle, which stated that there is a minimum associated with the total energy does not exist for the time-dependent systems because the energy is not a conserved quantity. A similar amount of energy, which is applied on the stationary principle, is defined as:

$$\mathrm{{\rm A}}\left[\psi \right]={\int }_{t_0}^{t_1}\mathit{dt}\left.⟨\psi \left(t\right)\right|i\frac{\partial }{\partial t}-\widehat{H}\left(t\right)\left|\psi \left(t\right)\right.⟩$$E44

where $\mathrm{{\rm A}}$ is called the action, and $\psi \left(t\right)$ is a function of the time-dependent poly electronic wave function.

5.2 The time-dependent approach of Kohn-Sham

The time-dependent approach [10] of Kohn-Sham shows an equation with partial derivatives for the effective system such as

$$i\frac{\partial }{\partial t}{\phi }_i\left(\overrightarrow{r},t\right)=H{\phi }_i\left(\overrightarrow{r},t\right)=\left[-\frac{1}{2}{\nabla }^2+{\widehat{V}}_{\mathit{eff}}\left(\rho \left(\overrightarrow{r},t\right)\right)\right]{\phi }_i\left(\overrightarrow{r},t\right)$$E45

$$\rho \left(\overrightarrow{r},t\right)=\sum _{i=1}^N{\left|{\phi }_i\left(\overrightarrow{r},t\right)\right|}^2$$E46

With ${\phi }_i\left(\overrightarrow{r},t\right)$ is a single electron time-dependent wave function of the Kohn-Sham. ${\widehat{V}}_{\mathit{eff}}\left(\rho \left(\overrightarrow{r},t\right)\right)$ defined in Eq. (46) contains in addition to the external field, a Hartree potential, and the exchange-correlation potential such as:

$${\widehat{V}}_{\mathit{eff}}\left(\rho \left(\overrightarrow{r},t\right)\right)={\widehat{V}}_{\mathit{ext}}\left(\overrightarrow{r},t\right)+\frac{\partial E_H\left(\overrightarrow{r},t\right)}{\partial \rho \left(\overrightarrow{r},t\right)}+{\widehat{V}}_{\mathit{xc}}\left(\rho \left(\overrightarrow{r},t\right)\right)$$E47

The exchange-correlation term ${\widehat{V}}_{\mathit{xc}}\left(\rho \left(\overrightarrow{r},t\right)\right)$ is unknown. There is a major difficulty to evaluate it because it does not depend only on the density $\rho \left(\overrightarrow{r},t\right)$, but also depends on the prior density at time $t$ in every point in space. The evolution of ${\widehat{V}}_{\mathit{xc}}\left(\rho \left(\overrightarrow{r},t\right)\right)$ therefore requires the introduction of the single fundamental approximation of TDDFT.

6. Linear response theory

The TDDFT approach allows one hand to study the perturbation of the system at time $t_0$, secondly to propagate this disturbance for a time $t>t_0$. The study of the evolution in the propagation of this disturbance leads to the production of the absorption spectrum. This theory [25, 26] is used in the space of frequencies ω rather than in the temporal space. The operation is performed by means of the Fourier transform.

Considering, the continuous function $f$ of the time variable. The Fourier transform of $f$ gives:

$$\mathit{TF}\left(f\right):\omega \mapsto f\left(\omega \right)=\frac{1}{\sqrt{2\pi }}{\int }_{-\infty }^{+\infty }f\left(t\right)e^{-\mathit{i\omega t}}\mathit{dt}$$E48

Considering a system in the ground state of electron density ${\rho }_0$ which is subject at $t=t_0$ to a low disturbance of the external potential $\partial {\widehat{V}}_{\mathit{ext}}$ the electron density obtained following this disturbance is, at first order given by:

$$\rho \left(\overrightarrow{r},\omega \right)={\rho }_0\left(\overrightarrow{r}\right)+\mathit{\delta \rho }\left(\overrightarrow{r},\omega \right)$$E49

In the interacting electrons system, the external potential is also written as the sum of the external potential calculated in the ground state, supplemented with a disruptive potential:

$${\widehat{V}}_{\mathit{ext}}\left(\overrightarrow{r},\omega \right)={\widehat{V}}_{\mathit{ext}}\left(\overrightarrow{r}\right)+\delta {\widehat{V}}_{\mathit{ext}}\left(\overrightarrow{r},\omega \right)$$E50

The perturbation $\mathit{\delta \rho }$ depends only on the potential $\delta V_{\mathit{ext}}$ and takes the expression:

$$\mathit{\delta \rho }\left(\overrightarrow{r},\omega \right)={\int }_{r\text{'}}\chi \left(\overrightarrow{r},{\overrightarrow{r}}^{\text{'}},\omega \right)\delta {\widehat{V}}_{\mathit{ext}}\left({\overrightarrow{r}}^{\text{'}},\omega \right)d{\overrightarrow{r}}^{\text{'}}$$E51

$\chi \left(\overrightarrow{r},{\overrightarrow{r}}^{\text{'}},\omega \right)$: represents the response function of the non-interacting electrons system. The evolution of $\chi$ is complex, but it can be simplified by using the Kohn-Sham approach that redefines $\mathit{\delta \rho }\left(\overrightarrow{r},\omega \right)$ according to the equation:

$$\mathit{\delta \rho }\left(\overrightarrow{r},\omega \right)={\int }_{r\text{'}}{\chi }_{\mathit{KS}}\left(\overrightarrow{r},{\overrightarrow{r}}^{\text{'}},\omega \right)\delta {\widehat{V}}_{\mathit{eff}}\left({\overrightarrow{r}}^{\text{'}},\omega \right)d{\overrightarrow{r}}^{\text{'}}$$E52

where ${\chi }_{\mathit{KS}}\left(\overrightarrow{r},{\overrightarrow{r}}^{\text{'}},\omega \right)$, according to Kohn-Sham response is easily calculated using the equation:

$${\chi }_{\mathit{KS}}\left(\overrightarrow{r},{\overrightarrow{r}}^{\text{'}},\omega \right)=\underset{\eta \to 0^{+}}{lim}\sum _{p=1}^{\infty }\sum _{q=1}^{\infty }\left(f_p-f_q\right)\frac{{\phi }_p^{\ast }\left(\overrightarrow{r}\right){\phi }_q\left(\overrightarrow{r}\right){\phi }_q^{\ast }\left({\overrightarrow{r}}^{\text{'}}\right){\phi }_p\left({\overrightarrow{r}}^{\text{'}}\right)}{\omega -{\omega }_{p\to q}+\mathit{i\eta }}$$E53

This expression (52) represents the various excitations between the occupied Kohn-Sham orbital’s ${\phi }_p$ and unoccupied ${\phi }_q$, for a total number N of orbital’s approaching infinity. The occupation of the orbital’s number p and q are respectively denoted $f_p$ and $f_q$. The frequency ${\omega }_{p\to q}={\epsilon }_p-{\epsilon }_q$; where ${\epsilon }_p$ and ${\epsilon }_q$ are the eigenvalues ​​respectively associated with wave functions ${\phi }_p$ and ${\phi }_q$. The infinitesimal positive number $\eta$ is introduced to reflect the causality of the system response.

Disruption of the effective potential $\delta {\widehat{V}}_{\mathit{eff}}$ (Kohn-Sham potential) can be written as the sum of three terms:

• the Coulomb term: ${\int }_{r\text{'}}\frac{\mathit{\delta \rho }\left({\overrightarrow{r}}^{\text{'}},\omega \right)}{\left|\overrightarrow{r}-{\overrightarrow{r}}^{\text{'}}\right|}d{\overrightarrow{r}}^{\text{'}}$

• the external term: $\delta {\widehat{V}}_{\mathit{ext}}\left(\overrightarrow{r},\omega \right)$

• the exchange-correlation term: $\delta {\widehat{V}}_{\mathit{xc}}\left(\overrightarrow{r},\omega \right)$

Such:

$$\delta {\widehat{V}}_{\mathit{eff}}\left(\overrightarrow{r},\omega \right)=\delta {\widehat{V}}_{\mathit{ext}}\left(\overrightarrow{r},\omega \right)+{\int }_{r\text{'}}\frac{\mathit{\delta \rho }\left({\overrightarrow{r}}^{\text{'}},\omega \right)}{\left|\overrightarrow{r}-{\overrightarrow{r}}^{\text{'}}\right|}d{\overrightarrow{r}}^{\text{'}}+\delta {\widehat{V}}_{\mathit{xc}}\left(\overrightarrow{r},\omega \right)$$E54

Following a Fourier transformation, it is customary to write disturbance of the exchange-correlation potential as:

$$\delta {\widehat{V}}_{\mathit{xc}}\left(\overrightarrow{r},\omega \right)=\mathit{TF}\left({\int }_{r\text{'}}{\int }_{t\text{'}}{f}_{\mathit{xc}}\left(r,t,r^{\text{'}},t^{\text{'}}\right)\mathit{\delta \rho }\left(r^{\text{'}},t^{\text{'}}\right){\mathit{dt}}^{\text{'}}{\mathit{dr}}^{\text{'}}\right)$$E55

Where

$$f_{\mathit{xc}}\left(r,t,r^{\text{'}},t^{\text{'}}\right)=\frac{\delta {\widehat{V}}_{\mathit{xc}}\left(\overrightarrow{r},\omega \right)}{\mathit{\delta \rho }\left(r^{\text{'}},t^{\text{'}}\right)},$$E56

$f_{\mathit{xc}}\left(r,t,r^{\text{'}},t^{\text{'}}\right)$ is the expression of exchange and correlation core.

In this step of the calculation, the development of the disturbance of the electron density (Eq. 51) through the expression of the effective potential disturbance (Eq. 52) gives:

$$\mathit{\delta \rho }\left(\overrightarrow{r},\omega \right)={\int }_{r\text{'}}{\chi }_{\mathit{KS}}\left(\overrightarrow{r},{\overrightarrow{r}}^{\text{'}},\omega \right)\delta {\widehat{V}}_{\mathit{ext}}\left({\overrightarrow{r}}^{\text{'}},\omega \right)d{\overrightarrow{r}}^{\text{'}}+{\int }_{r_1}{\int }_{r_2}{\chi }_{\mathit{KS}}\left(\overrightarrow{r},{\overrightarrow{r}}_1,\omega \right)f_{\mathit{Hxc}}\left({\overrightarrow{r}}_1,{\overrightarrow{r}}_2,\omega \right)\mathit{\delta \rho }\left({\overrightarrow{r}}_2,\omega \right){\mathit{dr}}_2{\mathit{dr}}_1$$E57

Where

$$f_{\mathit{Hxc}}\left(\overrightarrow{r},{\overrightarrow{r}}^{\text{'}},\omega \right)=\frac{1}{\left|\overrightarrow{r}-{\overrightarrow{r}}^{\text{'}}\right|}+f_{\mathit{xc}}\left(\overrightarrow{r},{\overrightarrow{r}}^{\text{'}},\omega \right)$$E58

$f_{\mathit{Hxc}}\left(\overrightarrow{r},{\overrightarrow{r}}^{\text{'}},\omega \right)$ is the expression of the Hartree exchange-correlation core.

Expression of the disturbance of the electron density is known both in the interacting electron system (Eq. 57) and non-interacting electrons system (Eq. 51). The equality of these two expressions leads to that of the response function of the interacting electron system:

$$\chi \left(\overrightarrow{r},{\overrightarrow{r}}^{\text{'}},\omega \right)={\chi }_{\mathit{KS}}\left(\overrightarrow{r},{\overrightarrow{r}}^{\text{'}},\omega \right)+{\int }_{r_1}{\int }_{r_2}\chi \left(\overrightarrow{r},{\overrightarrow{r}}_1,\omega \right)f_{\mathit{Hxc}}\left({\overrightarrow{r}}_1,{\overrightarrow{r}}_2,\omega \right){\chi }_{\mathit{KS}}\left({\overrightarrow{r}}_2,{\overrightarrow{r}}^{\text{'}},\omega \right){\mathit{dr}}_2{\mathit{dr}}_1$$E59

Equation of the susceptibility $\chi$ is important because this linear response has poles exactly to the energies of the electronic transitions of the system.

• When $f_{\mathit{Hxc}}\left({\overrightarrow{r}}_1,{\overrightarrow{r}}_2,\omega \right)=0$, the electronic transitions are exactly equal to those given by the Kohn-Sham system in the ground state. Their oscillator strengths are then given by the poles of ${\chi }_{\mathit{KS}}$.

• When $f_{\mathit{Hxc}}\left({\overrightarrow{r}}_1,{\overrightarrow{r}}_2,\omega \right)\ne 0$, the excitation energies are corrected by the exchange-correlation Hartree core, and the oscillator strengths of the transitions are given by the poles $\chi$.

7. Resolution of time-dependent problem

7.1 Adiabatic approximation

No analytical expression exists for the time-dependent exchange correlation potential. So, we work using the adiabatic approximation [26] by ignoring the dependence on previous electron densities of these two terms and use only the instantaneous electron density. Furthermore, when the time-dependent potential changes slowly, thus adiabatically, the system remains in its instantaneous ground state. The exchange-correlation potential developed in the DFT to describe the ground state of the system can then be translated to time-dependent systems. The potential and the adiabatic exchange-correlation core then take local forms in time given to equations (Eq. 60) and (Eq. 59). It is possible to note that in this approximation, the exchange-correlation core is independent of the frequency of the disturbance potential applied to the system.

$${\widehat{V}}_{\mathit{xc}}^{\mathit{\text{adia}}}\left(\overrightarrow{r},t\right)={\widehat{V}}_{\mathit{xc}}^{\mathit{DFT}}{\left[\rho \left(\overrightarrow{r}\right)\right]}_{\rho \left(\overrightarrow{r}\right)=\rho \left(\overrightarrow{r},t\right)}$$E60

$$f_{\mathit{xc}}^{\mathit{\text{adia}}}\left(\overrightarrow{r},t,r^{\text{'}},t^{\text{'}}\right)=\frac{\delta {\widehat{V}}_{\mathit{xc}}\left(\overrightarrow{r}\right)}{\mathit{\delta \rho }\left({\overrightarrow{r}}^{\text{'}}\right)}\delta \left(t-t^{\text{'}}\right)$$E61

The simplicity of this approximation allows to make after Fourier transform, the exchange-correlation core $f_{\mathit{xc}}^{\mathit{\text{adia}}}$ independent of the frequency $\omega$. Therefore, $f_{\mathit{Hxc}}$ also becomes independent of the frequency as:

$${\widehat{f}}_{\mathit{Hxc}}^{\mathit{\text{adia}}}\left(\overrightarrow{r},{\overrightarrow{r}}^{\text{'}},\omega \right)=f_{\mathit{Hxc}}^{\mathit{\text{adia}}}\left(\overrightarrow{r},{\overrightarrow{r}}^{\text{'}}\right)=\frac{1}{\left|\overrightarrow{r}-{\overrightarrow{r}}^{\text{'}}\right|}+f_{\mathit{xc}}^{\mathit{\text{adia}}}\left(\overrightarrow{r},{\overrightarrow{r}}^{\text{'}}\right)$$E62

In the adiabatic approximation, using functional developed for DFT retains weaknesses of these functional. For example, functional LDA and most functional GGA falsely represent the decay of exchange-correlation potential in neutral finite systems.

7.2 Equation with Eigenvalues

When the exchange-correlation core was independent of the frequency (in the adiabatic approximation case), then search for the poles of the response factor by solving a system of equations to the eigenvalues [22, 26] ​​of a matrix equation as the specific form:

$$\sum _{p\text{'}=1}^{\infty }\sum _{q\text{'}=1}^{\infty }{M}_{\mathit{pq},p\text{'}q\text{'}}\left[f_{\mathit{Hxc}}^{\mathit{\text{adia}}}\right]X_{p\text{'}q\text{'}}={\Omega }_{\mathit{pq}}^2{X}_{\mathit{pq}}$$E63

where $p$, $p^{\text{'}}$ et $q$, $q^{\text{'}}$ respectively represent indices of non-occupied and occupied orbitals for the orbital’s base with size tends to infinity. $M_{\mathit{pq},p\text{'}q\text{'}}$ is an element of the operator matrix writing such that:

$$M_{\mathit{pq},p\text{'}q\text{'}}\left[f_{\mathit{Hxc}}^{\mathit{\text{adia}}}\right]=2W_{\mathit{pq},p\text{'}q\text{'}}\left[f_{\mathit{Hxc}}^{\mathit{\text{adia}}}\right]\sqrt{{\omega }_{p\to q}{\omega }_{p\text{'}\to q\text{'}}}+{\omega }_{p\to q}^2\delta \left(\mathit{pq}\right)\left(p^{\text{'}}{q}^{\text{'}}\right),$$E64

The matrix element ${\omega }_{\mathit{pq},p\text{'}q\text{'}}$ is function as the Hartree exchange-correlation core, itself-independent of the frequency as:

$${\omega }_{\mathit{pq},p\text{'}q\text{'}}\left[f_{\mathit{Hxc}}^{\mathit{\text{abia}}}\right]={\int }_r{\int }_{r\text{'}}{\phi }_p^{\ast }\left(\overrightarrow{r}\right){\phi }_{p\text{'}}^{\ast }\left({\overrightarrow{r}}^{\text{'}}\right)f_{\mathit{Hxc}}^{\mathit{\text{adia}}}\left(\overrightarrow{r},{\overrightarrow{r}}^{\text{'}}\right){\phi }_q\left(\overrightarrow{r}\right){\phi }_{q\text{'}}\left({\overrightarrow{r}}^{\text{'}}\right)d\overrightarrow{r}d{\overrightarrow{r}}^{\text{'}}$$E65

Where $X_{\mathit{pq}}$ et ${\Omega }_{\mathit{pq}}$ are respectively the vectors and eigenvalues.

8. Conclusion

Our main objective is the determination of the electronic structure, so we have proposed the model of the Functional Density Theory dependent and independent of time which is widely used; it has the advantage of taking into account electron correlation directly in its formalism. To test the relevance of these theoretical calculations and the limit of validity of all these approximations used during the development of this theoretical model, a computer simulation model is necessary. Indeed, it will allow us to obtain the eigenvalues and eigenvectors. They provide enough information to access to the electronic structure and total energy. This theoretical investigation may allow to work toward the realistic modeling of the electronic structure of spherical and cylindrical nanoparticles, to provide effective potentials and orbitals that we can employ to calculate photoemission spectra, and may allow to improve the modeling of the ground-state electronic structure of metal surfaces within the framework of the streaked photoemission calculations, and to develop a theoretical framework for dielectric response of surfaces and nanostructures to external IR fields in order to model plasmonic electric fields [27, 28, 29, 30] enhancement near nanostructures systems.

Acknowledgments

We are grateful to Professor Uwe Thumm who hosted us for 3 months in his James R. Macdonald laboratory at the Kansas State University in the USA, and who offered us an opportunity to collaborate on this subject, as part of the Fulbright Grant Merit Award.