A theory for the stabilization of polar crystal surfaces by a liquid environment

Crystal morphology is a determining factor in the mechanical, rheological, and catalytic properties of ionic crystals.¹ As such, there is considerable effort to establish simple and inexpensive routes to control which facets of a crystal are exposed at its surfaces. In this article, we will concern ourselves exclusively with polar crystal faces, which we describe in more detail in Sec. II A. Examples of polar crystal surfaces include halite (1 1 1), wurtzite (0 0 0 1), and zinc blende (1 1 1). Polar surfaces have been studied extensively from a solid state physics perspective^2–4 owing to their use in microelectronic devices and the now well-established fabrication techniques of thin metal oxide films. Polar surfaces are also important for catalysis. For example, the {1 1 1} facets of MgO exhibit ultrahigh catalytic activity⁵ for Claisen–Schmidt condensation of benzaldehyde and acetophone, and the (0 0 0 1) face of ZnO is active in methanol⁶ decomposition.

What makes polar surfaces both challenging and interesting is that basic electrostatic arguments show that their surface energy diverges with increasing crystal thickness along a polar crystallographic direction;^2,3,7 this is the so-called “polar catastrophe.” The fact that polar crystal surfaces are observed, thus, implies a “polarity compensation mechanism” that overcomes the underlying instability. Known polarity compensation mechanisms are as follows:^3,4 (i) charge transfer via nonstoichiometric reconstruction, which amounts to an effective transfer of ions from one side of the crystal to the other; (ii) electronic reconstruction, where interfacial charges are modified by partial filling of electronic interface states; and (iii) adsorption of foreign atoms or ions that provide an appropriate amount of compensating interfacial charge. The theoretical rationalization of the polar catastrophe from energetic considerations of a crystal surface in contact with a low-density vapor phase lends itself most naturally to describing mechanisms (i) and (ii).

Polarity compensation through adsorption of ions from a surrounding liquid environment can be considered an extreme example of mechanism (iii), and it provides a potentially simple and cost-effective route to controlled exposure of polar crystal facets. For example, using the molten salt synthesis route, Xu et al.⁸ synthesized particles of ZnO, MgO, and Co₃O₄ that exclusively expose polar facets, and subsequent studies from Susman et al. suggested that K⁺ and Cl⁻ ions potentially stabilize the {1 1 1} facets of MgO,⁹ while a complex interplay between the molten salt and resulting crystal morphology is evident for NiO.¹⁰ To understand the processes that underlie polarity compensation by ion adsorption from a surrounding liquid, it will be useful to have a theory that treats the environment on an equal footing to the crystal. The purpose of this article is to develop such a theoretical framework.

A second motivation is to reiterate the central message of Tasker’s seminal article² on the theory of polar surfaces: to demonstrate that the polar catastrophe has a true physical origin and is not the result of lattice summation techniques. Such a (re)clarification is necessary as the use of molecular simulations—and consequently, lattice summation techniques—has exploded since Tasker’s original 1979 publication, and the technical nature of treating long-range electrostatic interactions under periodic boundary conditions can obfuscate physical intuition. Similar to Tasker, then, the theory presented here will not rely on periodic boundary conditions. Nonetheless, we will see implications for simulations that do employ periodic boundary conditions, and we make two observations that are useful from a practical viewpoint when modeling, with a slab geometry, the surfaces of polar crystals in contact with a conducting medium, such as an electrolyte:

1. Tinfoil Ewald approaches, briefly described in Sec. II B, are appropriate if we are interested in genuinely thin polar crystals.

2. An electric displacement field emerges as a boundary condition.

The second point clarifies the role of the electric displacement field imposed in previous simulation studies.^11–13 Specifically, the current work emphasizes that it is conceptually incorrect to associate this electric displacement field with the removal of spurious interactions between periodic replicas of the simulation cell and places the observation made in Ref. 13—that standard “slab corrections” lead to qualitatively incorrect results—on a firmer theoretical footing.

The rest of this article is organized as follows: In Sec. II A, we derive an expression for the equilibrium adsorbed surface charge density at halite (1 1 1) in contact with an electrolyte solution as a function of crystal thickness. We also show how the polar catastrophe manifests as a diverging free energy cost for fluctuations away from equilibrium as the crystal thickness increases. In Sec. II B, we show that the theoretical predictions are consistent with results from molecular dynamics simulations that employ a tinfoil Ewald approach to treat electrostatic interactions. In Sec. III, we discuss the emergence of the electric displacement field in the slab geometry and also outline a justification for the use of tinfoil Ewald approaches to model genuinely thin polar crystals. We summarize our findings in Sec. IV, before providing a brief summary of the simulation methods in Sec. V. In the  Appendix, we also provide a more formal justification for the use of tinfoil Ewald approaches.

