NMR Spin-Rotation Relaxation and Diffusion of Methane
Philip M. Singer D. Asthagiri Walter G. Chapman George J. Hirasaki


J. Chem. Phys.



ps41@rice.edu
dna6@rice.edu




Department of Chemical and Biomolecular Engineering, Rice University, 6100 Main St., Houston, TX 77005, USA






	
The translational-diffusion coefficient D_T and the spin-rotation contribution to the ^1H NMR relaxation time T_1J for methane (CH_4) are investigated using MD (molecular dynamics) simulations, over a wide range of densities ρ and temperatures T, spanning the liquid, supercritical, and gas phases. The simulated D_T agree well with measurements, without any adjustable parameters in the interpretation of the simulations. A minimization technique is developed to compute the angular-velocity for non-rigid spherical molecules, which is used to simulate the autocorrelation function G_J(t) for spin-rotation interactions. With increasing D_T (i.e. decreasing ρ), G_J(t) shows increasing deviations from the single-exponential decay predicted by the Langevin theory for hard spheres, and the deviations are quantified using inverse Laplace transforms of G_J(t). T_1J is derived from G_J(t) using the kinetic model “km" for gases (T_1J^km), and the diffusion model “dm" for liquids (T_1J^dm). T_1J^km shows better agreement with T_1 measurements at higher D_T, while T_1J^dm shows better agreement with T_1 measurements at lower D_T. T_1J^km is shown to dominate over the MD simulated ^1H-^1H dipole-dipole relaxation T_1RT at high D_T, while the opposite is found at low D_T. At high D_T, the simulated spin-rotation correlation-time τ_J agrees with the kinetic collision time τ_K for gases, from which a new relation 1/T_1J^km∝ D_T is inferred, without any adjustable parameters.



NMR Spin-Rotation Relaxation and Diffusion of Methane
    G. J. Hirasaki
    December 30, 2023
=====================================================





§ INTRODUCTION


Theoretical and experimental investigations into ^1H NMR (nuclear magnetic resonance) spin-rotation relaxation T_1J and translational-diffusion D_T of methane (CH_4), deutero-derivatives (CH_4-nD_n), and halide-derivatives (CH_4-n-mF_nCl_m) dates back over 50 years <cit.>. Likewise, the theoretical and experimental investigations into molecular-beam magnetic-resonance of methane and its derivatives dates back over 50 years <cit.>, from which came (among many other things) the coupling constants used to interpret T_1J.

More recently, the influence of dissolved methane on the T_1 and D_T of complex crude-oils and hydrocarbon mixtures has been investigated <cit.>, which reveal the important influence of the spin-rotation contribution from methane. Of particular significance is the well established spin-rotation component for methane in the fast-motion regime <cit.>:

    1/T_1J^A = A T^3/2/ρ,

where A is an empirically derived constant. Eq. <ref> states that the relaxation rate 1/T_1J increases with temperature T, i.e. T_1J decreases with T. This is in stark contrast to the ^1H-^1H dipole-dipole relaxation T_1RT which increases with T, for all hydrocarbons, including methane. Also of interest recently is the influence of pore confinement on the NMR response of methane <cit.>, which has practical applications for characterizing the light hydrocarbons in the organic nano-pores of kerogen and bitumen in organic-rich shale. One of the current mysteries is why the T_1S/T_2S ratio for surface-relaxation of methane in organic-shale is typically T_1S/T_2S≲ 2, while for higher-order alkanes it is typically higher T_1S/T_2S≳ 4. This has practical applications for separating the NMR response of light hydrocarbon from water in organic-shale, and for determining the hydrocarbon saturation in the organic-shale reservoir.

In order to properly characterize NMR relaxation of methane in the bulk and under nano-pore confinement, one must first separate the ^1H spin-rotation relaxation T_1J from ^1H-^1H dipole-dipole relaxation T_1RT. Traditionally this has been done by partially deuterating CH_4 to, for instance CHD_3, which dramatically reduces the dipole-dipole contribution. However, this has the drawback of turning a spherical molecule (CH_4) with one principle moment of inertia, into a symmetric-top molecule (CHD_3) with two principle moments of inertia. The theory of spin-rotation relaxation for symmetric-top molecules is much more complex than for spherical molecules, thereby making comparisons with measurements more complex.

MD (molecular dynamics) simulations provide an ideal tool for separating ^1H NMR relaxation mechanisms. As already shown for liquid-state n-alkanes in Ref. <cit.>, MD simulations can naturally separate intramolecular T_1R from intermolecular T_1T ^1H-^1H dipole-dipole relaxation, without deuteration, and without any adjustable parameters in the interpretation of the simulations. Such simulations yield unique insights into (a) the relative strengths of intramolecular versus intermolecular relaxation, (b) the influence of internal motions on the molecular dynamics of non-rigid molecules, and (c) the validity of traditional hard-sphere models <cit.> for different molecular geometries. 

In this report we simulate spin-rotation relaxation. In Section <ref>, we develop a technique to determine the autocorrelation function for angular-velocity G_J(t) of non-rigid spherical molecules. In Section <ref> we interpret G_J(t) to yield the spin-rotation relaxation using the kinetic model and the diffusion model. In Section <ref> we compare simulation versus measurement for spin-rotation relaxation and translational diffusion. In Section <ref> we compare and analyze the correlation times for the different relaxation mechanisms. In Section <ref> we propose a new kinetic model to account for Eq. <ref>.







§ METHODOLOGY
 



 §.§ Molecular simulation
 

