Beyond the Debye–Hückel limit: Toward a general theory for concentrated electrolytes

Previous theoretical treatments of electrolyte solutions have been largely confined to symmetric monovalent electrolyte systems that are known to support a KT. Analytical solutions of these simplified models predict the inverse length scale Q₀ will emerge precisely when the slope of a₀ vs concentration changes sign.^3,8–13 In addition, a₀ begins to deviate from the inverse DH length, κ_D, at concentrations below the KT, and a₀ becomes equal to κ_D at low concentrations. In what follows, we refer to this as a classic KT, where at a given concentration only one dominant screening length, a₀, exists.^8,11,14 Underscreening is defined as the concentration range where a₀ < κ_D (see Fig. 1). Since the seminal work of Kirkwood, modern views of the charge–charge correlations have been put forth.^{8,11–13,15–21} A goal of this study is to unify these approaches within a generalized formulation using a single function, Σ(Q) (see Fig. 2), and show the implications on KTs and breakdown of classic KTs for multivalent electrolytes. Here, we develop both the theory and numerical procedure that can predict Σ(Q), a correlation–interaction term beyond DH, for concentrated electrolytes.

An additional impetus for a more quantitative and conceptual approach to concentrated electrolytes has been motivated by two distinct experimental observations. First, surface force apparatus (SFA)^23–26 and atomic force microscopy (AFM) measurements by Hjalmarsson et al.²⁷ have reported large notional screening lengths in ionic fluids above KTs.^23–25 In contrast, other AFM experiments did not detect large screening lengths.^28,29 Although a large screening length has been obtained from surface techniques, it has been attributed to phenomena dominated by the bulk properties of electrolytes²⁵ and referred to as the “anomalous” screening length in the recent literature.^21,30 Although recent theoretical and simulation studies have reported on large screening lengths in Coulomb liquids,^21,31,32 a consensus on the precise molecular-scale phenomena responsible for anomalous screening remains an active area of research.^5,6,33,34

In contrast to the surface techniques, small-angle x-ray scattering (SAXS) measurements suggest a different interpretation of charge correlations under bulk homogeneous conditions. Specifically, the mapping of the SAXS experimental signal in the prepeak region (low-Q) to charge–charge correlations predicts underscreening and KTs.¹⁴ In the SAXS measurements, underscreening was shown to be consistent with standard liquid state theory estimates and support a picture where large, anomalous screening lengths do not dominate in concentrated bulk electrolytes.¹⁴ Therefore, it is imperative to understand and contrast the theoretical frameworks that can describe charge correlations beyond DH in Coulombic liquids.

Despite many theories being developed for electrolyte solutions to go beyond the DH limit,^{8,13,15,16,20,35,36} we will demonstrate that quantitative connections to the SAXS measurements require the incorporation of molecular detail that cannot be accessed through practical analytical treatments. Our results suggest that although current theories of electrolyte solutions^{8,13,15,16,20,35–37} are able to capture the phenomena of underscreening, they are unable to explain the experimentally observed concentration dependencies and the dominant inverse length scales a₀ and Q₀ that describe the emergent correlations in real multivalent electrolyte solutions.¹⁴ A necessary ingredient in a theory for concentrated electrolytes is the presence of a non-Coulombic short-range (SR) interaction among ions that competes with the slowly-varying LR Coulomb interaction. This SR interaction can be as simple as a hard sphere (HS) interaction used in numerical and simulation approaches of restricted primitive models (RPMs) of electrolytes.^21,38–40 Here, we generalize the SR interaction, u^SR(r), to include molecular details from condensed phase simulations using both classical and ab initio-based interaction potentials from quantum density functional theory (qDFT). We also present a theoretical framework that will allow us to differentiate between choices of molecular-informed SR interactions while preserving the conceptual simplicity of the “frustrated charge” models^16,41–43 (discussed below), providing a seamless connection to previous formulations of concentrated electrolytes.

The general framework upon which we develop our theory is schematically presented in Fig. 2, showing the progression of complexity from intrinsic to concentrated. In the intrinsic limit, i.e., an infinitely dilute dielectric continuum, two ions with charges q_α and q_β interact via the LR Coulomb interaction, $uLRrαβ=qαqβ/ϵrαβ$⁠, where ϵ is the static dielectric constant and r_αβ is the distance between ions. This gives rise to the well-known Coulomb interaction u(Q) in reciprocal Q-space (the left panel in Fig. 2). The effective response function for dilute systems in the DH limit, $χDHQ$⁠, can be derived by using the intrinsic interaction u(Q) in conjunction with linear response to get $χDHQ=κD2/(Q2+κD2)$ (the middle panel in Fig. 2).²² Moreover, $χDHQ$ is related to h_zz(r) via Fourier Transformation (FT), yielding the DH result of h_zz(r) ∼ exp(−κ_Dr)/r. In the concentrated limit (the right panel in Fig. 2), we show that the effective response function beyond the DH limit, χ(Q), can be formulated in terms of a single function $ΣQ$⁠. In what follows we show how $ΣQ$ is related to the charge–charge structure factor, S_zz(Q), defined by S_zz(Q) = ∑_αβz_αz_βS_αβ(Q). Here, $Sαβ(Q)=cαδαβ+cαcβnĥαβ(Q)$ is the partial structure factor, where the Greek subscripts refer only to the ions and $ĥαβ(Q)$ is the FT of h_αβ(r).⁴ 

