Simulating spin dynamics in organic solids under heteronuclear decoupling Solid

Although considerable progress has been made in simulating the dynamics of multiple coupled nuclear spins, predicting the evolution of nuclear magnetisation in the presence of radio-frequency decoupling remains challenging. We use exact numerical simulations of the spin dynamics under simultaneous magic-angle spinning and RF decoupling to determine the extent to which numerical simulations can be used to predict the experimental performance of heteronuclear decoupling for the CW, TPPM and XiX sequences, using the methylene group of glycine as a model system. The signal decay times are shown to be strongly dependent on the largest spin order simulated. Unexpectedly large differences are observed between the dynamics with and without spin echoes. Qualitative trends are well reproduced by modestly sized spin system simulations, and the effects of ﬁ nite spin-system size can, in favourable cases, be mitigated by extrapolation. Quantitative prediction of the behaviour in complex parameter spaces is found, however, to be very challenging, suggesting that there are signi ﬁ cant limits to the role of numerical simulations in RF decoupling problems, even when specialist techniques, such as state-space restriction, are used. & 2015 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Effective decoupling of the 1 H nuclear spins is essential for achieving high resolution 13 C NMR spectra from typical organic molecules. Such heteronuclear decoupling is particularly difficult in the solid state due to the strong dipolar interactions between the different magnetic nuclei, which are not averaged out by molecular motion as they are in the solution state. While considerable progress has been made in developing approaches to decoupling and understanding how they work [1,2], there is not a comprehensive theory that allows decoupling performance to be predicted.
In principle, exact numerical simulation of nuclear spin systems [3][4][5] ought to allow the prediction of decoupling performance. The multi-spin nature of the dipolar-coupled network in solid systems is not necessarily an obstacle; we have, for example, shown that simulations on modest numbers (9-10) of spins are sufficient to predict exactly the 1 H spin dynamics under magic-angle spinning (MAS) [6,7], and numerical solutions have proved invaluable in understanding how decoupling sequences work [8][9][10][11][12][13][14][15]. We show here, however, that predicting the performance of decoupling sequences, particularly in regions of interest around optimal conditions, would require very large numbers of spins to be included to obtain quantitative agreement with experiment. This number is larger than the practical limit for exact simulation, which is typically less than 12 spins (although this can be extended using artificial model geometries with additional symmetry [16]).
Recently a number of groups have demonstrated simulations of much larger numbers of coupled nuclear spins by restricting the size of the state space used for the simulations [17][18][19][20][21]. Different researchers have used slightly different methods for restricting the evolution of the spin system to coherences below a certain order, but it is argued that the success of such calculations relies on the populations of higher spin orders (i.e. the number of correlated spins involved in a coherence) remaining relatively small [22]. This is clearly the case in solution-state NMR, where high spin order coherences relax relatively quickly, and some promising results have also been obtained for simulations of 1 H spin-diffusion in powder samples under MAS [23][24][25]. It is not obvious, however, that state-space restriction is generally appropriate in the solid Contents lists available at ScienceDirect journal homepage: www.elsevier.com/locate/ssnmr state, where very high spin orders can be observed amongst 1 H nuclei [26][27][28], and so we also investigate the role of higher spin orders in heteronuclear decoupling.