The MD simulations were performed using NAMD <cit.> version 2.11. Methane was modeled using the CHARMM General Force field (CGenFF) <cit.>. The protocol for setting-up the initial simulation configuration was exactly as before <cit.>. As before, we created the initial simulation system by packing N copies of the molecule into a cube of volume L^3 using the Packmol program <cit.>. The volume was chosen such that the number density N/V corresponds to the experimentally determined number density at the specified temperatures listed in Table <ref>. The simulation approach for these systems using NAMD was as before <cit.>. 

The angular-velocity computation went as follows. From simulations, for each atom in a molecule, we have a site velocity v_i. Without loss of generality, we assume v_i is relative to the velocity of the center of mass. By the definition of angular velocity, we have 

    ω×r_i = v_i.

But a direct application of the above equation cannot be used to calculate ω because the matrix equation is singular. We therefore define 

    ℒ = ∑_i | v_i - ω×r_i|^2.

We minimize ℒ with respect to the x, y, z-components of ω, giving the following matrix equation:

    [ ∑_i z_i^2 + y_i^2      -∑_i x_i y_i        -∑ x_i z_i;      -∑_i x_i y_i ∑_i x_i^2 + z_i^2        -∑ y_i z_i;      -∑_i x_i z_i        -∑ y_i z_i ∑_i x_i^2 + y_i^2 ][ ω_x; ω_y; ω_z ]
    = 
    [ ∑_i (r_i ×v_i)_x; ∑_i (r_i ×v_i)_y; ∑_i (r_i ×v_i)_z ]

In tensor notation, we can write this more compactly as

    r_j ω_p r_q ϵ^ljkϵ_pqk = r_j v_k ϵ_jkl,

summed over all sites. 

We constructed the above matrix for each molecule and solved for ω_x, ω_y, ω_z. This calculation is repeated for all the frames in the trajectory. Then we calculated the autocorrelation of each component. Since the system is isotropic, all the component relaxations are the same. The above calculation is repeated for other molecules and the results averaged.




The translational-diffusion coefficient D_T was derived in a similar fashion to Ref. <cit.>. The linear slope of the mean-squared displacement was computed in the interval between 5 ps ↔ 10 ps. The periodic-boundary correction term <cit.> was derived using viscosity η and box size L, both listed in Table <ref>. The correction term resulted in a 5 ↔ 10 % boost in diffusion coefficient.



 §.§ Spin-rotation relaxation
 

The Hamiltonian for the spin-rotation interaction is given by the following <cit.>:

    ℋ_J = ħ∑_i=1^4  I_i· C_i· J.

I_i are the four ^1H nuclear spins on the methane molecule. J is the angular momentum of the molecule, which is related to the angular velocity 
ω by the following ħ J = I ω, where I = 5.33 × 10^-47 kg m^2 <cit.> is the moment of inertia for methane. The coupling tensor C has two principle components C_⊥ and C_∥, which are proportional to the magnetic field (per unit J) generated at a ^1H by rotations about an axis ⊥ (∥) to the C-H bond axis, respectively. The most practical formulation separates ℋ_J into a scalar component proportional to the average coupling constant C_a = (2C_⊥ + C_∥)/3, and a tensor component proportional to the diagonal coupling constant C_d = (C_⊥ - C_∥). The accepted values of the coupling constants, and the ones used here, are <cit.>: C_a/2π = 10.4 ± 0.1   kHz, C_d/2π = 18.5 ± 0.5   kHz (which are the experimentally determined quantities), from which the following can be inferred C_⊥/2π = 16.6 ± 0.3  kHz, C_∥/2π = -1.9 ± 0.2  kHz.

The expression for the autocorrelation function G_J^lk(t) used to determine T_1J for liquids <cit.> and gases <cit.> can be summarized as follows (in units of s^-2):

    G_J^lk(t) = I^2/2 ħ^2(C_a^2 < ω_l(t+τ)  ω_k(τ)  >_τ  + 
    29C_d^2 ∑_m,n,o,pα_l,m,n^k,o,p< Y_2^n(t+τ)  ω_m(t+τ)   Y_2^p(τ)  ω_o(τ)  >_τ).

The first term in Eq. <ref> is the scalar term, while the second term is the tensor term. The more general expression also includes a cross-term proportional to C_a C_d, however it reduces to zero for spherically-symmetric molecules such as methane. The terms ω_k(t) is the angular velocity around the k axis at time t, and Y_2^m(t) is the spherical harmonic of rank 2 and order m at time t. α_l,m,n^k,o,p is composed of two sets of 3-j symbols defined in <cit.>. 





The tensor term greatly simplifies provided the rotational motion Y_2^m(t) is independent of the angular velocity ω_k(t). In such cases the bracketed term separates as < Y_2^n(t+τ)    Y_2^p(τ)>_τ<ω_m(t+τ)  ω_o(τ) >_τ. Traditional theories state the following expressions for hard spheres:

    < ω_l(t+τ) ω_k(τ) >_τ   =  δ_l,k k T/Iexp(-t/τ_J),
    < Y_2^l*(t+τ)   Y_2^k(τ) >_τ   = δ_l,kexp(-t/τ_R).

The angular-velocity autocorrelation function in Eq. <ref> is derived from the Langevin model with “friction time" τ_J <cit.>. The orientation autocorrelation function in Eq. <ref> is the Deybe model with rotational correlation-time τ_R <cit.>. The orientation and angular-velocity autocorrelations can be separated provided either τ_J ≫τ_R (gases) or τ_J ≪τ_R (liquids) <cit.>. In the case of liquids (τ_J ≪τ_R), this has the effect of making < Y_2^l*(t+τ)   Y_2^k(τ) >_τ = δ_l,k in Eq. <ref>. Based on comparison with measurements, the relation < Y_2^l*(t+τ)   Y_2^k(τ) >_τ = δ_l,k is also used in the case of gases <cit.>. In such cases, Eq. <ref> simplifies to the following expression which is used in the MD simulations:

    G_J(t) = I^2/2 ħ^2(C_a^2 + α29C_d^2)  < ω_k(t+τ)  ω_k(τ)  >_τ

