Describing screening in dense ionic fluids with a charge-frustrated Ising model

The thermodynamic properties of ionic fluids are governed by long-range Coulomb interactions between ions¹ in addition to the short-range molecular interactions present in neutral liquids. Strong electrostatic correlations lead to counter-intuitive phenomena in dense ionic fluids such as charge inversion,^2–8 altered capacitance at electrode–fluid interfaces,^6,9–11 and the recently observed “anomalous screening” in surface force experiments.^12,13 These effects could be important in the self-assembly of a variety of biomolecules¹⁴ and soft materials.¹⁵ Electrostatic correlations can be particularly pronounced in molten salts and ionic liquids which comprise ions alone and no neutral solvent molecules. The novel properties of such purely ionic fluids make them useful for a variety of scientific and technological applications, such as energy storage,^16,17 as industrial lubricants,¹⁸ and of serving as media capable of supporting stable colloidal nanoparticles.^19,20

Pure ionic fluids are ideal model systems for the theoretical study of the statistical physics of strongly correlated electrostatics without the complicating ion-specific effects of hydration in aqueous solution.²¹ A theoretical description of dense ionic fluids must go beyond the classic Debye–Hückel (DH) theory, which is valid only for dilute electrolytes with weak inter-ionic correlations,²² or equivalently, small inverse Debye screening length (also known as the Debye constant) in relation to the inverse molecular size $κD≡(4πρq2)/(ϵkBT)≪σ−1$⁠, where ρ is the concentration of ions (per unit volume), q is the unit charge, ϵ is the dielectric constant of the electrolyte, k_BT is the thermal energy, and σ is the ion diameter. Indeed, recent surface force experiments using concentrated solutions of salts and ionic liquids measure screening lengths, 1/κ_s, well in excess of the DH prediction, 1/κ_D, and show non-monotonic dependence of κ_s on κ_D.^12,13,23 Especially surprising is the universal scaling collapse of κ_sσ when plotted against κ_Dσ, despite the use of a range of ion types, solvent types, and ion concentrations.^13,23 The particular scaling behavior in the dense ionic regime, $κs∼κD−2$⁠, is not predicted by existing theoretical results, suggesting the need to go beyond standard approaches in the field.

A variety of theoretical techniques have been used to extend the DH theory to the strong Coulomb coupling or high κ_D regime.²⁴ To take two examples, Attard uses a standard closure from the theory of liquids,²⁵ while Lee and Fisher generalize the DH theory by considering an oscillatory potential that intuitively arises from the preference of oppositely charged ions to arrange in alternating layers.²⁶ Both of these theories result in a regime at large κ_D where spatial correlations between ions cannot be ignored as they are in the DH theory.¹ Indeed, the manifestation of these correlations as oscillations in the charge density was predicted long ago by Kirkwood.²⁷ In this large κ_D regime, the charge correlation length can become much longer than the screening length predicted by DH theory, qualitatively similar to observations of anomalous screening in the aforementioned surface force experiments. More recent work based on both simulations and phenomenological theories reproduces this oscillatory, large κ_D regime.^3–8 However, none of these theoretical studies reproduces the universal scaling reported in Ref. 23.

Here, we use a model framework to investigate long length scale phenomena in ionic fluids: the Coulomb or charge-frustrated Ising (FI) model,^10,28,29 a lattice model which accounts for both the long-range Coulomb and the short-range molecular interactions present in ionic fluids. While many statistical mechanical formulations of ionic correlations treat ions as charged hard spheres within the minimal Restrictive Primitive Model (RPM),^25,30 we call attention here to the importance of short-range attractive interactions, such as dispersion (or van der Waals) forces. The short-range molecular interaction is included in the FI model as a nearest neighbor “ferromagnetic” interaction (similar molecular species tend to attract), and we show that it controls the crossover between the small and large-κ_D regimes. Intuitively, the length scale of the short-range interaction, l_c, competes with that of the electrostatic interaction, 1/κ_D, and when the two become similar, the DH theory breaks down.

