Charge compensation at the interface between the polar NaCl(111) surface and a NaCl aqueous solution

Periodic supercell models of electric double layers formed at the interface between a charged surface and an electrolyte are subject to serious finite size errors and require certain adjustments in the treatment of the long-range electrostatic interactions. In a previous publication (C. Zhang, M. Sprik, Phys. Rev. B 94, 245309 (2016)) we have shown how this can be achieved using finite field methods. The test system was the familiar simple point charge model of a NaCl aqueous solution confined between two oppositely charged walls. Here this method is extended to the interface between the (111) polar surface of a NaCl crystal and a high concentration NaCl aqueous solution. The crystal is kept completely rigid and the compensating charge screening the polarization can only be provided by the electrolyte. We verify that the excess electrolyte ionic charge at the interface conforms to the Tasker 1/2 rule for compensating charge in the theory of polar rocksalt (111) surfaces. The interface can be viewed as an electric double layer with a net charge. We define a generalized Helmholtz capacitance $C_\text{H}$ which can be computed by varying the applied electric field. We find $C_\text{H} = 8.23 \, \mu \mathrm{Fcm}^{-2}$, which should be compared to the $4.23 \, \mu \mathrm{Fcm}^{-2}$ for the (100) non-polar surface of the same NaCl crystal. This is rationalized by the observation that compensating ions shed their first solvation shell adsorbing as contact ions pairs on the polar surface.


I. INTRODUCTION
Crystals exposing a face bearing a net charge are intrinsically unstable if, in addition, the unit cell also has a net dipole moment perpendicular to the surface.Such a termination is referred to as a type III polar surface in the classification by Tasker 1 .Tasker explained the instability of type III surfaces by showing that the energy diverges with increasing thickness of the crystal (polar catastrophe).Yet, surfaces with type III orientations do occur in nature.The electrostatic instability is avoided by the accumulation of compensating charge which cancels the dipole moment. 2 Various compensating mechanisms have been observed or suggested.The dominant mechanisms, as reviewed by Noguera, 2 are a change of surface composition (nonstoichiometric reconstruction), adsorption of charged species, and electronic reconstruction (charging of gap or defect states).The review of Ref. 2 is restricted to oxide materials.It was updated in 2008 in collaboration with Goniakowski 3 and again in 2013, now also including nano objects, such as thin films. 4Thin films are of special interest because some structures can exhibit an unreconstructed polar surface.Because of their small size they can sustain the polarization field driving the instability in larger systems.
Polar surfaces are clearly a challenge for computational methods.The problem of how to calculate the total energy of model systems was first addressed as a theoretical exercise in the summation of electrostatic interactions arising from an array of point charges 1,5,6 which was followed later by studies involving electronic a) Electronic mail: tes36@cam.ac.uk structure calculations.The favourite model systems are: MgO (111)(rocksalt), free standing [7][8][9][10][11][12][13] or deposited on a metal substrate; 14 , NiO(111)(rocksalt), 9 and ZnO, either the (111) zinc blende 11 or (0001) (wurtzite) 15 surface.More complex surfaces that have been modelled include tetragonal ZrO 2 (110) 16,17 and α-Al 2 O 3 (0001) 18 .][21] Atomistic model systems necessarily have a slab geometry of limited width.Electric fields in the solid are permitted. 22Similar to the thin films of experiment, the finite slab width can therefore mask a polar catastrophe. 4n this contribution we will outline a finite field approach adjusting the electrostatics.The method is inspired by our work on the atomistic modelling of an electric double layer (EDL) formed by a solid in contact with an aqueous solution (electrolyte). 23Finite electric fields penetrating the solid are again a serious concern if the solid is an insulating mineral.The surface charge is of chemical origin: the result of exchange (adsorption or desorption) with an ionic solution (electrolyte).The surface charge is compensated by a zone of excess ionic charge on the electrolyte side of the interface.In a macroscopic, semiinfinite solid the net charge in the EDL is zero.However this basic rule of EDL theory can be violated if the system is represented by a slab with surfaces of opposite charge by creating an internal electric field.The residual electric field manifests as a finite, net charge in the EDL.The surface charge is not properly screened by the mobile ions of the solution, instead the slab acts as a nano-capacitor.
The nanocapacitor effect, while of interest by itself, must be regarded as a finite system size error if the aim is to simulate an interface between a semi-infinite insulator and an electrolyte.The magnitude of the error was investigated in Ref. 23 for the familiar simple point charge (SPC) model of an aqueous NaCl solution confined between two oppositely charged walls.The walls separating the solution from vacuum were only an atomic diameter thick.Since full 3D periodic boundary conditions (PBC) were applied, moving the supercell over half its length generates a different perspective of the same system.Now the vacuum is in the middle and the system can be viewed as parallel plate capacitor.Even for gaps as large as 20 to 100 Å, the deficit charge in the EDL turns out to be a significant fraction of the surface charge (up to 20% at 20 Å).
The results of Ref. 23 were interpreted in terms of a simple continuum (Stern) model of the two EDL's.Within this model the missing EDL charge was found to be linearly correlated with the electric field in the insulator (the vacuum space).This observation suggested that charge neutrality of the EDL can be restored by applying an external bias field of the correct magnitude to cancel the field in the insulator.We verified that this is indeed the case, applying a classical molecular dynamics (MD) version of the finite field methods developed by Stengel et al. [24][25][26] A further prediction of the continuum model was that the compensating field, termed the field of zero net charge (ZNC), is proportional to the inverse capacitance of the EDL.This relation was also confirmed by the atomistic simulation and used to compute the capacitance. 23n the present contribution the finite field methodology of Ref. 23 is extended to polar interfaces.The idea is again that application of an appropriate external electric field should cancel the internal field associated with the polarization.In this first application the configuration of the ions in the solid is strictly fixed.Only the ions and water molecules in the electrolyte move and should therefore screen the polarization.2][3] The Tasker rule makes a very precise statement about the compensating charge which should be relatively easy to verify in a simple classical point charge model as will be used here.