The resulting expression G_J(t) is independent of the direction k in G_J^kk(t), therefore the superscript is removed for clarity. The MD simulations compute all three k directions independently, and the average is then taken to improve the signal to noise ratio. The simulation results for G_J(t) are shown in Fig. <ref>(a), where both x and y axes have been normalized for better comparison of the functional form of the decay between the different states. The normalization also allows for comparison with the Langevin model in Eq. <ref>.


A significant parameter in the analysis is the autocorrelation
at t = 0, which is given by the following expression:

    G_J(0) = I^2/2 ħ^2(C_a^2 +  α29C_d^2)  <  ω^2_k(τ)  >_τ= 
    I kT/2 ħ^2(C_a^2 +  α29C_d^2) = 1/2Δω_α^2,

where the first equality is directly from Eq. <ref>. The second equality uses the time zero expression <  ω^2_k(τ) >_τ  = kT/I from Eq. <ref>, which was verified from MD simulations to be within ±1.5 % over the entire temperature range of interest. The third equality defines the second moment Δω_α^2 (i.e. strength) of the spin-rotation interaction <cit.>, given by:

    Δω_km^2    = I kT/ħ^2(C_a^2 +  445C_d^2)   [ k.m., α = 25], 
    Δω_dm^2    = I kT/ħ^2(C_a^2 +  29C_d^2)   [ d.m., α = 1].

The only free parameter is α, which is given by α = 1 in the diffusion model “dm" for liquids <cit.>, or by α = 2/5 in the kinetic model “km" for gases <cit.>. In describing the kinetic model, we have adopted the exact relation for the scalar term ħ^2<J(J+1)>/3= I kT, and for the tensor term we assume (without loss in accuracy) that <(2J-1)(2J+3)>/4<J(J+1)> ≃ 1 in the classical limit J ≫ 1 <cit.>. The reason for α = 2/5 in the kinetic model for gases is that the oscillatory terms Δ J ≠ 0 do not contribute (giving rise to × 1/5), and from the statistical independence of ω_k(t) and Y_2^m(t) (giving rise to × 2) <cit.>. 

The next quantity of interest is the spin-rotation correlation time τ_J determined from the MD simulations, which is determined from the integral of the normalized G_J(t) as such <cit.>:

    τ_J = 1/G_J(0)∫_0^∞G_J(t)  dt.

The NMR relaxation times are then derived from the spectral density J_J^α(ω) <cit.>:

    J_J^α(ω)    = 2∫_0^∞G_J(t)cos(ω t) dt, 
    1/T_1J^α   = 2 J_J^α(ω_0), 
    1/T_2J^α   = J_J^α(0) + J_J^α(ω_0),

where ω_0 = γ B_0 is the Larmor frequency for ^1H. Given the short correlation times τ_J ∼ 1 ps, and given typical Larmor frequencies ω_0/2π < 500 MHz, it is clearly the case that ω_0 τ_J≪ 1, i.e. the fast-motion regime applies. In such cases J_J^α(0) = J_J^α(ω_0)  = Δω_α^2 τ_J, and therefore T_1J^α = T_2J^α. Using the relation 1/T_1J^α = 2J_J^α(0) = 2 Δω_α^2  τ_J results in the final expressions:

    1/T_1J^km   = 2I kT/ħ^2(C_a^2 +  445C_d^2) τ_J, 
    1/T_1J^dm   = 2I kT/ħ^2(C_a^2 +  29C_d^2) τ_J.

Note that T_1J^dm is often expressed in terms of C_⊥ and C_∥ instead <cit.>, where 
(C_a^2 +  29C_d^2) = 13(2C_⊥^2 + C_∥^2).

In order to compare with measurements, we also include contributions from the ^1H-^1H dipole-dipole interactions, which separate into intramolecular T_1R and intermolecular T_1T relaxation. Details of the methodology behind the MD simulations of T_1R and T_1T can be found in Ref. <cit.>. The final expression for the total relaxation time is given by:

    1/T_1^α = 1/T_1R + 1/T_1T+ 1/T_1J^α = 
    10/3Δω_R^2  τ_R +  10/3Δω_T^2 τ_T +  2 Δω_α^2 τ_J,

which takes on two different values T_1^km and T_1^dm, depending on the model T_1J^km and T_1J^dm used. Δω_R,T^2 and τ_R,T are the second-moments and correlation times for intramolecular (R) and intermolecular (T), respectively. For comparison purposes, the square-root of the second-moments at 298 K and 349.1 bar (for instance) are 
Δω_R/2π≃ 23.7 kHz, Δω_T/2π≃ 7.5 kHz, and Δω_J^km/2π≃ 51.8 kHz.
We also define the total relaxation from ^1H-^1H dipole-dipole interactions T_1RT as such:

    1/T_1RT = 1/T_1R + 1/T_1T

Given that the fast-motion regime applies, all of the above results for longitudinal relaxation T_1 apply equally to transverse relaxation T_2. As such, the subscript 2 has been dropped everywhere for clarity. 



  §.§.§ Distribution in correlation times
 