Importantly, we develop a Gaussian field theory and formulas for $ΣQ$ that arise from a simple coupling between number and charge density (Sec. III). The theory presented in this perspective is general enough to describe charge correlation in realistic descriptions of electrolytes in addition to recovering the DH limit and the well-studied correlations of primitive models of electrolytes (the so-called RPM and their variants). An outcome of this research is a formulation that can be used to directly interpret the SAXS experimental data in terms of Σ(Q) and connect to the choice of the molecular representation that determines u^SR(r). Moreover, this work provides the necessary theoretical framework to connect all standard numerical and theoretical models of electrolytes and to demonstrate the importance of molecular details in providing robust formulations of charge–charge correlations that relate local molecular structure to screening in aqueous electrolytes.

This perspective is organized as follows: in Sec. II, we review important relevant theories of electrolyte solutions and discuss how they can be characterized via Σ(Q). We use simple RPM models to explain the general form Σ(Q) and demonstrate that Σ(Q) cannot be formulated as a simple power series in Q. Utilizing a Gaussian field theory in terms of the number density and charge density fields, we provide a theoretical formulation of the general form Σ(Q) in Sec. III. Importantly, it will be shown that Σ(Q) can be solely determined from SR direct correlations that can be used to quantitatively reproduce numerical simulations of various models. In Sec. IV, we present how to extract the pole structure from Σ(Q) and the resulting KT for real electrolyte solutions, predicting screening lengths to be about $∼1$ nm for concentrations for multivalent electrolytes above the KTs. In Sec. V, we conclude the study and present an outlook of the approach and future directions.

In this perspective, we will demonstrate that the entirety of correlations in concentrated electrolyte solutions can be described with a single function $ΣQ$⁠. $ΣQ$ can be numerically computed from simulation and can be determined theoretically under reasonable assumptions that are guided by liquid state theory. Importantly, our theoretical analysis suggests that the precise functional form of $ΣQ$ necessitates the coupling of both a number density field, which describes non-Coulombic SR packing, to a charge density field. The coupling of both charge and number density fields through a Gaussian field theory provides a useful framework to understand the origins of various correlations in $ΣQ$⁠. In addition, $ΣQ$ can be used to compare and contrast all standard theories of electrolyte solutions implicating the precise molecular nature of u^SR(r).

A formulation of the theory of concentrated electrolytes is a subject of broad scientific interest. The vast majority of theories to describe the collective effects of screening beyond the DH theory treat water as a continuum dielectric and find practical statistical mechanical relations between the local densities of the ions and the electrostatic potentials that satisfy the Poisson equations.^{4,8,9,11,12,17,20,37,44–48} In this regard, the electrostatic potentials have been obtained by the Green function techniques within the linearized Poisson–Boltzmann framework.^8,17 One can then construct electrostatic free energy functionals in terms of a single charge density, $ρZr$⁠, yielding a direct route to charge–charge correlations to determine a₀ and Q₀, yielding the classic KT shown in Fig. 1.^8,17

Another distinct approach is based on the modified Poisson–Boltzmann (MPB) closure that was introduced in Ref. 49 and successfully predicted classic KTs for 1:1 electrolytes. The advantage of the MPB formalism is that one can apply various hierarchies^{9,37,44,45,49,50} to access the concentration dependent radial distribution functions. A recent study by Outhwaite and Bhuiyan suggests that the present MPB formulation in a given hierarchy predicts classic KTs even for multivalent electrolyte solutions.³⁷ 

The liquid state theory with RPMs has also provided important insights into charge correlations in electrolytes.^{38–40,51,52} Here, the non-Coulombic u^SR(r) is given by a HS interaction that competes with the LR Coulomb interaction. These primitive models of electrolytes utilize the Orstein–Zernike (OZ) equation of liquid state theory in conjunction with the HyperNetted-Chain (HNC) or Percus–Yevick (PY) closures, which have shed light on deviations from DH theory.^38–40 Although these numerical methods are robust, there is little to be gleaned conceptually regarding the essential physics that drives charge correlations beyond DH.

Importantly, all of the above-mentioned approaches (in addition to other formulations similar in spirit^11,13,35), although seemingly disparate, yield substantially similar outcomes of the charge–charge correlations, screening, and the concomitant KTs for 1-1 electrolyte solutions that follow the picture provided in Fig. 1.^{8,12,13,17,20,37,45} One could build upon these approaches by formulating a perturbation or diagrammatic theory to connect to the framework of classical density functional theory. Rather, we show convincingly that the salient physics to explain the complexity of Σ(Q) can be formulated within a coupling of two Gaussian fields representing charge and number density^53,54 (see Sec. III). This Gaussian field theory provides a direct route to structure factors in addition to a direct connection to liquid state theory and both numerical and primitive models of electrolytes through a molecular-informed effective 2-body interaction⁵⁵ that will be discussed in detail below. In addition, the theory presented here is flexible enough to utilize input from any theory that provides a direct correlation function (say, from a classical density functional theory) to predict the form of Σ(Q) and provide a complete description of charge–charge correlations, the screening lengths, and the breakdown of classic KTs for multivalent electrolytes.

More insight into the behavior of the Σ(Q) function can be gleaned from the simple RPM models of electrolyte solutions whose analytical approaches can provide a conceptual picture but are usually limited to describing model electrolytes where all components are represented by identical sphere sizes. Examples include Kjellander’s dressed ion theories of RPM models in which the dressed ion has an effective charge determined by a diffuse charge associated with an ion.^11,20 Lee and Fisher provided an analytical expression for S_zz(Q) enforcing soft oscillatory modes in the charge density that are found in solutions to the FCI models.^{8,17,43,56,58} Xiao et al. followed Kjellander’s work based within the mean-spherical approximation (MSA) framework and assumed that the electric potential could be a linear combination of solutions from a DH theory.^35,36 Recently, Adar presented a theory of electrolyte solutions in which the two-body Coulomb interactions were modified.¹³ 