A. Point charge model for the NaCl-electrolyte interface
The main objective of this study is validation of the finite field method for electrochemical interfaces of polar surfaces.We opted therefore for the simplest of model systems, a slice of NaCl crystal of (111) orientation.The slab is terminated on one side by a Na + plane and the other side by a Cl − plane and is kept rigidly in the bulk rock salt geometry.The solid is in contact with a high concentration aqueous NaCl solution.Fig. 1a shows an instantaneous configuration sampled from the MD simulation.The classical force model is the same as used in The theory (Tasker rule) says that the surface charge density required to cancel the dipole generated by cleaving a rocksalt structure along a (111) plane is −σ 0 /2, where σ 0 is the surface charge density of the terminating plane.A crystal with this solid/electrolyte interface should be stable and have no net internal electric field.The procedure is therefore the same as in Ref. 23.The applied bias field is varied until the internal Maxwell field in the slab is on average zero.We computed the charge imbalance (space charge) in the electrolyte in contact with the crystal face and checked whether it has the theoretical value of −σ 0 /2.

B. Stern model for the polar surface-electrolyte interface
Anticipating our results we found that the Tasker half surface charge density rule for the unreconstructed rocksalt (111) termination is indeed satisfied for our model.This raises questions about the electrostatics of the polar surface/electrolyte interface.What are the differences compared to the regular charge neutral EDL induced by the surface charge of a non-polar dielectric solid?What is the capacitance or how can we even define a capacitance?We will address this problem by generalizing the continuum Stern model of Ref. 23.The model is pictured in Fig. 2. As before the electrolyte is partitioned in a proper ionic conductor ( = ∞) and two boundary layers of each of thickness l H .The dielectric constant in a boundary layer, to which we will continue to refer as a "Helmholtz layer", is H .The Maxwell field in the conducting region is strictly zero while it can be finite in a Helmholtz layer.This field is denoted by E H .The polarization of the electrolyte is represented by the surface charge density ±σ of the planes separating the conductor and the Helmholtz layer.With the surface charge of the solid fixed, the electrolyte surface charge σ is the central variable in the model.
To model the polar slab we follow Noguera, 2 representing the solid by a succession of planes with alternating surface charge density σ 0 .The planes are a distance R apart.There are n + 1 of these planes dividing the solid up into n slices, where n is an odd number.The dielectric constant is homogeneous throughout the solid and will again be indicated by d .The electric field is however not the same everywhere.Counting from the left in Fig. 2 the field E 1 in the first layer (i = 1) is different from the field E 2 in the next layer (i = 2).Because of the strict alteration of the surface charges on the planes the field in all odd numbered regions is E 1 and in the even numbered layers E 2 (note the convention of the field directions in Fig. 2).To represent a MD supercell the Stern model is periodically repeated in the normal direction.The length of the supercell is L.
In the finite field method of Ref. 23 the periodic MD system is subject to an electric field Ē using an extended Hamiltonian (see section II C).Note that Ē is not an external field E 0 but the average of the Maxwell field.The product V = − ĒL can therefore be directly interpreted as the potential difference across the MD cell.The corresponding cell potential V is the sum of all potential differences over the uniform layers making up the system and we can write with the field in the electrolyte set to zero (see Fig. 2).n Ēd R is (minus) the potential across the polar solid with average field Ēd given by The Maxwell equations for boundaries give the following dependencies: Exchanging the fields in the right-hand side of Eq. 1 for charge densities using Eqs.3, 4 and 5 we obtain an expression for σ in terms of σ 0 and Ē, where C d = d /(4πR) and C H = H /(4πl H ). The parameters C d and C H have the familiar form of the capacitance of a parallel plate capacitor and will be interpreted as such.
The surface charge σ 0 is fixed.The Maxwell field Ē is the only control parameter in the model.Ē can be varied until the average field Ēd within the crystal is zero.
In this state, referred to as the point of "Compensating Net Charge" (CNC), the two interfaces of the slab are decoupled, and so the infamous finite size error is removed.Setting Ēd = 0 in Eq. 1 and inserting into Eq.6 using Eq. 4 gives the compensating charge provided by the electrolyte: Indeed in the limit n → ∞, the prefactor of equation 7 tends towards 1 2 in agreement with the Tasker rule for a rocksalt(111) polar surface. 1,2The compensating charge is determined by the surface charge only and is independent of all other structural parameters.
The Ē = ĒCNC field plays the same role as the field of zero net charge ĒZNC of Ref. 23.At Ē = ĒZNC the field in the dielectric slab is zero restoring net charge neutrality to the EDL.Moreover the Stern model led to a simple relation between the capacitance of the EDL and ĒZNC (Eq.37 of Ref. 23) Pursuing the parallel with the "regular" EDL of Ref. 23 further, we can similarly express the capacitance of the Helmholtz layer in terms of ĒCNC .Combining Eq. 6 and 7, equivalent to imposing a potential − ĒCNC L, we find Recall that the factor 2 multiplying σ 0 in Eqs. 9 and 8 is there because −L Ē is the potential over a pair of EDL's in series each with a capacitance C H .The status of C H as the Helmholtz capacitance of the polar surface/electrolyte interface may look ambiguous because the response charge of the electrolyte is only half the surface charge.It is therefore of interest to compare the capacitance for the polar NaCl (111) surface as defined by Eq. 9 to the capacitance obtained by charging a non-polar surface.The obvious candidate is the (100) surface.The surface was charged by slightly enhancing the positively (negatively) charged ions compared to the counter charge.The expression for C H is given by Eq. 9 leaving out the factor (n + 1)/(2n) consistent with Eq. 8.

