Nuclear spin diffusion under fast magic-angle spinning in solid-state NMR

Solid-state nuclear spin diffusion is the coherent and reversible process through which spin order is transferred via dipolar couplings. With the recent increases in magic-angle spinning (MAS) frequencies and magnetic fields becoming routinely applied in solid-state nuclear magnetic resonance, understanding how the increased 1 H resolution obtained affects spin diffusion is necessary for interpretation of several common experiments. To investigate the coherent contributions to spin diffusion with fast MAS, we have developed a low-order correlation in Liouville space model based on the work of Dumez et al. (J. Chem. Phys. 33 , 224501, 2010). Specifically, we introduce a new method for basis set selection, which accounts for the resonance-offset dependence at fast MAS. Furthermore, we consider the necessity of including chemical shift, both isotropic and anisotropic, in the modeling of spin diffusion. Using this model, we explore how different experimental factors change the nature of spin diffusion. Then, we show case studies to exemplify the issues that arise in using spin diffusion techniques at fast spinning. We show that the efficiency of polarization transfer via spin diffusion occurring within a deuterated and 100% back-exchanged protein sample at 60 kHz MAS is almost entirely dependent on resonance offset. We additionally identify temperature-dependent magnetization transfer in beta-aspartyl L-alanine, which could be explained by the influence of an incoherent


I. INTRODUCTION
Spin diffusion is a reversible and coherent process through which spin order may be transferred via dipolar couplings in the solid state. In 1 H solid-state nuclear magnetic resonance (NMR) at slow magic-angle spinning (MAS) frequencies, spin diffusion occurs in a manner analogous to macroscopic diffusion owing to the nucleus's low chemical shift dispersion and strong dipolar couplings. The spatial diffusional nature of this transfer has led to it being applied to the study of systems from materials 1 and biomaterials, 2,3 to small molecules 4-6 and proteins. 7,8 Additionally, spin diffusion plays an important role in dynamic nuclear polarization NMR (DNP-NMR). [9][10][11][12][13][14][15] Recently, experimental methods relying on selective pulses to exploit the increased resolution at faster MAS and higher magnetic fields have been introduced. For example, selective pulses have found use in reducing t 1 noise, 16 in increasing the rate of experimental acquisition, 17 and in selectively investigating pharmaceuticals in the presence of excipients. 18 In addition, low power pulses with narrow bandwidth were used for implementing chemical exchange saturation transfer (CEST) in the solid state, 19,20 where spin diffusion may be an alternative mechanism to chemical exchange that needs to be considered. 21 Modeling spin diffusion transfer is of key importance to understanding the results of such experiments.
The spin dynamics at slow spinning frequencies have been shown to be adequately reconstructed using diffusion-based perturbation theory simulation approaches. 4,6,22,23 In these approaches, perturbation theory is used to derive rate expressions that are then used to model the system as a diffusive process. 24,25 It has even been shown that such models are able to solve crystal structures from known unit cell parameters to excellent precision. 22 However, with the increase in resolution obtained using higher MAS frequencies and higher magnetic fields, the assumption that spin diffusion may be treated in an entirely spatial manner begins to break down. [24][25][26] As an energy conserving process, it follows that spin diffusion between spins with dissimilar energy level separations (i.e., different chemical shifts) is only possible if interacting with a spin energy bath, such as that provided by a dense dipolar coupled proton spin network. The decrease in spectral overlap with higher spinning frequencies arises because these dipolar coupling networks are more effectively averaged, and this combined with the larger energy level separations at higher magnetic fields means that spin diffusion becomes strongly dependent on the resonance offset between two spins. The importance of this resonance offset dependence was recently demonstrated in the work of Agarwal, 27 where it was shown that, in proton spin diffusion spectra of L-Histidine ⋅ HCl ⋅ H 2 O, negative cross peaks may be observed arising due to the interaction of four spins simultaneously, where the difference in differences between pairs of spins lead to a n = 0 rotational resonance transfer with inverse sign.
This four spin interaction effect would not arise in perturbation theory-/diffusion-based approaches, with the exception of qualitative models intended explicitly to study the effect. 28 Indeed, the majority of such models published to date either include the resonance offset through an exponential or Gaussian approximation of a zero-quantum line shape 23 or exclude it entirely. 4,22 Computational calculations in which the spin evolution of the density matrix is simulated under the spin Hamiltonian would, in theory, accurately reconstruct the coherent spin dynamics. Unfortunately, owing to their exponential scaling (∝2 n , where n is the number of spins), such simulations are typically restricted to systems with fewer than 12 spins. 23,29 As a result, they are unable to accurately model spin diffusion for which interactions with many more spins must be considered.
One approach that has been used to remedy this scaling problem is the use of restricted basis sets. [30][31][32][33] In such approaches, the number of basis states for which the evolution must be considered are drastically reduced by omitting those that can be assumed to contribute negligibly to the evolution of the spin system. Restricted basis set methods have been shown to enable accurate simulation of spin systems containing thousands of interacting spins. 11 In the work of Dumez et al., 31 the low-order correlation in Liouville space (LCL) method was introduced, where only zero-quantum operators are considered, and product states are limited to those containing at most q interacting spins. Such an algorithm scales polynomially as n q and allows for the number of spins in simulations to be drastically increased. The LCL method was further developed by Perras and Pruski, 32 who introduced local restriction (LR-LCL), where only the N closest spins to each spin were considered to be interacting, resulting in a linear scaling algorithm (∝n × N q−1 ). Such an approach has been applied to modeling DNP in systems containing thousands of atoms. 11,34 LR-LCL simulations are, however, considered to be accurate up to only ∼40 kHz MAS frequencies. 11 This limitation arises due to the aforementioned increasing dependence on resonance offset and the chemical anisotropy. Though resonance offset was included in the LR-LCL model introduced in the work of Perras and Pruski, 32 to our understanding, it was included solely as an isotropic chemical shift.
We introduce a new method for basis set selection and further develop the LCL method to include both isotropic and anisotropic chemical shift (chemical shift anisotropy, CSA). We then consider the effects of the full chemical shift, MAS frequency, magnetic field, and dynamics on the evolution of the spin system. We show experimental results for the dipeptide β-aspartyl L-alanine (β-AspAla), for which we find agreement with the simulated trends, but additionally observe temperature-dependent behavior, which may be indicative of an incoherent 1 H-1 H homonuclear nuclear Overhauser effect. Finally, we show that the efficiency of polarization transfer via 1 H spin diffusion in a deuterated fully back-exchanged protein at 60 kHz MAS is dominated by resonance offset.