Methods
The decoupling performance is quantified experimentally by measuring the 'T 0 2 ' decay constant of the 13 C magnetisation, that is, the time constant for the decay of 13 C magnetisation where a 1801 pulse is applied on 13 C at the mid-point of the decay period [29,30]. The spin-echo refocuses decay due to inhomogeneous effects, such as B 0 inhomogeneity and magnetic susceptibility broadenings [31]. As a result, T 0 2 is much more sensitive to decoupling quality than T n 2 values obtained from measured linewidths; linewidths are often not particularly sensitive to changes in decoupling quality [30,32]. T 0 2 values continue to increase as the RF decoupling power is increased, tending towards the fundamental limit set by true T 2 relaxation [33], well after the limiting linewidth is observed. We use T c 2 here to refer to the coherent decay of 13 C magnetisation in the absence of a spinecho, to distinguish it from true (incoherent) T 2 relaxation and the overall time constant for magnetisation decay, T n 2 , which includes inhomogeneous contributions associated with the sample and any instrumental factors. The mechanisms for the decay of 13 C magnetisation are different for T c 2 and T 2 , but they both result in loss of the original coherence which cannot be readily refocussed. Decoherence is used here to refer to this magnetisation decay. In contrast to T 0 2 and T n 2 , T c 2 is not directly measurable, while both T 0 2 and T c 2 can, in principle, be directly observed in numerical simulations of magnetisation decay with or without a spin echo. It is important to note, however, that all the T 2 values used here are phenomenological quantities obtained by fitting experimental or simulated magnetisation decays to exponential functions. The absence of molecular tumbling in the solid state means that the magnetisation decays will generally be orientation dependent, and their powdered-averaged sum may not fit well to a single exponential.

Experimental
Experimental measurements of T 0 2 were performed on a polycrystalline sample of glycine-2-13 C, 15 N (99% 13 C, 98% 15 N) purchased from CortecNet. The sample was confirmed to be α-glycine based on the 13 C carbonyl peak at 176.5 ppm, which is sensitive to polymorphic changes [34,35]. As expected from the stability range of this form, 5-500 K [36] 13 C spectra of natural abundance samples under TPPM and XiX decoupling. The 13 C magnetisation was then measured after a spin-echo period, τ-π-τ, during which either CW, TPPM or XiX proton decoupling was applied, as shown in Fig. 1. For both the CW and TPPM experiments, the acquisition and recycle delay times were 30:77 ms and 4 s respectively, while under XiX decoupling they were 25:6 ms and 5 s respectively. Note that the TPPM pulse width, τ p , is generally parameterised below in terms of the corresponding nutation angle, θ ¼ τ p ν 1 3601.
The same 1 H decoupling was used in both spin-echo and acquisition periods. Although using a fixed decoupling sequence for acquisition would lead to reasonably consistent line-shapes in the acquired spectra, significant mismatches between spin-echo and acquisition decoupling were observed to distort fitted T 0 2 values via the orientation dependence of decoupling efficiency [38]; magnetisation that has been preserved by efficient decoupling during the spin-echo period may rapidly decohere under the acquisition decoupling, leading to an underestimate of intensity at longer spin-echo times and hence an underestimation of T 0 2 . The variation of the orientation dependence of T c 2 with decoupling parameters is illustrated in the Supplementary Information, Fig. S2. The 1 H decoupling nutation rate, ν 1 , was measured using the same sequence with a zero-length spin-echo period, incrementing the initial 1 H pulse width to acquire a 1 H nutation spectrum and taking the peak position as the nominal ν 1 .
Full decay curves were obtained at selected decoupling conditions by incrementing the evolution time, 2τ, linearly in 41 steps from zero to approximately twice the maximum expected T 0 2 . The free induction decays were zero-filled and Fourier transformed (without apodisation) using matNMR [39]. The decay of the methylene 13 C peak height as a function of 2τ was fitted to a decaying exponential to obtain T 0 2 using MATLAB s [40]. Where detailed parameter maps as a function of the parameters of a decoupling sequence were acquired, T 0 2 values were inferred from a pair of experiments at 2τ ¼0 and 2τ % T 0 2;max , assuming simple exponential decay of the peak height between these points. Discrepancies between the T 0 2 values obtained by this quick, but approximate, approach and those obtained from full decays were reduced by re-scaling the approximate values using a quadratic function fitted to approximate vs. accurate T 0 2 values at between three and five characteristic points in the parameter space. As illustrated in the SI, Fig. S1, this rescaling resulted in relatively modest changes in the T 0 2 values (up to 20% at maxima and 30% around minima), and allowed good T 0 2 estimates to be measured efficiently for a wide parameter space.

