Quadrupolar ²³Na⁺ NMR relaxation as a probe of subpicosecond collective dynamics in aqueous electrolyte solutions

The proper characterization and modeling of the solvation structure of alkaline cations (e.g., Li⁺, Na⁺, and K⁺) in aqueous solution are of considerable interest both in physiological systems^{1,2,3,4,5,6,7,8} and electrolytes used for electrochemical devices^{9,10,11,12,13}. Nuclear magnetic resonance (NMR) spectroscopy provides an excellent source of dynamic and structural information for a number of nuclear species, including ²³Na with a nuclear spin 3/2 and close to 100% natural abundance that produces the second strongest NMR signal after protons in biological tissues¹⁴. The NMR sensitivity of sodium is 9.2% of that of proton, while its typical concentration can be three, or more, orders of magnitude lower. Thus, in biological systems the sodium signal-to-noise ratio is 3000–12,000 times lower than that of ¹H⁴. Nonetheless, the longitudinal relaxation time T₁ of ²³Na (typically 40 ms and below) is short compared to that of ¹H (on the order of seconds)⁴, allowing for rapid averaging of the signals such that quantitative analysis is made possible within reasonable time scales¹⁵.

The shortness of the ²³Na NMR relaxation times is due to a fluctuating quadrupole interaction related to the changes in the solvation shell and the proximity of other ions¹⁶. The relaxation rate is determined from a combination of the strength of the electric field gradient (EFG) at the nucleus quantified by means of the quadrupolar coupling constant (QCC) C_Q and the characteristic correlation time τ_c with which the memory of fluctuations is lost. While the knowledge of C_Q and τ_c can potentially provide information about the hydration sphere structure^17,18 and useful dynamic properties (e.g., diffusion coefficients, viscosity, or conductivity^19,20,21), respectively, their unambiguous determination from the experimentally-measured rates in solution has remained essentially impossible²².

Different models have been suggested to rationalize quadrupolar relaxation using dielectric description^23,24, mode-coupling analysis²⁵, definite molecular processes (e.g., water reorientation^26,27,28 and collective symmetry-breaking fluctuations²⁹), or Brownian rotational diffusion^19,20,21. Ab initio^{30,31,32,33,34} and classical^{29,35,36,37,38,39,40,41,42} molecular dynamics (MD) simulations have been indispensable in assessing predictions of such theories, and invalidated the isotropic monoexponential character of the quadrupolar relaxation that is often assumed under a continuous-solvent description. A pronounced role of intermolecular cross-correlations on the relaxation was emphasized^17,29,38. While classical MD often relies on the Sternheimer approximation to incorporate electron cloud effects^39,42,43,44, ab initio methods provide the best accuracy of the computed EFG at the ion position^{30,32,45,46,47}. However, the associated high computational cost often impedes the long-time sampling of EFG fluctuations^30,31 and the accuracy of the correlation time estimates, even for aqueous ions at infinite dilution^30,32. Hence, uncertainties arising both in ab initio and classical MD-based approaches hindered the quantitative comparison with experimental NMR relaxation rates and systematic analysis of the quadrupolar relaxation mechanisms.

Here we show that applying ab initio calculations to parametrize C_Q in conjunction with classical MD to evaluate τ_c allows reaching a good agreement between the calculated and experimentally-obtained quadrupolar rates of ²³Na⁺ in electrolyte solutions over a broad range of salt concentrations and temperatures, thereby enabling a systematic analysis of the relaxation pathways and detailed modeling of the underlying dynamics. We find that the main effect of increased relaxivity is due to a lengthening of the correlation times, rather than a change of the average quadrupolar coupling constant. Counterintuitively, the latter varies mildly over the range of considered parameters, slightly decreasing with concentration and increasing with temperature. We conclude that, contrary to the commonly-assumed picture, rotational models based on the water dipole reorientation or Stokes-Einstein-Debye relation significantly overestimate the EFG correlation times. Rather, our results indicate that the EFG relaxation is mainly determined by the dynamics in the first two solvation shells around the solute and occurs over a time scale comparable to that of solution structural rearrangements. This work thus suggests that the subpicosecond collective dynamics of the liquid primarily drive the quadrupolar relaxation at the sodium-ion position, thereby offering insights into the quadrupolar relaxation mechanisms in electrolyte solutions.

Electron cloud contribution to electric field gradients

We perform density functional theory (DFT) calculations to determine quantum EFGs at the Na⁺ position in aqueous sodium chloride (NaCl) solutions at varying salt concentrations c = 1–5 molal (denoted with mol kg⁻¹ or m) at T = 25 °C (see Methods for details). The projector-augmented wave (PAW) method^45,47,48 is used to reconstruct the all-electron charge density at the nucleus. A configuration of a NaCl solution at 5 mol kg⁻¹ with converged charge densities is highlighted in Fig. 1a.

a Representative system configuration of a NaCl solution at 5 m. The gray opaque regions around water molecules show charge densities obtained with DFT PAW calculations (see Methods). b Component-wise comparison of ab initio EFGs, \({V}_{\alpha \beta }^{{{{{{{{\rm{AI}}}}}}}}}\), against classical EFGs, \({V}_{\alpha \beta }^{{{{{{{{\rm{ext}}}}}}}}}\), at the position of Na⁺ ions on the same set of configurations for different salt concentrations c. The solid lines indicate the fit for an effective Sternheimer factor γ_eff: \({V}_{\alpha \beta }^{{{{{{{{\rm{AI}}}}}}}}}=(1+{\gamma }_{{{{{{{{\rm{eff}}}}}}}}}){V}_{\alpha \beta }^{{{{{{{{\rm{ext}}}}}}}}}\). c Effective Sternheimer factors for Na⁺ obtained from the linear fit (yellow squares) or from the ratio \({(1+{\gamma }_{{{{{{{{\rm{eff}}}}}}}}}^{{\prime} })}^{2}=\langle {{{{{{{{\bf{V}}}}}}}}}_{{{{{{{{\rm{AI}}}}}}}}}^{2}\rangle /\langle {{{{{{{{\bf{V}}}}}}}}}_{{{{{{{{\rm{ext}}}}}}}}}^{2}\rangle\) (blue circles) at different c. d EFG variance at the Na⁺ position for different c as obtained directly with ab initio calculations (green triangles), or using the value of γ_eff (yellow squares) and \({\gamma }_{{{{{{{{\rm{eff}}}}}}}}}^{{\prime} }\) (blue circles) at infinite dilution. The error bars in c, d were calculated using bootstrapping.

