Molecular dynamics simulations of ¹H NMR relaxation in Gd³⁺–aqua

To model the Gd³⁺–aqua system, we use the AMOEBA polarizable force field to describe solvent water [55] and the ion [56]. The Gd³⁺ ion has an incomplete set of $4f$ orbitals, but this incomplete shell lies under the filled $5s$ and $5p$ orbitals. Thus ligand-field effects [57] are not expected to play a part in hydration and a spherical model of the ion ought to be adequate. However, polarization effects are expected to be important given the large charge on the cation.

Experimental NMR studies on Gd³⁺–aqua use concentrations of about 0.3 mM. This amounts to having one GdCl₃ molecule in a solvent bath in a cubic cell about 176 Å in size. Such large simulations are computationally rather demanding with the AMOEBA polarizable forcefield. Further, at such dilutions, the ions essentially “see” only the water around them, and since understanding the behavior of water around the ion is of first interest, here we study a single Gd³⁺ ion in a bath of 512 water molecules. (Note that within the Ewald formulation for electrostatic interactions, there is an implicit neutralizing background. This background does not impact the forces between the ion and the water molecules that are of interest here.) The partial molar volume of Gd³⁺ in water at 298.15 K has been estimated by Marcus (1985) to be $-59.6$ cc/mol. We use this to fix the length of the cubic simulation cell to $24.805$ Å. (From constant pressure simulations in a 2006 water system we find a Gd³⁺ partial molar volume of about -63 cc/mol, in good agreement with the value suggested by Marcus. However, in this work, we will use the value suggested by Marcus.)

All the simulations were performed using the OpenMM-7.5.1 package [58]. The van der Waals forces were switched to zero from 11 Å to 12 Å. The real space electrostatic interactions were cutoff at 9 Å and the long-range electrostatic interactions were accounted using the particle mesh Ewald method with a relative error tolerance of 10^-5 for the forces. In the polarization calculations the (mutually) induced dipoles were converged to 10^-5 Debyes.

The equations of motion were integrated using the “middle” leapfrog propagation algorithm with a time step of 1 fs coupled; this combination of method and time step provides excellent energy conservation in constant energy simulations. Exploratory work shows that the NMR relaxation in the Gd³⁺-aqua system is sensitive to the system temperature. So we additionally use with a Nosé-Hoover chain [59, 60] with three thermostats [61] to simulate the system at 298.15 K. The collision frequency of the thermostat was set to 100 fs to ensure canonical sampling. We carried out extensive tests with and without thermostats to ensure that the Nosé-Hoover thermostat does not affect dynamical properties. Our conclusions are consistent with an earlier study on good practices for calculating transport properties in simulations [62].

Initially, the system was equilibrated under $NVT$ conditions at 298.15 K for over 200 ps. We then propagated the system under the same $NVT$ conditions for 8 ns, saving frames every 0.1 ps for analysis. We used the last 6.5536 ns of simulation (equal to 65536($=2^{16}$) frames) for analysis. The mean temperature in the production phase was $298.15$ K with a standard error of the mean being 0.03 K.

Fig. 1 shows the ion-water radial distribution function. The location and magnitude of the peak is in good agreement with earlier studies founded on either ab initio [29] or empirical force field-based [31] simulations. Consistent with those studies, we find that the first sphere (i.e. inner sphere) contains between $q=8\leftrightarrow 9$ water molecules, with a mean of $q\simeq$ 8.5 consistent with published values [15, 17, 22, 21, 24, 31, 32].

To estimate the residence time of water molecules in the inner sphere, we need to keep track of the water molecule as it moves into/out of the inner sphere. To this end, we follow the residence time of a defined water molecule $w$ using an indicator function $\chi_{w}$ that equals 1 if the water molecule is in the inner sphere and zero otherwise. The inner sphere is defined as a sphere of radius $r\leq 3.3$ Å around the Gd³⁺ ion; this corresponds to the first hydration of the ion (Fig. 1). We perform this analysis for all the water molecules that visit the inner sphere at least once during the simulation. (We emphasize that the time here is discrete because configurations are saved only every 100 fs.)

