Quadrupolar 23Na+ NMR relaxation as a probe of subpicosecond collective dynamics in aqueous electrolyte solutions

Nuclear magnetic resonance relaxometry represents a powerful tool for extracting dynamic information. Yet, obtaining links to molecular motion is challenging for many ions that relax through the quadrupolar mechanism, which is mediated by electric field gradient fluctuations and lacks a detailed microscopic description. For sodium ions in aqueous electrolytes, we combine ab initio calculations to account for electron cloud effects with classical molecular dynamics to sample long-time fluctuations, and obtain relaxation rates in good agreement with experiments over broad concentration and temperature ranges. We demonstrate that quadrupolar nuclear relaxation is sensitive to subpicosecond dynamics not captured by previous models based on water reorientation or cluster rotation. While ions affect the overall water retardation, experimental trends are mainly explained by dynamics in the first two solvation shells of sodium, which contain mostly water. This work thus paves the way to the quantitative understanding of quadrupolar relaxation in electrolyte and bioelectrolyte systems.

In this section, we provide a more detailed derivation of the expressions used for calculating the quadrupolar relaxation rates in Section II A of the main text. As given by the formalism of Spiess et al. [1][2][3], the longitudinal 1/T 1 and transverse 1/T 2 quadrupolar relaxation rates of nuclei with spin I > 1/2 can be determined through where ℏ is the reduced Planck constant, eQ is the quadrupolar moment of the nucleus, and ω 0 = γB 0 is its Larmor frequency of the nucleus in the external magnetic field B 0 with γ denoting its gyromagnetic ratio. G Q 2,0 (ω) and G Q 2,±1 (ω) are spectral densities of the EFG at the nucleus position (see below). As the NMR relaxation of spin components may generally not be singleexponential [4], strictly speaking, the longitudinal and transverse quadrupolar relaxation rates above can be only defined under the following conditions [1]: (i) for the nuclear spin I = 1; (ii) provided that the spin subsystem can be described by a sufficiently high spin-temperature; (iii) in the extreme narrowing regime, for which the typical microscopic EFG correlation time τ c is much smaller in comparison to the Larmor period ω −1 0 of the nucleus, yielding ω 0 τ c ≪ 1. As we discuss in the main text, the latter condition is in fact fulfilled for 23 Na in aqueous solutions considered in this work and implies that the spectral densities in Eqs. (1) and (2) can be effectively evaluated at ω = 0.
The EFG spectral densities G Q 2,0 (ω) and G Q 2,±1 (ω) can be explicitly written as [1]: 2g Q 2,±1 (±ω) + 3g Q 2,0 (0), where g Q l,m are half-sided Fourier transforms of the autocorrelation functions of the spherical EFG components V l,m (l = 2, m = 0, ±1, ±2): with the brackets ⟨. . . ⟩ denoting the time averaging and V l,m for l = 2 given by [3] V 2,0 = 3 2 V zz (6) where V αβ (α, β = x, y, z) stand for the Cartesian components of the EFG tensor at the nucleus. The rotational invariance of the considered system implies that all autocorrelation functions ⟨V l,m (t)V * l,m (0)⟩ in Eq. (5) for different m become equal [1]. Combined with the extreme narrowing of the signal, we find from Eqs. (3) and (4) that G Q 2,0 (0) = G Q 2,±1 (0) and can be recast as follows: where ⟨V(t):V(0)⟩ = α,β ⟨V αβ (t)V αβ (0)⟩ with α, β = x, y, z. Accordingly, the two effective quadrupolar rates 1/T 1 and 1/T 2 also equalize and are given by The Eq. (10) above can be rewritten in terms of the variance of the EFG tensor ⟨V 2 ⟩ ≡ ⟨V(0):V(0)⟩ and the effective EFG correlation time τ c : yielding the final expression Similarly to the intramolecular quadrupolar relaxation, the value of V 2 can be recast in terms of the so-called quadrupolar coupling constant (QCC) C Q [7]. The latter is defined through the the three eigenvalues of the instantaneous EFG tensor at the ion position, V 33 , V 22 , and V 11 , that are ordered in a descending manner, |V 33 | ≥ |V 22 | ≥ |V 11 |. It can be shown that the ensemble average V 2 satisfies permitting us to define the value (eq) 2 ≡ V 2 33 and the asymmetry parameter η Q ≡ (V 22 − V 11 ) 2 / V 2 33 . Finally, the QCC can be defined as which enables to rewrite Eq. (12) as follows: It should be noted that for a purely intramolecular mechanism the quadrupolar relaxation is primarily driven by the molecule reorientation [7]. In such case, the components of the EFG tensor are fixed in a reference frame associated with the molecule that makes the values of eq and η Q discussed above constant. On the other hand, for an intermolecular mechanism, as for monoatomic quadrupolar ions considered here, the components of the EFG tensor at the nucleus position are not fixed in any reference frame, but are stochastic quantities, typically satisfying a Gaussian distribution [8]. Precisely because of the latter fact, in the present case eq and η Q are determined from the statistically averaged values of the eigenvalues.