C. Hamiltonian and finite electric field molecular dynamics
All simulations were performed under ambient conditions using a modified version of the GROMACS package. 27The water model is SPC/Extended. 28The SPC model of aqueous Na + and Cl − was taken from Ref. 29 and has been validated for high salt concentrations. 30,31(see also the review by Nezbeda et al. 32 ).Identical force field parameters were used for the interactions with the ions in the (rigid) NaCl crystal.The supercell cross-sectional area were 2.200 and 2.545 nm 2 for the (111) and (100) orientation respectively.The corresponding distance between the charged planes of the (111) slab, referred to as R in Fig. 2, is R = 1.63 Å.The MD cell lengths were adjusted to keep variations in solvent properties to a minimum.All supercells contained 20 aqueous NaCl ion pairs in ∼ 600 waters, evenly dispersed in the initial configuration.This particular concentration was chosen such that the electrolyte remains a good ionic conductor after the formation of EDLs has reduced the concentration (the total number of ions is fixed).The MD parameters were as follows: the NVT ensemble was employed, and a temperature of 298 K mainted by Nosé-Hoover thermostat with a coupling constant of 0.4 ps; 33 the timestep was 2 fs and the simulations were run for a total of 1 ns; the electrostatics were computed using 4 th order Particle Mesh Ewald (PME) summation with a Fourier spacing of 0.6 Å and a realspace cutoff of 6.5 Å. 34 The first 200 ps were discarded as equilibration time.The system was considered equilibrated if the electrostatic potential in the bulk electrolyte was flat.Ref 23 investigated the equilibration time for C H and the choice of 200 ps achieves this level of convergence in the present work.
The modification of the GROMACS package concerns the implementation of the constant field Hamiltonian method according Vanderbilt and co-workers. 24Forces and energies are derived from the extended Hamiltonian where H PBC (v) is the Hamiltonian as defined by the SPC model.v = (r N , p N ) stands for the collective momenta and position coordinates of the N particles in the system.The subscript PBC indicates that the electrostatic energies and forces are computed using standard Ewald summation ("tin foil" boundary conditions).Ω is the volume of the MD cell.P is the polarization perpendicular to the crystal slab with Ē the magnitude of the electric field.P is computed from the total dipole moment of the supercell including the contribution from the ions in the solution and the solid.Typical values of Ē applied to our model system are in the order of 1 V/nm.Finite electric field Hamiltonians of the form of Eq. 10 are familiar in classical MD. 35 As already pointed out by Yeh and Berkowitz, the electric field in the dipole coupling term −Ω ĒP , when combined with Ewald summation, must be interpreted as the average of the Maxwell field rather than the applied external field. 35n Ref. 24 this feature, specific to standard Ewald summation, is given a firm thermodynamic foundation (see also Ref. 25).A further important point is that including ions in the polarization makes P a multi-valued function depending on the supercell boundary.This is a central theme in the modern theory of polarization which applies to electronic polarization as well as classical point charge systems. 36This issue becomes critical for a supercell of the geometry of Fig. 1.The mobile ions can cross the cell boundaries.When this happens the ions must be followed out of the cell in the calculation of the polarization.In classical simulation, this definition of polarization was introduced in the context of the MD study of electrolytic solutions and is generally referred to as itinerant polarization 37 .Note that the forces derived from the Hamiltonian Eq. 10 are not affected by ambiguities in the definition of polarization (see Ref. 23 for a more detailed discussion).