A schematic of the system we are interested in is shown in Fig. 1(a). Here, a crystal exposing polar facets is immersed in an electrolyte environment. We assume that the crystal assumes a bulk truncated structure, i.e., nonstoichiometric reconstruction [mechanism (i)] has not taken place. Moreover, we assume an ionic model such that electronic reconstruction [mechanism (ii)] is not possible. With these approximations, the exposed crystal faces have surface charge density ±σ₀. We now suppose that ions from solution adsorb to the polar crystal faces with surface charge density ∓σ⁽ⁿ⁾. The reason for the “(n)” superscript will become clear shortly. To investigate these surfaces, we imagine taking a cut of the system far from the edges of the crystal, as indicated by the dotted box in Fig. 1(a), and ignore any edge effects. This amounts to assuming that the planes of charge span the entire plane orthogonal to the surface normal. This “slab geometry,” which is shown in greater detail in Fig. 1(b), is similar to the textbook parallel plate capacitor model,¹⁴ a fact we draw upon heavily in what follows. Despite relying on the parallel plate capacitor model extensively, we should acknowledge that neglecting edge effects in the slab geometry is not an innocent omission; in Sec. III, we will see that doing so ultimately gives rise to a boundary condition that is absent in the system of interest [Fig. 1(a)].

Also given in Fig. 1(b) is a more detailed description of the structure of the polar crystal. For simplicity, we consider halite (1 1 1)—relevant to, e.g., MgO, NiO, and NaCl—in which n + 1 crystal planes are separated by a distance R, where n is an odd integer, though results can be generalized to other geometries. Later, we will also perform “simulation experiments” in which we change the separation between the central planes in the crystal to R_c ≠ R. In the slab geometry, the underlying average microscopic charge density ρ_c(z) only varies in the direction normal to the surface, which we indicate by z. The total adsorbed surface charge density is then determined by

where “int−” and “int+” indicate that the integrations are performed over interfacial regions corresponding to the negatively and positively charged crystal surfaces, respectively.¹⁵ While Eq. (1) provides a direct means to obtain σ⁽ⁿ⁾ from average microscopic properties, deciding where to locate the planes of adsorbed surface charge density is less straightforward. Nonetheless, it seems reasonable to place the plane with charge density σ⁽ⁿ⁾ a distance ℓ₋ from the crystal’s negatively charged surface and that with −σ⁽ⁿ⁾ a distance ℓ₊ from the crystal’s positively charged surface, where ℓ₋ and ℓ₊ both have a clear microscopic interpretation, inasmuch as they can be related to features of ρ_c(z). In the following, we will, in fact, find an optimal apparent length scale entering the continuum model that has no such clear microscopic interpretation.

Let us denote the free energy per unit area of the system as f(σ⁽ⁿ⁾). Despite discussing the system in terms of uniform planes of charge, we have already acknowledged that the underlying system comprises atomic and molecular entities. It is natural, then, to partition f into “capacitive” and “noncapacitive” contributions,

The capacitive contribution, u_cap, captures the energy stored in the electric fields assuming that the system comprises uniform planes of charge, as shown in Fig. 1(b). The noncapacitive term, f_nc, captures contributions from everything else, e.g., non-electrostatic interactions, electrostatic interactions not captured by the simple capacitor model, solvent effects, and any entropic contributions.

To find u_cap, we view the crystal in the slab geometry as (n + 1)/2 parallel plate capacitors, between which the electric field is E = 4π(σ⁽ⁿ⁾ − σ₀), along with (n − 1)/2 capacitors between which the electric field is E = 4πσ⁽ⁿ⁾. For both sets of capacitors, the separation between plates is R (for the moment, we consider R_c = R). In addition, the electric field between an adsorbed plane and the surface of the crystal is also E = 4πσ⁽ⁿ⁾. Recalling that the energy density of an electric field is |E|²/8π, it follows that

where 2ℓ = l₊ + l₋. (Throughout this formulation, we work in a unit system in which 4πϵ₀ = 1, where ϵ₀ is the permittivity of free space.) At equilibrium, σ⁽ⁿ⁾ assumes a value that minimizes f,