Numerical simulations
To explore how the spin dynamics change as a function of spinsystem size, spin-systems containing different numbers of protons at increasing distance from a selected methylene C atom were created, based on the room temperature neutron structure of α- glycine (CSD refcode GLYCIN20 [41]). These spin systems are labelled as CH n , with n indicating the number of protons in the system. CASTEP version 6.0 [42] was used to optimise the 1 H positions in the unit cell using a planewave cut-off energy of 600 eV. Brillouin zone integrals used a minimum sampling density of 0:1˚A À 1 apart with the sampling grid offset by 0.25, 0.25, and 0.25 in fractional coordinates of the reciprocal lattice. The exchange-correlation functional was approximated at the generalised-gradient level, specifically that of Perdew, Burke and Ernzerhof (PBE) [43]. Ultrasoft pseudopotentials [44] consistent with the PBE approximation were generated by CASTEP on-the-fly. Shielding tensors were subsequently calculated using the GIPAW method [45][46][47].
The effects of dynamics on the dipolar and shielding tensors of protons of the NH 3 þ groups, which are in rapid exchange at ambient temperature, were accounted for by averaging the chemical shift and dipolar coupling tensors over the three H positions and diagonalising to obtain the new principal components and mean tensor orientation. Dipolar coupling tensors between the spins of the NH 3 þ were reconstructed by re-orienting the averaged dipolar tensor along the C-NH 3 þ bond vector and scaling by P 2 ð cos 901Þ ¼ 1=2. This task of combining shielding tensor information from CASTEP and dipolar couplings determined from the geometry was handled with in-house software (available with the pNMRsim simulation program [48]). The dynamics of 1 H coupled networks are strongly determined by the root-sum-square of the 1 H dipolar couplings, d rss , at a given site [6], and so the contributions of neglected protons outside the extracted 'cluster' of spins to d rss were compensated for by scaling the 1 H homonuclear dipolar couplings so that the d rss at one of the methylene 1 H sites (H5 in GLYCIN20) of the reduced spin-system matched that of the extended lattice. This d rss value converges to 27:8 kHz when sufficient unit cells are considered (the value for the other methylene proton, H4, is very similar, 27:3 kHz). Note that d rss for H5 without motional averaging is 30:2 kHz. The heteronuclear dipolar couplings were not scaled since the heteronuclear couplings between C α and non-methylene protons have a negligible effect on the heteronuclear d rss values. Spin systems with unscaled 1 H homonuclear couplings were also created and used in indicated cases. The 1 H chemical shift referencing was chosen to bring the methylene protons on resonance by subtracting the calculated chemical shielding values from 26:56 ppm. The 13 C chemical shift and the negligible J couplings were not included in the spin systems. The resulting 13 C, ( 1 H) n spin systems are given the labels CH n below.
Simulations of RF decoupling under magic-angle spinning were performed in Hilbert space with pNMRsim [48], using a minimum time-step for propagator calculation of 1 μs. The theoretical background to such simulations has been extensively described elsewhere [49][50][51][52]. The simulations started with a state of 13 C x magnetisation and measured the remaining x magnetisation as a function of the duration of the decoupling period to create a simulated freeinduction decay (FID) or spin-echo decay. In spin-echo simulations, an ideal refocusing π-pulse [53] was applied at the mid-point of the rotation-synchronised decay time. Unless otherwise indicated, powder averages were performed over all three Euler angles describing the crystallite orientation, using 150 orientations distributed over a hemisphere generated with the ZCW algorithm [54][55][56]. Where the cycle times of the RF pulse sequence and sample spinning are not too dissimilar, it is generally possible to find a common time base for both the timing of the RF pulse sequence and MAS period. For phasemodulated RF pulse sequences, this allows the evolution of the density matrix to be determined from a limited number of propagators evaluated over a single period of rotation [11], greatly reducing the simulation time, usually by an order of magnitude or more. T 2 relaxation can be safely omitted from these simulations by noting that at room temperature the relaxation of the C α site of glycine is in the extreme narrowing limit where T 1 ¼ T 1ρ ¼ T 2 ; relaxation time constants on the order of seconds have been observed experimentally [57], much longer than the maximum T 0 2 observed for this site [33]. Although the 1 H T 1 is somewhat shorter (about 0.5-1 s), this is also orders of magnitude longer than the time constants for decay of the 1 H magnetisation due to "spin diffusion". The fast dynamics of the methyl group is helpful in shortening T 1 without contributing significantly to T 1ρ . When comparing time constants for coherent decay from simulation, T c 2 , with experimental T 0 2 values, it is important to take into account the inhomogeneity of the RF (B 1 ) field in typical NMR probes. The incorporation of RF inhomogeneity into the simulations is discussed in the SI.
Calculation times for the 10-spin CH 9 system required on average about one hour per orientation. As an example, the results for CH 9 data set shown in Fig. 7 were acquired on an institutional HPC cluster with the calculations for each powder orientation run in parallel on separate processors and required about 33,000 CPU hours. In contrast, the corresponding calculation times per orientation for the CH 6 systems were about three seconds. This is consistent with calculation times scaling as OðN 3 Þ, where N is the size of the Hilbert space.  Dash-dot lines correspond to TPPM decoupling at θ¼ 198.51, ϕ¼ 171, ν 1 ¼ 105 kHz and solid lines to CW decoupling at ν 1 ¼ 300 kHz. The spin systems are CH 3 (black), CH 6 (cyan), CH 9 (magenta) and CH 9 using only isotropic components of 1 H shifts (brown). (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this paper.) decoupling parameters have been chosen to result in similar rates of decay. It can be seen in Fig. 2(a) that physically reasonable decoherence is observed only when 1 H CSA parameters are included. This contrasts to the observation that the decay of 1 H magnetisation under simple MAS is essentially independent of the CSA parameters [7], but is consistent with the behaviour being largely determined by the second-order cross-terms between the heteronuclear dipolar couplings and 1 H CSA tensors [58,59,13]. The dynamics beyond the first few milliseconds depend on the size of the spin system; both the monoexponentiality of the decays and the observed T c 2 tend to increase and converge to sizeindependent limits as the number of spins increases. The increase in T c 2 with increasing spin-system size is consistent with increased "self-decoupling" in a larger dipolar-coupled network [60]. These trends are similar for the two sequences. Fig. 2(b) shows the decay of 13 C magnetisation as a function of the spin-echo time 2τ, but with otherwise identical simulation conditions. The behaviour observed is markedly different, with incomplete decoherence of the starting magnetisation that tends towards a plateau rather than decaying to zero. Although this unphysical behaviour is reduced in the larger spin systems, the convergence is much slower than observed without the spin echo. Simulations with a spin-echo pulse of varying tip angle showed behaviour which evolved smoothly between the limits of no spinecho and a full 1801 refocusing pulse, providing some reassurance that the effects observed are not artifacts of an over-idealised simulation. The possible origins of this unexpected behaviour, which was also observed in test simulations of static samples, are discussed below.

