16. A potential paradox

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

In this section I highlight some subtleties in the meaning and definition of the microscopic potential $\phi$ and its relationships with its macroscopic counterpart ${\bphi}$ and the microscopic charge density $\rho$. As mentioned at the beginning of Sec. 15, the value of the MIP is believed to be positive [Sanchez and Ochando, 1985; Wilson et al., 1987; Wilson et al., 1988; Wilson et al., 1989; Pratt, 1992; Sokhan and Tildesley; Leung, 2010; Kathmann et al., 2011; Cendagorta and Ichiye, 2015; Blumenthal et al., 2017; Hörmann et al., 2019; Yesibolati et al., 2020; Madsen et al., 2021; Kathmann, 2021]. This contradicts my finding that it is zero. Therefore, to illustrate the subtleties, I use the example of Hans Bethe’s 1928 derivation of an approximate expression for the MIP, which is sometimes known as the Bethe potential (${\bphibethe}$), from several different perspectives. I begin, in Sec. 16.1, by outlining his derivation and line of reasoning.

My focus is on the ‘paradox’ referred to in the section title and I do not address the question(s) of most relevance to those using the Bethe potential, or one of its descendants, as a parameter in the analysis and/or interpretation of their calculations (e.g., theoretical electrochemistry) or experiments (e.g., electron holography). In most of these applications $\bphibethe$ is used as an approximation to the average potential experienced by an electron as it passes through the material. This quantity is likely to depend heavily on the electron’s energy as it enters the material and the time that it spends inside the material. Furthermore, one should not calculate it from the probability density ${n(\rvec)}$ of an electron being at ${\rvec}$, but on the conditional probability density ${n_c(\rvec_1|\rvec_2)}$ of there being an electron at ${\rvec_1}$ given that the probe electron is at ${\rvec_2}$.

As an illustration of the importance of basing calculations on ${n_c}$ rather than $n$, consider the example of a neutral atom meeting a stray electron in a vacuum. One might deduce from the atom’s electron density ${n(\rvec)}$ that they would not be attracted to one another; but by considering how the distribution of the atom’s electrons is changed by their interaction with the stray electron, one can quickly deduce that they do attract one another and that all singly-charged anions are stable in vacuum when the electrons’ temperature is sufficiently low.