where the prime indicates a partial derivative with respect to σ⁽ⁿ⁾. Equation (4) is a central result of this article. While noncapacitive contributions are irrelevant as n → ∞, their effects become increasingly important as the thickness of the crystal decreases.

Although formally exact, analysis of Eq. (4) is complicated by the presence of the noncapacitive contributions. To make exploratory progress, we simply postulate that

The equilibrium adsorbed surface charge density then reads

Within the phenomenological model specified by Eq. (5), one can view the effects of noncapacitive contributions as modifying the apparent length scale associated with ion adsorption, from ℓ to ℓ + a_nc, in the continuum representation of the system. A similar expression for σ^(n,eq) was derived previously by Hu¹⁶ for the slab geometry under periodic boundary conditions, within the framework of symmetry-preserving mean field theory.^17,18 Within that framework, an effective length scale is rationalized by considering the slowly varying components of ρ_c, whose structural features need not coincide with those of the full average microscopic charge density. We will postpone further discussion of the effective length scale ℓ + a_nc to Sec. II B and the  Appendix. Importantly, the physical interpretation of Eq. (6) is somewhat different from that of Ref. 16, where σ⁽ⁿ⁾ = σ^(n,eq) was attributed to the presence of a spurious electric field arising from the use of tinfoil Ewald sums. In contrast, in this study, Eq. (6) has been derived for a nonperiodic slab geometry; the driving force for ion adsorption is a significant reduction in the capacitive energy stored in the system upon ion adsorption.

To understand the origin of the polar catastrophe from this free energy perspective, we consider the reversible work required to change the adsorbed surface charge density by an amount δσ from its equilibrium value,

Neglecting terms higher than second order in δσ and combining with Eqs. (2), (3), and (5) give

Equation (7) is another key result of this study; it is analogous to the familiar diverging surface energy of an unreconstructed polar crystal in contact with vacuum.^2,3 We see that the reversible work required to change the adsorbed surface charge density from σ^(n,eq) by any finite amount diverges linearly with the crystal thickness nR.

Despite the clear analogy to the well-known diverging surface energy of polar crystal surfaces in vacuum, Eq. (7) contains more information. This can be seen more clearly by considering the behavior of the fluctuations of the total charge Q within a probe area A⁽ⁿ⁾ at the surface of the crystal. At equilibrium, the total adsorbed charge within such an area is Q^(n,eq) = A⁽ⁿ⁾σ^(n,eq). Equation (7) states that the area required to observe δQ = Q − Q^(n,eq) ≈ e with appreciable probability, i.e., βA⁽ⁿ⁾[f(σ^(n,eq) + δσ) − f(σ^(n,eq))] ≈ 1, is

where e is the unit of elementary charge and β = 1/k_BT (where k_B is Boltzmann’s constant and T is the temperature). As the thickness n of the crystal increases, ever increasing probe areas are required to observe appreciable fluctuations in the adsorbed surface charge. The theory presented here not only reveals the strong constraint placed on the average adsorbed surface charge density at polar crystal surfaces, but also that fluctuations at the interface are strongly suppressed.

Equation (6) provides an expression for the equilibrium adsorbed surface charge density at polar halite (1 1 1) surfaces in contact with an electrolyte solution, assuming an ionic model for the crystal. As such, Eq. (6) lends itself naturally to validation with molecular simulations. However, two comments are in order. The first concerns the treatment of long-range electrostatic interactions under periodic boundary conditions, for which lattice summation techniques are often employed. For Ewald summation, which underlies many of the most popular techniques used in molecular simulations, the expression for the Coulomb energy for a set of charges {q_i} with positions {r_i} reads

In Eq. (9), b is a lattice vector, and the prime indicates that the term i = j is omitted when b = 0. The volume of the simulation cell is Ω, and $ρ̃(k)$ is the Fourier transform of the charge density ρ(r) = ∑_iq_iδ(r − r_i). The potentials ϕ_SR(r) and ϕ_LR(r) are defined by a splitting of the Coulomb potential into short- and long-ranged contributions, 1/|r| = ϕ_SR(r) + ϕ_LR(r), and $ϕ̃LR$ is the Fourier transform of the long-ranged part. Most important for the current discussion is the final “surface term,” which depends upon the total dipole moment of the simulation cell, M, and the depolarization tensor, J. This surface term is determined by the specified summation order of the lattice sum, and for a spherical summation order, it is equal to 2π|M|²/(2ϵ′ + 1), where ϵ′ is the dielectric constant of a surrounding medium “at infinity.” Tinfoil Ewald approaches take this surrounding medium to be a perfect conductor, ϵ′ = ∞, and amounts to ignoring the surface term. The above is not intended as a detailed discussion of Ewald summation, for which there is extensive literature.^19–29 The reader is referred to Ref. 30 for a particularly clear discussion on the topic.