Fig. 2 shows the trace of the indicator function for a particular water molecule.

Note that during the approximately 4 ns window (from the $\approx 8$ ns trajectory), the water molecule makes several excursions out of the inner sphere before permanently leaving the inner sphere around 4.5 ns. The length of time that the water molecule spends continuously inside the inner sphere is denoted by $\tau_{i}$. As Fig. 2 shows, there can be several such islands of continuous occupation.

To test if the transient excursion is a bona fide escape from the inner sphere, we window average the data as follows: we consider a transient escape as a bona fide escape only if it persists for a defined number of consecutive time points. Fig. 2 shows the trace (blue curve) for such a “windowing” for a window length of 200 fs, i.e. two consecutive frames in the trajectory of configurations. As can be expected, windowing extends the length of time that the molecule is defined to be inside the inner sphere.

For all the water molecules that visit the inner sphere, we accumulate the set $\{\tau_{i}\}$ and then construct the histogram $h(\tau)$. Figure 3 shows the $h(\tau)$ for the raw data and the data with window length of 200 fs.

The $h(\tau)$ curve shows that there are many cases where water only transiently resides in the inner sphere. But we also find several water molecules that spend upwards of 0.5 ns continuously within the inner sphere. Specifically, for the data that has not been smoothed by “windowing”, we find three instances of water molecules continuously spending between 1 to 1.5 ns inside the inner sphere (and these also happen to be three distinct water molecules [data not shown]); with windowing using a 200 fs window, the upper limit is extended to just over 2 ns. Thus residence time, $\tau_{m}$, between $1\leftrightarrow$ 2 ns are predicted for a complete rejuvenation of the inner sphere population. This time-scale is in accord with the range of published ¹⁷O NMR data that suggest residence times $\tau_{m}\simeq 1.0\leftrightarrow 1.5$ ns [16, 63, 20, 64, 26].

The enhanced ¹H NMR relaxation rate $1/T_{1}$ for water is given by the average over the $N=1024$ ¹H nuclei in the $L=24.805$ Å box containing one paramagnetic Gd³⁺ ion:

	$\displaystyle\frac{1}{T_{1}}$	$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\frac{1}{T_{1i}},$		(1)
	$\displaystyle\frac{1}{T_{2}}$	$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\frac{1}{T_{2i}},$		(2)

where $T_{1i}$ is the $T_{1}$ relaxation for the $i$’th ¹H nucleus. The gadolinium molar concentration is given by $[Gd]=[H]/N$, where $[H]$ is the molar concentration of ¹H in the simulation box. Equivalently, $[H]=2[W]$ where $[W]=$ 55,705 mM is the molar concentration of H₂O at 25^∘C. This leads to the following expression for the relaxivity in units of (mM^-1s^-1):

	$\displaystyle r_{1}$	$\displaystyle=\frac{1}{[Gd]}\frac{1}{T_{1}}=\frac{1}{[H]}\sum_{i=1}^{N}\frac{1}{T_{1i}},$		(3)
	$\displaystyle r_{2}$	$\displaystyle=\frac{1}{[Gd]}\frac{1}{T_{2}}=\frac{1}{[H]}\sum_{i=1}^{N}\frac{1}{T_{2i}}.$		(4)

Note that $r_{1,2}$ are independent of $N$ (or box size $L$, equivalently). The “fast-exchange” regime ($T_{1,2}\gg\tau_{m}$) is assumed [11], and the chemical shift term in $r_{2}$ [10] is neglected for simplicity. Note also that the ¹H–¹H dipole-dipole relaxation [48] is not considered in these simulations.

The computation of $T_{1i,2i}$ originates from the Gd³⁺–¹H dipole-dipole autocorrelation function $G(t)$ [45, 5, 65, 66, 67, 68] shown in Fig. 4(a), where $t$ is the lag time. This autocorrelation is well suited for computation using MD simulations [69, 70]. Using the convention in the text by McConnell [66], $G_{i}(t)$ in units of (s^-2) is determined by:

$$G_{i}(t)=\frac{1}{4}\!\left(\frac{\mu_{0}}{4\pi}\right)^{2}\!\hbar^{2}\gamma_{I}^{2}\gamma_{S}^{2}S(S+1)\times\\ \left<\frac{(3\cos^{2}\!\theta_{i}\!(t+t^{\prime})-1)}{r_{i}^{3}\!\left(t+t^{\prime}\right)}\frac{(3\cos^{2}\!\theta_{i}\!(t^{\prime})-1)}{r_{i}^{3}\!(t^{\prime})}\right>_{\!\!t^{\prime}},$$

(5)

for the $i$’th ¹H nucleus. $\theta_{i}$ is the angle between the Gd³⁺–¹H vector ${\bf r}_{i}$ and the applied magnetic field ${\bf B}_{0}$. $\mu_{0}$ is the vacuum permeability, $\hbar$ is the reduced Planck constant. $\gamma_{I}/2\pi=42.58$ MHz/T is the nuclear gyro-magnetic ratio for ¹H with spin $I=1/2$, and $\gamma_{S}=658\,\gamma_{I}$ is the electron gyro-magnetic ratio for Gd³⁺ with spin $S=7/2$.

Note that in Eq. 5 we assume a spherically symmetric (i.e. isotropic) system, and therefore $G_{i}^{m}(t)=G_{i}(t)$ is independent of the order $m$, which amounts to saying that the direction of the applied magnetic field ${\bf B}_{0}=B_{0}\bf z$ is arbitrary. This assumption was verified in Ref. [31]. For simplicity, we therefore use the $m=0$ harmonic $Y_{2}^{0}(\theta,\phi)=\sqrt{5/16\pi}\,(3\cos^{2}\theta-1)$ for the MD simulations [48].

The second-moment $\Delta\omega_{i}^{2}$ (i.e. strength) of the dipole-dipole interaction is defined as such [67]:

$\displaystyle G_{i}(0)$

$\displaystyle=\frac{1}{3}\Delta\omega_{i}^{2}.$

(6)

Assuming the angular term in Eq. 5 is uncorrelated with the distance term at $t=0$, the relation $\left<(3\cos^{2}\!\theta_{i}\!(\tau)-1)^{2}\right>_{\tau}=4/5$ (which is independent of $i$) reduces the second moment to:

$\displaystyle\Delta\omega_{i}^{2}=\frac{3}{5}\!\left(\frac{\mu_{0}}{4\pi}\right)^{2}\!\hbar^{2}\gamma_{I}^{2}\gamma_{S}^{2}S(S+1)\left<\!\frac{1}{r_{i}^{6}(t^{\prime})}\!\right>_{\!\!t^{\prime}}.$

(7)

The next step is to take the (two-sided even) fast Fourier transform (FFT) of the $G_{i}(t)$ to obtain the spectral density function:

$\displaystyle J_{i}(\omega)$

$\displaystyle=2\int_{0}^{\infty}\!G_{i}(t)\cos(\omega t)\,dt.$

(8)

The relaxation rates are then determined for unlike spins [66]:

	$\displaystyle\frac{1}{T_{1i}}$	$\displaystyle=J_{i}(\omega_{0})+\frac{7}{3}J_{i}(\omega_{e}),$		(9)
	$\displaystyle\frac{1}{T_{2i}}$	$\displaystyle=\frac{2}{3}J_{i}(0)+\frac{1}{2}J_{i}(\omega_{0})+\frac{13}{6}J_{i}(\omega_{e}),$		(10)

assuming $\omega_{e}\gg\omega_{0}$, where $\omega_{0}=\gamma_{I}B_{0}=2\pi f_{0}$ is the ¹H NMR resonance frequency, and $\omega_{e}=658\,\omega_{0}$ is the electron resonance frequency.