DFT
The crystal structure of β-AspAla (CCDC: FUMTEM) 35 was geometry optimized by density functional theory (DFT) using CASTEP 16.1. 36-38 CASTEP implements DFT using a planewave basis set. The default CASTEP 16.1 ultrasoft pseudopotentials were used. The Perdew-Burke-Ernzerhof (PBE) implementation of the generalized gradient approximation was used as the exchange-correlation functional. 39 Plane waves up to 700 eV were used. The same cutoff energies were then used to determine magnetic resonance parameters using the Gauge-Including Projector Augmented Wave (GIPAW) method [40][41][42][43][44] under the same DFT conditions to determine the CSA tensors, which were then extracted using Magresview. 45

Spin diffusion simulations
We develop upon the excellent Tourbillon model introduced in the work of Dumez et al. 31 Specifically, we have further extended this model to use locally restricted basis sets and to implement the use of an unordered map for storing the density states, as introduced by Perras and Pruski. 32 Additionally, we implemented isotropic and anisotropic chemical shift evolution, the ability to output zero-quantum line shapes, and an implementation of our basis set selection method (see below). The code may be found at https://www.github.com/ThatPerson/Tourbillon_fastMAS. Spin diffusion simulations were run using a complete unit cell of β-AspAla (4 molecules in the unit cell and 12 1 H per molecule, i.e., n = 4 × 12 = 48 spins) using periodic boundary conditions. Unless otherwise stated, all simulations began with the inversion of all carboxylic acid protons (the site with the most separated 1 H NMR resonance), for which the interatomic distances and resonance offsets are shown in Table I. This system is illustrated in Fig. 1, where the spins of particular interest are color coded: Asp HA (orange), the spatially closest proton to COOH; Ala NH (green), the closest in chemical shift to COOH; and Ala CH 3 (lilac), for which there is particularly interesting spin evolution. In the case of the alanine CH 3 , we only considered one of the protons (labeled 12 in Table I) when plotting the trajectories as the evolution differs slightly between nonsymmetrically equivalent sites. The spins are numbered 1-48, where the protons are numbered sequentially for each individual β-AspAla molecule, i.e., molecule 1 is numbered 1-12, molecule 2 is numbered 13-24, etc. Experimental isotropic chemical shifts were used, 46 but the CSA tensors were calculated as described in a prior section. We note that these experimental isotropic chemical shifts differ from those given at lower spinning frequencies likely due to sample rotation heating. 47 The rotational motion of the NH 3 and CH 3 groups was considered by assuming averaging of the chemical shift tensors in the molecular frame prior to conversion into the  interaction frame for both of these sites, though no explicit averaging of dipolar couplings was considered unless indicated explicitly; the effect of dynamics on spin diffusion will be considered in Sec. IV C. Simulations were performed for t mix = 100 ms [see Fig. 2(a)] using a REPULSION-48 set of crystallite orientations 48 and a time step of 0.2 μs. Simulations were run using the University of Warwick Scientific Computing Research Technology Platform (SCRTP) High Performance Computing clusters, on nodes consisting of two Intel Xeon 24 core processors giving 48 cores per node. Parallelism was implemented with each crystallite running in an individual thread using OpenMP. 192 GB of RAM was present per node; however in the case of some larger models, high memory nodes were used with up to 1.5 TB of RAM. We additionally ran simulations using a REPULSION-128 set of crystallite orientations using the HPC Midlands Tier 2 High Performance Computing cluster Sulis, on nodes containing two AMD EPYC 7742 (Rome) 2.25 GHz 64 core processors, giving 128 cores and 512 GB of RAM per node. There was no additional benefit to using more crystallites (see Fig. S1).
B. Experimental solid-state NMR 1. Samples β-AspAla was purchased from Bachem (Switzerland) and packed as received into both a 1.3 mm zirconia rotor and a 1.6 mm zirconia rotor. The 1.6 mm rotor then had a plug of KBr inserted prior to cap insertion. Perdeuterated ( 2 H, 13 C, 15 N) GB1 with 100% back-exchange was prepared as described in Ref. 49 and packed into a 1.3 mm zirconia rotor, sealed with silicon glue.

Spectrometer and MAS probe
Experiments were performed at three fields: a 1 H Larmor frequency of 599.5 MHz with a Bruker Avance Neo console using a Bruker 1.3 mm HXY probe; a 1 H Larmor frequency of 850 MHz with a Bruker Avance Neo console using a Phoenix 1.6 mm HXY probe; and a 1 H Larmor frequency of 1 GHz with a Bruker Avance Neo console using a Bruker 1.3 mm HCN probe.