As shown in Fig. <ref>(a), G_J(t) deviates from single-exponential decay predicted by the Langevin model in Eq. <ref>. More specifically G_J(t) has a more “stretched" (i.e. multi-exponential) decay, which we quantify by inverting the following Laplace transform <cit.>: 

    G_J(t) = ∫ P_J(τ) exp(-t/τ) dτ

P_J(τ) (in units of s^-3) is the probability distribution function derived from the inversion. 

In the case of the Langevin sphere model, P_J(τ) is a delta-function at τ_J, i.e. P_J(τ) = G_J(0)  δ(τ - τ_J). However, as shown in Fig. <ref>(a), G_J(t) is always stretched (i.e. multi-exponential) to some degree, therefore P_J(τ) has a finite distribution. The decomposition of G_J(t) into a sum of exponential decays is common practice  <cit.>, where the more complex the molecule dynamics, the more terms are required. This justifies our general approach of decomposing G_J(t) into a “model free" sum of exponential decays in Eq. <ref>, for the purposes of quantifying the departure from the Langevin sphere model.

The resulting P_J(τ) distributions, shown in Fig. <ref>(b), were determined by using the discrete form of Eq. <ref>, using 100 logarithmically-spaced τ bins ranging from 0.03 ps ≤τ≤ 30 ps, and a fixed regularization parameter of 10^-2 <cit.>. Fig. <ref>(b) shows that the mean correlation-times τ_J get longer with increasing D_T, and that the width of the distributions get somewhat larger with increasing D_T. In all cases, the widths are much larger than the delta function prediction from the Langevin model in Eq. <ref>.

The widths of the P_J(τ) distributions were then quantified using the following:

    μ_J    = 1/G_J(0) ∫ P_J(τ) ln(τ)  dτ, 
    σ_J^2    = 1/G_J(0) ∫ P_J(τ) (ln(τ) - μ_J )^2 dτ.

σ_J is the standard deviation and μ_J is the mean of the variable ln(τ). A natural logarithm in τ is used as the variable since the underlying P_J(τ) distributions are discrete and evenly spaced in ln(τ). Also shown in Fig. <ref>(b) is the Langevin model which predicts σ_J = 0 (i.e. a delta function), or σ_J ≃ 0.038 due to regularization. Equivalent quantities were also derived for the intramolecular σ_R and intermolecular σ_T ^1H-^1H dipole-dipole interactions.



§ RESULTS
 



 §.§ Simulation versus measurement
 

The cross-plot of measured versus simulated translational-diffusion coefficient D_T in Fig. <ref>(a) indicates a strong correlation coefficient R^2  = 0.996, and an average absolute deviation of δ_abs = 9.2 %, where δ_abs is defined as:

    δ_abs = 1/N∑_i=1^N|Y_i - X_i/X_i|× 100.

Y_i are the simulated values, X_i are the measured quantities, and N is the number of points. The deviation is noticeable at the lowest temperature in the liquid phase (T<T_cr), which may be due to the proximity of the liquid-vapor phase transition, and/or temperatures are low enough that nuclear quantum effects may be important (further investigations are beyond the scope of this work).







The cross-plot of measured T_1 versus simulated T_1^α (defined in Eq. <ref>) total-relaxation time in Fig. <ref>(b) indicates strong correlation coefficients of R^2  = 0.982 for T_1^km, and R^2  = 0.954 for T_1^dm, while δ_abs varies as a function of measured T_1. At the higher D_T (i.e. lower ρ) end where the measured T_1 < 5 s, the deviation is lower for the kinetic model T_1^km (δ_abs = 11.0 %) than for the diffusion model T_1^dm (δ_abs = 15.8 %). This is expected given that the kinetic model is more appropriate for fluids at high D_T and low ρ. 

At the lower D_T (i.e. higher ρ) end where the measured T_1 > 5 s and T > T_cr (i.e. still in the supercritical phase), the deviation is lower for the diffusion model T_1^dm (δ_abs = 8.6 %) than for the kinetic model T_1^km (δ_abs = 19.7 %). This is expected given that the diffusion model is more appropriate for fluids at low D_T and high ρ. 

In the liquid phase where T < T_cr, ^1H-^1H dipole-dipole relaxation begins to dominate over the spin-rotation interaction, and the two interpretations T_1^km and T_1^dm become comparable. Nevertheless, the deviation is still lower for the diffusion model T_1^dm (δ_abs = 18.9 %) than for the kinetic model T_1^km (δ_abs = 29.1 %). The source of the deviation is from the dipole-dipole contribution T_1RT, more specifically from the dominating intermolecular contribution T_1T. This systematic deviation is potentially due to the proximity of the liquid-vapor phase transition, and/or temperatures are low enough that nuclear quantum effects may be important (further investigations are beyond the scope of this work).

Another potential systematic error are the experimental uncertainties in the coupling constants C_a and C_d <cit.> listed in Section <ref>. The maximum and minimum deviations reported in <cit.> result in a ±2.7 % uncertainty in T_1J^km, and a ±3.3 % uncertainty in T_1J^dm. These uncertainties are not insignificant in the above analysis, and should be taken into consideration when interpreting the simulations. There may also be uncertainty in the moment of inertia I for methane <cit.>.

Yet another potential systematic error in the measurements is the presence of dissolved oxygen <cit.>. O_2 is paramagnetic, which would tend to shorten the measured T_1 compared to simulations. Given the large values of measured T_1 ≤ 20 s, any trace amounts of oxygen could affect the results.





 §.§ Correlation times
 