Figure 3 depicts the numerically extracted Σ(Q) for symmetric 1-1 RPM electrolyte solutions with HS diameters of 4.25 Å at 6M from various theoretical treatments. Here, the HNC results are used as the standard reference, as they give consistent results with independent all-atom molecular dynamics simulations (see  Appendix C 4 and Fig. 12). Although the FCI and RPM are known to provide qualitatively similar screening and KTs, the form of Σ(Q) from the numerical models is distinct from that obtained from Eq. (5). Closer examination of Fig. 3 reveals that all primitive models depict a similar peak structure at low-Q qualitatively similar to realistic models of electrolytes that will be developed in Sec. II C.

To this end, making use of the MSA expression presented by Xiao et al.^35,36 for Σ(Q) (see  Appendix A), one observes better agreement, with no fitting parameters, to the numerical HNC predictions of the RPM results over a wider Q range (see Fig. 3). Similarly, the Σ(Q) from the approach of Lee and Fisher⁸ also shows reasonable agreement with the numerical RPM results over the entire Q-region. Noting that Adar et al.’s theory does not consider any coupling of number and charge density fields from their restricted Coulomb interactions,¹³ their approach yields a significant overestimation of Σ(Q) with respect to the numerical HNC results in the low Q region. When different variants of Adar’s model are considered, we find a wide range of results, including a significant underestimation of Σ(Q) in the low Q region (see  Appendix A).

Importantly, the inset of Fig. 3 suggests that a two parameter fit to a simple TS distribution (i.e., ion correlations predicted by the FCI) can capture the behavior of Σ(Q) only in a limited selected range of Q < 1.2 Å⁻¹. This provides a hint as to why these primitive models of electrolytes yield qualitatively similar screening lengths as a function of concentration in this important asymptotic regime. However, a quantitative fit over the entire low Q-range requires a Padé-type function [Eq. (G2)]. In general, it is clear that a simple polynomial in even powers of Q (⁠$∑i=04biQ2i$ with the same number of parameters as the Padé-type function) cannot adequately approximate Σ(Q) over a wide range of Q-space. A general functional form of Σ(Q) that justifies the aforementioned Padé approximant will be derived in Sec. III and constitutes the main theoretical result of this perspective [see Eqs. (18)–(20)].

Figure 4 shows the size dependencies of Σ(Q) for symmetric 1-1 RPM electrolytes at 4.5 M, again solved numerically. HS diameters ranging from 2 to 6 Å are shown. As can be seen, the magnitudes of Σ(Q) at Q < 0.5 Å⁻¹ are larger for larger ions, showing greater deviations from the DH limit, tracking with our understanding that small HS can be well represented by simple theory.^{17,44,49,59–63} These deviations are also reflected in the susceptibilities $χQ$ [see Fig. 13(a)], which also show larger amplitudes for larger ions. Figure 4 demonstrates that Σ(Q) has damped oscillatory behavior in Q-space with the oscillation wavelengths decreasing as the ion gets larger. This is due to packing and the emergence of extended structure in the electrolytes. Interestingly, Σ(Q) values for the largest studied ion size (a = 6 Å) can become negative around 2.3 and 3.3 Å⁻¹ in the Q < 4.0 Å⁻¹ region, resulting in significant changes in the slope of $χQ$ around 2.3 and 3.3 Å⁻¹ [see the inset of Fig. 13(a) in the  Appendix D].

In this section, we extract Σ(Q) for realistic molecular models of electrolyte solutions and demonstrate the sensitivity of Σ(Q) to the details of u^SR(r). For such realistic models of electrolyte solutions, the exact solutions of Σ(Q) in the low-Q region cannot be reliably estimated directly from all-atom molecular simulations with explicit water due to insufficient sampling required at the length scales relevant for the collective effects of screening (low-Q limit). Therefore, as alluded to above, in the low-Q limit, the collective effects of screening must be estimated from reduced theoretical methods. The interested reader is referred to  Appendix C for details. Here we present them concisely.

In previous studies, we and others have shown the efficacy of this approach for predicting a wide range of collective thermodynamic properties of concentrated electrolytes and comparing well with all-atom molecular treatments with explicit water.^14,64,66 For the remainder of the text, we refer to the use of $uSRr$ in conjunction with integral equations as the molecular-informed HNC (mHNC) approach to distinguish it from previous numerical procedures stated earlier for studies utilizing HS ions, which we will refer to as the HNC method. The intrinsic SR interactions, u^SR(r), for various models of SrCl₂ and NaCl electrolyte solutions are shown in Figs. 5 and 6, which are obtained from classical neural network potential (NNP) trained on qDFT^67,68 and qDFT representations of interaction (see  Appendix C for details).

As can be seen in the top panel of Fig. 5, there are significant differences in u^SR(r) based on the form of molecular interactions used to compute them. In particular, u^SR(r) potentials are shown for NNP, classical Kirkwood–Buff (KBFF), and classical Smith–Dang (SD) models between Na–Na, Na–Cl, and Cl–Cl ion pairs, where u^SR(r) is extracted from the infinite dilution PMFs after removing u^LR(r). Importantly, the peak positions of u^SR(r) and the NaCl ion-pair for all three models are in qualitative agreement. Quantitative differences in the like–like ion u^SR(r) obtained with NNP are softer than those obtained from the classical force fields (FFs). This is consistent with the current literature where such differences between classical FFs and NNP representations of interaction require the effects of explicit polarization that can be corrected by charge scaling to account for the electronic polarization.⁶⁹ 