1D 1 H selective spin diffusion experiments on β-AspAla
Spin diffusion experiments were performed using the 1.3 mm HXY probe operating in double resonance mode spinning at 55 kHz at a 1 H Larmor frequency of 599.5 MHz. 1D proton spin diffusion (PSD) experiments 50 were performed in a manner analogous to saturation transfer difference methods in the solution state, 51,52 using variable length "trains" of Gaussian inversion (180 ○ ) pulses for the saturation of the highest ppm resonance, the carboxylic acid resonance, at 13.1 ppm. Gaussian inversion pulses were optimized for a pulse length of 1.1 ms. Hard pulses were applied to 1 H with a nutation frequency of 100 kHz, corresponding to a 90 ○ pulse length of 2.5 μs. The pulse sequence for this is shown in Fig. 2(b), where it is compared with the effective pulse sequence for the simulations [ Fig. 2(a)]. Note that simulations were performed using a single inversion [ Fig. 2(a)] as the Gaussian pulse trains used experimentally have the potential to generate nonzero-quantum coherence, which cannot be represented using our model. On the other hand, experimentally saturation was used as this ensures a greater degree of magnetization transfer. The number of inversion pulses applied was varied linearly from 0 to 100, with eight coadded transients acquired per increment. Except where stated, a recycle delay of 2 s was used. The resulting data were then analyzed by taking the peak intensities for each peak. Spectra were referenced according to the resonance offset for the methyl group used for simulations at 1.4 ppm. 46 Temperatures were calibrated using the chemical shift dependence on temperature of the 79 Br resonance in KBr. 53 Further experiments were performed: 1D PSD experiments were performed using a 1.6 mm HXY probe at 40 kHz at a 1 H Larmor frequency of 850 MHz over a wider range of temperatures; these results are shown in the supplementary material (Figs. S7 and S8). A 2D PSD experiment was performed using a 1.3 mm HCN probe at 60 kHz at a 1 H Larmor frequency of 1 GHz.

