Stern and Diffuse Layer Interactions During Ionic Strength Cycling

Second harmonic generation amplitude and phase measurements are acquired in real time from fused silica:water interfaces that are subjected to ionic strength transitions conducted at pH 5.8. In conjunction with atomistic modeling, we identify correlations between structure in the Stern layer, encoded in the total second-order nonlinear susceptibility, chi(2)tot, and in the diffuse layer, encoded in the product of chi(2)tot and the total interfacial potential, phi(0)tot. chi(2)tot:phi(0)tot correlation plots indicate that the dynamics in the Stern and diffuse layers are decoupled from one another under some conditions (large change in ionic strength), while they change in lockstep under others (smaller change in ionic strength) as the ionic strength in the aqueous bulk solution varies. The quantitative structural and electrostatic information obtained also informs on the molecular origin of hysteresis in ionic strength cycling over fused silica. Atomistic simulations suggest a prominent role of contact ion pairs (as opposed to solvent-separated ion pairs) in the Stern layer. Those simulations also indicate that net water alignment is limited to the first 2 nm from the interface, even at 0 M ionic strength, highlighting water's polarization as an important contributor to nonlinear optical signal generation.

To provide an estimate of what atomic structures recapitulate the experimentally determined ($) values, we combine experiment and atomistic insights from molecular dynamics simulations to follow how interfacial structure and total potential vary as we transition an aqueous solution over fused silica between various concentrations of NaCl while maintaining the bulk solution pH at 5.8. We provide concrete evidence that the dynamics in the Stern and Diffuse layers are decoupled from one another under some conditions, while they are strongly coupled under other conditions that are readily identified. Furthermore, we obtain quantitative structural and electrostatic information that informs on the molecular origin of hysteresis in ionic strength cycling over fused silica. 38,39,70 In the experiments (see Methods and Supporting Information Section S1-S4), we apply heterodyne-detected second harmonic generation (HD-SHG) 9,41,42,71 to evaluate the following relationship describing the second-and third-order mixing with the interfacial potential: Here, the SHG amplitude, %&' , and phase, %&' , yield the total second-order response, !"! ($) . We use a potential that decays exponentially from the surface into bulk water. 42,46,49,54,72 In this case, the DC phase angle, 34,637 , is given by arctan( : 3 ), with ∆ : being the wavevector mismatch of the optical process (1.1 x 10 7 m -1 in our case) and lD being the Debye-Hückel screening length in the bulk ionic solution. 41 At low ionic strength, 3 ≈1 x 10 -7 m and 34,637 → ; $ , whereas at high ionic strength, 34,637 → 0. 41 We recently reported a newly identified term, 9 (,) , which is ~1.5 × /0!12 (,) , 9 where, for off-resonant SHG, /0!12 (,) = (9.6 ± 1.9) × 10 <$$ m $ V <$ from experiments 65 and 10.3 × 10 <$$ m $ V <$ from quantum mechanical calculations. 73 Results. Fig. 1A shows results for an experiment in which the salt concentration is first lowered quickly from 0.1 M to 10  changes approximately as much as in the previous scenario (Fig. 1B). Finally, jumping from 1 mM to 10 μM produces a comparably smaller change in the SHG amplitude relative to the jumps starting at 0.1 M, whereas the phase changes by about 20° (Fig. 1C).
Next, we obtained point estimates of the total surface potential and the second-order nonlinear susceptibility of the interface using the following expressions: 9 cos $ 1 34,637 2 (2b).
We calibrate the SHG response from a given aqueous:solid interface using a vertically aligned piece of z-cut a-quartz employed as an external IEEE phase reference standard 74,75 with GHFI,J ($) = 8 × 10 <K, m V -1 in place of the water while properly accounting for Fresnel coefficients. 43,62 This procedure yields a calibration and referencing ratio, C/R, of 3.6 × 10 <$$ $ <K in our spectrometer. 9 The results are shown in Figure 2. In the first five minutes, the 0.1 M to 10 µM jump results in a continuous change in the surface potential to increasingly negative values, as one would expect. In contrast, ($) first rises, reaches a maximum at ca. 3 minutes, and then decreases to a constant value at longer times. This time scale is comparable to what was recently reported from time-resolved X-ray reflectivity measurements of ion exchange between the inner and outer Helmholtz plane over mica: water interfaces. 15 The discontinuity in the ($) and the (0) !"! values that occurs at ca. 7 minutes is due to the in-line fast conductivity meter reaching its sensitivity limit (please see Supporting Information Section S2) at that time. Off-line measurements of the 10 µM eluent conductivity performed using a more sensitive conductivity probe that however requires a much larger sample volume show the conductivity at that point should arrive at approximately 2 µS cm -1 in the flow cell. Absent a reliable interpolation model, the discontinuity in conductivity results in a discontinuity in the DC phase angle used in eqn. 2 and the the ($) and the (0) !"! point estimates. The return jump is quick, with the changes in the ($) and (0) !"! values not resolvable using our existing time resolution of 12 seconds.
The 0.1 M to 1 mM jump shows the expected decrease in surface potential along with a monotonic increase in ($) until ca. double its starting value. The return jump results in quick return to the starting values, just like in the initial jump. Finally, the 1 mM to 10 µM jump results in the expected change in surface potential to more negative values, while ($) undergoes a brief small increase followed by a slightly smaller value than what is observed at the start. Hysteresis is not observed in any of the three ionic strength cycling experiments, for a given experiment i.e. the starting and final ($) and (0) !"! point estimates obtained upon completion of the return jump are comparable to one another within error (5 percent). 9 We find that the ($) and (0) !"! values determined for the 1 mM condition depend on whether we start the experiment at 1 mM or whether we jump to it from 0.1 M or from 10 µM.   3B). As the current consensus in the field is that ($) and (0) !"! report on structure and dynamics in the Stern and diffuse layers, respectively (vide supra), we conclude from Fig. 3A and 3B that the two do not necessarily act in concert, depending on the path by which one changes the ionic strength.
We then asked whether a surface in contact with low ionic strength water has a different interfacial structure depending on whether it had been pre-exposed to 0.1 M NaCl or not. Fig. 3C presents the point estimates of ($) and (0) !"! for a freshly cleaned a hemisphere that, after mounting it on the flow cell, was exposed first to 10 µM water held at pH 5.8 Then, the ionic strength was jumped to 0.1M at the same pH (the measured amplitude and phase data are shown in in Supporting Information Figure S7B). Fig. 3D displays the corresponding ($) : (0) !"! correlation plot for this ionic strength transition, along with the one in which the surface had been Discussion. The experimental results presented here quantify how interfacial structure and electrostatics at fused silica:water interfaces depend on how ionic strength is varied, even as the bulk solution pH is held constant. When in contact with low ionic strength aqueous solutions, the interfacial potential is shown to depend on whether or not the surface had been pre-exposed to salt, while the interfacial structure encoded in ($) does not. We therefore find that dynamic changes of interfacial structure in the Stern layer (reported by ($) ) and the diffuse layer (reported by the The jumps from high to low ionic strength should involve Na + ion desorption. We therefore constructed several plausible interfacial models having a varying proportion of SiOH:SiO -:Na + :H2O:NaCl (Fig. 4A). To simulate surfaces in contact with 0.1 M salt, we placed water containing 0.1 M NaCl in contact with a silica surface having SiOgroups to which Na + cations are coordinated as a direct contact ion pair (CIP). The 1 mM experiment is not feasible to simulate in general, so we employed instead an ion-free aqueous phase in contact with the same number of CIPs we modeled in the 0.1 M case. This model choice is motivated by the notion that the Debye length at 1 mM is around 10 nm, 85 resulting in a considerable number of ions at the interface. We also contrasted contact-with solvent-separated ion pairs (SSIPs). Finally, we simulated the 10 µM experiments using an ion-free aqueous phase in contact with the same negatively charged silica surface we modeled in the 0.1 M case that is, however, void of interfacial Na + cations (no ion pairs, "NIP"; charge neutrality is provided by countercharges deep inside the aqueous phase that are separated from the interface using a semipermeable boundary described in Supporting Information section S9). In the cases of the solvent-separated and no-ion pairs, semipermeable boundaries (please see Methods and Supporting Information) are included to prevent Na + ions from reaching the interface.
To further describe the CIP and SSIP cases, we present Fig. 4B-D, which show the Na + ion density and the SiO -⋯Na + distance probability densities as a function of SiO -⋯Na + distance from the surface, respectively. The results follow the expected trends, with the solvent-separated pair distance being further apart than that of the direct contact pairs. In addition to the Si-O -⋯Na + distance differences, the angle:distance correlation plots for the Si-O -⋯Na + ion pairs indicate a most probable Si-O -⋯Na + angle of 45° to 70° for the contact ion pairs, while the solvent-separated ion pairs show most probable Si-O -⋯Na + angles around 10°. In all three 0 M cases, the local water density remains the same.
We then computed  (Table I; to connect to the s-in/p-out polarization combination used in the experiments, we computed the zxx tensor element). Contributors to !"! ($) , such as the Na + , Cl -, and SiOgroups, are neglected. We instead limit our calculations to the contributions of the much more abundant SiOH groups and the water molecules, i.e. !"! ($) = L*MN ($) + /0!12 ($) .
The calculated /0!12 values (Table I, (Table I) 5). This figure also shows that the skewness, or asymmetry of the distribution (reported by the third moment), is close to zero at 2 nm distance and beyond. Despite the distinct water orientation distributions, the local water density remains the same in all three 0 M cases (Fig. 5). Taken together, the results indicate that the net alignment of the water molecules does not extend far into the bulk phase, even when the interface is void of any ions. We can then propose that the diffuse layer contribution to the SHG process, again reported by the (,) (0) !"! product, is likely due to the polarization of water molecules in the diffuse layer as opposed to their alignment. Table I shows that while the computational estimates qualitatively recapitulate the results obtained from the experiments, they are about 10 times smaller in !"! ($) when compared to the experiment, depending on charge density. Likewise, the total computed surface potential does not change by as much as it does in the experiment, even though it follows the expected trend (lower potential at higher [salt], or at lower change density). We attribute this mismatch to the lack of an electronic structure calculation in our all-classical MD trajectories, the simplifying assumptions in our largely rigid silica model, possible contributions to !"! ($) , from Na + , Cl -, and SiO -, and the use of a cubic crystalline bulk of inert atoms as opposed to amorphous bulk silica. Yet, the electrostatic potentials due to only the water molecules listed in Table I are comparable to those reported by Chen and Singer, 24 and span the range of the experimentally derived point estimates of (0) !"! .
Moreover, we find qualitative agreement between the experimental and molecular dynamics results regarding !"! ($) , as discussed next.
We first compare the !"! ($) and (0) !"! point estimates from our atomistic simulations and the experiments for the case shown in Fig. 3D (jump from ultrapure water to 0.1M [salt] on a fused silica surface that had vs had not been pre-exposed to 0.1M [NaCl], see "start" marks in Fig. 3D).
The latter case is modeled as "0M/NIP" in Table I, i.e. an ion-bare negatively charged silica surface in contact with pure water. The former case would be the "0M/CIP" case in Table I, i.e. a silica surface containing SiO -⋯Na + contact ion pairs, left behind from the prior salt exposure, that is in contact with pure water. Table I  We caution that unlike in our atomistic model, the interface probed in our experiments is unlikely to be entirely void of adsorbed ions at the lowest ionic strength examined here (10 µM).
Moreover, we caution that our idealized model neglects many aspects of the experiment, such as acid-base chemistry of the amphoteric SiOH groups, surface reconstructions, dissolved carbonate, protons, hydroxide ions, etc. In addition, the experiment employs fused silica, which is amorphous, whereas the model in the MD simulations is built upon a fixed lattice of inert atoms with a fictious mass. Finally, our water model choice (extended single point charge) is not polarizable. Despite these shortcomings, we are excited to lay out a path for connecting the structural and electrostatic information from the experiments to atomic structure. Likewise, the results presented here may serve for further analysis of existing atomistic models of silica:water interfaces, as !"! ($) and (0) !"! are straight-forward to obtain from already completed production runs or total energy/geometry optimization calculations.
Taken together, our combined computational and experimental approach opens a door to quantifying interfacial structure and electrostatics at charged, buried aqueous interfaces in real time. The approach allows one to devise experiments that tackle the chicken or the egg question as applied to charged systems other than oxide:water interfaces. For instance, one can now investigate if the net interfacial structure changes before or after (0) !"! if one applies an external stimulus, such as a electrostatic potential on an electrode. 88 Other avenues to explore concern corrosion 89 or transmembrane potentials. 90,91 Values of d !"! ($) /dt during the various stages of the process followed here can be in principle obtained from fits to coupled rate expressions describing intra-Stern (inner-and outer Helmholtz layer) and Stern-diffuse layer exchange dynamics, but this step requires reliable estimates for Na + and SiOsurface coverages for each time point, which have not yet been obtained for our experimental conditions. Future work includes speeding up the data acquisition time on our HD-SHG spectrometer and adding plug flow to our MD simulations to probe the dynamics of ion adsorption and desorption. Overall, studies on charged interfaces are thus likely to produce papers in The Journal of Physical Chemistry for another 125+ years.

Methods.
Our current HD-SHG spectrometer, described in detail in our previous work, records the SHG amplitude and phase every 12 seconds, which is a reasonable compromise given the maximum speed of our 0.1 M delay stage motor and acceptable signal-to-noise levels. The reference state conditions we applied were the following (see Supporting Information section S1 for the raw fringe data, fitting procedure, values of Esig and jsig, and referencing): the SHG amplitude obtained at 10 μΜ salt divided by the square root of the SHG signal intensity from the quartz calibration crystal put in place of water at 10 μΜ salt and the same pH, All molecular dynamics (MD) simulations were performed using the LAMMPS software package. 92 In all cases during both equilibration and production runs, the MD trajectories were integrated using the velocity-Verlet methods with a timestep of 1.5 fs. Rigid-body constraints for the water molecules and the terminal silanol groups were enforced using SHAKE. 93 The Nosé-Hoover thermostat (100 fs relaxation) and the Nosé-Hoover barostat (1000 fs relaxation) along xy directions were applied in all simulations to control the temperature (298.15 K) and the lateral pressure (1 atm). The simulation cell, measuring 2.77 nm by 2.88 nm in each direction, was periodically replicated in x and y coordinate space and filled using SPC/E water, 94 NaCl ions, 95 and model silica slabs. 71,96 Long-range contributions of Coulomb interactions were treated using a particle-particle particle-mesh method. 97 To prevent contributions from these interactions along z direction, a vacuum region was introduced on both sides of the simulation cell. 24,83,84 All the quantities reported here were averaged using simulation trajectories of 2-4 independent initial configurations over 30 ns after equilibration at least during 10 ns.
To generate the solid interfaces, two substrates were placed symmetrically in a simulation cell, with the center at z = 0 (Supporting Information Fig. S8). Each side was composed of four layers of neutral atoms previously described 71 terminated with silanol groups. 96  donors and acceptors to terminal groups. 48 The surface charge density of the silica surface is controlled by the number of deprotonated hydroxyl groups, ranging from 0 C m -2 (zero deprotonated silanol groups) to -0.08 C m -2 (4 deprotonated silanol groups).
A total of 3,000 SPC/E water molecules were inserted between two silica surfaces. The z positions of the silica substrates were adjusted to reduce pressure along the confinement to 1 atm.
The final distance between of oxygen atoms of the silanol groups in the z-direction was 11.54 nm.
After another short equilibration, NaCl ions were placed between the surfaces. Interactions between all atoms are described by Coulomb and Lennard-Jones (LJ) potentials without polarization. For cases involving a charged silica interface with no contact ion pairs, semipermeable boundaries were introduced to interact only with Na + ions, described by a truncated LJ potential. More details of the model and interaction potentials are described in Supporting Information.
We calculated the total mean electric potential, Φ -.
where is a running index for species, Y and Z are simulation box size along x and y coordinates, respectively, Δz (0.02 nm) is the grid size along z axis, and * and * are charge and z-position of th atom, respectively. The water contribution, Φ [\-(z), to the potential is calculated using a linear polarization 24,83,84 : where `i s the local field correction factor at frequency , and /0! ( ) is local number density of water molecule. The angle ( ) is calculated between z axis and a dipole vector of water molecule at z=z. The factor of -1 in front of cosine functions is included since the surface normal vector points from the aqueous region to the silica region. The values of the first hyperpolarizability are taken from Jansen et al. 99 as in Singer et al. 75 The susceptibility of water contribution is calculated as follows, by being integrated along z axis and normalized by the local field correction factor: where V represents the boundary between a SHG-active region and the bulk (a SHG-inactive region), and ) represents the boundary between ionic water and the silica surface. Here, V = 0 nm and ) = -5.77 nm, where oxygens of the silanol sites are placed.
Second, the contribution of SiOH to the susceptibility, χ a'bc ($) , in the same polarization is calculated following the same procedure for the χ [\-($) . To calculate , SiOH is placed in the zx plane. The OH bond is aligned along c axis, and b axis is orthogonal to the ca plane and parallel to y axis. The first hyperpolarizabilities in the non-resonant (NR) condition for the OH oscillator are calculated using the following relation, 47, 100, 101 : where β 0V] e is the first hyperpolarizability in the resonant condition, and β 0V] &-\-'B is its magnitude at the static limit. The harmonic-oscillator approximation is applied to calculate β 0V] e as follows: where bc is reduced mass of the OH oscillator, bc is the frequency of the oscillator, 0V (K) is ab element of the linear polarizability tensor, ] is c element of the dipole moment vector, and Γ bc is the dissipation from the environment. In this work, bc = 3000 cm <K , and Γ bc = (0.5 ps) <K . 102 Both derivatives of the linear polarizabilities and the dipole moment are obtained from Backus et al. 71 and Gaigeot et al. 103 Then, χ ()*+ (-) as a function of the fixed tilt angle ( / ) is given as follows: where a'bc is the number density of the hydroxyl groups at the silica surface. The azimuthal angle is found uniform so 〈cos $ 〉 = 〈sin $ 〉 = 0.5. As for χ [\-($) , the factor of -1 in front of cosine functions is included since the surface normal vector points from the aqueous region to the silica region, and the local field correction is not included.

Associated Content
Supporting Information: Heterodyne-detected SHG interference fringes, fitting procedures, referencing, SHG intensity data, flow-dependent data, additional computational details and methods.     Point estimates of χ (2) and Φ(0)tot as a function of time for 10 μM to 0.1 M jump on a fused silica substrate that had not (light circles) and had (dark circles) been pre-exposed to salt prior to jump.