In the rest of the paper, we first introduce the FI model and its simple, continuum mean field form which is sufficient to predict a crossover in regimes between small and large κ_D. The mean field theory is only valid when the ratio of Coulomb and ferromagnetic interaction strengths is small, and fails as this ratio is increased. We then present our Monte Carlo simulation results. The simulation results are quantitatively well described by the mean field theory in the limit of low Coulomb interaction strength and are in qualitative agreement with them in other regimes. The simulations allow us to comment on the screening behavior in regimes inaccessible by the mean field theory. The simulations and the mean field theory also elucidate how a short-ranged attractive interaction can modify the screening behavior of ionic fluids, such as the crossover to the strong Coulomb coupling regime as well as the scaling of the screening length with the Debye constant seen in simulations.

We study the Coulomb or charge-frustrated Ising model on a three dimensional (d = 3) simple cubic lattice with each site occupied by a positive or negative charge as a simple model for ionic fluids. Since the positive and negative ions in an ionic fluid are chemically different species, the differences in their size or van der Waals interactions may lead to a preferential attractive interaction between like ions.¹⁰ In this model, the charges interact through a nearest-neighbor ferromagnetic Ising interaction, representing short-range molecular attraction between like charges, as well as the Coulomb interaction. The corresponding Hamiltonian is

with N being the number of lattice sites, $qi=qri=±1$ being the instantaneous charge density at site i located at position r_i, Q > 0 being the Coulomb interaction strength, r_ij = |r_i − r_j|, and

where J > 0 governs the strength of the Ising interaction. The ensemble average of the charge density ⟨q_r⟩ → 0 in the bulk and the unit of length is the lattice length, or nearest neighbor distance, a.

We can use the static charge structure factor $Sq(k)=qkq−k$⁠, to extract a screening length, where q_k is the Fourier transform of the instantaneous charge density q_r. In the continuum limit of the mean field theory, k ≪ a⁻¹, the static charge structure factor has the form¹⁰ 

with T being the temperature and the Boltzmann constant, k_B, set to 1, and ρ = 1/a³ in this study. The Ising critical temperature is defined by $T¯cI≡2dJ$ (overbarred variables are continuum mean field results). Inverse Fourier transforming the structure factor gives the charge–charge correlation function, $Gq(r,r′)=qrqr′$⁠. The continuum S_q(k) in Eq. (3) corresponds to, for an isotropic fluid at large r, the real space charge correlations given by

with A being a normalization constant dependent on the parameters T, J, and Q; ω being the spatial oscillation frequency; θ being a phase factor fixed by the electroneutrality condition; and κ_s being the calculated screening constant corresponding to the decay of charge correlations. The latter may differ from the Debye inverse screening length, which for the FI model is identified with

The phases and regimes of the FI mean field theory are revealed by examining how the inverse length scales κ_s and ω vary while changing the parameters Q, J, and T. In the rest of the paper, we fix the value of Q and treat κ_D as a parameter. By varying κ_D at fixed Q, we access different temperature regimes.

Long-range modulated order characterizes the phase below the critical point,^29,31 and so the FI continuum mean field critical temperature is simply given by the temperature at which κ_s → 0 from positive values

In this work, we focus on the fluid-like regime above the critical point where there is no real long-range order (κ_s > 0). There are two regimes above the critical point which are differentiated by the value of ω: when T is very high, ω = 0, while at intermediate temperatures, ω > 0. The transition between these two regimes occurs at

At high temperatures, $T>T¯*$⁠, or equivalently, small κ_D, charge correlations decay exponentially. Furthermore, the screening constant tends to the Debye constant when temperature is very large, $T≫T¯*$⁠: κ_s → κ_D. This small κ_D regime corresponds to low Coulomb coupling and is equivalent to the Debye–Hückel theory. For large κ_D, obtained at low temperatures (equivalent to strong coupling), oscillations with frequency ω appear in the charge correlations, while the inverse decay length κ_s decreases with κ_D,

where l_c is the mean field FI correlation length and

with