For systems in the slab geometry under periodic boundary conditions employing tinfoil Ewald, where the slab is surrounded on either side by vacuum, it has long been recognized that if a net dipole M_z exists along z, then there is a finite electric field in the vacuum region. It is, therefore, common to employ “slab corrections” to mitigate possible spurious interactions between periodic images. In classical simulations of liquid–solid interfaces, the standard choice is the method of Yeh and Berkowitz,³¹ which includes a surface term $2πMz2/Ω$⁠. A similar expression exists in the surface science community (which typically uses periodic density functional theory calculations),^32,33 where it is known as the “dipole correction.” In Ref. 13, however, it was demonstrated that typical slab correction schemes impose vanishing adsorbed surface charge density; more generally, it was found that an electric displacement field directly determines the adsorbed surface charge, irrespective of the slab’s thickness. Such a scenario is clearly incompatible with σ^(n,eq) predicted by Eq. (6). In contrast, simulations performed with tinfoil Ewald approaches were found to give adsorbed surface charge densities that increase with crystal thickness. We will, therefore, compare the predictions of the theory presented above to simulations that use a tinfoil Ewald approach; a theoretical justification for doing so will be given in Sec. III, where it will also become clear why standard slab corrections fail.

The second comment concerns the parameters that determine σ^(n,eq). Here, we will consider crystals where the positions of the ions have been clamped and whose charges do not fluctuate; σ₀ and R are then straightforward input parameters. The extent to which agreement between the theory and simulation is observed is an indicator of how faithfully the theory presented in this article describes the ionic model of polar surfaces. Moreover, we will consider the specific example of NaCl (1 1 1) in contact with an aqueous NaCl solution. We then anticipate that ions in direct contact with the crystal will form a plane separated from the crystal by a distance ℓ ≈ R. The only remaining unknown parameter is a_nc, the length scale associated with noncapacitive contributions.

In a first validation step, we compare σ^(n,eq) obtained from molecular simulations of a concentrated aqueous NaCl solution in contact with halite (1 1 1), where a_nc is simply treated as a free fitting parameter. As seen in Fig. 2(a), Eq. (6) with ℓ = R = 1.628 Å and a_nc = −0.568 Å captures the trend seen in the molecular simulations very well. To determine the significance of noncapacitive contributions, it is instructive to compare the effective length scale ℓ + a_nc = 1.06 Å in the continuum model to real length scales in the system. To that end, Fig. 2(b) shows the number density profile of Na⁺ ions at halite (1 1 $1̄$⁠) with n = 5. As expected, the plane separated from halite (1 1 $1̄$⁠) by R coincides with the plane of Na⁺ directly adsorbed to the surface. In contrast, the plane at ℓ + a_nc from (1 1 $1̄$⁠) is located in a region that cannot be readily associated with ion adsorption; it does not appear that a_nc simply captures effects due to thermal fluctuations.

The catch-all nature of f_nc, and the ad hoc assertion of its quadratic form [Eq. (5)], makes direct physical interpretation of the value of a_nc challenging. Recent work suggests that dielectric boundaries for systems with water as a solvent are closely associated with water’s hydrogen density.^34,35 Also plotted in Fig. 2(b), therefore, is the number density profile of water’s hydrogen atoms. While some hydrogen density is found marginally closer to halite (1 1 $1̄$⁠) than that of the Na⁺ ions, it is difficult to draw any firm conclusions on a possible relationship between a_nc and the solvent. More likely is that noncapacitive contributions account for the rather drastic approximation of collapsing the entire adsorbed surface charge density into a single plane. As can be seen in Fig. 2(c), ρ_c(z) exhibits a pronounced structure beyond the first peak at ℓ = R. In Sec. III, we will see that the effective length scale ℓ + a_nc in the continuum model is that which equalizes the electrostatic potential in the solution on either side of the crystal, a condition that is automatically satisfied by molecular simulations of a crystal slab immersed in an electrolyte solution, in which tinfoil Ewald approaches are used (see the  Appendix).

