# 1 Introduction

## 1 Introduction

The dipole moment of any finite $`N`$–electron system in its ground state is a simple and well defined quantity. Given the many–body wavefunction $`\mathrm{\Psi }`$ and the corresponding single–particle density $`n(𝐫)`$ the electronic contribution to the dipole is:

$$𝐑=𝑑𝐫𝐫n(𝐫)=\mathrm{\Psi }|\widehat{𝐑}|\mathrm{\Psi },$$
(1)   
where $`\widehat{𝐑}=_{i=1}^N𝐫_i`$ (atomic Hartree units are adopted throughout). This looks very trivial, but we are exploiting here an essential fact: the ground wavefunction of any finite $`N`$–electron system is square–integrable and vanishes exponentially at infinity; the density vanishes exponentially as well.

Considering now a macroscopic solid, the related quantity is macroscopic polarization, which is a very essential concept in any phenomenological description of dielectric media : this quantity is ideally defined as the dipole of a macroscopic sample, divided by its volume. The point is that, when using Eq. (1), the integral is dominated by what happens at the surface of the sample: knowledge of the electronic distribution in the bulk region is not enough to unambiguosly determine the dipole. This looks like a paradox, since in the thermodynamic limit macroscopic polarization must be an intensive quantity, insensitive to surface effects.

Macroscopic polarization in the bulk region of the solid must be determined by what “happens” in the bulk as well. This is the case if one assumes a model of discrete and well separated dipoles, à la Clausius-Mossotti: but real dielectrics are very much different from such an extreme model. The valence electronic distribution is continuous, and often very delocalized (particularly in covalent dielectrics). Most textbooks attempt at explaining the polarization of a periodic crystal via the dipole moment of a unit cell, or something of the kind . These definitions are incorrect : according to the modern viewpoint, bulk macroscopic polarization is a physical observable completely independent from the periodic charge distribution of the polarized crystalline dielectric.

In condensed matter physics the standard way for getting rid of undesired surface effects is to adopt periodic Born-von Kármán boundary conditions (BvK). Indeed, the BvK choice is mandatory in order to introduce even the most elementary topics, such as the free–electron gas and its Fermi energy, or the Bloch theorem . Unfortunately, the adoption of BvK does not solve the polarization problem. In fact the dipole cannot be evaluated as in Eq. (1) when the wavefunction obeys BvK: the integrals are ill–defined due to the unbounded nature of the quantum–mechanical position operator.

For this reason macroscopic polarization remained a major challenge in electronic structure theory for many years. The breakthrough came in 1992, when polarization was defined in terms of the wavefunctions, not of the charge. This definition has an unambiguous thermodynamic limit, such that BvK and Bloch states can be used with no harm . In the following months a modern theory of macroscopic polarization in crystalline dielectrics has been completely established , thanks to a major advance due to R.D. King-Smith and D. Vanderbilt , who expressed polarization in terms of a Berry phase . A comprehensive account of the modern theory exists . Other less technical presentations are available as well ; for an oversimplified nontechnical outline see Ref. . First–principle calculations based on this theory have been performed for several crystalline materials, either within DFT (using various basis sets) or within HF (using LCAO basis sets) .

All of the above quoted work refers to a crystalline system within an independent–electron formulation: the single–particle orbitals have then the Bloch form, and the macroscopic polarization is evaluated as a Berry phase of them. The related, but substantially different, problem of macroscopic polarization in a correlated many–electron system was first solved by Ortíz and Martin in 1994 . However, according to them, polarization is defined—and computed —by means of a peculiar “ensemble average”, integrating over a set of different electronic ground states: this was much later (1998) shown to be unnecessary. In Ref. , in fact, a simpler viewpoint is taken: the polarization of a correlated solid is defined by means of a “pure state” expectation value, although a rather exotic kind of one. By the same token, it was also possible to define —and to compute —macroscopic polarization in noncrystalline systems.

In the present work we take advantage of the most recent developments for reconsidering the whole theory under a new light. At variance with previous presentations we are not going to introduce explicitly the Berry phase concept: in fact, within the formulation of Ref. , the Berry phase appears very much “in disguise”. Instead, a major role is played by a precursor work , apparently unrelated either to the polarization problem or to a Berry phase, which will be reexamined here and used to introduce the polarization theory.

