PHYSICAL REVIEW B 85, 195127 (2012)

Analytical optimization of spread and change of exponential decay of generalized Wannier functions
in one dimension

Denis R. Nacbar1 and Alexys Bruno-Alfonso2

1Programa de Pós-Graduação em Ciência e Tecnologia de Materiais, Unesp, 17033-360 Bauru, São Paulo, Brazil
2Departamento de Matemática, Faculdade de Ciências, Unesp, 17033-360 Bauru, São Paulo, Brazil
(Received 7 November 2011; revised manuscript received 22 March 2012; published 14 May 2012)

Generalized Wannier functions of a couple of bands in a one-dimensional crystal are investigated. A lower
bound for the global minimum of the total spread is obtained. Assumption of such a value being the minimum
leads to a first-order differential equation for the transformation matrix. Simple analytical solutions leading to real
generalized Wannier functions are presented for consecutive bands in a crystal with inversion symmetry. Results
are displayed for a particle in a diatomic Kronig-Penney potential. For the lowest couple of bands, calculated
single-band Wannier functions resemble orbitals of a diatomic molecule, whereas generalized Wannier functions
seem like orthogonalized atomic orbitals. The latter functions are neither symmetric nor antisymmetric and
display increased exponential decay. For the next two pairs of bands, Wannier functions retain their centers and
symmetries, and exponential decay does not increase. Results are also shown to be in agreement with solutions
of the eigenvalue problem of the band-projected position operator.

DOI: 10.1103/PhysRevB.85.195127 PACS number(s): 71.20.−b, 71.15.Dx, 31.15.aq

I. INTRODUCTION

During the centennial year of Gregory Hugh Wannier,
nearly a hundred publications have demonstrated that func-
tions named after him are key computational tools in several
research areas, such as insulators,1 magnetism,2,3 metals,4

photonics,5 polymers,6 semiconductors,7,8 superconductivity,9

and electronic transport.10 Of course, applications have
fueled basic research where Wannier functions (WFs) are
the subject matter,11 and textbooks devote some attention to
these functions.12–16 Since their introduction in 1937,17 the
properties of WFs have challenged the solid-state physics
community.18–20 This is mainly because of their lack of
uniqueness, which is a feature inherited from Bloch functions
(BFs). In fact, electronic BFs are common eigenfunctions of
the Hamiltonian and translation operators that are conveniently
chosen as periodic functions in the reciprocal space. Hence,
taking normalization into account, the complex argument of
BFs is determined up to a position-independent shift that
should be periodic in the wave vector. Moreover, due to
its periodicity, the Bloch function of each band may be
represented by a Fourier series where coefficients are the
Wannier functions of the band.

Contrasting the extended nature of BFs, Wannier functions
are localized in real space and their nonuniqueness opens
the possibility of optimizing the localization. The case of
an electron in a one-dimensional (1D) crystal with inversion
symmetry was investigated by Kohn in his outstanding paper
of 1959.21 He proved that the phase of the BFs may be chosen
in order to produce real, symmetric or antisymmetric, and
exponentially localized WFs. The existence of an additional lo-
calization in the form of a power law was demonstrated in 2001
by He and Vanderbilt.22 This was extended in 2007 to deal with
maximally localized WFs of the noncentrosymmetric case.23

In the latter work, localization is measured by the variance, a
criterion that was successfully used in 1997 by Marzari and
Vanderbilt,24 to calculate WFs of three-dimensional crystals.
Their approach leads with a generalization of WFs, in the
sense established a few years after Kohn’s paper by Blount,25

des Cloizeaux,26 and Eilenberger,27 and became the basis of
modern computational packages, like WANNIER90,28 devoted
to electronic-structure calculations. Moreover, the 1D version
of those generalized functions has been investigated by using
iterative and parallel-transport approaches.24,29

The aim of this work is to speed advances in the study
of generalized Wannier functions (GWFs), by performing the
analytical optimization of their total quadratic spread in the
case of two consecutive bands of a 1D crystal with inversion
symmetry. Several questions about the obtained GWFs are
addressed, such as the possibility of being real, whether
they may present increased exponential localization when
compared with WFs of each band, how the inversion symmetry
manifests on their shape, and what is their relation with
atomic orbitals. In Sec. II, we review the construction of WFs
and maximally localized GWFs. The problem is analytically
solved in Sec. III, and numerical results and discussions for
a diatomic Kronig-Penney model are given in Sec. IV. Final
considerations are presented in Sec. V.

II. BASIC CONCEPTS AND EQUATIONS

Let us first recall basic ideas on Wannier functions of a
1D crystal. We denote a Bloch function of band index j and
wave number k as ψj,k(x). This family of functions satisfies
the Schrödinger equation Ĥψj,k(x) = Ej (k)ψj,k(x), the Bloch
condition ψj,k(x + a) = eika ψj,k(x), where a is the lattice
period, and the orthonormalization condition 〈ψ∗

j,k ψj ′,k′ 〉 =
2πa−1δj,j ′δ(k − k′), for any k and k′ in the first Brillouin
zone, where brackets refer to integration over the whole x

axis. Moreover, since ψj,k+2π/a(x) = ψj,k(x), BFs are given
by the Fourier series

ψj,k(x) =
∑

n

wj,n(x) eikna, (1)

where the coefficient wj,n(x) is the Wannier function of the
j th band and nth cell. This set of functions is orthonormal
over the x axis and satisfies wj,n(x) = wj (x − na) with

195127-11098-0121/2012/85(19)/195127(7) ©2012 American Physical Society

http://dx.doi.org/10.1103/PhysRevB.85.195127


