Using molecular dynamics to quantify the electrical double layer and examine the potential for its direct observation in the in-situ TEM

Understanding the fundamental processes taking place at the electrode-electrolyte interface in batteries will play a key role in the development of next generation energy storage technologies. One of the most fundamental aspects of the electrode-electrolyte interface is the electrical double layer (EDL). Given the recent development of high spatial resolution in-situ electrochemical fluid cells for scanning transmission electron microscopy (STEM), there now exists the possibility that we can directly observe the formation and dynamics of the EDL. In this paper we predict electrolyte structure within the EDL using classical models and atomistic Molecular Dynamics (MD) simulations. Classical models are found to greatly differ from MD in predicted concentration profiles. It is thus suggested that MD must be used in order to accurately predict STEM images of the electrode-electrolyte interface. Using MD and image simulation together for a high contrast electrolyte (the high atomic number CsCl electrolyte), it is determined that, for a smooth interface, concentration profiles within the EDL should be visible experimentally. When normal experimental parameters such as rough interfaces and low-Z electrolytes (like those used in Li-ion batteries) are considered, observation of the EDL appears to be more difficult.


Background
Electrical double layer capacitors (EDLCs) are now being explored as energy storage devices for modern electronics such as memory-backup systems and electric vehicles [1][2][3]. These devices have power density capabilities that exceed conventional energy systems. Charge/discharge cycles can be performed very fast, and EDLCs can withstand more than 1 million operational cycles [1,3]. EDLC behavior is different from that found in conventional electrostatic capacitors due to the nature of the electrical double layer (EDL) that forms at the interface between charged electrode and electrolyte, resulting in ion space charge-based energy storage [1,[4][5][6][7][8][9]. Development of EDLC technology necessitates an understanding of EDL structure as a function of electrode and electrolyte (e.g. electrode surface defects can increase capacitance by breaking down electrolyte solvation shells [2]). EDLs influence capacitor properties including maximum operating voltage and energy density, ion conductivity, operational temperature, self-discharge rates, and safety [1][2][3]10,11].
EDL structure has not been directly observed in experiment, and theoretical techniques are often used to assess it. Classical theory can be used to approximately characterize the structure of the EDL [4][5][6]. As shown in Figure 1, classical models often assume a multilayer electrolyte structure at the electrode surface involving rigid counterion packing followed by a diffuse ion layer. Structure is then mathematically modeled with simple electrolyte concentration curves as a function of distance from the electrode surface. As suggested with Molecular Dynamics (MD) simulations, however, the EDL structure is instead a complicated function of atomistic environment [7][8][9]. The reason for this difference is that classical models generally do not handle fine details of electrolyte structure, not accounting for interactions at the electrode surface between solvent, solute, and electrode. As MD addresses these details, it can be considered preferable for modeling EDL structure.
As microscopic techniques develop for application to electrochemical systems, the possibility of direct observation of the EDL is rising [12]. Recent developments for in-situ transmission electron microscopy (TEM) have made it possible to analyze features of the electrodeelectrolyte interface. Thus far, in the application of insitu TEM to electrochemical problems, degradation products at the interface during battery cycling have been analyzed [13], and ionic concentration gradients have been observed [14]. Efforts are now being made to directly observe the EDL, and image simulation techniques can help identify the potential to observe it.
In this paper, we show that understanding structural features and image contrast features of the EDL generally requires an intensive simulation procedure combining MD and image simulation techniques. Through the use of MD, we recreate the interface of a 1.5 M CsCl solution (expected to have high Z-contrast) at a gold electrode under various electric field strengths. Then, using the MD-obtained structural information, million-atom fluid-stage models are constructed and used as input for image simulation. Image contrast features are found to be a function of the atomistic structure of the electrolyte and are thus influenced by the use of classical models or MD in structure generation. Observation of the cesiumbased EDL at a smooth electrode surface, under standard microscope operating conditions, is predicted to be possible at high field strength; however, for low field strength, low Z-contrast electrolyte, and/or high electrode roughness, observation of the EDL becomes more challenging.

