Efficient O„ N
2
… method to solve the Bethe-Salpeter equation
W. G. Schmidt,* S. Glutsch, P. H. Hahn, and F. Bechstedt
Institut fu ¨r Festko ¨rpertheorie und Theoretische Optik, Friedrich-Schiller-Universita ¨t Jena, Max-Wien-Platz 1, 07743 Jena, Germany
Received 3 October 2002; published 6 February 2003
We present a numerically efficient approach to solve the Bethe-Salpeter equation for the polarization func-
tion. Rather than from the usual eigenvalue representation, the macroscopic polarizability is obtained from the
solution of an initial-value problem. This reduces the computational effort considerably and allows for calcu-
lating excitonic and local-field effects in optical spectra of complex systems consisting of many atoms. As an
example we investigate the optical anisotropy of the monohydride Si(001)(2 1) surface. While excitonic
effects influence the surface optical properties considerably, the local-field effect induced changes are minimal.
DOI: 10.1103/PhysRevB.67.085307 PACS numbers: 73.20.At, 78.68.+m, 71.15.Qe
I. INTRODUCTION
Recent years have seen impressive methodological
progress in the accurate numerical modeling of optical prop-
erties from first principles see, e.g., Ref. 1. It has become
possible not only to calculate single-particle electronic exci-
tation energies accurately using the GW approximation
GWA, but also to solve the Bethe-Salpeter equation BSE
for pair excitations in order to account for excitonic and
local-field LF contributions to the optical response.
2–4
However, the large numerical effort required to solve the
BSE has restricted such calculations to the interaction of
relatively few electron-hole pairs. Therefore its application
has been limited to bulk semiconductors,
2,3,5
strongly local-
ized surface states,
6,7
small clusters,
8
or molecules.
9
At the same time, methods of optical spectroscopy are
rapidly gaining importance for materials characterization.
Techniques such as reflectance anisotropy spectroscopy
RAS have evolved from experimental methods to charac-
terize static surfaces to very powerful in situ diagnostic
probes which allow for the monitoring and controlling of
surface growth in real time and in challenging environments
such as in high pressures or under liquids.
10
In order to fully
exploit the potential of such methods, however, the accurate
theoretical modeling of the optical properties of large and
complex systems—such as surfaces—is required. A number
of technical improvements such as optimized schemes to cal-
culate the electron-hole interaction in reciprocal space
11
and
methodological developments which allow to obtain the po-
larization function from iterative schemes
12
have been sug-
gested in order to extend the applicability of the BSE to
larger and potentially more interesting systems.
In the present work we suggest an alternative approach to
solve the BSE. It is characterized by a O( N
2
) scaling of the
operation count with N being the number of electron-hole
pair states and allows for the accurate modeling of excitonic
and LF effects in systems consisting of comparatively many
atoms. After a brief description of the proposed methodology
and its test for bulk Si we demonstrate its applicability to
large systems by calculating the optical anisotropy of the
monohydride Si(001)(2 1) surface in a wide spectral
range.
II. METHODOLOGY
We start from first-principles pseudopotential calcula-
tions, using a massively parallel real-space finite-difference
implementation of the density-functional theory in local-
density approximation DFT-LDA.
13
A multigrid technique
is used for convergence acceleration. In order to include
electronic self-energy effects one needs to replace the local
exchange and correlation potential V
XC
( r) in the LDA by the
nonlocal and energy-dependent self-energy operator
( r, r'; E ) see, e.g., Refs. 14 and 15. For the calculation of
we use the GW approximation,
16,17
where the self-energy
operator is expressed as the convolution =iGW of the dy-
namically screened Coulomb potential W and the single-
particle propagator G. Since the calculation of surface optical
spectra involves a very large number of electronic states,
however, we introduce further approximations following the
schemes developed by Hybertsen and Louie
18
and Bechstedt
et al.:
19
the GW quasiparticle energies are obtained from the
DFT-LDA eigenvalues in a perturbative manner by
n
k
QP
=
n
k +
1
1 +
n , k
n , k
st
+
n , k
dyn
n
k -V
n , k
XC
,
1
where the self-energy operator has been divided into static
( st ) and dynamic ( dyn ) contributions. Indices at and V
XC
indicate diagonal matrix elements with the respective wave
functions.
n , k
is the linear coefficient in the expansion of
dyn
around the DFT-LDA eigenvalue
n
( k). The static part
can be further divided into two parts,
st
r,r' =
1
2
n , k
n , k
r
n , k
* r' W r,r';0 -v rÀr'
-
v , k
v , k
r
v , k
* r' W r,r';0 , 2
representing the Coulomb hole
COH
and the screened ex-
change
SEX
. The
n , k
are the DFT-LDA wave functions.
SEX
contains a sum over the occupied valence states v only.
The major bottleneck in the GW calculation is the computa-
tion of the screened interaction W. An extreme acceleration
can be achieved by using a model dielectric function, for
which several functional forms have been suggested. We use
the version suggested by Bechstedt et al.
19
PHYSICAL REVIEW B 67, 085307 2003
0163-1829/2003/678/0853077/$20.00 ©2003 The American Physical Society 67 085307-1