DENIS R. NACBAR AND ALEXYS BRUNO-ALFONSO PHYSICAL REVIEW B 85, 195127 (2012)

wj (x) = ψj,k(x), where the line above the expression indicates
the mean value over a Brillouin zone. We do not lose generality
by assuming that Wannier functions of each band are real;23

i.e., BFs satisfy ψj,−k(x) = ψ∗
j,k(x). More importantly, we

suppose such WFs present maximal localization, i.e., minimal
variance.24

To deal with Bloch functions of a group of J bands, we
consider a row J -vector �k whose j th component is ψj,k . The
orthonormalization of the BFs reads 〈�†

k �k′ 〉 = 2πa−1δ(k −
k′)I. Moreover, we introduce the J -vector Wannier function
as W = �k . Its normalization is written as 〈W †W 〉 = I.

The centers of the WFs are the diagonal elements of the
following matrix:

〈W †xW 〉 = 〈x �
†
k �k′ 〉 = X(k), (2)

with

X(k) = i

∫ a

0
u
†
ku̇k dx (3)

being the Berry connection.1 Here uk = exp(−ikx)�k is the
periodic part of �k and a dot over a variable refers to
differentiation with respect to k. More clearly, the center of
wj (x) is xj = 〈W †xW 〉j,j = Xj,j (k). It may be shown that
X(k) is Hermitian and satisfies X(k + 2π/a) = X(k).

To find the variance of WFs we define the matrix

〈W †x2W 〉 = 〈x2 �
†
k �k′ 〉 = Y (k), (4)

where

Y (k) =
∫ a

0
u̇
†
k u̇k dx (5)

is a Hermitian matrix satisfying Y (k + 2π/a) = Y (k). Hence,
the variance of wj (x) is given by

σ 2
j = 〈W †x2W 〉j,j − x2

j = Yj,j (k) − (Xj,j (k))2. (6)

One may construct a set of J multiband Bloch functions
ψ̃j,k(x) by performing linear combinations of the J pure
functions ψj,k(x). To do so, we introduce a J × J matrix
U (k) such that24

ψ̃j,k(x) =
J∑

j ′=1

ψj ′,k(x) Uj ′,j (k) (7)

applies for j = 1, . . . ,J . These are also called quasi-Bloch
functions, since they satisfy the Bloch condition but are not
eigenfunctions of the Hamiltonian operator.

By introducing the row J -vector �̃k with j component
ψ̃j,k , we are able to rewrite Eq. (7) as �̃k = �k U (k).
Since 〈�̃†

k �̃k′ 〉 = 2πa−1 δ(k − k′)U †(k)U (k), Bloch func-
tions remain orthonormalized provided U (k) is unitary, i.e.,
U †(k)U (k) = I. Moreover, the periodicity of �̃k in k, with
period 2π/a, leads toU (k + 2π/a) = U (k). Therefore, GWFs
are given by W̃ = �̃k and satisfy the orthonormalization
condition 〈W̃ † W̃ 〉 = I.

The center of the j th generalized Wannier function is x̃j =
〈W̃ †xW̃ 〉j,j = X̃j,j (k), where

X̃(k) = U †(k)X(k)U (k) + i U †(k)U̇(k) (8)

is a Hermitian matrix. Moreover, the variance of w̃j (x) is given
by

σ̃ 2
j = 〈W̃ †x2W̃ 〉j,j − x̃2

j = Ỹj,j (k) − (X̃j,j (k))2, (9)

where

Ỹ (k) = U †(k)Y (k)U (k) + U̇ †(k)U̇(k)

+ i[U †(k)X(k)U̇(k) − U̇ †(k)X(k)U (k)]

= U †(k)[Y (k) − X2(k)]U (k) + X̃†(k)X̃(k). (10)

To improve the quality of the generalized Wannier func-
tions, we minimize their total quadratic spread, i.e, the sum of
their variances,24 which is given by

� =
J∑

j=1

σ̃ 2
j = Tr[Ỹ (k)] −

J∑
j=1

x̃2
j . (11)

Since similar matrices have the same trace, this functional
takes the form

� = Tr[Y (k) − X2(k) + X̃†(k)X̃(k)] −
J∑

j=1

[X̃j,j (k)] 2

= �0 + 2
J∑

j=1

j−1∑
j ′=1

|X̃j,j ′ |2 +
J∑

j=1

[
X̃2

j,j − (X̃j,j ) 2
]
, (12)

where �0 = Tr[Y − X2] does not depend on U . The second
and third terms in the second line of Eq. (12) are non-negative.
Therefore, �0 is a lower bound for the total spread. Moreover,
if �0 is the minimum, then X̃(k) should be a real and constant
diagonal matrix. According to Eq. (8), the optimization
problem reduces to finding diagonal elements of X̃ such that a
unitary transformation having the periodicity of the reciprocal
lattice solves the first-order differential equation

U̇(k) = i[X(k)U (k) − U (k) X̃]. (13)

The resulting GWFs are real if the solution satisfies

U (−k) = U∗(k). (14)

We note that expressions in Eq. (12) resemble the previous
approach in Ref. 24, though we limit ourselves to the particular
case of a 1D crystal. After some algebraic work, it is possible
to prove that Eqs. (14), (19), and (20) in Ref. 24 correspond
to the first, second, and third terms, respectively, in the second
line of the present Eq. (12). In this latter version, however,
only matrix elements of finite matrices Y (k), X(k), and X̃(k)
are involved, the invariance of the first term is more apparent,
and the minimum variance is simply given in terms of Y (k)
and X(k). Moreover, for the case of a general 1D crystal,
it has been shown that eigenfunctions of the band-projected
position operator are maximally localized GWFs.24 Therefore,
diagonalization of that operator would lead to an optimal
matrix U . Instead, in next section, such a matrix is obtained
by solving Eq. (13), and analytic solutions are given for two
consecutive bands in a 1D crystal having inversion symmetry.
Later on, such solutions are compared with those of the
eigenvalue problem for the band-projected position operator.