A summary of the simulated correlation times is shown in Fig. <ref>(a). The spin-rotation correlation time τ_J shows a monotonic increase with increasing D_T/T (translational-diffusion divided by absolute temperature). At high D_T/T, i.e. in the supercritical/gas phase (T > T_cr), τ_J is consistent with the kinetic collision time τ_K defined as:

    D_T    = 1/3λv̅,   v̅ = √(8kT/π M),
    τ_K    = λ/v̅ = 3D_T/v̅^2 =  3π M/8k D_T/T,

where M = 2.664 × 10^-26 kg is the mass of the methane molecule, and v̅ is the mean thermal velocity. More specifically, the correlation between τ_J and τ_K is strong R^2 = 0.997 and the absolute deviation is δ_abs = 7.7 % in this region. Given the compelling observation τ_J ≃τ_K and the relation τ_K ∝ D_T/T in Eq. <ref> motivates using D_T/T for the x-axis in Fig. <ref>. The relation τ_J ≃τ_K is indicative of the “strong collision" regime <cit.>, and is used below to infer a new relation for T_1J^km. As expected, the relation τ_J ≃τ_K breaks down in the liquid phase (T < T_cr) where the diffusion model is more appropriate.



In the case of spherical molecules, the following simple theoretical predictions exist for relations between the intramolecular correlation time τ_R and spin-rotation correlation time τ_J:

    τ_R    = τ_J/2l+1 = τ_J/5    [ k.m.], 
    τ_R τ_J   = I/l(l+1)kT = I/6kT    [ d.m.].

In the kinetic model (T > T_cr) for the supercritical/gas phase, Eq. <ref> states that τ_R = τ_J/5 <cit.>, where l=2 is the rank of rotational diffusion tensor for hard-spheres. As shown in Fig. <ref>(a) for T > T_cr, the relation τ_R ∝τ_J is indeed found, however Eq. <ref> overestimates τ_R by a factor ≃ 2.5. In the diffusion model (T < T_cr) for liquids, Eq. <ref> <cit.> makes use of the Stokes-Einstein relation for hard-spheres, which predicts that rotational-diffusion D_R = 1/l(l+1)τ_R = 1/6τ_R and translational-diffusion D_T are related by D_T = 4a^2D_R/3 <cit.>, where a is the methane radius. At the lowest D_T value the relation τ_R ∝ 1/τ_J is indeed found, however Eq. <ref> underestimates τ_R by a factor ≃ 2.5. 

In the liquid phase (T < T_cr), the intermolecular correlation time τ_T decreases with increasing D_T, in general accordance with the Stokes-Einstein relation for hard-spheres τ_D = 52τ_T = 2a^2/D_T <cit.>. A transport radius of a ≃ 1.0  Å can be inferred at the lowest D_T value, which is consistent with the C-H internuclear distance of  a ≃ 1.093  Å deduced from the relation I = 83 m_H a^2 <cit.>, where m_H is the ^1H mass. At high D_T, i.e. in the supercritical/gas phase (T > T_cr), the Stokes-Einstein relation breaks down, and τ_T becomes independent of D_T.

A summary of the simulated standard deviations σ in correlation-times is shown in Fig. <ref>(b). For spin-rotation, σ_J increases by a factor of two from the lowest D_T to the highest D_T, indicating a larger distribution in correlation times for the gas phase. In the case of  intramolecular dipole-dipole, σ_R shows a low value in the liquid phase, and a possible increase with increasing D_T. In the supercritical/gas phase, G_R(t) shows signs of oscillations at early times t ≲ 0.3 ps, and then decays monotonically for t ≳ 0.3 ps. As such, σ_R is not computed for T>T_cr. Note that oscillations in G_R(t) for T>T_cr were previously predicted using the extended diffusion model <cit.>.

Meanwhile, σ_T decreases with increasing D_T, which is the opposite trend to σ_J and σ_R. At high D_T, σ_T plateaus to the σ_T ≃ 1.25 (note that Fig. <ref>(b) displays σ_T/4 for clarity), which is consistent with the inherent multi-exponential value of the Torrey hard-sphere model <cit.>. 

As stated in Eqs. <ref> and <ref>, the longer the spin-rotation correlation-time τ_J, the shorter the relaxation time T_1J, i.e. the more significant the relaxation mechanism. The same is true for the ^1H-^1H dipole-dipole mechanism <cit.>. It is therefore informative to compare the ratio of relaxation times between different mechanisms in order to determine which mechanism is dominant. Fig. <ref>(a) shows the ratio T_1RT/T_1J^km, where T_1RT is the total dipole-dipole relaxation defined in Eq. <ref>. A large value T_1RT/T_1J^km≫ 1 indicates that spin-rotation dominates over dipole-dipole, as found at high D_T in the supercritical/gas phase (T > T_cr). Meanwhile a small value T_1RT/T_1J^km≪ 1 indicates that dipole-dipole dominates over spin-rotation, as found at low D_T in the liquid phase (T < T_cr). 

A similar analysis can be made between the intramolecular and the intermolecular dipole-dipole interactions. A large value T_1T/T_1R≫ 1 indicates that intramolecular dominates over intermolecular, as found at high D_T in the supercritical/gas phase (T > T_cr). Meanwhile a small value T_1T/T_1R≪ 1 indicates that intermolecular dominates over intramolecular, as found at low D_T in the liquid phase (T < T_cr). The most likely reason for this is that τ_T ≫τ_R in the liquid phase (T < T_cr). Meanwhile in the supercritical/gas phase (T > T_cr), even though τ_T ≃τ_R, the intermolecular second-moment Δω_T^2 ∝ρ <cit.>, and therefore Δω_T^2 decreases with increasing D_T (i.e. decreasing ρ).







 §.§ New Kinetic model
 