Full size image

In classical MD, the electron cloud contribution to the EFG can be incorporated by means of the Sternheimer approximation^43,44, in which the full EFG at the nucleus V is proportional to the EFG created by the external charge distribution V_ext: V ≃ (1 + γ)V_ext. Here, the electron cloud polarization is included via the simple rescaling factor 1 + γ, with the Sternheimer (anti)shielding factor γ being typically large γ ≫ 1⁴⁴. To validate the Sternheimer approximation for present systems, we have compared ab initio, V_AI, against classical, V_ext, EFGs at the Na⁺ position, as determined on the same set of classically generated solution configurations (see Methods). Consistently with aqueous ions at infinite dilution^39,42, we find a strong correlation between V_AI and V_ext for all c = 1–5 mol kg⁻¹, as seen in Fig. 1b for the two extreme cases. The latter allows us to define effective Sternheimer factors γ_eff through the linear fit \({V}_{\alpha \beta }^{{{{{{{{\rm{AI}}}}}}}}}=(1+{\gamma }_{{{{{{{{\rm{eff}}}}}}}}}){V}_{\alpha \beta }^{{{{{{{{\rm{ext}}}}}}}}}\). As seen in Fig. 1c, the resulting γ_eff feature a small increase with c (<5% compared to the infinite dilution value γ_eff,0 = 10.54 ± 0.11⁴²) associated with the modifications of the ion’s solvation sphere (see Supplementary Note 5 and 6).

Despite the small changes of γ_eff with increasing c, the Sternheimer approximation for the EFG variance, \({(1+{\gamma }_{{{{{{{{\rm{eff}}}}}}}}})}^{2}\langle {{{{{{{{\bf{V}}}}}}}}}_{{{{{{{{\rm{ext}}}}}}}}}^{2}\rangle\), that is necessary for the NMR relaxation rate computation (Eq. (1) in Methods) underestimates the ab initio value \(\langle {{{{{{{{\bf{V}}}}}}}}}_{{{{{{{{\rm{AI}}}}}}}}}^{2}\rangle\) by more than 20% (highlighted in Fig. 1d using γ_eff,0). This again underlines the deficiencies of the Sternheimer approximation⁴² that does not take into account non-electrostatic electron cloud polarization effects, such as short-range repulsion^18,49,50. To improve upon the variance predictions, we formally define the Sternheimer factor \({\gamma }_{{{{{{{{\rm{eff}}}}}}}}}^{{\prime} }\) as \({\left(1+{\gamma }_{{{{{{{{\rm{eff}}}}}}}}}^{{\prime} }\right)}^{2}=\langle {{{{{{{{\bf{V}}}}}}}}}_{{{{{{{{\rm{AI}}}}}}}}}^{2}\rangle /\langle {{{{{{{{\bf{V}}}}}}}}}_{{{{{{{{\rm{ext}}}}}}}}}^{2}\rangle\) with state-dependent values of \(\langle {{{{{{{{\bf{V}}}}}}}}}_{{{{{{{{\rm{AI}}}}}}}}}^{2}\rangle\) and \(\langle {{{{{{{{\bf{V}}}}}}}}}_{{{{{{{{\rm{ext}}}}}}}}}^{2}\rangle\). Similarly, \({\gamma }_{{{{{{{{\rm{eff}}}}}}}}}^{{\prime} }\) slowly grows with c, yet starting from a markedly enhanced value of \({\gamma }_{{{{{{{{\rm{eff,0}}}}}}}}}^{{\prime} }=12.09\pm 0.14\) at infinite dilution (Fig. 1c). The EFG variance prediction \({(1+{\gamma }_{{{{{{{{\rm{eff,0}}}}}}}}}^{{\prime} })}^{2}\langle {{{{{{{{\bf{V}}}}}}}}}_{{{{{{{{\rm{ext}}}}}}}}}^{2}\rangle\) using \({\gamma }_{{{{{{{{\rm{eff,0}}}}}}}}}^{{\prime} }\) at infinite dilution is within 5% accuracy of \(\langle {{{{{{{{\bf{V}}}}}}}}}_{{{{{{{{\rm{AI}}}}}}}}}^{2}\rangle\) within the considered concentration range, a much better estimate in comparison to the simple Sternheimer approximation (Fig. 1d). While not capturing all condensed-phase effects that arise with increasing c, the estimate \({(1+{\gamma }_{{{{{{{{\rm{eff,0}}}}}}}}}^{{\prime} })}^{2}\langle {{{{{{{{\bf{V}}}}}}}}}_{{{{{{{{\rm{ext}}}}}}}}}^{2}\rangle\) provides a fair accuracy, reproduces the trend of \(\langle {{{{{{{{\bf{V}}}}}}}}}_{{{{{{{{\rm{AI}}}}}}}}}^{2}\rangle\) to decrease with the salt concentration (see Fig. 1d), and permits avoiding computationally expensive DFT calculations at multiple system state points of interest. As discussed below, in combination with the EFG relaxation dynamics captured at the classical level, this approach provides a good description of the quadrupolar ²³Na⁺ NMR relaxation rates in aqueous solutions.

Relaxation of electric field gradient fluctuations

We perform classical MD simulations employing the Madrid-2019 force field (FF)^51,52 to facilitate the long-time sampling of EFG fluctuations and to investigate the mechanisms behind the concentration and temperature behavior of the quadrupolar ²³Na⁺ NMR relaxation rate in aqueous sodium chloride, bromide (NaBr), and fluoride (NaF) solutions (see Methods for simulation details). Two facts give confidence in this approach: (i) a very strong correlation between the full and classical (external) EFGs (Fig. 1b), indicating that the dynamics of the former should be largely determined by that of the latter; (ii) while classical MD with rigid water molecules do not quantitatively reproduce the librational or hydrogen-bond stretching water dynamics that occur at very short times below ~ 50 fs⁵³, it is expected that these high-frequency motions do not significantly affect the dominating long-time (~1 ps) EFG relaxation mode (e.g., see ref. 40 and below).

