Charge correlations in dense ionic fluids give rise to novel effects such as long-range screening and colloidal stabilization which are not predicted by the classic Debye–Hückel theory. We show that a Coulomb or charge-frustrated Ising model, which accounts for both long-range Coulomb and short-range molecular interactions, simply describes some of these ionic correlations. In particular, we obtain, at a mean field level and in simulations, a non-monotonic dependence of the screening length on the temperature. Using a combination of simulations and mean field theories, we study how the correlations in the various regimes are affected by the strength of the short ranged interactions.

The thermodynamic properties of ionic fluids are governed by long-range Coulomb interactions between ions1 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 biomolecules14 and soft materials.15 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,18 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.21 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,22 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, kBT 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κD2, 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.24 To take two examples, Attard uses a standard closure from the theory of liquids,25 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.26 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.1 Indeed, the manifestation of these correlations as oscillations in the charge density was predicted long ago by Kirkwood.27 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, lc, 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.10 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

H=12iNjiNqiQrijJijqj,
(1)

with N being the number of lattice sites, qi=qri=±1 being the instantaneous charge density at site i located at position ri, Q > 0 being the Coulomb interaction strength, rij = |rirj|, and

Jij=Ji,j  nearest neighbors,0otherwise,
(2)

where J > 0 governs the strength of the Ising interaction. The ensemble average of the charge density ⟨qr⟩ → 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)=qkqk, to extract a screening length, where qk is the Fourier transform of the instantaneous charge density qr. In the continuum limit of the mean field theory, ka−1, the static charge structure factor has the form10 

ρ2Sq(k)/T=k2/a2Jk4+T2dJk2+4πρQ,
(3)

with T being the temperature and the Boltzmann constant, kB, set to 1, and ρ = 1/a3 in this study. The Ising critical temperature is defined by T¯cI2dJ (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 Sq(k) in Eq. (3) corresponds to, for an isotropic fluid at large r, the real space charge correlations given by

Gq(r)=A4πrexp(κsr)cos(ωr+θ),
(4)

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

κD4πρQT.
(5)

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

T¯cFI=T¯cI16πa2JρQ.
(6)

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

T¯*=T¯cI+16πa2JρQ.
(7)

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, TT¯*: κ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,

κ¯s=12lc¯TT¯cFI4a2J,T<T¯*,
(8)

where lc is the mean field FI correlation length and

ω¯=κ¯s*2κ¯s2,T<T¯*,
(9)

with

κ¯s*4πρQa2J1/4
(10)

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κDa1. The temperature dependence of κs in the “liquid-like” regime can be seen in Eq. (8) when T¯cFIT<T¯*. The mean field prediction for κs is plotted against κD in Fig. 1 for ρQ/J = 0.5/a2. 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 line7), 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)].

FIG. 1.

Mean field screening constant, κs, identified with the inverse decay length of charge correlations, displays non-monotonic trend as the Debye screening constant, κD=4πρQ/T, [Eq. (5)] is increased, plotted here for ρQ/J = 0.5/a2. The solid black line shows the predicted screening constant, κs, in the two regimes. Note the inverse dependence of κs on T in the two regimes [see Eq. (8)]. Near, but slightly above the regime change, the screening constant from simulation shows an apparent scaling κsaκDa1. The dashed line shows the Debye constant κD, and the dotted line shows the temperature scaling of the inverse Ising correlation length T/(a2J)1/lc. The dashed–dotted line is a second inverse length scale which goes as 1/lc for small κD; it merges with κs at the regime change κ¯D*, which also marks the peak in the screening constant, κ¯s*.

FIG. 1.