Spin diffusion in GB1 probed via a 2D 1 H-15 N Experiment
Experiments were performed using the 1.3 mm HXY probe operating in triple-resonance HCN mode spinning at 60 kHz. Cooling was applied to achieve a sample temperature of ∼300 K as calibrated using the difference in shift between 2,2-dimethyl-2-silapentane-5-sulfonate sodium salt (DSS) and H 2 O. 54 2D spin diffusion experiments were performed in a manner analogous to saturation transfer difference methods in the solution-state 51,52 or CEST experiments. 19 Low-powered (∼100 Hz) bandwidth saturation was applied for 500 ms, both on resonance with various N-H sites and off resonance. The pulse sequence for this is shown in Fig. 2(c). A total of 32 coadded transients were acquired for each of 128 t 1 free-induction decays (FIDs), using States-TPPI in F 1 for sign discrimination and a 2 s recycle delay. Except for the saturation pulse, the 1 H transmitter was placed at 2.46 ppm with a spectral width of 39.7 ppm. The 15 N transmitter was placed at 120 ppm, with a spectral width of 54.9 ppm. Hard pulses were applied to 1 H with a nutation frequency of 100 kHz, corresponding to a 90 ○ pulse length of 2.5 μs. On 15 N, hard pulses had a nutation frequency of 50 kHz, corresponding to a 90 ○ pulse length of 5 μs. Cross-polarization (CP) was applied with a 70:100% linear ramp 55 on 1 H, meeting the Hartmann-Hahn condition with 45 kHz on 1 H and 15 kHz on 15 N. [56][57][58] This was applied for 800 μs for the 1 H → 15 N CP and for 600 μs for the 15 N → 1 H CP. 100 ms MISSISSIPPI water saturation using 10 kHz 1 H irradiation was used. 59 Spectra were referenced according to DSS at 0 ppm, with 15 N referenced indirectly to liquid NH 3 at 0 ppm using the IUPAC recommended frequency ratio. 54,60 Low-powered (ν 1 = 10 kHz) WALTZ-16 heteronuclear decoupling was applied to the 1 H channel during 15 N evolution and again to the 15 N channel during 1 H acquisition. 61 C. Experimental processing β-aspartyl L-alanine experiments: Recorded experimental spectra were integrated using a python script using both nmrglue 62 and lmfit 63 to fit a polynomial background and Lorentzian line shapes to each peak (code can be found at https://github.com/ThatPerson/PeakInt1D). Following optimization of the offset on the first slice, the offsets were fixed and only the peak height and width were allowed to vary. Fitting of the peak width was performed to account for noise, noting that we would not expect variation in peak width with the duration of the saturation pulse.

III. THEORY
A. Low-order correlations in Liouville space theory A brief introduction to the LCL approach is given here, based on that presented in the work of Brüschweiler and Ernst, 66 Dumez et al., 31 and Perras and Pruski. 32 For clarity, a full list of symbols and associated descriptions is given in Table II.
In Liouville space, the density matrix for a set of n spin-1/2 nuclei may be represented as a 4 n vector, consisting of elementŝ whereÎi,r represent single spin operatorsÎi,r ∈ {Ê i,Îiz ,Îi+/ √ 2,Îi−/ √ 2}, and q represents the order of the product operator, that is, the number of nonidentity operators in the product operator, corresponding to the number of interacting spins.
The density matrix,σ(tmix) after a mixing time t mix , is then a combination of these elements, which may be propagated under the action of the Liouvillian superoperator,L, This propagation may be approximated by considering the Liouvillian to be piecewise time independent, where Δt is a small step such that P × Δt = tmix. The low-order correlation in Liouville (LCL) method relies on basis set reduction by excluding certain terms. 30,31,67,68 For the simulation of spin diffusion, the most obvious terms to omit are those that are not zero-quantum since spin diffusion is fundamentally a zeroquantum process. 24,25 Further basis set reduction can be performed as considered below in Sec. III B below.
Owing to the very large size of the full Liouvillian,L, we make use of the first-order Suzuki-Trotter algorithm 69 as was done in the work of both Dumez et al. 31 and Perras and Pruski. 32 This enables Length of simulation, corresponding to PSD mixing time (see Fig. 2

Br
Arbitrary product operator, defined as in Eq. (1).
the density matrix evolution to be calculated piecewise, without requiring storage of the full Liouvillian. That is, where the inner product represents the product over all states considered to be "neighbors" of spin i, noting that, in the case of the basis set methods used here, these are not necessarily the most spatially proximate sites (see Sec. III B). Higher-order Suzuki-Trotter algorithms were not applied, since these double the computational time and have minimal benefit for such simulations. 32 The exp (iLi, j (pΔt)Δt) terms in Eq. (5) are evaluated sequentially using the cog-wheel approach used in the work of Brüschweiler and Ernst. 66 Under a dipolar coupling, ωD,i, j (t), subspaces of the form {Î izBr ,Î jzBr ,Îi+Î j−Br ,Îi−Î j+Br } (whereBr does not include spins i or j and must be such that the total product operator formed is zero-quantum in nature) evolve under the effect of the following rotation:R where c = cos (ωD,i, j (pΔt)Δt) and s = sin (ωD,i, j (pΔt)Δt). An example subspace would be {Î 1zÎ3z ,Î 2zÎ3z ,Î 1+Î2−Î3z ,Î 1−Î2+Î3z }, whereBr =Î 3z , which would evolve under the rotationˆR 0,1,2 .
Under MAS, the dipolar coupling is calculated in a timedependent fashion as where Ω PC is the change in orientation going from the principal axis system of the dipolar interaction to the crystalline frame and Ω CR relate this crystalline frame to the rotor frame (that is the powder crystallites). Finally, Ω RL (t) is the time-dependent (under rotation) The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp conversion from the rotor frame to the laboratory frame. The dipolar coupling constant, R DD,i,j , is given as in units of rad s −1 . The evolution under chemical shift is treated in an analogous manner. The time-dependent chemical shift is calculated as with δ (l) where δ iso is the isotropic chemical shift (in ppm), δ aniso the anisotropic chemical shift (in ppm), and η the chemical shift anisotropy. 70,71 These are obtained from the diagonalized chemical shift tensor δ components, δxx, δyy, and δzz as 71 The evolution under chemical shift is then considered for any product operator of the formÎi±Br, as where B Hz/ppm is the magnetic field strength in units of Hz ppm −1 , e.g., 599.5 Hz ppm −1 for a 1 H Larmor frequency of 599.5 MHz, and b i±,r (t) is the contribution of the product operatorÎi±Br to the density matrix as defined in Eq. (2). Note that product operator stateŝ IizBr would not evolve under this term. Computationally, evolution under the dipolar and chemical shift interactions was implemented sequentially such that, for each time step, the dipolar evolution under Eqs. (6) and (7) was applied first, followed by the chemical shift evolution under Eq. (15).

B. Basis set selection
Ideally, we would like to consider a basis set in which we consider all spin interactions. Such a simulation, akin to a full density matrix treatment as implemented by SIMPSON 70,72,73 and SpinEvolution, 74 would scale exponentially (∝2 3n assuming Hilbert space simulations). 73 To reduce the complexity of the problem, the simulation can be approximated using a reduced basis set. In the work of Edwards et al., 75 three basis sets were introduced for the solution state with the aim to reduce the number of terms that must be included: (1) IK-0, in which only product operator states consisting of a given number of spins are considered; (2) IK-2, in which only product operator states containing spins that are close in proximity are included; and (2) IK-1, where both criteria are applied. In the solid state, the analogous restricted basis sets used in the work of Dumez et al. 31 correspond to an IK-0 basis set. The local restriction introduced by Perras and Pruski 32 further develops this into an IK-1 basis set.
The computational complexity of the locally restricted loworder correlations in Liouville space model is approximately given by where P represents the number of time steps that must be calculated, n is the number of spins in total, and Nmax is the number of spin neighbors considered per spin. In the case of the original LCL algorithm (i.e., Nmax → n − 1), the complexity goes as n q , that is, it will scale polynomially. On the other hand, introducing local restriction ensures that Nmax is independent of n. The complexity then scales linearly in n. The choice of Nmax is important for such linear scaling models because including too few spin pairs will lead to the simulation not adequately recreating the dynamics of a larger model, while including too many will not be an efficient use of computational time. In Fig. 3(a), it can be seen that at least 16 spin pairs per spin are required to recreate the spin dynamics, while in Fig. 3(b), increasing the average number of spin neighbors (note Navg; the difference between this and Nmax will be discussed below in the description of our computational implementation) leads to a polynomially greater computational time. We have therefore chosen to include on average 16 spin pairs per spin, in good agreement to the 15 found to be optimal by Perras and Pruski, 32 who found that ∼15 spin pairs per spin were able to recreate spin diffusion well in their system, a linear alkane.
In our modeling, we limited ourselves to product operator states containing at most q = 4 spins, an upper limit that has been shown to allow for accurate simulation when compared to exact simulation at slower spinning. 68 While a larger value of qmax may model the spin evolution better, 30,31 such a larger basis set is infeasible. The faster spinning frequencies modeled here mean that the rotor cycle time is shorter and, therefore, more steps are required to sample the rotation adequately for any given time, and additionally the slower rate of spin diffusion means that longer times must be simulated. A significantly larger absolute number of steps, P, therefore needs to be calculated. As such, given the exponential scaling nature of q, 30,31 [see Eq (16)] it is not feasible to simulate with larger product operators included.
A method for selecting which spin pairs contribute significantly is necessary to perform a basis set reduction in spin space. It has been shown that applying spatial restrictions (i.e., selecting the Nmax closest spins) to the choice of spin states to include is valid [30][31][32]68,75 when the evolution of the spin dynamics is determined predominantly by through-bond J-couplings (as in solution-state NMR) or dipolar couplings in the presence of a strongly coupled spin bath (as at slow spinning frequencies). However, this assumption begins to fail for solid-state NMR in the fast spinning regime. 26,27,46 At of Chemical Physics The simulated z-magnetization evolution rapidly converges, with a basis set consisting of on average 16 spin neighbors per spin approximating well a system containing an average of 24 spin neighbors. Only up to a spin system containing 24 spin neighbors per spin was considered as above this was not computationally feasible.
attainable fast spinning frequencies, the dipolar couplings are not averaged to zero as in the solution state but are averaged sufficiently to make the transfer of magnetization strongly truncated by resonance offset.
To consider this transition from slow to fast MAS, we introduce a method in which a short simulation is performed including all spin pairs in which we restrict the basis sets to product operators containing q ≤ 4 spins. We then calculate a score for each spin pair, based on the population of all spin operators that include both spins, where Wn is a weighting to ensure that all spin orders contribute equally and bp(t) refers to the population of spin states made up of the spins within the set p. Weighting is necessary, as there are n times as many four spin states as there are three spin states, and n times as many three spin states as two spin states. Note that this simple weighting does not a priori assume any dependence on the spinning frequency, which may emerge from the simulation. Figure 4 shows that the contributions of these states no longer significantly vary such that they would lead to different spin interactions giving a greater contribution after 0.1 ms; in our modeling, we ran these simulations for 0.2 ms for this selection procedure. Our selection method corresponds to ordering the spin pairs from the highest score to lowest score and selecting the n × Navg highest scoring spin pairs to include in the basis set, thereby including only the states that contribute the most to the spin evolution. The spin pairs that contribute the most are not necessarily those closest in proximity; Fig. 5 shows the magnitude of the score in Eq. (17) as a function of distance and resonance offset for a methyl group proton (spin-10 in Table I), where there is no clear dependence on either parameter. Figure S2 depicts the spatial arrangement of the spins included as interacting for this methyl group proton at νr = 20 kHz and ν 0 ( 1 H) = 600 MHz.
Note here that we use the average number of spin neighbors per spin (Navg) as opposed to a fixed number of spin neighbors per spin (Nmax). This is because when ordering the spin pairs by score, we  Fig. 2(a) including all spin pairs at νr = 60 kHz, ν 0 ( 1 H) = 600 MHz, beginning with inversion of the carboxylic acid proton. LCL simulations without restriction are run for 1000 steps (0.2 ms), with the weighted population of two, three, and four spin states plotted above. Spin pairs involving the COOH and the Ala NH, Asp HA, and Ala CH 3 in β-AspAla are shown in green (solid), orange (dashed), and lilac (dotted), respectively (see Fig. 1 and Table I). Beyond ∼0.1 ms, the contribution of the 2, 3, or 4 spin states no longer significantly changes such that the ordering of which spin interactions to include is likely fixed at this point.
find that some spins are generally more important to the evolution of other spins; for example, the carboxylic acid protons are inverted at the start of the simulations and as a result have the greatest magnetization gradient with their neighbors and therefore tend to have a much larger contribution to many other spins. This is seen in Fig. 6, where generally the carboxylic acid protons (gray) are found to contribute more to the spin evolution of far more other spins than is the case for the NH 3 protons (lilac). This is also evident in Fig. 5 where the biggest circles (highest scores) are at the greatest resonance offset. At faster magic-angle spinning frequencies, the relative spin pair scores change as shown in Fig. 5, leading to a different basis set selection, thereby making different spin pairs contribute more or less as in Fig. 6. For example, two spins that are spatially proximate but of Chemical Physics  well separated in resonance frequency will contribute more at lower spinning frequencies, but, with the truncation of resonance offset at faster spinning frequencies, they will consequently contribute less and therefore be less likely to be included. This has important consequences for the evolution of the spin system as will be discussed in the following sections.

A. Effect of resonance offset on spin diffusion
The chemical shift of a spin relates to how the local environment changes the effective magnetic field experienced by that spin and therefore the change in the difference in energy between the spin-up and spin-down states. Spin diffusion, being an energy conserving diffusion of spin order between sites, is therefore truncated by the offset between two resonances corresponding to two chemically distinct sites in the solid-state structure. Dipolar couplings to the large proton bath provide an external reservoir of energy to allow for spin diffusion to occur as apparent by the broad overlapping peaks in 1 H solid-state NMR spectra. In the limit of low magnetic fields and low spinning frequencies, the contribution from this spin bath is sufficient to ensure that the rapid magnetization transfer between different spins occurs in an approximately spatial manner (i.e., that the efficiency of the polarization transfer is approximately related to the distance between spins). 22,23 With faster magic-angle spinning and higher magnetic fields, however, the spatial nature of this transfer begins to break down. Figure 7 shows the effect of including both isotropic and anisotropic chemical shift in our spin diffusion model. It is notable that, in the absence of isotropic or anisotropic chemical shift (blue dotted lines), the reduction in the rate of spin diffusion is negligible as higher spinning frequencies are attained, as seen by comparing with the red dashed lines. This suggests that it is the truncation of Chemical Physics  Table I), after inversion of the COOH proton resonance [see pulse sequence in Fig. 2(a)].
by resonance offset that dominates the reduction in spin diffusion rather than this being due to the direct averaging of dipolar couplings. Note that, in the long-time limit, the curves converge at 83.3%. This occurs because, in the absence of any relaxation, the diffusion will eventually cause the system to equilibrate at the average of the initial magnetization. Noting the inversion of some of the spins, namely, the carboxylic acid protons, this average is then {44 × (+1) + 4 × (−1)}/48 = 0.833. Figure 7 further exemplifies the importance of including not only the isotropic component of the chemical shift tensor but also the anisotropic component. In the absence of anisotropic chemical shift, it appears that the rate of spin diffusion is further reduced as the anisotropic component of the chemical shift tensor does not commute with the homonuclear dipolar interactions and so aids the spin diffusion.
The increasing dependence of spin diffusion on resonance offset between sites may also be seen in Fig. 8. Here, the transfer of magnetization for the Ala NH (4.6 ppm offset from COOH), Asp HA (8.6 ppm offset), and an Ala CH 3 proton (11.7 ppm offset) is shown as a function of spinning frequency (10-150 kHz) and magnetic field ( 1 H Larmor frequencies of 100 MHz, 600 MHz, and 1 GHz). At 100 MHz, it is observed that there is little truncation due to offset even for high spinning frequencies (≥60 kHz), indicating little isotropic chemical shift resolution even at these higher spinning frequencies. This is exemplified by comparison between 20 kHz MAS at 1 GHz and 60 kHz MAS at 100 MHz, where the rate of transfer appears qualitatively similar. At 1 GHz, spinning frequencies greater than 100 kHz would be expected to severely truncate coherent spin diffusion.
Coherent inverse sign spin diffusion, as observed in L-Histidine ⋅ HCl ⋅ H 2 O at 60 kHz at 500 and 700 MHz, 27 is also seen. In these simulations, such spin diffusion is manifested by the magnetization associated with a site increasing, e.g., Mz > 1, as is evident for the inset axes on Fig. 8 (l, n, o, q, and r). This is expected to occur when four interacting spins, I, S, R, and P, meet a resonance condition of (ωI − ωS) − (ωR − ωP) = 0. In this system, we envisage that it is the Ala COOH (13.1 ppm), Ala NH (8.5 ppm), Ala HA (5.4 ppm), and Ala CH 3 (1.4 ppm) protons that approximately meet this four spin resonance condition, as 8.5 − 1.4 = 7.1 ≈ 13.1 − 5.4 = 7.7. It is found to be particularly relevant in this system for the methyl protons (lilac) at spinning frequencies >60 kHz at 1 GHz, and >100 kHz at 600 MHz; a more thorough analysis will be provided in Sec. IV C.

B. Relation to perturbation approaches
While the restricted basis set low-order correlation in the Liouville space model used here is orders of magnitude faster than simulating the complete evolution of the full density matrix, it is still far slower than perturbation theory diffusion equation-based approaches. As such, relating the resonance-offset dependence here to these methods is of interest to the broader applicability of such models at faster spinning frequencies. Kubo and McDowell (1988) give an equation for the dependence of the rate of spin diffusion on the dipolar couplings and the resonance offset (via the "zero-quantum line shape") under MAS [their Eq. (25)]. 24 Their equation is reproduced here, where KR y (ω) represents the zero-quantum line shape function and RDD,i,j is the dipolar coupling as defined in Eq. (9). This equation, or one of its equivalent equations, forms the basis of most diffusionbased simulations of spin diffusion. The zero-quantum line shape function arises from the evolution of the zero-quantum coherence (e.g., the zero-quantum transition operator, ZQT = I+S−). 76 In such models, the dependence of spin diffusion on both spinning frequency and resonance offset arises via this line shape function. As such, understanding the nature of this function has been of interest. 77,78 Typically, however, this has been done using convolutions of single-quantum line shapes 78 or through master equation-based modeling. 77 By comparison, here, we present simulations in Fig. 9 at ν 0 = 600 MHz, where we consider the evolution of the I+S− state for different spinning frequencies during the normal evolution of the system after inversion of the carboxylic acid proton resonance. In Sec. S2 of the supplementary material (see Figs. S4-S6), we present the resulting ZQ line shapes arising from higher-order terms that are zero-quantum in these same spins (e.g., I+S−Rz). We find that these higher-order zero-quantum operators show the same trend upon increasing MAS frequency, however, with increased broadening.
Specifically, the zero-quantum line shape of the first spinning sideband between an Asp NH 3 + proton and Ala NH protons of molecule 1 extracted from our model are shown in Fig. 9 for a range of spinning frequencies. At spinning frequencies ≤15 kHz, the line shape is broadened such that it is not possible to identify a single peak, which suggests that spin diffusion in this regime will have little dependence on the exact form of this function. As the spinning frequency increases, however, the line shapes narrow significantly. It is also noted that the intensity of the spinning sidebands is significantly reduced with higher MAS frequencies, reflecting the reduction in spin diffusion rate with spinning frequency. A result of this zero-quantum line narrowing is that we would expect spin diffusion to occur only within a well-defined range of resonance offsets, with relayed transfer slightly broadening this by allowing magnetization to transfer between rather more separated spins. Modeling this line narrowing in the manner utilized here may enable the development of perturbation theory approaches at faster MAS frequencies by accounting more directly for the change in ZQT line shape, rather than approximating this effect with a simple decaying exponential as performed in earlier work. 23

FIG. 9.
Simulated zero-quantum line shapes at the first spinning sideband as obtained from evolution of the I+S− ZQT coherence after inversion of the carboxylic acid proton resonance. Simulations were run with ν 0 ( 1 H) = 600 MHz and variable νr, for a full unit cell of β-AspAla (n = 48) with periodic boundary conditions. The interactions between spin-2 and spin-5 are shown where spin-2 and spin-5 are the NH 3 + and the Ala NH proton of molecule 1, respectively (see Table I).

C. Effect of dynamics on spin diffusion
Many of the systems of interest for characterization with spin diffusion-based techniques exhibit dynamics. 2,[79][80][81] Dynamic motions affect the coherent evolution in spin space by leading to averaging of dipolar couplings and CSA tensors. 82,83 Even in crystalline solid samples, it is known that significant motions may be present. [83][84][85][86][87][88] Of particular importance are rotational motions of axially symmetric groups, such as methyl groups, primary amines, and phenyl groups, which are typically fast on the NMR timescale and with limited energetic barriers. 89,90 FIG. 10. The simulated variation in spin diffusion evolution is shown for the Ala NH, Asp HA, and Ala CH 3 sites after inversion of the carboxylic acid resonance (see pulse sequence in Fig. 2(a) under different models of dynamical averaging. All dipolar couplings within each system were scaled by the parameter S to reflect overall dynamical motion. In (b), (d), (f), all dipolar interactions within a CH 3 or NH 3 group were additionally scaled by 0.5 to represent their axial rotation. Simulations were run with ν 0 ( 1 H) = 600 MHz and νr = 60 kHz, for a full unit cell of β-AspAla (n = 48) with periodic boundary conditions.
Motions occurring at frequencies coincident with various spin interactions can give rise to relaxation. In the case of interacting proton spins, motion could lead to incoherent cross relaxation and the nuclear Overhauser effect. Such an effect would interfere with the coherent spin evolution. A relaxation superoperator-based treatment 91 of such effects is beyond the scope of this paper, where we are solely interested in the coherent evolution of the system.
In our model, we introduce dynamics through two parameters. Overall motion is reflected by an order parameter S, which applies a linear scaling to all dipolar couplings between spins. C 3 rotation of the NH 3 and CH 3 groups was modeled by averaging the dipolar couplings between protons within the same group by P 2 (cos(θ)) (= 1 2 (3 cos 2 (θ) − 1)), where θ is the angle between the rotational axis and the interaction of interest. The angle is usually taken to be 90 ○ so that the averaging factor is 0.5. As stated previously, we also averaged the CSA tensors in the molecular frame between all protons within a given CH 3 or NH 3 group prior to rotation into the interaction frame. 77 The spin evolution as shown in Fig. 10 (60 kHz MAS, 600 MHz) appears relatively insensitive to the overall motion order parameter. Under CH 3 or NH 3 rotation, however, the dynamics significantly change. Given that six of the 12 protons in each molecule of β-AspAla are within these axially symmetric groups, this is unsurprising.
Of important note is the behavior of the inverse sign spin diffusion (Mz > 1) apparent for the methyl (lilac) group. As was noted by Agarwal, 27 this inverse transfer is aided by scaling the couplings involved. For this effect to occur, the four spin third-order Hamiltonian contribution must dominate over the two and three spin Hamiltonians. While at such fast spinning frequencies, the two spin terms are effectively completely averaged, relayed transfer may still come to dominate over the four spin interactions at long simulation times, hence giving the short lifetime of the effect observed in Fig. 10(e) in the absence of methyl rotation. When these relay interactions are further averaged by the methyl rotation, however, the four spin term can dominate for far longer as seen in Fig. 10 Table I]. In these experiments, variable length trains of selective Gaussian inversion pulses were applied selectively to the carboxylic acid proton peak [see Fig. 2(b)] and the corresponding intensity of all other sites recorded in a manner analogous to saturation transfer difference experiments conducted in the solution state. 51,92 Under these conditions, it is expected that the transfer of spin order will occur via spin diffusion with the saturated site acting as a magnetization sink. We chose to apply a saturation type experiment because, experimentally, we would experience influence from incoherent relaxation processes and so would expect a fast recovery of the COOH spin making analysis of the behavior more complex. On the other hand, under saturation, the system would evolve toward a "steady state" and, thus, the rate at which this occurs may be more readily quantified. Additionally, mis-set inversion pulses may give rise to transverse components that cannot be simulated in our  Fig. 2(b)], at temperatures of 279-320 K (calibrated using KBr). There is a significant difference in the spin diffusion observed under each temperature condition. Errors are shown as twice the standard error from voigt line shape fitting with lmfit. Simulated spin diffusion curves under the same conditions are shown in red. (g) A 2D 1 H-1 H spin diffusion spectrum recorded at 60 kHz MAS and a 1 H Larmor frequency of 1 GHz with a mixing time of 50 ms. model, and may not completely invert the site, meaning that the initial state would not be completely along −z as simulated. On the other hand, saturation would cause the magnetization of this site to approach 0.
While we additionally present a comparison with simulation in Figs. 11(a)-11(f), we note that the presence of incoherent autoand cross relaxation means that full agreement cannot be expected: Experimentally, the spin diffusion is a combination of incoherent and coherent processes, while our modeling here is purely coherent.
We observe that, under saturation, there is a strong resonance offset dependence with regard to the magnetization change under saturation. Figures 11(a)-11(f) are ordered with the site highest in chemical shift (closest in resonance offset to the carboxylic acid) at the top, with the chemical shift progressively getting lower (with greater resonance offset) on moving down the figure. Interestingly, we note that the Ala CH 3 group experiences significant inverse sign spin diffusion similar to that which has been seen in other systems, 27,28,93 with an enhancement of ∼20%. This is significantly greater than any inverse sign spin diffusion we observed within our simulation [the maximum seen being 1.82%, in Fig. 10(f) with S = 0.94]. Moreover, we find that this effect exhibits a significant temperature dependence, with an increase in the amount of inverse sign spin diffusion at higher temperatures. Section S3 of the supplementary material investigates this temperature dependence in more detail. This observation is in contrast to the similar effect observed for L-Histidine by Agarwal 27 and suggests that this effect may be an incoherent cross relaxation effect.
In Fig. 12, the calculated nuclear Overhauser effect (nOe) enhancement is shown as a function of temperature for this system. The temperature dependence of the timescale was approximated using an Arrhenius relation, 94 with the activation energy, Ea, assumed to be 17.4 kJ mol −1 and the timescale at 300 K, τ(300 K) = 200 ps, though we note that these values are approximated from the solution state 89 and are likely somewhat different in this system. 90,95 Approximating the spectral density as a single timescale Lorentzian function, the nuclear Overhauser effect is then calculated as with ω0 being the 1 H Larmor frequency. Note that the spectral density function in principle has an orientation dependence; however, here this prefactor is omitted since it would cancel out in the calculation of the nOe. At a 1 H Larmor frequency of ν 0 = 600 MHz, it is expected that a positive nOe would be observed and hence a positive enhancement in our signal over this temperature range. That this is an incoherent effect is additionally supported by the 2D 1 H-1 H spin diffusion MAS NMR spectrum of β-AspAla [ Fig. 11(g)], where we observe negative cross peaks to other sites (specifically, Ala NH, Ala HA), which would not arise in the case of the inverse sign coherent spin diffusion mechanism.
While both homonuclear 96 ( 11 B) and heteronuclear nOes 97-100 ( 1 H-13 C, 1 H-15 N) have been observed before in the solid state, we believe that this is one of the first observations of a 1 H-1 H solid state nOe. Indeed, the only example we were able to find following an extensive literature review is the suggestion that a 1 H-1 H homonuclear nOe is responsible for the difference between 1 H-1 H SD and 1 H-1 H fp-RFDR spectra occurring in a study of bone at 110 kHz MAS. 101 However, in light of the results of this study, we would suggest that what was observed in those spectra was not actually a nuclear Overhauser effect, rather it was the result of resonance offset truncation in spin diffusion. While this paper was under review, a study of hydrogen-π interactions with trapped water molecules was published showing a similar nOe effect involving a methyl group. 102 The presence of this temperature-dependent incoherent effect may suggest that at fast spinning frequencies and high magnetic fields, the spin terms responsible for coherent spin evolution are sufficiently truncated by resonance offset such that there may be a competing influence of incoherent effects.