We provide a systematic description of the ²³Na⁺ EFG relaxation at short, intermediate, and long times ranging from a few fs to tens of ps. Increasing salt concentration c or decreasing temperature T causes a profound slow-down of the EFG fluctuations at the ion position (Fig. 2). Due to a qualitative similarity of the EFG relaxation in the solutions considered, here we will focus on the case of NaCl; see Supplementary Information (SI) for NaBr and NaF. Figure 2a shows the autocorrelation functions (ACFs) of the classical EFG at the Na⁺ position, \({C}_{{{{{{{{\rm{EFG}}}}}}}}}(t)\equiv \left\langle {{{{{{{{\bf{V}}}}}}}}}_{{{{{{{{\rm{ext}}}}}}}}}(0):{{{{{{{{\bf{V}}}}}}}}}_{{{{{{{{\rm{ext}}}}}}}}}(t)\right\rangle\), as a function of c at T = 25 °C (see Supplementary Fig. 9 for other T). Similarly to a single Na⁺ in water^32,38,39,42, C_EFG(t) relaxes in two steps: (i) a rapid initial decay happening at t ≲ 0.2 ps that corresponds to ≈70% of the EFG decorrelation. This is in good agreement with the ab initio MD results for Na⁺ at infinite dilution³², highlighting the validity of the classical approach; (ii) a much slower secondary decay occurring in the picosecond regime. As seen in Fig. 2a, the increase in c leaves the initial fast decay practically unchanged, while causing a pronounced slow-down of the second decay mode. The latter is highlighted in the inset of Fig. 2a showing the EFG ACFs for t < 1 ps for different c at T = 25 °C (see also Supplementary Fig. 10). A qualitatively similar trend is found with decreasing temperature, as we show in Fig. 2b at c = 4 m and in Supplementary Fig. 9 for other c.

a Normalized autocorrelation functions C_EFG(t)/C_EFG(0) of the EFG at the position of a Na⁺ ion obtained using classical MD simulations for different salt concentrations c at T = 25 °C in aqueous NaCl solutions (c increases from bottom to top). Qualitatively similar trends are found for other concentrations and temperatures (see Supplementary Fig. 9). Insets in a highlight the short-time behavior of the ACFs for t < 1 ps and a typical system configuration at c = 4 m (Na⁺ and Cl⁻ ions are blue and green, respectively). b Temperature behavior of C_EFG(t) at c = 4 m (T decreases from bottom to top). c Long-time behavior of C_EFG(t) plotted on a double logarithmic scale at c = 4 m for different temperatures (the legend shown in b). The black solid and dashed lines highlight a power-law scaling ~t^α with α = −5/2 and α = −3/2, respectively. See also Supplementary Fig. 15 for C_EFG(t) multiplied by t^5/2 and t^3/2. Shaded regions in a–c indicate standard errors from multiple independent simulation runs.

Full size image

The form of the EFG ACF decay in Fig. 2 suggests a collective pathway behind the relaxation. After the initial fast decay that can be described with an exponential \(\sim {e}^{-t/{\tau }_{{{{{{{{\rm{f}}}}}}}}}}\) with τ_f ≈ 62 fs, we find a development of a much slower relaxation mode that profoundly depends on c and T. Compared to earlier results⁴⁰, our long-time sampling reveals that the slow part of the EFG ACF is not exponential, as clearly seen from the behavior of C_EFG(t) on a semi-logarithmic scale in Fig. 2a–b and as we show with explicit fits in Supplementary Note 9. Except at very long times, we find that the slow decay can be modeled either with a two-exponential (Supplementary Fig. 11) or a stretched exponential fit \(\sim {e}^{-{(t/{\tau }_{{{{{{{{\rm{s}}}}}}}}})}^{\beta }}\) with β = 0.67 ± 0.05 (Supplementary Figs. 12 and 13), which suggests a broad distribution of contributing relaxation modes (Supplementary Fig. 14). Although observed over a limited time range (up to a decade), we find that the long-time tail of the EFG ACFs is consistent with a power-law ~t^−5/2, as shown with C_EFG(t) on a log-log scale for c = 4 m in Fig. 2c and with t^5/2C_EFG(t) in Supplementary Fig. 15. Such a hydrodynamic tail was predicted by a mode-coupling theory of Bosse et al. for the EFG ACF in molten salts²⁵. It originates from the coupling between the ion motion and shear excitations in the liquid, a mechanism causing the well-known ~t^−3/2 tail of the velocity ACF⁵⁴. While sampling of the EFG fluctuations at even longer time scales is necessary to decisively confirm to presence of ~t^−5/2 regime, our results for Na⁺ in Fig. 2c suggest that its relative contribution may be marginal because the apparent onset of the algebraic decay occurs at times when the ACF has decayed considerably.

Quadrupolar relaxation rates

