12. Surface charge ($\bsigma$)

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

A charged surface or interface is not stable unless the electric potential from it is compensated by, for example, an oppositely charged surface or interface. The instability of isolated charged surfaces is due to the divergence of the electric potential (see Sec. 15). If a pristine isolated crystal surface is charged, and therefore unstable unless neutralized by a change in its composition with respect to the bulk, it is classified as polar. An important question, about which a great deal has been written [Tasker, 1979; Finnis, 1998; Noguera, 2000; Goniakowski et al., 2008; Stengel and Vanderbilt, 2009; Bristowe et al., 2011; Goniakowski and Noguera, 2011; Stengel, 2011; Noguera and Goniakowski, 2013; Bristowe et al., 2014; Goniakowski and Noguera, 2014; Goniakowski and Noguera, 2016] is how to determine whether a particular surface is polar or non-polar and to quantify its instability by calculating its surface charge.

12.1 Calculating surface charge from $\rho$

It does not seem difficult to intuit the meaning of the surface areal charge density $\bsigma$ when one first encounters the concept. However, as soon as one tries to define it, in order to calculate it, difficulties become apparent.

The obvious way to calculate $\bsigma$ is simply to integrate the volumetric charge density ${\rho(\rvec)}$ from a point above the surface to a point far beneath it. Assuming that the surface is perpendicular to the $x-$axis, and that $\varrho(x)$ is the average of ${\rho(\rvec)}$ over the $yz-$plane at $x$, the obvious definition of the $yz-$averaged areal density of excess surface charge is

\begin{align*} \sigma_s(x_b) \equiv \int_{x_L}^{x_b}\varrho(x)\,\dd{x} \tag{57} \end{align*}

Fig. 11 illustrates why this definition fails. It depicts the ‘surface’ of a one-dimensional crystal whose microscale charge distribution is a semi-infinite periodic array of alternating positive and negative point charges of magnitude one. For this simple case the integral in Eq. (57) becomes the sum ${\sigma_s(x_b)=1-1+1-1+\cdots}$. Its value is either zero or one, depending on the choice of $x_b$, and it continues to vary between these values ad infinitum as $x_b$ increases. Therefore ${\sigma_s(x_b)}$ is a microscopic function of ${x_b}$ and, as a result, Eq. (57) fails as a definition of $\bsigma$.

Finnis presented an elegant solution to this problem in [Finnis, 1998], and a generalization of his result to amorphous materials is derived by a different route in Sec. 9 and Appendix J. I quote and explain the more general result below. I then quote Finnis’s result for crystals, which is simpler and easier to relate to the example depicted in Fig. 11.

12.2 Macroscale surface charge

The problem of how to express $\bsigma$ in terms of $\rho$ is easy to resolve once it is realised that $\bsigma$ only has meaning at the macroscale. Microscopically, surfaces and interfaces are ill-defined entities because their widths are indeterminate: in the vicinity of a surface, both structure and composition differ from the bulk, in general, but they gradually become more bulk-like with depth. This gradual relaxation means that there is no clear boundary separating surface-like material from bulk-like material.

However, as Fig. 11 illustrates, even if a material could be terminated abruptly at a plane and prevented from changing its local structure (bond lengths and angles) or composition, such that surface structure and composition were identical to the bulk, the concept of a surface charge density simply does not apply at the microscale: The microstructure is defined on a simply connected subset of $\realthree$. One can define an areal charge density on any plane (e.g., ${\sigma(y,z;x)\equiv \rho(x,y,z)\dd{x}}$), but no special surface plane exists.

However, as explained in Sec. 8, the spatial resolution $\prectheo$ at the macroscale is unavoidably finite and all points separated by microscopic distances coincide at the macroscale. As a result, the surface region of indeterminate width is contracted to zero width by the homogenization transformation. It becomes a two dimensional manifold.

The mesoscale average $\brho$ of $\rho$ at any microscale point whose macroscale image is on this manifold differs significantly (by more than ${\precRho/2}$), in general, from its value elsewhere. To understand why, consider the material depicted in Fig. 12. There exist planes (e.g., Plane 4) on which the planar average of $\rho$ does not vanish. It follows that the three dimensional mesoscale average $\brho$ at any point (not just points on the charged planes) only vanishes as a result of cancellation of positive and negative contributions whose displacements from the point have components normal to those planes. If all material from one side of such a plane is removed to create a surface, this balance is disrupted and $\brho$ becomes finite, in general, at any point within a distance ${\prectheo/2}$ of the plane. The areal charge density $\bsigma$ at a point on the surface manifold is simply the integral of ${\brho}$ over the point’s preimage under the homogenization transformation. Therefore it is the integral of $\brho$ along on an interval of width ${\prectheo}$ on an axis normal to the surface. The macroscopic volumetric charge density $\Rho$ is simply the average of $\brho$ on this interval.

