7. Modern Theory of Polarization

https://doi.org/10.48550/arXiv.2403.13981

7.1 Introduction

The Modern Theory of Polarization (MTOP), which was developed in the 1990s by Resta and Vanderbilt and their collaborators, solved the problem of how to define the polarization current ${\Jconv}$ of an insulating material in terms of the evolving microstructure of its bulk, and provided a widely-applicable method of calculating it [Resta, 1993; King-Smith and Vanderbilt, 1993; Vanderbilt and King-Smith, 1993; Resta, 1994; Resta and Vanderbilt, 2007].

I emphasize that the MTOP defines $\Jconv$ in terms of the changing bulk microstructure because the alternative, namely defining it as ${\Jconv\equiv\bsigmadot}$ where ${\bsigma}$ is the surface charge density, would be of much less practical use: Calculations of $\Jconv$ would require knowledge of the time-dependent surface microstructure. Since the numbers of particles in surface regions tend to be very large, accurate calculations of $\Jconv$ would be intractable in many or most cases. However, the MTOP can be used to calculate the electronic contribution to ${\Jconv}$ from a simulation that only takes explicit account of the electrons in a small region whose microstructure is representative of the entire bulk. This makes it possible to calculate $\Jconv$ in a wide range of materials.

In subsections Sec. 7.3 and Sec. 7.5 I present illustrative derivations of the MTOP expressions for the electronic contribution to $\Jconv$; and in Sec. 13 I present an explanation of the MTOP definition of $\Jconv$ that is based on the relation ${\J=\bsigmadot}$ and on Finnis’s expression for the (macroscale) excess of charge on a surface, ${\bsigma[\rhop,\rhom]}$, in terms of the microscopic densities of positive ($\rhop$) and negative ($\rhom$) charge. These derivations make it clear that the MTOP definition of the electronic polarization current density ${\Jconvm}$ in an insulator whose electrons are in their ground state is basically exact.

7.2 How is the polarization current density $\Jconv$ defined?

Suppose that the microstructure of a material is in a stationary or quasi-stationary statistical state and that $\zeta$ is one or more parameters on which the stationary state has a continuous and continuously-differentiable dependence. Then, as discussed in Sec. 6.8, the polarization current density,

\begin{align*} \Jconv = \Jconvp+\Jconvm, \end{align*}

is the current that flows while the stationary-state charge density,

\begin{align*} \rho(\zeta)=\rhop(\zeta)+\rhom(\zeta), \end{align*}

is changing continuously in response to $\zeta$ changing quasistatically; where quasistatically means so slowly that deviations of the microstructure’s statistical state from its $\zeta$-dependent stationary state are negligible; and ${\Jconvp}$ and ${\Jconvm}$ are the contributions to $\Jconv$ from the charge densities of the nuclei (${\rhop}$) and electrons (${\rhom}$), respectively, changing.

A changing microscopic charge density implies the existence of a microscopic current density,

\begin{align*} \jconv(\rvec,t)=\jconvp(\rvec,t)+\jconvm(\rvec,t), \end{align*}

which satisfies the charge conservation equation,

\begin{align*} \div\jconv &= -\pdv{\rho}{t} = -\dot{\zeta}\pdv{\rho}{\zeta} \\ &= -\dot{\zeta} \pdv{\rhop}{\zeta} -\dot{\zeta} \pdv{\rhom}{\zeta} = \div\jconvp+\div\jconvm. \end{align*}

The macroscopic polarization current densities $\Jconv$, $\Jconvp$, and $\Jconvm$ are the mesoscale spatial averages of $\jconv$, $\jconvp$, and $\jconvm$, respectively, in the material’s bulk.

7.2.1 Contribution of the nuclei to $\Jconv$

Let us assume that ${\rho}$ has the form of Eq. (3). Therefore the nuclei can be treated as point charges and, if $V$ is a region of the bulk whose microstructure is representative of the bulk, their contribution to the polarization current density can be calculated from the classical expression,

\begin{align*} \Jconvp = \frac{e}{\abs{V}}\sum_{\Rvec_i\in V}Z_i\dv{\Rvec_i}{t} = \frac{e \dot{\zeta}}{\abs{V}}\sum_{\Rvec_i\in V}Z_i\dv{\Rvec_i}{\zeta}, \end{align*}

where ${Z_i}$ and ${\Rvec_i}$ are the ${i^\text{th}}$ nucleus’s atomic number and position, respectively, and ${\abs{V}}$ is the volume of $V$.

7.2.2 Contribution of the electrons to ${\Jconv}$

This subsection states one form of the MTOP expression for the electronic polarization current density ($\Jconvm$) that flows while $\zeta$ is changing.

The electrons’ charge density at position $\rvec$ is

\begin{align*} \rhom(\rvec;\zeta)=-e n(\rvec;\zeta), \end{align*}

where ${n(\zeta)}$ is the number density of electrons, or simply the ‘electron density’. Let us assume that the electron density has a ${\zeta}$-dependent decomposition,

\begin{align*} n(\rvec;\zeta)= -\rhom(\rvec;\zeta)/e= \sum_i n_i(\rvec;\zeta), \end{align*}

where each density ‘packet’ ${n_i(\zeta)}$ evolves smoothly with $\zeta$, without its integral changing and therefore without the charge,

\begin{align*} q_i =-e\int n_i(\rvec;\zeta)\ddpow{3}{r}, \end{align*}

that it carries changing. In Sec. 7.4 I will outline how such a decomposition can be found in practice.

Assuming that it can be found, the electronic contribution to $\Jconv$ is

\begin{align*} \Jconvm = \frac{1}{\abs{V}}\sum_{\rvecsub{i}\in V} q_i \dv{\rvecsub{i}}{t}= \frac{\dot{\zeta}}{\abs{V}} \dvone{\zeta} \left(\sum_{\rvecsub{i}\in V} q_i \rvecsub{i}\right), \tag{10} \end{align*}

where the sum is over all ${\rvecsub{i}}$ within a region $V$ of the bulk that is representative of the entire bulk, and

\begin{align*} \rvecsub{i}(\zeta) \equiv \frac{\intthree \rvec \,n_i(\rvec;\zeta)\ddpow{3}{r}}{\intthree n_i(\rvec';\zeta)\ddpow{3}{r'}}, \end{align*}

is the center of packet ${n_i}$. Eq. (10) is the MTOP definition of $\Jconvm$.

In practice, $\Jconvm$ is almost always calculated after the ground state electron density ${n(\zeta)}$ has been calculated with methods based on the density functional theory (DFT) of Hohenberg, Kohn, and Sham [Hohenberg and Kohn, 1964; van Leeuwen, 2003; Kohn and Sham, 1965; Dreizler and Gross, 1990]. At the end of Sec. 7.1 I described the MTOP definition of $\Jconvm$ as ‘exact’ for an insulator in its ground state. I meant two things by this. The first is that if $V$ is exactly representative of the bulk, Eq. (10) is an exact expression for ${\Jconv}$. The second is that a natural byproduct of the most widely used DFT-based methods of calculating the ground state density of an insulator is a decomposition of ${n(\zeta)}$ into packets whose centers are in $V$ and whose integrals are independent of $\zeta$. Eq. (10) cannot be used without such a decomposition.