The expressions for $T_{1i,2i}$ are then summed in Eqs. 3 and 4 to compute $r_{1}$ and $r_{2}$, respectively. We also define the following quantities summed over the $N=1024$ ¹H nuclei:

	$\displaystyle G(t)$	$\displaystyle=\sum\limits_{i=1}^{N}\!G_{i}(t),$		(11)
	$\displaystyle\Delta\omega^{2}$	$\displaystyle=\sum\limits_{i=1}^{N}\Delta\omega_{i}^{2}$		(12)
	$\displaystyle J(\omega)$	$\displaystyle=\sum\limits_{i=1}^{N}J_{i}(\omega).$		(13)

Note that the summed quantities $G(t)$, $\Delta\omega^{2}$, and $J(\omega)$ are independent of $N$ (or box size $L$, equivalently). The $G(t)$ simulation data is plotted in Fig 4(a), while the $J(\omega)$ simulation data is plotted in Fig 4(b).

We also define the average correlation time $\left<\tau\right>$ as the normalized integral [67]:

$\displaystyle\left<\tau\right>$

$\displaystyle=\frac{1}{G(0)}\int_{0}^{\infty}\!G(t)\,dt.$

(14)

The low frequency (i.e. extreme narrowing) limit $r_{1}(0)=r_{2}(0)=r_{1,2}(0)$ can then be expressed as:

$\displaystyle r_{1,2}(0)=\frac{1}{[H]}\frac{20}{9}\Delta\omega^{2}\!\left<\tau\right>.$

(15)

Note how $r_{1,2}(0)$ (for unlike Gd³⁺–¹H spin pairs) is a factor 2/3 less than the equivalent expression for like ¹H–¹H spin pairs [48].

We assume the “fast-exchange” regime in the above formulation, which is discussed in more detail in Section III.2. The fast-exchange regime can be inferred directly from measurements since $r_{1}$ increases with decreasing temperature [15]. Investigations are underway to extend the simulations to the slow-exchange regime.

The above analysis also assumes the electron spin is a point-dipole centered at the Gd³⁺ ion [4]. Given that the simulation agrees with measurements in the range $f_{0}\simeq$ 5 $\leftrightarrow$ 500 MHz, the point-dipole approximation is considered valid in the present case for ¹H [10].

The FFT result for $J(\omega)$ in Fig. 4(b) is sparse. Besides the $f$ = 0 data point, the lowest frequency FFT data point is given by the resolution $\Delta f=1/2t_{max}=$ 10³ MHz, where $t_{max}=500$ ps is the longest lag time in $G(t)$.

As an alternative to zero-padding the FFT, we see to model $G(t)$ in terms of molecular modes. To this end, we expand $G(t)$ as

$\displaystyle G(t)$

$\displaystyle=\int_{0}^{\infty}\!P(\tau)\exp\left(-\frac{t}{\tau}\right)d\tau,$

(16)

where $P(\tau)$ is the underlying distribution in molecular correlation times, $\tau$. We solve this Fredholm integral equation of the first kind to recover the $P(\tau)$ distribution. Since $G(t)$ is available only at discrete time intervals, the inversion is an ill-posed problem. We address this by using Tikhonov regularization [50, 71], with the vector P being one for which

$\displaystyle{\bf P}=\underset{{\bf P}\geq 0}{\mathrm{arg\,min}}\,\,||{\bf G}-K\,{\bf P}||^{2}+\alpha||{\bf P}||^{2}$

(17)

is a minimum. Here G is the column vector representation of the autocorrelation function $G(t)$, P is the column vector representation of the distribution function $P(\tau)$, $\alpha$ is the regularization parameter, and $K$ is the kernel matrix:

$\displaystyle K=K_{ij}=\exp\left(-\frac{t_{i}}{\tau_{j}}\right).$

(18)

The results for $P(\tau)$ are shown in Fig. 5, from which the following are determined:

	$\displaystyle\left<\tau\right>$	$\displaystyle=\frac{1}{G(0)}\int_{0}^{\infty}\!P(\tau)\,\tau\,d\tau,$		(19)
	$\displaystyle G(0)$	$\displaystyle=\int_{0}^{\infty}\!P(\tau)d\tau=\frac{1}{3}\Delta\omega^{2}.$		(20)

The spectral density $J(\omega)$ is then determined from the Fourier transform (Eq. 8) of $G(t)$ (Eq. 16):

