Ionic structure and decay length in highly concentrated confined electrolytes

The behavior of electrolyte ions in liquids confined between charged macromolecules regulates many processes in soft materials such as colloids, emulsions, polymeric membranes, and proteins.^1,2 The local ionic environment can modulate the effective interaction between the charged macromolecules, thus changing their assembly behavior.^3,4 Many energy storage applications^5–9 and separation process technologies^10–12 based on these materials rely on the complex organization and transport of ions near macromolecular surfaces. Therefore, investigating the self-assembly of ions in the nanoconfinement created by material surfaces has been the focus of many experimental and theoretical studies.^13–25 

From a theoretical standpoint, the surfaces are often modeled as planar interfaces considering the large size difference between the ions and the confining macromolecules. In many systems under dilute electrolyte conditions, the electric field associated with the charged interfaces decays exponentially with distance in good agreement with the behavior predicted by the Debye–Hückel theory²⁶ where the characteristic decay length scales inversely with the square root of the electrolyte concentration. In conditions where correlations between ions are strong, e.g., for highly concentrated electrolytes or multivalent electrolytes, an accurate description of the ionic structure requires departure from the simplest mean-field theories. Approaches such as density functional theory and integral equations where correlation effects are explicitly included take into account the important contributions of the fluctuations to the thermodynamic potential in such systems.^{18,20,22,27–34} Computationally, the effects of correlations on the ionic structure in confinement are studied using coarse-grained models that treat the ions as finite-sized particles and replace the molecular structure of the solvent and material surfaces with dielectric continua.^18,35–48 Some studies have also probed the effects of the solvent structure on the organization of ions near interfaces using explicit-solvent models.^23,49–51 These investigations have shown that strong ionic correlations can have profound effects on the ionic structure, which can significantly alter the screening of the charged surfaces and change the character of the effective interaction between them.

In addition to the explicit methods of calculation, several phenomenological approaches, including modified Poisson–Boltzmann theories and heuristic constructions, have been useful in developing an intuitive understanding of the effects of correlations and discussing the results of simulations describing the structural organization of ions.^52–56 In the context of electrolytes and their interactions with macroions or charged surfaces, these include the notions of ion condensation⁵⁷ and overcharging,⁵⁸ the identification of forces as induced by ion correlations,⁵⁹ and the description of ion environments in terms of Wigner cells⁶⁰ or ionic glasses.⁶¹ 

Recent experiments have directly measured the behavior of the decay length associated with the electric field originating from the charged surface across a wide range of electrolyte solutions at high concentrations.^62–68 Based on surface force measurements, Perkin and co-workers^62,63 showed that the decay of the electrostatic force between two charged surfaces in a highly concentrated electrolyte solution is much weaker than the Debye–Hückel prediction. A similar “underscreening” effect was observed by Gaddam and Ducker via an alternate approach involving the measurements of the surface excess of fluorescein.⁶⁴ These experiments show evidence for the electrostatic decay or screening length that greatly exceeds the theoretical Debye length and rises with an increase in concentration c for high c values. This observation has been made for a wide range of electrolytes from simple salt solutions (LiCl, NaCl, and CsCl in water) to pure ionic liquids to ionic liquids dissolved in solvents, signaling the universal nature of the behavior of highly concentrated electrolytes.

These experiments have inspired many recent theoretical investigations of highly concentrated electrolyte systems. Efforts have been made to describe the properties of electrolytes and the non-monotonic behavior of the associated decay length using both analytical approaches^{33,34,69–73} and computer simulations.^72,74–76 Most studies have focused on understanding the non-monotonic behavior of the decay length by computing pair correlation functions in bulk electrolyte systems.^70,74 Some computational studies have focused on measuring the screening behavior of electrolytes in confined electrolyte systems.^42,72,75,76

In an earlier publication,¹⁸ we investigated the structure of electrolyte ions in the confinement formed near uncharged planar interfaces using coarse-grained molecular dynamics (MD) simulations of restricted primitive models of electrolytes where cations and anions are assumed to have the same size. The distribution of ions in the confinement created by unpolarizable or polarizable interfaces was extracted for a fixed interfacial separation of 3 nm for different values of electrolyte concentration (<1M), ion valency, and dielectric mismatch at the interface. Ions were found to deplete from the interfaces, adsorb to the interfaces, and show density oscillations. We showed that the forces that govern these structural features are not always directly exerted by the interfaces, which were charge neutral in some cases, but arise from the thermal motion of the ions and the geometric constraints that the interfaces impose.

In this paper, we use MD simulations to perform a systematic study of the ionic structure of aqueous monovalent electrolyte solutions confined by two planar interfaces over a wide range of concentrations c ∈ (0.1, 2.5)M, interfacial separations h ∈ (5, 8) nm, surface charge densities σ_s ∈ (−0.005, −0.02) C/m², and counterion sizes d₊ ∈ (0.2–0.63) nm. Our focus is on understanding the behavior of ions in the confinement created by charged interfaces under high electrolyte concentrations.^62,64 To keep the computational costs tractable for investigations over a wide range of interfacial separations and solution conditions, the primitive model for the electrolyte system is adopted. The model system of cations and anions with diameters corresponding to hydrated diameters of Na and Cl ions is comprehensively studied. To complement the MD simulation results, we also develop a variational approach based on a free energy functional that captures the key features of the electrolyte system. We note that other energy-minimization based approaches have been developed recently, in particular, in Ref. 72, which explore the transition in the screening behavior as a function of increasing surface charge density. Our variational approach is specifically tailored to explore the transition as a function of the concentration.

In order to probe the non-monotonic behavior of the decay length observed for a diverse set of electrolytes and material surfaces,^62–64,66 the effects arising due to specific material surfaces, including the effects of surface polarization charges, are neglected in favor of examining the general interplay of electrostatics-driven ion accumulation (depletion) near charged interfaces and the effects of steric correlations arising due to finite ion size. Our previous paper showed that the thermal forces arising due to the steric correlations can often overwhelm the effects due to the surface polarization charges for monovalent electrolytes at high concentration (≳0.1M).¹⁸ 