195127-2



ANALYTICAL OPTIMIZATION OF SPREAD AND CHANGE . . . PHYSICAL REVIEW B 85, 195127 (2012)

III. EQUATIONS FOR TWO BANDS

Now we limit ourselves to the case of two bands, i.e., J = 2.
Since U (k) is unitary, we have det[U (k)] = exp[i τ (k)],
where τ (k) is a real function. Moreover, we obtain
U1,2 = −U ∗

2,1 exp(iτ ), U2,2 = U ∗
1,1 exp(iτ ), and | det[U ]| =

|U1,1|2 + |U2,1|2 = 1. Thus, there exist real functions θ (k),
φ1,1(k), and β(k), such that U1,1 = cos(θ ) exp(iφ1,1) and
U2,1 = sin(θ ) exp[i(φ1,1 + β)]. Since |U2,2| = |U1,1|, we
conveniently denote the phase difference between U2,2 and
U1,1 as α(k). Hence, matrix U may be written as

U = e
iτ
2

(
cos(θ ) − sin(θ ) e−iβ

sin(θ ) eiβ cos(θ )

)(
e− iα

2 0

0 e
iα
2

)
. (15)

We assume that the WFs of each single band have been
already optimized. Consequently,23 X1,1 = x1 and X2,2 = x2.
Also, since X is Hermitian, diagonal terms are real, while off-
diagonal ones satisfy X1,2 = X∗

2,1, with X2,1 = C exp(i γ ),
where C(k) and γ (k) are real functions. Hence, assumingU (k)
leads to a global minimum �0, we have X̃1,1 = x̃1, X̃2,2 = x̃2,
and X̃1,2 = X̃2,1 = 0. In this way, the sum of variances of the
maximally localized GWFs is

σ̃ 2
1 + σ̃ 2

2 = �0 = σ 2
1 + σ 2

2 − 2 C2. (16)

Moreover, Eq. (13) leads to a rather simple system of
differential equations, namely,

τ̇ = (x1 + x2) − (x̃1 + x̃2), (17)

θ̇ = C sin(β − γ ), (18)

α̇ = �x − �x̃ − 2C cos(β − γ ) tan(θ ), (19)
and

β̇ = �x + 2C cos(β − γ ) cot(2θ ), (20)

where �x = x2 − x1 and

�x̃ = x̃2 − x̃1

= �x − 2 C cos(β − γ ) tan(θ )

− a

2π
[α(π/a) − α(−π/a)]. (21)

The periodicity of U leads to conditions τ (k +
2π/a) = τ (k) + 2πμτ , θ (k + 2π/a) = (−1)μβ θ (k) + πμθ ,
α(k + 2π/a) = α(k) + 2πμα , and β(k + 2π/a) = β(k) +
πμβ , where μτ , μθ , μα , and μβ are integers, with μα −
μτ and μθ having the same parity. Moreover, Eq. (14)
implies τ (−k) = −τ (k) + 2πντ , θ (−k) = (−1)νβ θ (k) + πνθ ,
α(−k) = −α(k) + 2πνα , and β(−k) = −β(k) + πνβ , where
ντ , νθ , να , and νβ are integers, with να − ντ and νθ having the
same parity.

At this point we note that the expression in Eq. (15) for
transformation matrix U has the advantage of the differential
equation for τ (k) being decoupled. Solving Eq. (17) with the
corresponding boundary conditions leads to τ (k) = aμτk +
πντ , with μτ = [(x1 + x2) − (x̃1 + x̃2)]/a. Then, by taking
μτ = 0, the condition x̃1 + x̃2 = x1 + x2 holds; i.e., the sum
of centers of Wannier functions is conserved.24 Additionally,
by choosing ντ = 0, τ (k) identically vanishes. In other words,
the matrix U belongs to the special unitary group SU (2). Of

course, in agreement with these choices, μα and μθ should
have equal parities. The same thing applies to να and νθ .

Now we further restrict the discussion to the case of a 1D
crystal with inversion symmetry and two consecutive bands
labeled by j = 1 and j = 2. As shown in Appendix A, we
may take γ (k) = k �x + sπ/2, where s = 0 (s = 1) when the
parities of w1(x) and w2(x) are different (equal). Moreover,
C(k) is even (odd) when s = 0 (s = 1), and C(k) is periodic
(antiperiodic) over the reciprocal lattice, when n = 2�x/a is
an even (odd) integer. At this point we note that Bloch functions
at the edges of a band gap have different parities.30 Therefore,
the case s = 1 and even n does not occur for consecutive
bands.

Since γ̇ = �x, Eqs. (18)–(20) may be rewritten as

θ̇ = C sin(ξ ), (22)

α̇ = 2C cos(ξ ) tan(θ ) − 2C cos(ξ ) tan(θ ) + aμα, (23)

and

ξ̇ = 2C cos(ξ ) cot(2θ ), (24)

respectively, where ξ = β − γ . By combining Eqs. (22) and
(24), it is easily shown that cos(ξ ) sin(2θ ) = A, where A is a
constant satisfying −1 � A � 1. In what follows, only cases
A = −1 and A = 0 need to be considered.