The combination of EFG fluctuations captured at the classical level and consistent inclusion of the electron cloud contribution to the EFG enables reaching a good quantitative agreement between the calculated and experimentally-measured quadrupolar NMR relaxation rate for ²³Na⁺ in aqueous NaCl, as we compare in Fig. 3 with filled and open symbols, respectively. As seen in Eq. (1) in Methods, the quadrupolar relaxation rate is proportional to the product of the effective correlation time of EFG fluctuations, \({\tau }_{{{{{{{{\rm{c}}}}}}}}}={C}_{{{{{{{{\rm{EFG}}}}}}}}}^{-1}(0)\int\nolimits_{0}^{\infty }{{{{{{{\rm{d}}}}}}}}t\,{C}_{{{{{{{{\rm{EFG}}}}}}}}}(t)\), and the EFG variance, which we approximate as \(\langle {{{{{{{{\bf{V}}}}}}}}}^{2}\rangle={(1+{\gamma }_{{{{{{{{\rm{eff,0}}}}}}}}}^{{\prime} })}^{2}\langle {{{{{{{{\bf{V}}}}}}}}}_{{{{{{{{\rm{ext}}}}}}}}}^{2}\rangle\) with \({\gamma }_{{{{{{{{\rm{eff,0}}}}}}}}}^{{\prime} }=12.09\) and \(\langle {{{{{{{{\bf{V}}}}}}}}}_{{{{{{{{\rm{ext}}}}}}}}}^{2}\rangle={C}_{{{{{{{{\rm{EFG}}}}}}}}}(0)\). The integration of C_EFG(t) over tens of picoseconds is necessary to obtain well-converged correlation times τ_c (Supplementary Fig. 16), notably at high salt concentrations and low temperatures (Fig. 2). Finally, our estimates in Supplementary Note 3 for the dipole-dipole contribution to the ²³Na⁺ rate 1/T₁ due to interactions with the spins of ¹H, ²³Na, and ³⁵Cl are more than four orders of magnitude smaller compared to the quadrupolar contribution, indicating that the latter dominates the ²³Na NMR relaxation.

\({T}_{1}^{-1}\) of ²³Na⁺ in aqueous NaCl as a function of the salt concentration c as obtained in experiments (E, solid lines and open symbols) and simulations (S, filled symbols) for different temperatures. The error bars for simulation results are associated with the approximation for incorporating the electron cloud contribution to the EFG.

Full size image

The NMR relaxation rate 1/T₁ grows with increasing salt concentration c and with reducing temperature T (Fig. 3). Under the extreme narrowing condition, which is fulfilled for the considered cases (see Methods), 1/T₁ ∝ τ_c, thereby suggesting that the slowing down of EFG fluctuations (Fig. 2), as reflected in the augmented correlation time, determines the rate behavior. Experimentally, 1/T₁ rises by about 50% within the considered range of concentrations c = 0.17–5.1 m for temperatures T = 20–50 °C, in line with the previous results^22,28,55,56. At T = 30 °C, 1/T₁ increases from around 15.9 s⁻¹ at c = 0.17 m to 25.2 s⁻¹ at c = 5.1 m. With increasing T from 25 °C to 50 °C, 1/T₁ reduces by >25% for considered salt concentrations. In general, our computational results for 1/T₁ of ²³Na⁺ agree well with the experimental data, especially at lower salt concentrations c ≲ 2.5 m, reproducing both the concentration and temperature behavior. For higher salt concentrations, we find that 1/T₁ grows systematically faster with increasing c as compared to the experiments, yet the relative error remains less than 15% over the considered range of conditions. The latter difference is likely caused by the shortcomings of the employed FF in capturing dynamic properties of the solution for c ≳ 2 m⁵².

Microscopic parameters of the relaxation

We find that the slowing down of EFG fluctuations at the Na⁺ position primarily causes a marked increase in the quadrupolar NMR relaxation rate with increasing c and decreasing T (Fig. 3). In Fig. 4, we quantify the role of dynamic and static effects that are reflected in the changes of τ_c and 〈V²〉, respectively, with varying salt concentration and temperature, as obtained in MD simulations of aqueous NaCl (see Supplementary Fig. 20 for other electrolyte solutions). While τ_c increases by a factor of ~1.5–2.5 with increasing c and decreasing T within the considered range of parameters (Fig. 4a), the value of 〈V²〉 reduces concurrently by up to 10% (Fig. 4c), indicating that the augmented correlation times are mainly responsible for the rate behavior.

a Effective correlation time τ_c of the EFG fluctuations at the Na⁺ position as a function of salt concentration c for different temperatures T. b Relative contribution of the slow EFG relaxation mode to τ_c for different c and T. The contribution was estimated using the stretched exponential fit of the normalized EFG ACFs (Supplementary Note 9). c Variance of the total EFG at the ion position \(\langle {{{{{{{{\bf{V}}}}}}}}}^{2}\rangle={(1+{\gamma }_{{{{{{{{\rm{eff,0}}}}}}}}}^{{\prime} })}^{2}\langle {{{{{{{{\bf{V}}}}}}}}}_{{{{{{{{\rm{ext}}}}}}}}}^{2}\rangle\) as a function of c for different temperatures. The legend in c is the same as in a. d EFG variance, 〈V²〉, evaluated from water molecules and ions located within a different number of solvation shells around the central Na⁺ ion as a function of c at T = 25 °C. The standard error from multiple independent simulation runs is either explicitly shown or does not exceed the symbol size.

Full size image

For considered c and T, τ_c of Na⁺ is quite short and below 1 ps (Fig. 4a), a feature already pointed out in previous classical^35,38,39,42 and ab initio^30,32 MD studies at infinite dilution. At T = 25 °C, we find that τ_c increases from 0.41 ps at c ≈ 0.06 m to 0.65 ps at c = 4 m. Despite the rapid decorrelation of EFG ACFs for t ≲ 0.2 ps (Fig. 2), we find that the contribution of the slow relaxation process to τ_c yields >85% of its overall value and also grows with increasing c and decreasing T (Fig. 4b). The dominance of the slow non-exponential decay of EFG ACFs over the τ_c behavior again exemplifies the governing role of collective processes behind the quadrupolar Na⁺ relaxation.

While the EFG variance at the Na⁺ position is largely determined by the first solvation shell contribution (Fig. 4d), a quantitative understanding of the QCC is only achieved if we take into account point charges within a radius of r ≳ 8 Å around the central ion, approximately the length scale of pronounced ion-ion and ion-solvent correlations (Supplementary Figs. 3 and 4). Similarly, we find that the EFG relaxation dynamics is well captured by point charge contributions located within the first two solvation shells around the ion, whereas the EFG due to the first solvation shell relaxes much more slowly (Supplementary Fig. 8). The first two solvation shells of Na⁺ are predominantly populated by water molecules even at the highest c = 4 m considered (see Supplementary Note 5), suggesting that the solvent provides the largest contribution to the EFG at the Na⁺ position and that other ions mostly retard water dynamics.