The ionic structure is quantified by evaluating the ion number densities, net charge densities, integrated charges,^72,75,77,78 and characteristic decay lengths associated with the screening of the charged interface. The results show the presence of two distinct regimes of screening behavior as the concentration c is changed from 0.1M to 2.5M for a wide range of electrolyte systems generated by tuning the interfacial separation, surface charge density, and ion size. In the low concentration regime, e.g., c ≲ 1M for electrolytes with cations and anions of sizes corresponding to hydrated sodium and chloride ions, the integrated charge exhibits a smooth, monotonic decay to 0 as the distance from the interface is increased. In this regime, the associated decay length decreases sharply with increasing c. In stark contrast, for high concentrations (c ≳ 1M), the integrated charge exhibits a non-monotonic, oscillatory decay to 0 as the distance from the interface is increased, with a decay length that rises with increasing c. These changes in the screening behavior are attributed to the rise in the strength of the steric ion–ion correlations with increasing concentration, which produce dramatic changes in the ionic structure including the enhanced accumulation of counterions and co-ions near the interface and the non-monotonic behavior of the net charge density.

This paper is structured as follows: In Sec. II, we provide the details of the simulation models and methods. In Secs. III A and III B, we present the simulation results for the ionic density profiles and the integrated charge. Simulation results for the effects of changing the surface charge density and the ion size on the ionic structure and the screening behavior are given in Secs. III C and III D. In Sec. III E, we provide the details of the phenomenological model and the associated results of the variational approach. Finally, Sec. IV presents discussion and conclusions.

Our model system consists of an electrolyte solution confined within two charged, planar interfaces parallel to each other. Each interface is characterized with a uniform surface charge density σ_s < 0. The simulation cell is a rectangular box of dimensions l_x × l_y × h, with l_x = l_y. Periodic boundary conditions are employed in the x- and y-directions. l_x and l_y are chosen to be sufficiently large to avoid any artifacts due to the periodic boundary conditions. z = −h/2 and z = h/2 planes are chosen as the location of the interfaces, and h is defined as the interfacial separation. The coordinate system is defined such that z = 0 corresponds to the midpoint between the interfaces. The distribution of ions is examined for 5 nm ≤ h ≤ 8 nm. The uniform charge density on the planar interfaces is simulated by meshing each interface with a large number of points M and assigning each point the same charge $q=σslx2/M$⁠. For all systems, M = 2500 is used. We verified that similar results for the ionic distributions were obtained with larger values of M. We also note that more efficient simulations can be performed using a recently developed method that bypasses the need to explicitly include charges on the interfaces by modifying the Ewald summation technique.⁷⁹ 

Ions are modeled using the primitive model of electrolytes^16,18,37 with hydrated diameters. For example, a model electrolyte having cations and anions of sizes informed by the hydrated size of Na⁺ (0.474 nm) and Cl⁻ (0.627 nm)⁸⁰ is extensively studied. Water confined within the planar interfaces is modeled as an implicit solvent with a relative dielectric permittivity of 80.

The Hamiltonian of the confined electrolyte system is the sum of the total steric potential energy between ions, the total electrostatic energy between ions, and the total energy associated with the interaction between the ions and the two planar interfaces. Ion–ion steric interactions are modeled using the standard purely repulsive and shifted Lennard-Jones (LJ) potential.^18,35 For a pair of ions separated by distance r, this potential is given as

for r ≤ 2^1/6d, where d is the average diameter of the two ions, k_B is the Boltzmann constant, and T is the temperature. For r > 2^1/6d, U_LJ = 0. Steric interactions between an ion and the interface are modeled using a similar repulsive LJ potential with a different size parameter d that represents the distance of the closest approach of an ion to the interface. Electrostatic interactions between ions and between an ion and the interface are modeled using the Coulomb potential. The long-range of the Coulomb potential is properly treated using Ewald sums⁸¹ or the method of charged sheets.^18,82,83 Either approach gives similar results within statistical uncertainties.

Simulations are performed using LAMMPS⁸⁴ as well as codes developed in our laboratory.⁸⁵ For all systems, the dimensions of the simulation box in the unconfined x and y directions are taken to be l_x = l_y = 18 nm. Depending on the electrolyte concentration, interfacial separation and surface charge density, the number of particles in the main simulation cell varied between 96 and 7884. Note that in addition to electrolyte ions, counterions are included in the confinement to ensure electroneutrality. Counterions are modeled as positively charged ions of the same diameter and charge as electrolyte cations. All simulations are performed in an NVT ensemble at a temperature of T = 298 K that is maintained using a Nose–Hoover thermostat.⁸⁶ Each system is simulated for 1 ns to reach equilibrium with a time step of 1 fs. After equilibration, systems are simulated for another 9 ns and trajectory data for computing ionic distributions are collected every 0.1 ps.

Table I shows the key model parameters and the range over which they are varied in order to probe the ionic structure associated with confined monovalent electrolytes. The electrolyte concentration is defined as c = N/V, where N is the number of electrolyte cations or anions and V = l_x × l_y × h is the volume of the simulation box. In most investigations, the anions (co-ions) are modeled as particles of 0.627 nm diameter (associated with the hydrated size of a Cl⁻ ion).

The number density in the direction perpendicular to the interfaces for cations and anions is extracted using the trajectory data collected during the simulation. The net charge density ρ(z) at z is calculated as

where e is the electronic charge and n₊(z) and n₋(z) are the number density profiles of positively charged and negatively charged ions, respectively. To characterize the effects of ion depletion or accumulation near the interface on the screening of the surface charge, the integrated charge⁷⁷ S is extracted using the following equation:^42,75

Here, σ_s is the surface charge density and ρ is the charge density given by Eq. (2). Note that S(z) is defined using the left planar interface as the reference surface. For clarity, we plot S vs z + h/2 where the latter denotes the distance from the left interface. S(z) is well defined for −h/2 ≤ z ≤ 0, i.e., up to the center of the confined region. The first term σ_s in Eq. (3) is added to shift the integrated charge in order to satisfy an appropriate boundary condition at the interface: S(−h/2) = σ_s. This condition implies that as the distance from the interface approaches 0, the surface is not screened and S is simply the bare surface charge density σ_s.