It is logical, then, to define the surface charge as

\begin{align*} \bsigma = \int_{\mxl}^{\mxb}\Rho(\bx)\,\dbx . \tag{58} \end{align*}

This integral converges with respect to both of its limits because ${\Rho=0}$ in the bulk and in the vacuum above the surface. By substituting the definition of $\Rho$ as the mesoscale average $\brho$ of $\rho$ (Eq. (54)), it is shown in Sec. 9 that

\begin{align*} \bsigma \equiv \int_{x_L}^{x_b}\rho(x)\,\dd{x}-\bmonerho(x_b)\equiv\sigma_s(x_b) + \sigma_b(x_b), \tag{59} \end{align*}

where $x_b$ is any point deep below the surface and ${\bmonerho(x_b)\equiv -\sigma_b(x_b)}$ is defined as follows: A mesoscopic neighbourhood of $x_b$ is partitioned into a set of contiguous microscopic intervals ${\interval(\bar{x}_m,\Delta_m)}$ such that $x_b$ is at a boundary between two of these microintervals, and such that the integral of $\rho$ on each microinterval is zero. The second condition is always possible because $\rho$ fluctuates microscopically about zero. The dipole moment density of the $m^\text{th}$ interval is

\begin{align*} \monerho(\bar{x}_m,\Delta_m) \equiv \frac{1}{\Delta_m}\int_{-\Delta_m/2}^{\Delta_m/2}\rho(\bar{x}_m+u)\,u\,\dd{u}, \end{align*}

and $\bmonerho(x_b)$ is defined as the average of this quantity over all microintervals in the discrete and finite set that partitions the mesoscale neighbourhood of $x_b$. It is shown in Appendix J that the value of ${\bmonerho(x_b)}$ is the same for all sets of microintervals that satisfy the conditions stated above.

Although both ${\sigma_s(x_b)}$ and ${\sigma_b(x_b)}$ depend sensitively on $x_b$, their sum is independent of it. This is easy to see in the special case of a periodic bulk charge density, ${\rho(x+a)=\rho(x)}$. The points $x_b+ma$, where $m\in\mathbb{Z}$, can be chosen as the microinterval boundary points. All microintervals are identical, in this case, and Eq. (59) simplifies to Finnis’s result:

\begin{align*} \bsigma & = \int_{x_L}^{x_b}\rho(x)\,\dd{x}- \frac{1}{a} \int_0^a\rho(x_b+u)\,u\,\dd{u} \\ & = \sigma_s(x_b) - \mpp(x_b+a/2;a). \tag{60} \end{align*}

Note that the definitions of $\monerho$ and $\mpp$ are identical. I use $\mpp$ when it is useful to make clear that it is a dipole moment density. I use $\monerho$ when I favour consistency with Sec. 9 and with related quantities that will be introduced in Sec. 15.

Referring again to Fig. 11, and comparing the choices ${x_b=x_1}$ and ${x_b=x_2}$, we find that ${\sigma_s(x_1)=0}$ and ${\sigma_s(x_2)=1}$. If the minimum distance between positive and negative charges is denoted by $b$, then ${\sigma_b(x_1)= b/a}$ and ${\sigma_b(x_2)=-(a-b)/a}$. Therefore,

\begin{align*} \sigma_s(x_1)+\sigma_b(x_1)=\sigma_s(x_2)+\sigma_b(x_2)=\frac{b}{a}. \end{align*}

It may be illuminating to note that applying the homogenization transformation is a lot like taking a thermodynamic limit; and we can think of macroscopic quantities as thermodynamic quantities which, in a non-ergodic system, can only be defined on macroscopic length scales. Indeed, Finnis’s reasoning when deriving Eq. (60) differed slightly from the reasoning outlined above. He reasoned that one should average over an ensemble of terminating planes ($x_b$), in order to "reconcile the atomistic picture, in which [surface] excesses appear to oscillate on the atomic length scale as a function of the [surface] region size, with the thermodynamic picture." In the language that I have chosen to use here and in Sec. 8, he found the mesoscale average of the microscopic function $\sigma_s(x_b)$. The same result is found by substituting ${\Rho=\bar{\rho}}$ into Eq. (58) because, by changing the order of integration, the integral of a mesoscale average becomes the mesoscale average of an integral.