If the material is a crystal, $V$ can be a single primitive unit cell, $\Omega$. In an amorphous material $V$ must be large enough to sample all relevant features of the microstructure and, by expanding it to improve the sampling of the bulk, ${\Jconv}$ can be calculated to any desired precision.

**Figure 3.** Schematic of a polarized three dimensional material, which complements the schematic of the quasi-one dimensional material in Fig. 1. The material is charge neutral, and its bulk (shaded green) is macroscopically charge neutral. The material has a finite dipole moment, ${\mbd\hat{x}}$, and on the macroscale the surfaces with outward normals ${\hat{x}}$ and ${-\hat{x}}$ carry equal and opposite areal charge densities, $\bsigma$ and $-\bsigma$, respectively. See Sec. 7.3 for more details.

7.3 Illustrative derivation of Eq. (10)

The physical idea underpinning Eq. (10) is simple, and can be illustrated using the isolated nonconducting material depicted in Fig. 3, which we will assume to be a crystal for the purposes of this section.

The crystal’s dipole moment is $\mbd\hat{x}$ and the surfaces perpendicular to $\hat{x}$ carry equal and opposite excesses of charge, in the form of areal charge densities of ${\pm\bsigma}$. Between ${\xl}$ and ${\xr}$ the crystal’s shape and cross-sectional area ${\abs{\plane}}$ are uniform: The region ${\plane(x)}$ of each plane perpendicular to $\hat{x}$ that the crystal occupies is the same for all ${x\in[\xl,\xr]}$.

Although the width of each surface is zero at the macroscale, each macroscopically-planar surface corresponds to a region of finite thickness $\prectheo$ at the microscale, half of which is in vacuum and half of which is within the material. In Fig. 3 the two surface regions are coloured blue and red/pink, with lighter shades indicating the parts of them that are in vacuum. The parts that are within the material are indicated by both darker shades and the labels ${\surfaceL}$ and ${\surfaceR}$. The widths of the surface regions are exaggerated. In reality $\prectheo$ would be many orders of magnitude smaller than $\size$; therefore,

\begin{align*} \bulksize\equiv\abs{\bulk}=\size-\prectheo\approx\size. \end{align*}

To understand Eq. (10), let us choose the representative region of the bulk, $V$, to be an arbitrary unit cell in the crystal’s bulk, $\unitcell$, and let us derive the equation,

\begin{align*} \Jconvm\cdot\hat{x} = \left(\frac{1}{\volume}\sum_{\rvecsub{i}\in \unitcell} q_i \dv{\rvecsub{i}}{t}\right)\cdot\hat{x}, \tag{11} \end{align*}

which is the part of Eq. (10) that pertains to the $x$-component of $\Jconvm$. Note that, in this subsection, the unit cell ${\unitcell}$ is three dimensional and ${\volume}$ is its volume, but in Sec. 7.5, in later chapters, and in Appendices Appendix B, Appendix G, and Appendix H, each primitive unit cell $\unitcell$ is (quasi-) one dimensional and $\volume$ denotes its width.

Let us begin by assuming that the electron number density in the bulk can be expressed as a superposition of ${\NUnitcell}$ identical electron densities, whose integrals are independent of $\zeta$, such that each bulk primitive unit cell contains the center of exactly one of them. Let ${\nunitcell(\zeta)}$ denote the density packet whose center is in $\unitcell$, and which carries a total charge,

\begin{align*} q^-= -e\intthree \nunitcell(\rvec;\zeta)\ddpow{3}{r}. \end{align*}

Now let us choose an arbitrary point in unit cell $\unitcell$ to be the origin, which means that the density in the bulk is

\begin{align*} n(\rvec;\zeta) = \sum_{\Avec\in\bulk} \nunitcell(\rvec-\Avec;\zeta),\;\forall \rvec\in\bulk, \end{align*}

where the sum is over all distinct lattice vectors ${\Avec}$ such that the point displaced from the origin by $\Avec$ is in ${\bulk}$.

The charge neutrality of the crystal’s bulk implies that ${\rhocell(\zeta)}$, defined as the sum of ${\rhomcell(\zeta)=-e\,\nunitcell(\zeta)}$ and the delta charge distribution of the nuclei in the cell $\unitcell$, is a charge neutral distribution with a well defined (i.e., origin independent) dipole moment. The $x$-component of this dipole moment will be denoted by $\dipcell$.

When the sum of all $\Nunitcell$ bulk unit cells’ charge densities is subtracted from the charge density of the entire crystal, $\rho$, the result is

\begin{align*} \Drho(\rvec;\zeta)=\rho(\rvec;\zeta)- \sum_{\Avec\in\bulk} \rhocell(\rvec-\Avec;\zeta). \end{align*}

This vanishes everywhere in the bulk, but not everywhere in the surface regions, ${\surfaceL}$ and ${\surfaceR}$. Therefore if $\Drhobar(x;\zeta)$ denotes the average of ${\Drho(\rvec;\zeta)=\Drho(x,\svec;\zeta)}$ on the plane perpendicular to $\hat{x}$ at $x$, the areal charge density,

\begin{align*} \sigma_s(\zeta) \equiv -\int_{-\eta}^{\eta}\Drhobar(\xl+u;\zeta)\dd{u}= \int_{-\eta}^{\eta}\Drhobar(\xr+u;\zeta)\dd{u}, \end{align*}

is finite and independent of $\eta$ as long as ${\prectheo/2<\eta\ll\bulksize}$. The integrals of ${\Drhobar}$ across the left and right surfaces are equal in magnitude and opposite in sign because we have assumed that the material is charge-neutral overall. Note, however, that ${\bsigma}$ does not equal ${\sigma_s}$ unless ${\dipcell=0}$. This will be explained and discussed in Sec. 9 and Sec. 12.

The volume of the crystal can be expressed as

\begin{align*} \size\abs{\plane}=\left(\bulksize+\prectheo\right)\abs{\plane} =\NUnitcell\volume+\prectheo\abs{\plane}, \tag{12} \end{align*}

where the second term on the right hand side is much smaller than the first, because ${\NUnitcell\volume=\abs{\plane}\bulksize}$; which also means that

\begin{align*} \prectheo\abs{\plane}/(\NUnitcell\volume)=\prectheo/\bulksize \sim a/S < a/L \ll a/l \ll 1. \end{align*}

The $x$-component of the crystal’s dipole moment can be expressed as

\begin{align*} \mbd=\NUnitcell \dipcell + \sigma_s \abs{\plane}(\NUnitcell\volume/\abs{\plane}+\prectheo), \end{align*}

where the first term on the right hand side is the sum of the $x$-components of the dipole moments of the bulk unit cells; and the second term is the product of the net charge ${\sigma_s\abs{\plane}}$ in ${\surfaceR}$ and the distance between the surfaces. It follows that the average, over the entire crystal, of the $x$-component of the current density flowing within it is