Despite the lack of a simple physical interpretation, the importance of f_nc can be readily seen by plotting Eq. (6) with a_nc = 0, as shown in Fig. 2(a). Clearly, a model that only considers capacitive contributions in which the system is approximated as planes of uniform charge density poorly describes the simulation data, with predicted adsorbed surface charge densities much lower than that observed.

A possible cause for concern is that simply treating a_nc as a free fitting parameter amounts to nothing more than a convenient fix. Challenges in interpreting the value of a_nc notwithstanding, it is reasonable to assume that noncapacitive contributions largely arise from effects that are localized to the solution and interface regions. While changes to the structure of the crystal in regions far from the interface will affect u_cap, we do not anticipate significant impact on the coefficient a_nc. In a second validation step, we, therefore, perform a “simulation experiment” in which the separation between the crystal’s central planes is varied (i.e., R_c ≠ R). The manner in which σ^(n,eq) is affected by R_c ≠ R depends on whether the central planes bound a region with E = 4πσ⁽ⁿ⁾ [(n + 1)/2 is even] or E = 4π(σ⁽ⁿ⁾ − σ₀) [(n + 1)/2 is odd]. Following the same approach leading to Eq. (6) gives

Results obtained from simulations with R_c = 3R/4, R/2, and R/4 are presented in Fig. 3, along with the theoretical predictions described by Eqs. (10a) and (10b), in which no further fitting has been performed. That is, the same value a_nc = −0.568 Å has been used. The good agreement between the theoretical predictions and the simulation data lends support to the notion that fitting a_nc captures genuine noncapacitive effects that are localized to the solution and interface regions.

Results presented so far demonstrate consistency between the theory presented in Sec. II A for the equilibrium adsorbed surface charge density at halite (1 1 1), and molecular simulations that use a tinfoil Ewald approach. In Sec. III, we will discuss why it is appropriate to compare σ^(n,eq) for the system of interest [Fig. 1(a)] to simulations that employ the slab geometry under periodic boundary conditions. In doing so, we will also shed light on the role of the electric displacement field in the slab geometry.^{12,13,36–39}

Recall the procedure employed in Sec. II A to investigate the polar surfaces of a crystal in solution. First, a particle of finite size, with polar surfaces exposed, was immersed in an electrolyte, and ions from solution were imagined to adsorb to its surfaces. Then, we constructed the slab geometry by taking a cut of the system far from the edges of the crystal and considered the free energy per unit area. Why did we adopt this procedure, rather than starting immediately with the slab geometry? Consider again the system of interest, which, to avoid notational clutter, is shown again in Fig. 4(a). As the particle has finite size, the ions adsorbed to the positive and negative crystal faces originate from the same pool of ions in the surrounding electrolyte solution. The remaining electrolyte solution is electroneutral. As the electrolyte solution is conducting, the total electric field in its interior vanishes, E_e = 0. Now, consider the vapor region. We will assume that the vapor pressure is sufficiently low that this region can be approximated as vacuum. In the absence of an external field, the total electric field there too vanishes, E_v = 0. It then follows that the surface charge density at the electrolyte–vapor interface is zero, i.e., $σev=r̂⋅(Ev−Ee)=0$⁠, where $r̂$ is the outward normal to the electrolyte surface. This scenario is possible with the procedure employed in Sec. II A, as the electrolyte solution remains uncharged. Moreover, there is a clear physical pathway by which the ions have adsorbed from the solution to the surface.

Imagine that we instead start directly with the slab geometry. As the crystal spans the entire plane orthogonal to the surface normal, there are now two distinct regions of electrolyte. Upon ion adsorption, electroneutrality must be preserved in each region separately (the ions cannot pass through the crystal). As the electrolyte is a conductor, the remaining charge in solution must be located at the electrolyte–vapor boundary, |σ_ev| = σ⁽ⁿ⁾, as indicated in Fig. 4(b). As before, E_e = 0, but the electrolyte is now polarized, with uniform polarization $ẑ⋅Pe=σ(n)$⁠. The electric displacement field in the electrolyte then satisfies

In the vapor, we assume that atomic density is sufficiently low that P_v = 0. Matching the electric displacement field across the boundaries leads us to conclude that

This result for E_v in the slab geometry is most severe. Integrating the electrostatic energy density |E_v|²/8π over the macroscopic volume occupied by the vapor phase indicates a prohibitively expensive energy cost for σ⁽ⁿ⁾ ≠ 0. Another route to arrive at the conclusion σ⁽ⁿ⁾ = 0 when starting in the slab geometry is that, as the electric field external to the parallel plate capacitors vanishes, there is no driving force for ion adsorption.