Building on the observation that τ_J ≃τ_K in the supercritical/gas phase (T > T_cr) (see Fig. <ref>(a)), and the observation that spin-rotation dominates over dipole-dipole T_1RT/T_1J^km≫ 1 in this region (Fig. <ref>(b)), it is informative to infer a new relation for T_1J^km≃ T_1^km in the supercritical/gas phase. Assuming  τ_J = τ_K in Eq. <ref>, and using the definition for τ_K in Eq. <ref> results in the following prediction:

    1/T_1J^K = 3π I M/4 ħ^2(C_a^2 +  445C_d^2) D_T     [ k.m.].

The new prediction in the kinetic regime T_1J^K is plotted in Fig. <ref>(b) against the simulated T_1J^km and the measured D_T. In the supercritical/gas phase (T > T_cr), the correlation coefficient between T_1J^km and T_1J^K is found to be R^2 = 0.998, and the absolute deviation is found to be δ_abs = 7.7 %, without any adjustable parameters in the derivation of T_1J^K in Eq. <ref>. Moreover, there are no assumptions about an effective molecular radius. In the liquid phase (T < T_cr), T_1J^km starts to deviate from the kinetic model T_1J^K, as expected. Likewise, the total relaxation T_1^km deviates from T_1J^K due to the T_1RT contribution, as expected. 

Also shown in Fig. <ref>(b) is the empirical relation T_1J^A from Eq. <ref>, where A = 6.37 × 10^-6 <cit.> in units of temperature T (K), density ρ (g/cm^3), and T_1J^A (s). In the supercritical/gas phase (T > T_cr), the correlation between T_1J^A and T_1J^K is strong R^2 = 0.957, while the absolute deviation is δ_abs≃ 21.2 %, for possibly the same reasons as discussed in Section <ref>. Eq. <ref> and Eq. <ref> imply that the T and ρ dependence of T_1J comes entirely from D_T ∝ T^3/2/ρ. This is confirmed in the supercritical/gas phase (T > T_cr), where the correlation coefficient between measured D_T and measured T^3/2/ρ in Table <ref> is strong R^2 = 0.991.



§ CONCLUSIONS


We develop a minimization technique to compute the angular-velocity for non-rigid spherical molecules, which is used to simulate the autocorrelation function G_J(t) for the spin-rotation interaction of methane, over a wide range of densities ρ and temperatures T, spanning the liquid (T<T_cr), and supercritical/gas (T>T_cr) phases. The Langevin model predicts that G_J(t) should decay with a single-exponential function with correlation time τ_J. However, inverse Laplace transforms of G_J(t) indicate a distribution in correlation times τ, with a standard deviation σ_J (i.e. width) which increases with increasing D_T (i.e. decreasing ρ).

T_1J^α is derived from G_J(t) using the kinetic model “km" for gases (T_1J^km) <cit.>, and the diffusion model “dm" for liquids (T_1J^dm)  <cit.>. The total relaxation time T_1^km shows better agreement with measurements for T_1 < 5 s, with an absolute deviation of δ_abs = 11.0 % in this region. T_1^dm shows better agreement with measurements for T_1 > 5 s (and T>T_cr), with an absolute deviation of δ_abs = 8.6 % in this region. Uncertainties in the measured spin-rotation coupling-constants may contribute to these deviations. Meanwhile the simulated D_T agree well with measurements, without any adjustable parameters in the interpretation of the simulations.

MD simulations of the intramolecular and intermolecular ^1H-^1H dipole-dipole relaxation are computed at the same state-points. T_1J^km is shown to dominate over the total dipole-dipole relaxation T_1RT at high D_T, while the opposite is found at low D_T. 

The predicted relations between the intramolecular correlation-time τ_R and the spin-rotation correlation-time τ_J is tested, both in the liquid (T<T_cr) and the supercritical/gas phase (T<T_cr). At the highest D_T in the supercritical/gas phase, the relation τ_R = τ_J/5 <cit.> is found to hold within a factor ≃ 2.5. At the lowest D_T in the liquid phase, the relation τ_R τ_J = I/6kT <cit.> is also found to hold within a factor ≃ 2.5.

In the supercritical/gas phase (T>T_cr), τ_J is found to agree with the kinetic collision time τ_K, with an absolute deviation of δ_abs = 7.7 % in this region. Given this compelling finding, a new expression for the spin-rotation relaxation  1/T_1J^K∝ D_T is inferred without any adjustable parameters, and shows a strong correlation R^2 = 0.957 with the previously reported empirical finding 1/T_1J^A∝ T^3/2/ρ.



§ ACKNOWLEDGMENTS
 

