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

TL;DR

This study extends finite field electrostatic methods to analyze charge compensation at the polar NaCl(111) surface in contact with an aqueous solution, revealing how ionic screening adheres to theoretical rules and affects capacitance.

Contribution

It applies and verifies finite field methods to a polar crystal surface in electrolyte, demonstrating charge compensation mechanisms and calculating the Helmholtz capacitance.

Findings

Excess ionic charge conforms to Tasker 1/2 rule.

Helmholtz capacitance for NaCl(111) surface is 8.23 μF/cm².

Contact ion pairs form on the polar surface, affecting capacitance.

Abstract

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.

Figures

Equations

$$\bar{E}L=-2E_{\mathrm{H}}{l}_{\mathrm{H}}+n{\bar{E}}_{d}R,$$

$${\bar{E}}_{\text{d}}=-\frac{n-1}{2n}{E}_2+\frac{n+1}{2n}{E}_1.$$

$${ϵ}_{\mathrm{H}}{E}_{\mathrm{H}}+{ϵ}_{\text{d}}{E}_1$$

$${ϵ}_{\mathrm{H}}{E}_{\mathrm{H}}$$

$${ϵ}_{\text{d}}{E}_2+{ϵ}_{\text{d}}{E}_1$$

$$\sigma =\left(\frac{n+1}{2}\frac{{\sigma }_0}{C_{\text{d}}}-\bar{E}L\right){\left(\frac{2}{C_{\mathrm{H}}}+\frac{n}{C_{\text{d}}}\right)}^{-1},$$

$$\sigma =\frac{n+1}{2n}{\sigma }_0.$$

$$C_{\mathrm{H}}=\frac{2{\sigma }_0}{{\bar{E}}_{\mathrm{ZNC}}L}$$

$$C_{\mathrm{H}}=\frac{2{\sigma }_0}{{\bar{E}}_{\mathrm{CNC}}L}\frac{n+1}{2n}.$$

$$H_E\left(v,\bar{E}\right)=H_{\mathrm{PBC}}(v)-\Omega \bar{E}P(v),$$

Full text

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 TaskerTasker (1979). 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.Noguera (2000)

Various compensating mechanisms have been observed or suggested. The dominant mechanisms, as reviewed by Noguera,Noguera (2000) are a change of surface composition (non-stoichiometric 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 GoniakowskiGoniakowski, Finocchi, and Noguera (2008) and again in 2013, now also including nano objects, such as thin films.Noguera and Goniakowski (2013) Thin 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 chargesTasker (1979); Nosker, Mark, and Levine (1970); Wolf (1992) which was followed later by studies involving electronic structure calculations. The favourite model systems are: MgO (111)(rocksalt), free standingRefson et al. (1995); Goniakowski and Noguera (2002); Wander, Bush, and Harrison (2003); Finocchi et al. (2004); Goniakowski, Noguera, and Giordano (2007); Spagnoli, Allen, and Parker (2011); Gaddy et al. (2014) or deposited on a metal substrate;Goniakowski, Giordano, and Noguera (2010), NiO(111)(rocksalt),Wander, Bush, and Harrison (2003) and ZnO, either the (111) zinc blendeGoniakowski, Noguera, and Giordano (2007) or (0001) (wurtzite) Meyer and Marx (2003) surface. More complex surfaces that have been modelled include tetragonal ZrO2(110)Anez, Sierraalta, and Sautet (2007); Anez et al. (2009) and $\alpha$-$\mathrm{Al}_{2}\mathrm{O}_{3}(0001)$Ruberto, Yourdshahyan, and Lundqvist (2003). More recently ternary oxides have attracted attention, in particular polar terminations of LaTiO3 and interfaces with SrTiO3.Bristowe, Littlewood, and Artacho (2011); Chen et al. (2015); Krishnaswamy et al. (2015)

Atomistic model systems necessarily have a slab geometry of limited width. Electric fields in the solid are permitted.Bengtsson (1999) Similar to the thin films of experiment, the finite slab width can therefore mask a polar catastrophe.Noguera and Goniakowski (2013) In 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).Zhang and Sprik (2016a) Finite 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, semi-infinite 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.Stengel, Spaldin, and Vanderbilt (2009); Zhang and Sprik (2016b); Zhang, Hutter, and Sprik (2016) 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.Zhang and Sprik (2016a)

In 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. The question is thus whether the excess charge supplied by the electrolyte can play the role of the compensating charge density in the theory of polar surfaces.Tasker (1979); Noguera (2000); Goniakowski, Finocchi, and Noguera (2008) 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.