Supplementary Note 2: Force field parameters
As mentioned in the main text, aqueous NaCl solutions were simulated using the recently-developed Madrid-2019 [5,6] FF for aqueous electrolytes that employs scaled ionic charges. The electrolyte FF is based on the TIP4P/2005 water model [9]. The interparticle interaction energy V (r ij ) in the FF is given by a sum of Coulomb and Lennard-Jones (LJ) potentials: where q i and q j the charges of the i th and j th atoms, r ij is the distance between two particles, ϵ 0 is the vacuum permittivity, ϵ ij is the minimum of the LJ potential, and σ ij is the LJ diameter. In the TIP4P/2005 model [9], the rigid water molecules are comprised of four interaction sites: an oxygen atom with q O = 0 that is used as a LJ interaction center, two hydrogens with q H = 0.5564e,  1. The charges of the Na + cations and Clanions are q Na = +0.85e and q Cl = −0.85e, respectively, with e denoting the fundamental unit of charge.
The use of scaled ionic charges in the Madrid-2019 FF [5,6] enables not only an excellent description of the solution densities, but also a good representation of the dy-  Table 2.
Number of ion pairs Np and respective salt concentrations c in the simulation boxes with aqueous NaCl solutions. In each case, the system contained Nw = 1000 water molecules. namic properties of electrolytes, as given by the solution viscosity, self-diffusion coefficients of ions and water, as well as the quadrupolar NMR relaxation rates for a series of ions at infinite dilution [10]. As we show in Fig. 1, the NaCl densities obtained explicitly in N P T simulations at P = 1 bar are in excellent agreement with the experimental values within the considered range of temperatures and concentrations [11]. To compute equilibrium NaCl densities in the Madrid-2019 FF, we performed N P T simulations of length 1 ns using the MetalWalls open-source code [12] by coupling the systems to chains of Nosé-Hoover barostats and thermostats [13][14][15], both with a time constant of 1 ps. Finally, at a given system temperature T and salt concentration c, NaCl solutions were initialized and equilibrated at the resulting solution densities ρ(c, T ) that corresponds to the 1 bar isobar and the production runs were performed in the N V T ensemble. In each case, N w = 1000 water molecules and N p ion pairs were simulated in a box at a given density ρ(c, T ) (precise parameters are listed in Tab. 2). In addition to the NaCl systems, we have also simulated aqueous NaBr and NaF solutions using the extended Madrid-2019 FF parameters [6] (listed in Tab. 1) using the same protocol as described above. NaBr solutions were simulated at the same salt concentrations as in the NaCl case, whereas NaF salts were simulated at c = 0.06, 0.5, and 1 m because of their rather small solubility in water [6]. The simulations were performed at T = 25°C using equilibrium densities obtained in N P T simulations at 1 bar. The latter solution densities are also in excellent agreement with the experimental ones within the considered range of salt concentrations [6].

Supplementary Note 3: Impact of dipole-dipole couplings
In this section, we will estimate a potential effect of the dipole-dipole coupling mechanism on the longitudinal relaxation rate T −1 1 of 23 Na + ions. Besides the quadrupolar contribution to the rate discussed in Sec. 1 that is Supplementary Table 3. Estimates for the contribution T −1