$\displaystyle J(\omega)$

$\displaystyle=\int_{0}^{\infty}\!\frac{2\tau}{1+(\omega\tau)^{2}}P(\tau)d\tau\,,$

(21)

from which $T_{1,2}$ at any desired $f_{0}$ can be determined (Eqs. 9 and 10). More in-depth discussions of the above procedure, loosely termed as “inverse Laplace transform”, can be found in Refs. [54, 53, 52, 51, 72, 73] and the supplementary material in Refs. [50, 71]. We hasten to add that inverting Eq. 16 to recover $P$ is formally not a Laplace inversion [74], but this terminology is common in the literature. Possible alternatives to the above formulation are also discussed in Ref. [51].

The $P(\tau)$ is binned from $\tau_{min}=0.05$ ps to $\tau_{max}=500$ ps using 200 logarithmically spaced bins. In the present case of low-viscosity fluids ($\eta\simeq$ 1 cP), the choice of $\tau_{max}$ does not impact $J(\omega)$, and is therefore not a free parameter in the analysis in terms of molecular modes. The constant “div” in Fig. 5 is a “division” on a log-scale. More specifically, div $=\log_{10}(\tau_{i+1})-\log_{10}(\tau_{i})$ is independent of the bin index $i$, and ensures unit area when $P(\tau)$ is of unit height and a decade wide [71].

As shown in Fig. 4(a), the residual between the $G(t)$ data and the fit using molecular modes is not dominated by Gaussian noise. As such, the regularization parameter is fixed to $\alpha=$ 10^-1 in accordance with previous studies [50, 52, 53, 54, 51]. As shown in Fig. 4(b), we find that $\alpha=$ 10^-1 gives excellent agreement with the parameter-free $J(\omega)$ from FFT, which validates the analysis in terms of molecular modes. The results in Fig. 4(b) further emphasize the following advantages: (1) the expansion (Eq. 16) filters out the noise while still honoring the FFT data (including the $f=0$ data point), (2) Eq. 21 provides $J(\omega)$ for any desired $f$ value, and (3) the expansion in terms of molecular modes leads to physical insight into the distribution $P(\tau)$ of molecular correlation times $\tau$.

An independent computation of translational diffusion $D_{T}$ was performed from MD simulations. We calculate the mean square displacement $\left<\!\Delta r^{2}\right>$ of the water oxygen and Gd³⁺ ion as a function of the diffusion evolution time $t$ ($<$ 10 ps), and additionally average over a sample of 50 molecules to ensure adequate statistical convergence.

As shown in Fig.  6, at long-times ($t$) the slope of the linear diffusive regime gives the translational self-diffusion coefficient $D_{sim}$ according to Einstein’s relation:

$$D_{sim}=\frac{1}{6}\frac{\delta\!\left<\!\Delta r^{2}\right>}{\delta t}\,$$

(22)

where $D_{sim}$ is the diffusion coefficient obtained in the simulation using periodic boundary conditions in a cubic cell of length $L$. In the linear regression procedure, the early ballistic regime and part of the linear regime is excluded to obtain a robust estimate of $D_{sim}$.

Following Yeh and Hummer [75] (see also Dünweg and Kremer [76]), we obtain the diffusion coefficient for an infinite system $D_{T}$ from $D_{sim}$ using

$$D_{T}=D_{sim}+\frac{k_{B}T}{6\pi\eta}\frac{\xi}{L}$$

(23)

where $\eta$ is the shear viscosity and $\xi=2.837297$ [75] is the same quantity that arises in the calculation of the Madelung constant in electrostatics. (In the electrostatic analog of the hydrodynamic problem, $\xi/L$ is the potential at the charge site in a Wigner lattice.) For the system sizes considered in this study, $L\simeq 25$ Å, the correction factor constituted $\simeq 13\%$ of $D_{0}$, and $\simeq 16\%$ of $D_{W}$. The correction factor was not applied to $D_{Gd}$.