The subscripts of $\sigma_s$ and $\sigma_b$ are abbreviations of ‘surface’ and ‘bulk’, respectively. $\sigma_s$ includes all contributions from compositional differences between the surface and the bulk, including charged adsorbates, surpluses or deficits of electrons, charged impurities, and non-stoichiometry associated with reconstructions or coordination defects. On the other hand, $\sigma_b$, depends only on the charge density in the bulk and is independent of the surface composition. However, it is simplistic and wrong to view $\sigma_s$ as the contribution from extrinsic surface charges and $\sigma_b$ as the contribution from the bulk charge distribution. For example, it is always possible to choose $x_b$ such that ${\sigma_b(x_b)=0}$ and ${\bsigma=\sigma_s(x_b)}$. Therefore, as well as containing all extrinsic contributions to $\bsigma$, $\sigma_s$ can contain some, or all, of the contribution from the bulk. Choosing ${\sigma_b(x_b)=0}$ is equivalent to the "dipole-free unit cell" strategy used by Goniakowski et al. to deduce surface charge and stability [Goniakowski et al., 2008].

**Figure 12.** A lattice of point charges can notionally be cleaved along an infinite number of lattice planes. Here, four lattice planes are shown in cross section as dashed lines. Plane 4 intersects a lattice of positive point charges and is included to illustrate a point made in Sec. 12.2. Cleaving along Planes 1, 2, and 3, produce Surfaces 1, 2, and 3, respectively; Surfaces 1 and 2 are shown on the right. Deducing whether a surface is polar is as simple as constructing a unit cell of the crystal (shown in green), which has two lattice vectors that are parallel to the surface plane and one face at the surface, and calculating the dipole moment of this cell. The surface is polar if and only if the projection of the dipole moment onto the surface normal is non-zero. Surface 1 is non-polar but Surface 2 is polar. Surface 3 would be positively charged, despite Planes 2 and 3 being parallel and Surface 2 being negatively charged. Therefore, the orientation of the plane determines whether or not it is polar, but not the value of the surface charge density. One could also construct the dipole-free unit cell (navy boundary) and, if this cannot be constructed with one of its faces at the surface, the integral of the charge density from its uppermost face to the surface is the surface charge density, $\bsigma$.

12.3 Surface Stability

A pristine crystal surface is specified by the structure of the bulk crystal, the surface-plane orientation, and the surface termination. Consider the crystal depicted in Fig. 12 and the two examples given of surfaces of that crystal, which I’ll refer to as Surface 1 and Surface 2. These surfaces are defined by the planes (Plane 1 and Plane 2, respectively) at which one could imagine cleaving the perfect crystal. Each plane is defined by an orientation, which can be specified by the outward surface normal $\normal$, and by a position along an axis parallel to $\normal$. The importance of the relative displacements of parallel surface planes is illustrated by the fact that Surface 2 is negatively charged, whereas the surface created by cleaving at Plane 3 would be positively charged. In Fig. 12 I am considering an unphysical frozen surface where I have not allowed any relaxation from the bulk charge density Nevertheless, this unphysical surface suffices to allow the surface polarity to be quantified using Eq. (60).

If we could prepare the frozen surfaces depicted in Fig. 12, when we allowed them to relax in a vacuum they would relax and/or reconstruct. The surface might even melt. However, as long as the bulk crystal did not melt, or the surface didn’t banish ions or electrons from it (in the case of Surface 2 we would probably need to apply an electric field to prevent this), none of the structural change near the surface would have any impact on the surface charge, $\bsigma$, because the integral of the charge density between the surface and the crystalline bulk, which is the first term in Eq. (60), would be unchanged.

12.4 Interface charge

By treating an interface between two materials as a pair of adjoined surfaces it is straightforward to show that the interface charge is

\begin{align*} \bsigma = \int_{x_1}^{x_2} \rho(x)\,\dd{x} + \momone{x_1}-\momone{x_2} \end{align*}

where ${x_1}$ is an arbitrary point in the bulk of the material on one side of the interface and ${x_2>x_1}$ is an arbitrary point in the bulk of the material on the other side.