1X
to the longitudinal rate of 23 Na + that arises from the dipoledipole interaction between the spin of the sodium ion and the spin IX of X = 1 H, 23 Na + , and 35 Cl -. γX are the respective nuclear gyro-magnetic ratios and r is the location of the first peak of Na-X radial distribution functions (see Fig. 3). In all cases, we assumed that τ dd c ≈ 4 ps. expected to dominate the NMR relaxation for quadrupolar nuclei [7], additional additive contributions T −1 1X to T −1 1 might arise from the 23 Na + -X dipole-dipole spin interaction, where X = 1 H, 23 Na + , and 35 Cl -. In this case, the NMR relaxation for two unlike spins I and I X is determined from the dipole-dipole ACF G(t) [16,17]: where we assume that the spin I is relaxed by the spin I X , γ I and γ X are the two nuclear gyro-magnetic ratios, µ 0 is the vacuum permeability, ℏ is the Planck constant, r(t) is the distance between two spins at time t, θ(t) is the angle between the distance vector r(t) and the magnetic field direction at time t, and brackets ⟨· · · ⟩ τ denote ensemble averaging. Provided that the extreme narrowing condition holds, the longitudinal rate contribution T −1 1X is given by [16]: where G(0) is determined from Eq. (17) at t = 0, τ dd c is the correlation time of the dipole-dipole interaction that is evaluated as the integral of the normalized ACF G(t)/G(0), and the constant c X is equal to 3/2 in the case of like spins with I = I X and 1 otherwise [16]. To estimate G(0), let us assume that the angular part in Eq. (17) is uncorrelated to the radial one at t = 0, yielding ⟨ 3 cos 2 θ(τ ) − 1 2 ⟩ τ = 4/5 [16]. Thus, We estimate the value ⟨1/r 6 ⟩ from the first peak of the respective Na-X radial distribution functions (see Fig. 3 below) and assume that τ dd c ≈ 4 ps that corresponds to the time scale of intermolecular 1 H dipole-dipole relaxation in pure water at room temperature conditions [18].
The resulting values are listed in Tab. 3. The potential contribution to the NMR rate coming from the dipoledipole interaction with 1 H is the largest (≈ 2.7·10 −4 s -1 ), whereas those with 23 Na + and 35 Clare much smaller. Nevertheless, the overall dipole-dipole contribution is estimated to be orders of magnitude smaller than the quadrupolar rate (≈ 16 s -1 ) that we extract in our MD simulations at low salt concentrations.