Molecular simulation
The electrode-electrolyte interface is the main focus of this work, and its structural description must be carefully considered. In order to construct a full, atomic model of the interface for use in image simulation, it is necessary to characterize its structure through an analysis technique such as atomic packing density profiles (aka probability density distributions, or PDDs). Molecular Dynamics is performed here in order to determine the PDDs as a function of applied field.
Molecular Dynamics is carried out with an in-house code (certain publications were particularly helpful to the author in writing this code [15][16][17]) for 1.5 M aqueous CsCl solution interacting with a negatively charged gold electrode surface. Our approach is largely inspired by the wealth of electrode-electrolyte interface simulation studies [7][8][9][18][19][20][21][22][23][24][25]. Simulations are performed in the NpT ensemble at 298 K using the velocity-rescaling thermostat (coefficient = 0.1 ps) [26] and the barostat described below, a 1 fs timestep, a 10 ns (for field strength < 2.5 V/Å) or 30 ns (for field strength > 2.5 V/ Å) equilibration period, and a 10 ns production period. The system is two-dimensionally periodic with fixed cell dimensions of 3.1 × 3.1 nm in the interfacial plane. Along the remaining dimension, 4.1 nm of solution (containing 1074 water molecules, 30 anions, and 30 cations) is encapsulated between two walls that interact with the solution atoms.
One wall represents the cathode surface and is held fixed. In all cases except for zero-field simulation, it has an effective charge simulated through an applied electric field (i.e. electrode charge is not explicitly included; see electrostatics description below) [7], and this charge is equal and opposite to the charge of excess cations in solution (added by substituting water molecules). Obtaining structural features of the electrolyte solution at the cathode is the purpose of the simulation. The opposing wall, on the other hand, is effectively the anode, but only serves to encapsulate the system. No valuable information would be obtained by including anode charge in simulation, so this is not performed. Structural information within 1 nm of the anode surface is discarded, and we expect that the anode surface has no significant effect Beyond the first few atomic layers at the surface, the electrolyte has bulk-like structure. Note that concentration profiles calculated with molecular simulation are found to be significantly different.
on structural features across the remainder of the simulation box. The anode is an excellent candidate for use in a piston-like control of system pressure. After the first 10 ps of simulation (i.e. after initial relaxation), at every timestep, the anode is displaced by an amount proportional to the pressure acting upon its constant surface area. Exact displacements (in Å) are determined by multiplying the net two-body force acting on the wall (in eV/Å) by the simulation timestep (in fs) and a factor of 0.002. Note that a typical wall displacement per timestep is 0.0003 Å. External pressure effects are assumed to be insignificant and are not included. Pressure relaxation thus occurs by dynamic equilibration of the wall position, allowing the density at a given voltage to equilibrate. The chosen water model is the flexible SPC/Fw model [27]. Lennard-Jones (LJ) parameters for ion-ion and ionwater interactions are from the literature, and are used along with Lorentz-Berthelot mixing rules [28]. All wallatom interactions (i.e. cathode/anode-electrolyte interactions) are modeled by a 9-3 potential [7]. For the anode, a weak wall-atom interaction taken from a previous work is used for all non-hydrogen atoms [7]. Hydrogen atoms do no interact with the anode. For the cathode, parameters (see Table 1) are based on a DFT study on water monomer adsorption at the gold {111} surface [29]. Parameters are assigned in order to best reproduce the monomer binding energy (−0.13 eV), the metaloxygen bond length (0.302 nm), the tilt angle (~13°, near flat), and the water molecule's rotational barriers. We assume that these parameters are qualitatively sufficient for representing both the {111} and {100} facets of gold (see a comparative study that shows this to be reasonable for copper [30]; unfortunately, studies that consider both facets within the same work are rare which limits our ability to verify this assumption for gold [31]). The use of this simple model in simulation of the water-Au {111} interface gives generally good agreement with oxygen atom density profiles recently determined through dispersion-corrected ab initio molecular dynamics [32]. There is one exception to this agreement, however, as the ab initio simulations predict that the primary water layer at the gold surface is characterized by two density profile peaks as opposed to a single peak. Our model reproduces the strong peak from this study well, but it is too simple to reproduce the second, weak peak closer to the surface. Ion-wall interactions are determined in an approximate fashion by mixing united-atom O-wall parameters (i.e. the above O-wall parameters scaled as to have the full binding energy of −0.13 eV), the above mentioned water-water interaction, and ion-water interactions. Mixing rules are: ε ion-wall = ε water-wall * ε ion-water / ε water-water ; σ ion-wall = σ water-wall + (σ ion-water -σ waterwater ), where ε is the potential energy well depth and σ is the distance at zero energy. Atom-atom and wall-atom LJ interactions are cutoff in simulation at 1 nm, using the shifted-force form [33].
Shifted-force electrostatics has been shown to give reasonable accuracy, especially for structural details, if the cutoff is sufficiently large; force form is that in [33][34][35][36][37][38]34]. The electrostatic cutoff is set to 1.55 nm (i.e. one-half the box size along the periodic dimensions), which is assumed to describe well the electrostatic interactions taking place in a concentrated, 1.5 M CsCl solution (Debye length < 0.3 nm [9]). Thus, for this system, it is assumed that this electrostatic treatment is very accurate along the periodic dimensions of the cell (i.e. within the interfacial plane). Along these dimensions, the charge distribution is isotropic. Interface physics along the interface normal, however, may not be correctly treated [35,37]. Charge distribution is not isotropic along this direction, posing a problem for short-range electrostatic methods. Also of concern is that the electric field created by the cathode surface is applied along the interface normal as a long-range, distance-independent force: F field = q*E, where q is the atom charge and E is the uniform electric field [7]. If only shifted-force electrostatics were used, the shielding of that electric field would not even remotely be accurate (e.g. an ion 3 nm away from the cathode would not feel shielding ions and would incorrectly feel the full electric field). Therefore, we add a corrective treatment of electrostatics specifically along the direction of the interface normal that enables charge anisotropy and shielding effects to be handled accurately. The procedure for correction is outlined in Additional file 1.

