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.

## INTRODUCTION

The thermodynamic properties of ionic fluids are governed by long-range Coulomb interactions between ions^{1} 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^{14} 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 $\kappa D\u2261(4\pi \rho q2)/(\u03f5kBT)\u226a\sigma \u22121$, where *ρ* is the concentration of ions (per unit volume), *q* is the unit charge, *ϵ* is the dielectric constant of the electrolyte, *k*_{B}*T* 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, $\kappa s\u223c\kappa D\u22122$, 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, *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.

## MODEL

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

with *N* being the number of lattice sites, $qi=qri=\xb11$ 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)=qk\u2009q\u2212k$, 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*^{−1}, the static charge structure factor has the form^{10}

with *T* being the temperature and the Boltzmann constant, *k*_{B}, set to 1, and *ρ* = 1/*a*^{3} in this study. The Ising critical temperature is defined by $T\xafcI\u22612dJ$ (overbarred variables are continuum mean field results). Inverse Fourier transforming the structure factor gives the charge–charge correlation function, $Gq(r,r\u2032)=qr\u2009qr\u2032$. 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\xaf*$, or equivalently, small *κ*_{D}, charge correlations decay exponentially. Furthermore, the screening constant tends to the Debye constant when temperature is very large, $T\u226bT\xaf*$: *κ*_{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\xaf*$ [see the peak in Fig. 1, which occurs at the $\kappa \xafD*$ corresponding to $T\xaf*$, Eq. (7)]. Thus, in the FI mean field theory, $\kappa \xafD*$ 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 $\kappa sa\u223c\kappa Da\u22121$. The temperature dependence of *κ*_{s} in the “liquid-like” regime can be seen in Eq. (8) when $T\xafcFI\u226aT<T\xaf*$. The mean field prediction for *κ*_{s} is plotted against *κ*_{D} in Fig. 1 for *ρQ*/*J* = 0.5/*a*^{2}. 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^{7}), 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 $\omega \xaf$ 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\u22121$, given in Eq. (8). The larger of the two length scales approximately determines the effective screening length, $\kappa s\u22121$, 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, $\kappa \xafD*$. 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* [$\rho Q/J>d2/4\pi 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 $\kappa \xafs*$, Eq. (10), $\kappa \xafs*\u22121>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.

## SIMULATION

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, *G*_{q}(*r*), from simulations [see Fig. 2(a) for *ρQ*/*J* = 0.5/*a*^{2}]. For small *κ*_{D}, $\kappa D<\kappa D*$, the charge–charge correlation functions trend purely exponentially as predicted by the DH theory. For large *κ*_{D}, $\kappa D>\kappa 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*^{2} 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 $\kappa 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 $\kappa D\u22121$ 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*^{2}.

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^{35}

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*^{2}. 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 $\kappa D\u226a\kappa D*$ and $\kappa s\u223c\kappa D\u22121$ 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 $\kappa D*>\kappa \xafD*$. 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 $\kappa D*$ but saturates toward an asymptotic value as *κ*_{D} continues to increase, in line with the continuum mean field theory [$\omega \xaf$ 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, *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 $\kappa 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 $\kappa D\u226a\kappa D*$ and $\kappa s\u223c\kappa D\u22121$ 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 $\kappa 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.

## CONCLUSIONS

The recent experimental discovery of universal scaling of the screening length, $\kappa sa\u223c\kappa Da\u22122$, 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 approximations^{3,25} and a generalization of the Debye charging process,^{26} as well as a molecular dynamics simulation study of molten NaCl salt,^{4} suggest $\kappa sa\u223c\kappa Da\u22121/2$ for *κ*_{D} just above the peak $\kappa D*$. Considering additional effects such as the formation of Bjerrum ion pairs may modify the scaling to $\kappa sa\u223c\kappa Da\u22121$ 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 $\kappa sa\u223c\kappa Da\u22121$ 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 scaling^{23} 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,^{41} integral equations,^{25} field theories^{42} 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.

## ACKNOWLEDGMENTS

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).

### APPENDIX: METHODS

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* = 32*a* 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 ($\u2211jNqj=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 *k*^{2} 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.