Supplementary Note 4: Finite-size effects
To assess the finite-size effects on the EFG ACFs and the resulting NMR relaxation rates, we have simulated three systems at c = 4 m with N w = 125, 250 and 500 water molecules, in addition to the original one with N w = 1000, with 9, 18, 35, and 72 ion pairs, respectively, at T = 25°C. As seen in Fig. 2, the reduction of the system size does not influence the resulting ACF of the EFG at the ion position. Interestingly, neither we find any effect on the long-time tail of the ACF that is consistent with a hydrodynamic power law decay ∼t 5/2 [19]. The emergence of such algebraic decays is typically associated with a superposition of an infinite number of relaxation modes with wave vector k in a liquid [20]. In simulations with a box of size L and periodic boundary conditions, however, k is limited by the minimal value 2π/L, and infinitely small wave numbers k → 0 are not accessible. This causes an apparent exponential relaxation of the ACFs at long times that would otherwise be algebraic in a macroscopic system [20]. Nevertheless, within the considered time scales, we do not observe this effect on the EFG ACFs. To elucidate structural changes in aqueous NaCl solutions that arise with changing macroscopic parame-ters c and T in the Madrid-2019 FF, in Figs. 3 and 4 we show various radial distribution functions (RDFs) for T = 25°C for increasing the salt concentration c and at c = 1 m for increasing T , respectively. We characterize the solvation shell structure of Na + ion by means of the Na-O g Na−O (r) (Figs. 3a and 4a), Na-H g Na−H (r) (Figs. 3b and 4b), Na-Cl g Na−Cl (r) (Figs. 3c and 4c), and Na-Na g Na−Na (r) (Figs. 3d and 4d) RDFs. For completeness, we also show the Cl-O g Cl−O (r) (Figs. 3d and 4d) as well as Cl-Cl g Cl−Cl (r) (Figs. 3e and 4e) RDFs. We define the boundaries of the first and the second solvation shells of Na + from the first two minima of the Na-O RDF, g Na−O (r) (Figs. 3a and 4a). Despite a small drift of the first minimum position from around 3.15Å to 3.19Å in the range from c ≈ 0.06 m to c = 4 m, for the sake of simplicity we have considered the boundary of the first solvation shell to be at an intermediate value r 1s = 3.17Å in all cases. The boundary of the second solvation shell was considered to be at r 2s = 5.39Å.
The coordination numbers were obtained as where g(r) is the relevant RDF, n is the number density, and r c defines the integration limit. For instance, for the oxygen coordination number in the first solvation shell, g(r) = g Na−O (r), n is the number density of oxygen atoms in the simulation box, and r c = r 1s . For the first solvation shell of Na + at T = 25°C, we find that the oxygen coordination number decreases from around 5.55 to 5.46 for increasing c from 0.06 to 4 m. Concurrently, we find that the hydrogen coordination number decreases from 9.70 to 9.59. Within the two solvation shells of Na + , the number of oxygens decreases from 22.94 to 20.61 for increasing c from 0.06 to 4 m. Naturally, the latter behavior implies that water molecules in the first solvation shell are partly replaced by ions. We find that the latter are mostly due to Clanions, the number of which in the first solvation shell increases from about 0.01 to 0.08. In comparison, the number of Na + cations in the first solvation shell remains much lower (0.01) even at the highest concentration c = 4 m here. Furthermore, we find that at c = 4 m there is on average 0.8 cations and 1.9 anions within the second solvation shell of Na + . We find that with increasing c, there is an enhanced probability of observing more like-species in the solvation shells of Na + and Cl -, as evidenced by the increase in the peaks of g Na−Na (r) and g Cl−Cl (r) (Figs. 3d and 4f), whereas an opposite trend is seen for g Na−Cl (r), whose peaks decrease with increasing c. The latter behavior is consistent [21] with the trends found in the MD simulations of aqueous NaCl solutions that employed a deep neural network potential trained on the state-of-the-art DFT calculations using the strongly constrained and appropriately normed (SCAN) functional [22]. The temperature behavior of the RDFs at c = 1 m is shown in Fig. 4. For increasing T from 10°C to 50°C, the oxygen coordination number in the first solvation shell decreases from 5.55 to 5.45, whereas that of hydrogen from 9.80 to 9.43. Interestingly, we find that first peak of g Na−Cl (r) at 2.84Å rises with increasing T , corresponding to a small increase of the number of Clanions within the first solvation shell (from around 0.015 at T = 10°C to 0.027 at 50°C). The latter effect is, however, balanced by the decrease in the height of the second peak of g Na−Cl (r) at 4.6Å that leaves the overall number of Clanions within the two solvation shells (0.57) practically unchanged, see Fig. 4c. We classify different cation-anion solvation states by the minima of the Na-Cl RDF, g Na−Cl (r) (Figs. 3c and  4c). Given that r denotes the distance between a Na-Cl ion pair, we define the following three possibilities: (i) contact ion pairs (CIP) for r < 3.39Å; (ii) solvent-shared ion pairs (SIP) 3.39 < r < 5.76Å; (iii) solvent-separated ion pairs (SSIP) r > 7.9Å. The probability of observing the latter solvation states of Na + with increasing c at T = 25°C is shown in Fig. 5. We find that SIPs are the most likely states for c > 1 m, while there is a rather modest increase in the number of CIPs from few to around ∼10% in the range from 1 to 5 m. The SSIPs are very unlikely at salt concentrations above c = 3 m. By integrating g Na−Cl (r) up to 3.39Å, we find the number of CIPs to increase from around 0.01 at 0.5 m to 0.1 at 4 m at T = 25°C. In comparison to neural network MD based on the SCAN DFT calculations [21], we find that the computed number of CIPs is ∼2-3 times smaller. The latter is due to the fact that in our case the first peak of g Na−Cl (r) is much less pronounced than the second one (Figs. 3c and 4c), as compared to g Na−Cl (r) reported in Ref. [21]. In Fig. 6, we show the effective Sternheimer factors extracted on sodium states forming a CIP, SIP, and SSIP with the Clanion. The configuration-resolved γ eff were extracted from a set of at least 50 independent hydration shell configurations of Na + ions. We find that the ensemble-averaged values of the Sternheimer factor γ eff are consistent with those calculated from a subset of configurations forming the SIP, highlighted with red diamonds in Fig. 6, which constitutes the most probable Supplementary Figure 6. Configuration-resolved γ eff for Na + forming contact ion pairs (CIP, green triangles), solventshared ion pairs (SIP, red diamonds), and solvent-separated ion pairs (SSIP, violet pentagons) with the Clanion. The error bars were computed using bootstrapping.
solvation state of Na + cations for c > 1 mol·kg -1 (Fig. 5). The values of Sternheimer factors obtained on a subset of CIPs are somewhat enhanced in comparison to the ensemble average (green triangles in Fig. 6), whereas the ones found for SSIPs, are consistent with the infinite dilution value of γ eff (violet pentagons in Fig. 6).