Model construction
Virtual electron microscopy requires the construction of a detailed atomic model [39,40]. It is thus necessary for structural models of the amorphous silicon nitride window membrane, electrolyte solution, and gold electrode to be known. The combination of these components yields a million-atom fluid-stage model on which image simulation may be performed.
An amorphous silicon nitride (a-Si 3 N 4 ) window (~1 million atoms in size) at a target density of 3.14 g/cm 3 is created through a random, stoichiometric distribution of silicon and nitrogen atoms [41]. Atoms are positioned under the following spatial constraints that are intended to prevent unrealistic atomic potential overlap: d Si-N > 0.14 nm, d Si-Si > 0.23 nm, d N-N > 0.21 nm (values are based roughly on structural expectations [41], but are reduced for algorithm speed). The locally random nature of the atom distribution is assumed to be insignificant to image simulation results [42,43]. All image simulations use the same window. Representative electrolyte solution is derived from the PDDs determined with Molecular Dynamics in section 2.1. It is necessary to use a procedure by which the results for the small Molecular Dynamics box can be used to create a much larger structural sample. First, along the effective interface normal of the structural sample, the system is divided into 0.01 nm slices. Each slice corresponds to a bin of the PDDs. Next, atom counts of each PDD bin are scaled to account for the increase in size between the original simulation box and the large sample here constructed. Then, these desired atom counts are used to position a net number of atoms at random within each slice while following the spatial constraints: d H-H > 0.13 nm, d H-X > 0.13 nm, d X-X > 0.26 nm (where H=hydrogen atom, X=non-hydrogen atom). Through this procedure, an "atom soup" is created that captures the change in atom densities along the interface normal. This procedure is applied out to a distance of 1.5 nm along the interface normal as to fully generate the gold-solution interface. All of the remaining space (i.e. well beyond the interface) is filled with bulk solution-like atom densities, which are determined by taking the average atom densities between 1.5 and 2.5 nm in the PDD. As with the silicon nitride membrane construction, we assume that the loss of fine, local structure resulting from the random atom placement is inconsequential to image simulation.
Creating the gold electrode starts with cutting a large region from a bulk fcc gold lattice such that the resulting facets belong to the {001} plane family. The isotropic Debye-Waller factor, taken to be 0.016 nm 2 , is used to thermally displace the atoms of the gold electrode, creating a snapshot of one thermal phonon [44]. We are interested in how varying the electrical double layer structure (over different applied fields) changes image contrast features. Modeling additionally how gold's thermal phonons cause small contrast variations in the calculated image would complicate analysis and greatly increase computational effort while not providing valuable information. Therefore, a single snapshot of the gold electrode is used in all simulations as to reasonably simplify both computational effort and image analysis.
Once the complete fluid-stage model is built, it contains about 2,000,000 atoms and consists of a 10 × 10 × 49.7 nm silicon nitride membrane positioned 0.3 nm above a 5 × 10 × 100 nm electrolyte solution next to a gold electrode of equal dimensions. The (001) surface of the gold electrode is located along the silicon nitride membrane. It would be a trivial task to complete this fluid-stage holder model by attaching the remaining fluid layer and additional silicon nitride membrane to the bottom of the model. This step is unnecessary, however, as essentially the same information can be obtained in simulation without including atoms located below the specimen of interest. As discussed in section 2.3, it is simple to correct simulation results for the missing influence of these atoms.