On the one hand, A = −1 solves the problem pro-
vided θ (k) = −π/4, α̇ = 2(C − C) + aμα , ξ (k) = 0 [i.e.,
β(k) = γ (k)], μθ = [(−1)n − 1]/4, μβ = n = 2�x/a, νθ =
[(−1)s − 1]/4, and νβ = s. Therefore, this solution applies
for s = 0 and even n only, with μθ = νθ = νβ = 0. Moreover,
taking μα = να = 0, with the same parity as μθ and νθ ,
Eq. (23) is solved to obtain

α(k) = 2
∫ k

0
[C(κ) − C] dκ, (25)

where C is the mean value of C(k). Furthermore, �x̃ =
�x + 2C leads to x̃1 = x1 − C and x̃2 = x2 + C. The optimal
transformation matrix is

U = 1√
2

(
1 e−inka/2

−einka/2 1

) (
e− iα

2 0

0 e
iα
2

)
. (26)

On the other hand, A = 0 solves the problem for odd
values of n = 2�x/a, if θ̇ = −C, μθ = νθ = 0, α̇ = a μα ,
and ξ (k) = −π/2, i.e., β(k) = γ (k) − π/2, with μβ = n and
νβ = s − 1. Since μα and να should be even, we take μα =
να = 0, which leads to α(k) = 0 for every k. Moreover,

θ (k) =
∫ sπ/a

k

C(κ) dκ (27)

is odd (even) for s = 0 (s = 1), and the difference between
centers is preserved, i.e., �x̃ = �x. Since the same happens
with their sum, both centers remain untouched (x̃1 = x1 and
x̃2 = x2). The transformation matrix becomes

U =
(

cos(θ ) −i1−s sin(θ ) e−inka/2

is−1 sin(θ ) einka/2 cos(θ )

)
. (28)

It is worth noticing that we have obtained the forms of
matrix U (k) leading to a constant, real diagonal matrix X̃(k),

195127-3



DENIS R. NACBAR AND ALEXYS BRUNO-ALFONSO PHYSICAL REVIEW B 85, 195127 (2012)

for two consecutive bands, in the three symmetry cases of
interest. This guarantees �̃ reaches its global minimum, which
is given by Eq. (16).

IV. NUMERICAL RESULTS AND DISCUSSION

To be concrete, we deal with a diatomic version of the
Kronig-Penney model,31,32 where a particle of mass m has
potential energy given by

V (x) = − vh̄2

m a

∑
n

[δ(x − na + b/2) + δ(x − na − b/2)].

(29)

The index n runs from −∞ to +∞, v > 0 is the dimensionless
interaction strength, and b < a is an interatomic distance.
The energy bands are obtained from E = h̄2ε/(2ma2) and
cos(ka) = μ(ε), where

μ(ε) = (1 − v2/ε) cos(
√

ε) − 2v sin(
√

ε)/
√

ε

+ v2 cos[(1 − 2b/a)
√

ε]/ε. (30)

Since the potential presents inversion centers at integer multi-
ples of a/2, phases of Bloch functions are chosen appropriately
in order to obtain maximally localized Wannier functions for
each isolated band.30 The numerical results below correspond
to parameters v = 4 and b = 3a/8. Moreover, energies and
Bloch functions are calculated over a sample of 400 wave
vectors uniformly spaced in the first Brillouin zone.

Figures 1(a) and 1(b) display the WFs of the lower two
bands. Hence, the functions are labeled by superscript (1,2).
It is apparent that single-band Wannier functions resemble
molecular orbitals, with the first (second) band corresponding
to the bonding (antibonding) one. Both molecular orbitals
are centered at the origin of coordinates, i.e., x1 = x2 = 0,
and the first (second) one is even (odd), corresponding to
the case s = 0 and n = 0 in Sec. III. Moreover, the standard
deviation of the WF for the first (second) band is σ1 ≈ 0.2579a

(σ2 ≈ 0.3323a). These functions decay exponentially with
coefficients that may be given in terms of the quantities

FIG. 1. (Color online) Wannier functions of (a) first and
(b) second bands and (c,d) generalized Wannier functions for the
pair (1,2) of bands. Open circles indicate atomic positions.

hj = a−1 arccos[μ(εj )], where εj is the j th root of μ′(ε) = 0.
For parameters in the figure, we have h1 ≈ 0.9969 a−1 and
h2 ≈ 1.3542 a−1. Hence,21 the decay of the WF of the first
(second) band is given by h1 (the minimum between h1 and
h2, i.e., h1). This corresponds to the distance from the real
k axis to the closest branch points of εj (k). Of course, an
additional power-law decay with exponent 3/4 also occurs.22

Generalized Wannier functions for the lower two bands
are calculated by the transformation matrix in Eq. (26). Since
n = 0, it takes a simpler form, namely,

U (k) = 1√
2

(
e− iα(k)

2 e
iα(k)

2

−e− iα(k)
2 e

iα(k)
2

)
, (31)

where α(k), given by Eq. (25), is an odd function. Figures 1(c)
and 1(d) display results for the same parameters as Figs. 1(a)
and 1(b). The centers of these functions are x̃1 ≈ −0.21125a

and x̃2 ≈ 0.21125a. Moreover, they present identical standard
deviations, namely σ̃1 = σ̃2 ≈ 0.1516a. This means that the
total variance has decreased by about 74%. It is worth noting
that these GWFs are neither symmetric nor antisymmetric.
Indeed, crystal symmetry manifests in the sense that they are
mirror images of each other about line x = 0, i.e., w̃