S(z) provides a way to measure the screening of the surface charge on the interface by the presence of electrolyte ions at a given concentration. The distance z^* from the left interface beyond which the integrated charge S remains confined to values within a small neighborhood of 0 can be interpreted as a length-scale associated with the decay of the electric field originating from the charged surface. In the approach to this region, S can exhibit either a smooth, monotonic behavior or an oscillatory, non-monotonic behavior with gradually decreasing amplitude. In practice, z^* is extracted as the minimum distance from the left interface such that for all distances greater than z^*, |S| is confined to values within 5% of its value at the interface (|σ_s|). In the case of the oscillatory integrated charge profile, S may become 0 multiple times, but our definition ensures that the decay length corresponds to the distance beyond which the amplitude of the oscillations has decreased to an extent that the largest fluctuations in S do not exceed a small fraction of its value at the interface.

It is useful to note that under conditions where the screening behavior is well described by the Debye–Hückel theory, S exhibits a simple exponential decay and the distance at which S decays to less than 5% of its initial value is ≈3λ_D, where λ_D is the Debye length. Therefore, for low concentrations, our chosen characteristic length z^* is identified up to a factor with the standard Debye length. However, as our results show, z^* diverges substantially from the Debye length at high concentrations, allowing us to identify the onset of a new regime.

The concept of integrated charge has been used in many studies to quantify the effects of the ionic structure on the charged surface.^{42,72,73,75,77,78,87} Instead of Eq. (3), S is often defined by normalizing the integral (second term on the right-hand side of the equation) by σ_s such that S decays to 1 instead of 0.^72,77 We also note that in a recent study using classical density functional theory, the decay length was extracted by examining the behavior of ions near a charged surface and calculating the slope of log |ρ|.³⁴ 

We first present the results for the ionic distributions extracted using simulations of the model electrolyte system comprising cations and anions of 0.474 nm and 0.627 nm diameters, respectively. These ion sizes are informed by the hydrated diameters of Na and Cl ions in water.⁸⁰ Figure 1 presents the density profiles n₊(z) and n₋(z) of cations and anions confined between planar interfaces separated by 5 nm and characterized with a surface charge density σ_s = −0.01 C/m². Profiles are shown for different electrolyte concentrations in the range c ∈ 0.1M–2.5M. For each system, n₊(z) and n₋(z) are symmetric around z = 0, as expected. For c = 0.1M, cations (counterions) are accumulated near the planar interfaces, while anions (co-ions) are depleted near the interfaces. This behavior can be attributed to the electrostatic force exerted by the negatively charged surface. As c ≳ 0.5M, both cations and anions are accumulated near the interfaces as evidenced by the peaks of n₊(z) and n₋(z) near the interfaces. The accumulation of anions near the interfaces can be attributed to the thermal forces arising from the stronger steric correlations between anions that overcome the repulsive electrostatic force between the anions and the charged surface. The location of the peak of n₋(z) is farther from the interface compared to that of n₊(z) because of the larger size of the anions.

The accumulation of cations and anions near the interfaces continues to rise with increasing c. Figures 1(e) and 1(f) show that the peak of n₋ near the interface becomes higher than that of n₊ for c ≳ 2M. This effect can be attributed to the larger size of anions compared to cations and would disappear in size-symmetric electrolytes (see Fig. 16 in the  Appendix). The near-interface accumulation of ions of both species leads to the formation of structured layers that act as soft walls and induce secondary accumulation, thus nucleating more layers. The cascade of soft walls resulting from the external planar interface creates oscillations in the ionic structure for c = 2M and 2.5M. The changes in the distribution of confined cations and anions with increasing c for an interfacial separation of h = 5 nm are also observed for different h = 6 nm, 7 nm, 8 nm (Fig. 13 in the  Appendix shows the results for h = 7 nm).

We next examine the evolution in the ionic structure using the net charge density ρ(z) extracted from the associated number densities n₊(z) and n₋(z) by employing Eq. (2). Figure 2 shows ρ(z) for the same systems explored in Fig. 1. ρ is symmetric around the center of the confinement (z = 0). For the ease of exposition, we discuss the behavior of ρ near the left interface. In the immediate vicinity of the interface, ρ is positive and exhibits a peak close to the interface that can be attributed to the accumulation of cations (counterions). The height of this peak increases with increasing c from 0.1M to 2.5M. For c ≲ 0.5M, ρ exhibits a monotonic decay to 0 as z approaches 0 (inset of Fig. 2). For c ∼ 1M, a region close to the interface develops where ρ becomes negative, reaching a minimum of ≈−0.1eM before decaying to 0. This non-monotonic behavior can be attributed to the gradual increase in the accumulation of anions (co-ions) near the interface with increasing concentration.

As c rises to values ≳2M, a prominent ρ < 0 region develops adjacent to the first positive peak close to the interface and the minimum value of ρ drops sharply. The drop can be attributed to the strong accumulation of anions near the interface shown in Figs. 1(e) and 1(f). The excess net negative charge associated with the layer of anions induces a second positive peak adjacent to it, which, in turn, nucleates a second ρ < 0 region farther from the interface (seen more clearly for c = 2.5M in Fig. 2). The resulting oscillatory, non-monotonic behavior of ρ is associated with the structural modulations present in n₊ and n₋ shown in Fig. 1. The amplitude of the oscillations is reduced, and ρ gradually decays to 0 as z approaches the center of the confinement (z → 0).

It is useful to analyze the charge density ρ for different interfacial separations at a given electrolyte concentration. Figure 3 shows ρ as a function of the distance z + h/2 from the interface for h = 5 nm, 6 nm, 7 nm, 8 nm for c = 0.5M (a) and 2.0M (b). By symmetry considerations, the distance is measured relative to the left interface. For both low and high electrolyte concentrations, the data for ρ under different h = 5 nm, 6 nm, 7 nm, 8 nm fall on the same universal curve. For c = 0.5M, ρ exhibits an initial peak near the interface and then decays monotonically to 0 as z is increased. For higher c = 2M, ρ exhibits a non-monotonic, oscillatory behavior (with a prominent ρ < 0 region) before decaying to 0 for larger z + h/2. The distinct behavior of ρ under low and high c conditions is similar to the trend observed in Fig. 2.