A. Finding the field of compensating net charge (CNC)
The first step is locating the field of compensating net charge.In the state of CNC, Ēd = 0 by definition.The change ∆φ d in the electrostatic potential over the length of the crystal must therefore be zero.∆φ d was computed from the electrostatic potential profile φ(z).As is evident from Fig. 3(a) the potential profile, while flat in the electrolyte, shows a saw tooth pattern in the crystal reflecting the alternating charge of the (111) crystal planes.In contrast, the potential profile in the non-polar slab (Fig. 3(b)) is smooth because (100) planes are net neutral.In a CNC state ( Ēd = 0), from which the configuration of Fig. 3(b) was sampled, the potential is therefore constant (∆φ d = 0).
It is perhaps instructive to analyze the polar profile of Fig. 3(a) in somewhat more detail because its appearance may at first seem inconsistent with the condition of vanishing average internal electric field required at CNC. Referring back to Fig. 2 we note that the period of the modulation of the potential in Fig. 3(a) is 2R.However, the width of the slab is nR where n is an odd integer.The CNC field therefore aligns a maximum of the saw tooth at one face with a minimum at the opposite face as indicated by the red dashed line in the figure.This can only be achieved by canting the sequence of maxima (minima) which is the feature that stands out in Fig. 3(a).∆φ d was tracked against the applied field Ē to locate the point of CNC for increasing width of the crystal.The result is shown in Fig. 4(c).The width of the slab is represented by the number n of capacitors in series (see Fig. 2).The figure suggests that a minimum of n = 15 layers (16 planes) is needed to reach convergence.This corresponds to a slab width of 25 Å.In the search for the CNC it was observed that the response of charge and potential to changes in the field was linear.The exception is the terminal n=1 system.This is likely due to dielectric saturation at the high CNC field for n=1 as suggested by Fig. 4(c).Therefore, the values for n=1 are obtained by extrapolation from the low-field, linear regime.
Slab width is of course always a critical parameter in periodic models of interfaces.Solid NaCl/aqueous NaCl interfaces are popular model systems and the question of size dependence has been investigated in great detail in calculations of solubility from the direct equilibrium between solid and solution 32 .Using an interaction model identical to the one used here, Espinoza and coworkers found in a careful study 38 that the solubility is essentially converged for slabs of a width larger than 40 Å.This number exceeds but is still comparable to the 25 Å inferred from Fig. 2. Considering however, that the orientation of the slab used in Ref. 38 is the nonpolar (100) surface, the similarity between the finite size effects is perhaps somewhat surprising.In this context we should also point out that the size effect for our polar interface seems to be captured to great accuracy by an analytical expression (Eq.7) specified by only a single parameter, the surface charge density σ 0 of the polar (111) surface for which there is no direct counter part in the dissolution of a (100) interface.On the other hand the consistency in finite size effects can also be interpreted as an indication that surface charge and polarity play a role in the dissolution of non-polar interfaces.