Simulation results
Analogous decay curves under XiX decoupling are shown in Fig. 3. As for CW and TPPM, the T c 2 tends to increase with increasing spin-system size, with the larger spin-systems showing clearly monoexponential decays. In this case, neglecting the CSA of 1 H shift tensors (dashed lines) has a much smaller effect, which is consistent with cross-terms between the heteronuclear and the homonuclear dipolar couplings being the limiting factor for XiX decoupling away from resonance conditions [61]. The spin-echo decays under XiX, Fig. 3(b), show similar unphysical incomplete decoherence, which is gradually reduced as more spins are included in the simulation, as previously observed for CW and TPPM.
The failure of the simulated spin-echo curves to decay to zero with increasing echo time is investigated in Fig. 4 using decay curves for CH 4 and CH 7 spin systems under CW decoupling for three non-special crystallite orientations. The first half, up to the π pulse at τ ¼ 25 ms, shows the 13 C magnetisation decaying towards zero as expected, although there is still significant oscillatory behaviour in the curves for individual orientations. These oscillations, associated with the finite size of the spin systems, are a function of the crystallite orientation, and so are effectively disguised by powder averaging, i.e. the success of the small spinsystem simulations in reproducing realistic magnetisation decays should not be overplayed. As can be seen from the individual curves, the π pulse has the effect of largely reversing the oscillatory evolution, and the final amplitude, which corresponds to the 2τ ¼ 50 ms point in the spin-echo simulations shown in Fig. 2(b), is almost completely refocused.
Further insight is provided by analysing the evolution of the density matrix during the spin-echo simulations. Fig. 5 shows the distribution of the matrix norm between the different 1 H spin orders as a function of time for a sample crystallite orientation under CW decoupling. The 1 H spin order for a given element of the density matrix, 〈ij σj j〉, corresponds to the number of 1 H spins that need to be flipped to convert the bra i to the match the ket j. Fig. 5(a) shows that the populated states initially all have a 1 H spin order of zero (corresponding to the starting state of pure 13 C magnetisation), but the higher spin orders are quickly populated due to "spin-diffusion" driven by the homonuclear 1 H couplings. The distribution of magnetisation immediately before the π pulse is close to statistical (i.e. reflecting the relative numbers of states of given order), but this evolution is mostly reversed by the π pulse. This behaviour is much more marked in the artificial system, Fig. 5(b), in which homonuclear couplings have been removed. Although the initial decay of spin order 0, due to the cross-term between the heteronuclear dipolar coupling and 1 H CSA, is the same in (a) and (b), the subsequent evolution shows much stronger oscillations, with the higher spin orders being populated much more slowly, and the starting 13 C magnetisation is fully refocussed. It is clear that the homonuclear couplings have a significant impact on the dynamics, but the underlying behaviour is still driven by the dipolar/CSA cross-term.
The unphysical refocusing of the magnetisation can only be an artifact of the finite size of the spin system, and has a significant impact on the utility of simulations based on finite spin-systems. Such simulations can only be guaranteed to reflect experimental observations at short timescales before the higher spin orders are populated and the "phase space" of the simulations has been filled. In the real spin-system, the density matrix norm can spread indefinitely into an infinite phase space, and a refocusing π pulse is unable to reverse this evolution. The essentially coherent nature of the evolution is less obvious in the absence of a refocusing pulse and the small system simulations are surprisingly effective [7], particularly when the evolution is averaged over multiple orientations.
The role of high spin orders in the spin dynamics raises the question of the extent to which decoupling performance is influenced by the parameters of distant spins. Fig. 6 shows simulated T c 2 values as a function of the TPPM pulse width using three model CH 8 spin systems: one (solid line) the normal glycinederived CH 8 system, the second (dash-dot line) with only isotropic chemical shifts on distant protons (further than the two methylene 1 H spins bonded to the 13 C), and the third (dashed line) without heteronuclear dipolar couplings to those distant protons. When T c 2 is small, i.e. the magnetisation decays quickly, the parameters of the remote spins have negligible impact. Around regions of peak decoupling, however, very different decay rates are observed. Unsurprisingly, given the mode of action of TPPM, including the CSA and heteronuclear dipolar couplings to the remote spins (solid line) results in poorer performance. Very similar effects of distant-spin parameters are observed across XiX parameter maps, except these are more sensitive to distant heteronuclear couplings rather than 1 H CSAs, see Fig. S3 in the SI. Although these distant heteronuclear couplings contribute negligibly to the total d rss , their impact on XiX decoherence times is significant. This makes it difficult to make quantitative predictions (a) CH 4 (b) CH 7  of peak decoupling since the dynamics clearly depend on both a large number of spin parameters and having a large numbers of spins i.e. it is not sufficient to 'rescale' the parameters of a small spin-system simulation.