We next analyze the integrated charge S(z) as a function of the distance z + h/2 from the left planar interface for different electrolyte concentrations. S is extracted using Eq. (3) by employing the results for the net charge density ρ(z) shown in Fig. 2 for an interfacial separation of h = 5 nm. Figure 4 shows S(z) vs z + h/2 for the same systems studied in Figs. 1 and 2. For all concentrations c ∈ (0.1, 2.5)M, S = σ_s = −0.01 C/m² when z + h/2 ≲ d₊/2. In other words, the surface charge on the interface is unscreened up to a distance of the closest approach of the cation. For all c, S(z) exhibits an initial rise as the distance z + h/2 is increased, implying the reduction in the integrated charge due to the screening of the surface by the electrolyte ions. For c ≲ 0.5M, S(z) exhibits a monotonic rise before decaying to 0 as the distance z + h/2 from the interface becomes large. Under these conditions, S ≤ 0 for all z, showing no evidence for any charge inversion. For c = 1M, S increases sharply with z, reaching a value of 0 at a distance z + h/2 ≈ 0.3 nm away from the interface. S continues to rise further, peaking at a value of ≈0.003 C/m² (≈0.3|σ_s|) before decaying to 0. The presence of the S(z) ≳ 0 region for z + h/2 ≳ 0.3 nm suggests a weak charge inversion.

For c ≳ 1.5M, S continues to rise rapidly to values beyond 0, reaching much higher peaks (e.g., >0.03 C/m² for 2.5M), signaling greater charge inversion. More importantly, instead of decaying smoothly down to 0 (like for c = 1M), S exhibits a non-monotonic, oscillatory behavior and becomes negative with increasing z + h/2, exhibiting a minimum before decaying to 0. The oscillations are stronger for c = 2.5M, where a clear secondary peak ≈1 nm away from the interface is observed. Similar results for S(z) are obtained for electrolytes confined under different interfacial separations h = 6 nm, 7 nm, 8 nm. This finding is expected, given the variation in ρ as a function of z + h/2 for different h (Fig. 3). Figure 14 in the  Appendix shows S for h = 7 nm.

Using S(z), we extract the decay length z^* following the procedure described in Sec. II A. Recall that z^* is the minimum distance from the left interface such that for all distances greater than z^*, |S| is confined to values within 5% of its value at the left interface. z^* acts as a measure of the extent of the screening of the charged surface by the electrolyte ions. The availability of data for S(z) vs z + h/2 for many h values enables an accurate assessment for z^*, providing also the associated statistical uncertainties. z^* is computed by taking the average of the decay lengths obtained for different h = 5 nm, 6 nm, 7 nm, 8 nm. The inset of Fig. 4 shows z^* vs c. The error bars are computed by extracting the standard deviation associated with the z^* data obtained for different h = 5 nm, 6 nm, 7 nm, 8 nm. We find that z^* decreases sharply as c is increased up to ≈1M. However, for c ≳ 1M, z^* increases slightly with increasing c. We note that the variation in the decay length vs c for c ≳ 1M exhibits greater noise, which can be attributed to the extraction process described in Sec. II A that examines the region where |S| is confined to values within 5% of its value at the interface (−0.0005 C/m² < S < 0.0005 C/m²). In some cases, at high concentrations, the oscillatory nature of S produces variations in amplitude close to these small threshold values defining the region, inducing greater sensitivity and uncertainty in the determination of z^*.

The results shown thus far have been obtained for planar interfaces characterized with a surface charge density σ_s = −0.01 C/m². It is useful to assess how the distinct ionic structure observed under low and high electrolyte concentrations is affected by changing σ_s. We performed simulations of the same model electrolyte system confined by interfaces for two different values of σ_s = −0.005 C/m², −0.02 C/m². Changing σ_s alters the distributions of ions in the confinement; however, the distinct features associated with the ionic structure separating the two regimes of low c (≲1M) and high c (>1M) persist for different σ_s values. Figure 5 shows the distribution of cations and anions confined by surfaces characterized with σ_s = −0.02 C/m² for c ∈ (0.1, 2.5)M. We find that the larger negative charge on the interfaces leads to enhanced accumulation of positive ions and stronger depletion of anions compared to the σ_s = −0.01 C/m² case. For example, unlike the behavior observed in Fig. 1(b), anions do not exhibit a clear accumulation near the interfaces for c = 0.5M. Furthermore, even at the highest c = 2.5M, the peak associated with the anions does not significantly exceed the cation peak, as was the case for σ_s = −0.01 C/m² for c ≳ 2M [Figs. 1(e) and 1(f)].

Figure 6 shows the integrated charge S(z) as a function of the distance z + h/2 from the left interface characterized with σ_s = −0.02 C/m² for different c ∈ (0.1, 2.5)M. Quantitative differences are observed in S(z) compared to the results for the system with σ_s = −0.01 C/m² (Fig. 4). However, the overall S behavior at low c is qualitatively distinct from the behavior at high c, consistent with the smaller σ_s results. When the concentration is low (c ≲ 1M), S(z) ≤ 0 and exhibits a monotonous rise before decaying to 0 for large z + h/2. The screening behavior for high c ≳ 1.5M is entirely different. Similar to S(z) for σ_s = −0.01 C/m² (Fig. 4), S rises rapidly to values >0 and then exhibits a non-monotonic, oscillatory decay to 0 for large z + h/2. Similar results are observed for the case of σ_s = −0.005 C/m². Figure 6 (inset) shows z^* as a function of electrolyte concentration c for different σ_s values. We find that z^* vs c for σ_s = −0.005 C/m², −0.02 C/m² exhibits a similar behavior as the result for σ_s = −0.01 C/m². z^* decreases as c rises up to ∼1M. For c > 1.0M, z^* rises with increasing c.

We now investigate the effects of changing the ion size on the distributions of ions confined by planar interfaces for different electrolyte concentrations. Our focus is on assessing how the changes in the steric correlations influence the observed distinct regimes of screening behavior under low and high c conditions. The diameter of the anion is fixed at d₋ = 0.627 nm. We perform simulations of different electrolyte systems generated by changing cation diameter to d₊ = 0.209 nm, 0.313 nm, and 0.627 nm. Other parameters are the same as in the model electrolyte system studied in Sec. III A. Simulations are performed for interfaces characterized with a surface charge density σ_s = −0.01 C/m² and for different interfacial separations h = 5 nm, 6 nm, 7 nm, 8 nm. We also simulate a size-symmetric electrolyte with smaller-sized ions of diameters d₊ = d₋ = 0.313 nm as a reference system to assess the contributions of stronger steric effects.