16.1 Mean inner potential (Bethe's fallacy)

Bethe assumed that the charge densities of materials are not too dissimilar from a superposition of atomic charge densities. For a crystal with one spherically-symmetric atom in its primitive unit cell $\unitcell$, the expression he derived is

\begin{align*} \bphibethe = \frac{2\pi e}{3 \epsilon_0\volume}\int_0^\infty n(r) r^4 \dd{r} > 0, \tag{88} \end{align*}

where ${\volume}$ is the volume of $\unitcell$; and ${n(r)}$ is the number of electrons per unit volume in each atom’s electron cloud at a distance $r$ from its nucleus. Bethe deduced from Eq. (88) that $\bphi$ is positive, which contradicts my finding that it is zero. I will now rederive Eq. (88) via a more explicitly-careful mathematical route than Bethe chose to present, but using his starting point and physical reasoning. His starting point was the expression

\begin{align*} \phir(r) \equiv \frac{1}{\epsilon_0 r}\int_0^r \rho(u) u^2\dd{u} + \frac{1}{\epsilon_0}\int_r^\infty \rho(u) u\dd{u}, \tag{89} \end{align*}

for the electric potential ${\phir}$ at a distance $r$ from the center of an isolated spherically-symmetric charge distribution ${\rho(r)}$. Eq. (89) can be derived from Gauss’s law by assuming that the electric field inherits spherical symmetry from $\rho$ and by expressing the potential at a distance $r$ from the nucleus, $\phi_r(r)$, as the integral of the spherically-symmetric field from an infinitely-distant point to one whose distance from the nucleus is $r$, along an axis passing through the nucleus.

Bethe reasoned that the average potential in the crystal is the potential emanating from one atom, integrated over all points in space, and divided by the volume per atom, i.e., the ${R\to \infty}$ limit of

\begin{align*} \phibar(R) & = \frac{4\pi}{\volume}\int_0^R \phir(r) r^2 \dd{r}. \tag{90} \end{align*}

His reason for integrating $\phir$ over all space, but dividing by the volume of only one unit cell, was that the potential in each cell has contributions from atoms in all other cells and, either by symmetry, or when averaged over all other cells, the sum of the contributions of the atom in a given cell to the averages of the electrostatic potential in all other cells must equal the sum of the contributions of atoms in all other cells to the average electrostatic potential in the given cell.

Bethe expressed each atom’s charge density (charge per unit volume at a distance $r$ from its center) as ${\rho(r)=\rhop(r)+\rhom(r)}$, where ${\rhop}$ is the density of nuclear charge and ${\rhom(r) = -e \, n(r)}$ is the density of electron charge. He expressed ${\rhop}$ as a delta distribution, i.e.,

\begin{align*} \rho(r) = Z\, e\, \delta(r) - e\,n(r), \tag{91} \end{align*}

but I will keep it more general for now. Substituting Eq. (89) into Eq. (90) and simplifying leads to

\begin{align*} \phibar(R) & = -\frac{2\pi}{3\epsilon_0\volume}\int_0^R \rho(r) r^4 \dd{r} \\ &+ \frac{\frac{4}{3}\pi R^3}{\volume}\left[\frac{3}{2}\frac{\kappa Q(R)}{R} + \frac{1}{\epsilon_0}\int_R^\infty \rho(r) r\dd{r}\right] \tag{92} \end{align*}

where

\begin{align*} Q(R)\equiv 4\pi\int_0^R \rho(r) r^2 \dd{r} \end{align*}

is the net charge in a sphere of radius $R$. When $R$ is chosen large enough that the total charge outside this sphere is negligible, only the first term on the right hand side remains, i.e.,

\begin{align*} \phibar(R) & = -\frac{2\pi}{3\epsilon_0\volume}\int_0^R \rho(r) r^4 \dd{r} \end{align*}

Since ${\rho(r)<0}$ when ${r}$ exceeds the spatial extent of ${\rho_+}$, ${\phibar(R)}$ is positive and we recover Eq. (88) in the limit of large $R$ if ${\rhop}$ is localized at a point.

If, instead, we assume that ${\rhop}$ has a finite width and denote the total nuclear charge by ${Q_+}$, we can use the atom’s overall charge neutrality, i.e.,

\begin{align*} Q_+\equiv 4\pi\int_0^\infty u^2\rhop(u)\dd{u} = -4\pi\int_0^\infty u^2\rhom(u)\dd{u}, \end{align*}

to express $\bphibethe$ as

\begin{align*} \bphibethe& = \frac{Q_+}{6\epsilon_0\volume}\left(s_-^2-s_+^2\right) \tag{93} \end{align*}

where $s_+^2$ and $s_-^2$ are the mean squared distances of positive and negative charges, respectively, from the atom’s center, i.e.,

\begin{align*} s_\pm^2 & \equiv \pm \frac{4\pi}{Q_+}\int_0^\infty \rhopm(u) u^4 \dd{u} \end{align*}

From Eq. (93) it seems clear that the MIP being positive is a consequence of the electrons being more delocalized than the nuclei (${s_-^2>s_+^2}$). For example, if we assume that ${\rhop}$ and ${\rhom}$ are constant within concentric spheres of radii ${r_+}$ and ${r_-}$, respectively, and zero outside them, then Bethe’s derivation would lead to

\begin{align*} \bphibethe = \frac{Q_+}{10\epsilon_0\volume}\left(r_-^2 - r_+^2\right) > 0. \end{align*}

It seems reasonable to interpret this expression as follows: Gauss’ Law implies that ${\phir(r)}$ vanishes if ${r>r_-}$ because the net charge within a sphere of radius $r_-$ centered at the atom’s center is zero. The potential is positive at the atom’s center and it decreases monotonically to its value of zero at ${r=r_-}$. Therefore the potential is positive in a sphere of radius ${r_-}$ and zero everywhere else. It seems obvious, then, that ${\bphi}$ must be positive.

I will now show, by illustration, why this obviously-right result must be wrong. Then I will explain the flaw in Bethe’s reasoning and show that a more careful treatment of the problem leads to the conclusion that ${\bphi}$ is zero. I illustrate the flaw from several different perspectives to highlight some of the many pitfalls that exist when working with the electric potential. Readers who have already spotted the flaw, or who don’t like whodunnits, might want to skip the illustrations and proceed directly to Sec. 16.2.

**Figure 16.** If a crystal is built from a superposition of atoms (Crystal 1), Bethe's derivation leads to the result ${\bphibethe_1>0}$. Crystal 2, shown schematically above, is built from a charge-neutral building block comprised of an electron cloud surrounded by fractions of nuclei (bottom). Bethe's derivation leads to the result ${\bphibethe_2<0}$ for Crystal 2. Both crystals are identical in the bulk (yellow region), but differ at their surfaces; however the surfaces of both crystals are charge-neutral (${\bsigma=0}$). Crystal 2 can be transformed into Crystal 1 by cancelling each of the outermost ${+Ze/4}$ charges, where $Z$ is the atomic number, with a charge of ${-Ze/4}$ at the same position and adding a charge of ${+Ze/4}$ to the centers of the outermost electron clouds. If we assume this to be equivalent to adding a dipole moment density (dipole moment per unit length, in this 2-d example) of ${-\left(Ze/4\right) \hat{n}}$ to each surface, where ${\hat{n}}$ is the surface's outward unit normal vector, it is easy to calculate the effect of such a layer on the MIP for the 3-d analogue of the crystal shown. It would shift the average potential below the surface relative to its value in the vacuum above the surface by exactly ${\bphibethe_1-\bphibethe_2}$, thereby making the values of the MIP calculated for Crystal 1 by two different routes equal. However, at the macroscale, adding this dipole moment density should have no effect on the value of ${\bphi}$ because a microscopic distance is indistinguishable from a distance of zero at the macroscale (see Sec. 8). Therefore an isolated dipole moment shrinks to a point of no net charge under the homogenization transformation and the addition of a plane of such dipoles to a surface does not change either the net charge ${\bsigma}$ at the surface or the macroscopic charge density ${\Rho}$ at any point inside or outside the material.

16.1.1 Existence of a flaw - Illustration 1

Bethe chose to build his material from a spherically-symmetric charge density, with a localized distribution of positive charge at its center and relatively delocalized distribution of negative charge surrounding it. However, just as there is no ‘right’ way to partition the charge density of a material for the purpose of defining its average dipole moment density (see Sec. 5.0.1), there is no right way to partition it for the purpose of calculating its average potential. It is no less justified to build a crystal from a superposition of charge densities of the form

\begin{align*} \rho(\rvec) = -e\,n(r) + \frac{Q_+}{N_A}\sum_{i: |\Rvec_i|=A} \delta(\rvec-\Rvec_i) \tag{94} \end{align*}

where, as before, the origin is chosen to coincide with a nucleus; ${r=\abs{\rvec}}$; each ${\Rvec_i}$ is a lattice vector and therefore a relative displacement of two nuclei (for simplicity I assume that ${T\to 0}$); and the sum is over all ${\Rvec_i}$’s of a given finite magnitude ${A}$, of which there are a total of ${N_A=\sum_{i:|\Rvec_i|=A} 1}$.

An example of a building block of this form is shown schematically for a 2-d crystal in Fig. 16. It has negative charge at its center, positive charge further away, and it is charge neutral overall, meaning that the net flux of the electric field through any surface that encloses it is zero, as it is for an atom. The flux from an atom is zero at all points on a surface enclosing it, whereas the flux from the charge density of Eq. (94) is finite almost everywhere on a surface enclosing it, but with regions of the surface where it is positive and regions where it is negative. Nevertheless, Gauss’s law implies that the net potential outside the surface from charge within it is zero, as it is for an atom.

Using the same physical reasoning with which Bethe deduced that ${\bphibethe>0}$ for a material built from atoms, it can be shown that ${\bphibethe<0}$ in a material built from this charge distribution, because ${s_+=A > s_-}$. Furthermore, the magnitude of ${\bphibethe}$ depends on the value of $A$. For example, consider a simple cubic crystal with lattice spacing ${a}$ and let us build it from Eq. (94) with the choice ${A=a}$ (${\implies N_A=6}$). I will refer to the crystal built in this way as Crystal 2, I will refer to the crystal built by Bethe from atoms as Crystal 1, and I will denote their MIPs, as derived using Bethe’s approach, by ${\bphibethe_2}$ and ${\bphibethe_1}$, respectively. Then if, following Bethe, we assume the nuclei to be localized at a point, we find that

\begin{align*} \bphibethe_2=\bphibethe_1-Q_+ \frac{a^2}{6\epsilon_0\volume}. \end{align*}

Crystal 1 and Crystal 2 are identical in the bulk; they differ only near surfaces. However, because all surfaces of each crystal are charge-neutral (${\bsigma=0}$), and because ${\Rho=0}$ in the bulk in each case, the macroscale theory would not be internally consistent if ${\bphi>0}$ in one case and ${\bphi<0}$ in the other. If the value of ${\bphi}$ is defined it must be the same in each case because the macrostructures of the two crystals are identical: electrostatically, each one is indistinguishable from empty space.

Crystal 1 and Crystal 2 could be made identical by adding a pair of charges of opposite signs, of magnitudes ${Q_+/6}$, and separated by a distance $a$, to each surface unit cell of one of the crystals. For example, to make Crystal 2 into Crystal 1 we would have to add a charge of ${Q_+/6}$ to the center of each of the electron clouds closest to its surface and a charge of ${-Q_+/6}$ at a displacement ${a\,\hat{n}}$ from the first charge, where ${\hat{n}}$ is the surface’s outward unit normal. Let us temporarily assume that, from the perspective of a point whose depth below the surface is much greater than $a$, this is equivalent to adding to the surface an approximately-uniform areal density of dipole moments,

\begin{align*} \sigma_\mpp=-\left(Q_+ a/6 a^2\right) \hat{n} = -\left(Q_+ a^2/6\volume\right)\hat{n}. \end{align*}

Then the result would be an upward shift of the potential in the crystal relative to the vacuum near the surface of

\begin{align*} \Delta\phi = \sigma_\mpp/\epsilon_0=Q_+ a^2/6\epsilon_0\volume. \end{align*}

This cancels the difference between ${\bphibethe_2}$ and ${\bphibethe_1}$. Therefore, although the values of ${\bphibethe}$ derived by Bethe’s method differ, it appears possible to correct the difference between them by changing the surface structure of either crystal to make the two crystals identical.

Unfortunately, although we have corrected the difference between the two values of ${\bphibethe}$, this does not solve our problem. We still do not have any reason to prefer one building block over another; therefore we do not have any reason to prefer one of the resulting crystal surface structures over the other. We have two derivations, which appear equally valid, and from which we deduce two different values of the MIP. This appears to imply that there is a flaw in the construction that Bethe used for his derivation.

Bethe did not involve surfaces in his derivation because, when ${\bsigma=0}$, he regarded the MIP as a property of the bulk. However, the fact that ${\bphibethe_1}$ and ${\bphibethe_2}$ differ suggests that the MIP is a surface property, which can be changed by adding or removing equal amounts of positive and negative charge at each surface. This is problematic if we wish to identify ${\bphibethe}$ as the macroscopic potential ${\bphi}$ because, as mentioned above, the addition of a layer of microscopic dipoles to a surface should not change ${\bphi}$ in an internally-consistent linear macroscale theory. This is because, at the macroscale, a microscopic distance is equivalent to ($\Lequiv$) a distance of zero (see Sec. 8 and Sec. 9), and because a layer of microscopic dipoles, ${qa\hat{n}}$, is equivalent to two layers with equal and opposite charges per unit area that are separated by a distance $a$. Since ${a\Lequiv 0}$ these two layers are equivalent, at the macroscale, to a single charge-neutral layer, which would not change ${\bphi}$ inside the crystal. Therefore if Bethe’s derivation is right, and if ${\Rho}$ is a linear spatial average of $\rho$, the MIP cannot be identified as the macroscopic potential because that would be tantamount to saying that the same macroscale distribution of charge can give rise to different values of ${\bphi}$. This would imply that, even when $\Rho$ and the macroscale boundary conditions are known, the value of $\bphi$ cannot be calculated; its value depends, in some way, on certain microscopic details of $\rho$ that are lost by the ${\rho\mapsto\Rho}$ homogenization transformation.

If we can assume that ${\phi}$ is a linear functional ${\phi[\rho]}$ of $\rho$, and that $\Rho$ is a linear spatial average $\expval{\rho}$ of $\rho$, then the linearity of both operations implies that ${\phi[\Rho]= \phi[\expval{\rho}]=\expval{\phi[\rho]}}$. Therefore, if there exist two microscopic charge densities, $\rho$ and ${\rho+\Delta\rho}$ , with the same macroscopic charge density ${\Rho}$ (${\implies\expval{\Delta\rho}=0}$) and different macroscopic potentials, $\bphi$ and ${\bphi+\Dbphi}$, then ${\bphi}$ does not equal ${\phi[\Rho]}$ and is a nonlinear functional ${\bphi[\rho]}$ of $\rho$. Linearity would imply that

\begin{align*} \Dbphi\equiv\bphi[\rho+\Delta\rho]-\bphi[\rho]=\bphi[\Delta\rho]\neq 0; \end{align*}

and it would also imply that the mesoscale spatial average,

\begin{align*} \expval{\Dbphi}=\expval{\bphi[\Delta\rho]}=\bphi[\expval{\Delta\rho}], \end{align*}

is zero because ${\expval{\Delta\rho}=0}$. Therefore ${\Delta\bphi}$ would be a harmonic function (${\mathbf{\laplacian}\,\Dbphi=0}$) that fluctuates microscopically about zero. It would follow that ${\Delta\bphi}$ and $\bphi$ are microscopic quantities, not a macroscopic ones. On the other hand, nonlinearity of ${\bphi[\rho]}$ implies that a material’s macroscopic potential can depend on its history; for example, if the material’s microstructure ${\rho=\rho_1+\Delta\rho_1=\rho_2+\Delta\rho_2}$ is built by superimposing the charge densities ${\rho_1}$ and ${\Delta\rho_1}$, its macroscopic potential would differ, in general, from its value if it was built by superimposing the two different charge densities, ${\rho_2}$ and ${\Delta\rho_2}$. There are many problems that arise if we are tempted to assume that ${\bphi}$ depends on microscopic details of $\rho$ that are washed away by the homogenization transformation; I have only mentioned a few of them.

Returning to the example of Crystal 1 and Crystal 2: if we rigidly shift the MIP of Crystal 2 by ${\Delta \phi>0}$ by coating its surfaces with a layer of dipoles, the same layer would shift the average potential in the vacuum just outside the crystal by ${-\Delta \phi}$. Outside Crystal 1, $\phi$ appears to be zero because the field from each atom is zero. Outside Crystal 2, the average of $\phi$ appears to be zero because the average electric field emanating from each building block is zero. Adding a dipole layer to Crystal 2 to turn it into Crystal 1 appears to shift the mean vacuum potential (MVP) in a layer surrounding the crystal up, while shifting $\bphibethe$ down by the same amount. In the limit of large distance from the crystal the potential vanishes, because the crystal is charge neutral overall, but it does not begin to decrease in magnitude significantly until the distance to the closest surface is comparable to one or more of the surface’s linear dimensions; therefore the MVP is shifted by ${-\Delta \phi}$ in a macroscopic layer of vacuum surrounding the crystal. So if Bethe’s derivation was correct, and if a macroscopic layer of microscopic dipoles could shift the average potential in a macroscopic region, the MVP would be zero in a macroscopic layer of vacuum surrounding the crystal both before and after it had been shifted by a finite amount ${-\Delta \phi}$! Clearly, this is absurd.

Now consider Fig. 17, which uses the concept of a dipole layer to illustrate one argument for why the MIP is positive. The crystal in question is identical to Crystal 1, so this construction appears to validate Bethe’s result. However, there is no justification for dividing the material’s microstructure into the blue and pink layers. If, for example, we combined each adjacent pair of pink and blue layers into a single charge-neutral layer, we would find that the MIP vanishes. Therefore, as with the construction Bethe used in his derivation, two equally-justified ways to partition and spatially-average the microstructure leads to two different values of the MIP.

The superposition principle, on which Bethe’s derivation and much of electromagnetic theory are based, allows us to do the following: Let us partition the space $\materialregion$ occupied by an electron cloud into $M$ partitions of volume ${\materialregionvolume/M}$ and let us divide the nucleus’s charge into $M$ ‘pieces’, such that for each partition there is a piece of nucleus with the same magnitude of charge. Now, after taking the large-$M$ limit, let us displace the pieces to the partitions so that each one becomes charge neutral. The atom’s spherical symmetry initially ensures that, after displacing all of the pieces, it is again spherically symmetric. After this redistribution of charge, the MIP must be zero because the nucleus and every partition have become charge neutral. If we view the displacement of nuclear charge as the superposition of the negative of an atom’s charge density on each atom, this makes sense. We have simply superimposed a crystal’s charge density and its negative, so of course the MIP of the superposition vanishes. However, we could also view the displacement of each piece of nucleus as the placement of a negative charge at the nucleus and a positive charge in the partition. Placing a dipole at a point in space changes the potential everywhere but, by symmetry, it cannot change the spatial average of the potential. Therefore placing all of these dipoles inside the crystal should not change the MIP. It appears that the superposition principle does not apply.

**Figure 17.** One can think of each of the coloured layers parallel to the surface as planes of charge, with the layers coloured pink carrying charge densities of ${\sigma>0}$ and those coloured blue carrying charge densities of ${-\sigma}$. Alternatively, one can think of each blue layer and the red layer next to it as a plane of dipoles. The plane of dipoles closest to the surface is pointing into the surface, causing an upward step in $\phi$, which in this example denotes the planar average of the microscopic potential, relative to its value of ${\phi_{\text{vac}}}$ in the vacuum close to the surface (dotted green line). The next layer causes a downward step of $\phi$, and so on. The average potential is clearly positive if ${\phi_{\text{vac}}}$ is regarded as the zero of potential. It would be negative if the layer closest to the vacuum was positive instead of negative, as in Fig. 16, and if the value of ${\phi_{\text{vac}}}$ was again set equal to zero, despite its value differing from its value above the surface of the crystal pictured above.

16.1.2 Existence of a flaw - Illustration 2

Another way to see that there must be problem with Bethe’s result is to treat electrons as point particles instead of expressing ${\rhom}$ as a smooth and delocalized density. Using a line of reasoning similar to Bethe’s we could say that the total potential in each unit cell from all electrons and nuclei outside the cell is approximately equal to the sum of the total potentials emanating from the particles inside the cell. Then we could calculate the total potential from each point particle at all points within a distance $R$ of it, add together the total potentials from all particles within each unit cell of the crystal, and take the limit of large $R$ to get the total potential emanating from each unit cell. This total would vanish because the potential emanating from a point charge ${Ze}$ is the negative of ${Z}$ times the total potential emanating from a point charge ${-e}$. Cancellation is obvious when ${Z=1}$ (hydrogen), but Bethe’s construction does not exclude this case, either explicitly or implicitly. Therefore his derivation leads to the conclusion that the magnitude of the spatial average of the potential from a proton is greater than the magnitude of the spatial average of the potential from an electron.

This suggests that the problem in Bethe’s derivation might be related to his use of a continuous charge density for electrons and a (discrete) delta distribution of charge for nuclei. Usually this form of $\rho$ is regarded as a time average of the true charge distribution: electrons are whizzing around the more massive nuclei so fast that their charge, when observed on a timescale of about ${10^{-16}-10^{-15}}$ seconds, appears to be smeared into a continuous charge density. This timescale is too short for nuclei to move significantly, but long enough for each electron to trace out a very long trajectory. Nevertheless, because the integral of

\begin{align*} \phi\left(\uvec \right) \equiv \overbrace{\sum_{i\in\text{nuclei}}\frac{\kappa Z e}{\abs{\uvec-\Rvec_i}}}^{\displaystyle \phi_n(\uvec)} \overbrace{- \sum_{j\in\text{electrons}}\frac{\vphantom{Z}\kappa e}{\vphantom{\Rvec}\abs{\uvec-\rvec_j}}}^{\displaystyle \phi_e(\uvec)} \end{align*}

over all space vanishes, which implies that its average over all space, ${\expval{\phi(\uvec)}_{\uvec}}$, vanishes, one might expect the spatial average ${\expval{\expval{\phi(\uvec)}}_{\uvec}}$ of its expectation value with respect to the material’s wavefunction $\Psi$,

\begin{align*} \expval{\phi(\uvec)}\equiv\frac{\expvaltwo{\phi(\uvec)}{\Psi}}{\braket{\Psi}{\Psi}}, \end{align*}

to also vanish. However, it is easy to reduce this expectation value and its average over all $\uvec$ to the forms

\begin{align*} \expval{\phi(\uvec)} = &\frac{\expvaltwo{\phi(\uvec)}{\Psi}}{\braket{\Psi}{\Psi}} =\phi_n(\uvec) -\kappa\,e\, \int \frac{n(\rvec)}{\abs{\uvec-\rvec}}\ddpow{3}{r}, \\ \expval{\expval{\phi(\uvec)}}_{\uvec} & = \expval{\phi_n}_{\uvec} + \expval{-\kappa\,e\, \int \frac{n(\rvec)}{\abs{\uvec-\rvec}}\ddpow{3}{r}}_{\uvec}, \tag{95} \end{align*}

where ${n(\rvec)}$ is the probability density that there is an electron at ${\rvec}$. At first glance, Eq. (95) might appear to validate Bethe’s approach, because it appears to be the spatial average of the potential from a charge distribution whose form becomes equivalent to the one he used (Eq. (91)) when there is spherical symmetry. It would be very strange if it were equivalent: it would mean that ${\expval{\expval{\phi(\uvec)}}_{\uvec}}$ is finite but that the expectation value ${\expval{\expval{\phi(\uvec)}_{\uvec}}}$ vanishes. Therefore, it would mean that its value is changed simply by changing the order of integration such that the integral with respect to ${\uvec}$ is performed first. We need to understand this better - both physically and mathematically.

16.2 Flaw in Bethe's reasoning

The flaw in Bethe’s derivation is that he calculated the electrons’ contribution to the potential from a volumetric density of negative charge, ${\rhom(\rvec) = -e\,n(\rvec)}$, which is defined at all points in space. Then, because the electron density is spread over a greater volume than the nuclear density on femtosecond time scales, the spatial average of the potential does not vanish.

To understand why it does not vanish, let us again consider the potential ${\phi_r(r;\eta)}$ at a distance $r$ from the center of a spherically-symmetric nonpositive or nonnegative charge density ${\rho(u;\eta)}$, where the width of ${\rho}$ is proportional to the value of $\eta$, which is a parameter. Let us denote the integral of $\rho$ within a sphere of radius $r$ by ${Q(r;\eta)}$ and its integral over all space by

\begin{align*} Q_\infty\equiv\lim_{r\to \infty}Q(r;\eta)= \lim_{\eta\to 0}Q(r;\eta). \end{align*}

Using Eq. (89) we can express the magnitude of the potential at a distance ${r}$ from the center of $\rho$ as

\begin{align*} \abs{\phi_r(r;\eta)} & = \abs{\frac{\kappa Q(r;\eta)}{r} + \frac{1}{\epsilon_0}\int_r^\infty \rho(u;\eta) u\dd{u}} \\ & < \abs{\frac{\kappa Q(r;\eta)}{r} + \frac{1}{\epsilon_0}\int_r^\infty \rho(u;\eta) u\left(\frac{u}{r}\right)\dd{u} } \\ \therefore \abs{\phi_r(r;\eta)} & < \frac{\kappa\abs{Q_\infty}}{r}. \end{align*}

At very large distances (${r\gg \eta}$) the potential is approximately equal to ${\kappa Q_\infty/r}$, and it gets closer to this value as $r$ increases. At short distances (${r \sim \eta}$) the magnitude of the potential is significantly smaller than ${\kappa \abs{Q_\infty}/r}$ and the ratio ${\phi_r(r;\eta)/\left(\kappa Q_\infty/r\right)}$ gets smaller as $r$ decreases.

Now let us consider the average of $\phi_r$ over all points within a fixed distance $R$ of the center of $\rho$ as the value of $\eta$ changes. When ${\eta/R}$ is very small, the magnitude of $\phi_r$ at almost all points is approximately ${\kappa \abs{Q_\infty}/r}$ and it is only significantly smaller than that value in a volume fraction ${\sim \left(\eta/R\right)^3}$ of the sphere. In the limit ${\eta/R\to 0}$, the average potential is the same as it would be if $\rho$ was the delta distribution of a point charge. However, as $\eta$ increases, the magnitude of the average potential in the sphere of radius $R$ decreases because the fraction of the volume occupied by points at which $\abs{\phi_r}$ is significantly less than ${\kappa \abs{Q_\infty}/r}$ increases. Therefore, the average potential in the sphere reduces in magnitude as $\eta$ increases.

This is why the potential from the electrons does not cancel the potential from the nuclei in Bethe’s derivation: the value of ${\eta}$ is finite for the electrons, but vanishingly small for nuclei, which makes the magnitude of the average potential from the nuclei greater than that from the electron cloud.

I will now explain, from three different perspectives, what is wrong with Bethe’s derivation and with his use of the charge density in Eq. (91).

16.2.1 Perspective 1

Bethe’s use of a continuous electron charge density ${\rhom(\rvec)=-e\,n(\rvec)}$ suggests that he interpreted it as the time average of the electrons’ instantaneous delta distribution of charge. However ${n(\rvec)}$ is not a time average of the electrons’ positions, it is a probability density that an electron (any one of them) is at position $\rvec$ at any precisely-specified time. The average, over a time interval ${\interval(t,\Delta t)}$, of the delta distribution of a set of moving particles is not a volumetric charge density, but a set of linear charge densities defined only along the segments of the trajectories followed during ${\interval(t,\Delta t)}$. Therefore the time average of the electrons’ delta distribution is a set of linear charge densities defined on a set of curves and there is no charge at points that do not lie along these curves.

This means that the set of points at which the time-average of the electrons’ charge distribution is nonzero is a set whose measure in $\realthree$ is zero, regardless of the magnitude of ${\Delta t}$; therefore electrons are not more delocalized than nuclei because both occupy zero volume. Using the fact that ${Ze/r}$ is cancelled by ${Z\times -e/r}$, it is easy to show that the potential from the true time average of the electrons’ charge distribution exactly cancels the potential from the nuclei.

16.2.2 Perspective 2

Although Eq. (95) appears to be the spatial average of the potential from a set of point nuclei and a continuous density of negative charge ${\rho_-(\rvec)=-e\,n(\rvec)}$, it is not. The quantity

\begin{align*} \expval{\phi_e(\uvec)} = -\kappa e \int \frac{n(\rvec)}{\abs{\uvec-\rvec}}\ddpow{3}{r} \tag{96} \end{align*}

is not the potential at $\uvec$ from ${\rhom(\rvec)=-e\,n(\rvec)}$ because Coulomb’s law does not hold for the charge from a continuous charge density within a region of infinitesimal volume. To see that it does not hold, consider a spherical Gaussian surface of radius $\Delta r$ centered at the point $\rvec$. The charge within the sphere is ${\rho(\rvec)\times \frac{4}{3}\pi(\Delta r)^3}$, but the electric field, $\me$, on its surface is not directed radially outward from $\rvec$ because it has a contribution from the nucleus. This would not present a problem for point charges because one could choose ${\Delta r}$ to be arbitrarily small without changing the amount of charge enclosed by it; and, as ${\Delta r}$ got smaller, the direction of the field passing through the surface would become arbitrarily close to radially outward.

However this reasoning does not apply to a continuous charge density because the magnitude of the charge enclosed by the surface scales like ${\left(\Delta r\right)^3}$ in the small ${\Delta r}$ limit. A correct application of Gauss’s law for a spherically-symmetric charge density leads to Eq. (89), which does not give the same result as Eq. (96). Eq. (96) is the correct expression for the expectation value of the potential at ${\uvec}$ from the electrons, but it is not the correct expression for the potential from charge density ${\rhom}$; when there is spherical symmetry, the latter is Eq. (89).

Eq. (95) can be expressed in the slightly more general form

\begin{align*} \expval{\expval{\phi(\uvec)}}_{\uvec} & = \kappa \expval{\int \frac{\rho(\rvec)}{\abs{\uvec-\rvec}}\ddpow{3}{r}}_{\uvec} \end{align*}

where ${\rho=\rhop+\rhom}$ and $\rhop$ is the distribution of the nuclei. Assuming that $\rho$ has spherical symmetry, and choosing its center as the origin, the integral of ${\expval{\phi(\uvec)}}$ over all points within a distance $R$ of its center is

\begin{align*} \int_{\abs{\uvec}<R} \expval{\phi(\uvec)} \ddpow{3}{\uvec} & = \kappa \int_{\abs{\uvec}<R} \ddpow{3}{u}\int_{\realthree} \ddpow{3}{r}\left(\frac{\rho(\rvec)}{\abs{\uvec-\rvec}}\right) \\ & = \kappa \int_{\realthree} \ddpow{3}{r} \rho(\rvec)\int_{\abs{\uvec}<R} \frac{1}{\abs{\uvec-\rvec}}\ddpow{3}{u} . \end{align*}

In the limit of large $R$ the right hand side becomes arbitrarily close to

\begin{align*} \kappa \left(\int_{\realthree} \rho(\rvec) \ddpow{3}{r} \right) \left(\int_0^R 4 \pi u^2 \times \frac{1}{u} \dd{u}\right), \end{align*}

which vanishes because ${\int_{\realthree} \rho(\rvec)\ddpow{3}{r}}$ vanishes. Therefore the total potential from the atom, divided by the volume of a primitive unit cell, which is the quantity calculated by Bethe, is zero. Furthermore, to deduce that it is zero I have not assumed anything about the degrees to which the nuclear density and the electron density ${n(\rvec)}$ are localized.

16.2.3 Perspective 3

The relations ${\me=-\grad \phi}$ and ${\rho/\epsilon_0 = -\laplacian\phi}$ are preserved by the homogenization transformation because $\grad$, $\laplacian$, and the spatial average are all linear operations, which commute when they are applied in a mutually-consistent manner. For example, if $\expval{\;}_x$ and $\expval{\;}_y$ denote averages along the $x$-axis and $y$-axis, respectively, then

\begin{align*} \expval{\rho}_x/\epsilon_0 = -\expval{\laplacian\phi}_x=-\laplacian\expval{\phi}_x. \end{align*}

However, it cannot generally be true that ${\expval{\rho}_x/\epsilon_0 = -\laplacian\expval{\phi}_y}$. Calculating the average potential, along an axis normal to a surface, from the average of the charge density on planes parallel to the surface does not, in general, lead to a meaningful result.

This is why it is not physically reasonable to calculate the MIP from the average of the green curve in Fig. 17 and it is why the value of the MIP deduced by averaging the charge distribution in layers parallel to the surface depends on the choice of the layers’ positions and thicknesses. For example, if, instead of the division into the pink (P) and blue (B) layers depicted in Fig. 17, each layer was chosen to be a layer of atoms, the green curve would be flat because each layer would be charge neutral. A layer of atoms comprises two pink layers and two blue layers in the order BPPB, so I will call it a BPPB-layer. One could choose the first layer at the surface to be a negatively charged B-layer and all others to be charge-neutral PPBB-layers. In that case the green curve would be a straight line with a negative slope; therefore, not only would the MIP appear to be negative, there would be a macroscopic electric field in the material emanating from the plane of negative charge. If the first layer is a BPP-layer and subsequent layers alternate between B- and BPP-layers, the MIP appears to be negative, with the potential as a function of depth resembling a skewed version of the negative of the green curve in Fig. 17.

The electron charge density used by Bethe in his derivation can be regarded as the result of performing the following sequence of temporal and spatial averages: first, the time average of the electron charge distribution is calculated to give a set of curves carrying linear charge densities; next, this set of linear charge densities are turned into a volumetric charge density by averaging in the radial direction over a small width ${\dd{r}}$; finally, the resulting distribution is given spherical symmetry by setting the charge density at distance $r$ equal to its spatial average on the spherical surface of radius $r$. The resulting charge density is then used to calculate the potential as a function of position along the radial direction, i.e., along an axis that, at its point of intersection with the surface on which the final spatial average is taken, is perpendicular to this surface. There is no reason to expect this procedure to produce results that are any more meaningful than those derived by partitioning the surface in Fig. 17 into artificial layers using an arbitrary, unjustified, and mutually-inconsistent sequence of partial spatiotemporal averages.

I avoided the problems with Bethe’s derivation in Sec. 15 by not making any assumptions about the microscopic charge density, except that it is mathematically smooth; and by only averaging along the $x$-axis; and by using the same mesoscopic interval width for all spatial averages. Because I used a general form of $\rho$, the derivations of Sec. 15 apply to charge distributions that are arbitrarily close to delta distributions; and Coulomb’s law can be applied to a smooth charge distribution in this limit.