Finally, let me just mention a latest development, not to be discussed in the present work, where ideas spawned from the polarization theory—and more specifically from Ref. —are used to investigate wavefunction localization .

## 2 The “electron in broth”

Adopting a given choice for the boundary conditions is tantamount to defining the Hilbert space where our solutions of Schrödinger’s equation live. For the sake of simplicity, I am presenting the basic concept by means of the one–dimensional case. For a single–particle wavefunction BvK reads $`\psi (x+L)=\psi (x)`$, where $`L`$ is the imposed periodicity, chosen to be large with respect to atomic dimensions. Notice that lattice periodicity is not assumed, and BvK applies to disordered systems as well.

By definition, an operator maps any vector of the given Hilbert space into another vector belonging to the same space: the multiplicative position operator $`x`$ is therefore not a legitimate operator when BvK are adopted for the state vectors, since $`x\psi (x)`$ is not a periodic function whenever $`\psi (x)`$ is such. It is then obvious why Eq. (1) cannot be used in condensed matter theory. Of course, any periodic function of $`x`$ is a legitimate multiplicative operator in the Hilbert space: this is the case e.g. of the nuclear potential acting on the electrons.

Before switching to the polarization problem, it is expedient to discuss an important precursor work, apparently unrelated to the polarization problem, where nonetheless the expectation value of the position operator plays the key role.

Some years ago, A. Selloni et al.  addressed the properties of electrons dissolved in molten salts at high dilution, in a paper which at the time was commonly nicknamed the “electron in broth”. The physical problem was studied by means of a mixed quantum–classical simulation, where a lone electron was adiabatically moving in a molten salt (the “broth”) at finite temperature. The simulation cell contained 32 cations, 31 anions, and a single electron. KCl was the original case study, which therefore addressed the liquid state analogue of an F center; other systems were studied afterwards . The motion of the ions was assumed as completely classical, and the Newton equations of motion were integrated by means of standard molecular dynamics (MD) techniques, though the ionic motion was coupled to the quantum degree of freedom of the electron. The electronic ground wavefunction was determined solving the time–dependent Schrödinger’s equation at each MD time step. As usual in MD simulations, periodic boundary conditions were adopted for the classical ionic motion. Ideally, the ionic motion occurs in a simulation cell which is surrounded by periodic replicas: inter-cell interactions are accounted for, thus avoiding surface effects. Analogously, the electronic wavefunction is chosen in the work of Selloni et al. to obey BvK over the simulation cell, and therefore features periodic replicas as well. A plot of such an electronic distribution, in a schematic one–dimensional analogue, is given in Fig. 1.

One of the main properties investigated in Ref.  was the electronic diffusion, where the thermal ionic motion is the driving agent (within the adiabatic approximation). In order to perform this study, one has to identify first of all where the “center” of the electronic distribution is. Intuitively, the distribution in Fig. 1 appears to have a “center”, which however is defined only modulo the replica periodicity, and furthermore cannot be evaluated simply as in Eq. (1), i.e. $`x=𝑑xx|\psi (x)|^2`$, precisely because of BvK. Selloni et al. solved the problem by means of a very elegant and far–reaching formula, presented below. The work of Ref.  can be regarded as the many–body generalization of it.

## 3 The main formula: One electron

According to Refs. , the key quantity for dealing with the position operator within BvK is the dimensionless complex number $`𝔷`$, defined as:

$$𝔷=\psi |\mathrm{e}^{i\frac{2\pi }{L}x}|\psi =_0^L𝑑x\mathrm{e}^{i\frac{2\pi }{L}x}|\psi (x)|^2,$$
(2)   
whose modulus is no larger than 1. The most general electron density, such as the one depicted in Fig. 1, can always be written as a superposition of a function $`n_{\mathrm{loc}}(x)`$, normalized over $`(\mathrm{},\mathrm{})`$, and of its periodic replicas:

$$|\psi (x)|^2=\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}n_{\mathrm{loc}}(xx_0mL).$$
(3)   
Both $`x_0`$ and $`n_{\mathrm{loc}}(x)`$ have a large arbitrariness: we restrict it a little bit by imposing that $`x_0`$ is the center of the distribution, in the sense that $`_{\mathrm{}}^{\mathrm{}}𝑑xxn_{\mathrm{loc}}(x)=0`$. Using Eq. (3), $`𝔷`$ can be expressed in terms of the Fourier transform of $`n_{\mathrm{loc}}`$ as:

$$𝔷=\mathrm{e}^{i\frac{2\pi }{L}x_0}\stackrel{~}{n}_{\mathrm{loc}}(\frac{2\pi }{L}).$$
(4)   
If the electron is localized in a region of space much smaller than $`L`$, its Fourier transform is smooth over reciprocal distances of the order of $`L^1`$ and can be expanded as:

$$\stackrel{~}{n}_{\mathrm{loc}}(\frac{2\pi }{L})=1\frac{1}{2}\left(\frac{2\pi }{L}\right)^2_{\mathrm{}}^{\mathrm{}}𝑑xx^2n_{\mathrm{loc}}(x)+𝒪(L^3).$$
(5)   
A very natural definition of the center of a localized periodic distribution $`|\psi (x)|^2`$ is therefore provided by the phase of $`𝔷`$ as:

$$x=\frac{L}{2\pi }\text{Im ln}𝔷,$$
(6)   
which is in fact the formula first proposed by Selloni et al. . The expectation value $`x`$ is defined modulo $`L`$, as expected since $`|\psi (x)|^2`$ is BvK periodic.

The above expressions imply $`xx_0\text{mod }L`$; in the special case where $`n_{\mathrm{loc}}(x)`$ can be taken as an even (centrosymmetric) function, its Fourier transform is real and Eq. (4) yields indeed $`xx_0\text{mod }L`$. In the case of extreme delocalization we have instead $`|\psi (x)|^2=1/L`$ and $`𝔷=0`$: hence the center of the distribution $`x`$, according to Eq. (6), is ill–defined. For a more general delocalized state, we expect that $`𝔷`$ goes to zero at large $`L`$ .

We have therefore arrived at a definition of $`x`$ within BvK which has many of the desirable features we were looking for: nonetheless, there is a property that is even more important, and which we are going to demonstrate now. Suppose the potential which the electron moves in has a slow time dependence—as was the case in Ref. —and we wish to follow the adiabatic evolution of the electronic state $`|\psi `$. If we call $`|\phi _j`$ the instantaneous eigenstates at time $`t`$, the lowest order adiabatic evolution of the ground–state density matrix is :

$$|\psi \psi ||\phi _0\phi _0|+i\underset{j0}{}\left(|\phi _j\frac{\phi _j|\dot{\phi }_0}{ϵ_jϵ_0}\phi _0|\text{Hc}\right),$$
(7)   
where the phases have been chosen in order to make $`|\phi _0`$ orthogonal to its time derivative $`|\dot{\phi }_0`$. The macroscopic electrical current flowing through the system at time $`t`$ is therefore:

$$j=\frac{1}{L}\psi |p|\psi \frac{i}{L}\underset{j0}{}\frac{\phi _0|p|\phi _j\phi _j|\dot{\phi }_0}{ϵ_jϵ_0}+\text{cc}.$$
(8)   
It is then rather straightforward to prove that $`j`$ to lowest order in $`1/L`$ equals $`(1/L)dx/dt`$, where $`x`$ is evaluated using in Eq. (6) the instantaneous ground eigenstate:

$$j\frac{1}{2\pi }\text{Im }\frac{d}{dt}\text{ln }\phi _0|\mathrm{e}^{i\frac{2\pi }{L}x}|\phi _0.$$
(9)   
This finding proves the value of the “electron–in–broth” formula, Eqs. (2) and (6) in studying electron transport .

## 4 The main formula: Many electrons

So much about the one–electron problem: we are now going to consider a finite density of electrons in the periodic box. To start with, irrelevant spin variables will be neglected: for the sake of notation simplicity, I will first illustrate the main concepts on a system of “spinless electrons”.

Even for a system of independent electrons, our approach takes a simple and compact form if a many–body formulation is adopted. BvK then imposes periodicity in each electronic variable separately:

$$\mathrm{\Psi }_0(x_1,\mathrm{},x_i,\mathrm{},x_N)=\mathrm{\Psi }_0(x_1,\mathrm{},x_i+L,\mathrm{},x_N).$$
(10)   
Our interest is indeed in studying a bulk system: $`N`$ electrons in a segment of length $`L`$, where eventually the thermodynamic limit is taken: $`L\mathrm{}`$, $`N\mathrm{}`$, and $`N/L=n_0`$ constant. We also assume the ground state nondegenerate, and we deal with insulating systems only: this means that the gap between the ground eigenvalue and the excited ones remains finite for $`L\mathrm{}`$.