Figure 7 shows the density profiles of cations and anions near the left interface for systems with cations of diameters d₊ = 0.209 nm, 0.313 nm, 0.474 nm, and 0.627 nm under low c = 0.1M (top row) and high c = 2M conditions (bottom row). For c = 0.1M, as d₊ increases, cations accumulate farther from the interface due to the increase in the excluded volume [Fig. 7(a)]. In addition, their peak density increases slightly with increasing d₊. The distribution of anions is weakly affected by changes in the cation size [Fig. 7(b)]. In all cases, anions deplete from the interfaces, and the depletion is slightly stronger for larger d₊. This can be attributed to a more effective screening of the surface charge by smaller-sized counterions as they can accumulate closer to the interface, which decreases the electrostatic force with which the anions are repelled by the interface.

Changing cation size has more dramatic effects for c = 2M. Here, again due to the increase in the excluded volume, cations accumulate farther from the interface as d₊ increases [Fig. 7(c)]. However, the peak density of cations rises much more rapidly with cation size (from ≈2.8M at d₊ = 0.209 nm to ≈7.5M at d₊ = 0.627 nm). This can be attributed to an overall increase in the ion–ion steric correlations, which push the ions toward the interface (similar to the case where the ion density near the interface rises because of increasing c). More importantly, the distribution of anions is strongly affected by the changes in cation size. Contrary to the trend observed for c = 0.1M, the peak anion density increases significantly with increasing d₊ at c = 2M [Fig. 7(d)]. Furthermore, as d₊ increases, the density of anions exhibits more modulations. In Fig. 7(d), all anions have the same size and all electrolytes are at the same concentration (2M). The dramatic changes in the ionic structure seen here arise due to the strong steric correlations between cations and anions under high c conditions.

It is instructive to examine the changes in the charge density ρ vs concentration c for electrolytes with different cation sizes. Figures 8(a) and 8(b) show ρ near the left interface for electrolytes with cations of diameters d₊ = 0.209 nm and d₊ = 0.627 nm, respectively. The results are shown for electrolytes confined within interfaces separated by h = 5 nm. ρ is derived using the cation and anion distributions associated with these systems (Figs. 15 and 16 in the  Appendix) for c ∈ (0.1, 2.5)M. The overall behavior of ρ vs c is similar to that observed in Fig. 2 for the electrolyte with cations of size d₊ = 0.474 nm. ρ exhibits a transition from a monotonic behavior for sufficiently low c to a non-monotonic, oscillatory behavior for high c. For the electrolyte system with smaller-sized cations, this transition occurs at a lower c (∼0.5M) compared to the case of electrolytes with cations of size d₊ = 0.474 nm. The modulations in the ionic structure at high c for the electrolyte system with cations of size d₊ = 0.627 nm are found to be relatively milder, which can be attributed to the size-symmetric nature of this electrolyte. Similar evolution in ρ with changes in c is observed for systems confined under different interfacial separations h = 6 nm, 7 nm, 8 nm.

We now examine how the variations in the ionic distributions resulting from changing d₊ affect the screening of the surface charge by the electrolyte solution. Figure 9 shows the integrated charge S(z) vs the distance z + h/2 from the left interface for electrolytes confined by the interfaces at c = 2M and 0.1M (inset). The interfacial separation is h = 7 nm. Open symbols correspond to results for electrolytes with anions of fixed size d₋ = 0.627 nm and cations of different sizes d₊ = 0.209 nm, 0.313 nm, 0.474 nm, and 0.627 nm. Closed symbols are the results for the size-symmetric electrolyte system with ions of diameter d₊ = d₋ = 0.313 nm.

Figure 9 (inset) shows that when the concentration is low (c = 0.1M), S ≤ 0 and exhibits a monotonous rise toward saturation to 0 as z + h/2 increases for all systems. The behavior at high c = 2M is entirely different (except for the case of the size-symmetric electrolyte with ions of 0.313 nm diameter; see below). S(z) is non-monotonous and exhibits charge inversion for all d₊. For d₊ = 0.209 nm, 0.313 nm, and 0.474 nm, S rises rapidly to values >0, exhibits a peak, and then decays to 0 in an oscillatory fashion. The maximum peak of S decreases as d₊ increases from 0.209 nm to 0.474 nm; however, its location remains the same (at ≈0.4 nm). For d₊ = 0.627 nm, S rises much more slowly to values >0 and exhibits a smaller and broader peak, before decaying to 0 with relatively milder oscillations. This can be attributed to the size-symmetric nature of this electrolyte (d₊ = d₋ = 0.627 nm). Overall, these results indicate that the distinct nature of the ionic structure and the screening behavior for low and high c conditions are present over a wide range of cation sizes for electrolytes with anions of 0.627 nm size.

For the size-symmetric electrolyte with smaller-sized ions of diameter d₊ = d₋ = 0.313 nm, S(z) behavior under low and high c conditions is similar. We find that even at c = 2M, S exhibits a monotonic rise with increasing z + h/2 before decaying to 0 for large z + h/2. For this system, there is no evidence of two distinct regimes of different ionic structures and screening behavior within c ∈ (0.1, 2.5)M. We expect that the non-monotonic behavior in S for this particular system appears at large c > 2.5M where steric correlations become strong enough to produce structured layers of ions near the interfaces.

Figure 10 shows the decay length z^* as a function of the concentration c for the different systems shown in Fig. 9. z^* is computed as an average of the results of the decay lengths extracted using S(z) data for different interfacial separations h ∈ (5, 8) nm. We ensure that data from sufficiently large h are employed to enable a meaningful extraction of the decay length, in particular, when the latter is large. z^* decreases sharply for all cases as c is increased up to 0.5M. For electrolytes with anions of diameter d₋ = 0.627 nm and cations of diameter d₊ ∈ (0.209, 0.627) nm, z^* vs c behavior is non-monotonic. For these systems, as c ≳ 1M, z^* does not decrease with increasing c but rises slightly. The concentration where the behavior switches from sharp fall to mild rise exhibits a dependence on the ion size. For example, the crossover c for the electrolyte with d₊ = 0.627 nm is ≈0.5M compared to ≈1M for the system with cations of 0.474 nm size. In striking contrast, for the size-symmetric electrolyte with smaller ions of 0.313 nm diameter, z^* exhibits a monotonic decrease with increasing c for c ∈ (0.1, 2.5)M. We hypothesize that for this system, the crossover concentration is >2.5M.

