### The Bethe-Salpeter equation solvers

In this section we will describe two different methods to derive the macroscopic dielectric function including excitonic effects:

"Diagonalization",

"Haydock recursive algorithm". Before presenting such methods we should define the macroscopic dielectric function which, according to

Adler and

Wiser, can be directly related to the inverse of the microscopic dielectric matrix in the following way

It is important to note that the optical absorption spectrum, which is the quantity we can compare to the experimental data is given by the imaginary part of the last expression.

In order evaluate the macroscopic dielectric function using the "Diagonalization algorithm", we have to write the expression defining the macroscopic dielectric function in a different way. In particular, as it has been shown in

Onida and

Albrecth, it can be rewritten as follows

where, without entering in the details of the calculations

and

Further, in order to describe the problem in terms of electron-hole pairs, it is useful to introduce the four-point polarization

which corresponds to the old polarization if x1=x2 and x3=x4. Thus the macroscopic dielectric function will be given by

where we have introduced the Fourier transform of the four-point polarization. Teh first step we are going to perform, in order to evaluate this expression, is the projection of the four-point polarization in energy space using a basis of LDA Bloch functions. The idea to describe the problem in single particles energy space is suggested by the fact that we expect that only a "little" number of excitations will contribute to each excitation, so the computational effort to evaluate the excitonic effects will be reduced. This practically means that we will perform the following transformation

From now on we will omit the frequency dependence if not necessary in the understanding of the formulation. As it has been shown in

BSK the four-point polarization satisfies a Dyson like equation: the Bethe-Salpeter equation. Usually this equation is written for L, which is defined as follows

In practice, in the basis of the single particle Bloch wavefunctions, we have to solve

in this expression L0 represents the independent particle polarization, which has the following expression

while

*K* is the kernel describing the particle-hole interaction effects. In the rest of this section we will try to describe these effects and show how they crucially modify the optical spectra. Now, in order to solve the Bethe-Salpeter equation, we define

so we will have for L

It can be shown with some algebraic calculation, which can be found in

Albrecth, that P can be written as follows

where we have defined the excitonic Hamiltonian

This two-particle hamiltonian has two contributions: the first, containing the energy difference between a conduction and a valence state, is analogous to the free particle energy in the single particle hamiltonian, while the second, which is equivalent to the self-energy, describes the electron-hole interaction. This means that we have reduced the problem to the inversion of the excitonic hamiltonian. Indeed

where

Thus we have to invert a matrix for each absorption frequency. This is computationally very expensive, so we need some kind of "trick". In practice we use a spectral representation for the inversion, using the relation

which holds for a system of eigenvectors and eigenvalues of a general non-Hermitian matrix defined by

and where S is the overlap matrix (of the in gneral non orthogonal eigenstates) of the excitonic hamiltonian.

More explicitly we have to solve the

*effective two particle Schrodinger equation* for the excitonic system by:

The eigenvalues, which represent the excitonic states energies, and the eigenvector, which are the excitonic wavefunctions, will be determined by the

**diagonalization** of the excitonic hamiltonian. Once we have determined the eigenvalues and eigenvectors of the excitonic hamiltonian, we get for L

The presence of the difference between two Fermi functions in the last equation means that we can restrict to the subspace in which the couples (n1, n2) refere to a conduction and valence state (c,v) or (v,c). Within this subspace the two-particle excitonic hamiltonian will have the form (immaginate che sia un array)

where the

*resonant part* is Hermitian

while the

*coupling part* is symmetric

We also denote the part in the lower right part of the array describing the excitonic hamiltonian as

*anti-resonant*. It contains negative transition energies, thus only contributing to negative frequencies.
Note that when running the code you have the possibilty to choose to take into account the resonant part or the coupling part of the excitonic hamiltonian:

BSKmod. In the code one has even the possibility to take into account the exchange part of the excitonic hamiltonian, the meaning of this term will be soon clarified when the analytical expression of the single terms in the excitonic Hamiltonian will be presented. Indeed concerning the

**resonant** part we have

the diagonal part is trivial, as it contains the difference between single particle valence and conduction states

while the unscreened short-range exchange term, can be easily evaluated and results to be

where V represents the volume of the crystal. Finally the Coulomb term can be written as follows

where we have defined the

*symmetrized inverse dielectric* function

these expressions completely descibe the

*resonant* contribution to the excitonic hamiltonian. It should be noted that some terms in the sums diverge if G, G prime and eventually q are equal to zero, this means that some tricks have to be used in order to cure these divergencies, altough it is beyond the scope of this introductory description of the Bethe-Salpeter to present these details, the interested reader can refer to

Albrecth. Now we have to describe the

**Coupling** part of the excitonic Hamiltonian, this term can be split as follows

the unscreened short-range exchange term is

the Coulomb term gives