We observe that 〈V²〉 is reduced in bipyramidal complexes with octahedral symmetry, coordinated by six water molecules, yet only by 10% compared to the ensemble average (Supplementary Note 7). The contribution of the first solvation shell to the EFG variance features an increase with c (Fig. 4d), correlated with the fact that the six-coordinated state becomes less likely with increasing the salt concentration (Supplementary Fig. 7). Our consistently calculated QCC for ²³Na⁺ in aqueous NaCl is in the range between 19 × 10⁶ and 20.6 × 10⁶ rad × s⁻¹ for considered c and T (Supplementary Fig. 17), a value ~3−4 times larger than previous estimates based on the assumption that the EFG primarily decorrelates by translational and reorientational water dynamics with τ_c ≈ 3 − 7 ps^{19,20,21,22,57}. We thus conclude that the aforementioned modes of motion provide only a minor contribution to the observed relaxation.

Assessment of the relaxation models

We utilize the information available in experiments and molecular simulations in Fig. 5 to shed light on the mechanisms behind the quadrupolar relaxation. First, we focus on the possibility to model the EFG correlation time τ_c using the commonly used Stokes–Einstein–Debye (SED) relation \({\tau }_{{{{{{{{\rm{c}}}}}}}}}^{{{{{{{{\rm{SED}}}}}}}}}=4\pi \eta {r}_{0}^{3}/3{k}_{{{{{{{{\rm{B}}}}}}}}}T\), where η is the dynamic viscosity of the solution, r₀ is the sodium’s hydrodynamic (Stokes) radius, and k_B is the Boltzmann constant. Within the SED picture, the EFG relaxation at the Na⁺ position is governed by the Brownian rotational diffusion, likely to be related with collective reorientations of ion-water solvation complexes⁵⁵. While the SED model assumptions are not expected to hold down to the molecular scale^58,59, we systematically explore \({\tau }_{{{{{{{{\rm{c}}}}}}}}}^{{{{{{{{\rm{SED}}}}}}}}}\) in relation to τ_c, as it is often exploited to rationalize quadrupolar relaxation dynamics of ²³Na^+ 19,20,21.

a Na⁺ diffusion coefficient as a function of salt concentration c for different temperatures T in experiments (open symbols, legend shown in a) and simulations (filled symbols, legend shown in b). b Dynamic solution viscosity as a function of c for different T in experiments (open symbols, legend shown in a) and simulations (filled symbols, legend shown in b). Experimental viscosities were taken from ref. 60. Solid lines in a, b are polynomial fits of the experimental data. c Stress tensor ACFs C_stress(t) normalized by the system volume V and k_BT for increasing C at T = 25 °C. The black dotted line shows the stretched exponential fit of the long-time decay at c = 4 m. Inset: the time scale of solution structural relaxation τ_struct, as extracted from the long-time decay of C_stress(t) (see main text), as a function of c at multiple temperatures. T decreases from top to bottom, and the legend is shown in f, g. d Stokes-Einstein-Debye time plotted versus the EFG correlation time \({\tau }_{{{{{{{{\rm{c}}}}}}}}}^{\exp }\), as extracted from experimental data for different c and T. e \({\tau }_{{{{{{{{\rm{c}}}}}}}}}^{\exp }\) as a function of η/k_BT for different temperatures. \({r}_{0}^{{{{{{{{\rm{eff}}}}}}}}}\) is the effective hydrodynamic radius of a Na⁺ ion extracted under assumption that \({\tau }_{{{{{{{{\rm{c}}}}}}}}}^{\exp }\) can be modeled by a SED relation (the gray line shows the best fit). d, e Share the same legend. f Mean water dipole reorientation time τ_dip plotted versus the EFG correlation time τ_c, as extracted in simulations for different c and T. g τ_struct plotted versus τ_c for different c and T. The gray line indicates the linear dependence, τ_struct = τ_c. a–c, f, g The standard error from independent simulation runs are either explicitly shown or do not exceed the symbol size.

Full size image

We use the translational Stokes-Einstein relation D = k_BT/6πηr₀ to determine the concentration- and temperature-dependent values of the Stokes radius from the experimental Na⁺ diffusion coefficients (Fig. 5a) and highly accurate NaCl viscosity values provided by Kestin et al.⁶⁰ (Fig. 5b). D and η calculated in our MD simulations (see Methods) are in good agreement with the experiments, especially for c ≲ 2 m, capturing both the concentration and temperature behavior (compare filled and open symbols in Fig. 5a, b). The viscosity η in MD was obtained via the Green-Kubo formula using the stress ACFs (Fig. 5c), as detailed in Eqs. (3) and (4) in Methods. In Fig. 5d, we compare \({\tau }_{{{{{{{{\rm{c}}}}}}}}}^{{{{{{{{\rm{SED}}}}}}}}}\) calculated from state-dependent Stokes radii r₀(c, T) against the effective EFG correlation time \({\tau }_{{{{{{{{\rm{c}}}}}}}}}^{\exp }\) obtained from the experimental NMR relaxation rates and the sodium’s QCC from simulations (Supplementary Fig. 17) rather than those from previous estimates^{19,20,21,22,57}. r₀ assumes values between 1.5 and 2.0 Å for considered parameters (Supplementary Fig. 19). While both \({\tau }_{{{{{{{{\rm{c}}}}}}}}}^{{{{{{{{\rm{SED}}}}}}}}}\) and \({\tau }_{{{{{{{{\rm{c}}}}}}}}}^{\exp }\) generally lengthen with increasing c and decreasing T, \({\tau }_{{{{{{{{\rm{c}}}}}}}}}^{{{{{{{{\rm{SED}}}}}}}}}\) exceeds \({\tau }_{{{{{{{{\rm{c}}}}}}}}}^{\exp }\) by a factor of 8–17 (Fig. 5d). Similar results are obtained in our simulations (Supplementary Fig. 19). At T = 25 °C, \({\tau }_{{{{{{{{\rm{c}}}}}}}}}^{{{{{{{{\rm{SED}}}}}}}}}\) increases from around 7.5 to 8.2 ps, larger by more than one order of magnitude than \({\tau }_{{{{{{{{\rm{c}}}}}}}}}^{\exp }\) that grows from 0.44 to 0.66 ps for increasing c from 0.17 to 5.1 m. Thus, we conclude that the EFG correlation times cannot be understood on the basis of the SED relation parameterized using the translational hydrodynamic radius of sodium ions r₀ = k_BT/6πηD.