As a sanity check, let us imagine a plane perpendicular to ${\hat{x}}$ at position $x_b$ in the bulk of a mesocopically-uniform material. This plane can be viewed as an interface between two perfectly-aligned identical materials. The charge density at this imaginary interface is

\begin{align*} \bsigma = \int_{x_b-u}^{x_b+u}\rho(x)\,\dd{x} + \momone{x_b-u}-\momone{x_b+u} \end{align*}

where $u$ can have any value since all points on either side of $x_b$ are in the bulk. For the sake of brevity, I choose the limit ${u\to 0^+}$, which leads immediately to ${\bsigma=0}$, as should be the case. It is not too difficult to prove that $\bsigma=0$ for an arbitrary value of $u$ in a non-periodic system.

12.5 Consistency with the standard model of macroscale electrostatics

Standard treatments of macroscale electrostatics tend to distinguish between a free charge density ${\Rhofree=\div{\D}}$ and a bound charge density ${\Rhobound=-\div{\pp}}$, where

\begin{align*} \Rho = \epsilon_0\div{\E} = \div{\D} -\div{\pp} = \Rhobound+\Rhofree. \end{align*}

Substituting into Eq. (58) gives

\begin{align*} \bsigma & = \overbrace{\int_{\mx_L}^{\mxb} \Rhofree(\mx)\, \dbx}^{\displaystyle \bsigmafree} \overbrace{-\int_{\mx_L}^{\mxb} \div\pp \, \dbx}^{\displaystyle \bsigmabound} \\ & = \int_{\mx_L}^{\mxb} \Rhofree(\mx)\, \dbx + \pp\cdot\vec{\hat{n}} \end{align*}

and consistency with Eq. (60) requires that

\begin{align*} \int_{\mx_L}^{\mxb} \Rhofree(\mx) \dbx & + \pp\cdot\vec{\hat{n}} \\ = & \int_{x_L}^{x_b} \rho(x) \dd{x} + \mbp\left(x_b+a/2,a\right)\cdot\vec{\hat{n}}. \tag{61} \end{align*}

Both terms on the right hand side depend sensitively on $x_b$. However, because $x_b$ is arbitrary, both terms on the left hand side must be independent of it if physical meaning can be attributed to them independently of one other. One way to resolve this is by defining $\Rhofree$ and $\pp$ to be the mesoscale averages of $\rho$ and $\mbp$, respectively. That is, ${\Rhofree \equiv \bar{\rho} = \Rho}$ and ${\pp \equiv \bar{\mbp} = 0}$ (see Eq. (56)), which implies that ${\Rhobound = 0}$ and ${\bsigmabound \equiv \pp\cdot\vec{\hat{n}} =0}$. These definitions preserve consistency between the standard model of electrostatics in dielectrics and the apparently-reasonable definitions of $\Rho$ and $\bsigma$ presented herein; namely, $\Rho$ is the mesoscale average of $\rho$ and $\bsigma$ is its integral across a surface. However, achieving consistency in this way entails discarding several quantities from the standard model of electrostatics: $\D$, $\pp$, and $\Rhobound$ all vanish and ${\Rhofree}$ is simply ${\Rho}$. The macroscopic Maxwell equations are now identical in form to their microscopic counterparts because averaging commutes with differentiation; for example,

\begin{align*} \epsilon_0 \div\me(x) = \rho(x) \Rightarrow \epsilon_0 \div\E(\mx) = \Rho(\mx). \end{align*}

12.5.1 Quantization of $\pp$

It is important to consider carefully whether or not the less drastic option of keeping quantities $\pp$, $\D$, $\Rhobound$, and $\Rhofree$ within the macroscale theory is logical or viable. In this section I assume that $\pp$ remains an element of the theory and I show that consistency with the definitions of $\Rho$ (Eq. (54)) and $\bsigma$ (Eq. (58)) requires it to be quantized.

The fact that $\pp$ would be quantized for a classical crystal in the same way as it is quantized in the MTOP appears to have been appreciated from the beginning [Vanderbilt and King-Smith, 1993; Vanderbilt, 2018]. I explain it here for completeness, and also to emphasize it, because it is easy to misinterpret the term quantum of polarization as referring to something quantum mechanical.