\begin{align*} \Jconv\cdot\hat{x}& =\bsigmadot\cdot\hat{x} =\mbpdot\cdot\hat{x} =\frac{\mbddot}{\NUnitcell\volume+\prectheo\abs{\plane}} \\ & =\frac{\NUnitcell \dipcelldot + \dot{\sigma}_s (\NUnitcell\volume+\prectheo\abs{\plane})}{\NUnitcell\volume+\prectheo\abs{\plane}} \\ &= \frac{\dipcelldot}{\volume} + \dot{\sigma}_s + \order{a/S} = \frac{\dipcelldot}{\volume} + \order{a/S}, \tag{13} \end{align*}

where ${\mbp\equiv \mbd/(\NUnitcell\volume+\prectheo\abs{\plane})}$ is the $x$-component of the crystal’s dipole moment divided by the crystal’s volume; and $\dot{\sigma}_s$ vanishes because we have assumed that there is no conduction current and that the crystal is isolated.

The term of order ${a/S}$ on the right hand side of Eq. (13) is negligible. Therefore,

\begin{align*} \Jconv&\cdot\hat{x} = \frac{\dipcelldot}{\volume} \\ &=\frac{1}{\volume}\dvone{t}\left(-e \intthree \rvec \nunitcell(\rvec;\zeta)\ddpow{3}{r} + e\sum_{\Rvecsub{i}\in\Omega} Z_i\Rvecsub{i}\right)\cdot\hat{x} \\ &=\underbrace{\vphantom{\sum_{\Rvecsub{i}\in\Omega}}\frac{q^-}{\volume} \dv{\rvecsub{\unitcell}}{t}\cdot\hat{x}}_{\displaystyle \Jconvm\cdot\hat{x}} + \underbrace{\frac{e}{\volume} \sum_{\Rvecsub{i}\in\Omega} Z_i\dv{\Rvecsub{i}}{t}\cdot\hat{x}}_{\displaystyle \Jconvp\cdot\hat{x}}. \end{align*}

The electronic contribution to the polarization current is

\begin{align*} \Jconvm\cdot\hat{x} &= \frac{1}{\volume} q^- \dv{\rvecsub{\unitcell}}{t}, \tag{14} \end{align*}

which has the same form as Eq. (11) when only one density packet has its center in $\Omega$.

The derivation of Eq. (14) neglected contributions to ${\mbddot}$ from redistributions of charge within the surface regions ${\surfaceL}$ and ${\surfaceR}$, because changes of surface polarization should not be counted as contributing to the bulk polarization current, and because these surface contributions vanish in the limit ${\prectheo/\bulksize\to 0}$.

7.3.1 Merging and subdividing density packets

Now let us suppose that density packet ${\nunitcell(\zeta)}$ can be expressed as

\begin{align*} \nunitcell(\rvec;\zeta) = \sum_{i} n_i(\rvec;\zeta), \end{align*}

where packet ${n_i(\zeta)}$ is centered at ${\rvecsub{i}(\zeta)}$ and carries a charge,

\begin{align*} q_i^- = -e\intthree n_i(\rvec;\zeta)\ddpow{3}{r}, \end{align*}

that is independent of $\zeta$. Then,

\begin{align*} q^- \rvecsub{0}(\zeta) &= -e\intthree \rvec\,\nunitcell(\rvec;\zeta)\ddpow{3}{r} \\ &= -e\sum_i \intthree \rvec\,n_i(\rvec;\zeta)\ddpow{3}{r} = \sum_i q^-_i \,\rvecsub{i}(\zeta). \end{align*}

Therefore Eq. (14) can be expressed as

\begin{align*} \Jconvm\cdot\hat{x} &= \left(\frac{1}{\volume}\sum_i q_i^-\dv{\rvecsub{i}}{t}\right)\cdot\hat{x}. \tag{15} \end{align*}

This demonstrates that decomposing density packets into smaller packets does not change ${\Jconv}$ as long as the smaller packets’ integrals are independent of ${\zeta}$. It can also be shown that ${\Jconv}$ is invariant under mergers of multiple representative regions $V$ into larger representative regions, and under mergers of multiple $\zeta$-independent density packets into larger density packets.

Note that Eq. (15) is identical to Eq. (11) if the centers of the packets ${n_i(\zeta)}$ are all in cell $\unitcell$. However if, for example, the center ${\rvecsub{j}}$ of ${n_j(\zeta)}$ is not in $\unitcell$, there exists an image ${\nunitcell(\rvec-\Avec;\zeta)}$ of ${\nunitcell(\rvec;\zeta)}$, under a translation by a lattice vector $\Avec$, whose counterpart ${n_j(\rvec-\Avec;\zeta)}$ of ${n_j(\rvec;\zeta)}$ has its center in $\unitcell$. Therefore, Eq. (11) could be derived by repeating the derivation that led to Eq. (14) with the definition of ${\rhomcell(\zeta)}$ changed to

\begin{align*} \rhomcell(\rvec;\zeta) \equiv -e\sum_{\Avec\in\bulk}\sum_{\Avec+\rvecsub{i}\in\unitcell} n_i(\rvec-\Avec;\zeta), \end{align*}

i.e., by defining it as the sum of all of the smaller charge density packets whose centers are in $\unitcell$.