Inspired by other studies,^62,74 an attempt is made to collapse the data from different systems by scaling z^* with λ_D and plotting it against the scaled concentration represented as d₋/λ_D, where $λD=1/8πlBc$ is the Debye length. The inset of Fig. 10 shows the log–log plot of z^*/λ_D vs d₋/λ_D. Two distinct regimes are observed. For d₋/λ_D ≲ 2 (small c), the scaled decay length is roughly constant with a value of ≈3. In other words, the screening of the charged surface follows the behavior predicted by the Debye–Hückel theory under these conditions. Recall that the decay length, as noted in Sec. II A, is ≈3λ_D for conditions where the Debye–Hückel description is applicable. For d₋/λ_D ≳ 2 (large c), z^*/λ_D rises with increasing c. We find that this rise exhibits a power-law behavior $z*/λD=(d−/λD)n$ with exponent n ∼ 1.5, although we realize that this scaling is observed over a comparatively limited range of data.^62,74

Simulation results point to the emergence of structured layers of ions near the interface and associated distinct screening behavior when the electrolyte concentration c is increased. Figures 11(a) and 11(b) show representative simulation snapshots of ions near the left interface for electrolyte concentrations of 0.5M and 2.5M, respectively. The images shown correspond to the electrolyte system of cations (red) and anions (blue) of sizes of 0.474 nm and 0.627 nm, respectively, confined by interfaces characterized with a surface charge density of σ_s = −0.02 C/m². The leftmost images in (a) and (b) show ions whose centers are within a distance of ≈0.6 nm from the left interface. The next three images from left to right represent the features in the different layers extracted by taking thin slices of the leftmost images. Slice I shows ions within ≈0.3 nm distance from the interface. Slice II shows ions within a thin layer of ≈0.1 nm width, bordering slice I to its left. Slice III shows ions within a layer of ≈0.2 nm width, bordering slice II to its left.

Slice I in either system is dominated by counterions; the number of counterions is greater for the high c system (b) compared to the low c system (a). Slice II associated with the high c system exhibits a much larger population of co-ions (anions) compared to cations in contrast with the low c case where a similar number of cations and anions are observed. Images associated with slice III show that farther from the interface, regions have a mixed population of cations and anions, with the high c system exhibiting a more packed arrangement of ions compared to the sparse layer of ions for the low c system. These snapshots are consistent with the quantitative plots shown above capturing the differences in the ionic structure under low and high c conditions.

It is possible to gain insight into the simulation results using a simple mean free energy model where key features of the electrolyte system emerge naturally. A variational scheme is employed, where two distinct types of behavior are considered. The first type admits a Poisson–Boltzmann description based on a regular solution model entropy for the ions in solution. A second option considers that ions near the charged interfaces appear in highly structured distributions, forming soft aggregates [similar to the visualizations shown in Fig. 11(b)].

We consider a semi-infinite electrolyte with a single charged planar interface. Near the interface, the system can acquire a spatially varying net charge density ρ. This charge density, along with external charges, produces a mean field potential ϕ. The contribution to the energy of the system is then the integral of ρϕ/2. Schematically, the free energy can be written as an integral over the sum of the mean field electrostatic energy and a functional density f_ℓ describing the local behavior of the ions,

The functional f_ℓ has different forms according to the location considered. Taking the interface to be at z = 0, we define a near-interface region, 0 < z ≤ Ld, extending L ion diameters into the bulk. For simplicity, we have considered cations and anions of the same diameter d. In this region, the functional is f_ℓ = f_s, where f_s is constructed using a description of the ions assembled into soft aggregates. Beyond this region, where z > Ld, we use f_ℓ = f_b. This second functional f_b assumes a regular solution model for the entropy of the ions.

In the variational scheme we consider, the near-interface region is not determined from the outset and can have zero thickness, in which case the whole system is described by the bulk properties. The thickness of this layer is an outcome of the analysis of the model and depends on the bulk concentration and other parameters. Both explicit functionals f_s and f_b depend on the mean number density n = (n₊ + n₋)/2 and the charge density ρ. These quantities vary only along the z direction. In the bulk, the mean density takes the value n₀. In both bulk and near-interface regions, the fields that appear in the free energy are required to satisfy several conditions. The potential should decay to zero in the bulk. The potential must satisfy the Poisson equation for a source equal to the mean-field charge. The electric field at the interface is determined by the external charges. Additionally, the near-interface and bulk regions have equal chemical potentials for the species. All these requirements can be implemented explicitly in the variational functional by means of Lagrange multipliers. For expediency in this presentation, we omit these details.

In the bulk region, the free energy density has the regular solution model form,

where j runs over the ion species and solvent and v is a thermal volume for the solution. The functional can be evaluated and expanded around a uniform bulk state and can be shown to lead to a Poisson–Boltzmann description of the bulk response to external fields.

Next, we construct the functional for the near-interface region. Here, ions are assumed to be found in aggregates of size N_c that are considered to act as fundamental particles. The entropy of the aggregates within the solvent is described by a regular solution model. In addition, within the aggregates, the ions interact with their neighbors and produce an effective cohesive energy per particle, E_c. This energy arises mainly from electrostatic interactions. We note, however, that the presence of the interface imposes important steric constraints that facilitate the organization of the ionic aggregates. Thus, this average free energy per particle still contains information about the entropy of the aggregate. For simplicity, all ions within the near-interface region are assumed to be present in these aggregates. In the near-interface region, the model explores the possibility of a number density n_w, with a step function profile uniform through the region, that might be different from the bulk value n₀. The associated free energy density is

The second term in Eq. (6) is the entropy of a mixture of only two species: the aggregates with number density n_c = 2n_w/N_c and the solvent.