The slab geometry, despite widespread use in both theoretical and computer simulation studies, does not correspond to a physical reality. Consider the electrostatic potential in the regions of the electrolyte indicated by the dotted and dotted-dashed circles in Fig. 4(a), which we denote V_a and V_b, respectively. As the electrolyte is conducting, it follows immediately that V_a = V_b. In the slab geometry, the crystal provides an insulating layer between the two distinct regions of electrolyte. For σ⁽ⁿ⁾ = 0, the electrostatic potential difference across the crystal is

Even for n = 1, this electrostatic potential difference is substantial: V_a − V_b ≈ 20 V. In contrast, for σ⁽ⁿ⁾ = σ^(n,eq) [Eq. (6)],

where, as discussed in the  Appendix, we have treated the adsorbed planes as if they are separated from the crystal’s surfaces by the effective distance ℓ + a_nc. The slab geometry fundamentally breaks the system: for σ⁽ⁿ⁾ = 0, E_v = 0 is accurately captured, but there exists an erroneous potential difference across the crystal; for σ⁽ⁿ⁾ = σ^(n,eq), the potential difference across the crystal is correct, but at the cost of introducing a finite field in the vapor phase, $ẑ⋅Ev=4πσ(n,eq)$⁠. This is a consequence of the finite electric displacement field in the electrolyte, $ẑ⋅De=4πσ(n,eq)$ [Eq. (11)], that emerges in the slab geometry.

What are the implications for molecular simulations? Let us consider cases where periodic boundary conditions are used in all three Cartesian directions, restricting ourselves to cases where periodic replicas of the crystal slab sandwich the electrolyte solution such that there is no vapor region, as shown in Fig. 4(c). Similar to the nonperiodic slab geometry, Fig. 4(b), the electrolyte is polarized, $ẑ⋅Pe=σ(n)$⁠, and there is a finite electric displacement field in the electrolyte, $ẑ⋅De=4πσ(n)$⁠. For simulations that employ tinfoil Ewald approaches to compute electrostatic interactions, there is no electrostatic potential difference across the simulation cell. As E_e = 0, it then follows that V_a − V_b = 0. (In the  Appendix, it is shown explicitly that σ^(n,eq) is the adsorbed surface charge density that enforces E_e = 0 when using tinfoil Ewald approaches.) By not including a vapor region between periodic images, the electrostatic energy is well behaved; as D_e · E_e = 0, we can extend the electrolyte region between periodic replicas without incurring any cost in electrostatic energy. Simulations using tinfoil Ewald sums, thus, appear to give a faithful representation of polar surfaces in solution—in the sense of taking a cut of the system of interest as shown in Fig. 1(a), albeit with a finite polarization and electric displacement field in the electrolyte, the effects of which, however, seem largely benign.

Implicit in the preceding discussion is that the polarization is treated in an itinerant fashion—the ions are included in its definition—which has been shown^40,41 to satisfy key statistical mechanical properties for electrolyte solutions, such as the Stillinger–Lovett sum rules.^42,43 Any discussion on the multivaluedness of polarization under periodic boundary conditions has also been neglected, along with implications of whether the crystal or the electrolyte straddles the boundaries of the simulation cell. The reader is referred to Refs. 12, 36, and 44 for a discussion of these issues in the context of molecular dynamics simulations.

As already mentioned, Ref. 13 demonstrated that imposing an electric displacement field directly determines the adsorbed surface charge density. Theoretical discussion in that work, however, was brief, and what little there was relied heavily upon the framework provided by the “finite field approach” of Zhang and Sprik.^36,39,45 The picture to emerge from the finite field approach is that of a set of “virtual electrodes” that either controls the electrostatic potential difference across the simulation cell (“constant E”) or controls the electric displacement field at the cell boundaries (“constant D”). As the finite field approach was formulated in the context of periodic boundary conditions, it can be challenging to disentangle its deeper significance from issues concerning conditional convergence of lattice sums and unwanted interactions between periodic images. Moreover, the picture of virtual electrodes can give the impression that one is forcing the system adopt a particular adsorbed surface charge density with an arbitrary field. In this article, the status of the electric displacement field as a control variable has been clarified: it should be viewed as matching an electrostatic boundary condition that emerges in the slab geometry—even in the absence of periodic boundary conditions—rather than an ad hoc external electric field that is applied to the system. In Ref. 13, the surfaces of macroscopically thick crystals in contact with an electrolyte solution were effectively simulated by imposing $ẑ⋅De=4πσ(∞,eq)$ using the finite field approach. Similar control of the adsorbed surface charge with an electric displacement field has also been observed to work away from equilibrium.⁴⁶ 