B. Compensating electrolyte charge and capacitance
In the simple Stern model of Fig. 2, the charge induced in the electrolyte is represented as a surface charge density σ at the sharp interface between the electrolyte and the Helmholtz layer.In the atomistic model system of Fig. 1, the interface is more diffuse.We have estimated the compensating charge as the integral of the excess charge density ("space charge") of the electrolyte in contact with the surface.This procedure is outlined in more detail in Ref. 23.The excess charge per unit area is identified with the σ of the Stern model.The results are plotted in Fig. 4(a) where we have represented the estimates of σ as total surface charges Aσ with A the MD cell cross section.
The surface charge of a NaCl (111) plane in our model system is Aσ 0 = 16e.The theoretical compensating charge (Tasker rule) is therefore 8e.Fig. 4(a) shows that the compensating charge Aσ approaches 8e when the number of layers of polar solid gets large enough.Therefore, our simulation confirms that the charge imbalance is in accord with the theoretical value.Indeed, figure 4(b) shows excellent agreement with Eq. 7 derived from the continuum model Fig. 2.This also gives us confidence that Eq. 9 can be used to extract a value for the Helmholtz capacitance from simulations.Indeed, Fig. 4(d) shows the desired linear (n+1)/2n dependence of Eq. 9.The slopes gives an estimation of C H of 8.23 µFcm −2 for the polar (111) surface.
8.2 µFcm −2 is a relatively modest capacitance, not much larger than the 4.4 µFcm −2 we found in Ref. 23  for the EDL formed by a 1.4 M aqueous NaCl electrolyte confined between two walls of uniform opposite surface charge densities.While our procedure of a fractional increase of the charge of all ions of one species (Na + or Cl − ) is nonphysical, we argue that the electrostatics of an EDL formed by this system is similar to that of the "regular" EDL of Ref. 23. Accordingly, the Helmholtz capacitance can be obtained from Eq. 8 or equivalently from Eq. 9 by leaving out the (n+1)/2n factor.The EDL of the non-polar surface is charge balanced for vanishing internal field.The CNC and ZNC are interchangeable in this case.The result of the applied voltage at CNC as a function of the surface charge for the non-polar (100) surface is shown in Fig. 5.The slope yields the Helmholtz capacitance C H with a value 4.23 µFcm −2 , effectively identical to the capacitance obtained in Ref. 23.