(1,2)
2 (x) =

w̃
(1,2)
1 (−x). As shown in Appendix B, this follows from choices

leading to Eq. (31) and gauge choices21 ψ1,k(−x) = ψ1,−k(x)
and ψ2,k(−x) = −ψ2,−k(x). We also note that, because of the
diagonal character of X̃(k), the condition 〈w̃(1,2)

1 |x|w̃(1,2)
2 〉 = 0

should apply. The symmetry relation between the GWFs also
leads to this result, namely,〈

w̃
(1,2)
1

∣∣x∣∣w̃(1,2)
2

〉 =
∫ ∞

−∞
x w̃

(1,2)
1 (x)w̃(1,2)

2 (x) dx

=
∫ ∞

−∞
x w̃

(1,2)
1 (x)w̃(1,2)

1 (−x) dx = 0, (32)

since w̃
(1,2)
1 (x)w̃(1,2)

1 (−x) is an even function.
Regarding the asymptotic behavior of GWFs, symmetry

allows us to focus on w̃
(1,2)
1 (x) only. We notice that both

the coefficient of exponential decay and the exponent of
the power-law decay may be extracted from the list of
probabilities for the particle to be found within each unit
cell.22,23 For the functions in the figure, we have obtained
identical parameters for x → −∞ and x → ∞. Namely, the
coefficient of exponential localization is noticeably larger for
these generalized WFs than for functions in Figs. 1(a) and 1(b),
and it is found to be h2 ≈ 1.3542. This is due to cancellation
between terms in the linear combination of ψ1,k and ψ2,k at
their common branch points, while farther branch points of
ψ2,k survive. The 3/4 power-law decay also occurs.

When a minimization procedure of total variance is applied
to the pair (2,3) of bands, the results are quite different.
This is because the symmetries of WFs of second and third
bands are different, while the distance between their centers
is an odd multiple of half period. In fact, it is apparent in
Figs. 2(a) and 2(b) that x1 = 0 and x2 = a/2. Moreover, the
variance and coefficient of exponential localization for the
second (third) band are σ1 ≈ 0.3323a and h1 = 0.9968a−1

(σ2 ≈ 0.3468a and h3 = 0.7994a−1), respectively. This case
is characterized by parameters s = 1 and n = 1 of Sec. III.
Hence, the transformation in Eq. (28) is used to obtain

195127-4



ANALYTICAL OPTIMIZATION OF SPREAD AND CHANGE . . . PHYSICAL REVIEW B 85, 195127 (2012)

FIG. 2. (Color online) Wannier functions of (a) second and
(b) third bands and (c,d) generalized Wannier functions for the pair
(2,3) of bands. Open circles indicate atomic positions.

the generalized functions shown in Figs. 2(c) and 2(d). By
comparing WFs and GWFs, we observe that the centers and
symmetries are preserved, while the total spread diminishes.
On the one hand, the properties w̃

(2,3)
1 (−x) = −w̃

(2,3)
1 (x) and

w̃
(2,3)
2 (a − x) = w̃

(2,3)
2 (x) are shown in Appendix B to follow

from Eq. (28) and gauge choices21 ψ2,k(−x) = −ψ2,−k(x) and
ψ3,k(a − x) = ψ3,−k(x). On the other hand, the variances of
w

(2,3)
1 (x) and w

(2,3)
2 (x) are σ̃1 ≈ 0.3157a and σ̃2 ≈ 0.2108a,

respectively, corresponding to a decrease in total variance of
about 37%. Furthermore, both GWFs decay exponentially with
the same coefficient, namely, h3 = 0.7994 a−1. To understand
this, it may be noticed that the first, second, and third gaps
with coefficients h1, h2, and h3 are relevant for the pair (2,3)
of bands. Therefore, even if cancellation at branch points of the
second gap were to occur, min(h1,h3) = h3 would dominate
the decay.

The results are similar for the pair (3,4) of bands. As
shown in Figs. 3(a) and 3(b), now single-band WFs have

FIG. 3. (Color online) Wannier functions of (a) third and
(b) fourth bands and (c,d) generalized Wannier functions for the pair
(3,4) of bands. Open circles indicate atomic positions.

the same parity but their centers are separated by half
a period. Therefore, the transformation in Eq. (28) leads
to the maximally localized GWFs displayed in Figs. 3(c)
and 3(d). Clearly, centers and symmetries are preserved
again. Namely, the conditions w̃

(3,4)
1 (a − x) = w̃

(3,4)
1 (x) and

w̃
(3,4)
2 (−x) = w̃

(3,4)
2 (x) hold. These properties may be derived

from Eq. (28) and gauge choices21 ψ3,k(a − x) = ψ3,−k(x) and
ψ4,k(−x) = ψ4,−k(x) (see Appendix B). Moreover, standard
deviations of WFs (σ1 = 0.3468a, σ2 = 1.1982a) and GWFs
(σ̃1 = 0.6915a, σ̃2 = 0.9618a), corresponding to a decrease in
total spread of about 10%. It is interesting to note in Fig. 3 that
one function has become broader [Fig. 3(a)→Fig. 3(c)] while
the other has narrowed [Fig. 3(b)→Fig. 3(d)]. At the same
time, the first one has slowed down its exponential decrease.
In fact, the coefficient of decay of both GWFs is found to be
h4 ≈ 0.0545 a−1.

For the sake of completeness, we now compare the results
presented above with those obtained by diagonalization of
the position operator x. Conveniently, this operator can be
represented in terms of Wannier functions wj,n(x), with j =
1, . . . ,J and n running over the integers. The matrix elements
of x in this basis are given by