For a simple monovalent electrolyte, NaCl, the bottom panels of Fig. 5 show the implications of different $uSRr$ for the three aforementioned models on the longitudinal susceptibility, $χQ$⁠, and correlation–interaction term beyond DH, Σ(Q). The results show significant sensitivities of response functions of electrolyte solutions to minor differences in $uSRr$ interactions. Interestingly, over the whole concentration range, the general behaviors of $χQ$ and Σ(Q) for NNP and SD models are in good agreement with each other, while $χQ$ and Σ(Q) for the KBFF model are significantly different from those obtained from the NNP and SD models. The general behavior of Σ(Q) at high concentrations for the KBFF model at Q < 0.2 produces a positive Σ(Q) in contrast with the negative values of Σ(Q) obtained with NNP and classical SD models.

Again, for the divalent electrolyte, SrCl₂, the top panel of Fig. 6 denotes the infinite dilute PMFs extracted from both classical and qDFT molecular simulations for all of the pair interactions. A similar trend is observed with the like–like repulsion being significantly softer for qDFT. The left bottom panel of Fig. 6 shows the resulting longitudinal susceptibility functions from our mHNC calculations at concentrations ranging from 0.01 to 4.0 M. For both qDFT and classical FFs, the longitudinal susceptibility functions at low concentrations decrease as Q increases, while at high concentrations both show maxima, which appear at relatively smaller Q values for the qDFT $uSRr$⁠.

Figures 5 and 6 demonstrate clear, qualitative differences in the characterization of correlation from extracting $ΣQ$ through our mHNC procedure. These differences result in quantitative measurable differences in the collective behavior of concentrated electrolytes and will be addressed later in the paper. Indeed, it has been demonstrated that the choices in molecular interactions provide qualitative and quantitative differences in our understanding of the intrinsic properties of ions, namely local structure and single ion solvation free energies.^63,66,70,71 In this perspective we demonstrate how collective properties of concentrated electrolyte solutions can be quantitatively understood via Σ(Q) and the choice of $uSRr$⁠.

In Sec. II B 2, we developed a picture where the charge–charge correlations are characterized by Σ(Q) for various models that could not, in general, be well approximated by a power series in Q, and a Padé expansion of the form of Eq. (G2) was necessary to find near quantitative fits over the entire low-Q region studied here. This section provides a phenomenological formulation of Σ(Q) that captures the complexity of the Padé approximation. The theory involves an implicit averaging over the microscopic degrees of freedom to obtain a functional (an effective Hamiltonian) at length scales relevant for describing the collective effects of screening in electrolyte solutions.⁷² To this end, a general formulation of a Gaussian field theory coupled to an auxiliary field to demonstrate emergent phenomena in the liquid state has a long history.^53,73,74 Although providing immense theoretical insight, field theoretical formulations are often difficult to implement for general molecular systems and require further simplifications.^15,75 Below, we present a formulation that is both numerically and theoretically tractable for concentrated electrolytes, providing a molecular-scale picture of underscreening phenomena in Coulombic systems.

We follow previous studies on anisotropic next-nearest neighbor Ising (ANNNI) models that use Hubbard–Stratonovich transformations to develop a Gaussian field theory based on a single number density-based order parameter to obtain the correlation across the phase diagram representing a multicomponent complex fluid.^62,76–78 Going beyond a Gaussian field theory in the ANNNI model via diagrammatic analysis provided a route to compute corrections to mean field theory near a critical phase boundary where taking the Q → 0 limit is justified.⁶² The continuous forms of both the FCI model and the ANNNI model are similar. However, in this application it is required that we are deep in the isotropic phase well above any critical point of the model to capture electrostatic screening.¹⁶ Therefore, we seek different theoretical tools to describe the complexity of $ΣQ$ for the phenomena of screening. Motivated by previous works, we propose a simple Gaussian phenomenological model with two fields, namely, ρ^z, representing the charge density, and ψ, representing an auxiliary field for the number density.^15,53,73,74 We will demonstrate that this formulation provides a direct route to the charge–charge scattering function, S_zz(Q), and allows us to examine the precise details of u^SR(r) on the collective property of screening.

Our key insights are based on two assumptions. First, A_zz(Q), A_NN(Q), and A_Nz(Q) are assumed to capture the Q dependence of SR interactions in the kernels, noting that the $κD2Q2$ term captures the LR Coulombic interactions in the electrolyte solutions and b is a phenomenological coefficient with dimensions of length, which can originate from non-local coupling of density field gradients. Interestingly, the density and charge fields are fully decoupled in symmetric charged systems when A_Nz(Q) = 0. In the weak coupling limit for asymmetric electrolyte solutions, one may also assume A_Nz(Q) to be a constant.

We have now provided explicit expressions in Eq. (20) to directly compute $SzzQ$ from the combination of molecular simulations and estimates of $cαβSR(Q)$ through either theory or numerical approaches⁸² as outlined above and will be explicitly shown later. The salient point of Eq. (20) is that we have directly implicated the SR direct correlation function, $cαβSR(Q)$⁠, to contain all corrections beyond the DH limit. For real electrolyte solutions, a formal connection between u^SR(r) and the terms that appear in Eq. (20) remains unclear and the subject of future research. Nevertheless, below we demonstrate the veracity of these assumptions and compare the $ΣQ$ between closed expressions from analytical solutions of model electrolyte solutions directly through Eqs. (19) and (20) to the numerically extracted solutions based on Eq. (3).

Figure 7 shows that the MSA auxiliary-field results obtained from the analytical expressions of SR direct correlation functions for symmetric 1-1 electrolytes are based on Eqs. (21)–(B5). The results are in quantitative agreement with the independent MSA analytical results³⁵ [Eq. (A3)], demonstrating the accuracy of the auxiliary-field formulation in conjunction with Eq. (20) for these primitive models of electrolyte solutions.