As before I consider the surface of a pristine perfect crystal whose outward unit normal is $\normal$ and whose structure and composition have not been allowed to change after all the material on one side of the surface plane was removed. I choose the crystal’s primitive lattice vectors ${(\avec_1, \avec_2, \avec_3)}$ such that ${\avec_1\cdot\normal>0}$, ${\avec_2\cdot\normal=\avec_3\cdot\normal=0}$, and ${\normal\cdot\left(\avec_2\times\avec_3\right)=A_\Omega>0}$. The volume of the bulk crystal’s unit cell is ${\Omega\equiv\avec_1\cdot\left(\avec_2\times\avec_3\right)=\abs{\avec_1} A_\Omega}$. The surface, which is a two dimensional lattice, has primitive lattice vectors ${(\avec_2,\avec_3)}$ and the area of its primitive unit cell is ${A_\Omega}$. Given the surface normal $\normal$, this choice of primitive unit cell of the bulk crystal allows all surface terminations to be identified by a single parameter $\alpha$, which is the position along $\normal$ at which the uppermost primitive cell is sliced to form the surface. For example, surfaces formed by cleaving at Planes 2 and 3 of Fig. 12 differ only by their values of $\alpha$. In this simple case the value of $\alpha$ determines only whether the uppermost plane is a plane of cations or a plane of anions. In more general cases the electron density would also be divided; however, it would be unphysical to remove fractions of electrons by removing all density above the termination plane, so I assume that the integral of the density that remains in the uppermost cell is rounded up to a whole number. How this density is distributed has no bearing on the arguments to follow.

The excess bound charge at the surface of the crystal is ${\bsigmabound=\pp_\perp\equiv\pp\cdot\normal}$. Within the standard model of electrostatics $\pp$ is a bulk quantity and so it must be independent of surface termination $\alpha$. Therefore $\bsigmabound$ must be the same for all surfaces whose outward normal is $\normal$. However, as discussed in Sec. 12.3, and as Fig. 12 illustrates, ${\bsigma=\bsigmabound+\bsigmafree}$ is not the same for all values of $\alpha$. One could choose to include all of the $\alpha$-dependence of $\bsigma$ in $\bsigmafree$, which would leave $\bsigmabound$ independent of surface termination. However, this is not the approach taken within the MTOP [Vanderbilt and King-Smith, 1993; Stengel, 2011]. The MTOP assumes the standard convention that $\Rhofree$ and $\bsigmafree$ only contain contributions from charges that are not intrinsic to the material [Ashcroft and Mermin, 1976; Jackson, 1998]. As a consequence of preserving this old convention, $\pp$ must be quantized [Vanderbilt and King-Smith, 1993; Vanderbilt, 2018]. I now prove this.

There are no extrinsic charges in the idealized surfaces constructed; therefore ${\bsigmafree=0}$ and ${\bsigma=\bsigmabound=\pp_\perp}$. Now, because $\bsigma$ is known and single-valued, and because it can be changed to the value it would have for any other value of $\alpha$ by adding/removing the same numbers and types of particles (nuclei and electrons) to/from each unit cell of the surface lattice, $\bsigmabound$ must be multivalued. Its set of values must be the set of values of $\bsigma$ for every possible choice of surface termination, $\alpha$. These values differ by integer multiples of ${e/A_\Omega}$. Therefore, $\pp_\perp$ is quantized such that if ${\Delta\pp}$ is the difference between two values of $\pp$ that are consistent with Eq. (54) and Eq. (58), then

\begin{align*} \Delta\pp_\perp \equiv \Delta\pp\cdot\normal &= \frac{m_1\, e}{\abs{\avec_2\times\avec_3}} = \frac{m_1\, e}{\avec_1\cdot\left(\avec_2\times\avec_3\right)}\, \avec_1\cdot\normal \\ \implies \Delta\pp & = \frac{m_1\,e}{\Omega}\,\avec_1 + \Delta\pp_{\parallel, 2}\,\avec_2 + \Delta\pp_{\parallel, 3}\,\avec_3 \end{align*}

where ${\Delta\pp_\parallel \equiv \Delta\pp_{\parallel, 2}\,\avec_2 + \Delta\pp_{\parallel, 3}\,\avec_3 =\Delta\pp - \Delta\pp_\perp\,\normal}$, and $m_1$ is an integer. By considering surfaces perpendicular to ${\avec_3\times\avec_1}$ and ${\avec_1\times\avec_{2}}$, the same logic would lead me to the following general expression for the difference between any two values of polarization that are consistent with Eq. (54) and Eq. (58).

\begin{align*} \Delta\pp = \frac{e\A}{\Omega} \tag{62} \end{align*}

