Unified microscopic theory of equilibrium thermodynamics and ion association in aqueous and non-aqueous electrolytes with explicit hard-core size

Within the framework of a functional integral formalism incorporating ionic charge and hard-core (HC) interactions on an equal footing, we formulate a unified theory of equilibrium thermodynamics and ion association in charged solutions. Via comparison with recent Monte-Carlo (MC) simulation results (J. Forsman et al., PCCP 26, 19921 (2024)), it is shown that our approach is able to predict with quantitative precision the pair distributions of monovalent ions with the typical hydrated sizes d = 3.0 A and 4.0 A up to the molar concentration ni = 2.0 M. Moreover, comparison with additional simulation data from the literature indicates that within the characteristic regime of ionic packing fraction eta <0.1, the formalism can accurately account for the ion size dependence of the excess energy and pressure from d = 14.3 A down to d = 1.6 A. Via the adjustment of the hydration radius, the theory can also reproduce the non-monotonic salt dependence of the experimentally measured osmotic coefficients of various aqueous and non-aqueous solutions. In accordance with AFM experiments involving non-aqueous electrolytes, the underlying sharp competition between the opposite charge attraction and the excluded volume constraint is shown to limit the occurrence of substantial ionic pair formation to the submolar concentration regime ni <50 mM; at larger concentrations, HC repulsion hinders ion association and results in the quasi-saturation of the pair fraction curves.

💡 Research Summary

The paper presents a unified microscopic theory that simultaneously describes equilibrium thermodynamics and ion association in both aqueous and non‑aqueous electrolytes, explicitly incorporating ionic charge and hard‑core (HC) repulsion on an equal footing. Starting from a minimal model of p ion species, each represented as a charged sphere of diameter d, the authors formulate the grand‑canonical partition function using a functional integral representation. The total interaction is split into Coulombic and HC contributions, with the Coulomb potential further decomposed into short‑range (v_s) and long‑range (v_l) parts via a variational splitting length σ. This splitting is motivated by the need to treat long‑range electrostatics within a weak‑coupling Gaussian approximation while handling short‑range, strongly coupled interactions (including HC exclusion) with a virial‑type treatment.

A key technical advance is the derivation of self‑consistent Schwinger‑Dyson (SD) identities that relate the two‑point correlation functions G_γ(r) of the fluctuating fields (γ = s, l, h) to the physical parameters of the electrolyte. By imposing the variational condition ∂Ω/∂σ = 0 on the grand potential Ω, the optimal σ is obtained, ensuring that the splitting adapts to the specific ion size, concentration, and dielectric environment. The authors also provide an analytical expression for the Gaussian‑level correlation function, dramatically reducing the number of Fourier transforms required and accelerating numerical implementation.

The resulting Self‑Consistent Debye‑Hückel (SCDH) formalism—an upgraded version of the previously introduced Cumulative‑Corrected DH (CCDH) approach—captures both the excess internal energy and pressure with high accuracy across a broad parameter space. Validation against extensive Monte‑Carlo (MC) data, including recent simulations by Forsman et al. (PCCP 2024), shows that the theory reproduces radial distribution functions g_ij(r) for monovalent ions with hydrated radii d = 3.0 Å and 4.0 Å up to molar concentrations n_i ≈ 2 M, with errors typically below 5 %. Moreover, in the regime of low packing fraction (η ≲ 0.1), the theory accurately predicts the ion‑size dependence of excess thermodynamic quantities for diameters ranging from 1.6 Å to 14.3 Å.

Ion association is addressed through an explicit pair‑fraction model derived from the partition function. The analysis reveals a sharp competition between opposite‑charge attraction and HC exclusion: substantial ion pairing occurs only below ~50 mM, where electrostatic attraction dominates; at higher concentrations, the excluded‑volume constraint suppresses further pairing, leading to a quasi‑saturation of the pair‑fraction curves. This behavior aligns with atomic‑force microscopy observations in non‑aqueous electrolytes.

The authors also demonstrate that by adjusting the effective hydration radius, the SCDH framework reproduces the experimentally observed non‑monotonic dependence of osmotic coefficients on salt concentration for a variety of solvents. This highlights the theory’s capability to incorporate solvent dielectric properties (ε_s) and specific ion solvation effects without resorting to ad‑hoc parameters.

Computationally, the upgraded SCDH scheme reduces the number of required Fourier transforms by analytically evaluating the Gaussian correlation function, achieving speed‑ups of an order of magnitude compared with earlier implementations. This efficiency enables systematic scans over ion size, concentration, and solvent permittivity, facilitating direct comparison with experimental data.

Limitations are acknowledged: the variational σ is derived within a mean‑field Gaussian framework, so at very high packing fractions (η > 0.1) or for multivalent ions, higher‑order correlations may become important. Nonetheless, within its validated domain, the theory bridges the gap between simple Debye‑Hückel models and computationally intensive simulations, offering a versatile, analytically tractable tool for predicting thermodynamic and structural properties of concentrated electrolytes.