In addition, one can use $cαβSR(Q)$ from the numerical integral equation (HNC) framework and use it as an input to Eqs. (19) and (20) to obtain the auxiliary-field formulation results based on an accurate u^SR(r), which is denoted as Σ^Aux(Q). On the other hand, one can denote the $ΣQ$ from Eq. (3) as $ΣHNCQ$⁠, where the non-Coulombic SR interactions are treated as HS. Figure 7 shows agreement between Σ^HNC(Q) and Σ^Aux(Q) for symmetric electrolyte solutions.

Importantly, the self-consistency of the formulation can be shown for 1-1 asymmetric RPMs in which the cross term A_Nz(Q) in Eq. (19) plays a crucial role in resulting in the excellent agreement of Σ^Aux(Q) and Σ^HNC(Q). This is also shown in Fig. 7. Therefore, the approximations made for deriving Eq. (19) appear to be robust. Furthermore, it is clear from Fig. 7 that, unlike symmetric RPMs, Σ(Q) in various regions of Q-space can be negative for asymmetric RPMs, which shows the significance of asymmetry in the behavior of the Σ(Q) function and the decrease of the electrostatic response in those ranges of Q-space. It is worth noting that the behavior of Σ(Q) can be due to the interplay between several factors such as asymmetry, ion size, charge valency, and details of SR interactions, and a full understanding of the negative behavior of Σ(Q) is an interesting future study.

In short, these results demonstrate that the dominant contributions to the Γ(Q) matrix in Eq. (16) can be well represented with the SR terms in Eq. (20) for both symmetric and asymmetric RPM electrolytes.

Here, we examine how Σ^Aux(Q) performs for realistic molecular descriptions of electrolytes such as NaCl and SrCl₂, whose effective SR interactions are obtained from NNP or qDFT all-atom molecular dynamics simulations that include explicit water. Using $uSRr$ constructed from NNP or qDFT, we construct $χQ$ and Σ^mHNC(Q) from the mHNC procedure and Σ^Aux(Q) for NaCl and SrCl₂ aqueous electrolyte solutions. Figure 8 demonstrates excellent agreement between Σ^mHNC(Q) and Σ^Aux(Q) for NaCl, where the b value in Eq. (16) is set to 1 Å. Although there are deviations between Σ^mHNC(Q) and Σ^Aux(Q) for SrCl₂, the shapes of both curves are in good agreement. It should be noted that different procedures have been used to extract u^SR(r) for NaCl and SrCl₂.⁶⁷ Therefore, the remarkable agreement of the NaCl mHNC with the theoretically predicted Σ^Aux(Q) could be attributed to u^SR(r) computed via inversion of the Poisson–Boltzmann equation, providing a superior reference of interaction.⁶⁷ According to the formulation presented in this work, if the exact u^SR(r) is used in Eq. (19), one would expect to get the exact numerical Σ^mHNC(Q) from Eq. (3). The sources of disagreement and the dependence on the procedure to compute u^SR(r) will be left for a future study.⁸³ 

Advances in machine learning, particularly equivariant NNPs, are leading to substantial reductions in the computational cost of computing u^SR(r) through both directly accelerating all-atom simulations and through coarse-graining.^67,68 This means that u^SR(r) will soon be determined with increasingly more accurate levels of qDFT theory⁸⁴ and for a wider range of systems, making methods of interpreting this information increasingly valuable. In addition to computing $ΣQ$⁠, u^SR(r) is also critically important for determining key thermodynamic quantities such as osmotic coefficients/activity coefficients through various routes such as the virial equation.⁶⁴ Reproducing these quantities will be an important task for future theories of electrolyte solutions.

We have now established that Σ(Q) can be obtained within the mHNC framework described earlier. In this section, we provide connections to SAXS experiments and determine the dominant poles by establishing the function form of $ΣQ$ as a ratio of rational polynomials in Q going beyond the standard TS form.⁵⁷ As we show, we find this form of $ΣQ$ supports the observation from the SAXS experiment and confirms that KTs are more general than the classic KT. The dominant poles can also be obtained in the r-space (see  Appendix F), and we show that one can get consistent poles both in Q-space and r-space. Determining the pole structure, we define generalized KTs as the decoupling of the onset of oscillations and underscreening, where the length scale Q₀ emerges at concentrations before the maximum in a₀. We show that NaCl follows classical KTs, while SrCl₂ and CaCl₂ electrolytes show generalized KTs. In general, in both types of KTs, we show that at concentrations above KTs, the screening lengths predicted from this approach are $∼1$ nm.

Having obtained Σ^mHNC(Q), for the molecular models of SrCl₂ solutions, we use Σ^fit(Q) to extract the poles by fitting to Eq. (G2) in the Q-region from 0.04 to 2.3 Å⁻¹. The poles obtained from Eq. (23) are verified by comparing them to the TS single pole fits obtained in r-space from h_zz(r) defined in  Appendix F.

Figure 9 compares the pole structure extracted from h_zz(r) and Σ(Q) when they are fitted to a damped harmonic oscillator function [Eq. (F2)] in r-space and Eq. (G2). As can be seen, the pole structure is in excellent agreement. The salient result is that the explicit coupling between the charge and number density fields is required to provide excellent fits to the scattering function S_zz(Q) over a large range of Q-space.