Image simulation
Aberration-corrected HAADF-STEM image simulations are performed with the QSTEM image simulation software (V2.22) [40]. Our simulations follow a similar procedure as in our previous work [45]. Microscope parameters are based on the 200 kV JEOL 2100F AC-STEM microscope. These parameters are C s = 0.01 mm, C c = 1.8 mm, ΔE = 1.5 eV, α = 30 mrad, and HAADF inner detector angle = 100 mrad. This inner detector angle is roughly chosen as to be large enough to ensure thermal diffuse scattering is the dominant signal contribution and small enough to ensure sufficient signal strength. The outer detector angle, limited by the maximum scattering angle in the simulation, is 167.2 mrad. Additional, approximate parameters are C 5 = 3.2 mm and source size = 0.1 nm. The defocus plane is set to the Scherzer defocus of 6.1 nm below the top atomic plane of the gold electrode. We assume a brightness of 1.5 × 10 10 A/cm 2 sr and pixel dwell time of 2 μs. Slice thickness is set to a low value of 0.111 nm to better include the effects of electron beam interaction with the gold atoms; a highly quantitative study may desire even thinner slices. The simulation window size is set to 6.0 × 6.0 nm (i.e. 1200 × 1200 pixels with 0.005 nm point sampling). Thermal diffuse scattering is included in the calculated image, but only one structure is simulated (see section 2.2). We note that image simulation does not include the influence of the applied electric field on beam propagation. Based on the thickness of the electrode and applied field strengths utilized in this study, we generally expect the beam scattering due to the field to be less than 1 mrad, an insignificant effect compared to general resolution concerns.
Statistical analysis of simulated images is performed by taking a pixel-averaged line scan along the interface normal (40 pixels averaged per data point). A correction factor must be applied to these line scans in order to account for the fluid and window membrane located beneath the electrode interface which have not been included in our image simulation model. In the simulations of our previous work [45] as well as in electrode-electrolyte simulations in which these atoms are instead included, we have considered how~400-500 nm of fluid and a 50 nm silicon nitride window located below the specimen of interest influence image features. It has been found that the effect on the image is to reduce the specimen signal-to-noise (SNR) by a factor of~1.2. Thus, we apply this SNR reduction factor to all results as to statistically include the influence of the postspecimen environment without arguably wasting computational effort. Note that longer fluid path lengths would require the application of larger correction factors (thus resulting in smaller SNR values).