Supplementary Note 7: Effect of the local solvation shell structure on the electric field gradient variance
The changes in the EFG variance ⟨V 2 ⟩, which determines the strength of the quadrupolar relaxation, are typically associated with the modification of the hydration shell structure and its symmetry [23][24][25]. For example, Versmold [23] demonstrated that the EFG at the position of a solute embedded into a dipolar liquid profoundly depends on the symmetry of the solvation shell and can be significantly reduced due to mutual cancellations for highly symmetric arrangements (e.g., in the case of tetrahedral or octahedral solvation shell structures). Furthermore, as the EFG is a rather short-range quantity that is proportional to r −3 and r −4 at the distance r away from a point charge and dipole, respectively, the dominant contribution to ⟨V 2 ⟩ is due to molecules located in the immediate neighborhood of the solute [8].
To elucidate the effect of the solvation shell structure of 23 Na + on the EFG variance at its position, in Fig. 7 we consider the EFG variance ⟨V 2 (N w )⟩ averaged over configurations with a different number of water molecules N w in its first solvation shell for different concentrations and temperatures. The sodium ion has a quite small ionic radius of 1.02Å that results in strong interactions with neighboring water molecules and the development of multiple coordination shell structures [21,26]. As seen in Figs. 7a and 7b and as discussed in Refs. [21,26], the latter are dominated by the square pyramidal and tri-angular bipyramidal complexes coordinated by N w = 5 water molecules, as well as by the 6-coordinated square bipyramidal state. The much less likely tetrahedral structures with N w = 4 are also present. In agreement with recent machine learning based simulations parameterized on the strongly constrained and appropriately normed functional [21,22], we find that the hydration shell state with N w = 6 is the most likely at low salt concentrations, yet is becoming somewhat less probable with increasing c at the expense of tetrahedral structures with N w = 4 as well as complexes with N w = 5 (Fig. 7a). The latter rearrangements are due to ions that are more likely to penetrate into the first two sodium's hydration shells with increasing c (see Figs. 3). Qualitatively similar trends are seen with decreasing T (Fig. 7b). Consistently at different concentrations and temperatures, we find that the octahedral hydration complexes formed by 6 water molecules in the first solvation shell feature a reduced value of the EFG variance at the Na + position relative to the ensemble average ⟨V 2 ⟩ (Figs. 7c and 7d). While consistent with the Versmold's symmetry argument [23], the magnitude of the reduction at hand does not exceed 10% in comparison to ⟨V 2 ⟩. Furthermore, the EFG variance computed for less coordinated polyhedra with N w = 4 and 5 are about 10% larger than ⟨V 2 ⟩. This highlights that the changes of the QCC in different Na + hydration complexes are rather small within the considered range of salt concentrations and temperatures. Nonetheless, the trend of the variance reduction in sodium's octahedral hydration complexes can be understood on the basis of EFG cancellations in highly symmetric environments [23].