Mean field screening constant, κs, identified with the inverse decay length of charge correlations, displays non-monotonic trend as the Debye screening constant, κD=4πρQ/T, [Eq. (5)] is increased, plotted here for ρQ/J = 0.5/a2. The solid black line shows the predicted screening constant, κs, in the two regimes. Note the inverse dependence of κs on T in the two regimes [see Eq. (8)]. Near, but slightly above the regime change, the screening constant from simulation shows an apparent scaling κsaκDa1. The dashed line shows the Debye constant κD, and the dotted line shows the temperature scaling of the inverse Ising correlation length T/(a2J)1/lc. The dashed–dotted line is a second inverse length scale which goes as 1/lc for small κD; it merges with κs at the regime change κ¯D*, which also marks the peak in the screening constant, κ¯s*.

Close modal

The correlation length associated with short-range Ising interactions, lc, 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, lc1, given in Eq. (8). The larger of the two length scales approximately determines the effective screening length, κs1, 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.32 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,32 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.29 We simulate a wide range of temperatures and extract the charge–charge correlation function, Gq(r), from simulations [see Fig. 2(a) for ρQ/J = 0.5/a2]. 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|Gq(r)|, which has the form of a decaying exponential [mean field, or large r, form of Gq(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/a2 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 κD1 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/a2.

FIG. 2.

Spatial correlations in the FI model for various inverse Debye screening lengths, κD, for the parameter ρQ/J = 0.5/a2. (a) Absolute values of charge–charge correlation functions, r|Gq(r)|, plotted on log-linear scale for various Debye constants. For κDκD*, the correlations decay purely exponentially as shown in the bottom two plots, while oscillations appear when κDκD*, see the top two plots. The solid black lines correspond to the envelope of these functions from which κs can be extracted. The dotted black line is the DH prediction for the decay of correlations. (b) Structure factors scaled by temperature, Sq(k)/T, for various κD. j is an integer in [0, L). For small k, the structure factors scale as k2 (solid black line). For κDκD*, Sq(k) plateaus when k becomes large, but as κD increases, oscillations appear. The peak at k ∼ 1 shifts toward larger k with increasing Debye constant. The largest k value peak corresponds to the lattice length a.

FIG. 2.

Spatial correlations in the FI model for various inverse Debye screening lengths, κD, for the parameter ρQ/J = 0.5/a2. (a) Absolute values of charge–charge correlation functions, r|Gq(r)|, plotted on log-linear scale for various Debye constants. For κDκD*, the correlations decay purely exponentially as shown in the bottom two plots, while oscillations appear when κDκD*, see the top two plots. The solid black lines correspond to the envelope of these functions from which κs can be extracted. The dotted black line is the DH prediction for the decay of correlations. (b) Structure factors scaled by temperature, Sq(k)/T, for various κD. j is an integer in [0, L). For small k, the structure factors scale as k2 (solid black line). For κDκD*, Sq(k) plateaus when k becomes large, but as κD increases, oscillations appear. The peak at k ∼ 1 shifts toward larger k with increasing Debye constant. The largest k value peak corresponds to the lattice length a.

Close modal
FIG. 3.

Screening constant, κs, for different extraction methods and oscillation frequency, all from simulation for ρQ/J = 0.5/a2 and compared with theory. Solid and dashed black lines show mean field theory prediction for screening length and oscillation frequency, respectively. Blue dots show screening constant extracted from envelope fits of charge–charge correlation functions [method shown in Fig. 2(a)]. Red triangles show screening constant, while green squares show oscillation frequency extracted from small-k course of simulation Sq(k) (see section titled Simulation). The length scales from Sq(k) fits consistently overestimate the length scale in the small κD regime and underestimate it in the large κD regime.

FIG. 3.

Screening constant, κs, for different extraction methods and oscillation frequency, all from simulation for ρQ/J = 0.5/a2 and compared with theory. Solid and dashed black lines show mean field theory prediction for screening length and oscillation frequency, respectively. Blue dots show screening constant extracted from envelope fits of charge–charge correlation functions [method shown in Fig. 2(a)]. Red triangles show screening constant, while green squares show oscillation frequency extracted from small-k course of simulation Sq(k) (see section titled Simulation). The length scales from Sq(k) fits consistently overestimate the length scale in the small κD regime and underestimate it in the large κD regime.

Close modal
FIG. 4.

Screening constant, κs, displays non-monotonic trend as κD is increased, shown here for ρQ/J = 1/a2. Solid black line is continuum mean field theory prediction. Blue dots are screening constants extracted from the envelope of charge–charge correlation functions, Gq(r), in simulations. The effect of the negative T¯cFI is visible in the slight positive curvature of the mean field prediction when κD>κ¯D*. Near, but slightly above the regime change, simulation κsaκDa1.

FIG. 4.

Screening constant, κs, displays non-monotonic trend as κD is increased, shown here for ρQ/J = 1/a2. Solid black line is continuum mean field theory prediction. Blue dots are screening constants extracted from the envelope of charge–charge correlation functions, Gq(r), in simulations. The effect of the negative T¯cFI is visible in the slight positive curvature of the mean field prediction when κD>κ¯D*. Near, but slightly above the regime change, simulation κsaκDa1.

Close modal

Fitting the envelope of the charge–charge correlation function, Gq(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 definition35 

Sq(k)=1Nj,lqjqlexp2πiLkrjrl,
(11)

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

k0=ω+iκs,
(12)

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

The values of κs extracted from simulation using the large wavelength Sq(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κD1 just above the regime changeover, is the same between the two methods. When κD is small, the Sq(k) fits underpredict the screening constant. Relative to the mean field, the Sq(k) fits also predict κD*>κ¯D*. In the large κD regime, the Sq(k) fits overpredict the screening constant. The Sq(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 Sq(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.32 The short ranged ferromagnetic Ising interaction, described by J, causes spins which are alike to cluster, leading to a length scale, lc, 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 Sq(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κD1 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.

FIG. 5.

The screening constant, κs, against κD for different ρQ/J ratios. We extract κs here using the small-k course of Sq(k) discussed in the section titled Simulation. Increasing ρQ/J shifts κD* to the right, also increasing the maximum screening constant, κs*. Near but slightly above the regime change, simulation κsaκDa1 for each ρQ/J [the dotted lines show the scaling T/JκDa1 for each parameter set].

FIG. 5.

The screening constant, κs, against κD for different ρQ/J ratios. We extract κs here using the small-k course of Sq(k) discussed in the section titled Simulation. Increasing ρQ/J shifts κD* to the right, also increasing the maximum screening constant, κs*. Near but slightly above the regime change, simulation κsaκDa1 for each ρQ/J [the dotted lines show the scaling T/JκDa1 for each parameter set].

Close modal

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 RPM36,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.

FIG. 6.

Screening constant, κs, displays non-monotonic trend as κD [Eq. (5)] is increased, shown here for ρQ = 1/a2 and J = 0. Dashed black line is the Debye constant, κD, which is also the prediction of the continuum mean field theory presented in the section titled Model when J = 0. Blue dots are screening constants extracted from the envelope of simulation charge–charge correlation functions, Gq(r). Note that the domain and range of this plot differ from previous κs vs κD plots in this paper.

FIG. 6.

Screening constant, κs, displays non-monotonic trend as κD [Eq. (5)] is increased, shown here for ρQ = 1/a2 and J = 0. Dashed black line is the Debye constant, κD, which is also the prediction of the continuum mean field theory presented in the section titled Model when J = 0. Blue dots are screening constants extracted from the envelope of simulation charge–charge correlation functions, Gq(r). Note that the domain and range of this plot differ from previous κs vs κD plots in this paper.

Close modal

The recent experimental discovery of universal scaling of the screening length, κsaκDa2, in concentrated electrolytes and ionic liquids has rekindled theoretical interest in the large κD or strong Coulomb coupling regime.23 Past theoretical work based on the RPM of electrolytes using closure relations such as hypernetted chain approximations3,25 and a generalization of the Debye charging process,26 as well as a molecular dynamics simulation study of molten NaCl salt,4 suggest κsaκDa1/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κDa1 within a Poisson–Boltzmann framework.38 

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κDa1 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,11 or creating asymmetry in the charge carriers, either in magnitude or shape.39 These possibilities will be explored in future work. We also note that while the experimental universal scaling23 and much previous theoretical work3,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,41 integral equations,25 field theories42 or their hybrids,43 and molecular simulations,4 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.28 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.29 Monte Carlo move random numbers are generated using the PCG pseudo-random number generator.46 Lattice trajectories were visualized using VMD.47 

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.30 A formulation of the SL2 condition is that the charge structure factor tends to zero as k2 for small k.25 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.48 We confirm that condition as well.

1.
2.
A. Y.
Grosberg
,
T. T.
Nguyen
, and
B. I.
Shklovskii
,
Rev. Mod. Phys.
74
,
329
(
2002
).
3.
J.
Ennis
,
R.
Kjellander
, and
D.
Mitchell
,
J. Chem. Phys.
102
(
2
),
975
(
1995
).
4.
P.
Keblinski
,
J.
Eggebrecht
,
D.
Wolf
, and
S. R.
Phillpot
,
J. Chem. Phys.
113
,
282
(
2000
).
5.
J.
Janeček
and
R. R.
Netz
,
J. Chem. Phys.
130
,
074502
(
2009
).
6.
M. Z.
Bazant
,
B. D.
Storey
, and
A. A.
Kornyshev
,
Phys. Rev. Lett.
106
,
046102
(
2011
).
7.
D.
Frydel
,
J. Chem. Phys.
145
,
184703
(
2016
).
8.
N.
Gavish
,
D.
Elad
, and
A.
Yochelis
,
J. Phys. Chem. Lett.
9
,
36
(
2018
).
9.
D. J.
Bozym
,
B.
Uralcan
,
D. T.
Limmer
,
M. A.
Pope
,
N. J.
Szamreta
,
P. G.
Debenedetti
, and
I. A.
Aksay
,
J. Phys. Chem. Lett.
6
,
2644
(
2015
).
10.
D. T.
Limmer
,
Phys. Rev. Lett.
115
,
256102
(
2015
).
11.
B.
Uralcan
,
I. A.
Aksay
,
P. G.
Debenedetti
, and
D. T.
Limmer
,
J. Phys. Chem. Lett.
7
,
2333
(
2016
).
12.
M. A.
Gebbie
,
M.
Valtiner
,
X.
Banquy
,
E. T.
Fox
,
W. A.
Henderson
, and
J. N.
Israelachvili
,
Proc. Natl. Acad. Sci. U. S. A.
110
,
9674
(
2013
).
13.
A. M.
Smith
,
A. A.
Lee
, and
S.
Perkin
,
J. Phys. Chem. Lett.
7
,
2157
(
2016
).
14.
W. M.
Gelbart
,
R. F.
Bruinsma
,
P. A.
Pincus
, and
V. A.
Parsegian
,
Phys. Today
53
(
9
),
38
(
2000
).
15.
Y.
Li
,
M.
Girard
,
M.
Shen
,
J. A.
Millan
, and
M.
Olvera de la Cruz
,
Proc. Natl. Acad. Sci. U. S. A.
114
,
11838
(
2017
).
16.
17.
H.
Riazi
,
S.
Mesgari
,
N. A.
Ahmed
, and
R. A.
Taylor
,
Int. J. Heat Mass Transfer
94
,
254
(
2016
).
18.
K.
Gkagkas
,
V.
Ponnuchamy
,
M.
Dai
, and
I.
Stankovi
, in
43rd Leeds–Lyon Symposium on Tribology 2016
[
Tribol. Int.
113
,
83
(
2017
)].
19.
H.
Zhang
,
K.
Dasbiswas
,
N. B.
Ludwig
,
G.
Han
,
B.
Lee
,
S.
Vaikuntanathan
, and
D. V.
Talapin
,
Nature
542
,
328
(
2017
).
20.
V.
Srivastava
,
W.
Liu
,
E. M.
Janke
,
V.
Kamysbayev
,
A. S.
Filatov
,
C.-J.
Sun
,
B.
Lee
,
T.
Rajh
,
R. D.
Schaller
, and
D. V.
Talapin
,
Nano Lett.
17
,
2094
(
2017
).
21.
D.
Ben-Yaakov
,
D.
Andelman
,
R.
Podgornik
, and
D.
Harries
,
Curr. Opin. Colloid Interface Sci.
16
,
542
(
2011
).
22.
P.
Debye
and
E.
Hückel
,
Phys. Z.
24
,
185
(
1923
).
23.
A. A.
Lee
,
C. S.
Perez-Martinez
,
A. M.
Smith
, and
S.
Perkin
,
Phys. Rev. Lett.
119
,
026002
(
2017
).
24.
L. M.
Varela
,
M.
Garcia
, and
V.
Mosquera
,
Phys. Rep.
382
,
1
(
2003
).
25.
26.
B. P.
Lee
and
M. E.
Fisher
,
Europhys. Lett.
39
,
611
(
1997
).
27.
J. G.
Kirkwood
,
J. Chem. Phys.
2
,
767
(
1934
).
28.
D.
Wu
,
D.
Chandler
, and
B.
Smit
,
J. Phys. Chem.
96
,
4077
(
1992
).
29.
M.
Grousson
,
G.
Tarjus
, and
P.
Viot
,
Phys. Rev. E
64
,
036109
(
2001
).
30.
F. H.
Stillinger
, Jr.
and
R.
Lovett
,
J. Chem. Phys.
48
,
3858
(
1968
).
31.
M.
Seul
and
D.
Andelman
,
Science
267
,
476
(
1995
).
32.
M.
Grousson
,
G.
Tarjus
, and
P.
Viot
,
Phys. Rev. E
62
,
7781
(
2000
).
33.
S. A.
Brazovskiǐ
,
Sov. J. Exp. Theor. Phys.
41
,
85
(
1975
).
34.
G. H.
Fredrickson
and
K.
Binder
,
J. Chem. Phys.
91
,
7265
(
1989
).
35.
J. P.
Hansen
and
I. R.
McDonald
,
Theory of Simple Liquids
(
Elsevier
,
2013
).
36.
A. Z.
Panagiotopoulos
and
S. K.
Kumar
,
Phys. Rev. Lett.
83
,
2981
(
1999
).
37.
A.
Diehl
and
A. Z.
Panagiotopoulos
,
Phys. Rev. E
71
,
046118
(
2005
).
38.
J.
Zwanikken
and
R.
van Roij
,
J. Phys.: Condens. Matter
21
,
424102
(
2009
).
39.
M.
Ding
,
Y.
Liang
,
B.-S.
Lu
, and
X.
Xing
,
J. Stat. Phys.
165
,
970
(
2016
).
40.
R. L.
de Carvalho
and
R.
Evans
,
Mol. Phys.
83
,
619
(
1994
).
41.
J. N.
Israelachvili
,
Intermolecular and Surface Forces
(
Academic Press
,
2011
).
42.
R. R.
Netz
and
H.
Orland
,
Europhys. Lett.
45
,
726
(
1999
).
43.
C. E.
Sing
,
J. W.
Zwanikken
, and
M. O.
de la Cruz
,
Phys. Rev. Lett.
111
,
168303
(
2013
).
44.
45.
D.
Frenkel
and
B.
Smit
,
Understanding Molecular Simulations
(
Academic Press
,
San Diego
,
2002
).
46.
M. E.
O’Neill
, “
PCG: A family of simple fast space-efficient statistically good algorithms for random number generation
,” Technical Report No. HMC-CS-2014-0905,
2014
.
47.
W.
Humphrey
,
A.
Dalke
, and
K.
Schulten
,
J. Mol. Graphics
14
,
33
(
1996
).