Additional models for image simulation
A simplified and low-effort approach for interface model construction (i.e. for image simulation) would be to use basic analytical models, instead of MD, to describe atomic density profiles. In order to assess this approach, we consider two analytical models for the cesium ion distribution and compare image simulation results with MD-based results (believed to be more accurate). Interface models are built using Gouy-Chapman (GC) theory [4,5] or a modified Helmholtz (MH) theory [6]. The Helmholtz model has been found to be closer than GC to our MD results (see section 3.1), and we here slightly modify the use of the Helmholtz model in hopes of achieving good agreement with MD-based results. GC and MH model construction generally follows the same procedure as in Section 2.2 with the exception that probability density distributions are determined analytically. For both models, water density, based on MD results, is constant at 31.2 molecules/nm 3 . At distances less than 3 Å from the gold surface, there are no oxygen or hydrogen atoms present. Chloride density is assumed to be constant at 0.895/nm 3 . For the basic GC model, its density is also set to zero at < 3 Å, and for the MH model (which assumes that the neutral surface repels ions, see Figure 2), its density is set to zero at < 5.9 Å. Cesium density follows the same rules as chloride. In addition, however, at nonzero voltage, cesium density is increased as a function of distance from the gold surface in accordance with the theoretical model considered. For GC, cesium density is increased according to an exponential function (Debye length = 2.5 Å, as CsCl concentration is 1.5 M) such that the total charge of cesium ions added is equal to the amount of opposite charge on the electrode surface. Chloride repulsion effects are assumed to be insignificant. For MH, a sharp, linear decay of cesium density is added over 3 to 3.4 Å from the gold surface such that the total amount of cesium ions added is equal to the amount observed at the gold surface (i.e. within the first atomic layer) in MD. Chloride ions are added, in the same manner, at 2 Å distance from the cesium layer as to balance excess charge created by cesium overshielding. We note that this modified approach uses knowledge from MD, but the cesium packing density is a parameter that could easily be varied without MD. We wish to see if the use of just MD cesium packing densities is sufficient to give image simulation results that match those obtained with detailed MD.
Influence of surface roughness on image contrast has also been considered in additional models. The method of introducing roughness is arbitrary and qualitative, as a highly detailed, accurate description of surface roughness is well beyond the scope of this study. The following procedure is meant to create a surface that is intermediate between flat and truly random. After the construction procedure of Section 2.2, the gold surface in contact with the solution is divided into 100 × 100 spatial bins. Half of these bins are chosen at random, and solution atoms, within the projection of these bins along the interface normal, are moved into the surface by a random distance (per bin) between 0 and 10 Å along that normal. For each affected bin, gold atoms within the bin's projection that are within the random distance minus 1.5 Å from the gold surface plane are removed as to make space for the solution atoms. Note that, for comparative purposes, the exact same surface roughening (which is randomly generated once) is used for each constructed model. surface charge densities ranging from 0 to −2.08 e/ nm 2 ). Applying a negative field creates an electric double layer of solute ions and repels solvent molecules from the electrode surface. These structural effects are in turn significant to predicted features of the EDL image.

Results and discussion
Counterion packing at the electrode surface shows complex features dependent on the electric field regime. Integrated cesium ion densities (determined by integrating the ion population within 5 Å of the electrode surface; see Figure 2) are, in order of increasing electric field strength, 0.01, 0.04, 0.39, 1.20, 2.40, 3.02, and 2.91 ions/nm 2 . Cesium ion/absolute electrode surface charge ratios go as 0, 0.40, 1.26, 1.92, 2.31, 1.93, and 1.4. Three regimes can be identified. The first, at low electric field, shows little cesium packing at the electrode surface due to the field strength being insufficient to pull cesium ions out of their solvation shell. The second, at intermediate electric field, shows so much cesium packing that overshielding of electrode charge occurs, even resulting in a cesium ion/negative electrode surface charge ratio greater than 2 at −1.88 V/Å. The third, at high electric field, shows saturation of cesium packing, and this is best demonstrated by the decrease (as opposed to expected increase) in cesium ion density when changing the field from −2.82 to −3.77 V/Å. Also, an increase in field strength results in a more rigid cesium layer at a closer distance to the electrode surface (peak maxima range from 3.4 to 2.8 Å).
Many other structural features of the EDL can be seen in the PDD. Due to strong layering effects at the gold surface, there is a region of excluded-volume that repels cesium ions and has a field-dependent width. Beyond this region, particularly at field strengths lower than 1.13 V/Å, there is generally one cesium peak found between 6 and 8 Å. Considering the anion PDDs, between 4 and Figure 3 Oxygen PDD vs. electric field. Interfacial water structure forms at a density and location dependent on the electric field. Beyond 12 Å (not shown), density oscillations are insignificant and water density approaches that of the bulk. 7 Å, one or two chloride peaks are found. Overshielding by the cesium layer at the surface is charge balanced by chloride ions located in this region. Nothing else in the ion PDDs significantly deviates from bulk-like electrolyte behavior.
Application of the electric field results in reduced solvent density and distorted solvent structure (see Figure 3). For all fields, the oxygen atom PDD is characterized by oscillating density features out to 12 Å, physically due to the interaction of water with the flat surface. Integrated water densities at the surface (determined by integrating within 5 Å) are, in order of increasing electric field strength, 11.3, 11.1, 10.3, 8.8, 6.3, 5.4, and 6.4/nm 2 . Increasing field strength also dampens the second water peak at 6 Å. At zero-field, water molecules lie nearly flat along the metal surface. Increasing field strength reorients hydrogen atoms towards the metal surface, restructuring the hydrogen bond network. Note that, for this study, the Z-contrast of the cesium ion is much greater than that of oxygen. For other ions having lower Zcontrast, such changes in solvent density and structure can be expected to possibly dominate image contrast.