x(j,n),(j ′,n′) = 〈wj,n|x|wj ′,n′ 〉 = Xj,j ′,n′−n + na δj,j ′ δn,n′ ,

(33)

with

Xj,j ′,n = exp(−inak) Xj,j ′ (k) (34)

being the nth Fourier coefficient of Xj,j ′ (k). Clearly, the
matrix is Hermitian, and its eigenvalues are the centers
of GWFs with maximal localization.24 It is worth noting
that the spectrum of x is periodic, with period a, and the
corresponding eigenvectors are the coordinates of GWFs
in the basis set of single-band WFs. Therefore, it suffices
to accurately determine those eigenvalues between −a/2
and a/2. With this in mind, approximate eigenvalues may
be calculated by truncating the basis, namely, by taking
n = −N,1 − N, . . . , − 1,0,1, . . . ,N , for a sufficiently large
integer N . The truncated matrix is of dimension J (2N + 1) ×
J (2N + 1), and accuracy may be improved by increasing
the value of N . This is apparent in Fig. 4, where the two
relevant eigenvalues for the pair (1,2) of bands are shown as a
function of N . In fact, the lower [upper] eigenvalue in Fig. 4(a)

FIG. 4. (Color online) The two eigenvalues, within the range
−a/2 < x � a/2, of the position operator, when projected on
the set of Wannier functions wj,n(x), with j = 1 or 2 and n =
−N,1 − N, . . . ,N . Dots denote the lower (upper) eigenvalue, with
N = 0,1, . . . ,6. Horizontal dashed lines are at the centers of the
generalized Wannier functions in Fig. 1.

195127-5



DENIS R. NACBAR AND ALEXYS BRUNO-ALFONSO PHYSICAL REVIEW B 85, 195127 (2012)

[Fig. 4(b)] rapidly converges to the center of w̃1(x) [w̃2(x)]
as N increases. It should also be remarked that eigenvalues
for N = 0 coincide with the corresponding limit values. This
is in agreement with the analytical development in Sec. II,
where x̃1 and x̃2 have been shown to be x1 − C and x2 + C,
respectively, with x1 = x2 = 0 and C = x(1,0),(2,0) = x(2,0),(1,0).
We have also checked the convergence to GWFs of linear
combinations of the J (2N + 1) WFs, with the coefficients
being the corresponding eigenvectors. For the considered
parameters of the 1D crystal, the scalar product between the
analytically obtained GWFs and the linear combinations is
above 0.99, for all values of N . And, apparently, convergence
to 1 is quite rapid. Nevertheless, good agreement between
the exponentially decaying tails may require large values
of N . A similar analysis may be done for other pairs of
bands.

V. CONCLUSIONS

The total variance of generalized Wannier functions of two
consecutive bands in a 1D crystal with inversion symmetry
has been minimized by following a straightforward analytical
approach. This has led to real generalized Wannier functions of
maximal localization. Different expressions have been given
for the cases where the distance between the centers of
single-band Wannier functions is either an even or an odd
multiple of half a period of the crystal. Numerical results were
presented for the lowest four bands of a particle subject to
a diatomic Kronig-Penney potential. For the lowest couple of
bands, calculated single-band (generalized) Wannier functions
resemble molecular (orthogonalized atomic) orbitals, and the
latter display increased exponential localization. Moreover,
the centers change and the generalized functions are neither
symmetric nor antisymmetric, but they are related by an
inversion-symmetry transformation. Instead, for the second
and third pairs of consecutive bands, Wannier functions con-
serve both center and inversion symmetry, while exponential
decay does not increase. Changes in the coefficient of such
decay have been explained in terms of branch points of Bloch
functions of each band.

The presented explicit expressions of the matrix that
transforms the Bloch function into quasi-Bloch functions
noticeably simplify computation of generalized Wannier
functions, making them available to a broader community.
Moreover, convergence of truncation and iterative procedures
may be further tested. This has been done for the eigenvalue
problem of the band-projected position operator. At the same
time, interesting properties of symmetry and decay have been
explained. Calculated functions should find applications in
the study of defect-induced localized states with energy lying
within the gap separating the considered pair of bands, but
future work should deal with generalized Wannier functions
of several bands. We expect similar efforts will be devoted to
one-dimensional crystals lacking inversion symmetry, as well
as two- and three-dimensional crystals.

ACKNOWLEDGMENTS

We thank the Brazilian Agencies CNPq and FAPESP for
financial support.

APPENDIX A: OFF-DIAGONAL TERM X2,1(k)

In this section, properties of real functions C(k) and γ (k)
satisfying X2,1(k) = C(k) exp[i γ (k)] are demonstrated. From
Eq. (3), we have

X2,1(k) = i

∫ a

0
u∗

2,k(x)u̇1,k(x) dx

=
∫ a

0
ψ∗

2,k(x)[i ψ̇1,k(x) + x ψ1,k(x)] dx. (A1)

If wj (x) obeys wj (x) = sj wj (2xj − x), then it is even (odd)
about x = xj , when sj = 1 (sj = −1). The correspond-
ing Bloch function satisfies ψj,k(x) = sj ψj,−k(2xj − x) =
sj ψ∗

j,k(2xj − x). Taking the Bloch condition and orthogonal-
ity of Bloch functions into account, we obtain

X2,1(k) = s1s2

∫ 2x1

2x1−a

ψ2,k(2x2 − x)[i ψ̇∗
1,k(2x1 − x)

+ x ψ∗
1,k(2x1 − x)] dx

= s1s2 e2ik �x

∫ a