C. Electric double layer structure
The near 2 to 1 ratio for the capacitance of the polar (111) and non-polar (100) surfaces may not be a coincidence.A possible explanation for this behaviour can be seen in the double layer structure of Fig. 6.For the (111) case, the ions of the double layer are adsorbed to the surface at a distance of ∼1.5 Å, shedding their inner solvation shell.The formation of contact ion pairs is particularly noteworthy for the Na + ions, which are traditionally thought to have an almost unbreakable solvation shell.In the (100) case, at the lowest charge of 2e, there is no driving force for such a dehydration, and the first peak in the density occurs at the larger distance of 2.8 Å, followed by a small secondary peak further out.As the surface charge is increased, the profile shifts towards the surface.Even for a (100) surface charge as high as the compensating charge for the (111) polar surface (8e), the structural difference between EDLs is apparent.The capacitance of a compact (Helmholtz) EDL is inversely proportional to its width, therefore, the 2 to 1 ratio in C H between the surfaces is roughly consistent with the relative positions of the counter ions in Fig. 6.

IV. CONCLUSION AND OUTLOOK
Our calculations show that the excess charge in a high concentration NaCl aqueous solution adjacent to a rigid NaCl (111) surface complies to the Tasker rule for a polar surface of this geometry (half the surface charge).This was to be expected.The Tasker rule is based on general considerations involving geometry and electrostatics.It would have been a surprise if this rule would not hold for polar surface electrolyte interfaces.The results of the present study are therefore intended as further validation of our finite field method for the simulation of electric double layers under full periodic boundary conditions.This method was introduced in Ref. 23 and applied to a conventional electric double layer for which the electrolyte counter charge equals the opposite of the surface charge (zero net charge).Reproducing the half charge predicted by the Tasker rule for a model polar surface seemed to us a separate challenge which was taken up in the calculation reported here.
The key function of the applied external field was to compensate for the internal electric field in polar model slabs.Slabs of a width accessible to molecular simulation can sustain an internal field leading to violation of the Tasker rule for polar surfaces of semi-infinite crystals.The fact that the electrolyte is an ionic conductor is essential.This enforces a zero macroscopic field in the bulk region of the electrolyte whatever the magnitude of the polarization of the slab or the applied field.The external field can therefore be adjusted to cancel the internal field in the solid without inducing a field in the electrolyte.
The model system in this feasibility study was deliberately kept as simple as possible.In particular the structure of the solid slab was constrained to be rigid.This meant that a number of interesting issues could not be addressed.One most important question is the competition with non-stoichiometric reconstructions observed for vacuum surfaces.One of the candidate structures for the rocksalt (111) surface is the so-called octopolar reconstruction suggested by Wolf 6 (see also Noguera 2,3 ).Dissolution of NaCl is a facile process interfering with reconstruction and the stability of polar surfaces is usually studied for more robust rocksalt crystals such as MgO 12 and NiO 39 .
However, under saturation the NaCl surface is in equilibrium with its aqueous solution and it has been proven to be feasible to stabilize a NaCl(111)/NaCl(aq) interface under appropriate thermodynamic conditions 40 .The force field used in the present study (the Joung-Cheetham model of Ref. 29) is also suitable to carry out a simulation of a solid/solution equilibrium.After years of hard work there seems to be now a consensus in the computational literature regarding the solubility of this model (see the review of Ref. 32).This is below the experimental value, but above the 2 M concentration used in the present study for the solution phase.This ensures that there is no net tendency of the crystal to grow.These observations apply to non-polar surfaces (see however the results of Ref. 38 for a spherical piece of crystal).It would therefore be of interest to repeat these calculations for a polar interface.
Finally we return to the capacitance we computed in section III B for the polar surface electrolyte interface.The possibility of experimental realization of an unreconstructed polar surface electrolyte interface raises the question of the interpretation of this capacitance and how it might be observed.While the structure of the "double layer" should be a least in principle verifiable by experiment, the status of the corresponding capacitance is less clear.Capacitance is well defined for a nanoslab and should also be accessible to experiment.It would be the capacitance of the nanocapacitor formed by the solid slab with the electrolyte on either side acting as the two electrodes.But how to determine the capacitance of the interface of a non-conducting (insulator) semi-infinite crystal with a polar termination?The capacitance probably enters experimentally-measurable quantities only indirectly, such as in adsorption (complexation) energies of charged species.Maybe the double layer will also show electrokinetic signatures, such as a finite zeta-potential.This rather technical report is not the place to address these questions.This will have to be resolved in future investigations.
In conclusion we reiterate the observation by Noguera that the value of the compensating charge is a result of long range electrostatics and should remain the same even if the electronic structure is taken into account 2 .We therefore anticipate that this method can enable us to study the interaction of an electrolyte with more complex and realistic polar surfaces, possibly even treated by electronic structure calculation.