EDL contrast from analytical models and their shortcomings
Image simulation is performed for models of the electrode-electrolyte interface contained within an insitu TEM holder, with models built through either Gouy-Chapman (GC), modified Helmholtz (MH), or MD density distributions (see Figure 4 for example images). MD is found to predict better observability of the EDL than GC and MH (see Figure 5). Contrast peak position, magnitude, and width greatly differ between these three approaches.
Image contrast predicted by the GC model is diffuse, just as would be expected considering the model's emphasis on a diffuse layer of counterions. In addition to having different descriptions of counterion distribution, MH and GC ignore low-density regions (and other structural features) related to the solvent. Therefore, the overall mass distribution is possibly incorrect when using these models. The classical models also lack fine detail of the field-dependent solute structure. These shortcomings explain differences in predicted contrast features between classical models and MD. As can be concluded from such model-dependent predictions, obtaining correct contrast features requires either the use of accurate molecular simulation or more detailed analytical models that include detailed EDL structural properties.

EDL contrast from MD
Image simulation is performed for in-situ electrodeelectrolyte interfaces, using MD to describe interface structure. EDL structure-contrast relationships are determined, and the observability of the EDL is assessed. Additionally, analysis of EDL contrast is shown to be especially difficult when low Z-contrast electrolyte is instead considered.
Integrated line scans for image contrast, determined by subtracting zero-field signal from applied-field signal, show sharp contrast peaks at 3 Å corresponding to the shielding layer of Cs + ions (see Figure 6). Contrast maxima intensity is nearly a linear function of cesium packing density (i.e.~300 electron counts per cesium ion/nm 2 ). Based on the complicated structure of chloride and water at the gold surface present in the interface model, it can be concluded that such a simple, linear relationship is only likely for a high Z-contrast cation. In changing the field from −2.82 to −3.77 V/Å, EDL contrast slightly decreases. This is due to the saturation of the double layer in this field range (i.e. a lack of increase in cesium content prevents any contrast increase). Signal-to-noise (SNR) of the maximum contrast in the line scan is, in order of increasing field strength, 0.2, 1.1, 2.9, 5.7, 7.8, and 7.2 (per pixel). SNR is a near-linear function of cesium packing density, with the R 2 value of the linear fit being 0.997. The SNR values are ideal in that they assume Poisson noise in each image. It has also been assumed that no time-dependent ion diffusion effects, which would complicate image interpretation, occur. Rose criterion is thus met, in this ideal treatment, for field strengths higher than 1.7 V/Å. At weaker fields, observation of the EDL will be hindered by noise. There is an additional shoulder peak at 6 Å due to EDL structure, but SNR is always low for this peak. No other contrast features are found.
When Cs + ions are replaced with Li + ions (by pure substitution to the structure used in image simulation; MD not performed), the structure-contrast relationship is transformed (see Figure 7). A contrast fringe pattern appears near the surface, but the relationship between fringe contrast and field is not clear. Also, SNR is low Figure 6 Integrated contrast line scan vs. electric field. Two contrast features due to EDL structure are found at 3 and 6 Å. Peak maxima, at 3 Å, scale near-linearly with cesium packing density. The gold surface is located at the origin. MD-determined results are shown. for these peaks. At a field of −1.88 V/Å, summation of all contrast over the line scan gives a value that is only 2% of that found for the equivalent Cs + case. Interpreting contrast, for lithium chloride, appears to be much less straightforward.