0
ψ2,k(x)[i ψ̇∗

1,k(x)

+ (2x1 − x) ψ∗
1,k(x)] dx

= −s1s2 e2ik �x X∗
2,1(k). (A2)

Therefore, we may write

e−ik �x X2,1(k) = (−1)s [e−ik �x X2,1(k)]∗, (A3)

where (−1)s = −s1s2, and we take s = 0 (s = 1) when the
parities of w1(x) and w2(x) are different (coincide). In the
first (second) case, e−ik �x X2,1(k) is real [pure imaginary]
and may be written as C(k) [i C(k)]. These results may be
summarized as X2,1(k) = C(k) exp[i γ (k)], where C(k) is real
and γ (k) = k �x + s π/2.

From Eq. (A1) we may derive inversion-symmetry prop-
erties of C(k). In fact, since X2,1(−k) = X∗

2,1(k), we obtain
C(−k) = (−1)s C(k). This means C(k) is even (odd) for s = 0
(s = 1). Furthermore, since X2,1(k) has the periodicity of the
reciprocal lattice, C(k) satisfies C(k + 2π/a) = (−1)n C(k),
with n = 2�x/a; i.e., it is periodic (antiperiodic) when n is
even (odd).

APPENDIX B: SYMMETRY AND GENERALIZED
WANNIER FUNCTIONS

To obtain a symmetry relation between the generalized
Wannier functions of a couple of bands having Wannier func-
tions centered at x = 0 with the first (second) one being an even
(odd) function, we take into account conditions ψ1,k(−x) =
ψ1,−k(x) and ψ2,k(−x) = −ψ2,−k(x). From Eq. (31), we have

ψ̃
(1,2)
1,k (x) = e− i α(k)

2√
2

[ψ1,k(x) − ψ2,k(x)] (B1)

and

ψ̃
(1,2)
2,k (x) = e

i α(k)
2√
2

[ψ1,k(x) + ψ2,k(x)]. (B2)

195127-6



ANALYTICAL OPTIMIZATION OF SPREAD AND CHANGE . . . PHYSICAL REVIEW B 85, 195127 (2012)

Moreover, according to Eq. (25), α(k) is an odd function, and
we obtain

ψ̃
(1,2)
2,k (x) = e

−i α(−k)
2√
2

[ψ1,−k(−x) − ψ2,−k(−x)] = ψ̃
(1,2)
1,−k(−x).

(B3)

Taking the average over the first Brillouin zone, this leads to
w̃

(1,2)
2 (x) = w̃

(1,2)
1 (−x).

Now, let us try to demonstrate that centers and symmetries
of w̃1(x) and w̃2(x) coincide with those of w1(x) and w2(x),
respectively, when n = 2�x/a is odd. To be concrete, we con-
sider the pair (2,3) in Sec. IV, where ψ2,k(−x) = −ψ2,−k(x)
and ψ3,k(a − x) = ψ3,−k(x). We note that s = 0 (opposite
symmetries), θ (k) is odd, and n = 1. Then, from Eq. (28),
we have

ψ̃
(2,3)
1,k (x) = cos (θ (k)) ψ2,k(x) − i e

ika
2 sin (θ (k)) ψ3,k(x)

(B4)

and

ψ̃
(2,3)
2,k (x) = −i e− ika

2 sin (θ (k)) ψ2,k(x) + cos (θ (k)) ψ3,k(x).

(B5)

According to the Bloch condition, this leads to

ψ̃
(2,3)
1,k (−x) = −cos (−θ (−k)) ψ2,−k(x)

− i e
−ika

2 sin(−θ (−k)) ψ3,−k(x)

= − cos (θ (−k)) ψ2,−k(x)

+ i e
−ika

2 sin (θ (−k)) ψ3,−k(x)

= −ψ̃
(2,3)
1,−k(x) (B6)

and

ψ̃
(2,3)
2,k (a − x) = i e

ika
2 sin (−θ (−k)) ψ2,−k(x)

+ cos (−θ (−k)) ψ3,−k(x)

= −i e
ika
2 sin (θ (−k)) ψ2,−k(x)

+ cos (θ (−k)) ψ3,−k(x)

= ψ̃
(2,3)
2,−k(x). (B7)

Hence, after averaging over the first Brillouin zone, the results
are w̃

(2,3)
1 (−x) = −w̃

(2,3)
1 (x) and w̃

(2,3)
2 (a − x) = w̃

(2,3)
2 (x).

The conservation of symmetry for pair (3,4) in Sec. IV may
be proved in the same way. In doing so, the conditions ψ3,k(a −
x) = ψ3,−k(x) and ψ4,k(−x) = ψ4,−k(x) should be taken into
account. Moreover, we note that s = 1 (same symmetry), θ (k)
is even, and n = −1.

1R. Resta, Eur. Phys. J. B 79, 121 (2011).
2E. Krüger and H. Strunk, J. Supercond. Novel Mag. 24, 2103
(2011).

3T. Berlijn, D. Volja, and W. Ku, Phys. Rev. Lett. 106, 077005 (2011).
4T. Bazhirov, J. Noffsinger, and M. L. Cohen, Phys. Rev. B 84,
125122 (2011).

5T. Stollenwerk, R. Frank, A. Lubatsch, O. Zaitsev, S. Zhukovsky,
D. Chigrin, and J. Kroha, Appl. Phys. B 105, 163 (2011).

6M. Lorenz, D. Usvyat, and M. Schütz, J. Chem. Phys. 134, 094101
(2011).

7D. Benson, O. F. Sankey, and U. Häussermann, Phys. Rev. B 84,
125211 (2011).