Once we have diagonalized the excitonic Hamiltonian we can evaluate the macroscopic dielectric function. Indeed placing the expression of L (which is simply iX) in the equation which defines the macroscopic dielectric function, we obtain

Note that we have added a small imaginary constant to the frequency omega in order to smooth the curve by a Lorentzian shape and to shift the poles away from the real axis. If we consider only the resonant part of the excitonic Hamiltonian we obtain a simpler formula for the imaginary part of the dielectric function

The method just presented requires to perform the diagonalization of the excitonic hamiltonian in order to determine the macroscopic dielectric function. Such procedure, requiring the storage in memory of at least half of the excitonic hamiltonian, can not be applyed to the case of very large systems such as surfaces where the number of k-points and transitions needed can become huge. Thus a new procedure has to be introduced, indeed the Haydock recursive method allows to determine the dielectric function within an iterative scheme, where, in each step, a matrix-vector product has to be performed. This means that we don't need to store in memory the full excitonic hamiltonian but we need to know only a matrix line each time. A complete and detailed analysis of this method can be found in

Haydock,

Benedict and

Marsili. As a first step to apply this algorithm, we should rewrite the dielectric function in a different way, in particular it should have the following structure

where

*O* reprents a generic observable, while H is the hamiltonian of the system. This means that we have to derive a new expression of the dielectric function. To derive it we suppose to apply an electromagnetic field to the system. In the

*Coulomb gauge*, neglecting non linear effects, we will find that the interaction of such field with the electrons is described by the following Hamiltonian

where

**A** is the vector potential and

**v** is the velocity operator for each electron. If we further suppose that the elctromagnetic radiation is monocromatic, with frequency omega, and has an amplitude A0, from the Fermi golden rule, we will find that the transition rate from the initial state |i> to a generic final state |f> will have this form

where

*e* is the polarization vector. Following the derivation presented in

Albrecth we can connect the transition rates to the imaginary part of the dielectric function through the following expression

where |0> is the many body ground state. It is also necessary to write the velocity operator in a more convenient way (

Del Sole) as

so the imaginary part of the dielectric function will result as follows

This equation can be written in a different way, closer to the result we want to obtain, considering the following relation:

where

*P* represents the principal part of the integral. The dielectric function will be given by

Finally this is the expression we were looking for, indeed if we set

we recover the definition of

*I* given at the beginning of this section. Now we are going to show how it can be evaluated in terms of a continued fraction by an

**iterative procedure**. As a first step in our iterative scheme we set

*P*, which is defined by

equal to |1> and we build an orthonormal basis using the following scheme

it is straightforward to see that <1|2>=0, and that <1|H|2>=b1. We set a2=<2|H|2> and define |3> as

where we have implicitely defined b2. Starting from these three states it is possible to give an iterative definition of the rest of the basis:

It is easy to show this is an orthonormal set and that, in this basis, the Hamiltonian matrix has a tridiagonal form in terms of the coefficients (ai,bi). It is easy to show that

*I* can now be written as a continued fraction

This fraction can be terminated in Yambo by switching on the BSHayTer flag.
flag.

To summarize, in order to obtain the spectra, one has to perform the following steps:

- define the starting state |P>=|1>
- iteratively apply the Hamiltonian to the state |P>. In this way the orthonormal set of states {|i>} will be generated and the coefficients {ai,bi} can be stored
- build the continued fraction, and compute the imaginary part
- check the convergence of the spectra with the number of iterations used.

**NEW**
The Haydock procedure as described hereabove holds only for Hermitian
Hamiltonians and is therefore used in combination with the Tamm-Dancoff
approximation in which the coupling part of the Hamiltonian, rending the
matrix non-Hermitian, is neglected.
When the coupling is considered, Yambo uses instead the approach recently
proposed by Gruning et al that generalized the Haydock procedure to the full BSE Hamiltonian.

#### References

- S. L. Adler, Phys. Rev.
**126**, 413 (1962)
- N. Wiser, Phys. Rev.
**129**, 72 (1963)
- G. Onida, L. Reining, and A. Rubio, Rev. Mod. Phys.
**48**, 601 (2002)
- S. Albrecht, Ph.D. thesis, Ecole Polytechnique, Palaiseau (Paris), 1999
- R. Haydock, Comput. Phys. Comm.
**20**, 11 (1980)
- L. X. Benedict and E. L. Shirley, Phys. Rev. B
**59**, 5441 (1999)
- M. Marsili, Ph.D. thesis, Universita di Roma "Tor Vergata", 2006
- R. Del Sole and R. Ghirlanda, Phys. Rev. B
**48**, 11789 (1993)
- Myrta Gruning, Andrea Marini and
Xavier Gonze, Nano Letters
**9** (2009), pp 2820-2824.
Link to article