being the maximum screening constant, achieved at $T¯*$ [see the peak in Fig. 1, which occurs at the $κ¯D*$ corresponding to $T¯*$⁠, Eq. (7)]. Thus, in the FI mean field theory, $κ¯D*$ describes the transition between a DH-like regime with “gas-like” charge correlations and a second regime with “liquid-like” charge correlations, where κ_s has inverse dependence on temperature as in the DH regime $κsa∼κDa−1$⁠. The temperature dependence of κ_s in the “liquid-like” regime can be seen in Eq. (8) when $T¯cFI≪T<T¯*$⁠. The mean field prediction for κ_s is plotted against κ_D in Fig. 1 for ρQ/J = 0.5/a². The analogy with gas and liquid-like correlations is useful intuitively (and has been noted by others in connection with the so-called Fisher-Widom line⁷), but one important difference here is that the oscillation frequency is not fixed by the ion size and can instead vary significantly for different κ_D [see $ω¯$ given in Eq. (9)].

The correlation length associated with short-range Ising interactions, l_c, defines a molecular length scale in addition to the lattice size, a. In Fig. 1, we plot the inverse length-scales associated with the competing interactions of the FI model, namely, the Debye constant, κ_D, originating in Coulomb interactions, and the inverse FI correlation length, $lc−1$⁠, given in Eq. (8). The larger of the two length scales approximately determines the effective screening length, $κs−1$⁠, found within the FI model. The regime change of screening lengths in ionic fluids may then be understood in terms of these two competing length scales that are equal near the crossover point, $κ¯D*$⁠. At small κ_D, the correlations between ions are dominated by electrostatics, while at large κ_D, the short-range Ising correlations dominate. Importantly, even in the regime dominated by short-range interactions, electrostatics still plays a vital role, placing constraints on the system which appear as electroneutrality and higher moment conditions.^25,30

At large ρQ/J [$ρQ/J>d2/4πa2$], the continuum mean field theory breaks down, as noted by Grousson and Viot.³² One way the breakdown in the theory can be seen is through the FI critical temperature, Eq. (6), which becomes unphysically negative for large ρQ/J. The regime of validity can also be cast in terms of $κ¯s*$⁠, Eq. (10), $κ¯s*−1>a/d$ for validity. This form makes clear that the breakdown occurs when the minimum screening length for the system becomes similar to the lattice cell size. Grousson and Viot offer a correction by explicit treatment of the lattice,³² neglected here, and another route to improve the theory might be a more careful treatment of the finite size of ions. A third method to go beyond mean field theory, the incorporation of fluctuations, was considered as the correlation length is strongly renormalized near the critical temperature.^33,34 However, because the regimes we study are at temperatures far above criticality, the mean field results are not changed qualitatively. We use simulations of the FI model to investigate screening lengths and crossovers in the regime where the mean field theory breaks down.

We perform Monte Carlo simulations of the FI model to investigate its screening length behavior. We study parameter ranges strictly above the FI critical point.²⁹ We simulate a wide range of temperatures and extract the charge–charge correlation function, G_q(r), from simulations [see Fig. 2(a) for ρQ/J = 0.5/a²]. For small κ_D, $κD<κD*$⁠, the charge–charge correlation functions trend purely exponentially as predicted by the DH theory. For large κ_D, $κD>κD*$⁠, oscillations develop. By fitting the envelope of r|G_q(r)|, which has the form of a decaying exponential [mean field, or large r, form of G_q(r) shown in Eq. (4)], we can find the screening constant for a given κ_D. We plot the trending of the screening constant with κ_D for ρQ/J = 0.5/a² in blue dots in Fig. 3. For small κ_D, agreement between the DH theory, the continuum FI mean field theory, and the FI simulation is excellent. As κ_D increases beyond $κD*$⁠, estimates of the screening constant from both simulations and mean field theory begin to fall, with mean field scaling as in Eq. (8) and simulation scaling similarly, roughly as $κD−1$ near the screening constant peak. Overall, the agreement between the continuum mean field theory and simulation is excellent for small ρQ/J. The mean field theory is still reasonable at moderate ρQ/J, for example, see Fig. 4, where ρQ/J = 1/a².