FIG. 1 :
FIG. 1: MD snapshots of the polar (111) (top) and non-polar (100) (bottom) surfaces of a rigid crystalline NaCl slab interfaced with a high concentration NaCl aqueous solution.Na + ions are depicted in blue, Cl − ions in yellow.The non-polar (100) surface is given an artificial net charge of 8e matching the surface charge 16e of the polar(111) surface.A finite electric field is applied cancelling the internal field of the slab.With the polarization field removed the excess charge in the electrolyte near the surface should compensate the slab polarization.

FIG. 2 :
FIG. 2: Schematic drawing of the continuum model of a polar surface-electrolyte system under periodic boundary conditions (PBC).The (absolute) surface charge density of a polar surface is σ 0 and the compensating charge induced in the electrolyte solution is σ.The solid slab is separated from the electrolyte on both sides by Helmholtz layers.The dielectric constants of the Helmholtz layers and polar solid are H and d respectively.The box size is L, the width of Helmholtz layer is l H and the thickness of a layer in polar solid is R.The arrows indicate the convention for the sign of the uniform electric fields in the Helmholtz layers and crystal segments.

FIG. 3 :
FIG. 3: (a) Potential profile (potential vs distance) for the NaCl (111)-electrolyte system for n=23 at ĒCNC L = 1.2 V.The first and last crystal peak can be seen to occur at the same potential conform the requirement for CNC.(b) Potential profile for the NaCl (100)-electrolyte system with surface charge 8e at ĒCNC L = 1.8 V.

FIG. 4 :
FIG. 4: (a) The effective surface charge of the electrolyte at the point of CNC for increasing width n of the slab (see Fig. 2).It is represented as the difference between the total fixed surface charge Aσ 0 of the crystal and the response charge Aσ from the electrolyte solution, where A is the area of the polar plane in the supercell.Aσ 0 = 16e in our model.(b)Plot of A(σ 0 − σ) vs (n + 1)/2n.The slope of −Aσ 0 demonstrates that Eq. 7 is obeyed.(c) Plot of ĒCNC L vs n.(d) Plot of ĒCNC L vs (n + 1)/2n.The slope of 2σ 0 /C H allows us to extract the value of the Helmholtz capacitance C H .

FIG. 5 :
FIG. 5: CNC Potential ĒCNC L (is potential of zero net EDL charge) as a function of the artificially enhanced surface charge for the NaCl (100)-electrolyte system.The slope of the linear fitting (red solid line) is linked to the Helmholtz capacitance C H according to Eq. 8.

FIG. 6 :
FIG.6: Density plot of compensating ions (excess charge) from the electrolytic solution at points of CNC.d is defined as the distance of ions from the relevant surface of NaCl crystals.The solid black line is the density profile for the (111) polar surface with surface charge 16e.The solid red line is for the (100) non-polar surface with artificial charge Aσ 0 = 8e.The long-dashed orange, short-dashed green, and dotted blue lines are the (100) surface with Aσ 0 values of 6e, 4e, and 2e respectively.