We start defining the one–dimensional analogue of $`\widehat{𝐑}`$, namely, the multiplicative operator $`\widehat{X}=_{i=1}^Nx_i`$, and the complex number

$$𝔷_N=\mathrm{\Psi }|\mathrm{e}^{i\frac{2\pi }{L}\widehat{X}}|\mathrm{\Psi }.$$
(11)   
It is obvious that the operator $`\widehat{X}`$ is ill–defined in our Hilbert space, while its complex exponential appearing in Eq. (11) is well defined. The main result of Ref.  is that the ground–state expectation value of the position operator is given by the analogue of Eq. (6), namely:

$$X=\frac{L}{2\pi }\text{Im ln }𝔷_N,$$
(12)   
a quantity defined modulo $`L`$ as above.

The right–hand side of Eq. (12) is not simply the expectation value of an operator: the given form, as the imaginary part of a logarithm, is indeed essential. Furthermore, its main ingredient is the expectation value of the multiplicative operator $`\mathrm{e}^{i\frac{2\pi }{L}\widehat{X}}`$: it is important to realize that this is a genuine many–body operator. In general, one defines an operator to be one–body whenever it is the sum of $`N`$ identical operators, acting on each electronic coordinate separately: for instance, the $`\widehat{X}`$ operator is such. In order to express the expectation value of a one–body operator the full many–body wavefunction is not needed: knowledge of the one–body reduced density matrix $`\rho `$ is enough: I stress that, instead, the expectation value of $`\mathrm{e}^{i\frac{2\pi }{L}\widehat{X}}`$ over a correlated wavefunction cannot be expressed in terms of $`\rho `$, and knowledge of the $`N`$-electron wavefunction is explicitly needed. In the special case of a single–determinant, the $`N`$-particle wavefunction is uniquely determined by the one–body reduced density matrix $`\rho `$ (which is the projector over the set of the occupied single–particle orbitals): therefore the expectation value $`X`$, Eq. (12), is uniquely determined by $`\rho `$. But this is peculiar to uncorrelated wavefunctions only: this case is discussed in detail below.

As in the one–body case, whenever the many–body Hamiltonian is slowly varying in time, the macroscopic electrical current flowing through the system is given by

$$J=\frac{1}{L}\frac{d}{dt}X,$$
(13)   
where $`X`$ is evaluated using in Eq. (11) the instantaneous ground eigenstate of the Hamiltonian at time $`t`$: this result is proved in Ref. . Considering now the limit of a large system, $`X`$ is an extensive quantity: the macroscopic current $`J`$, Eq. (13), goes therefore to a well defined thermodynamic limit.

We stress that nowhere in our presentation have we assumed crystalline periodicity. Therefore our definition of $`X`$ is very general: it applies to any condensed system, either ordered or disordered, either independent–electron or correlated.

## 5 Macroscopic polarization

In the Introduction, we have discussed what polarization is not, by outlining some incorrect definitions  and their problems . We have not stated yet what polarization really is: to this aim, a few experimental facts are worth recalling. The absolute polarization of a crystal in a given state has never been measured as a bulk property, independent of sample termination. Instead, well known bulk properties are derivatives of the polarization with respect to suitable perturbations: permittivity, pyroelectricity, piezoelectricity, dynamical charges, In one important case—namely, ferroelectricity—the relevant bulk property is inferred from the measurement of a finite difference (polarization reversal). In all cases, the derivative or the difference in the polarization is typically accessed via the measurement of a macroscopic current. For instance, to measure the piezoelectric effect, the sample is typically strained along the piezoelectric axis while being shorted out with a capacitor (see Fig. 2).

The theory discussed here only concerns phenomena where the macroscopic polarization is induced by a source other than an electric field. Even in this case, the polarization may (or may not) be accompanied by a field, depending on the boundary conditions chosen for the macroscopic sample. The theory addresses polarization differences in zero field: this concerns therefore lattice dynamics, piezoelectricity (as in the ideal experiment sketched in Fig. 2), and ferroelectricity. Notably, the theory reported here does not address the problem of evaluating the dielectric constant: this can be done using alternative approaches, such as the well–established linear–response theory , or other more innovative theories .