II Model and method

II.1 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 Ref. 23 (for further detail see section II.3).

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 $-\sigma_{0}/2$, where $\sigma_{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 $-\sigma_{0}/2$.

II.2 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. Zhang and Sprik, 2016a. The model is pictured in Fig. 2. As before the electrolyte is partitioned in a proper ionic conductor ($\epsilon=\infty$) and two boundary layers of each of thickness $l_{\mathrm{H}}$. The dielectric constant in a boundary layer, to which we will continue to refer as a “Helmholtz layer”, is $\epsilon_{\mathrm{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_{\mathrm{H}}$. The polarization of the electrolyte is represented by the surface charge density $\pm\sigma$ of the planes separating the conductor and the Helmholtz layer. With the surface charge of the solid fixed, the electrolyte surface charge $\sigma$ is the central variable in the model.

To model the polar slab we follow Noguera,Noguera (2000) representing the solid by a succession of planes with alternating surface charge density $\sigma_{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 $\epsilon_{\text{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. Zhang and Sprik, 2016a the periodic MD system is subject to an electric field $\bar{E}$ using an extended Hamiltonian (see section II.3). Note that $\bar{E}$ is not an external field $E_{0}$ but the average of the Maxwell field. The product $V=-\bar{E}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

[TABLE]

with the field in the electrolyte set to zero (see Fig. 2). $n\bar{E}_{d}R$ is (minus) the potential across the polar solid with average field $\bar{E}_{d}$ given by

[TABLE]

The Maxwell equations for boundaries give the following dependencies:

[TABLE]

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 $\sigma$ in terms of $\sigma_{0}$ and $\bar{E}$,

[TABLE]

where $C_{\text{d}}=\epsilon_{\text{d}}/(4\pi R)$ and $C_{\mathrm{H}}=\epsilon_{\mathrm{H}}/(4\pi l_{\mathrm{H}})$. The parameters $C_{\text{d}}$ and $C_{\mathrm{H}}$ have the familiar form of the capacitance of a parallel plate capacitor and will be interpreted as such.

The surface charge $\sigma_{0}$ is fixed. The Maxwell field $\bar{E}$ is the only control parameter in the model. $\bar{E}$ can be varied until the average field $\bar{E}_{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 $\bar{E}_{d}=0$ in Eq. 1 and inserting into Eq. 6 using Eq. 4 gives the compensating charge provided by the electrolyte:

[TABLE]

Indeed in the limit $n\rightarrow\infty$, the prefactor of equation 7 tends towards $\frac{1}{2}$ in agreement with the Tasker rule for a rocksalt(111) polar surface.Tasker (1979); Noguera (2000) The compensating charge is determined by the surface charge only and is independent of all other structural parameters.

The $\bar{E}=\bar{E}_{\mathrm{CNC}}$ field plays the same role as the field of zero net charge $\bar{E}_{\mathrm{ZNC}}$ of Ref. Zhang and Sprik, 2016a. At $\bar{E}=\bar{E}_{\mathrm{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 $\bar{E}_{\mathrm{ZNC}}$ (Eq. 37 of Ref. Zhang and Sprik, 2016a)

[TABLE]

Pursuing the parallel with the “regular” EDL of Ref. Zhang and Sprik, 2016a further, we can similarly express the capacitance of the Helmholtz layer in terms of $\bar{E}_{\mathrm{CNC}}$. Combining Eq. 6 and 7, equivalent to imposing a potential $-\bar{E}_{\mathrm{CNC}}L$, we find

[TABLE]

Recall that the factor 2 multiplying $\sigma_{0}$ in Eqs. 9 and 8 is there because $-L\bar{E}$ is the potential over a pair of EDL’s in series each with a capacitance $C_{\mathrm{H}}$.

The status of $C_{\mathrm{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_{\mathrm{H}}$ is given by Eq. 9 leaving out the factor $(n+1)/(2n)$ consistent with Eq. 8.

II.3 Hamiltonian and finite electric field molecular dynamics

All simulations were performed under ambient conditions using a modified version of the GROMACS package.Hess et al. (2008) The water model is SPC/Extended.Berendsen, Grigera, and Straatsma (1987) The SPC model of aqueous Na+ and Cl- was taken from Ref. Joung and T E Cheatham, 2008 and has been validated for high salt concentrations.Zhang et al. (2010); Zhang and Carloni (2012) (see also the review by Nezbeda et al. Nezbeda, Moučka, and Smith (2016)). 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 nm2 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 $\sim$ 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;Martyna, Klein, and Tuckerman (1992) the timestep was 2 fs and the simulations were run for a total of 1 ns; the electrostatics were computed using 4th order Particle Mesh Ewald (PME) summation with a Fourier spacing of 0.6 Å and a real-space cutoff of 6.5 Å.Darden, York, and Pedersen (1993) 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 Zhang and Sprik, 2016a investigated the equilibration time for $C_{\mathrm{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.Stengel, Spaldin, and Vanderbilt (2009) Forces and energies are derived from the extended Hamiltonian

[TABLE]

where $H_{\textrm{PBC}}(v)$ is the Hamiltonian as defined by the SPC model. $v=(\mathbf{r}^{N},\mathbf{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). $\mathit{\Omega}$ is the volume of the MD cell. $P$ is the polarization perpendicular to the crystal slab with $\bar{E}$ 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 $\bar{E}$ 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.Yeh and Berkowitz (1999) As already pointed out by Yeh and Berkowitz, the electric field in the dipole coupling term $-\mathit{\Omega}\bar{E}P$, when combined with Ewald summation, must be interpreted as the average of the Maxwell field rather than the applied external field.Yeh and Berkowitz (1999) In Ref. Stengel, Spaldin, and Vanderbilt, 2009 this feature, specific to standard Ewald summation, is given a firm thermodynamic foundation (see also Ref. Zhang and Sprik, 2016b). 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.Resta and Vanderbilt (2007) This 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 polarizationCaillol (1994). Note that the forces derived from the Hamiltonian Eq. 10 are not affected by ambiguities in the definition of polarization (see Ref. Zhang and Sprik, 2016a for a more detailed discussion).

III Results

III.1 Finding the field of compensating net charge (CNC)

The first step is locating the field of compensating net charge. In the state of CNC, $\bar{E}_{\text{d}}=0$ by definition. The change $\Delta\phi_{\text{d}}$ in the electrostatic potential over the length of the crystal must therefore be zero. $\Delta\phi_{\text{d}}$ was computed from the electrostatic potential profile $\phi(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 ($\bar{E}_{\text{d}}=0$), from which the configuration of Fig. 3(b) was sampled, the potential is therefore constant ($\Delta\phi_{\text{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).

$\Delta\phi_{\text{d}}$ was tracked against the applied field $\bar{E}$ 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 solutionNezbeda, Moučka, and Smith (2016). Using an interaction model identical to the one used here, Espinoza and coworkers found in a careful studyEspinosa et al. (2016) 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 $\sigma_{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.

III.2 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 $\sigma$ 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. Zhang and Sprik, 2016a. The excess charge per unit area is identified with the $\sigma$ of the Stern model. The results are plotted in Fig. 4(a) where we have represented the estimates of $\sigma$ as total surface charges $A\sigma$ with $A$ the MD cell cross section.

The surface charge of a NaCl (111) plane in our model system is $A\sigma_{0}=16e$. The theoretical compensating charge (Tasker rule) is therefore $8e$. Fig. 4(a) shows that the compensating charge $A\sigma$ 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_{\text{H}}$ of 8.23 $\mu$Fcm*-2* for the polar (111) surface.

8.2 $\mu$Fcm*-2* is a relatively modest capacitance, not much larger than the 4.4 $\mu$Fcm*-2* we found in Ref. Zhang and Sprik, 2016a 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. Zhang and Sprik, 2016a. 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_{\text{H}}$ with a value 4.23 $\mu$Fcm*-2*, effectively identical to the capacitance obtained in Ref. Zhang and Sprik, 2016a.

III.3 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 $\sim$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_{\mathrm{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 WolfWolf (1992) (see also NogueraNoguera (2000); Goniakowski, Finocchi, and Noguera (2008)). 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 MgOSpagnoli, Allen, and Parker (2011) and NiOFreund (1997).

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 conditionsRadenović et al. (2006). The force field used in the present study (the Joung-Cheetham model of Ref. Joung and T E Cheatham, 2008) 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.2 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 accountNoguera (2000). 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.

Acknowledgements.

TS is supported by a departmental studentship (No. RG84040) sponsored by the Engineering and Sciences Research Council (EPSRC) of the United Kingdom. The Research fellowship (No. ZH 477/1-1) provided by the German Research Foundation (DFG) for CZ is gratefully acknowledged.