7.4 Decomposing particle number densities into `packets'

It should be clear from the derivation above that the partitioning of the electron density into packets (e.g., $\nunitcell$) does not constitute an approximation; and that, if the density can be decomposed into ${\zeta}$-dependent packets with $\zeta$-independent integrals, and if region $V$ represents the bulk microstructure exactly, Eq. (10) is an exact expression for $\Jconv$.

7.4.1 Applying the MTOP to electrons in their ground state

The derivation above does not assume that the material is an insulator or that the electrons are in their ground state, but both of these assumptions are made when the MTOP is used to calculate $\Jconvm$. It is also assumed that the electrons’ external potential ${\vext(\zeta)}$ is changing smoothly and quasistatically with $\zeta$. The external potential ${\vext(\rvec;\zeta)\in\realone}$ is the energy of interaction of an electron at $\rvec$ with the quasistatic nuclei and any externally-applied quasistatic fields [Kaxiras and Joannopoulos, 2019; van Leeuwen, 2003].

In most applications of the MTOP, the ground state density is calculated using the density functional theory (DFT) of Hohenberg, Kohn, and Sham [Hohenberg and Kohn, 1964; van Leeuwen, 2003; Kohn and Sham, 1965; Dreizler and Gross, 1990]. In most DFT calculations the ground state density is calculated as

\begin{align*} n(\rvec;\zeta) = \sum_{i=1}^\infty f_i\abs{\psi_i(\rvec;\zeta)}^2, \tag{16} \end{align*}

where the infinite orthonormal set,

\begin{align*} \densityset\big(\zeta;\hamfrakr\big) \equiv \bigg\{\psi_i(\zeta)\in&\lebesgue(\realone^3): \;1\leq i\leq \infty, \\ \hamfrakr\psi_i=\epsilon_i\psi_i,\;&\int_{\realone^3}\psi^*_i(\rvec;\zeta)\psi_j(\rvec;\zeta)\ddpow{3}{r}=\delta_{ij}\bigg\}, \end{align*}

is the set of eigenfunctions of the Hamiltonian,

\begin{align*} \hamfrakr(\zeta):\lebesgue(\realone^3)\to\lebesgue(\realone^3), \end{align*}

of a fictitious system of non-interacting electrons that has the same ground state density ${n(\zeta)}$ as the real system of mutually-repulsive electrons; and the eigenfunctions’ ground state ‘occupation numbers’, $f_i$, satisfy,

\begin{align*} \sum_{i=1}^\infty f_i = \Nelec, \end{align*}

where ${\Nelec}$ is the total number of electrons.

Let us assume that the eigenfunctions $\psi_i$ are indexed in order of increasing eigenvalue, $\epsilon_i$, i.e., such that,

\begin{align*} i\leq j \iff \epsilon_i\leq\epsilon_j. \end{align*}

Then, in an insulator, Eq. (16) simplifies to

\begin{align*} n(\rvec;\zeta) = 2\sum_{i=1}^{\Nelec/2} \abs{\psi_i(\rvec;\zeta)}^2 = \sum_{i=1}^{\Nelec/2} n_i(\rvec;\zeta), \tag{17} \end{align*}

where it is assumed that each eigenstate is either empty (${f_i=0}$) or occupied by two electrons (${f_i=2}$) with opposite spins, and that occupation numbers do not change as $\zeta$ changes. Therefore the integral of each density ‘packet’ $n_i(\zeta)$ is two.

The assumption that $\zeta$ and ${\vext(\zeta)}$ change quasistatically and the assumption that the material has a finite band gap, ${\bandgap>0}$, are necessary for the assumption that electrons are in their ground state while the polarization current is flowing to be valid. The assumption that the material is an insulator, rather than a semiconductor with a very small band gap, ${0<\bandgap\approx 0}$, is necessary for the assumption that electrons are in their ground state to be approximately true at all temperatures that are low relative to ${\bandgap/k_B}$ (e.g., ${1\,\text{eV}/k_B\approx 11605\,\text{K}}$); and a large gap ensures that the charge of ${-2e}$ carried by each smoothly-evolving packet of density, ${\abs{\psi_i(\zeta)}^2}$, does not change as its center moves.

It is known that any insulating material’s ground state electron density ${n(\rvec;\zeta)}$ can be approximated arbitrarily closely by Eq. (17) for some $1$-electron Hamiltonian of the form,

\begin{align*} \hamfrakr(\rvec;\zeta) = \hat{t} + \vext(\rvec;\zeta) + \Delta v(\rvec;\zeta), \end{align*}

where ${\hat{t}}$ is the $1$-electron kinetic energy operator; and ${\Delta v(\zeta)}$ adjusts ${\vext(\zeta)}$ to account exactly, or arbitrarily closely to exactly, for treating the electrons as if they do not interact with one another [Kohn and Sham, 1965; van Leeuwen, 2003; Dreizler and Gross, 1990]. A ground state electron density that is also the ground state density of a fictitious set of noninteracting electrons, which are confined by a potential ${v=\vext+\Delta v}$, it is said to be noninteracting $v$-representable.

Not only is it known that in principle every ground state electron density is arbitrarily close to a noninteracting $v$-representable density, it is known how to approximate ${\hamfrakr(\rvec;\zeta)}$ in practice [Hohenberg and Kohn, 1964; van Leeuwen, 2003; Kohn and Sham, 1965; Dreizler and Gross, 1990]. Therefore a set ${\{n_i(\zeta)\}}$ of $\zeta$-dependent density packets with $\zeta$-independent integrals can be calculated using one of many widely available and widely used software packages, such as Quantum ESPRESSO [Giannozzi et al., 2017] or ABINIT [Gonze et al., 2020], both of which calculate ${\Jconvm}$ if fed the appropriate keywords and parameters. These methods and codes approximate the ground state density, rather than calculating it exactly, but the MTOP method of calculating $\Jconv$ from set ${\{n_i(\zeta)\}}$ does not, in principle, introduce further approximations: Eq. (10) is an exact expression if region $V$ represents the bulk microstructure exactly.

Most calculations of set ${\{n_i(\zeta)\}}$ use Born-von Kármán boundary conditions, which will be discussed in Sec. 7.5, and which does constitute an approximation, albeit one that seems likely to be accurate.

7.4.2 Decomposing an arbitrary stationary state number density

Although the integral of each density packet tends to be one or two in standard applications of the MTOP to electrons, it is clear from the derivation of Eq. (11) above, and it will be clear from the derivation of Eq. (10) in Sec. 12.6 and Sec. 13, that any method of partitioning a stationary state number density ${n(\zeta)}$ into $\zeta$-dependent packets of fixed charge would give the same result, regardless of whether or not the packets’ charges are integer multiples of ${e}$.

Furthermore, neither the derivation above nor the derivation in Sec. 12.6 and Sec. 13, which is based on the homogenization theory of Sec. 8 and Sec. 9, require the particles whose number density is ${n(\zeta)}$ to be electrons, or require ${n(\zeta)}$ to be a ground state number density. There is only an implicit assumption that ${n(\zeta)}$ is the number density of a set of particles that are in a stationary statistical state.

To my knowledge, it is not currently known how to decompose an arbitrary stationary-state number density into $\zeta$-dependent packets with $\zeta$-dependent integrals; and it is not even known which densities admit such a decomposition in principle. However the literature on non-interacting $v$-representability of electron number densities does clarify two important issues, albeit not necessarily explicitly [van Leeuwen, 2003; Engel and Dreizler, 2011; Corso, 2025; Sutter et al., 2024; Penz et al., 2023]: First, there exists a broad range of particle number densities for which such a decomposition is possible in principle; and second, it is not necessary for a number density to be a property of a quantum mechanical stationary state for it to admit such a decomposition.

Therefore, if the number densities of other kinds of charged particles in other kinds of physical systems can be decomposed in practice, Eq. (10) is more generally applicable than its developers claimed it to be in at least some of their presentations of it [Resta, 1993; King-Smith and Vanderbilt, 1993; Vanderbilt and King-Smith, 1993; Resta, 1994; Resta and Vanderbilt, 2007].

7.4.3 Transforming between density decompositions

Here, and in the remainder of this work, particles’ spins, and the spin-degeneracies of eigenstates, will be ignored. It will be assumed that an insulator in its electronic ground state has ${\Nelec}$ singly-occupied eigenstates ${\ket{\psi_i(\zeta)}}$ of the operator,

\begin{align*} \hamfrak(\zeta)\equiv \intthree\ddpow{3}{r} \hamfrakr(\rvec;\zeta)\dyad{\rvec}, \end{align*}

which span a ${\Nelec-}$dimensional Hilbert space $\hilbert_\zeta$ that changes smoothly with $\zeta$. This means that, as $\zeta$ changes, ${\hilbert_\zeta}$ rotates within the infinite-dimensional Hilbert-Lebesgue space, ${\lebesgue(\realone^{\dimension \Nelec})}$, of which it is a subspace, where ${\dimension}$ is the dimensionality of the system. We will return to this concept in Sec. 13.2, and it is illustrated schematically for ${\hilbert_\zeta\equiv\Span\left\{\ket{\varphi_1},\ket{\varphi_2},\ket{\varphi_3}\right\}}$ in Fig. 13.

The charge density packet ${-e\abs{\psi_i(\rvec;\zeta)}^2}$ of each occupied eigenfunction ${\psi_i(\rvec;\zeta)\equiv\braket{\rvec}{\psi_i(\zeta)}}$ of ${\hamfrakr(\rvec;\zeta)}$ carries a charge of ${-e}$ and contributes to the polarization current via the motion of its center ${\rvecsub{i}(\zeta)}$ caused by ${\hilbert_\zeta}$’s rotation.

The remainder of Sec. 7 will focus on the case of a (quasi-)one dimensional material. Therefore ${\psi_i(\zeta)=\psi_i(x;\zeta)}$ is an eigenfunction of

\begin{align*} \hamfrakx:\lebesgue(\realone)\to\lebesgue(\realone), \end{align*}

and the center of $\psi_i$ will be denoted by $x_i$ instead of $\rvecsub{i}$. In Sec. 7.5, it will be assumed that the bulk of the material is in a torus, and ${\varphi_i(\zeta)}$ will denote the counterpart ${\psi_i(\zeta)}$ in that context. In Sec. 14, Appendix G, and Appendix H the counterpart of ${\hamfrakx(\zeta)}$ in $\onetorus$ will be denoted by ${\hamsmallx(\zeta)}$.

As will be discussed in Sec. 13 and in the caption of Fig. 13, it is the space ${\hilbert_\zeta}$ that determines ${n(\zeta)}$, not its basis of eigenstates. Therefore we are free to use a different basis to calculate ${\Jconvm}$; and it is common to transform the set ${\{\psi_i(\zeta)\}_{i=1}^{\Nelec}}$ to one for which the density packets ${\{\abs{\psi_i(\zeta)}^2\}_{i=1}^{\Nelec}}$ are highly localized. However, this step is unnecessary, in principle: the centers of the (delocalized) eigenstates can be used directly in the equation,

\begin{align*} \Jconvm = -\frac{e}{\abs{V}}\sum_{x_i\in V} \dot{x}_i = -\frac{e\dot{\zeta}}{\abs{V}}\sum_{x_i\in V}\dv{x_i}{\zeta} \tag{18} \end{align*}

which is the one dimensional version of Eq. (10) with ${q_i=-e}$, and where $V$ represents a one dimensional region of the bulk of width ${\abs{V}}$.

An alternative is to Fourier transform the eigenfunctions with respect to their position arguments, and to compute each $x_i$ as the expectation value of the Fourier transformed position operator ${\hat{x}^\ksuper\equiv i\!\partial/\partial k}$. For example, the center of ${n_i(x)=\abs{\psi_i(x)}^2}$, when ${\braket{\psi_i}{\psi_i}=1}$, is

\begin{align*} x_i &\equiv \int\dd{x} n_i(x) x = \int \dd{x} \psi_i^*(x)\left(x \psi_i(x)\right) \\ & = \fourierconst^2\!\int \dd{x} \int \dd{k} e^{ikx}\fourierspace{x \psi_i}(k) \int\dd{q} \fts{\psi}_i^*(k+q) e^{-i(k+q)x}\\ & = \fourierconst^2\!\int\dd{k} i\pdv{\fts{\psi}_i(k)}{k}\int\dd{q}\left(\int\dd{x} e^{-iqx}\right) \fts{\psi}_i^*(k+q) \\ & = \int\dd{k} \fts{\psi}_i^*(k)\left(i\pdvone{k}\right)\fts{\psi}_i(k) = \int\dd{k} \fts{\psi}_i^*(k)\,\hat{x}^\ksuper\,\fts{\psi}_i(k), \end{align*}

where ${\fourierconst\equiv 1/\sqrt{2\pi}}$. Despite my use of suggestive notation, there is no quantum mechanical content in these lines of mathematics. Any classical probability density ${n_i(x)}$ can be expressed as the square modulus of a Hilbert space vector,

\begin{align*} \psi_i(x)\equiv \sqrt{n_i(x)}e^{i\theta(x)}\in\lebesgue(\realone). \end{align*}

The phase ${\theta(x)}$ is often irrelevant and arbitrary when studying the stationary states of an isolated system, but it becomes important when dividing a system into subsystems (e.g., nuclei and electrons or surface and bulk) or when studying the time dependence of a nonequilibrium probability density.

7.5 The bulk subsystem

In theoretical solid state physics it is common to study a material’s bulk subsystem, meaning the material without its surfaces, by placing it in a torus. The bulk of a 1-d material would be represented in a 1-torus $\onetorus$, whose circumference is the width, ${\bulksize}$, of the material’s bulk. As discussed in Appendix G, representing bulk wavefunctions in $\onetorus$ is equivalent to enforcing $\bulksize$-periodic Born-von Kármán boundary conditions on wavefunctions in ${\realone}$ [Born and von Kármán, 1912].

In textbooks on solid state physics, Born-von Kármán boundary conditions are usually discussed in the context of crystalline materials [Cohen and Louie, 2016; Ashcroft and Mermin, 1976; Ibach and Lüth, 2012; Kittel, 2004], but in computational research they are also commonly used for amorphous materials [Trave et al., 2002; Pasquarello et al., 1992; Tangney and Scandolo, 2002]. The bulks of amorphous materials can be represented in $\onetorus$ by ensuring that $\bulksize$ is large enough to contain a representative sample of the material, where ‘representative’ means that the sample’s microstructure shares, or approximates, whatever statistical characteristics of the true material are relevant to the properties of it that are being studied.

The purpose of this subsection is to derive the MTOP expression for ${\Jconvm}$ in an insulator in terms of a set ${\varphiset}$ of smoothly-evolving functions, which are $\bulksize$-periodic and therefore respect the Born-von Kármán boundary conditions.

7.5.1 Calculating $\Jconvm$ in $\onetorus$

This subsection derives an expression for ${\Jconvm}$ in the limit ${\bulksize\to\infty}$ by assuming the existence of a set of functions,

\begin{align*} \varphiset\equiv \left\{\varphi_i(\zeta)\in\lebesgue(\onetorus): 1\leq i \leq \Nelecbulk \right\}, \end{align*}

whose elements vary smoothly with $\zeta$ while reproducing the bulk electron density ${n(\zeta)}$ via

\begin{align*} n(x;\zeta)=\sum_{i=1}^{\Nelecbulk} \abs{\varphi_i(x;\zeta)}^2 ,\;\;\forall x\in\onetorus; \end{align*}

and while satisfying the orthonormality condition,

\begin{align*} \braketT{\varphi_i(\zeta)}{\varphi_j(\zeta)}\equiv \int_\onetorus\varphi_i^*(x;\zeta)\varphi_j(x;\zeta)\dd{x}=\delta_{ij}, \tag{19} \end{align*}

for all ${i,j \in\{1\cdots\Nelecbulk\}}$, where ${\int_\onetorus \equiv \int_{x_0}^{x_0+\bulksize}}$ for any ${x_0\in\onetorus}$.

In a crystal the elements of ${\{\varphi_i(\zeta)\}}$ are usually Bloch functions, which are discussed in Appendix G. However, before we assume that ${n(\zeta)}$ is $\volume$-periodic, we will derive an expression for $\Jconvm$ in the more general case of a density ${n(\zeta)}$ that is ${\bulksize}$-periodic but not $\volume$-periodic.

The $\bulksize$-periodicity of $\onetorus$ allows each element of $\varphiset$ to be expressed as a Fourier series,

\begin{align*} \varphi_i(x;\zeta) \equiv \sum_{g\in \reciplattg} \ftsvarphi_i(g;\zeta) e^{igx}, \tag{20} \end{align*}

where $\reciplattg$ is the set of all wavevectors that are integer multiples of ${\hbulksize\equiv 2\pi/\bulksize}$ (see Appendix B) and

\begin{align*} \ftsvarphi(g;\zeta) \equiv \frac{1}{\bulksize}\int_\onetorus \varphi_i(x;\zeta) e^{-ig x}\dd{x}. \end{align*}

It follows from Eq. (19) that the normalization of ${\ftsvarphi_i}$ is

\begin{align*} \sum_g \abs{\ftsvarphi_i(g;\zeta)}^2=\frac{1}{\bulksize}, \end{align*}

where ${\sum_g}$ is an abbreviation of ${\sum_{g\in\reciplattg}}$.

Eq. (20) can be used to express the center of ${\varphi_i(\zeta)}$ as

\begin{align*} x_i(\zeta) &\equiv \int_\onetorus x \,\abs{\varphi_i(x;\zeta)}^2 \dd{x} \\ & = \sum_{g}\ftsvarphi^*_i(g;\zeta)\sum_{g'}\ftsvarphi_i(g+g';\zeta)\int_\onetorus x e^{ig'x}\dd{x}. \tag{21} \end{align*}

Now let us assume that ${x_i(\zeta)}$ is being calculated in the ${\bulksize\to \infty}$ limit, which is the limit in which $\hbulksize$, the smallest difference between elements of $\reciplattg$, vanishes. In this limit,

\begin{align*} \int_\onetorus x e^{i g' x}\dd{x} &= -\frac{i}{\hbulksize}\int_\onetorus \left(e^{i(g'+\hbulksize)x}-e^{ig'x}\right)\dd{x}. \tag{22} \end{align*}

The integrals ${\int_\onetorus e^{i(g'+\hbulksize)x}\dd{x}}$ and ${\int_\onetorus e^{ig'x}\dd{x}}$ vanish for all values of ${g'}$ in $\reciplattg$ except ${g'=-\hbulksize}$ and ${g'=0}$, respectively, in which cases their values are ${\int_\onetorus\dd{x}=\bulksize}$. Therefore when Eq. (22) is substituted into Eq. (21), it simplifies to

\begin{align*} x_i(\zeta) & = -i\sum_{g}\ftsvarphi^*_i(g;\zeta)\left(\frac{\ftsvarphi_i(g-\hbulksize;\zeta)-\ftsvarphi(g;\zeta)}{\hbulksize}\right) \\ & = i\sum_g \ftsvarphi^*_i(g;\zeta)\partial_g\ftsvarphi_i(g;\zeta). \end{align*}

In the limit ${\bulksize\to \infty}$ the sum over $g$ can be expressed as an integral over the set of all real wavenumbers, ${\hrealone}$; i.e.,

\begin{align*} \sum_g=\frac{1}{\hbulksize}\sum_g \hbulksize \to \frac{\bulksize}{2\pi}\int_{\hrealone} \dd{g}. \end{align*}

Therefore,

\begin{align*} x_i(\zeta) & = i\frac{\bulksize}{2\pi}\int_{\hrealone} \ftsvarphi^*_i(g;\zeta)\partial_g\ftsvarphi_i(g;\zeta) \dd{g} \\ \implies \dv{x_i}{\zeta} & = i\frac{\bulksize}{2\pi}\int_{\hrealone}\bigg[\partial_\zeta\ftsvarphi^*_i(g;\zeta)\partial_g\ftsvarphi_i(g;\zeta) \\ &\qquad\qquad\qquad + \ftsvarphi^*_i(g;\zeta)\partial_\zeta\partial_g\ftsvarphi_i(g;\zeta)\bigg]\dd{g} \end{align*}

If the second term in the integrand is integrated by parts, after swapping the order of the derivatives with respect to $\zeta$ and $g$, the boundary term vanishes because ${\lim_{g\to \infty} \varphi_i(\pm g;\zeta)=0}$. Therefore we are left with

\begin{align*} \dv{x_i}{\zeta} & = -\frac{\bulksize}{2\pi} \times 2\im\left\{\int_{\hrealone} \partial_\zeta\ftsvarphi^*_i(g;\zeta)\partial_g\ftsvarphi_i(g;\zeta)\dd{g}\right\}. \tag{23} \end{align*}

By substituting this expression into Eq. (18) and choosing $V$ to be the entire bulk (i.e., ${[x_0,x_0+\bulksize]}$, where ${x_0\in\onetorus}$ is arbitrary), we find that

\begin{align*} \Jconvm & = \frac{e\dot{\zeta}}{\pi} \sum_i \im \left\{ \int_{\hrealone} \partial_\zeta\ftsvarphi^*_i(g;\zeta)\partial_g\ftsvarphi_i(g;\zeta)\dd{g}\right\}. \tag{24} \end{align*}

7.5.2 Calculating $\Jconv$ in $\onetorus$ for a crystal

Now let us assume that the material whose bulk is represented in $\onetorus$ is a crystal that comprises ${\NUnitcell}$ unit cells of widths $\volume$. Then ${\bulksize=\NUnitcell\volume}$ and ${\Nelecbulk=\NUnitcell\Neleccell}$, where ${\Neleccell}$ is the number of electrons per unit cell.

Since ${n(\zeta)}$ is ${\volume}$-periodic, $\varphiset$ can be chosen to be closed under translations by lattice vectors. In other words, if $\tvarphi(x;\zeta)$ is any element of $\varphiset$, and if $m$ is any integer, then the image ${\tvarphi(x-m\volume;\zeta)}$ of ${\tvarphi(x;\zeta)}$ under translation by the lattice vector ${m\volume}$ is also an element of ${\varphiset}$. If ${m}$ is an integer multiple of $\Nunitcell$, translating ${\tvarphi(\zeta)}$ by ${m\volume}$ maps it onto itself; otherwise it maps it onto an identical function whose center is in a different unit cell.

It follows that $\varphiset$ is the union,

\begin{align*} \varphiset = \bigcup_{m=0}^{\Nunitcell-1} \left\{\varphi_{i}(x-m\volume;\zeta):1\leq i \leq\Neleccell \right\}, \tag{25} \end{align*}

of ${\Nunitcell}$ subsets, where the indices have been chosen such that the center of every element of

\begin{align*} \varphisetcell\equiv \left\{\varphi_i(\zeta):1\leq i \leq \Neleccell\right\}\subset\varphiset \end{align*}

is in a particular unit cell, $\unitcell$. Therefore the centers of all elements of each subset in the union on the right hand side of Eq. (25) are in the same cell; and any one of the subsets can be transformed into any other one by translating all of its elements by the same lattice vector.

Eq. (25) implies that subset ${\varphisetcell}$ can represent set $\varphiset$ exactly, and the electronic structure of the bulk of the crystal exactly. Therefore the sum in Eq. (24) can be expressed as ${\Nunitcell}$ times the sum over the $\Neleccell$ elements of ${\varphisetcell}$.

It is shown in Appendix G that each element of ${\varphisetcell}$ is a Bloch function ${\bloch_{\alpha k}(\zeta)}$, which is associated with a particular element $k$ of the subset ${\BZ}$ of ${\reciplattg}$ whose elements are in the first Brillouin zone (see Appendix B). Therefore, let us replace index ${i}$ with ${k\alpha}$, where $\alpha$ distinguishes between different elements of ${\varphisetcell}$ that are associated with the same element of $\BZ$. Let us also make all dependences on $\zeta$ implicit, instead of denoting them implicitly.

Then Eq. (24) can be expressed as

\begin{align*} \Jconvm & = \frac{e\dot{\zeta}}{\pi} \sum_{\alpha k} \im \left\{ \int_{\hrealone} \partial_\zeta\ftsbloch^*_{\alpha k}(g)\partial_g\ftsbloch_{\alpha k}(g)\dd{g}\right\}, \tag{26} \end{align*}

where ${\ftsbloch_{\alpha k}}$ is the Fourier transform of Bloch function ${\bloch_{\alpha k}}$, and where $\sum_i$ has been replaced by ${ \Nunitcell\sum_{\alpha=1}^{\Neleccell}\sum_{k\in\BZ} }$. The factor ${\Nunitcell}$ disappears because $V$ is now a single unit cell, which is ${1/\Nunitcell}$ times the width of the entire bulk.

Bloch functions have the form

\begin{align*} \bloch_{\alpha k}(x)=e^{ikx}\pbloch_{\alpha k}(x), \tag{27} \end{align*}

where ${\pbloch_{\alpha k}}$ is $\volume$-periodic. This implies that

\begin{align*} \ftsvarphi_i(g)=\fts{\bloch}_{\alpha k}(g) = \fts{\pbloch}_{\alpha k}(g-k), \end{align*}

where ${\fts{\pbloch}_{\alpha k}}$ is the Fourier transform of ${\pbloch_{\alpha k}}$. The $\volume$-periodicity of ${u_{\alpha k}}$ implies that its Fourier transform ${\ftspbloch_{\alpha k}(g)}$ vanishes unless ${g}$ is a reciprocal lattice vector, which means that ${\fts{\bloch}_{\alpha k}(g)}$ vanishes unless ${g-k}$ is a reciprocal lattice vector. Therefore if ${\int_{\hrealone}\dd{g}=\hbulksize\sum_g}$ is expressed as ${\sum_{G}\int_{\BZ}\dd{k'}}$, where ${\sum_G=\sum_{G\in\reciplatt}}$ is a sum over all reciprocal lattice vectors and ${\int_{\BZ}\dd{k'}= \hbulksize\sum_{k'\in\BZ}}$, we find that the integral in Eq. (26) can be expressed as

\begin{align*} \int_{\hrealone} & \partial_\zeta\ftsbloch^*_{\alpha k}(g)\partial_g\ftsbloch_{\alpha k}(g)\dd{g} \\ &= \sum_G\int_{\BZ}\dd{k'} \partial_\zeta\ftsbloch^*_{\alpha k}(G+k')\partial_{k'}\ftsbloch_{\alpha k}(G+k)\dd{g}. \end{align*}

This vanishes unless ${k=k'}$ because ${\ftsbloch_{\alpha k}(G+k')}$ vanishes unless ${G+k'-k}$ is a reciprocal lattice vector, and because the largest difference between two elements of $\BZ$ is less than the magnitude ${\hreciplatt}$ of the smallest finite reciprocal lattice vectors. Therefore Eq. (26) can be expressed as

\begin{align*} \Jconvm & = \frac{e\dot{\zeta}}{\pi} \sum_{\alpha=1}^{\Neleccell} \im \left\{ \int_{\BZ} \dd{k}\sum_{G} \partial_\zeta\fts{u}^*_{\alpha k}(G) \partial_k\fts{u}_{\alpha k}(G)\right\}. \end{align*}

Now, since ${\partial_\zeta u_{\alpha k}}$ and ${\partial_k u_{\alpha k}}$ are ${\volume}$-periodic, we can express them as the Fourier series

\begin{align*} \partial_\zeta u_{\alpha k}(x) & = \sum_G \partial_\zeta\fts{u}_{\alpha k}(G)e^{iGx} \end{align*}

and

\begin{align*} \partial_k u_{\alpha k}(x) &= \sum_G \partial_k\fts{u}_{\alpha k}(G)e^{iGx}, \end{align*}

respectively. By substituting these series, it is straightforward to show that

\begin{align*} \braketT{\partial_\zeta u_{\alpha k}}{\partial_k u_{\alpha k}} & =\int_\onetorus \partial_\zeta u^*_{\alpha k}(x) \partial_k u_{\alpha k}(x)\dd{x} \\ & =\sum_G \partial_\zeta\fts{u}^*_{\alpha k}(G) \partial_k\fts{u}_{\alpha k}(G) \end{align*}

Therefore the polarization current can be expressed as

\begin{align*} \Jconvm & = \frac{e\dot{\zeta}}{\pi} \sum_{\alpha=1}^{\Neleccell} \im \bigg\{ \int_{\BZ} \dd{k} \braketT{\partial_\zeta u_{\alpha k}}{\partial_k u_{\alpha k}} \bigg\}. \tag{28} \end{align*}

Finally, note that the normalization of each periodic Bloch function ${u_{\alpha k}}$ is the same as its non-periodic counterpart, ${b_{\alpha k}}$, i.e.,

\begin{align*} \braketT{\pbloch_{\alpha k}}{\pbloch_{\alpha k}}=\braketT{\bloch_{\alpha k}}{\bloch_{\alpha k}} = 1. \end{align*}

However it is more common to express $\Jconvm$ in terms of periodic Bloch functions ${\tpbloch_{\alpha k}\equiv \sqrt{\Nunitcell}\pbloch_{\alpha k}}$ that are normalized on a single unit cell, i.e.,

\begin{align*} \braketT{\tpbloch_{\alpha k }}{\tpbloch_{\alpha k}}&=\Nunitcell \\ \implies \braketcell{\tpbloch_{\alpha k}}{\tpbloch_{\alpha k}}&\equiv \int_\unitcell\dd{x}\abs{\tpbloch_{\alpha k}(x)}^2 = 1. \end{align*}

Then,

\begin{align*} \Jconvm & = \frac{e\dot{\zeta}}{\pi} \sum_{\alpha=1}^{\Neleccell} \im \bigg\{ \int_{\BZ} \dd{k} \braketcell{\partial_\zeta \tpbloch_{\alpha k}}{\partial_k \tpbloch_{\alpha k}} \bigg\}, \tag{29} \end{align*}

which is one of the most commonly-quoted forms of the MTOP definition of $\Jconvm$.

I emphasize, again, that quantum mechanics has not been used in this derivation, or anywhere in Sec. 7. It was assumed that the density could be expressed as ${n(x;\zeta)=\sum_{i=1}^{M}\abs{\varphi_i(x;\zeta)}^2}$, where the functions ${\varphi_i}$ are mutually orthogonal and normalized to one, and ${M=\Nelecbulk}$. However they could also have different normalizations, and ${M}$ could be bigger or smaller than $\Nelecbulk$, as long as ${n(x;\zeta)}$ was normalized to $\Nelecbulk$ and the normalization of each $\varphi_i$ did not vary with $\zeta$. Such a representation of ${n(\zeta)}$ is possible for at least some classically-generated probability densities. For example, when the nuclei are treated as classical particles, their delta distribution can be represented in this way.

7.6 Interpretation of the MTOP

In Sec. 7.3 and Sec. 7.5 I have derived the MTOP expression for $\Jconv$ via routes that are very different to the one originally used to derive it; and I have shown that it is more widely applicable than previously believed. Previous derivations have used quantum mechanical perturbation theory, which is nothing more than a Taylor expansion of a microstate about a stationary state (i.e., an eigenstate) of a Hermitian operator. The applicability of such an expansion is not restricted to quantum mechanical systems.

I have shown that the MTOP can be applied to a system of classical identical particles in an evolving steady state $n(x;\zeta)$ that can be expressed as a sum of densities from orthogonal functions $\varphi_i$ whose normalizations do not change as $\zeta$ changes. Mathematically, this means that the MTOP is applicable whenever $n$ can be expressed as a sum of contributions from eigenstates of a self-adjoint one-particle operator, such as the solutions of any Sturm-Liouville equation, and when the subset of eigenstates contributing to $n$ does not change suddenly as $\zeta$ changes: Each contributing eigenstate must evolve smoothly with $\zeta$, while preserving its normalization.

An alternative derivation of Eq. (10), based on the homogenization theory proposed in Sec. 8 and Sec. 9, is presented and discussed in Sec. 12.6 and Sec. 13.

Although I have concluded that the MTOP definition of polarization current is correct, and exact, I am not satisfied that the MTOP solves the problem discussed in Sec. 5 and I do not agree with one interpretation of the MTOP; namely, that polarization is a fundamentally quantum mechanical phenomenon [Resta, 1993; Resta, 1994] and a property of the phase of a material’s wavefunction.

The wavefunction of an isolated material at equilibrium does not have a relevant phase. It is real, apart from an arbitrary and irrelevant constant phase factor. Therefore if $\pp$ were exclusively a property of a wavefunction’s phase, $\pp$ could only exist in an isolated material when it was being observed. This particular strain of quantum weirdness may be palatable to some, but it should be noted that delving more deeply into it leads to difficult philosophical questions.

Note that, by applying the kinetic energy operator to a wavefunction ${\Psi = \sqrt{\pdf}e^{i\theta}}$ whose phase ${\theta(\{x_i\})}$ is nontrivial, the phase can be shown to be equivalent to a potential ${\propto \sum_i\abs{\partial \theta/\partial x_i}^2}$. Although this potential is positive, potentials are only defined up to a constant applied to the entire physical system of interest. This means that any positive or negative potential is either equivalent to, or arbitrarily close to, the combination of a phase and this irrelevant constant. Therefore, when studying the bulk subsystem, the interaction between the bulk and surface subsystems could be expressed as a phase. Then ${\Jconv=\bsigmadot}$ would manifest as a time-dependence of this phase, which is equivalent to a changing interaction with the surface. However this is not the rationale used to relate $\Jconv$ to a phase in the MTOP.

Furthermore, it has sometimes been implied or stated that $\pp$ ‘is’, rather than ‘can be’, a property of the phase and that this means that it does not have an analogue within classical physics. This interpretation of the mathematics cannot be correct because the only assumption that I made in Sec. 7.5 to derive Eq. (26), which is identical to the MTOP expression for $\Jconv$, is that the electron density $n$ can be represented by a sum of contributions from smoothly-evolving mutually-orthogonal functions in an ${\lebesgue}$ Hilbert-Lebesgue space.

The ‘quantum’ MTOP definition of $\pp$ was developed to solve the problem discussed in Sec. 5 and the problem of how to calculate the currents demanded by asymmetry (Sec. 6.8). However, I did not invoke quantum mechanics to explain those problems because they are problems that exist within classical physics. This implies that, despite the development of the MTOP, there was more to understand because classical mechanics and classical statistical mechanics have never been derived from quantum mechanics. The correspondence principle is conjecture: Quantum mechanics is not known to be more general than classical statistical mechanics; it is assumed to be.

Classical physics is arguably no less internally-consistent than quantum physics. The problem with classical physics is not internal inconsistency but that, in its current state of development, it does not agree with experimental observations of very small or sensitive systems. Therefore, with the exception of a disagreement with experiment, any physical problem that can be stated within the classical realm, and which does not expose an internal inconsistency in classical physics, must have a solution within the classical realm.

The well-defined problem of how to calculate the polarization current from the changing charge densities in Fig. 2 is an excellent example. There exist classical processes that result in the equilibrium distribution of charges (e.g., ions in solution) changing. In most cases the charges are either too numerous or too sensitive to the act of observation to allow the net flow of charge to be calculated within classical deterministic mechanics, and one must turn to statistical mechanics. How can the current resulting from an adiabatic evolution of a classical equilibrium statistical state be calculated? The MTOP, as it was originally derived and interpreted, does not answer this question directly. But if a classically evolving charge density can be decomposed into packets whose integrals are constant, it answers it indirectly: $\Jconv$ can be calculated from Eq. (10).

My approach to deriving Eq. (10) in Sec. 12.3 and Sec. 13 is conceptually very simple. In Sec. 8 and Sec. 9 I use a systematic and unbiased approach to structure homogenization to derive expressions for interfacial excesses and changes in macroscopic fields across interfaces. From these I derive an expression for $\bsigma$ in terms of the microstructure (a generalization of Finnis’s expression), which leads to an expression for ${\J\equiv\bsigmadot}$ that is slightly more general than the MTOP expression, but otherwise equivalent to it.

I do not find any reason to retain the concept of a $\pp$ field, and it appears to violate macroscale symmetry. However, as I will discuss in Sec. 12.5.1, if one chooses to retain it, and if one also demands that ${\pp}$ determines the excess charge on a pristine surface of a perfect crystal, which is devoid of extrinsic charges, then, even within classical physics, $\pp$ must be quantized in the same way as the MTOP prescribes.