Figure 10 presents the pole structure of SrCl₂, which is extracted from Σ^fit(Q) for both qDFT and KBFF. They are compared with the pole structure extracted from the experimental SAXS measurements over the entire concentration range. The details about the extraction of the pole structure for the experimental SAXS measurements are presented in  Appendix F, and supplementary details are available in Ref. 14. Interestingly, the a₀ values for the qDFT over the whole concentration range are, in general, in better agreement with the SAXS a₀ values. On the other hand, the Q₀ values obtained from the KBFF seem to better represent the SAXS Q₀ values. The origin of these differences can be traced back to $uSRr$ extracted from both models. Overall, the pole structure for SrCl₂ deviates from the classic KT. In particular, the emergence of Q₀ occurs before the slope of a₀ changes sign. This demonstrates that classic KTs break for multivalent electrolytes. In addition, it is clear from Fig. 10 that experimental and theoretical screening lengths for SrCl₂ at concentrations above KTs are predicted to be less than ∼1 nm.

Figure 11 depicts the pole structures for the NNP, KBFF, and SD NaCl models when they are extracted from the fits of Eq. (F2) to the numerical h_zz(r) from the mHNC method. Note that Σ(Q) for SD and NNP models are generally in good agreement, which translates to their root structures. Interestingly, for NaCl, the classic KT is established for all forms of interactions studied. It is worth mentioning that finding the dominant poles from Σ(Q) is numerically challenging due to the existence of multiple poles predicted in previous studies.^47,48,87 Although the numerical procedure used for extracting the dominant poles for SrCl₂ from Σ(Q) is consistent with those obtained in r-space [see Eq. (F2)], our preliminary results suggest that using a similar numerical approach for NaCl does not provide satisfactory agreement with the r-space approach. Therefore, a full analysis of pole structure for NaCl will be the subject of further examination in conjunction with the analysis of new SAXS experimental data. Nevertheless, our theoretical predictions based on the dominant pole analyses show that the screening lengths are less than $∼1$ nm at concentrations above KTs.

We have recently established the connection between the SAXS experimental measurements and charge–charge correlations in electrolyte solutions.¹⁴ Although previous theories of electrolyte solutions^8,13,20,37 provide a qualitative picture of the collective effects of screening, in their present form, they cannot predict the complex pole structure reported in Fig. 10 verified by SAXS measurements. In this perspective, we present a Gaussian field theory for $ΣQ$ to explain these recent SAXS measurements.¹⁴ The theory involves coupling between the number and charge density and provides simple formulas that can be utilized to study the complexities of molecular representations of concentrated electrolytes. In addition, we have demonstrated that $ΣQ$ can be used to unify the theories of electrolytes (see Fig. 3) and provide a conceptual understanding of the origins of charge-charge correlations in concentrated electrolytes.

To this end, our formulation provides a direct and simple route to $SzzQ$ [Eq. (18)] and provides a general formula $ΣQ$ [Eq. (19)] that can be cast in terms of SR interactions through SR direct correlation functions, which can be estimated by liquid state theory [Eq. (20)]. An exact $ΣQ$ can be extracted from simulation and, therefore, constitutes a theory for concentrated electrolytes. According to the formulation presented in this work, if the exact SR interactions are used in Eq. (19), one would expect to get the exact numerical Σ^mHNC(Q) from Eq. (3). We have shown in this work that $ΣQ$ can be well represented through a molecular-informed $uSRr$ within the liquid state framework. This demonstrates the importance of the coupling between dielectric theory and PMF that properly integrates over solvent fluctuations to obtain a molecular picture of the collective properties of electrolyte solutions.^63,88 The choice of $uSRr$ can influence the screening length scales over the entire concentration range, as shown in pole structure and, therefore, the details of the KT shown in Figs. 10 and 11. The agreement between the Gaussian (auxiliary) field theory presented in Sec. III with the mHNC results in Fig. 8 is remarkable. To the extent that this agreement relies on a more rigorous determination of $uSRr$ via the direct subtraction of the Coulomb term is an open question. Indeed, Ref. 67 provides a more sophisticated route to $uSRr$ and appears to give substantially better agreement (see Fig. 8). Future studies will focus on the self-consistent procedures such as those suggested by Outhwaite et al.^{9,37,44,45,50} to construct more rigorous representations of $uSRr$⁠.

The importance of treating the number density explicitly can be discerned through the form of $ΣQ$ obtained from treating only the charge density field, namely the FCI. Although the FCI can describe the qualitative physics of underscreening, it cannot reproduce the function form needed to reproduce molecular simulation and experimental data over a broad Q range. Therefore, the number density and charge density fields and their coupling are essential for the advanced theories of electrolyte solutions. With this in mind, we discuss the sensitivities of $ΣQ$ to the symmetry and asymmetry in $uSRr$ and show the significance of the coupling of number density and charge density in theory. We observe that $ΣQ$ for the studied symmetric RPM models in the low Q-region (Q < 2 Å⁻¹) remains positive, but $ΣQ$ can become negative for the asymmetric RPMs (see Fig. 7). This is further confirmed for the asymmetric classical and ab initio-based models of 1-1 and 1-2 electrolytes (see Figs. 5 and 6). It would be an interesting future study to fully understand the complex behavior of Σ(Q) in real electrolyte solutions and its role in predictions of the thermodynamic properties such as activity and osmotic coefficients.⁸⁹ 

We also demonstrate that Σ(Q) can be used to determine the screening lengths and KTs. We show that KTs can go beyond classic KTs in which the emergence of the order in LR charge–charge correlations, Q₀, is coupled with known changes in the behavior of underscreening (maximum in a₀). In this regard, we generalize KTs for multivalent electrolytes while our numerical analyses show minimal deviations from the classic KTs for monovalent electrolytes. Moreover, although the TS form seems to be robust as a single pole description of underscreening and the breakdown of DH, this work provides a formulation that can differentiate between different forms of molecular interaction and, therefore, brings us one step closer to having a molecular rationalization for underscreening in electrolytes. Finally, above the KTs, our dominant pole analyses of bulk electrolytes reveal screening lengths comparable to the size of ion hydration spheres. This is in contrast with the unusually large screening lengths reported in recent surface technique studies.^25,26,30