We prepared a Gd³⁺–aqua solution in de-ionized water at $[Gd]$ = 0.3 mM and measured $T_{1,meas}$ at a controlled temperature of 25^∘C, using static fields at $f_{0}=2.3$ MHz with an Oxford Instruments GeoSpec2, at $f_{0}=20$ MHz with a Bruker Minispec, and at $f_{0}=500$ MHz with a Bruker 500 MHz Spectrometer. The measured relaxivity $r_{1}$ was determined as such:

$\displaystyle r_{1}=\frac{1}{[Gd]}\left(\!\frac{1}{T_{1,meas}}-\frac{1}{T_{1,bulk}}\!\right).$

(24)

where $T_{1,bulk}=$ 3.13 s was found for bulk water (not de-oxygenated [77]) at 2.3 MHz and 500 MHz. Field cycling $r_{1}$ data at $[Gd]$ = 1 mM and 25^∘C were taken from Luchinat et al. [32] (supplementary material) using a Stelar SpinMaster 1T. The field cycling results agreed with our measurements at $f_{0}=$ 2.3 MHz and 20 MHz (within $\simeq$ 5%), while our $f_{0}=500$ MHz significantly extends the frequency range of the measurements.

The results of simulated relaxivity $r_{1}$ (Eq. 3) are shown in Fig. 7, alongside corresponding measurements of Gd³⁺–aqua solution at 25^∘C using field cycling by Luchinat et al. [32] (supplementary material), and static fields from this work. A cross-plot of simulated versus measured $r_{1}$ results are also shown in Fig. 8. The simulation is within $\simeq$ 8% of measurements in the range $f_{0}\simeq$ 5 $\leftrightarrow$ 500 MHz. Given that there are no adjustable parameters in the interpretation of the simulations, this agreement validates our simulations at high frequencies. We note that a similar degree of agreement (within $\simeq$ 7%) was found in previous studies of liquid alkanes and water [48].

For convenience, Table 1 lists the average correlation time $\left<\tau\right>$ (Eq. 19), the square-root of the second moment (i.e. strength) of $\Delta\omega$ (Eq. 20), and the residence time $\tau_{m}$ (Fig. 3). Note that these three quantities are model free.

The SBM inner-sphere model assumes a rigid Gd³⁺–¹H pair undergoing rotational diffusion, leading to the following mono-exponential decay in the autocorrelation function:

$\displaystyle G_{SBM}(t)=G_{SBM}(0)\exp\left(-\frac{t}{\tau_{R}}\right).$

(25)

This functional form is identical to the BPP model [45] which is based on the Debye model, where the rotational-diffusion correlation time $\tau_{R}$ is defined as the average time it takes the rigid pair to rotate by 1 radian. $G_{SBM}(t)$ is plotted in Fig. 4(a) assuming $\tau_{R}=\left<\tau\right>$ = 30 ps.

As shown in Fig. 4(a), the mono-exponential decay in $G_{SBM}(t)$ is not consistent with the multi-exponential (i.e., stretched) decay in $G(t)$. Equivalently, $J(\omega)$ in Fig. 4(b) does not follow the $f^{-2}$ power-law behavior at large $f$. This is expected since the simulations implicitly include both inner-sphere and outer-sphere (see below) contributions. Currently the simulations do not separate inner-sphere from outer-sphere contributions, therefore the simulations do not clarify the origin of the multi-exponential decay in $G(t)$.

It is also possible that the inner-sphere $G(t)$ itself has a multi-exponential decay of its own. The functional form of $G_{SBM}(t)$ is based on the Debye model, which was previously shown to be inaccurate when used in the BPP model for liquid alkanes and water [48]. It would therefore not be surprising if the inner-sphere $G(t)$ was also multi-exponential in nature.

Assuming that inner-sphere relaxation dominates, and therefore that the multi-exponential decay in $G(t)$ is due to inner-sphere dynamics alone, the correlation time $\left<\tau\right>\simeq$ 30 ps is consistent with published values from the SBM inner-sphere model, where a range of $\tau_{R}\simeq 23\leftrightarrow 45$ ps is reported [15, 17, 22, 21, 24, 31, 32]. Note that $\left<\tau\right>\simeq$ 30 ps is a factor $\simeq 10$ larger than that of bulk water ($\tau_{R}\simeq 2.7$ ps [48]), which is expected given the hindered rotational motion of the Gd³⁺–aqua complex.