Surface roughness and contrast effects
Image simulation is performed in the approach of Section 3.3 with the inclusion of surface roughness in the simulation model. Roughness causes the major contrast peak to essentially split into several weaker peaks (see Figure 8). At applied fields of −0.56, −1.13, and −1.88 V/ Å, the SNR of the major peak is consistently reduced bỹ 60%. This loss of major peak signal is likely due to (1) redistribution of cesium ions across the surface and (2) loss of atomic order at the interface, which reduces electron channeling effects. Minor peaks, for the tested fields, have SNRs less than unity. Therefore, surface roughness results in images with only low-contrast features, limiting our ability to observe EDL structure.
A more realistic roughness model would have to additionally include complicated structural features caused by the defects at the rough surface, and this would make assessment of the structure-contrast relationship more  difficult. Therefore, surface roughness can be expected to also reduce our ability to understand images of EDL structure.

Conclusions
Optimizing performance of electrical double layer capacitors requires understanding of the double layer structure that defines their operation, and being able to image the double layer with an electron microscope will yield that understanding. Image interpretation, nonetheless, will not be straightforward, and we have demonstrated simulation methodology that allows for more accurate interpretation of contrast features. Molecular Dynamics and image simulation have been used to characterize the electrical double layer of 1.5 M CsCl solution at a gold electrode, under various applied fields. By comparing MD and image simulation, the following has been deduced: (1)Electric fields induce various structural changes at the electrode-electrolyte interface. Cation density is a non-linear function of electric field (and therefore non-linear with electrode surface charge) due to effects such as overpacking and excluded-volume. At high field strength, cation saturation of the interface prevents an increase in cation density. Also, chloride density generally correlates with the extent of cesium overpacking. Solvent is generally repelled from the interface by the field with the remaining solvent molecules adapting new orientations. These structural features are not well captured by simple theoretical models of the electrical double layer. (2)Simple theoretical models of the electrical double layer, when used to perform image simulation, do not give the same image contrast features as those obtained through molecular simulation. Simple models tested in this paper actually predict poorer image quality than that found through MD. We believe MD results to generally be more accurate than simple model results as MD is able to capture fine structural detail. Therefore, in order to accurately perform image simulation on interfacial systems, molecular modeling or more detailed analytical models must first be used. (3)When imaging a high Z-contrast electrolyte under high field strength and ideal conditions (i.e. a smooth electrode surface and only Poisson noise), EDL contrast features can be observed and interpreted. SNR is found to be a near-linear function of cesium packing density (i.e.~2.5 SNR per ion/nm 2 ).
In general, however, noise greatly limits observation capabilities. Low Z-contrast electrolytes show both poor SNR and contrast features that are difficult to relate to interface structure. (4)Surface roughness greatly reduces SNR of the EDL, making EDL observation and interpretation much more difficult. If the electrical double layer is to be imaged, it is important to use an electrode with low surface roughness.
From all of these findings, it can be concluded that direct observation of the EDL is most likely possible for a high Z-contrast electrolyte under high electric field at a smooth electrode surface. Achieving observation of the EDL will lead to a deeper understanding of the operation of electrical devices, but concerns put forth in this paper must be carefully considered in order to accomplish it.

Competing interests
The one competing interest in this work is that the co-author NB is the editor-in-chief of this journal. We assume the review process will address any concerns of bias sufficiently and will assert the merit of this work. Aside from that, the authors declare that they have no competing interests.
Authors' contributions DW designed and carried out the methodology used in this work by combining approaches of past simulations with creativity as to meet this paper's specific demands. LM provided criticism, helped with project design, and has been working tirelessly to produce experimental validation of the simulations. HH created the graphic figure shown in Figure 1. RF ensured simulation method reliability was well tested and filtered out undesirable ideas proposed for methodology. JE provided criticism and future direction, particularly with matters of comparison between theory and experiment. NB advised this work from start to finish. All authors approve of this manuscript.