Fitting the envelope of the charge–charge correlation function, G_q(r), works well to extract the screening constant except when the screening constant is large. In principle, the oscillation frequency can also be extracted by fitting a decaying oscillatory function, such as Eq. (4), to simulation data directly. However, due to constraints arising from the finite nature of the lattice, length scales extracted from such a fitting procedure can be error prone particularly in regimes where the length scale is comparable with the lattice size. We instead extract the oscillation frequency by first computing the charge–charge structure factor from simulation. We use the standard definition³⁵ 

from which S_q(k) can be easily computed; see Fig. 2(b) for some S_q(k) from simulation with ρQ/J = 0.5/a². We then fit the large wavelength or small-k region of S_q(k) using the inverse quartic form of the mean field expression in Eq. (3). As mentioned in the section titled Model, S_q(k) contains information about the length scales of the system, which can be extracted from the pole of the structure factor,

with κ_s and ω being the length scales appearing in the charge–charge correlation function, Eq. (4). Thus, fitting the small-k form to simulation S_q(k) allows us to extract estimates of both κ_s and ω from simulation.

The values of κ_s extracted from simulation using the large wavelength S_q(k) fits exhibit the same qualitative trends as those extracted from charge–charge correlation fits, see Fig. 3. Importantly, the scaling of the two regimes, κ_s ∼ κ_D when $κD≪κD*$ and $κs∼κD−1$ just above the regime changeover, is the same between the two methods. When κ_D is small, the S_q(k) fits underpredict the screening constant. Relative to the mean field, the S_q(k) fits also predict $κD*>κ¯D*$⁠. In the large κ_D regime, the S_q(k) fits overpredict the screening constant. The S_q(k) fit inverse length scales are essentially shifted to the right with respect to mean field and charge–charge correlation fits but capture the qualitative features.

Given the qualitative agreement between values of κ_s estimated from direct simulations and from the fitting method described above, it is reasonable to speculate that the oscillation frequencies extracted via S_q(k) small-k fits will capture the qualitative trends exhibited by the simulations. We compare the oscillation frequencies and screening constants extracted from the structure factor fits, to mean field predictions in Fig. 3. The oscillation frequency grows rapidly as κ_D increases past $κD*$ but saturates toward an asymptotic value as κ_D continues to increase, in line with the continuum mean field theory [$ω¯$ given in Eq. (9)].

We also simulate a range of ratios ρQ/J to extend our results beyond the continuum mean field theory which is only strictly valid for small ρQ/J.³² The short ranged ferromagnetic Ising interaction, described by J, causes spins which are alike to cluster, leading to a length scale, l_c, which acts as a molecular length scale aside from the lattice length, a. As recognized some time ago in the context of RPM models,^25,30 it is the frustration between a short-range length scale and the Coulomb length scale that results in non-DH behavior. While RPM models have a fixed molecular length scale, the hard sphere size, the FI model can potentially afford tunability of the molecular length scale, as J can be varied.

In Fig. 5, we plot the screening constant trending, extracted from large wavelength fits of the simulation S_q(k), for different ρQ/J ratios. We see that $κD*$ changes as ρQ/J is varied, but the same qualitative trends hold for all ρQ/J examined here. Namely, there are two regimes, one governed by the Debye constant, and the other governed by the inverse Ising correlation length analogous to the mean field prediction in Eqs. (5) and (8). The scaling of κ_s in the two regimes remains unchanged—κ_s ∼ κ_D when $κD≪κD*$ and $κs∼κD−1$ just after the regime changeover—despite changing the ratio ρQ/J. Thus, the two distinct regimes are robust even beyond the validity of the continuum mean field theory; within the range of parameters studied here, increasing ρQ/J monotonically increases $κD*$⁠. The division between the DH and overscreened regimes can thus be controlled by tuning J, as predicted in Eq. (7) and borne out in simulations in Fig. 5.

Finally, we consider the limiting case that exists when varying ρQ/J, namely, when J → 0. That limit allows us to make some connection with previous work on the lattice RPM^36,37 whose short-range interaction is purely repulsive. We find that two regimes occur in simulation for J = 0, just as in the J > 0 case, see Fig. 6. Note that the simple FI continuum mean field theory fails in this regime, predicting that the J = 0 case is identical to the Debye–Hückel theory for all values of κ_D. The simulation lattice plays a role directly analogous to the RPM hard sphere interaction, providing a sense of finite size to each ion.