It is important to emphasize that the cohesive energy E_c though mostly of electrostatic origin is not captured by the mean field contribution proportional to ρϕ. This is the case, for example, in ionic crystalline solids where, when averaged over large length scales, both charge density and field are zero. Instead, the cohesive energy is well described by a Madelung constant, which evaluates the sum of Coulomb interactions of the whole system with a representative member of the lattice. This idea has been used to model macroions collapsed with their counterions and to describe the adsorption of ions to charged walls.^59–61 The precise value of the cohesive energy is difficult to calculate, but it should be nearly independent of the bulk concentration. Thus, we use the assumption that E_c is a constant and set its value to obtain agreement with the observed simulation results. It is important to note, however, that we assume that this term is only present near the interface.

Once the functional is specified, standard variational arguments indicate that the best approximation to the properties of the system is obtained when a selection of ion number density and charge distribution minimizes the free energy. For this, it is necessary to explore all possible ion distributions. However, these configurations are limited by steric effects. To make calculations feasible, as well as to incorporate features of the observed behavior in simulations, we consider only step-wise charge distributions, with steps corresponding to positions (k − 1/2)d away from the interface, where k is an integer. That is, we assume that within the aggregates, there are well-defined layers of ions. This is consistent with simulation results (Figs. 1, 2, 4, and 11) showing enhanced occupation of these near-interface regions by ions. The charge density ρ can only take values from −2en_w to +2en_w and can be expressed as a fraction of the mean ion number density ρ = p_k2en_w, where −1 ≤ p_k ≤ 1 and e is the elementary charge. p_k refers to the charge fraction associated with the layer number k. Thus, the possible charge distributions considered are

For a given bulk particle density, and a prescribed surface charge density, the free energy functional can now be minimized with respect to the charge fraction p_k in the layers, the number density n_w near the interface, and the total number of structured layers L. This minimization is carried out under the constraints noted above. Finally, while the number of layers L is variable in principle, we only consider cases where these layers exhibit non-negligible charge density. We do not consider possibilities that extend the region of aggregates indefinitely into the bulk. This is in accordance with the observation that the external planar interface has an ordering effect that is not likely to extend into the bulk. As observed in our simulations, the net charge density oscillates but rapidly decreases to zero. An ordering transition in the bulk of the electrolyte is expected to appear at high ion concentrations. The order observed at the interface appears, however, at lower concentrations aided by the steric effects of the interface and the presence of the electric field due to the surface charge.

Analysis of the model produces the following results. For a given value of the surface charge density, as the bulk concentration increases, a transition is observed from a fluctuation-dominated state without a distinct near-interface behavior to a state exhibiting structured layers of finite thickness near the interface. The crossover concentration depends on the surface charge density but is not highly sensitive to this value, indicating the important role of the layering favored by the steric effects due to the interface. A concrete example is shown in Fig. 12 where the integrated charge and the associated charge densities are provided. The parameters used roughly correspond to simulation conditions for size-symmetric electrolyte ions with diameter d equal to the Bjerrum length l_B in water and σ_s = −0.01 C/m². The example assumes an aggregate size of N_c = 30, but the results are insensitive to this selection and are nearly identical for values above this size. Similarly, different values of the surface charge density, comparable to those used in simulations, lead to qualitatively similar results.

The onset of aggregate formation corresponds to concentrations of 1M. To observe this transition, the required cohesive energy is E_c = −0.8k_BT. We find that the transition is sharp and a near-interface region of structured layers of thickness L = 3 is preferred above the concentration threshold (crossover concentration). In Fig. 12, this region extends up to z ≈ 2.1 nm. Within the near-interface region, the ionic density is higher than in the bulk, in contrast with the Poisson–Boltzmann case where it remains the same, and increases slowly with bulk concentration. Above the transition, the charge density is effectively independent of the bulk concentration. This density, in all cases, is oscillating. The first layer reverses the charge of the interface, the second layer of charge again inverts the net charge, and so does the third. This oscillatory nature of the charge density makes the integrated charge non-monotonic. As in this regime the charge density is effectively independent of the bulk concentration, the screening length becomes a constant, z^* = Ld, i.e., three ion diameters in the model considered. We emphasize that this behavior is in sharp contrast with the low concentration regime where the integrated charge has a monotonic behavior, and the resulting screening length is shortened with increasing concentration.

Using molecular dynamics simulations of the primitive model of electrolytes, we performed a systematic study of the ionic structure of aqueous monovalent electrolyte solutions confined by two planar interfaces over a wide range of electrolyte concentrations c ∈ (0.1, 2.5)M, interfacial separations h ∈ (5, 8) nm, surface charge densities σ_s ∈ (−0.005, −0.02) C/m², and counterion sizes d₊ ∈ (0.2–0.63) nm. Our focus was on understanding the behavior of ions in the nanoconfinement created by interfaces under high electrolyte concentrations and how the ionic structure influences the screening of the charged interfaces. The ionic structure was quantified by evaluating the density profiles of ions, net charge densities, integrated charges, and decay lengths associated with the screening of the charged interface.

The results show the presence of two distinct regimes of screening behavior as the concentration is changed from 0.1M to 2.5M for electrolytes with cations and anions of sizes corresponding to hydrated sodium and chloride ions. For low c ≲ 1M, the integrated charge exhibits a monotonic decay to 0 with a decay length that decreases sharply with increasing c. On the other hand, for high c ≳ 1M, the integrated charge has a non-monotonic, oscillatory behavior, signaling charge inversion and formation of structured layers near the interfaces. The decay length under these conditions rises with increasing c.

The changes in the screening behavior for electrolytes are observed over a wide range of systems generated by tuning the interfacial separation, surface charge density, and size of the ions. The distinct regimes of the screening behavior are attributed to the dramatic changes in the ionic structure with an increase in c including the enhanced accumulation of both counterions and co-ions near the interface and the non-monotonic behavior of the net charge density. Both these changes are driven by the rise in the strength of the steric ion–ion correlations. When the ion size is reduced to d₊ = d₋ = 0.313 nm, the screening behavior for electrolytes within the concentration range c ∈ (0.1, 2.5)M does not exhibit any significant deviations from the Debye–Hückel scaling of the decay length as a function of concentration.