To end this discussion, it is worth considering standard slab correction schemes again,^31–33 which typically introduce a vacuum region between periodic replicas along z and remove the electric field in the vacuum region by enforcing $ẑ⋅Dv=0$⁠. While such a line of attack seems reasonable from the perspective of removing interactions between periodic images, it presupposes that the true physical scenario corresponds to an isolated slab in contact with the liquid [e.g., Fig. 4(b)]. Slabs that completely span the plane orthogonal to the surface normal, however, do not exist; this fact cannot be ignored when modeling the surfaces of polar crystals.

A theoretical framework to describe polarity compensation arising from the solution environment has been presented in which the free energy of the system has been separated into “capacitive” and “noncapacitive” contributions. For crystals that are thin along a polar crystallographic direction, noncapacitive contributions were found to play a significant role. As the thickness of the crystal increases, the capacitive contribution dominates, and the adsorbed surface charge density at equilibrium agrees with known expressions for the required charge imbalance at polar surfaces.^2,3

The theory presented here only attempts to capture the net adsorbed charge at each polar crystal surface and says nothing about the complex structural details of the solution near the interface. A major source of the noncapacitive contributions likely arises from the drastic approximation of assuming that all adsorbed charge lies in a single plane. In addition, the theory assumes an ionic model for the crystal. While the severity of such an approximation will depend on the system under investigation,³ it may nonetheless provide a useful starting point for more sophisticated theoretical approaches.

Like previous theoretical treatments for polar crystal surfaces in contact with vacuum,^2,3 the presented theory relies heavily on modeling the crystal as a stack of parallel plate capacitors. It was shown that in such a slab geometry, an electric displacement field arises naturally in the solution. This electric displacement field is a consequence of the artificial change of topology imposed on the system in which the crystal blocks the passage of ions, resulting in a uniformly polarized electrolyte. Our analysis reveals that this is a consequence of the slab geometry itself, rather than the lattice summation techniques that are typically used to treat electrostatic interactions under periodic boundary conditions. We argue that tinfoil Ewald approaches, which impose zero electropotential difference in the electrolyte either side of the crystal, are appropriate when attempting to model the polar surfaces of crystals in solution, if the crystal is genuinely thin. In contrast, if one is interested in crystals that are macroscopic in extent, the electric displacement field can be chosen to enforce σ^(∞,eq), as demonstrated in Ref. 13.

Finally, our theoretical prediction [Eq. (7)] for the fluctuations in the adsorbed surface charge density is intriguing and suggests that dynamics at the interface may change significantly as the thickness of the crystal increases. Static energy calculations already reveal that the behavior of thin films in contact with vacuum can be complex, with structural relaxations that we have not considered in this work permitting uncompensated polarity to a certain extent.⁴⁷ The morphology of crystals is generally determined by the relative growth rates along different crystallographic directions.¹ Considering the potential interplay between structural relaxations, and the thickness dependence of both the equilibrium adsorbed surface charge density and its fluctuations, one can easily imagine that the growth mechanisms of a polar surface from either a melt or supersaturated solution will be a complex affair. In addition, we have implicitly assumed $A≫nR$ throughout, where $A$ is the surface area of the exposed polar facets in the crystal [Fig. 1(a)], whereas the polar catastrophe has nontrivial dependence on crystal morphology.^7,48 Future work will focus on disentangling these different contributions to the growth mechanisms of polar crystal surfaces.

The simulation methodology closely resembles that of Ref. 13, and indeed, for the R_c = R system, we reused trajectories for n = 3, 5, 7, 11, and 23. All simulations used the extended simple point charge water model (SPC/E),⁴⁹ whose geometry was constrained using the RATTLE algorithm⁵⁰ and the Joung–Cheatham NaCl force field.⁵¹ Dynamics were propagated using the velocity Verlet algorithm with a time step of 2 fs. The temperature was maintained at 298 K with a Nosé–Hoover chain^52,53 with a 0.2 ps damping constant. The particle–particle particle–mesh Ewald method was used to account for long-ranged interactions,⁵⁴ with parameters chosen such that the root mean square error in the forces were a factor 10⁵ smaller than the force between two unit charges separated by a distance of 1 Å.⁵⁵ A cutoff of 10 Å was used for non-electrostatic interactions.