This is further illustrated in Fig. 5e showing \({\tau }_{{{{{{{{\rm{c}}}}}}}}}^{\exp }\) plotted against η/k_BT for various temperatures. While a Stokes-Einstein-like relation holds for \({\tau }_{{{{{{{{\rm{c}}}}}}}}}^{\exp }\), that is a strong correlation \({\tau }_{{{{{{{{\rm{c}}}}}}}}}^{\exp }\propto \eta /{k}_{{{{{{{{\rm{B}}}}}}}}}T\) exists for the considered range of parameters, the effective Stokes radius \({r}_{0}^{{{{{{{{\rm{eff}}}}}}}}}=0.69\) Å that would correspond to the EFG correlation time \({\tau }_{{{{{{{{\rm{c}}}}}}}}}^{\exp }\) within the SED model is clearly unphysical and smaller than the ionic radius 1.02 Å. \({r}_{0}^{{{{{{{{\rm{eff}}}}}}}}}\) was obtained from the fit \({\tau }_{{{{{{{{\rm{c}}}}}}}}}^{\exp }=4\pi \eta {\left[{r}_{0}^{{{{{{{{\rm{eff}}}}}}}}}\right]}^{3}/3{k}_{{{{{{{{\rm{B}}}}}}}}}T+{\tau }_{{{{{{{{\rm{0}}}}}}}}}^{{{{{{{{\rm{eff}}}}}}}}}\) with an additional intercept \({\tau }_{{{{{{{{\rm{0}}}}}}}}}^{{{{{{{{\rm{eff}}}}}}}}}=0.11\) ps needed for the best data representation⁵⁹. We obtain a similar value of \({r}_{0}^{{{{{{{{\rm{eff}}}}}}}}}\approx 0.68\) Å from our MD simulations (Supplementary Fig. 19). Therefore, the validity of the relation \({\tau }_{{{{{{{{\rm{c}}}}}}}}}^{\exp }\propto \eta /{k}_{{{{{{{{\rm{B}}}}}}}}}T\) explains the correlation between \({\tau }_{{{{{{{{\rm{c}}}}}}}}}^{\exp }\) and D⁻¹ reported in refs. 19,20,21, rather than simplified assumptions of the rotational Brownian diffusion that yield much larger estimates of \({\tau }_{{{{{{{{\rm{c}}}}}}}}}^{\exp }\) (Fig. 5d).

We now return to microscopic time scales of molecular motion in relation to that EFG fluctuations. The average water dipole reorientation time \({\tau }_{{{{{{{{\rm{dip}}}}}}}}}=\int\nolimits_{0}^{\infty }{{{{{{{\rm{d}}}}}}}}t\,\langle {P}_{1}[{{{{{{{\bf{u}}}}}}}}(t)\cdot {{{{{{{\bf{u}}}}}}}}(0)]\rangle\) assumed to drive the quadrupolar relaxation within the Hertz model^17,26,27,38 is 11-14 times larger compared to τ_c, as extracted in our simulations (Fig. 5f and Supplementary Fig. 18). Above, u is a unit vector pointing along the HOH bisector of a water molecule and P₁(x) = x is the first Legendre polynomial. This indicates that the single molecule reorientation with neglected intermolecule cross-correlations cannot explain the EFG relaxation dynamics. Yet, as seen in Fig. 5f, both τ_dip and τ_c increase similarly with increasing c and decreasing T, suggesting that the overall deceleration of the electrolyte dynamics, marked by an enhanced viscosity, impacts in a similar way both the motions that drive water reorientation as well as those that cause the EFG relaxation at the ion position.

To illustrate the relationship between these effects, in our MD simulations we extract a typical time scale of solution structural relaxation τ_struct using the stress tensor ACFs (Fig. 5c). While the short-time behavior of C_stress(t) corresponding to elastic, vibrational contributions features little changes with varying c and T⁶¹, its long-time tail slows down with increasing c and decreasing T, indicating an overall deceleration of the viscous dynamics of the liquid. We find that the long-time tail can be modeled well using a stretched exponential decay, \(\sim {e}^{-{(t/{\tau }_{{{{{{{{\rm{K}}}}}}}}})}^{{\beta }_{{{{{{{{\rm{K}}}}}}}}}}}\), with β_K ≈ 0.61 ± 0.04 (consistent with earlier simulations of pure water⁶¹ and time-resolved spectroscopy experiments⁶²). τ_K is in the range between 0.13 and 0.70 ps for considered parameters. The mean structural relaxation time \({\tau }_{{{{{{{{\rm{struct}}}}}}}}}={\tau }_{{{{{{{{\rm{K}}}}}}}}}{\beta }_{{{{{{{{\rm{K}}}}}}}}}^{-1}{{\Gamma }}({\beta }_{{{{{{{{\rm{K}}}}}}}}}^{-1})\), defined through the integral of the stretched exponential expression, is strongly correlated and comparable to the subpicosecond EFG correlation time τ_c (Fig. 5g). While the stress and EFG tensors are not directly related to each other, both quantities are inherently collective, that is the relaxation of their fluctuations is mainly driven by many-body correlations, and features a similar stretched decay for t ≳ 0.4 ps. All these observations suggest that the fast collective dynamics of the liquid that drive its structural rearrangements are also responsible for the quadrupolar NMR relaxation.