Can decoupling performance be predicted?
These results demonstrate that decay of nuclear spin magnetisation under heteronuclear decoupling can only be effectively simulated using small spin-systems in a relatively narrow set of circumstances, for example, simulations of magnetisation decays under poor TPPM decoupling, Fig. 6, while simulating T 0 2 decays is even more challenging due to unphysical refocusing effects. As discussed in the introduction, state-space restriction methods [17,18,21,20,25] have recently been developed that allow much larger numbers of spins to be simulated, at the expense of neglecting higher spin orders. This is effective in solution-state NMR where the spin-spin couplings are relatively weak and highorder coherences relax quickly [22] the high order coherences do not build up sufficient population to have much impact on the spin dynamics. However, in general, this approach is not suitable for solids due to the strong couplings and slower relaxation. For glycine in the solid state we can infer that T 2 ¼ T 1 ¼ 5 s, as discussed previously. Given that 12-spins are about the limit for Hilbert-space simulation, we can estimate the maximum spinspin couplings that would allow accurate simulation of the spin dynamics for this relaxation rate. Following the procedure of Ref. [22], the couplings would have to be no more than 2.1 Hz if the relaxation rates are proportional to the coherence order, n, as in solution-state. This value is at least an order of magnitude smaller than the effective couplings observed in typical organic solids under MAS, as shown by the time constants for decay of the 13 C magnetisation of 1-2 ms observed in Fig. 5. We assume that the scaling of the relaxation rate with coherence order is the same in molecular solids and liquids based on experimental observations that relaxation in glycine is in the fast-motion limit [57], and noting that the dynamics in the solid state are not dissimilar to those in the solution state but just without translational degrees of freedom that are irrelevant to relaxation (with plastic crystals such as adamantane being an extreme example). If, however, we assume as Karabanov et al. [22] that the relaxation rates are proportional to ffiffiffi n p , then even smaller spin-spin couplings of less than 0.7 Hz would be needed for a realistic 12-spin simulation. See the SI for further discussion of relaxation in solids.
Nevertheless, restricted state-space simulations have had some success for solid-state simulations in the absence of RF decoupling [23,24]. Similarly "effective field" approaches have been used [60] to describe the effects of decay of coherence into the 1 H bath, but at the expense of introducing additional empirical parameters. We can somewhat crudely mimic these Liouville space simulations within the confines of a classical Hilbert space calculation by repeatedly nulling the amplitudes of higher spin orders in the density matrix, e.g. every rotor period. This might approximate the behaviour of an extended spin-system, with the magnetisation quickly evolving to higher spin orders and never returning. However, it was found (results not shown) that periodically nulling just the highest spin order had little effect on the unphysical refocusing. Nulling more spin orders progressively reduced the amount of refocusing, but, by effectively introducing relaxation to the decoherence, this artificially shortened the decay and so failed to reproduce the long-term evolution.
It is evident from previous literature [32,59,13] and the results above that simulations involving relatively few spins, 2-4, can reproduce the qualitative performance of decoupling sequences with respect to their parameters. However, the quantitative values of decay rates obtained from simulation depend strongly on the size of the spin system and the various parameters involved. Fig. 7(a) shows a cross-section of the TPPM parameter space at a fixed phase excursion as a function of spin system size. Larger than about CH 4 , the patterns are more-or-less consistent, which suggests it might be possible to "rescale" the results from calculations on a small system onto those obtained from much costlier multispin calculations. This is illustrated in Fig. 7(b) and (c), which plot the fitted T c 2 values from simulations of one spin-system against corresponding ones obtained using a different spinsystem. Although the decoherence times observed in CH 2 bear little correlation to those in CH 9 , there is a close-to-linear relationship between the results obtained from CH 6 and CH 9 systems. Given that each additional spin reduces the calculation efficiency by close to an order of magnitude, there seems little value in performing a calculation on a CH 9 system if very similar results can be obtained by rescaling (by a factor of γ¼1.58) results obtained on simulations of a CH 6 system. In principle, more complex analytical "transfer functions" could be used to perform this rescaling, but simple linear functions are adequate here. Fig. 8(a) show the results of plotting the correlations illustrated in Fig. 7(b) and (c), and extracting a transfer function gradient, γ, for different pairs of spin systems. Trivially, the gradient tends towards unity as the difference in the number of spins between the two systems decreases to zero. Simulations using fewer than six spins correlate poorly with larger simulations. On the other hand, when more distant 1 H spins are added, the transfer function gradients and correlation coefficients both steadily tend towards unity. The CH 5 system corresponds to the methylene carbon plus all the hydrogen atoms of a single glycine molecule, but it is not safe to assume that this explains the strong vs. erratic degree of correlation for spin systems that are larger vs. smaller than five 1 H spins. The number of spins necessary to start seeing simulation results representative of bulk behaviour is expected to depend on the decoupling performance in that region of the parameter space, with smaller numbers of spins needed to reproduce the T c 2 decay in regions of poor decoupling. The dashed line in Fig. 8(a) shows the mapping onto the CH 9 results where the homonuclear d rss has not been scaled to match the limit of an infinitely large system. Although the scaling factors are slightly different, the trends are the same, and so there is little advantage, in this case, of scaling the strength of the homonuclear couplings to match the extended lattice. It is worth noting that the transfer gradient increases as the homonuclear coupling network is strengthened in going from the unscaled to the scaled d rss , whereas it decreases when the size of the spin system increases. This suggests that it is the increased size of the Hilbert space, rather than the couplings themselves, that are responsible for the trends with increasing spin system size. The comparison of transfer function gradients for XiX, Fig. 8(b), shows the same overall trends as for TPPM, although there is a noticeable even-odd alternation of the gradient values, which is slowly damped as the size of the space increases. This has some similarities with the observation by Halse et al. of very strong even-odd alternation when varying the maximum allowed spin order in simulations of 1 H spin diffusion in a static solid [25]. Fig. 9(a) compares fitted T c 2 values obtained from the five largest spin-system simulations of Fig. 7(a), and the experimental T 0 2 recorded under the same conditions, as described in the Experimental section. The good agreement between the shape of the simulated and experimental parameter maps is reflected in Fig. 9(c), which shows the mapping between the simulated T c 2 decay constants to corresponding experimental T 0 2 values. Note how the very strong linear correlation between experimental T 0 2 and simulated T c 2 weakens as T c 2 increases beyond about 8 ms. Comparison with Fig. 5(a) suggests that this corresponds to the point where the finite size of the spin system has an increasing impact. The fact that the slope at smaller T c 2 is close to unity, however, is of limited significance; other slices through the parameter surfaces show the same trends, but with different slopes observed.
The comparison as a function of TPPM pulse phase, along a line passing through several local optima close to the theoretical optimum θ ¼ 1801= cos ðϕÞ [13], Fig. 9(b) and (d), is much poorer.
Both sets of transfer functions in Fig. 9(c) and (d) show the same pattern of convergence towards the experimental results, with largest deviations being observed for the longest decoherence times. The much weaker correlation observed in Fig. 9(d) is associated with the significantly greater complexity of this crosssection through the parameter surface; it includes regions of very poor decoupling, which are not well characterised by the choice of sampling point for the T 0 2 decays, combined with regions of good performance which require large spin-systems to be accurately described. Note that the corrections for RF inhomogeneity have little impact on the depth of the minima around ϕ ¼101 and 201, and so the discrepancies in this region are more likely to be associated with inaccurate measurements of very short T 0 2 values rather than deficiencies of the calculations. Similarly poor correlations when multiple resonance conditions are involved are also observed for XiX decoupling (Fig. S5 in the Supplementary  Information).
It is worth considering whether the disagreement in Fig. 9(d) is related to the differences between the experiments and the simulations. Firstly the experiments measure T 0 2 , using a spinecho to refocus inhomogeneous contributions to the decay rate (which are not present in simulation), but are compared to simulated T c 2 values, since T 0 2 simulations are complicated by unphysical behaviour, cf. Fig. 2. All the experimental evidence suggests that the positions of T c 2 and T 0 2 optima are the same [62]. This was confirmed in correlation plots using experimental peak height as the metric of experimental performance (Figs. S6 and S7 in the Supplementary Information), which show the same behaviour observed in Fig. 9. The simulations also assume that the initial 13 C magnetisation is uniform across all crystallite orientations. Simulations of the orientation-dependence of the CP efficiency, however, show that 1 H spin-diffusion during the contact times used (1.2 ms and 2.7 ms) results in close to uniform excitation. The final point of difference is that the simulations sample the total remaining 13 C magnetisation at the end of the decay period, but before acquisition, while the experimental measurement is based on the peak height of the methylene signal. As discussed in the Experimental section, using the same decoupling during spin-echo and acquisition periods minimises cases where magnetisation that has been retained during the decay period rapidly dephases and is not observed in acquisition period.