The recent experimental discovery of universal scaling of the screening length, $κsa∼κDa−2$⁠, in concentrated electrolytes and ionic liquids has rekindled theoretical interest in the large κ_D or strong Coulomb coupling regime.²³ Past theoretical work based on the RPM of electrolytes using closure relations such as hypernetted chain approximations^3,25 and a generalization of the Debye charging process,²⁶ as well as a molecular dynamics simulation study of molten NaCl salt,⁴ suggest $κsa∼κDa−1/2$ for κ_D just above the peak $κD*$⁠. Considering additional effects such as the formation of Bjerrum ion pairs may modify the scaling to $κsa∼κDa−1$ within a Poisson–Boltzmann framework.³⁸ 

In this work, we focus on the properties of the FI model well above its critical point and find that it captures important features required to model the correlations of bulk ionic fluids. From simulations of the FI model, we find that $κsa∼κDa−1$ in the strong Coulomb coupling regime. The introduction of short length scale fluctuations affects only the temperature at which the crossover from the DH to the oscillatory regime occurs and leaves the scaling behavior unchanged. This scaling is different from the universal scaling experimentally observed in Ref. 23. However, it may be possible to alter the scaling of the FI model in the overscreened regime via simple modifications such as the introduction of defects in the lattice,¹¹ or creating asymmetry in the charge carriers, either in magnitude or shape.³⁹ These possibilities will be explored in future work. We also note that while the experimental universal scaling²³ and much previous theoretical work^3,25,40 place an emphasis on the ion size as a determining factor for the strong coupling regime, the ion size is not as simple to interpret in the FI model and appears to some extent through the Ising coupling J.

In conclusion, the FI model complements other theoretical techniques commonly used to describe ionic fluids, such as mean-field Poisson–Boltzmann theories,⁴¹ integral equations,²⁵ field theories⁴² or their hybrids,⁴³ and molecular simulations,⁴ and has the merit of reproducing the essential features of ionic correlations relatively simply. The FI model may be generalized to model surfaces and solvents in ionic fluids—which are systems of great current experimental interest.^19,23 Overall, the Coulomb-frustrated Ising model is an attractive framework for the study of long-range non-DH correlations in ionic fluids due to its simplicity and its capture of broad qualitative trends.

We acknowledge helpful conversations with Tom Witten. N.B.L. was primarily supported by the University of Chicago Materials Research Science and Engineering Center, which is funded by the National Science Foundation under Award No. DMR-1420709. S.V. acknowledges support from the Sloan Foundation. We gratefully acknowledge the University of Chicago Research Computing Center for computer time and technical assistance. D.V.T. also acknowledges support from the National Science Foundation (NSF Award No. DMR-1611371).

The Coulomb interaction is implemented using the Ewald summation technique.^44,45 The long-range part is precomputed at the start of a run since the separation between all lattice sites is fixed.²⁸ We use periodic boundary conditions in all three dimensions. Our simulation box has sides of length L = 32a with a = 1 being the lattice cell length. The lattice is initialized with an equal number of positive and negative charges. We use cluster moves which preserve the net charge of the system (⁠$∑jNqj=0$⁠) and greatly reduce the autocorrelation times at low temperatures, improving efficiency.²⁹ Monte Carlo move random numbers are generated using the PCG pseudo-random number generator.⁴⁶ Lattice trajectories were visualized using VMD.⁴⁷ 

We use fundamental requirements for statistical mechanical electrostatic systems as a check for our simulations. The Stillinger–Lovett second-moment (SL2) condition constrains the long-length scale fluctuations of a Coulomb system.³⁰ A formulation of the SL2 condition is that the charge structure factor tends to zero as k² for small k.²⁵ We have demonstrated that our simulation produces the required trend, see in particular Fig. 2(b). In addition, the high-T energy scaling of a Coulomb system must reduce to that of the Debye–Hückel theory U ∼ −T^−1/2.⁴⁸ We confirm that condition as well.