We have shown that the multiscale methodology combining DFT PAW calculations to parameterize the QCC and classical MD simulations to sample long-time EFG fluctuations enables an accurate description of the quadrupolar NMR relaxation rates of ²³Na⁺ in aqueous electrolyte solutions over a broad range of salt concentrations and temperatures. The resulting NMR relaxation rates are in very good agreement with the experimental data, especially at low salt concentrations, as validated in aqueous NaCl at multiple system state points. We find that the growth of the relaxation rate \({T}_{1}^{-1}\) with increasing c and decreasing T is primarily due to the slowing down of the EFG fluctuations reflected in the augmented EFG correlation time τ_c, while the concurrent changes in the QCC are rather small. The availability of dynamic information over a broad range of system parameters enabled us to have a consistent discussion concerning the quadrupolar relaxation models. We have demonstrated that the commonly-assumed rotational relaxation models based on either the water dipole reorientation^26,27,28 or on the Stokes-Einstein-Debye relation^19,20,21 overestimate the consistently-determined τ_c by at least an order of magnitude. This disagreement is understandable as these models restrict the relaxation description to one- or two-body correlations, oversimplifying the collective dynamics of the intermolecular EFG at the ion position^17,29,38. The quantitative interpretation of the EFG correlation times in terms of such simple isotropic models should therefore be used with caution. Instead, our results indicate that the Na⁺ EFG relaxation is largely determined by the dynamics in the two first solvation shells of the ion and occurs over a subpicosecond time scale comparable to that of solution structural rearrangements τ_struct, as extracted from the relaxation of the stress tensor. This again invalidates a continuous-solvent hydrodynamic description assuming that τ_c ≫ τ_struct.

Rather than directly probing single ion diffusion or single water molecule reorientation, our results suggest that the quadrupolar NMR relaxometry of ²³Na⁺ may be used as a complementary tool to analyze electrolyte dynamics in the THz ___domain, as relevant for emerging sodium-ion battery technologies^9,10,11. As quadrupolar relaxation is largely determined by the processes in the immediate vicinity of the solute, it can provide supplementary information on the fast, collective, molecular motions in ionic solvation cages that have been associated with the high-frequency dielectric response⁶³, solvation dynamics⁶⁴, as well as structural relaxation⁶¹ in aqueous electrolytes. The ability to capture the NMR relaxation rates by means of classical MD allows elucidating the quadrupolar relaxation mechanisms that occur in multicomponent systems, such as concentrated aqueous solutions of multiple salts⁶⁵, mixtures of salts with glycerol⁶⁶, or polyelectrolytes⁶⁷, where the relaxation dynamics may be influenced by environment heterogeneity, interface formation, microphase separation, or ion binding to polyelectrolyte chains. Future work could also focus on developing mesoscopic approaches that would allow a quantitative description of quadrupolar relaxation in complex biological-type compartments, characterized either by slow-motion conditions with dynamics within intracellular and extracellular spaces in biological tissues that may include structural anisotropy with residual quadrupolar coupling and distribution of correlation times (as, for example, in connective tissue where sodium ions are surrounded by a collagen matrix⁴). Such models would further provide a foundation for the interpretation of magnetic resonance imaging contrast mechanisms that are sensitive to the quadrupolar interaction^68,69,70,71.

Quadrupolar NMR relaxation rates

The quadrupolar mechanism dominates the relaxation of nuclei with spin I > 1/2 and is due to the coupling between their quadrupolar moment eQ with the EFG tensor V at the nucleus position¹⁶. While the NMR relaxation of spin components is generally bi-exponential for ²³Na with I = 3/2^22,72, it is possible to define effective longitudinal and transverse quadrupolar relaxation rates, 1/T₁ and 1/T₂, respectively, provided that the “fast motion” or “extreme narrowing” regime holds^30,73,74. In this case, the characteristic EFG correlation time τ_c is much smaller than the Larmor period \({\omega }_{0}^{-1}\) of the nucleus, ω₀τ_c ≪ 1. The latter can be shown to be fulfilled for all cases considered in this work as the relevant correlation times for ²³Na in electrolyte solutions are below 100 ps, and the magnetic field used in the experiments is 11.7 T that corresponds to \({\omega }_{0}^{-1}\approx 7.6\) ns. As we show in more detail in the SI, the two quadrupolar relaxation rates become equal in the fast motion regime and, combined with the rotational invariance of the system, can be expressed as³⁸

where ℏ is the reduced Planck constant, τ_c is an effective correlation time of EFG fluctuations

where \(\left\langle {{{{{{{\bf{V}}}}}}}}(0):{{{{{{{\bf{V}}}}}}}}(t)\right\rangle={\sum }_{\alpha,\beta }\langle {V}_{\alpha \beta }(0){V}_{\alpha \beta }(t)\rangle\) with α, β = x, y, z and the brackets 〈…〉 denoting an ensemble average, and \(\langle {{{{{{{{\bf{V}}}}}}}}}^{2}\rangle \equiv \left\langle {{{{{{{\bf{V}}}}}}}}(0):{{{{{{{\bf{V}}}}}}}}(0)\right\rangle\) is the EFG variance at the ion position. For a ²³Na nucleus with I = 3/2 and Q = 104 × 10⁻³¹ m^2 75, the rate constant 1/T₁ can be recast as \(1/{T}_{1}={C}_{Q}^{2}{\tau }_{{{{{{{{\rm{c}}}}}}}}}/10\) with the quadrupolar coupling constant (QCC) defined as \({C}_{Q}^{2}=\frac{2}{3}{\left(\frac{eQ}{\hslash }\right)}^{2}\langle {{{{{{{{\bf{V}}}}}}}}}^{2}\rangle\)^16,30. Finally, Eq. (1), which follows from linear response theory, allows calculating the quadrupolar spin-lattice relaxation rate 1/T₁ from the EFG fluctuations in equilibrium MD simulations without an imposed magnetic field.

NMR experiments

Solution samples with nine different NaCl concentrations were prepared by mixing x mg of NaCl in (y–x) mg of deionized water in a beaker, with x = 0.1, 0.2, 0.5, 0.8, 1.1, 1.4, 1.7, 2.0, 2.3 mg and y = 10 mg, to make solutions of concentrations 0.173, 0.349, 0.901, 1.488, 2.115, 2.786, 3.505, 4.278, and 5.111 mol kg⁻¹ to 5 mm NMR tubes (sample volume = 0.5 mL). All mass measurements were performed on a Mettler Toledo ME204E balance with a resolution of 0.1 mg. The solution at 26% weight corresponds to NaCl saturation in water at 20 °C⁷⁶.