The authors would like to acknowledge J. Chun, Y. Levin, D. T. Limmer, D. V. Matyushov, and C. W. Outhwaite for their useful discussions. In particular, the authors acknowledge D. T. Limmer for insightful discussions and sharing his perspective on the two-field formulation. The authors were supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Chemical Sciences, Geosciences, and Biosciences Division, Condensed Phase and Interfacial Molecular Science programs Grant Nos. FWP 16249 and FWP 16248. Pacific Northwest National Laboratory is operated by Battelle for the U.S. DOE under Contract No. DE-AC05-76RL01830. Computing resources were generously allocated by PNNL’s Institutional Computing program. This research also used resources of the National Energy Research Scientific Computing Center (NERSC), a Department of Energy Office of Science User Facility using NERSC Award No. BES-ERCAP-26557. We acknowledge Paula J. Seeger for a thorough copy edit of the paper.

The authors have no conflicts to disclose.

Mohammadhasan Dinpajooh: Conceptualization (lead); Data curation (equal); Formal analysis (lead); Investigation (lead); Methodology (equal); Software (equal); Validation (lead); Visualization (equal); Writing – original draft (lead); Writing – review & editing (lead). Nadia N. Intan: Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal). Timothy T. Duignan: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal). Elisa Biasin: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Resources (equal); Visualization (equal); Writing – review & editing (equal). John L. Fulton: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Resources (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal). Shawn M. Kathmann: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal). Gregory K. Schenter: Conceptualization (lead); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (lead); Methodology (equal); Project administration (lead); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal). Christopher J. Mundy: Conceptualization (lead); Data curation (equal); Formal analysis (equal); Funding acquisition (lead); Investigation (lead); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (lead); Writing – review & editing (lead).

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

In this section, the $ΣQ$ is presented for 1-1 electrolyte solutions for other common theories/models.

Independent of the above-mentioned expressions, one can get the $ΣQ$ term in the MSA as presented in Eq. (A3).

In this section, we provide the details of numerical methods used in this perspective.

Ab initio Molecular Dynamics (AIMD) simulations for the SrCl₂ electrolyte solution in the infinite dilution limit were performed using the Quickstep module of the CP2K software.⁹⁰ The valence electrons were treated explicitly at the DFT level employing the revised version of the Perdew–Burke–Ernzerhof (rev-PBE) functional⁹¹ and the double-ζ MOLOPT basis set (DZVP–MOLOPT–SR–GTH)⁹² with a density cutoff of 400 Ha. The core electrons on all-atoms were represented by Goedecker–Teter–Hutter (GTH) pseudopotentials.⁹³ The Grimme dispersion correction (DFT-D3)⁹⁴ was employed to account for LR dispersion interactions. The canonical ensemble (NVT ensemble) was employed with the Nos$é$–Hoover thermostat, which was used to maintain the average temperature at 300 K. The initial configurations were generated by increasing the distance between ion pairs at a rate of 10⁻³ Å/fs to separate the ion pairs, and 50 windows with a separation of 0.1 Å were generated. For each window, a restraint of 80 kcal/Ų was applied to keep the ions separated within a fixed distance. A time step of 0.5 fs was used to generate a 20 ps trajectory, of which the first 2.5 ps of simulation time was dedicated to system equilibration and the rest for data collection.

The NaCl short range PMFs were originally presented in Ref. 95. In brief, short AIMD simulations of 2.5M NaCl in solution were run with CP2K using GPW and the strongly constrained and appropriately normed (SCAN) functional.⁹⁶ An equivariant neural network potential (NNP) was then trained to reproduce the forces and energies of this simulation using the NequIP software package.⁹⁷ Since the NNP only has access to local information (SR, $<5$ Å) the LR electrostatic ion–ion interactions were removed from the forces and energies prior to training using a dielectrically screened Coulomb interaction.^98,99 They were then added back in during all NNP simulations. The NNP enables much longer and larger simulations than standard AIMD, allowing for the ion–ion RDFs to then be converged. The modified Poisson–Boltzmann (PB) equation was then inverted to find the best SR PMF force that reproduces the RDF at the given concentration.⁶⁷ 

Here, a brief description of classical all-atom MD simulations is presented. The simulations were performed using the LAMMPS software program¹⁰⁰ at 298 K and a pressure of 1 atm. The SPC/E water model was used, and for ions, the Smith–Dang (SD) and Kirkwood–Buff force fields (KBFFs) were used¹⁰¹ and presented in Table I. The periodic boundary conditions were applied in 3-dimensions and simulations were performed in the isothermal-isobaric ensemble (NPT) using the Nos$é$–Hoover thermostat and barostat. The standard velocity-Verlet time integrator was used with a time step of 2 fs, and the SHAKE procedure was used to conserve the intramolecular constraints. The Lennard-Jones (LJ) and real-space part of the Coulombic interactions were truncated at 10 Å with additional switching/shifting functions that ramped the energy smoothly to zero between 10 and 14 Å. The conducting metal (tinfoil) boundary conditions were used to treat the electrostatic interactions, and the particle–particle particle–mesh solver was used.