where ${\A = m_1\,\avec_1+m_2\,\avec_2+m_3\,\avec_3}$ is an arbitrary lattice vector and ${m_1, m_2, m_3 \in\mathbb{Z}}$. This is identical to the quantization of $\pp$ deduced within the MTOP, but this derivation makes clear that the quantization of $\pp$ is not a quantum mechanical quantization because I have not invoked quantum mechanics to deduce it. It is a consequence of adopting the apparently arbitrary and unnecessary convention that, in the absence of extrinsic charges, $\bsigmabound$ is the total surface charge density, and of preserving the macroscale theory’s internal consistency while shoehorning $\pp$ into it.

12.6 Mapping to a set of localized charge packets

In this section I present a result that is pivotal for understanding the relationship between this work, which is founded on a systematic approach to structure homogenization, and the MTOP’s definition of polarization current, which is founded on quantum mechanical perturbation theory.

Let us express the charge density as the sum, ${\rho(x)=\sum_i\rho_i(x)}$, of a set of localized charge packets, ${\{\rho_i\}}$, where each $\rho_i$ is either nonpositive or nonnegative. The total charge in the $i^\text{th}$ packet is $q_i$ and its center of charge is $x_i$. That is,

\begin{align*} q_i \equiv \int_{-\infty}^\infty \rho_i(x) \dd{x}, \quad x_i \equiv \frac{1}{q_i}\int_{-\infty}^\infty \rho_i(x)\, x \,\dd{x}. \end{align*}

I assume that each ${\rho_i}$ can be chosen such that it is negligible outside interval ${\interval(x_i,\prectheo)}$. For a system of nuclei and electrons whose charge density is given by Eq. (4), the nonnegative packets are ${\rho_i(x)\equiv Z_i\,e\,\delta(x-X_i)}$ and the nonpositive packets are ${\rho_i(x)\equiv-e\,n_i(x)}$, where ${n(x)\equiv\sum_{i}n_i(x)}$ is the electron number density partitioned into a set of packets, ${\{n_i\}}$.

The localization transformation ${\rho_i(x)\to q_i\,\delta(x-x_i)}$ conserves charge and preserves $\rho_i$’s center of charge. Therefore, the transformation of ${\rho}$ into the discrete distribution of point charges ${\rho^q(x)\equiv\sum_i q_i\, \delta(x-x_i)}$ is an isotropic spatial redistribution of charge. By ‘isotropic’ I mean that it does not change the center of charge of either $\rho$ or ${\trho\equiv\sum_{i\in\iset}\rho_i}$, where $\iset$ is any subset of the set of packet indices. The equitable movement of charge in both directions cannot change the macrostructure if charge is only moved across distances smaller than ${\prectheo=\dbx}$. Therefore, the macroscale counterpart ${\Rho^q}$ of ${\rho^q}$ cannot differ from $\Rho$. From this fact, and from Eq. (59), it follows that

\begin{align*} \bsigma & = \int_{\mxl}^{\mxb}\Rho^q(\bx)\,\dbx = \sum_{i:x_i<x_b} q_i - \bmonerhoq(x_b) \tag{63} \end{align*}

Note that the derivation of Eq. (48) in Sec. 9 assumed that space could be partitioned into microintervals whose net charges were equal. This is always possible for a continuous charge density, but it is not possible for the $\rho^q$’s resulting from every possible partitioning of $\rho$ into packets because the integral of $\rho^q$ on each microinterval is a sum of point charges. If the magnitudes of these charges are irregular it is not possible, in general, to partition space such that each interval’s average charge density is precisely the same. It seems likely that a more general derivation, which does not require each microinterval to have the same average charge density, is discoverable.

In the bulk of a crystal with periodicity $a$, the charge packets can be chosen such that, for any packet $\rho_i$, whose center is $x_i$, there are identical packets with centers at ${x_i+m a}$ for all ${m\in\mathbb{Z}}$. In this case, the surface charge is

\begin{align*} \bsigma = \sum_{i:x_i<x_b} q_i - \frac{1}{a}\sum_{i:x_i\in(x_b,x_b+a)} q_i\,x_i, \tag{64} \end{align*}

where I am assuming that ${\sum_{i:x_i\in(x_b,x_b+a)} q_i=0}$. The time derivative of Eq. (64) is the current $\J$. If the conduction current vanishes, it is the polarization current $\Jconv$.