8J. S. Lee, G. A. H. Schober, M. S. Bahramy, H. Murakawa,
Y. Onose, R. Arita, N. Nagaosa, and Y. Tokura, Phys. Rev. Lett.
107, 117401 (2011).

9M. Casula, M. Calandra, G. Profeta, and F. Mauri, Phys. Rev. Lett.
107, 137006 (2011).

10M. Shelley, N. Poilvert, A. A. Mostofi, and N. Marzari, Comput.
Phys. Commun. 182, 2174 (2011).

11M. Chen, G.-C. Guo, and L. He, J. Phys.: Condens. Matter 23,
325501 (2011).

12N. W. Ashcroft and N. D. Mermin, Solid State Physics (Brooks
Cole, Belmont, 1976).

13O. Madelung, Introduction to Solid-State Physics (Springer, Berlin,
1978).

14R. M. Martin, Basic Electronic Structure (Cambridge University
Press, Cambridge, 2004).

15Conceptual Foundations of Materials, edited by S. G. Louie and
M. L. Cohen (Elsevier, Amsterdam, 2006).

16Condensed Matter Physics, 2nd ed., edited by M. P. Marder (Wiley,
Hoboken, 2010).

17G. H. Wannier, Phys. Rev. 52, 191 (1937).
18G. Nenciu, Commun. Math. Phys. 91, 81 (1983).
19C. Brouder, G. Panati, M. Calandra, C. Mourougane, and

N. Marzari, Phys. Rev. Lett. 98, 046402 (2007).
20H. D. Cornean, A. Nenciu, and G. Nenciu, J. Phys. A 41, 125202

(2008).
21W. Kohn, Phys. Rev. 115, 809 (1959).
22L. He and D. Vanderbilt, Phys. Rev. Lett. 86, 5341 (2001).
23A. Bruno-Alfonso and D. R. Nacbar, Phys. Rev. B 75, 115428

(2007).
24N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847 (1997).
25E. I. Blount, Solid State Phys. 13, 305 (1962).
26J. D. Cloizeaux, Phys. Rev. 135, A685 (1964); 135, A698 (1964).
27G. Eilenberger, Z. Phys. A 180, 43 (1964).
28A. A. Mostofi, J. R. Yates, Y.-S. Lee, I. Souza, D. Vanderbilt, and

N. Marzari, Comput. Phys. Commun. 178, 685 (2008).
29J. Bhattacharjee and U. V. Waghmare, Phys. Rev. B 71, 045106

(2005).
30A. Bruno-Alfonso and G.-Q. Hai, J. Phys.: Condens. Matter 15,

6701 (2003).
31F. B. Pedersen, G. T. Einevoll, and P. C. Hemmer, Phys. Rev. B 44,

5470 (1991).
32M. J. Rave and W. C. Kerr, Eur. Phys. J. B 45, 473 (2005).

195127-7

http://dx.doi.org/10.1140/epjb/e2010-10874-4
http://dx.doi.org/10.1007/s10948-011-1177-6
http://dx.doi.org/10.1007/s10948-011-1177-6
http://dx.doi.org/10.1103/PhysRevLett.106.077005
http://dx.doi.org/10.1103/PhysRevB.84.125122
http://dx.doi.org/10.1103/PhysRevB.84.125122
http://dx.doi.org/10.1007/s00340-011-4729-3
http://dx.doi.org/10.1063/1.3554209
http://dx.doi.org/10.1063/1.3554209
http://dx.doi.org/10.1103/PhysRevB.84.125211
http://dx.doi.org/10.1103/PhysRevB.84.125211
http://dx.doi.org/10.1103/PhysRevLett.107.117401
http://dx.doi.org/10.1103/PhysRevLett.107.117401
http://dx.doi.org/10.1103/PhysRevLett.107.137006
http://dx.doi.org/10.1103/PhysRevLett.107.137006
http://dx.doi.org/10.1016/j.cpc.2011.05.017
http://dx.doi.org/10.1016/j.cpc.2011.05.017
http://dx.doi.org/10.1088/0953-8984/23/32/325501
http://dx.doi.org/10.1088/0953-8984/23/32/325501
http://dx.doi.org/10.1103/PhysRev.52.191
http://dx.doi.org/10.1007/BF01206052
http://dx.doi.org/10.1103/PhysRevLett.98.046402
http://dx.doi.org/10.1088/1751-8113/41/12/125202
http://dx.doi.org/10.1088/1751-8113/41/12/125202
http://dx.doi.org/10.1103/PhysRev.115.809
http://dx.doi.org/10.1103/PhysRevLett.86.5341
http://dx.doi.org/10.1103/PhysRevB.75.115428
http://dx.doi.org/10.1103/PhysRevB.75.115428
http://dx.doi.org/10.1103/PhysRevB.56.12847
http://dx.doi.org/10.1016/S0081-1947(08)60459-2
http://dx.doi.org/10.1103/PhysRev.135.A685
http://dx.doi.org/10.1103/PhysRev.135.A698
http://dx.doi.org/10.1016/j.cpc.2007.11.016
http://dx.doi.org/10.1103/PhysRevB.71.045106
http://dx.doi.org/10.1103/PhysRevB.71.045106
http://dx.doi.org/10.1088/0953-8984/15/40/008
http://dx.doi.org/10.1088/0953-8984/15/40/008
http://dx.doi.org/10.1103/PhysRevB.44.5470
http://dx.doi.org/10.1103/PhysRevB.44.5470
http://dx.doi.org/10.1140/epjb/e2005-00204-6