The electrolyte comprised 600 water molecules and 20 NaCl ion pairs. The crystal consisted of alternating layers of Na⁺ and Cl⁻ ions, separated by R = 1.628 Å, and each layer comprised 16 ions. The lateral dimensions of the simulation cell were L_x = 15.952 Å and L_y = 13.815 Å along x and y, respectively. In the slab geometry with n = 3 and R_c = R, the length of the simulation cell along z was L = 94.841 Å, and L was increased with n accordingly. For instance, for n = 5, L was increased by 2R. Each simulation was 10 ns long post equilibration. The Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) simulation package was used throughout.⁵⁶ 

I am grateful to Rob Jack, Michiel Sprik, and Tom Sayer for many enlightening discussions on this topic and Rick Remsing for reading a draft of the manuscript. This work was supported by a Royal Society University Research Fellowship (Grant No. URF\R1\211144).

The author has no conflicts to disclose.

Stephen J. Cox: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Funding acquisition (lead); Investigation (lead); Methodology (lead); Writing – original draft (lead); Writing – review & editing (lead).

The data that support the findings of this study are openly available at https://doi.org/10.17863/CAM.83994.

In the main article, σ^(n,eq) was derived for a nonperiodic slab geometry [Eqs. (4) and (6)]. In this appendix, an expression for σ^(n,eq) is derived by considering the average molecular charge distribution under periodic boundary conditions. The derivation strongly hints that noncapacitive contributions arise from approximating σ^(n,eq) as being confined to a single plane close to the crystal’s surface (Fig. 1). The derivation largely resembles that of Ref. 16, but differs in that we do not consider a smeared average molecular charge distribution. Moreover, we take the present derivation as evidence that tinfoil Ewald sums provide a reliable estimate of σ^(n,eq), rather than ion adsorption resulting from a spurious treatment of electrostatic interactions.

When the average molecular charge density ρ_c only varies along z, the average electrostatic potential under tinfoil Ewald periodic boundary conditions is^57–59 

with

where L is the length of the simulation cell along z. For simplicity, let us first consider the case n = 1. In the simulations discussed in this article, the crystal ions do not move; approximating these as uniformly charged planes, we have

where $ρc(soln)$ is the average molecular charge distribution arising from the electrolyte solution (both solvent and ions). Let us assume that the crystal is centered at z = 0. In the region R/2 < z ≤ L/2, we then have

and an average electric field

At equilibrium, the electric field in the electrolyte vanishes, and, thus,

where z is understood to represent a point in the bulk electrolyte. The equilibrium electrolyte charge distribution satisfies Eq. (A6), no matter how complex its behavior. For a general $ρc(soln,eq)$⁠, however, Eq. (A6) is not straightforward to analyze.

To simplify matters, we introduce a model for $ρc(soln)$ whereby the adsorbed surface charge at each surface is confined to a single plane,

The relationship between $ℓ̄$ and the true average molecular charge distribution is left unspecified, and it should not necessarily be considered a physically meaningful length scale. For this model, the average electrostatic potential for $(R+2ℓ̄)/2<z≤L/2$ is

while the electric field is

Again enforcing the equilibrium condition, $ẑ⋅Ee=0$⁠, we find

which agrees with Eq. (6), provided that $ℓ̄=ℓ+anc$⁠. Comparing Eqs. (6), (A6), and (A10) suggests that noncapacitive contributions largely account for the approximation of treating the entire adsorbed surface charge density as if it were confined to a single plane.

Following the same procedure for general n, we find

To arrive at this result, we recall that n + 1 is an even integer, such that iⁿ⁺¹ = (−1)^(n+1)/2. At equilibrium,

Again, we find agreement with Eq. (6), provided that $ℓ̄=ℓ+anc$⁠. As we derived σ^(n,eq) by enforcing $ẑ⋅Ee=0$⁠, and as V(−L/2) = V(L/2) [Eqs. (A1) and (A2)], it immediately follows that V_a = V_b in the bulk electrolyte either side of the slab (Fig. 4; see also Fig. 4 of Ref. 13). Note that σ^(n,eq) is independent of L. In the limit L → ∞, any value of σ⁽ⁿ⁾ satisfies $ẑ⋅Ee=0$⁠, which reflects the fact that the electric field external to a set of infinite parallel plate capacitors vanishes. For σ⁽ⁿ⁾ ≠ σ^(n,eq), however, V_a ≠ V_b, and the system becomes increasingly unstable as n increases [Eq. (7)].