Continuing with the assumption that inner-sphere relaxation dominates, one can approximate the following inner-sphere quantities:

	$\displaystyle\frac{1}{r_{in}^{6}}$	$\displaystyle\simeq\frac{1}{2q}\sum_{i=1}^{N}\left<\!\frac{1}{r_{i}^{6}(t^{\prime})}\!\right>_{\!\!t^{\prime}}$		(26)
	$\displaystyle\Delta\omega^{2}_{in}$	$\displaystyle\simeq\frac{1}{2q}\sum_{i=1}^{N}\Delta\omega_{i}^{2}$		(27)
	$\displaystyle\frac{1}{T_{1in}}$	$\displaystyle\simeq\frac{1}{2q}\sum_{i=1}^{N}\frac{1}{T_{1i}}.$		(28)

The expressions are averaged over the $2q$ ¹H nuclei in the inner sphere, where $q=8.5$ is the H₂O inner-sphere coordination number determined from $n(r)$ (Fig. 1) at $r\simeq 3.5$ Å. Note again that these equations neglect the outer-sphere contribution, implying that $\Delta\omega_{in}^{2}$ is an upper bound, while $r_{in}$ and $T_{1in}(0)$ are lower bounds.

According to Table 1, the resulting $r_{in}\simeq$ 2.97 Å is consistent with published values from the inner-sphere SBM model, where a range of $r_{in}\simeq 3.0\leftrightarrow 3.2$ Å is reported [15, 17, 22, 21, 24, 31, 32]. The product $\Delta\omega_{in}\left<\tau\right>\simeq 0.002$ indicates that the Redfield-Bloch-Wagness condition ($\Delta\omega\left<\tau\right>\ll 1$) is satisfied [65, 66, 67], which justifies the relaxivity analysis used here. The fast-exchange regime [11] can also be verified by noting that $T_{1in}(0)\simeq$ 4.4 $\mu$s at $f_{0}=0$ is three orders of magnitude larger than $\tau_{m}\simeq 1$ ns, i.e. $(T_{1in}+\tau_{m})\simeq T_{1in}$ can be assumed.

We now discuss the outer-sphere contribution to relaxivity, although it is generally believed (though not proven) to be smaller than inner-sphere relaxivity [15, 18, 23]. The outer-sphere contribution is expected to follow the Hwang-Freed (HF) model for the relative translational diffusion between Gd³⁺ and H₂O assuming two force-free hard-spheres [6]:

$$G_{HF}(t)=G_{HF}(0)\frac{54}{\pi}\int\limits_{0}^{\infty}\frac{x^{2}}{81+9x^{2}-2x^{4}+x^{6}}\times\\ \exp\left(-x^{2}\frac{t}{\tau_{D}}\right)dx.$$

(30)

The translational-diffusion correlation time $\tau_{D}$ is defined as the average time it takes the molecule to diffuse by one hard-core diameter $d$. $G_{HF}(t)$ is a multi-exponential decay by nature, and therefore one expects the total $G(t)$ to be stretched, the extent of which depends on the relative contributions of outer-sphere to inner-sphere.

The correlation time $\tau_{D}$ can be predicted as such:

$\displaystyle\tau_{D}$

$\displaystyle=\frac{d^{2}}{D_{W}+D_{Gd}}=\frac{9}{4}\tau_{T},$

(31)

where the simulated diffusion coefficients of H₂O ($D_{W}$) and Gd³⁺ ($D_{Gd}$) in Gd³⁺–aqua are taken from Fig. 6, the results of which are listed in Table 2. The hard-core diameter $d$ is taken from the local maximum at $r\simeq 4.6$ Å in the pair-correlation function $g(r)$ in Fig. 1 (which is attributed to outer-sphere water). Note that the resulting $\tau_{D}\simeq 99$ ps is a factor $\simeq 10$ larger than that of bulk water ($\tau_{D}\simeq 9.0$ ps [48]), which is expected given the larger hard-core distance $d\simeq 4.6$ Å than bulk water ($d\simeq$ 2.0 Å [48]).