Example 2: Spin diffusion in 100% back-exchanged perdeuterated GB1
Since we have observed strong resonance offset dependence for spin diffusion at fast MAS in a protonated sample, we also wanted to investigate the extent of this phenomenon in a system with a more dilute proton network. For perdeuterated GB1 with 100% back-exchange, the combination of perdeuteration along with fast magic-angle spinning (at 60 kHz) leads to excellent resolution with proton linewidths smaller than 50 Hz. Selective saturation (∼100 Hz bandwidth, 500 ms) was applied to the LYS 10 N-H site, and the ratio of the intensities of all other resolved resonances under saturation on resonance and off resonance was recorded as a function of the resonance offset from the saturation transmitter. We observe a strong dependence of magnetization change on resonance offset [ Fig. 13(a)], yet negligible to no dependence on distance [ Fig. 13(b)]. Thus, under these conditions, considering the intensity of peaks as even a very rough indicator of spatial proximity may lead to wrong conclusions without careful consideration of the chemical shift. Additional discussion of this effect, including saturation of other sites within perdeuterated GB1, is given in Sec. S4 of the supplementary material.

V. CONCLUSIONS AND OUTLOOK
As faster magic-angle spinning and higher magnetic fields become routine for accessing higher resolution in 1 H solid-state NMR, understanding the behavior of spin diffusion becomes more important for designing experimental methods in which spin diffusion is utilized or where it is a competing effect to, e.g., chemical exchange like in CEST. Specifically, as 1 H sites become better resolved, the transfer of spin order between them via spin diffusion becomes increasingly dependent on their resonance offset. We showed that the inclusion of chemical shift evolution, and specifically chemical shift anisotropy, is of key importance for simulating how these systems evolve. We also investigated the effect of dynamics and found it to play an important role in the coherent evolution of the system. Finally, we observed that experimentally the situation is likely further complicated by the presence of incoherent cross relaxation effects.
While here we have only considered direct proton to proton spin diffusion, we note that the consequences and features observed likely affect other experiments. For example, CHHC 23,103 experiments as commonly applied at slower spinning rates will suffer the same resonance offset dependence under fast MAS as the proton spin bath relied upon for transfer becomes increasingly more resolved. This analysis highlights the importance of application of active recoupling under fast magic-angle spinning conditions both to collect spatial structural constraints [104][105][106][107][108][109][110] and in order to overwhelm the presence of incoherent cross relaxationbased nOe phenomena. On the other hand, better understanding of proton spin diffusion at fast MAS has implications for our ability to measure site-specific 1 H relaxation and quantification of 1 H CEST in the solid state.

SUPPLEMENTARY MATERIAL
See the supplementary material for additional discussion and data relating to the experimental phenomena observed here, as well as discussion of other considerations with regard to the model presented here, such as why the evolution of the spin system is marginally different for crystallographically equivalent sites, and on the number of crystallites that must be included.