Supplementary Note 8: Solvation shell resolved electric field gradient relaxation
In this section, we discuss the origin of the EFG fluctuations in terms of the microscopic environment of the Na + ion. In Supplementary Fig. 8, we show the EFG ACFs at T = 25°C and multiple concentrations, as computed from all point charges in the system using Ewald summation, from water molecules and ions residing in the first solvation shell only, and from water molecules and ions within the first two solvation shells. We find that the EFG ACF due to the first solvation shell decays much more slowly than the full EFG ACF, in particular with increasing salt concentration. In contrast, the EFG ACF due to point charges in the first two solvation shells provides a good description of the relaxation with the corresponding ACF capturing well both the shortand long-time dynamics (compare blue and black lines in Fig. 8). This indicates that the processes in the first two solvation shells are mainly responsible for the quadrupolar relaxation dynamics. In this section, we present additional results concerning the relaxation dynamics of the EFG at the position of Na + ions. In particular, in Fig. 9 we highlight the EFG ACFs for different temperatures and concentrations that supplement the results shown in the main text. Systematically for both the dilute and concentrated regimes, we find that the EFG ACFs have a characteristic form that consists of two steps: (i) fast initial decorrelation at short times (t ≲ 0.2 ps); (ii) a slower decay on the ps time scale. While the initial fast decay barely changes with c and T , as seen in Fig. 10 that shows the EFG ACFs for t < 1.5 ps, the main change in the EFG relaxation dynamics is reflected on the second decay of the ACFs that profoundly slows down with decreasing T and with increasing salt concentration c. The latter is particularly evident at high salt concentrations (see Fig. 9d-f).
The effective EFG correlation time τ c (11) is obtained directly from MD data by integrating the normalized ACFs. The running integrals of C EFG (t)/C EFG (0) are shown in Fig. 16. Due to a pronounced slow-down of EFG fluctuations with increasing c and decreasing T , the integration up to ∼50 ps may be necessary (e.g., at T = 10°C and c = 4 m) to obtain well-converged values of τ c . This again highlights the importance of long-time sampling in the description of the NMR relaxation rates in concentrated electrolyte solutions.
To assess the EFG relaxation behavior, we attempted to fit the normalized EFG ACFs, C EFG (t)/C EFG (0), using different functional relations. Given a two-step shape of the relaxation decay, we first modelled the EFG ACFs using a sum of two exponentials: where τ f is the time scale of the initial, fast process that occurs at t ≲ 0.2 ps, τ 1 is the time scale of the second, slow process, and α is its relative weight. The resulting fits using Eq. (21) are shown with dotted lines in Fig. 9 for different salt concentrations and temperatures (the average R 2 score is 0.995 over the ensemble of fitted curves). Furthermore, in Fig. 10 the same fits and EFG ACFs are plotted for t < 1.5 ps. Evidently, for both dilute salt concentrations (Fig. 9a) and in the concentrated regime (Fig. 9b-f) the two exponential fit (21) is only valid for t ≲ 2−3 ps, whereas the consecutive decay at longer times is slower. The latter highlights the nonexponential character of the relaxation and points to a collective pathway behind the EFG relaxation. Quantitatively, by performing a generic three-parameter fit in Eq. To improve the representation of EFG ACFs at long times, we have included an additional second slow process with a time scale τ 2 and weight α 2 in the fitting procedure: To reduce the number of independent fitting parameters in Eq. (22), we reused the values of τ f obtained in the fitting with Eq. (21). As seen in Fig. 11, the three exponential fit (21) provides a quite good representation of the relaxation dynamics of EFG fluctuations for most of the decay (more than two orders of magnitude relative to the initial value of the ACF, roughly up to ∼10 ps). The associated average R 2 score is 0.999 over the ensemble of fitted curves. Yet, the presence of even slower relaxation modes is evident in the concentrated solutions at very long times (Fig. 11d-f). For the considered c and T , the relative weight α 1 attains values 0.15-0.21, whereas α 2 around 0.1-0.14. Both time constants τ 1 and τ 2 grow with increasing c and decreasing T . For instance, at T = 25°C, τ 1 increases from 0.65 to 1.1 ps, whereas τ 2 from 2.21 to 3.95 ps. Since the fit (22) provides a good description of the bulk part of the decay, the effective EFG correlation times obtained by integrating the model expression (22), τ c, three exp. = (1 − α 1 − α 2 )τ f + α 1 τ 1 + α 2 τ 2 , is in good agreement with τ c obtained by directly integrating the EFG ACFs from MD.
We find that the EFG ACFs can be equally well represented by a simpler fit that is composed of a single exponential to model the fast initial decay and a stretched exponential to model to consecutive, slower part of the relaxation: where α s is the relative weight of the stretched exponential function, τ s is its time scale, and β is the stretching exponent. While the expression (23) generally consists of four independent fitting parameters, certain approximation can be applied: (i) we reuse the values of τ f obtained in the fitting of Eq. (21); (ii) we obtain the value of β = 0.67 ± 0.05 by fitting the tail of the EFG ACF at the lowest concentration considered, c = 0.06 m, and fix it in the following optimizations. This reduces the fit (23) to just two independent parameters, α s and τ s . The resulting fits shown with dotted lines are in good agreement with the measured EFG ACFs, as seen in Fig. 12 for multiple temperatures and concentrations. The average R 2 score is 0.998. Similarly to the three-exponential fit (22), the model employing the stretched exponential function (23) decays faster at long times in comparison to the EFG ACFs obtained in MD. This might indicate that the long-time EFG relaxation decay is slower than both exponential and stretched exponential relations. In fact, in the main text we show that the long-time tail of the EFG ACFs is consistent with a power-law ∼t −5/2 that was predicted within a mode coupling theory [19]. In Fig. 13, we show the concentration and temperature dependence of the optimized values of τ s and α s (Fig. 13a  and 13b, respectively) of the fit in Eq. (23). The parameter τ s mimics the behavior of τ c , yet adopts values that are about twice as high. α s , somewhat decreasing with increasing c and T , features values in range between 0.35 and 0.39. By integrating the expression (23), we find the effective, predicted EFG correlation time with Γ(x) denoting the gamma function. The values of τ c,stretched are in excellent agreement with the τ c results obtained by directly integrating the ACFs from MD. Finally, we find that the stretched exponential in Eq. (23) provides a dominant contribution to τ c . The latter is highlighted in Fig. 13d where we show β −1 Γ(β −1 )α s τ s scaled with τ c . It is evident that the second, slow part of the EFG ACF decay constitutes more than 85% of τ c , and its contribution increases with increasing c and decreasing T . The stretched exponential decay ∼e −[t/τs] β can be interpreted as a continuous superposition of exponential modes with a different rate constant k [27]: where H(k) is the distribution function that quantities the relative weight of each k-mode, satisfying +∞ 0 dk H(k) = 1. Formally, H β (k) can be obtained as the inverse Laplace transform of e −[t/τs] β . Analytical, closed-form expressions were derived for certain rational β, e.g. 1/3 and 2/3 [27]. Moreover, there exists a convergent series that can be used to determine H β (k) in the general case [27]. Nevertheless, the latter series features an oscillating behavior that complicates its use in practical numerical calculations [27]. Thus, to understand the behavior of H(k) in our case of the EFG ACFs (23), we have employed the approximate numerical expressions of Berberan et al. [27]: Both for increasing c at a fixed T (Fig. 14a) and for decreasing T at a fixed c (Fig. 14c), H(k) features a rather broad distribution with a single peak, whose maximum position shifts towards larger k. The inverse of the peak position of H(k), k −1 max , (Fig. 14b and 14d) generally follows the behavior of τ c and adopts values between 1 and 5 ps, typically 2-3 times larger than the actual value of the effective EFG correlation time τ c .
In Fig. 15, we additionally highlight the long-time behavior of the EFG ACFs by multiplying it by t 5/2 and t 3/2 . As seen in Fig. 15a, t 5/2 C EFG (t) saturates at long times, suggesting an algebraic decay ∼t −5/2 that is consistent with the mode coupling theory of Bosse et al. [19].
Supplementary Figure 13. Parameters of the fit in Eq. (23). a Time scale of the stretched exponential τs as a function of the salt concentration c for different temperatures T . b Weight of the stretched exponential contribution αs in the fit (23). c Effective EFG correlation time τc as obtained directly in MD (filled symbols) and by integrating the fitting expression (23) (solid lines). d Contribution of the stretched exponential αsτsβ −1 Γ(β −1 ) to τc. The standard error from multiple independent simulation runs is either explicitly shown or does not exceed the symbol size.