The bulk quantity of interest, to be compared with experimental measurements, is the polarization difference between two states of the given solid, connected by an adiabiatic transformation of the Hamiltonian. The electronic term in this difference is:

$$\mathrm{\Delta }P=_0^{\mathrm{\Delta }t}𝑑tJ(t),$$
(14)   
where $`J(t)`$ is the current flowing through the sample while the potential is adiabatically varied, i.e. precisely the quantity discussed in the previous Section, Eq. (13). Notice that in the adiabatic limit $`\mathrm{\Delta }t`$ goes to infinity and $`J(t)`$ goes to zero, while Eq. (14) yields a finite value, which only depends on the initial and final states. We may therefore write:

$$P=\underset{L\mathrm{}}{lim}\frac{1}{2\pi }\text{Im ln }\mathrm{\Psi }|\mathrm{e}^{i\frac{2\pi }{L}\widehat{X}}|\mathrm{\Psi },$$
(15)   
where it is understood that Eq. (15) is to be used twice, with the final and with the initial ground states, in order to evaluate the quantity of interest $`\mathrm{\Delta }P`$. Notice that $`L\mathrm{}`$ in Eq. (15) is a rather unconventional limit, since the exponential operator goes formally to the identity, but the size of the system and the number of electrons in the wavefunction increase with $`L`$.

## 6 The case of independent electrons

We now specialize to an uncorrelated system of independent electrons, whose $`N`$-electron wavefunction $`|\mathrm{\Psi }`$ is a Slater determinant. As discussed above, in this case the expectation value $`X`$, Eq. (12), is uniquely determined by the one–body density matrix. However, the formulation is simpler when expressing $`X`$ and the resulting polarization $`P`$ directly in terms of the orbitals.

We restore explicit spin variables from now on. Suppose $`N`$ is even, and $`|\mathrm{\Psi }`$ is a singlet. The Slater determinant has thus the form:

$$|\mathrm{\Psi }=\frac{1}{\sqrt{N!}}|\phi _1\overline{\phi }_1\phi _2\overline{\phi }_2\mathrm{}\phi _{N/2}\overline{\phi }_{N/2}|,$$
(16)   
where $`\phi _i`$ are the single-particle orbitals. It is then expedient to define

$$|\stackrel{~}{\mathrm{\Psi }}=\mathrm{e}^{i\frac{2\pi }{L}\widehat{X}}|\mathrm{\Psi }:$$
(17)   
even $`|\stackrel{~}{\mathrm{\Psi }}`$ is indeed a Slater determinant, where each orbital $`\phi _i(x)`$ of $`|\mathrm{\Psi }`$ is multiplied by the plane wave $`\mathrm{e}^{i\frac{2\pi }{L}x}`$. According to a well known theorem, the overlap amongst two determinants is equal to the determinant of the overlap matrix amongst the orbitals. We therefore define the matrix (of size $`N/2\times N/2`$):

$$S_{ij}=\phi _i|\mathrm{e}^{i\frac{2\pi }{L}x}|\phi _j=_0^L𝑑x\phi _i^{}(x)\mathrm{e}^{i\frac{2\pi }{L}x}\phi _j(x),$$
(18)   
in terms of which we easily get

$$P=\frac{1}{2\pi }\text{Im ln }\mathrm{\Psi }|\stackrel{~}{\mathrm{\Psi }}=\frac{1}{\pi }\text{Im ln det }S,$$
(19)   
where the factor of 2 accounts for double spin occupancy, and the expression becomes accurate in the limit of a large system.

The expression of Eq. (19) goes under the name of “single–point Berry phase” (almost an oxymoron!), and was first proposed by the present author in a volume of lecture notes . Since then, its three-dimensional generalization has been used in a series of DFT calculations for noncrystalline systems , and has been scrutinized in some detail in Ref. .

The case of a crystalline system of independent electrons is the one which historically has been solved first , though along a very different logical path  than adopted here. I am going to outline how the present formalism leads to the earlier results.

For the sake of simplifying notations, I am going to consider the case of an insulator having only one completely occupied band. The single–particle orbitals may be chosen in the Bloch form: indeed, the canonical ones must have the Bloch form. The many–body wavefunction, Eq. (16), becomes in the crystalline case:

$$|\mathrm{\Psi }=\frac{1}{\sqrt{N!}}|\psi _{q_1}\overline{\psi }_{q_1}\psi _{q_2}\overline{\psi }_{q_2}\mathrm{}\psi _{q_{N/2}}\overline{\psi }_{q_{N/2}}|.$$
(20)   
The lattice constant is $`a=2L/N`$, and the Bloch vectors entering Eq. (20) are equally spaced in the reciprocal cell $`(0,2\pi /a]`$:

$$q_s=\frac{4\pi }{Na}s,s=1,2,\mathrm{},N/2.$$
(21)   
Owing to the orthogonality properties of the Bloch functions, the overlap matrix elements in Eq. (18) vanish except when $`q_s^{}=q_s2\pi /L`$, that is $`s^{}=s1`$: therefore the determinant in Eq. (19) factors as

$$P=\frac{1}{\pi }\text{Im ln }\underset{s=1}{\overset{N/2}{}}\psi _{q_s}|\mathrm{e}^{i\frac{2\pi }{L}x}|\psi _{q_{s1}},$$
(22)   
where $`\psi _{q_0}(x)\psi _{q_{N/2}}(x)`$ is implicitly understood (so–called periodic gauge).

The expression in Eq. (22) is precisely the one first proposed by King–Smith and Vanderbilt  as a discretized form for the Berry phase: in fact, the multi–band three–dimensional generalization of Eq. (22) is the standard formula  implemented in first–principle calculations  of macroscopic polarization in crystalline dielectrics, either within DFT or HF.

## 7 Critical rethinking of DFT

The modern viewpoint about macroscopic polarization has even spawned a critical rethinking of density–functional theory in extended systems. The debate started in 1995 with a paper by Gonze, Ghosez, and Godby , and continues these days . The treatment of polarization provided in the present work allows discussing the issue in a very simple way.

The celebrated Hohenberg–Kohn (HK) theorem, upon which DFT is founded , states that there exists a universal functional $`F[n]`$, which determines the exact ground–state energy and all other ground–state properties of the system. The main hypotheses are that the magnetic field is vanishing, the ground state is non degenerate, and—most important to the present purposes—the system is finite, with a square–integrable ground wavefunction. This is precisely the key point when dealing with an extended system.

Ideally, it is possible to refer to a macroscopic but finite system: polarization is then by definition the dipole divided by the volume, and Eq. (1) safely applies. As a consequence, the polarization of the real interacting system is identical to the one of the fictitious noninteracting Kohn–Sham (KS) system. Unfortunately, such polarization depends on the charge distribution both in the bulk and at the surface: according to Ref. , this fact implies a possible “ultranonlocality” in the KS potential: two systems having the same density in the bulk region may have qualitatively different KS potentials in the same region, and different polarizations as well. The issue can be formulated in a more transparent way by recasting it in the language of the present work.

As discussed in the Introduction, condensed matter theory invariably works in a different way: one adopts BvK since the very beginning, and the system has no surface by construction. The original HK theorem was formulated for the case where the Schrödinger’s equation is solved imposing square–integrable boundary conditions, but the same theorem holds within BvK, with an identical proof, for a finite $`N`$-electron system. We can then define, even within BvK, the fictitious KS system of noninteracting electrons, having the same density as the interacting one: if the system is crystalline, the KS orbitals have the Bloch form.

Now the question becomes: do the interacting system and the corresponding KS noninteracting one have the same polarization? The answer is actually “no”: in fact, within BvK, polarization is not a function of the density, not even of the one–body density matrix, as stressed above. Numerical evidence of the fact that the two polarizations are not equal has been given .

Other important implications concern the occurrence of macroscopic electric fields within DFT: we refer to the original literature  about this issue, while here we limit ourselves to just remarking an important point. The potential (both one–body and two–body) within the Schrödinger’s equation must be BvK periodic, otherwise the Hamiltonian is an ill–defined operator in the Hilbert space. The periodicity of the potential is tantamount to enforcing a vanishing macroscopic electric field: therefore some ad–hoc strategies must be devised in order to cope with nonzero electric fields, as is indeed done in Refs. .

## Acknowledgments

Invaluable discussions with R.M. Martin, G. Ortíz, and D. Vanderbilt are gratefully acknowledged. Part of this work has been performed while attending the 1998 workshop “Physics of Insulators” at the Aspen Center for Physics. Partially supported by ONR grant N00014-96-1-0689.