NMR experiments were performed on an 11.7 T NMR Bruker Avance I spectrometer operating at 132.3 MHz for ²³Na, using a 5 mm double resonance broadband probe. The test tubes with the solutions were placed inside the spectrometer where the sample temperature could be controlled using gas flow and a temperature sensor providing precise, stable, and reliable temperature regulation. After each desired temperature reached steady state, a standard free induction decay was acquired followed by a longitudinal relaxation time T₁ mapping sequence, and a diffusion pulse sequence. At each temperature, the tuning and matching was checked. The duration of the 90 ° pulse was 9.6 μs, whereas that for the 180 ° pulse was 19.6 μs. A standard inversion-recovery pulse sequence was used to acquire T₁ with 32 logarithmically spaced steps. The delay was varied from 1 ms to 400 ms for ²³Na. Diffusion coefficients were measured using a Pulsed-Gradient-Spin-Echo in 32 steps with a maximum b-value of 2200 s mm⁻². The maximum diffusion gradient was 1 T m⁻¹ and the duration was 4 ms.

Molecular dynamics simulations

Aqueous sodium chloride (NaCl), bromide (NaBr), and fluoride (NaF) solutions were simulated using classical MD employing the Madrid-2019 FF⁵¹ that is based on the TIP4P/2005 water model⁷⁷ and uses scaled charges of +0.85e and −0.85e (e is the fundamental unit of charge) for Na⁺ cations and Cl⁻, Br⁻, and F⁻ anions, respectively. The FF parameters are listed in the SI. The scaled ionic charges aim at taking into account the electronic contribution to the dielectric constant at high frequencies in a mean-field fashion⁷⁸. At a moderate computational cost in comparison to fully polarizable models, the EFG relaxation within the Madrid-2019 FF⁵¹ has recently been shown to accurately describe the quadrupolar NMR relaxation rates of alkali metal ions at infinite dilution⁴², in particular, that of Na⁺. Solutions comprised of N = 1000 water molecules and N_p ion pairs were initialized at different salt concentrations c between 0.06 m (N_p = 1) and 4 m (N_p = 72) in a cubic box at the equilibrium solution density ρ(c, T) obtained in NPT simulations at P = 1 bar. The densities are in excellent agreement with the experimental ones, as discussed in Supplementary Note 2.

The equilibrated electrolyte systems were then simulated in the NVT ensemble. Both NPT and NVT simulation runs were carried out in the open-source MetalWalls package on graphics processing units⁷⁹ with electrostatic interactions computed with Ewald summation⁸⁰ and a short-range cutoff of 1.24 nm. The constant temperature was maintained using the Nose-Hoover chains thermostat with a time constant of 1 ps. System temperatures in range from 10 °C to 50 °C were considered. The equations of motion were integrated using the velocity Verlet algorithm and an integration time step of 1 fs. The effective rigidity of water molecules was imposed with the help of the RATTLE algorithm with a precision of 10⁻⁹. For each (c, T) state point, at least five independent runs of length 5 ns were performed to measure the EFG at the ion positions (sampled every 50 fs). Full Ewald summation expressions⁸⁰ were used in the computation of the EFGs, as recently implemented in MetalWalls⁴². For the considered system parameters, the relaxation of EFG fluctuations was found not to be affected by the finite box size, as we discuss in Supplementary Note 4.

Ab initio calculations

To determine EFGs with the electron cloud contribution, smaller systems containing 55 water molecules and N_p = 1, 2, 3, 4, and 5 NaCl ions pairs, corresponding to the salt concentrations c = 1, 2, 3, 4, 5 mol kg⁻¹, were simulated in the same way as the larger ones using the Madrid-2019 FF. In a single NVT simulation run at T = 25 °C, 2000 configurations were sampled with a period of 10 ps, and were later used in DFT-based EFG calculations with periodic boundary conditions in the Quantum Espresso (QE) package⁸¹. No additional geometry optimization of the configurations was performed in the DFT calculations. The pseudopotential-based projector-augmented wave (PAW) method^45,47,48 was used to reconstruct the all-electron charge density in the vicinity of the nucleus using the QE-GIPAW package⁸². The self-consistent electron densities were calculated using the PBE functional⁸³, a kinetic energy cutoff of 80 Ry, and norm-conserving pseudopotentials of the GIPAW package⁸⁴. In the case of Na⁺ ions, the EFGs obtained with the PBE functional were shown to be in good agreement³² with those obtained with the hybrid PBE0 functional⁸⁵.

Dynamical properties of electrolyte solutions

The shear viscosity of aqueous electrolyte solutions was obtained using the Green-Kubo relation⁸⁶:

with V being the system volume and k_B standing for the Boltzmann constant. The stress tensor ACF C_stress(t) was computed as⁸⁶

where α, β run over the three Cartesian components and P_αβ is the traceless symmetrized part of the stress tensor σ_αβ: \({P}_{\alpha \beta }=\frac{1}{2}\left({\sigma }_{\alpha \beta }+{\sigma }_{\beta \alpha }\right)-\frac{1}{3}{\delta }_{\alpha \beta }{\sum }_{\gamma }{\sigma }_{\gamma \gamma }\). For each salt concentration, the viscosity was measured over more than five independent simulation runs of length 5 ns with the stress tensor sampled every integration time step (1 fs).

The Na⁺ diffusion coefficients were extracted from the long-time limit of the ion’s mean-square displacement:

where N_p is the number of sodium ions in the system, r_i(t) is the position of the i-th ion at time t, and the brackets 〈⋯〉 stand for ensemble averaging. The obtained diffusion coefficients were corrected for finite-size effects using the Yeh-Hummer relation⁸⁷:

with the diffusion coefficient D_∞ corresponding to a macroscopic system, D being obtained in a cubic simulation box with side length L, and ξ ≈ 2.837297. The calculated values of viscosity η in Eq. (3) were used for evaluating D_∞ in Eq. (6). The finite-size correction term corresponded to 17–22% of the measured value D.