Concluding remarks
We have tested the fundamental limitations of exact simulations of the heteronuclear decoherence times T c 2 (directly observed FID) and T 0 2 (under a spin-echo) in the solid state under simultaneous MAS and RF decoupling. Counterintuitively, a spin-echo pulse on the observed nucleus is found to largely refocus the decay of 13 C magnetisation. Increasing the size of the spin system reduces the degree of refocusing, but it is not possible to reproduce experimental T 0 2 decays using the largest number of coupled spins that can be simulated exactly. Attempts to eliminate this unphysical refocusing by suppressing higher order coherences (mimicking simulations in which the state space is restricted to lower spin orders) were unsuccessful. Turning this around, however, we predict that liquid-crystalline materials, where the spin systems are limited to the liquid crystal molecules, would show this refocusing effect in spin-echo experiments.
In contrast, physically realistic T c 2 decays are obtained in simulation, but the nuclear spin decoherence around regions of good decoupling is found to be very sensitive to the size of the spin system and the details of the various spin-system parameters, including 1 H CSA parameters and/or dipolar couplings of distant spins. This makes quantitative prediction of peak decoupling performance particularly problematic. High spin-order coherences are rapidly populated, and while the experimental behaviour can be qualitatively reproduced by mimicking the effects of 1 H spin diffusion [1], it seems unlikely that simulations within a small Hilbert or Liouville space can quantitatively reproduce experimental behaviour without invoking adjustable and purely empirical parameters. This contrasts to earlier work that found that modest numbers (9-10) of 1 H spins were sufficient to predict the evolution of 1 H magnetisation under MAS [7] without the need for adjustable parameters. The difference is presumably related to the shorter timescales for the decoherence of the 1 H spin order compared to the much slower decay of 13 C magnetisation under RF decoupling.
Investigating the dependence of T c 2 decays on spin-system size, we find that the results for simulations in systems of 6-7 spins can often be usefully mapped on to those from larger, costlier calculations. Further increasing the size of the system leads to a monotonic rescaling of the decoherence times, and allows longer decoherence times to be observed at decoupling optima, but without significantly altering the map of performance vs. decoupling parameters. This confirms that an important determining factor in decoherence times is the size of "spin space" into which the magnetisation can spread. However, these 'transfer functions' tend to vary between different parts of the parameter space and between decoupling sequences, limiting their applicability to predict quantitatively the performance of arbitrary sequences. This also implies that state-space restriction techniques will similarly struggle to reproduce the quantitative behaviour of many-spin systems under RF decoupling.
Mapping of simulation results onto experimental behaviour is significantly less successful than mapping between simulation results in differently sized spin systems. As seen in Fig. 9, simulations correlate very well with experimental behaviour in smoothly varying regions of the parameter space (subject to the limitations of a finite simulation space discussed above). In regions where the parameter space is rapidly changing, however, there is not a straightforward mapping between experimental results and simulation. Considering Fig. 6, it is reasonable to suppose that the behaviour in rapidly changing regions of the parameter space is highly dependent on both long-and short-range NMR parameters as well as instrumental factors, such as RF inhomogeneity. In the   case of sequences such as XiX, the optimum decoupling conditions result from the interaction of multiple resonance conditions; quantitative prediction of decoupling performance in this situation is likely to be extremely challenging.
On the other hand, the optimum decoupling parameters are robustly reproduced using a modest number of spins, and the structure of the parameter maps is not sensitive to the parameters of remote spins. This is consistent with the routine experimental practice of optimising decoupling on a set-up system and using the same decoupling sequence parameters for the sample under study; it would not be possible to optimise experiments in this way if the positions of decoupling optima were strongly dependent on the parameters of multiple spins. Similarly, it is quite practical to perform multi-variable parametrisations of the performance of a decoupling sequence using a 6-7 spin system under different conditions by exploiting the efficient simulation techniques used here. While quantitative prediction of the performance of different local optima will remain a major challenge, this can be avoided in experimental practice by direct optimisation of the decoupling parameters [30].