Supplementary Note 10: Quadrupolar coupling constant
In Fig. 17, we show the concentration and temperature dependence of the QCC C Q that is related to the EFG variance ⟨V 2 ⟩ through Eq. (14) as follows: C Q decreases with increasing c and decreasing T in range from around 19×10 6 to 20.6×10 6 rad·s −1 . At T = 25°C and low c, C Q = 19.87 × 10 6 rad·s −1 .
Supplementary Figure 17. QCC given by Eq. (14) for 23 Na + ions in aqueous NaCl solutions as a function of the salt concentration c for different temperatures indicated in the legend.

Supplementary Note 11: Water dipole reorientation
We quantify the mean time scale τ dip of water dipole reorientation by means of the integral where P 1 (x) = x is the first Legendre polynomial and u(t) is a unit vector at time t that points in the direction of the dipole of a TIP4P/2005 water molecule (i.e., from the virtual M-site towards the oxygen atom). The resulting concentration and temperature dependence of τ dip is shown in Fig. 18.

Supplementary Note 12: Quadrupolar relaxation within a Stokes-Einstein-Debye model
In Fig. 19a and 19b, we extract the hydrodynamics radius r 0 of Na + ions in aqueous NaCl solutions by means are more than an order of magnitude larger than the respective EFG correlation time τ c . In agreement with the experimental results in Fig. 5, we find that τ c in simulations satisfies a generalized Stokes-Einstein relation, τ c ∝ η/k B T (Fig. 19d). The associated effective hydrodynamic radius that can be obtained using a fit τ c = 4πη r eff The EFG relaxation of 23 Na + in NaBr, NaF, and NaCl is further discussed in Fig. 20 for different salt concentrations at T = 25°C. The dynamics of EFG fluctuations for 23 Na + in NaBr and NaF are quite similar to the case of NaCl discussed in detail in the main text, as we show with the normalized EFG ACF in NaBr for different salt concentrations in Fig. 20a and by comparing the EFG ACFs in NaCl, NaBr, and NaF at c = 1 m in Fig. 20b. Quite small differences in the slow part of the EFG ACFs (Fig. 20b) cause a discrepancy in the effective correlation time τ c (Fig. 20c). Whereas the magnitude of the EFG variance at the ion position ⟨V 2 ⟩ features a small decrease with increasing c (Fig. 20d), the lengthening of the effective correlation time τ c (Fig. 20c) results in the increase of the quadrupolar relaxation rate (Fig. 20e) in all cases. Supplementary Figure 20. EFG relaxation at the position of Na + ions in different electrolyte solutions at T = 25°C. a Normalized EFG ACFs at the position of Na + ions in aqueous NaBr solutions for different salt concentrations. b The EFG ACF for Na + in aqueous NaBr, NaF, and NaCl solutions for c = 1 m and T = 25°C. In a and b, the shaded regions indicate the standard error from multiple independent simulation runs. The effective EFG correlation time τc (c) and the EFG variance (d) as a function of c in NaBr, NaF, and NaCl solutions. e The resulting quadrupolar NMR relaxation rates for 23 Na + in aqueous NaBr, NaF, and NaCl solutions as a function of the salt concentration c. In c, d, and e, the standard error from multiple independent simulation runs is either explicitly shown or does not exceed the symbol size.