For the potential of mean force calculations in the infinite dilution limit, the harmonic biasing/umbrella sampling method was used to enforce moving restraints, where the restraints were introduced on the collective variable ξ, which is the distance vector between the ion centers. A force constant of 100 kcal/mol/Ų was used, where the centers of the harmonic restraints were changed during the simulations (usually from faraway distances to closer distances) in discrete stages, where each stage consisted of 30 000 time steps. The initial configuration was equilibrated for 10⁵ MD steps, and a grid spacing of 0.1 Å was used to discretize the distance between the ions. The thermodynamics integration methods available in the colvars package¹⁰⁰ were used to get the free energies, which included the Jacobian term 2k_BT  ln(r) to the output free energies.

The OZ equations were then solved using the hyper-netted chain (HNC) closure. This numerical procedure for real electrolyte solutions is called the molecular-informed HNC (mHNC) method, which builds upon the previous studies by Rasaiah,⁵¹ Vrbka et al.,⁶⁴ and others.⁵ We used a modified version of the SunlightHNC code (PNNL-SunlightHNC),⁸² where the inputs to the PNNL-SunlightHNC code were $uSRr$ between ions, ion charges, and ion concentrations that were described in Sec. II C.

In a previous work, we demonstrated that the mHNC results are in good agreement with all-atom MD simulation results for both SR and LR correlations (see Figs. S5 and 4 in Ref. 14) of ErCl₃ electrolytes at various concentrations. In this section, we provide numerical results showing that when the mHNC method is applied to simple RPMs, the SR and LR correlations are, not surprisingly, in good agreement with those obtained from molecular simulations. The notation mHNC is used here because the all-atom molecular-informed SR interaction is utilized.

The RPMs are modeled using a 1 − tanh[(r − A)/0.05] function with A = 4.25 Å in the molecular simulations. The HS simulations were performed using the CP2K software program⁹⁰ for 60 ns. The simulation box length was set to 66.4644 Å, consisting of 1061 positive and 1061 negative HSs. A time step of 1 fs was used, and the dielectric constant of water was set to 80. The periodic boundary conditions were applied in 3-dimensions and the MD simulations were performed for 60 ns in the canonical ensemble (NVT) using the Nos$é$–Hoover thermostat and the standard velocity-Verlet time integrator. The pair correlations were obtained by analyzing 60 000 MD frames. In parallel, making use of the HS interactions, the mHNC calculations were performed under the same conditions as the explicit MD of HS with a diameter of 4.25 Å.

Figure 12(a) demonstrates good agreement between all-atom MD and mHNC results for SR pair correlation functions, while Fig. 12(b) shows that the all-atom MD and mHNC results are in good agreement up to r ∼ 20 Å, noting that the LR behaviors of pair correlations cannot be captured by direct simulation. We argue that for the purpose of this perspective, the accuracy of the mHNC method is plausible.

In this section, the longitudinal susceptibility functions, $χQ$⁠, of the symmetric 1-1 RPM electrolyte solutions at 4.5 M for various ion sizes are reported in Fig. 13(a) and compared with the DH limit with the DH S_zz(Q) given by $Q2/(κD2+Q2)$⁠. As can be seen, as the ion size increases, a peak in $χQ$ starts to appear around 0.5–0.6 Å⁻¹ for ion sizes with a > 3 Å. This indicates the significance of local packing on the behavior of $χQ$⁠. Interestingly, from 2.3 to 3.3 Å⁻¹, where Σ(Q) starts to become negative for an ion size with a = 6 Å (see Fig. 4 in the main text), the slope of $χQ$ changes significantly as illustrated in the inset in Fig. 13(a).

Here, we also show our numerical results for CaCl₂ electrolyte solutions. The top panel of Fig. 14 shows $uSRr$ for the CaCl₂ electrolyte solution, while the bottom panels show the longitudinal dielectric susceptibility, χ(Q), and $ΣQ$⁠. Similar to SrCl₂, the KBFF model shows that peaks in $χQ$ (at higher concentrations) appear at larger Q values when compared to the qDFT model. In addition, the KBFF model shows negative behavior of Σ(Q) in Q < 0.2 Å⁻¹, while in the qDFT model, $χQ$ remains positive.

Figure 15 shows the root structures of CaCl₂ electrolyte solution when they are extracted from the charge–charge correlation functions. While the pole structure for both models has similar features, the Q₀ inverse lengths for the qDFT model turn out to be smaller than the ones in the KBFF model at higher concentrations. The screening lengths, on the other hand, are larger for the qDFT model than the ones obtained from the KBFF model.

Finally, assuming that S_zz(Q) has a form of a TS function presented in Eq. (F1), a similar recipe has been used to fit the numerical mHNC Σ(Q) results [obtained from $χQ$ per definitions in Eq. (3)] in Fig. 17.

A careful reader may ask how auxiliary-field formalism can improve our understanding of the effects of SR interactions on charge–charge correlations. In fact, this is a deep question whose examinations require rigorous and full analyses. In this work, preliminary analyses are presented in Sec. IV. In this section, we reveal a related aspect to the correlations obtained from the SAXS measurements. We present the distributions resulting from the auxiliary formulation and compare them with the traditional TS distributions used in the analysis of SAXS spectra. In previous studies, we discussed the significance of TS distributions in the analysis of SAXS spectra and established that the pole structure of multivalent electrolytes for SAXS signals and charge-charge correlations are consistent.¹⁴ However, the poles from SAXS were obtained from the fits of SAXS signals to the TS distributions, which are very challenging due to the sensitivity of the analyses to the Q range chosen for fitting. In  Appendix F, we briefly present the details of the analyses. Here, we use numerical mHNC results to demonstrate that auxiliary field formulation can provide a unique avenue to obtain distributions that have implications for the analyses of SAXS spectra, which is the subject of future study.