The value $\tau_{D}\simeq 99$ ps is compared to the distribution $P(\tau)$ in Fig. 5. More specifically, the translational correlation time $\tau_{T}\simeq 44$ ps is plotted in Fig. 5, where the relation $\tau_{D}=9/4\,\tau_{T}$ and the origin of the factor 9/4 is explained in [67, 50]. $\tau_{T}$ lies within the $P(\tau)$ distribution, indicating that the outer-sphere may contribute to relaxivity. Further investigations are underway to separate inner-sphere from outer-sphere contributions in the simulations, without assuming any models.

We note that $P(\tau)$ in Fig. 5 has a small contribution at short $\tau\simeq 10^{-1}$ ps, which is a result of the sharp drop in $G(t)$ over the initial $t\simeq 0.2$ ps (Fig. 4(a)). This molecular mode is also present for intramolecular relaxation in liquid alkanes and water, while it is absent for intermolecular relaxation. In the case of alkanes, the ubiquitous intramolecular mode at $\tau\simeq 10^{-1}$ ps is attributed to the fast spinning methyl end-groups [50]. Investigations are underway to better understand the origin of this mode in Gd³⁺–aqua, which can perhaps be explained using a two rotational-diffusion model such as found in bulk water [78]. The origin of other modes in $P(\tau)$ at $\tau\simeq 10^{0}$ ps and $\simeq 10^{1}$ ps are also being investigated.

Finally, we note that the $r_{1}$ dispersion in Fig. 7 results in a mild increase in the ratio $T_{1}/T_{2}=r_{2}/r_{1}\simeq 7/6$ above $f_{0}\gtrsim 10$ MHz, until $f_{0}\gtrsim 6$ GHz where $T_{1}/T_{2}$ increases further. Combining Eqs. 9 and 10 with $J(\omega_{e})=0$ (i.e., slow-motion regime) and $J(\omega_{0})=J(0)$ (i.e., fast-motion regime) accounts for the factor $T_{1}/T_{2}=$ 7/6 within the frequency range $f_{0}\simeq$ 10 MHz $\leftrightarrow$ 6 GHz. The ratio $T_{1}/T_{2}=7/6$ was also used to explain water-saturated sandstones [38].

At low-frequencies ($f_{0}\lesssim 0.5$ MHz), the difference between $r_{1}$ measurements ($\simeq$ 29.7 mM^-1s^-1) and simulation ($\simeq$ 34.6 mM^-1s^-1) can be reconciled by taking the zero-field electron-spin relaxation time $T_{e0}=T_{1e}(0)=T_{2e}(0)$ into account. Assuming that the correlation times $P(\tau)$ are uncorrelated with the electron-spin relaxation times, the following expression results [3]:

	$\displaystyle r_{1,2}^{\prime}(0)$	$\displaystyle=\frac{1}{[H]}\frac{20}{9}\Delta\omega^{2}\!\left<\tau^{\prime}\right>,$		(32)
	$\displaystyle\frac{1}{\left<\tau^{\prime}\right>}$	$\displaystyle=\frac{1}{\left<\tau\right>}+\frac{1}{T_{e0}}.$		(33)

This is equivalent to introducing an exponential decay term $\exp(-t/T_{e0})$ inside the FFT integral of Eq. 8. The fitted value of $T_{e0}\simeq 180$ ps is determined by matching $r_{1}^{\prime}(0)$ to the average low-frequency ($f_{0}\lesssim 0.5$ MHz) measurement (Fig. 7). The resulting $T_{e0}\simeq 180$ ps is consistent with the published range of $T_{e0}\simeq 96\leftrightarrow 160$ ps [15, 17, 22, 21, 24, 31, 32]. Investigations are underway to incorporate $T_{1e}(\omega_{e})$ and $T_{2e}(\omega_{e})$ dispersion [4, 24, 26, 31] for predicting $r_{1}$ at low frequencies.