This work was funded by the Rice University Consortium on Processes in Porous Media, and the American Chemical Society Petroleum Research Fund [ACS-PRF-58859-ND6]. We gratefully acknowledge the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy [DE-AC02-05CH11231], for HPC time and support. We also gratefully acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin (URL: http://www.tacc.utexas.edu) for providing HPC resources, and Zeliang Chen for his assistance.






10
	url<#>1urlprefixURL href#1#2#2 #1#1gutowsky:prl1961
	H. S. Gutowsky, I. J. Lawrenson, K. Shimomura, Nuclear magnetic spin-lattice
	relaxation by spin-rotational interactions, Physical Review Letters 6 (7)
	(1961) 349–351.
	
	hubbard:pr1963
	P. S. Hubbard, Theory of nuclear magnetic relaxation by spin-rotational
	interactions in liquids, Phys. Rev. 131 (3) (1963) 1155–1165.
	
	bloom:cjp1967a
	M. Bloom, R. Dorothy, Determination of the spin-rotation interaction constants
	in CH_4 by means of nuclear spin relaxation measurements, Can. J. Phys.
	45 (1967) 3411–3413.
	
	bloom:cjp1967b
	M. Bloom, F. Bridges, W. N. Hardy, Nuclear spin relaxation in gaseous methane
	and its deuterated modifications, Can. J. Phys. 45 (1967) 3533–3554.
	
	mcclung:JCP1969
	R. E. D. McClung, Rotational diffusion of spherical-top molecules in liquids,
	J. Chem. Phys. 51 (1969) 3842–3852.
	
	dong:cjp1970
	R. Y. Dong, M. Bloom, Determination of spin-rotation interaction constants in
	fluorinated methane molecules by means of nuclear spin relaxation
	measurements, Can. J. Phys. 48 (1970) 793–804.
	
	dawson:aiche1970
	R. P. Dawson, F. Khoury, R. Kobayashi, Self-diffusion measured in methane by
	pulsed nuclear magnetic resonance, AIChE J. 16 (5) (1970) 725–729.
	
	mcclung:JCP1971
	R. E. D. McClung, Rotational diffusion of spherical-top molecules in liquids.
	II. reorientation in liquid methane and solutions of methane in liquid
	noble gases, J. Chem. Phys. 55 (7) (1971) 3459–3467.
	
	gerritsma:ph1971a
	C. J. Gerritsma, N. J. Trappeniers, Proton-spin-lattice relaxation and self
	diffusion in methanes, I. spin-echo spectrometer and preparation of the
	methane samples, Physica 51 (1971) 365–380.
	
	gerritsma:ph1971b
	C. J. Gerritsma, P. H. Oosting, N. J. Trappeniers, Proton-spin-lattice
	relaxation and self diffusion in methanes, II. experimental results for
	proton-spin-lattice relaxation times, Physica 51 (1971) 381–394.
	
	oosting:ph1971a
	P. H. Oosting, N. J. Trappeniers, Proton-spin-lattice relaxation and self
	diffusion in methanes, III. interpretation of proton-spin-lattice
	experiments, Physica 51 (1971) 395–417.
	
	oosting:ph1971b
	P. H. Oosting, N. J. Trappeniers, Proton-spin-lattice relaxation and self
	diffusion in methanes, IV. self-diffusion in methane, Physica 51 (1971)
	418–431.
	
	rajan:jmr1974
	S. Rajan, K. Lalita, S. V. Babu, Nuclear spin-lattice relaxation in
	CH_4-inert gas mixtures, J. Magn. Reson. 16 (1974) 115–129.
	
	mcclung:AMRIP1977
	R. E. D. McClung, On the extended rotational diffusion model for molecular
	reorientation in fluids, Adv. Mol. Relax. Interact. Processes 10 (1977)
	83–171.
	
	mcclung:emr2007
	R. E. D. McClung, Spin-rotation relaxation theory, eMagRes DOI:
	10.1002/9780470034590.emrstm0524.
	
	cederberg:pr1964
	J. W. Cederberg, C. H. Anderson, N. F. Ramsey, Rotational magnetic moments,
	Phys. Rev. 136 (4A) (1964) 960–961.
	
	anderson:pr1966
	C. H. Anderson, N. F. Ramsey, Magnetic resonance molecular-beam spectra of
	methane, Phys. Rev. 149 (1) (1966) 14–24.
	
	yi:pr1968
	P.-N. Yi, I. Ozier, C. H. Anderson, Theory of nuclear hyperfine interations in
	spherical-top molecules, Phys. Rev. 165 (1) (1968) 92–109.
	
	ozier:pr1968
	I. Ozier, L. M. Crapo, S. S. Lee, Nuclear radio-frequency spectra of a series
	of tetrahedral molecules, Phys. Rev. 172 (1968) 63–82.
	
	yi:jcp1971
	P.-N. Yi, I. Ozier, C. H. Anderson, Low-field hyperfine spectrum of CH_4,
	J. Chem. Phys. 55 (11) (1971) 5215–5227.
	
	zhang:spwla2002
	Y. Zhang, G. J. Hirasaki, W. V. House, R. Kobayashi, Oil and gas NMR
	properties: the light and heavy ends, Soc. Petrophys. Well Log Analysts
	(2002) SPWLA–2002–HHH.
	
	lo:SPE2002
	S.-W. Lo, G. J. Hirasaki, W. V. House, R. Kobayashi, Mixing rules and
	correlations of NMR relaxation time with viscosity, diffusivity, and
	gas/oil ratio of methane/hydrocarbon mixtures, Soc. Petrol. Eng. J. 7 (1)
	(2002) 24–34.
	
	hurlimann:petro2009
	M. D. Hürlimann, D. E. Freed, L. J. Zielinski, Y.-Q. Song, G. Leu,
	C. Straley, C. C. Minh, A. Boyd, Hydrocarbon composition from NMR diffusion
	and relaxation data, Petrophysics 50 (2) (2009) 116–129.
	
	yang:petro2012
	Z. Yang, G. J. Hirasaki, M. Appel, D. A. Reed, Viscosity evaluation for NMR
	well logging of live heavy oils, Petrophysics 53 (1) (2012) 22–37.
	
	straley:spwla1997
	C. Straley, An experimental investigation of methane in rock materials, Soc.
	Petrophys. Well Log Analysts (1997) SPWLA–1997–AA.
	
	sigal:petro2011
	R. F. Sigal, E. Odusina, Laboratory NMR measurements on methane saturated
	barnett shale samples, Petrophysics 52 (1) (2011) 32–49.
	
	kausik:spe2011
	R. Kausik, C.-C. Minh, L. Zielinski, B. Vissapragada, R. Akkurt, Y.-Q. Song,
	C. Liu, S. Jones, E. Blair, Characterization of gas dynamics in kerogen
	nanopores by NMR, Soc. Petrol. Eng. SPE-147198-MS.
	
	tinni:spe2014
	A. Tinni, E. Odusina, I. Sulucarnain, C. Sondergeld, C. Rai, NMR response of
	brine, oil, and methane in organic rich shales, Soc. Petrol. Eng. (2014)
	SPE–168971–MS.
	
	wang:ef2014
	H.-J. Wang, A. Mutina, R. Kausik, High-field nuclear magnetic resonance
	observation of gas shale fracturing by methane gas, Energy Fuels 28 (2014)
	3638–3644.
	
	papaioannou:PRA2015
	A. Papaioannou, R. Kausik, Methane storage in nanoporous media as observed via
	high-field NMR relaxometry, ‎Phys. Rev. Applied 4 (2) (2015) 024018,
	1–11.
	
	sigal:spej2015
	R. F. Sigal, Pore-size distribution for organic-shale-reservoir rocks from
	nuclear-magnetic-resonance spectra combined with adsorption measurements,
	Soc. Petrol. Eng. J. 20 (4) (2015) 824–830.
	
	kausik:petro2016
	R. Kausik, K. Fellah, E. Rylander, P. M. Singer, R. E. Lewis, S. M. Sinclair,
	NMR relaxometry in shale and implications on logging, Petrophysics 57 (4)
	(2016) 339–350.
	
	valori:ef2017
	A. Valori, S. V. den Berg, F. Ali, W. Abdallah, Permeability estimation from
	NMR time dependent methane saturation monitoring in shales, Energy Fuels 31
	(2017) 5913–5925.
	
	tinni:petro2018
	A. Tinni, C. Sondergeld, C. Rai, New perspectives on the effects of gas
	adsorption on storage and production of natural gas from shale formations,
	Petrophysics 59 (1) (2018) 99–104.
	
	singer:jmr2017
	P. M. Singer, D. Asthagiri, W. G. Chapman, G. J. Hirasaki, Molecular dynamics
	simulations of NMR relaxation and diffusion of bulk hydrocarbons and water,
	J. Magn. Reson. 277 (2017) 15–24.
	
	bloembergen:pr1948
	N. Bloembergen, E. M. Purcell, R. V. Pound, Relaxation effects in nuclear
	magnetic resonance absorption, Phys. Rev. 73 (7) (1948) 679–712.
	
	torrey:pr1953
	H. C. Torrey, Nuclear spin relaxation by translational diffusion, Phys. Rev. 92
	(4) (1953) 962–969.
	
	namd
	J. C. Phillips, R. Braun, W. Wang, E. Tajkhorshid, E. Villa, C. Chipot,
	R. Skeel, L. Kale, K. Schulten, Scalable molecular dynamics with NAMD, J.
	Comput. Chem. 26 (2005) 1781–1802.
	
	cgenff
	K. Vanommeslaeghe, E. Hatcher, C. Acharya, S. Kundu, S. Zhong, J. Shim,
	E. Darian, O. Guvench, P. Lopes, I. Vorobyov, A. D. M. Jr., CHARMM general
	force field: A force field for drug-like molecules compatible with the
	CHARMM all-atom additive biological force field, J. Comput. Chem. 31 (2010)
	671–690.
	
	packmol
	L. Martinez, R. Andrade, E. G. Birgin, J. M. Martinez, Packmol: A package for
	building initial configuration for molecular dynamics simulations, J. Comput.
	Chem. 30 (2009) 2157–2164.
	
	yeh:jpcb2004
	I.-C. Yeh, G. Hummer, System-size dependence of diffusion coefficients and
	viscosities from molecular dynamics simulations with periodic boundary
	conditions, J. Phys. Chem. B 108 (2004) 15873–15879.
	
	kremer:jcp93
	B. Dünweg, K. Kremer, Molecular dynamics simulation of a polymer chain in
	solution, J. Chem. Phys. 99 (1993) 6983–6997.
	
	herzberg:book
	G. Herzberg, Infrared and Raman Spectra of Polyatomic Molecules, D. Van
	Nostrand Company, Inc. , New York, 1945.
	
	mcconnell:book
	J. McConnell, The Theory of Nuclear Magnetic Relaxation in Liquids, Cambridge
	University Press, 1987.
	
	cowan:book
	B. Cowan, Nuclear Magnetic Resonance and Relaxation, Cambridge University
	Press, 1997.
	
	venkataramanan:ieee2002
	L. Venkataramanan, Y.-Q. Song, M. D. Hürlimann, Solving fredholm integrals
	of the first kind with tensor product structure in 2 and 2.5 dimensions, IEEE
	Transactions on Signal Processing 50 (5) (2002) 1017–1026.
	
	song:jmr2002
	Y.-Q. Song, L. Venkataramanan, M. D. Hürlimann, M. Flaum, P. Frulla,
	C. Straley, T_1-T_2 correlation spectra obtained using fast
	two-dimensional laplace inversion, J. Magn. Reson. 154 (2002) 261–268.
	
	beckmann:prep1988
	P. A. Beckmann, Spectral densities and nuclear spin relaxation in solids, Phys.
	Rep. 171 (3) (1988) 85–128.
	
	shikhov:amr2016
	I. Shikhov, C. Arns, Temperature-dependent oxygen effect on NMR D-T_2
	relaxation-diffusion correlation of n-alkanes, Appl. Magn. Reson. 47 (12)
	(2016) 1391–1408.
	