Our results directly probe the effect of increasing the electrolyte concentration on the screening of the charged surfaces. Recent studies have examined these effects by extracting the correlation length associated with the charge–charge pair correlation functions in bulk electrolyte solutions.⁷⁴ These studies have found that the correlation length rises with concentration for sufficiently high c as a power law. The exact scaling relation between the correlation length and concentration remains unclear; different approaches including liquid-state theories, atomistic simulations, and experimental studies place the power-law exponent n to be between 1 and 3.^63,74 Our finding of n ≈ 1.5 derived by examining the decay of the integrated charge near a charged surface screened by electrolyte ions is within this range.

We realize that although our results are accurate for the primitive electrolyte model, the model system itself assumes a homogeneous structureless solvent. The results obtained with the primitive model of ions of fixed, hydrated sizes ignore the important effects of the partial or complete shedding of the ion-hydration layers near the interface. These effects can modify the ion depletion and accumulation behavior near the interface and alter the associated screening of the charged surface. Work on quantifying the effects of the solvent structure is under way, and initial results of atomistic simulations of confined aqueous NaCl solutions show a screening behavior qualitatively similar to that obtained with the implicit-solvent model including the transition of the integrated charge from monotonic decay to non-monotonic decay as c is increased. A comprehensive study based on simulations of explicit-solvent models is needed to distill the contributions of solvent effects vs ionic correlations toward the evolution of structural features and associated changes in the screening behavior.

We note that theoretical studies using both implicit-solvent and explicit-solvent descriptions have predicted the rise in the decay length with increasing c for high c conditions.^30,74 However, the two models can yield distinct quantitative scaling behavior, e.g., the power-law exponent associated with the rise in the decay length can be different for a two-component model (cations and anions) compared to a three-component model (cations, anions, and solvent particles).⁷⁴ 

Our model system assumed unpolarizable material surfaces. The dielectric permittivity of the solvent was also assumed to be the same for all electrolyte concentrations. In our earlier work, we found the effects of polarization charges to be generally weaker for monovalent electrolytes at concentrations ≳0.1M.¹⁸ Furthermore, we performed simulations of monovalent electrolyte ions in water (dielectric permittivity of ɛ = 80) confined between two uncharged, polarizable surfaces separated by 5 nm and characterized with material permittivity ϵ_m = 2 using methods outlined in previous papers.^18,55,88 These simulations showed that surface polarization charges change the ionic structure minimally for 0.5M and the effects are much more suppressed for 2.0M. The effects of polarization charges are expected to be further overwhelmed by the free surface charges in the presence of charged interfaces. Furthermore, we note that our results for the primitive model system can offer pathways to quantify those additional effects arising due to polarizable surfaces.

Finally, we note that the deviation in the scaling behavior of decay length vs concentration can occur even at lower concentrations when electrostatic coupling is higher (e.g., with multivalent ions or under conditions of low temperatures or solvents with low dielectric permittivity).³³ Our results do not describe these effects.

To complement the MD simulation results, we constructed a minimal phenomenological model for the system that explores a subset of possible ion number and charge densities. The model can be considered as a simplified version of the more detailed density functional approaches^34,72 that nevertheless highlights the key features of the system. Two features not included in the model are the smooth decay of the number density into bulk and the detailed structure of the species within the near-interface region. We have not explicitly modeled the steric effects of the interface and their decay into the bulk that would lead to smoother behavior of the number density. A more detailed study of such effects has been carried out in Ref. 72 where the effects of finite size of ions are incorporated via modifications to the Poisson–Boltzmann equation. This approach successfully recovers the presence of layered structures. In our work, this feature is directly implemented as one of the possible states of the system, as an alternative to a fluctuation-dominated solution state. The description of the distributions of ions of different species is also not directly considered, which limits our results to the size-symmetric electrolyte.

The model implements the observation that, near the charged interface, both simulations and the simplest interpretations of experimental results⁶² indicate the presence of strongly correlated structures. Beyond articulating this observation, the model ties together other aspects of the system behavior. The structured region near the surface is described by a free energy functional that scales differently than at the bulk. As a result, we obtain number density profiles higher near the interface. This is consistent with observations of increased ion density in Figs. 1(d)–1(f), 5(c), and 5(d), which are not found at lower concentrations and are not predicted by Poisson–Boltzmann descriptions of charge accumulation at charged surfaces. For high concentrations, when the structured state of the model is preferred, the system can still select among charge distributions that are consistent with a layered structure. As shown in Fig. 12, this density is always oscillating and within the model largely independent of the bulk concentration, leading to constant decay lengths in this regime. The simulations show smoothly changing charge densities that nevertheless preserve peak locations consistent with layering and decay lengths that increase slightly at high concentrations.

This research was supported by the National Science Foundation through Award Nos. 1720625 and DMR-1753182. Simulations were performed using the Big Red II and III supercomputing systems.

The data that support the findings of this study are available within the article or from the corresponding author upon reasonable request.

Figure 13 shows the ionic density profiles for the electrolyte system discussed in Sec. III A when the ions are confined between two negatively charged interfaces separated by h = 7 nm at 0.1M (a), 0.5M (b), 2.0M (c), and 2.5M (d). The results are similar to the density profiles shown in Fig. 1 for the case where the interfacial separation is 5 nm. At 0.1M, cations accumulate near the interfaces, while anions are depleted near the interfaces. However, for concentration c ≳ 2M, both cations and anions exhibit enhanced accumulation near the interfaces.

Figure 14 shows the result for integrated charge S(z) vs z + h/2 for the electrolyte system considered in Sec. III B confined between two charged interfaces with a separation of h = 7 nm. The results are similar to the integrated charge S for h = 5 nm (Fig. 4).

Figures 15 and 16 show the ionic density profiles corresponding to the electrolyte system for which the net charge density profiles are shown in Fig. 8 in the main text. Figure 15 shows the ionic distributions for the system with cations of diameter d₊ = 0.209 nm and anions of diameter d₋ = 0.627 nm. Figure 16 shows the results for the same system when cations of diameter d₊ = 0.627 nm are considered. The anion peak density for c ≳ 1.0M for the system with cations of 0.209 nm size is higher than the cation peak density near the interface. On the other hand, under similar conditions, the anion peak density near the interface is smaller compared to the cation peak density for the size-symmetric electrolyte system with cations of 0.627 nm size.