Nuclear magnetic relaxation by the dipolar EMOR mechanism: Multi-spin systems

In aqueous systems with immobilized macromolecules, including biological tissues, the longitudinal spin relaxation of water protons is primarily induced by exchange-mediated orientational randomization (EMOR) of intraand intermolecular magnetic dipole-dipole couplings. Starting from the stochastic Liouville equation, we have previously developed a rigorous EMOR relaxation theory for dipole-coupled two-spin and three-spin systems. Here, we extend the stochastic Liouville theory to four-spin systems and use these exact results as a guide for constructing an approximate multi-spin theory, valid for spin systems of arbitrary size. This so-called generalized stochastic Redfield equation (GSRE) theory includes the effects of longitudinal-transverse cross-mode relaxation, which gives rise to an inverted step in the relaxation dispersion profile, and coherent spin mode transfer among solid-like spins, which may be regarded as generalized spin diffusion. The GSRE theory is compared to an existing theory, based on the extended Solomon equations, which does not incorporate these phenomena. Relaxation dispersion profiles are computed from the GSRE theory for systems of up to 16 protons, taken from protein crystal structures. These profiles span the range from the motional narrowing limit, where the coherent mode transfer plays a major role, to the ultra-slow motion limit, where the zero-field rate is closely related to the strong-collision limit of the dipolar relaxation rate. Although a quantitative analysis of experimental data is beyond the scope of this work, it is clear from the magnitude of the predicted relaxation rate and the shape of the relaxation dispersion profile that the dipolar EMOR mechanism is the principal cause of water-1H low-field longitudinal relaxation in aqueous systems of immobilized macromolecules, including soft biological tissues. The relaxation theory developed here therefore provides a basis for molecular-level interpretation of endogenous soft-tissue image contrast obtained by the emerging low-field magnetic resonance imaging techniques. © 2017 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). [http://dx.doi.org/10.1063/1.4991687]


I. INTRODUCTION
During the past two decades, a rigorous molecular theory has been developed for nuclear magnetic relaxation induced by exchange-modulated electric quadrupole [1][2][3] or magnetic dipole [4][5][6][7] couplings in aqueous systems with immobilized macromolecules. Such a theory is needed, for example, to interpret water 1 H or 2 H field-cycling (FC) 8-10 relaxation dispersion data in biophysical studies of water-protein interactions 3,[11][12][13][14][15][16][17] and to connect endogeneous soft-tissue image contrast to molecular-level phenomena in magnetic resonance imaging (MRI) analysis using emergent fast FC 18 or SQUIDdetected 19 low-field MRI techniques. For biophysical applications, the 2 H nuclide is the NMR probe of choice on account of the more straight-forward analysis of relaxation data involving the single-spin quadrupolar mechanism. But medical MRI applications invariably utilize the 1 H resonance, thereby forcing us to confront the more challenging multi-spin dipolar relaxation problem. a) bertil.halle@bpc.lu.se The present work completes a series of four papers devoted to the theory of dipolar relaxation by the mechanism of exchange-mediated orientational randomization (EMOR). In the three preceding papers, here referred to as Paper I, 5 II, 6 and III, 7 we treated two-spin and three-spin systems. Here, we make the leap to multi-spin systems, comprising one or two labile spins exchanging with an isotropic bulk phase and dipole-coupled to an arbitrary number of nonlabile spins in a solid-like environment.
In the EMOR mechanism, exchange plays a dual role: it transfers magnetization between the anisotropic macromolecular sites and the isotropic bulk phase and it drives relaxation by modulating the dipole couplings of the labile spin(s). The conventional Bloch-Wangsness-Redfield (BWR) perturbation theory of nuclear spin relaxation 20 therefore breaks down when, as is usually the case, fast-exchange conditions do not prevail. We have therefore developed a non-perturbative relaxation theory, based on the stochastic Liouville equation (SLE), 21,22 and valid without restrictions on exchange rate, dipole couplings, and magnetic field strength. [5][6][7] To better understand the rich relaxation behavior exhibited by the dipolar EMOR model, we have also developed a perturbation theory, based on the stochastic Redfield equation (SRE), that is, the BWR master equation supplemented with exchange terms, which, however, is valid only in the restricted motional-narrowing (RMN) regime, where the exchange rate exceeds all (fluctuating and static) dipole couplings.
For larger spin systems, the SLE theory becomes computationally intractable. We therefore turn to the more computationally efficient and physically transparent SRE theory, extending its validity beyond the RMN regime by certain physically inspired, but essentially ad hoc, modifications. These modifications were calibrated against the exact SLE solution of the EMOR model for four-spin systems, which is also presented here. The resulting generalized SRE (GSRE) theory, which is applicable in the full parameter space of the EMOR model, is then applied to spin systems comprising up to 16 protons extracted from crystal structures of two globular proteins. The GSRE theory is also compared to a previously proposed approximate theory, 4 based on the extended (by exchange terms) Solomon equations 23 (ESE). Unlike the GSRE and SLE theories, the ESE theory does not take into account longitudinal-transverse cross-mode relaxation or coherent transfer of magnetization (and higher spin modes) induced by the static dipole couplings.
The outline of this paper is as follows. In Sec. II, we briefly define the multi-spin dipolar EMOR model; more details can be found in Paper II. 6 In Sec. III, we develop the analytical multi-spin SRE theory, making extensive use of symmetry selection rules, and, in Sec. IV, we re-derive the multi-spin ESE theory, which, in the RMN regime, turns out to be a special case of the SRE theory. The exact solution of the SLE for four-spin systems, described in Sec. V, is used in Sec. VI to construct and assess the multi-spin GSRE theory for spin systems of increasing size. Finally, in Sec. VII, we discuss some limitations of the GSRE theory and how these can be overcome by suitable generalizations. The results presented in this paper are based on a substantial amount of analytical and numerical work, described in more detail in the 14 appendices of the supplementary material.

Spin system and EMOR mechanism
We consider an immobilized protein molecule containing one or two labile protons dipole-coupled to m nonlabile protein protons and exchanging with the protons in the surrounding bulk water phase. The case of a single labile proton, with spin I, might represent a hydroxyl or carboxyl proton in an amino acid side-chain. The case of two labile protons, with spins I and S, represents a water molecule transiently buried in a cavity inside the protein. For spins associated with nonlabile protons, we use the generic label P or lower-case Greek letters µ, ν, κ, λ, . . . .
To identify different spin systems and exchange cases, we use the notation "(spins in state A)-(spins in state B)," where state A comprises protons associated with the protein and state B includes the bulk water protons. The two generic multi-spin exchange cases treated here are thus denoted as IP m − I and ISP m −IS. The two-spin cases IS −IS and IP −I were treated in Papers I and II, 5,6 respectively, whereas the three-spin cases IP 2 −I and ISP −IS (as well as ISP −ISP) were considered in Paper III. 7 In the following, we present exact solutions for the four-spin cases IP 3 −I and ISP 2 −IS and approximate solutions for the generic multi-spin cases IP m −I and ISP m −IS without restriction on the number m of nonlabile protons.
The physical system that we have in mind is an aqueous protein gel or a soft biological tissue, where most or all protein molecules are effectively immobilized by chemical cross-links or non-covalent interactions. The system contains a large number of randomly oriented protein molecules, each of which is referred to as a site and identified by the generic label α. These sites, which make up state A, are identical in all respects except for their orientation (with respect to a lab-fixed frame), which conforms to an isotropic orientational distribution.
The labile protons undergo chemical or, in the case of internal water molecules, physical exchange with bulk water protons, and their mean survival time in the A sites is denoted by τ A . Detailed balance 24 then requires that P B τ A = P A τ B , where P A = 1 − P B is the fraction of all labile protons (including bulk water) that, at any given time, reside in A sites. In practice, 4,12,15 P A 1 and the following development will therefore focus on this dilute regime.
Exchange not only transfers magnetization between the A and B states and quenches multi-spin modes by fragmenting the spin system; 25 it is also the motion that induces spin relaxation by stochastically modulating all dipole couplings involving the labile spin(s) residing in A sites. In this dynamic model, which we refer to as exchange-mediated orientational randomization (EMOR), the orientation of all internuclear vectors involving the labile spin(s) is instantaneously randomized upon exchange. This simple model is justified if τ A is long compared to the time required for orientational randomization when the labile spin(s) has been transferred to state B. This is the case for chemical exchange of labile macromolecular protons with bulk water as well as for physical exchange of internal water molecules with bulk water. 5,26 In the interest of clarity, we have here described the model in terms relevant for the most important applications. However, the theory developed in the following is more general. For example, the macromolecules need not be proteins, the bulk phase need not be water, and the spin-1/2 nuclei need not be protons.

Spin Hamiltonian
From an NMR point of view, the A state has the peculiar property of combining solid-like and liquid-like features. Whereas I − µ, S − µ and I −S dipole couplings are averaged to zero by exchange, thereby inducing spin-lattice relaxation, µ − ν dipole couplings are static and therefore induce coherent evolution of the spin system. The dipolar spin Hamiltonian (in angular frequency units) for site α is where the first sum runs over all mutual dipole couplings among the m + 1 or m + 2 spins in site α. The irreducible spherical tensor operators (ISTOs) T 2 M (X) are normalized in the two-spin Liouville space of the spins involved in dipole coupling X, as in Table S1 of Paper II. 6 The argument of the rank-2 Wigner functions D 2 M0 (Ω α X ) is the Euler angles Ω α X that specify the orientation of the internuclear vector r X in site α with respect to the lab-fixed frame. 27 As in Papers I-III, [5][6][7] we incorporate a factor 3/2 in the definition of the dipole coupling: ω D,X ≡ (3/2) [µ 0 /(4 π)] γ 2 /r 3 X . As shown in Paper III, 7 the effect of proton chemical shifts is negligibly small under the conditions of interest. The spins can therefore be taken to be isochronous, so the Zeeman spin Hamiltonian is the same in all sites, where the S z term is present only for the ISP m −IS case. The Zeeman Hamiltonian in state B is the same as in state A, but without the sum over nonlabile spins. We omit the rapidly fluctuating dipolar Hamiltonian in state B, which gives rise to a small and frequency-independent relaxation contribution that, if so desired, can be added to the final expression for the relaxation rate. 5

B. Integral relaxation rate in the dilute regime
Our primary objective here is to calculate the frequencydependent longitudinal relaxation rate, R 1 (ω 0 ), of the abundant water proton spins. Under most conditions of interest, the system is in the dilute regime and relaxation is therefore strictly exponential. The so-called integral longitudinal relaxation rate (ILRR),R dil 1 (ω 0 ), that is, most readily obtained from theory, 1,6 is then identical to the observable R 1 (ω 0 ). For this reason, and in the interest of notational economy, we shall omit the caret as well as the "dil" superscript (indicating dilute regime conditions) used in previous papers [5][6][7] in this series. On the other hand, we will indicate with a subscript I or IS whether the relaxation rate pertains to the IP m −I or ISP m −IS case.
The general formalism for relaxation by the dipolar EMOR mechanism presented in Paper II 6 is valid for arbitrarily large spin systems. The relaxation rate in the dilute regime and with observation of the labile spin(s) is given, for the two exchange cases considered here, by R 1,IS = 2 P A τ A (U 11 + U 12 + U 21 + U 22 ) , where U np ≡ (n|U|p) is an element of the Liouville-space supermatrix U and the basis operators B n ≡ |n are taken to be ordered so that B 1 ∼ I z and B 2 ∼ S z . [The quantity here denoted by U corresponds to P A /τ A times the resolvent supermatrix U BB (s) at Laplace variable s = 0, as defined in Paper II. 6 ] According to Paper II, 6 in the dilute regime, where the exchange topology supermatrix T connects all B state basis operators (which only involve labile spins) with the corresponding operators in the larger A state basis, and T is the transpose of T. Furthermore, where the angular brackets indicate an isotropic orientational average (over all sites) and the supermatrix Λ α , associated with a particular site α, takes different forms depending on whether the subsequent theoretical development is based on the general stochastic Liouville equation (SLE) or its limiting form, the stochastic Redfield equation (SRE). 6,7 As seen from Eq. (3), for both exchange cases, the relaxation rate is strictly proportional to P A in the dilute regime.
In the multi-spin SLE theory, valid in the full parameter space of the EMOR model, Λ α is given by 6 where L α D and L Z are the Liouvillian supermatrices corresponding to the dipolar and Zeeman spin Hamiltonians in Eqs. (1) and (2), and the supermatrix K differs from the identity matrix only in that diagonal elements corresponding to basis operators that only involve nonlabile spins are zero. In the multi-spin SRE theory, valid in the restricted motional narrowing (RMN) regime, Λ α is instead given by 7 featuring the relaxation supermatrix R α and the coherent mode transfer supermatrix ∆ α , the explicit forms of which are given in Sec. III A. The supermatrix G A is isotropically averaged and therefore shares the cylindrical symmetry of the Zeeman Hamiltonian. The Wigner-Eckart theorem 28 then implies that G A is Q-block-diagonal in the ISTO basis for state A. The supermatrix T G A T is therefore Q-block-diagonal in the B state basis, as is the Zeeman Liouvillian supermatrix L Z in Eq. (4). Since the block structure is retained under matrix inversion, it follows from Eq. (4) that also U is Q-block-diagonal. As seen from Eq. (3), we only need elements from the Q = 0 block of U, so we can disregard the other Q-blocks. For the IP m −I case, the Q = 0 block is one-dimensional and [L Z ] Q=0 = L Z,11 = 0. For the ISP m −IS case, the Q = 0 block is 5-dimensional, but, if the labile spins I and S are isochronous in state B, we have [L Z ] Q=0 = 0. For both exchange cases, L Z can therefore be omitted from Eq. (4), whereby where T 0 connects B state Q = 0 operators with the corresponding A state operators. For exchange case IP m −I, so Eq. (8) yields where g np ≡ (n|G A |p). Combination of Eqs. (3a) and (10) then yields This is an exact result for the multi-spin IP m − I case in the dilute regime.
For exchange case ISP m − IS, the Q = 0 subspace is 5dimensional. However, to an excellent approximation, we need only retain the 2 × 2 block of G A spanned by the first two basis operators B 1 and B 2 , which have odd spin inversion conjugation parity. 29,30 As shown in Appendix A of the supplementary material, the relaxation rate then becomes R 1,IS = P A τ A 2 (1 − g 11 − g 22 + g 11 g 22 − g 12 g 21 ) (2 − g 11 − g 22 + g 12 + g 21 ) . (12) The test calculations reported in Appendix A of the supplementary material demonstrate that, in practice, this may be regarded as an exact result.

A. Relaxation and coherent mode transfer
The starting point for the multi-spin SRE theory is the restricted BWR master equation, which we derive in Appendix B of the supplementary material. As shown in Paper II, 6 this master equation leads to the exact results of Sec. II, notably Eqs. (5), (11), and (12) with Λ α given by Eq. (7). For isochronous spins, the relaxation supermatrix R α appearing in Eq. (7) is given explicitly by 7 The coherent mode transfer supermatrix ∆ α , which is present for m ≥ 2, is obtained as the obvious multi-spin generalization of the ISP−I result given in Paper III, 7 In Eq. (13), X and Y refer to dipole couplings that involve at least one labile spin, so they are randomized by the exchange. In Eq. (14), X refers to static dipole couplings between pairs of nonlabile spins. For the EMOR relaxation mechanism, the spectral density function (SDF) in Eq. (13) is of the form J(n ω 0 ) = τ A 1 + (n ω 0 τ A ) 2 . (15) As shown in Paper III, 7 the odd SDF has no effect for the ISP−IS case and a completely negligible effect for the IP 2 −I case. In Eq. (15), we therefore retain only the even part of the SDF. The angular functions in Eq. (13) are defined as To compute the Wigner functions D 2 * M0 (Ω α X ) in Eqs. (14) and (16), we need a convention for parametrizing the internuclear geometry. The Euler angles Ω α X = (ψ α X , ϑ α X , −) specify the orientation of the dipole vector r X with respect to the lab-fixed frame. First, we need a convention for the direction of the dipole vector r X . For fluctuating dipole vectors, r X is taken to point from the labile spin (I or S) to the nonlabile spin (µ) or, for the intramolecular water dipole coupling, from spin I to spin S. For static dipole couplings between nonlabile spins µ and ν, r X is taken to point from µ to ν with µ < ν.
The total number of (static or fluctuating) dipole couplings is m(m + 1)/2 for the IP m −I case and 1 + m(m + 3)/2 for the ISP m − IS case. A description that, for each site α, specifies the two angles ψ α X and ϑ α X for each of those dipole couplings is redundant. To obtain a more parsimonious description, we perform the transformation from the lab-fixed frame L, with the z L axis along the B 0 field, to the dipole-fixed frame X, with the z X axis along the dipole vector r X , via an intermediate sitefixed frame D. The first transformation step, which rotates the L frame into coincidence with the D frame of site α, is specified by the X-independent Euler angles Ω α LD = (φ α , θ α , ϕ α ). The second transformation step is specified by the site-independent Euler angles Ω DX = (γ X , β X , −), that is, the spherical polar angles that define the orientation of the dipole vector r X with respect to the D frame. In terms of the two successive transformations L → D → X, the Wigner functions can be expressed as In Appendix C of the supplementary material, we describe the conventions used to define the D frame and the Euler angles β X and γ X . The isotropic orientational average in Eq. (5) involves the Euler angles Ω α LD = (φ α , θ α , ϕ α ). For the angles φ α and θ α , which determine the orientation of the z D axis, we use Lebedev quadrature 31,32 with N L points on the unit sphere. For the angle ϕ α , which determines the orientation of the x D and y D axes, we use a uniform grid with N ϕ points. Unless otherwise noted, we use N Ω = N L × N ϕ = 14 × 5 = 70 Euler angle sets, but when extreme accuracy is desired we use N Ω = 50 × 15 = 750 sets.

B. Symmetry rules
The real-valued elements of the coefficient supermatrices in Eqs. (13) and (14) where the last equality follows from the cyclic permutation invariance of the trace. Whereas the ISTOs T 2 M (X) are normalized in the two-spin Liouville space corresponding to the two spins involved in the dipole coupling X, the basis operators B n are normalized in the multi-spin Liouville space, that is, Tr{B † n B p } = δ np with the trace including all m + 1 or m + 2 spins. The numerical prefactors in Eqs. (13) and (14) are consistent with these conventions.
The coefficients C XY MM ,np and D X M,np exhibit symmetries that allow us to formulate the SRE theory in a small part of the full (4 m+2 1)-dimensional spin Liouville space. The Liouville subspace notation used to describe these symmetries is defined in Appendix D of the supplementary material. For example, the subspaces L and N are spanned by basis operators involving only labile or nonlabile spins, respectively. The one-or two-dimensional subspace LZ is spanned by basis operators proportional to I z or S z , and the 27-dimensional subspace N 3 (µνκ) is spanned by basis operators involving only the three indicated nonlabile spins. These subspaces may be further decomposed into subspaces spanned exclusively by basis operators with odd (antisymmetric = a) or even (symmetric = s) spin inversion conjugation (SIC) parity. 29,30 For example, N = NA ⊕ NS. The odd and even subspaces can then be decomposed into subspaces with basis operators involving a specific number of distinct spins, for example, NA = N 1 ⊕ N 3 ⊕ N 5 if m = 6.
By using the general properties of the commutator and trace operations, we can establish useful symmetry rules for the coefficients C XY MM ,np and D X M,np in Eqs. (18) and (19). Here we merely state these rules; the proofs can be found in Appendix E of the supplementary material. The rules can be expressed most succinctly in terms of the relevant blocks of the supermatrices C XY MM or D X M . Symmetry Rule I. Labile single-spin modes are relaxationcoupled to nonlabile single-spin modes but not to other nonlabile odd-parity modes, and the relaxation coupling only involves self-correlations, Symmetry Rule II. Relaxation does not couple nonlabilespin modes involving different number of spins, Symmetry Rule III. Within the nonlabile k-spin subspace N k , relaxation only couples modes that involve the same set of k spins, C XY MM , N k (µνκλ··· )N k (µ ν κ λ ··· ) = δ µνκλ··· , µ ν κ λ ··· C XY MM , N k (µνκλ··· )N k (µνκλ··· ) .
The supermatrices C XY MM , N k N k and R α N k N k are therefore blockdiagonal in spins (k = 1), spin pairs (k = 2), etc. Furthermore, relaxation coupling between single-spin modes only involves self-correlations.
Symmetry Rule IV. Relaxation supermatrix elements within the single-spin subspace only involve self-correlations, Symmetry Rule V. Coherent dipolar evolution does not couple directly to labile-spin modes, Symmetry Rule VI. Coherent dipolar evolution does not couple nonlabile-spin modes with mixed modes, Together, symmetry rules V and VI show that coherent dipolar evolution can only couple nonlabile-spin modes to other nonlabile-spin modes.
Symmetry Rule VII. Coherent dipolar evolution couples nonlabile single-spin modes with nonlabile two-spin modes, but not with any other nonlabile k-spin modes, Symmetry Rule VIII. Coherent dipolar evolution only couples nonlabile spin modes differing by one spin. Thus, for For k = 1, only the plus sign applies and Eq. (27) reduces to Eq. (26). Symmetry Rule IX. Coherent dipolar evolution couples kspin modes with (k + 1)-spin modes only if k spins are shared between the modes. The mode coupling is induced by dipole couplings between the non-shared spin and each of the shared spins.

C. Integral relaxation rate
We have now defined all quantities in Eq. (7) that make up the supermatrix Λ α , which must then be inverted and isotropically averaged to obtain the elements of the supermatrix G A in Eq. (5) needed to calculate the relaxation rate from Eq. (11) or (12). Even though the dimension of the supermatrix Λ α is prohibitively large for realistic spin systems, a computationally manageable theory can be obtained by using the SIC symmetry of the supermatrices R α and ∆ α and the symmetry rules in Sec. III B. We shall thus generalize the three-spin treatment of Paper III 7 to the multi-spin cases IP m −I and ISP m −IS, without any limitation on the number m of nonlabile spins. This will be done exactly for m ≤ 3, whereas, for m ≥ 4, we neglect four-spin and higher spin modes. Here we only present an outline of the derivation of the multi-spin SRE theory; the full derivation can be found in Appendix F of the supplementary material.
Provided that the multi-spin ISTO basis for state A is ordered with the odd-parity operators before the even-parity operators, the SIC symmetry ensures that the relaxation supermatrix R α is block-diagonal whereas the coherent mode transfer supermatrix ∆ α is anti-block-diagonal. To make use of this symmetry, we partition Λ α into blocks associated with the antisymmetric (A) and symmetric (S) subspaces. We only need the AA block of the inverse (Λ α ) −1 . Expanding this inverse to first order in R α τ A , that is, to second order in ω D τ A , consistent with our use of the RMN approximation in Eq. (13), and performing the isotropic average in Eq. (5), the relevant 1 × 1 (IP m −I case) or 2 × 2 (ISP m −IS case) block of the supermatrix G A can be expressed in terms of the auto-spin relaxation supermatrix R α LZ,LZ and the cross-relaxation supermatrix Γ α LZ,LZ for site α, Using symmetry rules I and V, we obtain for the crossrelaxation supermatrix, The supermatrices X α N 1 N 1 and Y α N 1 N 1 describe the coherent transfer among single-spin modes associated with different nonlabile spins. This transfer proceeds via two-spin modes (described by X α N 1 N 1 ) and also three-spin and higher (up to m-spin) modes (described by Y α N 1 N 1 ), which simultaneously undergo relaxation induced by the fluctuating dipole couplings between the associated nonlabile spins and the labile spin(s).
The physical significance of the different parts of Eq. (29) may be appreciated with reference to Fig. 1, showing the pathways of dissipative (relaxation) and coherent (evolution under the static dipolar Hamiltonian) spin mode transfers for the IP 3 − I case. The 1 × 9 supermatrix R α LZ,N 1 describes the cross-spin relaxation between the labile spin I and the three nonlabile spins µ, ν, and κ in site α, including crossmode relaxation between longitudinal I z magnetization and transverse nonlabile-spin magnetizations. These pathways are indicated by dashed arrows in the lower part of Fig. 1. The 9 × 9 supermatrix R α , which is block-diagonal in spin (with three 3 × 3 blocks along the diagonal), describes the auto-spin relaxation of the single-nonlabile-spin modes (longitudinal and transverse magnetizations), as indicated by the three solid arrows pointing out from the corners of the triangle in the upper part of Fig. 1. The 9 × 9 supermatrix X α N 1 N 1 describes the coherent interconversion of single-nonlabile-spin modes via two-spin modes, corresponding to two consecutive solid wavy lines in Fig. 1. Finally, the 9 × 9 supermatrix Y α N 1 N 1 FIG. 1. Spin mode transfer pathways for the IP 3 − I case, showing B ↔ α exchange, auto-spin relaxation (solid arrows), and cross-spin relaxation (dashed arrows) in site α. In the upper part of the figure, wavy lines indicate the coherent transfer of among the single-spin modes of the three nonlabile spins in site α via two-spin modes (solid) or via the three-spin mode (dashed), as well as auto-relaxation (arrows) of all these modes. describes the coherent interconversion of single-nonlabilespin modes via two-spin and three-spin modes, for example, µ → µκ → µνκ → νκ → ν, involving two solid and two dashed wavy lines in Fig. 1. Simultaneously, the two-spin and three-spin modes undergo relaxation, as indicated by the thicker arrows in Fig. 1.
Expressions for the 3m × 3m supermatrices X α N 1 N 1 and Y α N 1 N 1 in terms of the relevant blocks of the relaxation and coherent mode transfer supermatrices R α and ∆ α can be obtained by using the symmetry rules of Sec. III B and expanding to second order in ω D τ A (consistent with the RMN approximation). For X α The supermatrix Y α N 1 N 1 involves coherent spin mode pathways of the general form (31) The mixing of the three-spin modes involves modes from twospin up to m-spin. To obtain a computationally tractable theory, we neglect k-spin modes with k ≥ 4. With this three-spin mode (3SM) approximation, we obtain Expressions suitable for machine computation of the relaxation rates R 1,I and R 1,IS are presented in Appendix G of the supplementary material. These expressions are based on Eqs. (11), (12), (28)- (30), and (32), but also make use of special properties, such as block-diagonality and Hermiticity, of the relaxation and mode transfer supermatrices. As required, these expressions reduce to the results presented in Paper III 7 for three-spin systems.
Explicit expressions for all relaxation supermatrices and coherent mode transfer supermatrices appearing in Eqs. (28)- (30) and (32) are given in Appendix H of the supplementary material. Whereas the single-spin mode relaxation supermatrices R α LZ,LZ , R α LZ,N 1 , and R α N 1 N 1 only involve self-correlations, that is, correlations of the dipole coupling between the same two spins at two different time points, the two-spin and threespin mode relaxation supermatrices R α N 2 N 2 and R α , appearing in Eqs. (30) and (32), also involve distinct correlations. Distinct correlations can be of two kinds: three-spin correlations, where one spin is shared between the two dipole couplings, and four-spin correlations, where four distinct spins are involved in the correlation. For the multi-spin EMOR cases considered here, a four-spin correlation must involve dipole couplings X = I µ and Y = Sν. However, it follows from Eq. (18) that such correlations do not contribute. The multispin SRE theory thus only involves two-spin and three-spin correlations.
The theory simplifies considerably if all static dipole couplings between nonlabile spins are neglected. Then the coherent mode transfer supermatrices ∆ α Computation of the cross-relaxation supermatrix Γ α LZ,LZ in Eq. (29) then only requires taking the inverse of the 3 × 3 blocks of the block-diagonal single-spin relaxation supermatrix R α

IV. MULTI-SPIN ESE THEORY
A simpler and less rigorous theoretical approach to the multi-spin EMOR problem for the IP m −I and ISP m −IS cases has been described. 4 The starting point for that analysis was the multi-spin extended Solomon equations (ESEs), describing the auto-spin and cross-spin relaxation of the longitudinal magnetizations of the labile and nonlabile spins and the exchange of the former between the anisotropic α sites and the isotropic B state. Being restricted to longitudinal magnetizations, the ESE theory cannot describe cross-mode relaxation and the inverted relaxation dispersion associated with it. 6,7 Moreover, because the orientation-dependent relaxation rates in the anisotropic α sites are replaced at the outset by isotropically averaged rates, the ESE theory disregards the dependence of the relaxation rate on the relative orientation of the dipole vectors, a dependence which, at least for small spin systems, is significant. 6,7 Finally, by neglecting static dipole couplings between nonlabile spins, the ESE theory omits all effects of the coherent mode transfer.
On the other hand, the ESE theory is more general than the SRE theory because it does not assume fast exchange (that is, ω D τ A need not be small), although it describes relaxation in the α sites with the aid of BWR theory (which presupposes that ω D τ A 1). The ESE theory is thus a hybrid approach, lacking internal consistency. Nevertheless, because the ESE theory is computationally expedient and at least approximately valid for arbitrarily slow exchange, it is of interest to compare it to the more rigorous, albeit more computationally demanding, theories developed here. Before doing so, we shall use the rigorous SRE framework established in Sec. III to re-derive the ESE theory. Apart from correcting a minor error in the original formulation, 4 this exercise allows us to precisely identify the approximations implicit in the ESE theory. The full derivation is relegated to Appendix I of the supplementary material; here we merely quote the final result.
According to the multi-spin ESE theory, the relaxation rate can be expressed in the familiar 33 dilute-regime form with the effective local longitudinal relaxation rate given by 4 with the SDF as in Eq. (15). Equations (33) and (34) are valid for both exchange cases, but with different meanings for the effective, frequency-dependent, dipole coupling Ω D . For the IP m −I case, with the cumulative dipole coupling ω D,I defined through For the ISP m − IS case, Ω D has a more complicated form (Appendix I of the supplementary material).

V. FOUR-SPIN SLE THEORY
The general theoretical framework for the exact multispin SLE theory, which is also valid outside the RMN regime, was established in Paper II. 6 Unfortunately, implementation of the SLE theory is not practically feasible for large spin systems because of the huge number of matrix elements that must be evaluated (by symbolic computer algebra) and because the required numerical computations (inversion of very large matrices at each of many Euler angle sets for each frequency point in the dispersion profile) become too timeconsuming.
In Paper III, 7 we implemented the SLE theory for the three-spin cases IP 2 − I and ISP − IS. Here, we describe the implementation of the SLE theory for the four-spin cases IP 3 −I and ISP 2 −IS. The motivation for this significant undertaking is to establish a benchmark and guide for the development of the more computationally efficient and physically transparent generalized multi-spin SRE theory described in Sec. VI. The three-spin cases are less suitable for this purpose, since they involve at most a single static dipole coupling. In contrast, the IP 3 − I case involves three static dipole couplings and it can therefore be expected to capture the essential features of larger multi-spin systems. As an additional benefit, the fourspin SLE solution can be compared to the SRE solution in the RMN regime, where both theories are exact, thereby providing a valuable check on the extensive algebra and computer coding underlying these calculations.
The 256 ISTOs spanning the four-spin Liouville space were generated by symbolic computer algebra using three successive angular momentum couplings, as described in Appendix J of the supplementary material. A complete list of these basis operators can be found in Table S3 of the supplementary material. The 65 025 elements of the supermatrix representation L α D in this basis of the dipolar Liouvillian corresponding to the dipolar Hamiltonian in Eq. (1) were generated as described in Appendix K of the supplementary material, again using symbolic computer algebra. The supermatrix representation L Z of the Zeeman Liouvillian corresponding to the isochronous Zeeman Hamiltonian in Eq. (2) is diagonal in the ISTO basis, with elements Q ω 0 . The supermatrix Λ α is then computed from Eq. (6), numerically inverted and isotropically averaged to yield the supermatrix G A in Eq. (5). Finally, the relaxation rate is calculated from Eq. (11) or (12).
Several checks were performed on the four-spin SLE solution. First, we demonstrated that it agrees quantitatively with the previously presented two-spin 6 and three-spin 7 SLE solutions in the special cases where one or two nonlabile spins are located far away from the labile spin(s). Second, we demonstrated that the SLE and SRE theories agree quantitatively in the RMN regime, where both theories are exact. This agreement is shown in Fig. 2 for the IP m −I and ISP m −IS cases. Here, we also include the corresponding results for two-spin and three-spin systems to highlight the substantial qualitative differences in the shape of the dispersion profile as more nonlabile spins are incorporated in the spin system.

VI. MULTI-SPIN GSRE THEORY
In this section, we develop a generalized SRE (GSRE) theory that reduces to the SRE theory of Sec. III in the RMN regime but also remains approximately valid outside the RMN regime, that is, for arbitrarily long τ A values. We consider first the multi-spin case IP m −I for m ≥ 3 and we then describe the few additional modifications needed for the ISP m −IS case.
The degree of success of the GSRE theory will be judged by comparison with results obtained from the exact SLE theory for the four-spin cases IP 3 − I and ISP 2 − IS. These test calculations were performed on 25 four-proton fragments (20 side-chain hydroxyl or carboxyl protons and five internal water molecules) extracted from crystal structures of the globular proteins, ubiquitin and BPTI, which have been studied extensively by water 1 H, 2 H, and 17 O MRD. 3,[34][35][36][37] The six dipole coupling constants ω D,X for each of the 25 fragments are listed in Tables S4 and S5 of Appendix L of the supplementary material. The nuclear coordinates from these structures were also used to compute the Euler angles β X and γ X (Appendix C of the supplementary material). In Sec. VI C, we use the GSRE theory to examine the scaling and convergence of the dispersion profile with an increasing number m of nonlabile spins. For these calculations, we used protein fragments with up to m = 15 nonlabile protons, selected in order of increasing distance from the labile spin(s). All calculations reported in this paper use the experimentally relevant value P A = 1 × 10 −3 . 3,12 However, as seen from Eqs. (11) and (12), as long as we are in the dilute regime, R 1 (ω 0 ) is strictly proportional to P A , so the value chosen is unimportant.

A. Exchange case IP m −I
For the two-spin cases IP−I and IS−IS, 5,6 the exact SLE theory leads to simple analytical results for the zero-field (ZF) relaxation rate R 1,I (0). For the IS−IS case, this result suggests that the validity of the SRE theory can be extended beyond the MN regime simply by replacing the SDF with a generalized SDF (GSDF), where the caret distinguishes the GSDF from the regular SDF in Eq. (15). For the IS−IS case, the resulting GSRE theory accurately reproduces the exact (SLE) dispersion profile R 1,I (ω 0 ) over the full τ A range and it is exact in the MN and ZF regimes. 5 A similar GSDF accurately describes the exact R 1,I (ω 0 ) profile for the quadrupolar spin-1 EMOR case. 2 For the IP−I case, the exact ZF rate R 1,I (0) is reproduced with a GSDF of the same form but with a factor 1/5 multiplying the (ω D τ A ) 2 term. 6 No attempt was made to reproduce the full dispersion profile with a GSRE theory for the IP −I case or for the three-spin cases treated in Paper III. 7 Although exact ZF results are not available for IP m −I cases with m > 1, physical considerations suggest a GSDF similar to that in Eq. (37), but with a different numerical coefficient in front of the (ω D τ A ) 2 term. Accordingly, we base the multi-spin GSRE theory on the GSDF, To assess the accuracy of this ansatz, we first consider R 1,I (0) for the IP m − I case without static dipole couplings (SDCs), which will be dealt with subsequently, whereby Eqs. (G.1)-(G.5) of the supplementary material and (38) yield Trial calculations of R 1,I (0) for the IP 3 − I case without SDCs show that the exact SLE result is reproduced to within a few percent over the full τ A range (Fig. 3, and Figs. S2 and S3 of Appendix M of the supplementary material), provided that the coefficient ζ I is allowed to vary from 1/3 near the MN regime to 2/15 in the ultraslow-motion (USM) regime, where (ω D,I τ A ) 2 1, according to The ZF rate predicted by Eqs. (39) and (40) coincides with the exact SLE result, not only in the MN limit, but also in the USM limit, where This is a particularly simple instance of the exact result, R 1,I (0) = C P A /τ A with a numerical constant C of order unity, which holds generally for the dilute EMOR model in the USM limit. Importantly, this result is valid for the dipolar 5-7 EMOR model for any m and with or without SDCs, as well as for the quadrupolar 1,2 EMOR model. For the two-spin dipolar EMOR model (the IP−I and IS−IS cases) and for the isomorphic axially symmetric quadrupolar I = 1 case (or any integral I), 1 An intuitive understanding of Eq. (41) can be obtained from the classical vector picture. The labile spin (or an ensemble of such spins) precesses (uninterrupted for many periods in the USM regime) with frequency ω D,I in the local magnetic field ω D,I /γ I produced by the m nonlabile spins in a particular site α, independent of the SDCs. A spin undergoing such local precession contributes a certain amount to the macroscopic longitudinal magnetization along any labfixed axis. After the spin leaves site α, it looses all correlation with the site orientation so that the next visited site can have any orientation drawn from an isotropic distribution. Consequently, the magnetization is randomized on the time scale τ A and we can regard 1/τ A as the local relaxation rate in the ZF USM limit. Note that the SDCs play no role here, consistent with our SLE results. This local relaxation rate is closely related to the so-called dipolar relaxation rate, for which a semi-phenomenological "strong-collision" theory (invoking the spin temperature concept) yields R 1,D = C D /τ c , where τ c is the correlation time for some local rotational motion and the numerical coefficient C D depends on what fraction of the dipolar energy is modulated by this motion. 38,39 The USM result in Eq. (41) can also be obtained from the ESE theory by taking the slow-exchange limit of Eq. (33), sincer τ A 1 in the USM limit. This is also true if the local relaxation is induced by a different motion than exchange (as in the EMOR model). Making use of the detailed balance condition, Eq. (41) can also be expressed as R 1,I (0) = 1/τ B , implying that spin I is relaxed as soon as it reaches any A site. However, in contrast to the GSRE theory, the ESE theory grossly overestimates R 1,I (0) at intermediate τ A values, and it is not even correct in the MN limit (Fig. 3, Figs. S2 and S3 of the supplementary material).
We now consider the full dispersion profile R 1,I (ω 0 ), but still without SDCs. The GSRE theory is based on Eqs. (G.1)-(G.5) of the supplementary material-but without the last two terms in Eq. (G.5) of the supplementary material, which vanish in the absence of SDCs-and the GSDF in Eq. (38) is used to calculate all relaxation supermatrix elements. This GSDF not only predicts R 1,I (0) accurately; it also describes how the frequency of the primary dispersion varies from ω 0 ≈ 1/τ A in the MN regime to ω 0 ≈ ζ 1/2 I ω D,I in the USM regime.
However, for ω D,I τ A 1, the GSRE profile also exhibits an inverted secondary dispersion step at ω 0 ≈ 1/τ A that is not present in the exact SLE profile (see Fig. 4 and Figs. S4-S7 of Appendix M of the supplementary material). The origin of this spurious GSRE dispersion is most readily appreciated in the USM regime. Well below the primary dispersion step, that is, for ω 0 ω D,I , the GSDF reduces toĴ(n ω 0 ) = J(0) = 15/(2 ω 2 D,I τ A ) so the elements of the auto-spin relaxation matrix R α N 1 (µ)N 1 (µ) in Eq. (G.5) of the supplementary material are of order 1/τ A . When ω 0 1/τ A , the term i ω 0 Q 1 in Eq. (G.5) of the supplementary material then cancels the off-diagonal cross-mode elements of R α N 1 (µ)N 1 (µ) , thus giving rise to the secondary dispersion step. To eliminate this spurious feature, we replace the Larmor frequency in the i ω 0 Q 1 term of Eq. (G.5) of the supplementary material, but not in the GSDF, by the renormalized frequency We refer to this modification of the SRE theory as suppression of nonsecular decoupling (SND). As seen from Fig. 4 (and Figs. S4-S7 of the supplementary material), a GSRE theory that features the GSDF in Eq. (38) as well as the SND modification in Eq. (42) reproduces the exact SLE profile quite well also outside the MN regime. The numerical factor of 2 in Eq. (42) was not rigorously optimized, but without it the GSRE profiles for τ A = 10 −4 s exhibit a small residual ND hump (more prominent for m > 3). In the USM limit (τ A = 10 −2 s), the ZF regime extends all the way up to the primary dispersion in the GSRE profile, whereas the SLE profile exhibits some fine structure in the low-frequency (LF) regime (1/τ A ω 0 ω D,I ). However, the difference between the two profiles is rather small (Figs. S6 and S7 of the supplementary material).
So far in this section, we have ignored the SDCs. What is their effect on the R 1,I (ω 0 ) profile? Judging from exact SLE calculations for the IP 3 − I case ( Fig. 5 and Fig. S8 of the supplementary material), the SDCs greatly increase the ZF rate R 1,I (0) in the RMN regime, but as τ A becomes longer their effect is diminished and, as we have seen, in the USM limit they have no effect at all on R 1,I (0). Moreover, the frequency of the primary dispersion is hardly affected by SDCs at any τ A value ( Fig. 5 and Fig. S8 of the supplementary material).
In contrast, the GSRE theory, as developed up to this point, predicts that the SDCs greatly increase R 1,I (ω 0 ) for all τ A values, also in the USM regime (Fig. 6). The origin of this spurious effect is the assumption in the derivation of the SRE theory that not only the fluctuating dipole couplings but also the SDCs are "motionally narrowed" in the sense that ω D,µν τ A 1. Formally, the spurious R 1,I (0) increase can be understood in terms of Eqs. (G.1)-(G.5) of the supplementary material. The ZF rate R 1,I (0) can be increased by decreasing the cross-relaxation rate Γ II zz (0) in either of the two ways: by increasing the nonlabile auto-spin relaxation rates or by increasing the rate of coherent mode transfer among the nonlabile spins, described by R α N 1 (µ)N 1 (µ) and [X α , respectively, in Eq. (G.5) of the supplementary material. These two phenomena are akin to incoherent and coherent spin decoupling, respectively. However, like the SRE theory itself, this picture is not valid outside the RMN regime. Indeed, according to the exact SLE theory, if the SDCs are artificially increased so that ω D,µν τ A 1, while keeping ω D,I τ A 1, then the SDCs actually decrease R 1,I (ω 0 ).
Outside the RMN regime, where ω D,µν τ A 1, we must therefore introduce a SDC correction in the GSRE theory. The finding that the (true) SDC effect vanishes in the USM limit suggests a renormalization such that the individual SDCs vanish when ω D,µν τ A 1. Accordingly, we replace the SDCs appearing in the coherent mode transfer supermatrices ∆ α N 1 N 2 and ∆ α N 2 N 3 in Eqs. (30) and (32) of the GSRE theory by the renormalized SDCs, FIG. 6. R 1,I (ω 0 ) dispersion profiles with τ A = 10 −6 s (a), 10 4 s (b), and 10 2 s (c) for labile protons in Thr-7 (left column) and Asp-39 (right column) side-chains of ubiquitin coupled to three nonlabile protons, computed from the SLE theory with SDCs (red, solid) and without SDCs (black, dasheddotted) and from the full GSRE theory (blue, dashed) and from GSRE theory without SND (green, dashed) or without SND and without SDC renormalization (magenta, dashed-dotted).
where the numerical coefficient was chosen to optimize the agreement with the exact SLE profiles for the IP 3 − I case.
In addition, the SND correction in Eq. (42) is also applied to the i ω 0 Q 2 and i L Z, N 3 N 3 terms of Eqs. (G.6)-(G.9) of the supplementary material to suppress nonsecular decoupling in the nonlabile two-spin mode and three-spin mode relaxation supermatrices. This additional SND correction has only a small effect for m = 3, but it is required for convergence of the dispersion profile with increasing m (Sec. VI C). As seen from Fig. 6 (and Figs. S9-S14 of Appendix M of the supplementary material), the GSRE theory, incorporating the renormalized quantities defined in Eqs. (38), (42), and (43), reproduces the exact R 1,I (ω 0 ) dispersion profile for the IP 3 −I case with good to excellent accuracy in the full parameter space of the EMOR model. In particular, the GSRE prediction of the R 1,I (ω 0 ) profile is virtually exact up to τ A ≈ 10 −6 s, R 1,I (0) is exact in the USM limit, and the dispersion frequency is close to the exact one under all conditions.

B. Exchange case ISP m −IS
The ISP m −IS case features three types of dipole couplings (I−µ, S−µ, and I−S) and two types of Euler angles (Ω Iµ and Ω Sµ ).
It is therefore not possible to derive a simple result analogous to Eq. (39). Nevertheless, we postulate a GSDF analogous to that in Eq. (38), but accounting for the three different types of dipole couplings, The function ζ I is given by Eq. (40) and ζ S by the analogous expression with ω D,I replaced by ω D,S . The function ζ IS is taken to be of the same form, but with different numerical coefficients, This choice of GSDF ensures that R 1,IS (0) = P A /τ A in the USM limit. In the (unrealistic) special case that ω D,IS (ω D,I , ω D,S ), the ISP 2 −IS case reduces to the IS −IS case, for which we have the exact ZF result ζ IS = 1, yielding the exact ZF USM result R 1,IS (0) = (2/3) P A /τ A . 5 Even though the I−S dipole coupling is dominant in the examined ISP 2 −IS cases (Table S4 of the supplementary material), the SLE theory yields R 1,IS (0) = P A /τ A exactly in the USM limit, just as for the IP 3 −I case. We assume that this result holds also for the ISP m −IS case with m > 2. To get the correct m-scaling in the GSRE theory, we must use a numerical coefficient in Eq. (45) that yields R 1,IS (0) = P A /τ A . With the coefficient 0.6425 in Eq. (45), the GSRE theory for the ISP 2 −IS case reproduces this USM result to better than 0.5% for all five internal water molecules.
As for the IP 3 −I case, the GSRE theory for the ISP 2 −IS case without SDC reproduces the exact SLE result for the ZF rate R 1,IS (0) to within a few percent over the full τ A range when Eqs. (44) and (45) are used for the GSDF (Fig. S15 of the supplementary material). In the MN regime, where the GSDF in Eq. (44) reduces to the regular SDF in Eq. (15), the (G)SRE profile agrees quantitatively with the SLE profile. This is true up to τ A ≈ 10 −6 s (Fig. S16 of the supplementary material). The dispersion shape at τ A = 10 −6 s exhibits a small (because of the dominance of ω D,IS ) inverted secondary dispersion at ω 0 ≈ ω 2 D,I τ A ≈ ω 2 D,S τ A (Fig. S16 of the supplementary material). The origin of this secondary dispersion step is nonsecular decoupling of longitudinal-transverse cross-mode relaxation in the nonlabile auto-spin relaxation supermatrix R α N 1 (µ)N 1 (µ) (ω 0 ), which does not involve ω D,IS , effected by the term i ω 0 Q 1 in Eq. (G.5) of the supplementary material.
In the USM regime, the primary dispersion step in the SLE and GSRE profiles occurs at ω 0 ≈ [0.6425 ω 2 D,IS + (2/15) ω 2 D,I + ω 2 D,S ] 1/2 rather than at ω 0 ≈ ω D,I (Fig. S17 of the supplementary material). However, the GSRE profile also exhibits a small inverted secondary dispersion step at ω 0 ≈ 1/τ A that is not present in the exact SLE profile (Fig. S17 of the supplementary material). To eliminate this spurious feature, we introduce a SND correction analogous to that in Eq. (42), As seen from Figs. S16 and S17 of the supplementary material, the GSRE profile based on the GSDF in Eq. (44) and the SND correction in Eq. (46) agrees rather well with the exact SLE profile without SDC, even though the fine structure in the USM regime is not captured by the GSRE theory. For the five ISP 2 −IS cases examined here, the single SDC is relatively weak (Table S4 of  coupling, the agreement is excellent also near the R 1,IS (0) maximum. This is a welcome result since internal water molecules often have τ A values in the range 10 6 -10 5 s. 3,16,17 Similarly to the IP 3 − I case, the GSRE theory is much more accurate than the ESE theory near the R 1,IS (0) maximum.
Applying also the SND correction in Eq. material, as for the IP 3 −I case, we find that the GSRE theory reproduces the exact R 1,IS (ω 0 ) dispersion profile for the ISP 2 − IS case with excellent accuracy in the full parameter space of the EMOR model ( Fig. 7 and Figs. S20 and S21 of the supplementary material).

C. Spin system scaling
Having calibrated the GSRE theory with the aid of exact SLE results for four-spin systems, we now examine the predictions of the GSRE theory for larger spin systems. We recall that the GSRE theory invokes the 3SM approximation for m ≥ 4 (Sec. III C). Specifically, we shall determine how, and at what rate, the dispersion profile converges as the number m of nonlabile spins increases. For the analysis of the IP m −I case, we focus on the labile protons in the side-chains of Thr-7, Thr-22, and Asp-39 in ubiquitin (Table S5 of the supplementary material), successively adding nonlabile protons up to m = 15 in order of increasing distance from the labile proton. Figure 8 shows how the cumulative dipole coupling squared, ω 2 D,I , increases with m for these labile protons.
We first consider the ZF rate R 1,I (0) (Fig. 9). At τ A = 10 −6 s, we are near the RMN regime, where, in the absence of SDCs, Eq. (39) predicts that R 1,I (0) scales with m as ω 2 D,I , so, at m = 15, R 1,I (0) has attained ∼95% of its converged value at m = 100 (as inferred from Fig. 8). In the presence of SDCs, R 1,I (0) converges even faster, albeit less smoothly, coming within 5% of convergence already at m ≈ 7. At τ A = 10 −4 s, we are approaching the USM limit, where R 1,I (0) = P A /τ A = 10 s 1 . Among the three examined labile protons, Thr-7 has the largest cumulative dipole coupling (Fig. 8) and is therefore closest to the USM limit. In the absence of SDCs, R 1,I (0) converges rapidly (∼95% convergence at m = 6) to a value ∼20% below the strict USM limit. In the presence of SDCs, R 1,I (0) converges more slowly to a value ∼30% above the true USM limit. At τ A = 10 −4 s, the SDC effect is modest for m = 3 ( Fig.  6), but, for larger m, the SDCs increase R 1,I (0) substantially ( Fig. 9). We now consider the convergence with increasing m of the entire R 1,I (ω 0 ) dispersion profile (Fig. 10) proportional to ω D,I (Figs. S26 and S27 of the supplementary material), and therefore converges at modest m values (as can be inferred from Fig. 8).
As expected on account of the dominant intramolecular I−S coupling (Appendix L of the supplementary material), both the ZF rate R 1,IS (0) and the shape of the R 1,IS (ω 0 ) dispersion profile converge faster with m for the ISP m −IS case than for the IP m −I case (Figs. S28 and S29 of the supplementary material).

VII. CONCLUDING REMARKS
This is the last one in a series of four papers developing the theory of longitudinal relaxation by the dipolar EMOR mechanism in systems comprising one or two spins exchanging with an isotropic bulk phase and dipole-coupled to a solid-like collection of m non-exchanging spins. The previous studies 5-7 provided exact solutions for two-spin and three-spin systems. Here, we have extended the exact theory to four spins and presented an approximate treatment-the GSRE theoryapplicable to arbitrarily large spin systems. For our Matlab implementation of the multi-spin GSRE theory, the computing time scales as m 9 . Fortunately, the relaxation rate converges rapidly with m; in practice, it is rarely necessary to go beyond m = 10, corresponding to a distance of 3.54.5 Å from the labile proton(s).
The multi-spin GSRE theory takes into account the transfer among the m nonlabile spins of magnetization, as well as of two-spin and three-spin modes, induced by the SDCs between these spins. This coherent process decisively influences the relaxation of the labile spin in the RMN regime, but the effect is attenuated as exchange becomes slower. In the more familiar so-called spin diffusion process, magnetization can be coherently transferred over large distances to a relaxation sink, which then drives the relaxation of the entire dipole-coupled spin system. 40 Spin diffusion from nonlabile to labile spins can greatly enhance the relaxation of the nonlabile spins, if their magnetization can be selectively detected. But when the magnetization of the labile spins is observed, as in a field-cycling experiment in the dilute regime, the magnetization starts at the relaxation sink ( Fig. 1), so long-range transfer of magnetization (and higher spin modes) is unimportant. In the context of the EMOR mechanism, the labile spin is the relaxation sink.
Whereas we have obtained exact solutions for small spin systems, up to and including the IP 3 − I and ISP 2 − IS cases, the multi-spin GSRE theory contains several approximations. The 3SM approximation neglects the coherent transfer involving four-spin and higher modes and is therefore invoked only for m ≥ 4. Given the limited range of coherent transfer in the EMOR context (see above) and the minor effect of SDCs outside the RMN regime, we do not expect the 3SM approximation to be a significant source of error. The GSRE theory was obtained by introducing three modifications designed to extend the validity of the SRE theory beyond the RMN regime: the GSDF and the renormalization of ω 0 and SDCs. Although inspired by theoretical considerations, these ad hoc modifications represent a trade-off between simplicity and accuracy. A more accurate GSRE theory could presumably be obtained at the expense of greater formal complexity. For example, to keep the theory simple, we have chosen to use a single GSDF, even though the multi-spin dipolar Hamiltonian has multiple eigenfrequencies. The non-axially-symmetric quadrupolar spin-1 EMOR theory involves three GSDFs, corresponding to the three NQR frequencies. 2 Since the dipolar Hamiltonian for a system of three or more spins also lacks axial symmetry, a GSRE theory featuring a single GSDF cannot be expected to achieve the same quantitative accuracy as for a two-spin system. 5 In 2006, one of us presented a more approximate relaxation theory for the dipolar EMOR mechanism in multi-spin systems in the dilute regime. 4 In Appendix I of supplementary material, we re-derived this so-called ESE theory in a way that clearly identifies the approximations involved. Although less rigorous than the GSRE theory, the ESE theory is less computationally demanding. It is therefore of interest to directly compare the predictions of the ESE and GSRE theories. This is done in Fig. 11 for a labile carboxyl side-chain proton and for an internal water molecule. For each case, we have chosen an experimentally relevant τ A value and included a sufficient number (m = 7 or 8) of nonlabile spins to approach convergence. For the carboxyl proton, the ESE profile differs surprisingly little from the GSRE profile, although the ESE theory fails to reproduce the low-field inverted dispersion step associated with nonsecular decoupling of longitudinal-transverse cross-mode relaxation. However, the near agreement in this case is largely a fortuitous result of error cancellation. Because the ESE theory does not include SDCs, it is more revealing to compare it with the GSRE profile calculated without SDCs. As seen from Fig. 11, the two theories then differ substantially. Also for the internal water case in Fig. 11, where SDCs play only a minor role, the ESE theory is seen to overestimate the relaxation rate significantly. The ESE theory was proposed with large spin systems in mind; for the two-spin and threespin systems treated in Papers II and III, 6,7 it is not a viable alternative.
There can be little doubt that the EMOR mechanism is the principal cause of water-1 H low-field longitudinal relaxation in aqueous systems of immobilized macromolecules, [12][13][14][15] including soft biological tissues. Nevertheless, the GSRE theory presented here ignores certain complications that should be addressed before the theory is used to quantitatively analyze experimental data. The two most obvious complications are as follows.
In the foregoing, we have tacitly assumed that the macromolecule is rigid, so exchange is the only motion capable of inducing spin relaxation. Conformational fluctuations tend to be fast (compared to τ A ) and, under typical conditions, their direct effect appears as an additive contribution to the labilespin relaxation rate. 2,5 The only significant effect of internal motion on the EMOR contribution is that all dipole couplings are multiplied by orientational order parameters. For simplicity, this has not been done for the sample calculations reported here, except for the internal water molecules (Appendix L of the supplementary material).
In the analyzed spin systems extracted from protein structures, all spins, except for the central labile spin(s), were treated as nonlabile, that is, non-exchanging. This assumption is justified for most of the cases analyzed here, where the ten nearest protons are bound to carbon or amide nitrogens. However, to correctly describe special cases with two adjacent polar side-chains or a cluster of internal water molecules, the theory needs to be modified. This complication can be approximately handled by associating an effective mean survival time τ A,Iµ = 1/(1/τ A,I + 1/τ A,µ ) to each I − µ dipole coupling. This simple amendment might suffice for labile protons in sidechain hydroxyl or carboxyl groups, as considered here. For the three protons in the ammonium group of a lysine side-chain, a more rigorous analysis is warranted, explicitly incorporating all three labile spins in the dipolar Hamiltonian and in the exchange topology supermatrix. However, at neutral or acidic pH, nitrogen-bound protons tend to exchange too slowly to contribute significantly by the EMOR mechanism. 12,37 Barring such complications, the total water-1 H relaxation rate in the dilute regime can be obtained as a population-weighted sum over all individual labile protons and internal water molecules. 2,5

SUPPLEMENTARY MATERIAL
See supplementary material for accuracy assessment of Eq. (12) (Appendix A); derivation of the restricted BWR equation (Appendix B); Euler angle conventions (Appendix C); Liouville subspace notation (Appendix D); proofs of symmetry rules (Appendix E); derivation of multi-spin SRE theory (Appendix F); explicit SRE results (Appendix G); relaxation and mode transfer matrix elements (Appendix H); derivation of ESE theory (Appendix I); table of four-spin ISTO basis operators (Appendix J); Liouvillian supermatrix elements (Appendix K); protein spin systems (Appendix L); accuracy

Supplemental Material
Nuclear magnetic relaxation by the dipolar EMOR mechanism: Multi-spin systems where q 3 , q 4 and q 5 are the sequence numbers in the A state multi-spin basis that correspond to basis operators 3, 4 and 5, respectively, in the B state two-spin basis, as defined in Table S2 of Paper III. 1 Combination of Eqs. (8) and (A.1) yields The five basis operators spanning the two-spin Q = 0 subspace are of two kinds: B 1 and B 2 are single-spin operators (proportional to I z and S z , respectively) with odd spin inversion conjugation (SIC) parity, 2, 3 whereas B q 3 , B q 4 and B q 5 are two-spin operators (involving both I-spin and S-spin operators) with even SIC parity. The SIC symmetry of the basis operators can be used to establish selection rules for the matrix elements g np , but only if the superoperator G A has definite SIC parity. In SLE theory this is not the case, because the superoperators L Z and L α D have different SIC parity. In terms of the supermatrices appearing in Λ α in Eq. (6), L Z is block-diagonal (or diagonal for isochronous spins) with respect to SIC parity whereas L D is anti-block-diagonal. Therefore, G A is neither block-diagonal nor anti-block-diagonal. In SRE theory, Λ α in Eq. (7) involves the supermatrices R α and L Z , which are block-diagonal with respect to SIC parity, and (for m ≥ 2) the supermatrix ∆ α , which is anti-block-diagonal. Therefore, just as in SLE theory, G A is neither block-diagonal nor anti-block-diagonal. Consequently, the full 5 × 5 matrix in Eq. (A.2) must be retained in an exact (SLE or SRE) treatment of the exchange case ISP m −IS.
In the absence of static dipole couplings, as for the three-spin ISP − IS case, the SRE Λ α matrix in Eq. (7) only involves block-diagonal (or diagonal) matrices, so also G A is block-diagonal with respect to SIC parity. In the MN regime, where SRE and SLE are equivalent, also the SLE-derived G A must be block-diagonal with respect to SIC parity. However, outside the MN regime, the SLE-derived G A is no longer blockdiagonal. Consequently, two conditions must be satisfied simultaneously for the mixed odd/even blocks in Eq. (A.2) to vanish: (1) there are no static dipole couplings, and (2) the fluctuating dipole couplings are in the MN regime (ω D τ A 1).

S3
In the zero-field (ZF) regime, the Wigner-Eckart theorem 4,5 implies that G A is blockdiagonal in the rank index K, as well as in the projection index Q (see Paper II 6 ). Apart from B 1 ∼ I z and B 2 ∼ S z , only one more of the five Q = 0 two-spin basis operators are of rank K = 1, namely B 4 ∼ i (I × S) · e z (with basis ordering as in Table S2 of Paper III 1 ). Consequently, in the ZF regime, Eq. (A.2) reduces to As we shall see, to a very good approximation, we can neglect the odd/even blocks in Eq. (A.2) so that U Q=0 becomes block-diagonal. In view of Eq. (3b), we only need to be concerned with the odd parity block,

Combination of Eqs. (3b) and (A.4) then yields
To assess the accuracy of this approximation, we use SLE theory for the ISP 2 −IS case to calculate For the ISP −IS case, where condition (1) is satisfied automatically, we can examine the consequences of violating condition (2) by comparing the ILRR computed from the standard SLE theory, corresponding to the full 5 × 5 matrix in Eq. (A.2), with the ILRR computed from a modified SLE theory where only the 2 × 2 odd-parity block is retained in Eq. (A.4). Within the MN regime, of course, the two methods yield the same ILRR. Outside the MN regime they do differ, but not by more than a few % and typically by less than 1 % over the full parameter space.

APPENDIX B: RESTRICTED BWR EQUATION
For completeness, we present here the derivation of the restricted BWR master equation, valid in the presence of static as well as fluctuating dipole couplings. This derivation closely follows the conventional one for the case with only fluctuating dipole couplings. 7 For clarity, we focus in this subsection on the exchange case IP m −I. However, the final result is valid also for the exchange case ISP m −IS.
For the exchange case IP m −I, where the spin system comprises a labile spin I and m nonlabile spins, the mutual dipole couplings are of two kinds: static (between two nonlabile spins) and fluctuating (involving the labile spin I). For site α, the corresponding Liouvillians are denoted L α D,0 and L α D (t), respectively. In a semiclassical description, the stochastic time-dependence of the latter is associated with the exchange of the labile spin I between an anisotropic site α and the isotropic bulk state. In addition, the Zeeman Liouvillian L Z describes the interaction of the m + 1 isochronous spins with the external magnetic field.
The evolution of the stochastic spin density operator for site α, ρ α (t), under these Liouvillians is governed by the Liouville -von Neumann equation, 7 The BWR master equation is obtained from a perturbation expansion, truncated after second order in the small quantity ω D τ A . Here, ω D is the typical magnitude of a fluctuating dipole coupling and τ A is the correlation time. Before performing this expansion, we must remove the coherent time-dependence from the density operator by transforming to the interaction representation according to In this interaction picture, the Liouville -von Neumann equation (B.1) becomes This differential equation is converted into an integral equation by a formal integration, which, when substituted into the right-hand side of Eq. (B.4), yields the exact result The time dependence of the fluctuating dipolar Liouvillian, L α D (t), is stochastic. Therefore, the time dependence of the spin density operator, ρ α (t), in Eq. (B.6) is also (partly) stochastic. The spin density operator related to the observable magnetization is obtained by taking an ensemble average, The ensemble may be thought of as a collection of trajectories, in each of which the labile spin I exchanges (but at different time points) from site α, thereby randomizing all dipole couplings that it is involved in. Neither the Zeeman Liouvillian L Z nor the static dipolar Liouvillian L D,0 is affected by this ensemble averaging.
Taking the ensemble average of Eq. (B.6), we obtain Ensemble averaging does not affect the initial spin density operator so the average in the first term on the right-hand side of Eq. (B.8) can be written where Eq. (B.2c) was used. The fluctuating dipolar Liouvillian is where c is a normalization constant. The ensemble average only involves the Wigner function, and with the EMOR propagator where τ A is the mean survival time of the labile spin I in site α. Consequently, Up to this point, the treatment is exact. We now introduce the motional narrowing (MN) approximation, stating that the time scale t for substantial variation in the observable of interest, that is, the longitudinal magnetization of the labile spin I, and the associated density matrix elements is much longer than the time scale τ A of fluctuations in the dipole couplings of the labile spin. Several simplifications follow from this assumption. First, since the exponential factor in Eq. (B.13) is 1 and, therefore, we can drop the first term in Eq. (B.8). Physically, this simplification results from the orientation Ω α Iµ (t) of the I −µ dipole couplings being completely randomized on the time scale t τ A . As a second consequence of the MN approximation, the integrand in Eq. (B.8) can be simplified as (B.14) Inserting this expression in Eq. (B.8), changing integration variable according to t = t − τ , and extending the upper integration limit to ∞ (a third consequence of the MN approximation), we find The master equation (B.15) in the interaction representation can be transformed back to the Schrödinger representation with the aid of Eq. (B.2), yielding Finally, we make use of the stationarity of the stochastic process to set t = 0 in the time correlation function in the integrand, whereby So far, we have only used the MN approximation, which is a condition on the dipole couplings with spin I: ω D,Iµ τ A 1. In the multi-spin SRE theory, we also impose a similar restriction on the static dipole couplings: ω D,µν τ A 1. As a consequence of the latter assumption, we can omit L α D,0 from Eq. (B.19) so that This simplification is possible because the essential contribution to the integral in Eq. (B.18) comes from τ values on the order of τ A or less, so that L α D,0 τ L α D,0 τ A 1. This argument only relies on the smallness of the numerical factor ω D,µν τ and is not invalidated by the fact that L Z and L α D,0 do not commute. This can be shown explicitly by Taylor expanding the exponentials in Eq. (B.19).
The master equation (B.17) predicts that the density operator σ α (t) evolves towards zero rather than towards the equilibrium density operator σ α eq . As usual, this deficiency is corrected by the ad hoc replacement of σ α (t) by the difference ∆σ α (t) ≡ σ α (t) − σ α eq in the last (relaxation) term of Eq. (B.17): Since σ α eq is independent of time, we can replace σ α (t) by ∆σ α (t) also on the left-hand side of Eq. (B.21). The same substitution can be made in the coherent term, because where β = 1/(k B T ) and Z = Tr{exp(−β H α 0 )} and last equality follows by Taylor expanding the exponential operator. We thus arrive at the restricted BWR master equation

APPENDIX C: EULER ANGLE CONVENTIONS
Here we delineate the conventions used to define the D frame and the Euler angles β X and γ X that specify the relative orientation of the dipole coupling vectors r X in any A site.

S10
For m ≥ 2, there are also dipole couplings between nonlabile spins. With the D frame defined as described above, Eq. (17) yields for these static dipole couplings where µ < ν and µν = 12 in Eq. (C.6b). Because the spins I, 1 and 2 all lie in the x α D −z α D plane, Eq. (C.6a) is of the same form as Eq. (C.1b) (that is, γ I2 = γ 12 = 0). Let r α µν be the vector from µ to ν with µ < ν and let e α µν = r α µν /r α µν be the corresponding unit vector. In analogy with Eq. (C.2), In analogy with Eq. (C.3), we have for µν = 12 To define the internuclear geometry of the IP m spin system with m ≥ 2, we need to specify m(m + 1)/2 distances r X and m(m + 1) − 4 internal angles (β X , γ X ). For m = 1, only one distance (r IP ) and no angles are needed.
For the ISP m −IS case, the vector r α IS pointing from spin I to spin S in site α defines the positive z α D axis, so that e α IS = e α z D and β IS = 0. The second vector needed to define the orientation of the D frame is taken to be r α I1 from spin I to the nearest nonlabile spin µ = 1. This vector is taken to lie in the x α D − z α D half-plane with x α D > 0, so that e α I1 · e α x D > 0 and γ I1 = γ S1 = 0. For the 2m dipole vectors emanating from spin I or spin S, Eq. (17) then yields with µ = 2, . . . , m in Eqs. (C.9d,e). In place of Eq. (C.2), we now have The azimuthal angles γ Iµ and γ Sµ are obtained from Eqs. (C.3) and (C.5) and their S-spin analogs (with I replaced by S everywhere). For the static dipole couplings between nonlabile spins (present when m ≥ 2), Eq. (C.6b) applies for all µν such that µ < ν (including µν = 12, since the nonlabile spin In place of Eq. (C.7), we have (for all µν) cos β µν = e α µν · e α IS , (C. 12) and the azimuthal angle γ µν is obtained from Eqs. (C.8).
To define the internuclear geometry of the ISP m spin system with m ≥ 1, we need to specify 1 + m(m + 3)/2 distances r X and m(m + 3) − 2 internal angles (β X , γ X ). For m = 0, only one distance (r IS ) and no angles are needed.
A slightly different convention was used for the three-spin system treated in Paper III, 1 where β I and β S were defined as the inner angles at spin I or S, respectively, in the triangle formed by the three spins I, S and P . For the IP 2 −I case, the correspondence is β I = β I2 and β S = π − β 12 . For the ISP −IS case, the correspondence is β I = β I1 and β S = π − β S1 . S12

APPENDIX D: LIOUVILLE SUBSPACE NOTATION
Here we define the notation for the spin Liouville subspaces used in the derivation of the SRE theory. For the IP m −I and ISP m −IS cases, the full spin Liouville space W can be written as the direct sum W = L ⊕ N ⊕ U, where L is the subspace spanned by basis operators involving only labile spins (I, or I and S), N is the subspace spanned by basis operators involving only nonlabile spins, and U is the remaining subspace, spanned by basis operators involving both labile and nonlabile spins (at least one of each).
Each of these subspaces, as well as the full Liouville space, may also be decomposed into subspaces spanned exclusively by basis operators with odd (antisymmetric = a) or even (symmetric = s) spin inversion conjugation (SIC) parity, that is, involving singlespin operators associated with an odd or even number of distinct spins, respectively. For example, W = A ⊕ S and N = NA ⊕ NS. The odd and even subspaces can be further decomposed into subspaces with basis operators involving a specific number of distinct spins. For example, if m = 6, then where NA is the subspace of odd nonlabile basis operators involving more than one spin and NS is the subspace of even nonlabile basis operators involving more than two spins. These generic subspaces may be further decomposed into subspaces involving specific spins. For example, Note that the order of the spins is irrelevant here: νµ and µν refer to the same subspace. For a spin system with m nonlabile spins, the number of distinct k-spin subspaces , that is, the number of distinct ways of choosing a subset of k spins, irrespective of their order, from a set of m spins. Each distinct k-spin subspace N k (µνκ · · · ) is spanned by 3 k basis operators, so it comprises 3 k k-spin modes. In other words, the dimension of N k (µνκ · · · ) is 3 k . The dimension of the total k-spin subspace N k , comprising all distinct combinations of k spins, For m = 3, for example, there are three single-spin subspaces - and N 2 (νκ) -each of dimension 9, and one three-spin subspace N 3 (µνκ) of dimension 27. The subscripts n and p denote basis operators B n and B p . To identify a basis operator belonging to a particular subspace, we give the subspace symbol as an argument. For example, the basis operator n(L 1 ) contains one single-spin operator associated with an unspecified labile spin, whereas n(L 1 (I)) = n(I) is one of the three basis operators (proportional to I z , I + or I − ) that span the I-spin Liouville space. The one-or twodimensional subspace spanned by basis operators proportional to I z or S z is denoted by LZ and the associated basis operators are n(LZ). Similarly, the basis operator n(N 1 ) contains one single-spin operator associated with an unspecified nonlabile spin, whereas n(N 1 (µ)) = n(µ) is one of the three basis operators (proportional to µ z , µ + or µ − ) that span the µ-spin Liouville space. Finally, the basis operator n(N 2 ) contains a product of two single-spin operators (or a sum of such products), each associated with a different nonlabile spin, whereas n(N 2 (µν)) = n(µν) is one of the nine basis operators that span the two-spin Liouville subspace associated with the specific nonlabile spins µ and ν (regardless of their order). The extension of these conventions to the nonlabile three-spin Liouville subspace N 3 is obvious.

APPENDIX E: PROOFS OF SYMMETRY RULES
Employing the subspace notation defined in Appendix D, we shall now prove nine symmetry rules for the coefficients C XY M M ,np and D X M,np in Eqs. (18) and (19). To do this, we shall make use of two general results. First, because two operators associated with different spins (such as I z and µ + ) necessarily commute, it follows that the commutator of two operators, each of which is a product of single-spin operators associated with distinct spins, vanishes if the two product operators have no spin in common. Second, since Tr I {I z } = Tr I {I + } = Tr I {I − } = 0, it follows that the multi-spin trace of a product operator, as in Eqs. (18) and (19), vanishes if the product operator contains at least one lone operator, that is, if any spin occurs only once in the product operator. Each symmetry rule will be derived from Eqs. (18) and (19) for basis operators B n and B p belonging to certain subspaces, but the symmetry rule will be expressed more succinctly in terms of the corresponding block of C XY M M or D X M . Symmetry Rule I. Labile single-spin modes are relaxation-coupled to nonlabile singlespin modes but not to other nonlabile odd-parity modes (involving 3 or 5 or ... nonlabile spins) and the relaxation coupling only involves self-correlations.
Proof. Assume that B n = n(L 1 ) = n(I). For the other possibility, B n = n(S), the following arguments are the same except that I and S are interchanged everywhere. The first commutator in Eq. (18) then vanishes unless the dipole coupling X involves spin I, that is, if X = Iµ or X = IS. Because B p = p(N k ) only involves nonlabile-spin operators, the second commutator in Eq. (18) vanishes unless the dipole coupling Y involves one nonlabile spin, that is, if Y = Iν or Y = Sν. To avoid having a lone I-spin or S-spin operator, which would cause the trace in Eq. (18) to vanish, we must have X = Iµ and Y = Iν. By assumption, B p is a product of an odd number of single-spin operators associated with distinct nonlabile spins (or a sum of such products). If B p involves 3 or more distinct nonlabile spins, the trace will vanish because its argument involves at least one lone nonlabile spin. For the trace to be nonzero, B p must involve only one nonlabile spin, which must also occur in the dipole couplings X and Y . In other words, X = Y = Iµ and B p = p(µ). Q.E.D.
Symmetry Rule II. Relaxation does not couple nonlabile-spin modes involving different numbers of spins.
Proof. Since B n and B p only involve nonlabile spins, the dipole couplings must be X = Iµ and Y = Iν or X = Sµ and Y = Sν and the nonlabile spins µ and ν must be contained in B n and B p , respectively. The trace over nonlabile spins then involves a product of k mutually distinct nonlabile-spin operators times a product of l mutually distinct nonlabilespin operators. If k = l, at least one lone nonlabile-spin operator occurs, causing the trace to vanish. Q.E.D.
Symmetry Rule III. Within the nonlabile k-spin subspace N k , relaxation only couples modes that involve the same set of k spins.
The matrices C XY M M , N k N k and R α N k N k are therefore block-diagonal in spins (k = 1), spin pairs (k = 2), etc. Furthermore, relaxation coupling between single-spin modes only involves self-correlations. Thus, Proof. For the trace to be nonzero, it is not only necessary that B n and B p involve the same number of nonlabile spins; actually, the identity of the nonlabile spins must be the same in B n and B p . Only then will all spins be pairwise matched so the trace can be nonzero. Moreover, for k = 1, the identifications X = Iµ and Y = Iµ , established in the proof of Eq. (E.2), imply when µ = µ , as dictated by Eq. (E.4a), that X = Y , so only self-correlations contribute to R α Symmetry Rule IV. Relaxation supermatrix elements within the single-spin subspace only involve self-correlations.
Proof. According to Eq. (18), C XY M M ,np can be nonzero only if B n and X share one spin (otherwise the first commutator vanishes) and if B p and Y share one spin (otherwise the second commutator vanishes). The argument of the trace in Eq. (18) is then a product of four single-spin operators associated with the four spins that are involved in the two dipole couplings X and Y . Because the (partial) trace of a lone single-spin operator vanishes, the four operators must consist of two operators associated with one spin and two operators associated with another spin. Since a dipole coupling involves two different spins, it follows that C XY M M ,np can be nonzero only if X = Y . Q.E.D.
Symmetry Rule V. Coherent dipolar evolution does not couple directly to labile-spin modes.
Proof. Since X = µν, the commutator in Eq. (19) vanishes if either B n or B p only involves labile-spin operators. Q.E.D.
Symmetry Rule VI. Coherent dipolar evolution does not couple nonlabile-spin modes with mixed modes.
Together, symmetry rules V and VI show that coherent dipolar evolution can only couple nonlabile-spin modes to other nonlabile-spin modes.
Proof. The trace in Eq. (19) vanishes because it involves at least one lone labile-spin operator, contributed by the basis operator from the U subspace. Q.E.D.
Symmetry Rule VII. Coherent dipolar evolution couples nonlabile single-spin modes with nonlabile two-spin modes, but not with any other nonlabile k-spin modes.
Proof. For B n = n(µ) or n(ν) and X = µν, the commutator yields µν. The only way to match these two operators without introducing any new lone-spin operators is if B p is a two-spin operator. In fact, it is necessary that B p = p(µν). Q.E.D.
Symmetry Rule VIII. Coherent dipolar evolution only couples nonlabile spin modes differing by one spin. Thus, for k ≥ 2, For k = 1, only the plus sign applies and Eq. (E.9) reduces to Eq. (E.8).
Proof. The commutator in the second form of Eq. (19) can be nonzero only if the basis operator B n = n(N k ) involves either or both of the spins in the dipole coupling X = µν.
In the former case, where B n = n(µ · · · ) does not involve spin ν or B n = n(ν · · · ) does not involve spin µ, the commutator yields µν · · · . The trace in Eq. (19) can then be nonzero only if B p = p(µν · · · ), that is, if B p contains the same spins as B n plus the second spin (here, ν) involved in the dipole coupling X. Hence, l = k + 1. In the latter case, where B n = n(µν · · · ) involves both spins in X, the commutator yields products of k + 1 operators associated with k spins, such as µ a µ b ν · · · and µν a ν b · · · . Since one spin is represented by two operators, the trace in Eq. (19) can be nonzero only if B p = p(µ · · · ) not involving spin ν or if B p = p(ν · · · ) not involving spin µ, that is, if B p contains the same spins as B n except one of the spins involved in the dipole coupling X. Hence, Symmetry Rule IX. Coherent dipolar evolution couples k-spin modes with (k + 1)-spin modes only if k spins are shared between the modes. The mode coupling is induced by dipole couplings between the non-shared spin and each of the shared spins. Proof. Implicit in the proofs of symmetry rules VII and VIII.
Special cases of symmetry rule IX are and (E.11)

APPENDIX F: MULTI-SPIN SRE THEORY
Here we present the full derivation of the multi-spin SRE theory. Specifically, we obtain the supermatrix G A , the elements of which determine the ILRR according to Eq. (11) or (12). First, we exploit spin inversion conjugation (SIC) symmetry. If the (4 m+2 − 1)dimensional multi-spin ISTO basis for state A is ordered with the odd-parity operators before the even-parity operators, the relaxation supermatrix R α is block-diagonal whereas the coherent mode transfer supermatrix ∆ α is anti-block-diagonal. Furthermore, K and L Z (for isochronous spins) are diagonal. We can therefore partition Λ α in Eq. (7) into blocks associated with the anti-symmetric (A) and symmetric (S) subspaces: To calculate the ILRR, we only need elements from the AA block of the inverse (Λ α ) −1 , We decompose the anti-symmetric subspace as A = L A + NA, where L = L + U. Noting that K NA,NA = 0, we can partition Eq. (F.2) as But we only need the L A, L A block of (Λ α ) −1 AA , This result is ill-suited for computations, because the L A subspace is very large for multi-spin systems. However, Eq. (F.5) can be simplified by implementing the RMN approximation. We thus expand the inverse in Eq. (F.5) to first order in R α τ A , that is, to second order in ω D τ A , consistent with our use in Eq. (13) of a relaxation supermatrix derived from BWR theory under the assumption that (ω D τ A ) 2 1 (for fluctuating and static dipole couplings). Performing also the isotropic orientational average, we then obtain from Eqs. (5) and (F.5), The L A subspace comprises all odd-parity basis operators that contain at least one labile-spin operator (I or S). But to obtain the ILRR, we only need matrix elements of G A L A,L A in the LZ subspace. In this subspace, It follows from Eqs. (14), (25) and According to Eqs. (13) and (20), R α LZ,NA is nonzero only in the nonlabile single-spin subspace N 1 , so Eq. (F.9) becomes Decomposing the NA subspace as N 1 + NA , we can write the inverse in Eq. (F.10) as According to Eqs. (F.10) and (F.11), We shall now use symmetry arguments and the RMN approximation to simplify the supermatrices X α , associated with spins µ and ν, can be expressed on a generic form involving matrix elements of the same form for all spins. To show this, we use Eq. (F.3) to write We now decompose the symmetric subspace as S = L S + NS, where L = L + U as before.
According to Eqs. (14), (26) and (27), Combination of Eqs. (F.23) -(F.25) yields where the upper limit of the sums is the largest integer contained in (m − 1)/2 or in m/2, respectively. Thus, for example, We decompose the even-parity subspace as S = NS + L S, with NS = N 2 + N 4 + . . . comprising the even parity subspaces containing only nonlabile spins and the remaining subspace L S = L + U including even-parity mixed labile/nonlabile subspaces (U 2 , U 4 , etc) and, if the spin system contains two labile spins, also the labile two-spin subspace The (NS, NS) block of the inverse of this matrix is We now apply the MN approximation again by truncating the expansion of (1+R α L S,L S τ A + i L Z,L S,L S τ A ) −1 after second order in ω D τ A , The last two terms are of order (ω D τ A ) 4 and (ω D τ A ) 6 , respectively, and they can therefore be neglected in the MN regime. Consequently, According to Eqs. (13) and (21), the supermatrix R α NS,NS is block-diagonal with matrices R α N k N k on the diagonal. Furthermore, L Z,NS,NS is purely diagonal. The matrix which is substituted into Eq. (F.26) to give, for k = 3, 5, . . . , Note that the second sum in Eq. (F.26), starting with k = 2, does not contribute to X α N 1 ,N k . In the same way, we obtain, for k = 3, 5, . . . , Now consider the matrix M α NA ,NA ≡ R α NA ,NA + X α NA ,NA + i L Z, NA ,NA . The Zeeman Liouvillian L Z, NA ,NA is completely diagonal, whereas R α NA ,NA is block-diagonal according to Eqs. (13) and (21). However, as we shall now show, X α NA ,NA is not block-diagonal. In analogy with Eq. (F.23), Using the selection rule for ∆ α S,NA in Eq. (F.25) and the corresponding relation for ∆ α NA ,S , we obtain for m ≥ 3 which is combined with Eq. (F.40) to give The supermatrix X α NA ,NA thus has a block-tridiagonal structure, as does the supermatrix M α NA ,NA in Eq. (F.38). To obtain the N 3 N 3 block of the inverse (M α NA ,NA ) −1 we must invert the full matrix (not just the diagonal blocks). The block-dimension of M α NA ,NA is [(m − 1)/2], that is, the smallest integer contained in (m − 1)/2. For m = 3 or m = 4, the supermatrix M α NA ,NA thus contains a single block, For m = 5 or m = 6, the supermatrix M α NA ,NA has a 2 × 2 block structure, where X α N 3 N 3 has the same form as for m = 4 (this form is valid for all m ≥ 4), Similarly, we obtain for m = 4, and, using also Eq. (F.48), for m = 5, and for µ = ν, (F.53) The supermatrix X α N 3 (µνκ)N 3 (µνκ) appearing in Eqs. (F.52) and (F.53) is given by (F.54) Here we have also introduced the 9 × 9 diagonal matrix Q 2 from Eq. (F. 19) and the analogous 27 × 27 matrix provided the 27 three-spin basis operators are ordered according to Q index, as in Table  S2 in Appendix H.
From the general result in Eq. (F.38), we see that, for any m ≥ 3, the coherent spin mode pathways are of the general form The mixing of the three-spin modes is brought about the matrix The block-tridiagonal matrix X α NA ,NA involves modes from two-spin up to m-spin. We now introduce the three-spin mode (3SM) approximation by neglecting k-spin modes with k ≥ 4. The 3SM approximation is thus needed only for m ≥ 4. The matrix M α NA ,NA then contains only one block, so that In the 3SM approximation, Eq. (F.49) applies also to m > 3, where the primed sums over spins are restricted to ensure that a k-spin subspace contains k distinct spins. For diagonal blocks (µ = ν), Equation (F.62) shows that the 27 × 27 blocks of X α N 3 N 3 vanish unless the two involved three-spin subspaces have at least two spins in common. For the nonzero off-diagonal blocks, connecting distinct subspaces sharing exactly two spins, only one of the three terms in Eq. (F.62) contributes. For the diagonal blocks, where the two subspaces are identical (sharing all three spins), all three terms in Eq. (F.62) contribute.
The effect on the ILRR of coherent spin mode transfer induced by static dipole couplings can be identified by computing the ILRR from the multi-spin SRE theory with and without the X α N 1 N 1 and Y α N 1 N 1 supermatrices, which vanish in the absence of static dipole couplings. Omitting these supermatrices from Eq. (F.13), we obtain The multi-spin SRE theory developed here reduces correctly to the results derived in Paper III for three-spin systems. 1 For the ISP − IS case, there are no static dipole couplings, so Eq. (F.63) applies. Setting m = 1 and N 1 (µ) = P , we find in agreement with Paper III. 1 For the IP 2 − I case, m = 2 and Y α N 1 N 1 = 0. From Eq. (F.13), we then obtain, in agreement with Paper III, 1

APPENDIX G: EXPLICIT SRE RESULTS
In this Appendix, we summarize the steps involved in computing the ILRR from the multi-spin SRE theory in the 3SM approximation (only used for m ≥ 4). In particular, we present explicit expressions for the special cases with m = 1, 2, 3 or 4, as well as general expressions, valid for any m, suitable for machine computation.
1. Exchange case IP m −I According to Eqs. (11) and (F.8), the ILRR for the IP m −I case is given exactly by The auto-relaxation rate R II zz is obtained from Eqs. (H.5) and (H.7a) as where the cumulative fluctuating dipole coupling ω D,I is defined as The cross-relaxation rate Γ II zz is obtained from Eqs. (F.13), (H.12a) and (H.13a) as .

S34
For m = 3, with a single three-spin subspace, the supermatrix X α N 3 N 3 consists of a single 27 × 27 block. For m = 4, with four three-spin subspaces, X α N 3 N 3 has 4 2 = 16 blocks given by Eq. (G.9). Each of the four diagonal blocks contains three terms, corresponding to the three two-spin subspaces contained in each three-spin subspace,

(G.22d)
For the off-diagonal blocks, only one of the three terms in Eq. (G.9) contributes, because there can only be one two-spin subspace contained in each of two different threespin subspaces (which must share two spins). Some of the 12 off-diagonal blocks for m = 4 are (the others are readily obtained by analogy), We now consider the general case, with m ≥ 3. The ILRR for the IP m −I case is given by Eqs. (G.1) -(G.5). In Eq. (G.4), we need to invert the 3m × 3m matrix V α N 1 N 1 , the constituent 3 × 3 blocks of which are given by Eq. (G.5). To obtain V α N 1 N 1 , we need to calculate m diagonal blocks R α N 1 (µ)N 1 (µ) , m 2 blocks X α N 1 (µ)N 1 (ν) , and m 2 blocks Y α N 1 (µ)N 1 (ν) . The matrices X α N 1 N 1 and Y α N 1 N 1 are neither symmetric, nor Hermitian. According to Eq. (G.6a), the off-diagonal (µ = ν) 3 × 3 blocks X α N 1 (µ)N 1 (ν) contain a single term, which involves the two-spin subspace N 2 (µν). The diagonal blocks X α N 1 (µ)N 1 (µ) contain m − 1 terms, corresponding to the two-spin subspaces involving spin µ and one of the remaining m − 1 spins. We can thus express Eq. (G.6b) as where N 2 (a µ ) denotes one of the m − 1 two-spin subspaces that includes spin µ. The 3 × 3 blocks Y α N 1 (µ)N 1 (ν) are given by Eqs. (G.7) -(G.9). Equations (G.7) and (G.8) may be expressed concisely as where N 2 (b µ ) denotes one of the two two-spin subspaces that contain spin µ and one of the other spins in the three-spin subspace This matrix inversion is the most computationally demanding step in the multi-spin SRE-3SM theory. For m = 3, 4, 5 and 6, the dimension of this matrix is 27, 108, 270 and 540, respectively. For large m, the dimension grows as m 3 . The three-spin-mode relaxation supermatrix R α N 3 N 3 is block-diagonal, with the 27 × 27 blocks computed as described in Appendix H. The three-spin Zeeman supermatrix L Z, N 3 N 3 is also block-diagonal, and all the blocks are equal to ω 0 Q 3 , with the 27 × 27 diagonal matrix Q 3 given by Eq. (F.55) if the three-spin basis operators are ordered as in Table S2 of Appendix H.

S36
The 27 × 27 blocks of the 3 → 2 → 3-spin mode transfer matrix X α N 3 N 3 are given by Eq. (G.9). The off-diagonal blocks may be expressed as where N 3 (A) and N 3 (B) are two different three-spin subspaces. The delta function δ c(A),c(B) signifies that X α N 3 (A)N 3 (B) vanishes unless these subspaces share two spins, in which case the two-spin subspace N 2 (c) comprises the two shared spins. For a given three-spin subspace N 3 (A), how many other three-spin subspaces N 3 (B) share exactly two spins with According to Eq. (G.9), the diagonal blocks of X α N 3 N 3 contain three terms, where c denotes one of the three pairs contained in A.

S38
Before considering the general case, we shall examine the special cases with m = 0 − 4. For m = 0, that is, the symmetric exchange case IS −IS, there are no nonlabile spins so the cross-relaxation rates in Eq. (G.31) vanish and the longitudinal auto-mode rates are given by Eq. (G.32) with ω D,I = ω D,S = 0.

(G.34d)
For m = 2, that is, the exchange case ISP 2 −IS, there is a single static dipole coupling, so Y α N 1 N 1 = 0 and the four 3×3 blocks of X α N 1 N 1 are given by Eq. (G.11), the only difference form the IP 2 −I case being that the two-spin relaxation matrix R α N 2 (12)N 2 (12) now involves S −µ as well as I −µ dipole couplings.
For m = 3 and m = 4, that is, the exchange cases ISP 3 −IS and ISP 4 −IS, X α and Y α N 1 N 1 are given by the expressions in Eqs. (G.12) -(G.23), the only difference form the IP 3 −I and IP 4 −I cases being that the two-spin and three-spin relaxation matrices R α N 2 N 2 and R α N 3 N 3 now involves S −µ as well as I −µ dipole couplings. In the general case, with m ≥ 3, the expressions for X α N 1 (µ)N 1 (ν) and Y α N 1 (µ)N 1 (ν) in Eqs. (G.24) -(G.29) are valid for the ISP m −IS case as well as for the IP m −I case. Again, the only difference is that the two-spin and three-spin auto-relaxation matrices R α N 2 N 2 and R α N 3 N 3 now involve S −µ as well as I −µ dipole couplings.

S39
APPENDIX H: R α AND ∆ α MATRIX ELEMENTS

Relaxation of single-spin modes
In this Appendix, we specify in explicit form all the matrix elements required in the multispin SRE theory, beginning with relaxation matrix elements between single-spin basis operators. According to Eq. (23), such matrix elements only involve self-correlations, that is, X = Y in Eq. (13). It also follows from Eq. (18) that for auto-spin rates, where n and p refer to the same (labile or nonlabile) spin, all (exchange-modulated) dipole couplings involving that spin contribute to the rate. Conversely, for cross-spin rates, where n and p refer to distinct (labile or nonlabile) spins, only the dipole coupling between those two spins contributes to the rate. Since single-spin relaxation matrices only involve self-correlations, the angular functions appearing in Eq. (13) and defined in Eq. (16) take the simpler form , we obtain the symmetry relations We consider first the isotropically averaged relaxation matrix in the labile single-spin longitudinal subspace LZ, which appears in Eq. (28), (H.6c)

S40
Here we have introduced the generic isotropically averaged longitudinal auto-spin and cross-spin relaxation rates associated with dipole coupling X, with the SDF J(n ω 0 ) given by Eq. (15). Next, we consider orientation-dependent local relaxation matrices in single-spin subspaces. These can all be expressed in terms of the generic auto-spin relaxation matrix associated with the dipole coupling X and the generic cross-spin relaxation matrix associated with the dipole coupling X A comment about indexing is in order here. In single-spin mode relaxation matrix elements, such as R IS zz or R νS +z , the left (row) basis operator is indicated by the left superscript and subscript, and the right (column) basis operator is indicated by the right superscript and subscript. In generic auto-spin and cross-spin rates, such as ρ Iµ z− or σ Sν +z , the subscripts are the same as in the corresponding relaxation matrix element (indicating row and column position in the matrix), but the superscript now refers to a dipole coupling X. For cross-spin rates, X involves the same two spins as the left and right basis vectors, but for auto-spin rates this is not always the case since several dipole couplings may contribute to a given auto-spin rate. In the generic rates, where the superscript refers to a dipole coupling X, the order of the two spins is irrelevant and we use the order corresponding to the convention adopted for the direction of the dipole vectors (Appendix C).
In terms of the angular functions in Eq. (H.1) and the SDF in Eq. (15), the elements of the generic auto-spin relaxation matrix are obtained from Eqs. (13) and (18) as Similarly, the elements of the generic cross-spin relaxation matrix are obtained as As indicated in Eqs. (H.9) and (H.11), the generic cross-spin relaxation matrix σ X is Hermitian. The generic auto-spin relaxation matrix ρ X is Hermitian only at ω 0 = 0.
The basis operators B n are normalized in the two-spin (µν) Liouville space.

Relaxation of two-spin modes
The 9 × 9 nonlabile two-spin relaxation matrices R α N 2 (µν)N 2 (µν) appearing in Eqs. (30) and (32) can be expressed as The elements ρ XY np of the 9 × 9 two-spin relaxation matrix ρ XY are all auto-(spin pair) rates, and either auto-mode (n = p) or cross-mode (n = p) rates. The nine two-spin modes are ordered as in Table S1. As noted in Sec. III.C of the main text, terms like ρ Iµ, Sν , which involve four-spin correlations, do not appear in the SRE theory.
Because the two-spin modes involve two different nonlabile spins (µ and ν), the twospin rates involve two different dipole couplings (X = Iµ and Y = Iν) for a given labile spin. In general, the two-spin rates therefore have contributions from self-correlations (s) as well as from distinct correlations (d), corresponding to the diagonal (X = Y ) and off-diagonal (X = Y ) terms, respectively, in the double sum of Eq. (13). Accordingly, we decompose the two-spin rates as where ρ d np (XY ) is the sum of the two terms in Eq. (13) that differ by X ↔ Y interchange. Consequently, the ordering of X and Y in ρ d np (XY ) is irrelevant. Within the Q-blocks (where Q n = Q p ; see Table S1), ρ XY np only has contributions from M = −M terms in Eq. (13). The Q-blocks are Hermitian, so ρ XY pn = (ρ XY np ) * , implying that ρ XY nn is real. The quantities ρ s nn (X) and ρ s nn (Y ) are identical, except for the X ↔ Y interchange. This is also true for off-diagonal elements ρ s np (X) and ρ s np (Y ) if the basis operators B n and B p have the same parity under spin permutation (Table S1). If the parity is different, ρ s np (X) and ρ s np (Y ) differ in sign (apart from the X ↔ Y interchange). Off-diagonal elements ρ XY np within the Q-blocks are real if B n and B p have the same spin permutation parity, and complex (for np = 45 and 67) or purely imaginary (for np = 13 and 23) otherwise. From Table S1, it is seen that basis operators with odd rank have odd parity under spin permutation.
The generic two-spin rates within the 3 × 3 Q = 0 block are obtained from Eqs. (13) and (18) as It may be noted that Eq. (16) For the 2 × 2 Q = 1 and Q = −1 blocks, we obtain in the same way Finally, for the Q = ±2 rates, we obtain ρ s 88 (X) = 1 9 ω 2 D,X 5 F 00 (X) Outside the Q-blocks, the Wigner-Eckart theorem requires that ρ XY np = 0 whenever |Q n − Q p | > 2, the rank of the ISTO T 2 M (X) in the dipole coupling. With the basis ordering of Table S1, the ten vanishing elements are np = 49, 94, 59, 95, 68, 86, 78, 87, 89, and 98. The remaining 81 − 19 − 10 = 52 nonzero off-Q elements have a more complicated structure since it is now the M = −M terms in Eq. (13) that contribute. Also, there is no simple relation between ρ XY np and ρ XY pn . As examples, we show only four of the 52 off-Q elements. First, for np = 14 or 41, with one basis operator from the Q = 0 subspace and the other from the Q = 1 subspace, we obtain from Eqs. (13) and (18), (H.31c) and ρ s 41 (X) = 1 18 ω 2 D,X − 3 Similarly, we obtain for np = 58 or 85, with one basis operator from the Q = 1 subspace and the other from the Q = 2 subspace,

Relaxation of three-spin modes
The 27 × 27 nonlabile three-spin relaxation matrices R α N 3 (µνκ)N 3 (µνκ) in the block-diagonal supermatrices R α N 3 N 3 appearing in Eqs. (F.60) and (F.61) can be expressed as The elements ρ XY Z np of the 27 × 27 three-spin relaxation matrix ρ XY Z are all auto-(spin triplet) rates, and either auto-mode (n = p) or cross-mode (n = p) rates. The 27 threespin modes are ordered as in Table S2.
Because the three-spin modes involve three different nonlabile spins (µ, ν and κ), the three-spin rates involve three different dipole couplings (X = Iµ, Y = Iν and Z = Iκ) for a given labile spin. In general, the three-spin rates therefore have contributions from self-correlations (s) as well as from distinct correlations (d), corresponding to diagonal S47 and off-diagonal terms, respectively, in the double sum of Eq. (13). Accordingly, we decompose the three-spin rates as where, for example, ρ d np (XY ) is the sum of the two terms in Eq. (13) that differ by X ↔ Y interchange.
a Parity of B n under spin permutation µ ↔ ν. b These basis operators have full spin permutation symmetry. c The basis operators B n are normalized in the three-spin (µνκ) Liouville space.

S48
The generic three-spin relaxation matrix ρ XY Z contains 27 2 = 729 elements. Among these, 141 elements are within the Q-blocks (where Q n = Q p ; see Table S2). Outside the Q-blocks, the Wigner-Eckart theorem requires that ρ XY np = 0 whenever |Q n − Q p | > 2. Consequently, there are 156 vanishing elements and 729 − 141 − 156 = 432 nonzero elements outside the Q-blocks.
The sequential coupling scheme used to construct the basis operators does not, in general, produce operators with definite parity under spin permutation. While all our basis operators have definite (odd or even) parity under µ ↔ ν interchange (see Tables S1 and S2), definite permutation parity involving other spins is only present in special cases. Among the 27 three-spin basis operators in Table S2, only eight have definite parity under all possible spin permutations: B 1 is odd whereas B 7 , B 13 , B 19 , B 22 , B 25 , B 26 and B 27 are even. Because our three-spin basis is not fully symmetry-adapted, the relations found among the elements of ρ XY only apply to ρ s np (X), ρ s np (Y ) and ρ d np (XY ), but not in general to ρ s np (Z), ρ d np (XZ) and ρ d np (Y Z), where the dipole coupling Z = Iκ involves spin κ.
To illustrate these considerations, we shall examine a few elements of the three-spin relaxation matrix ρ XY Z , calculated from Eqs. (13) and (18). Consider first matrix elements between fully permutation symmetric (odd or even) basis operators, In addition ρ XY Z 26,26 = ρ XY Z 27,27 . For these four diagonal elements between fully permutation invariant basis operators, As anticipated, these diagonal elements are real, and there is complete symmetry among the dipole couplings X, Y and Z.
In contrast, consider two diagonal matrix elements involving basis operators that only have µ ↔ ν symmetry, Also these diagonal elements are real, but now the components involving dipole coupling Z deviate from those involving only X and Y . Although all diagonal elements ρ XY Z nn are real, the Q-blocks of ρ XY Z are not Hermitian.
The matrix δ α N 1 (ν)N 2 (µν) is identical to δ α N 1 (µ)N 2 (µν) in Eq. (H.46), except for a sign reversal of all elements that involve a two-spin mode with odd permutation parity (and odd rank K). The parity is given in Table S1 and it is also reproduced on the line above the matrix in Eq. (H.46), showing that elements in columns 3, 4 and 6 are sign inverted in δ α N 1 (ν)N 2 (µν) . The two-spin to three-spin transfer matrices enter into Y α N 1 (µ)N 1 (ν) directly in Eqs. (F.60) and (F.61) and via X α N 3 N 3 in Eq. (F.62). The three distinct transfer matrices can be expressed in terms of six generic matrices as With the two-spin basis ordered as in Table S1 and the three-spin basis as in Table S2, Eqs. (14) and (19) yield the result in Eq. (H.48) for the dimensionless 9 × 27 matrix δ α N 2 (µν)N 3 (µνκ) (µκ) with D M ≡ D 2 * M 0 (Ω α µκ ), as given by Eq. (C.6) or (C.11). In Eq. (H.48), the numbering of the two-spin and three-spin basis operators is given to the left of and above the matrix, along with the µ ↔ ν permutation parity. The Wigner-Eckart theorem requires that δ α np = 0 whenever |Q n − Q p | > 2, the rank of the ISTO T 2 M (X) in Eq. (19). The 42 elements that vanish for this reason are shown in boldface in Eq. (H.48).
The matrix δ α N 2 (µν)N 3 (µνκ) (νκ) in Eq. (H.47a) is identical to δ α N 2 (µν)N 3 (µνκ) (µκ), except for a sign reversal of all elements that connect two-spin and three spin modes with different µ ↔ ν permutation parity. (Also the Euler angles are Ω α νκ instead of Ω α µκ .) The other four generic matrices in Eq. (H.47) are similar to the one shown in Eq. (H.48). However, the identity (apart from a sign reversal when the basis operators have different parity) of the pair of matrices applies to Eqs. (H.47b) and (H.47c) only for the eight columns associated with three-spin operators with full permutation symmetry (Table S2).

APPENDIX I: MULTI-SPIN ESE THEORY
In this Appendix, we use the SRE framework and well-defined approximations to derive the multi-spin ESE theory for the IP m −I and ISP m −IS cases. Our starting point is the multi-spin SRE framework as defined by Eqs. (5), (7), (11) - (19) of the main text. The RMN condition ω D τ A 1 will not be assumed except in so far that the relaxation supermatrix in Eq. (13) is based on BWR theory. Specifically, we will not implement the RMN approximation via series expansions, as we did in the derivation of G A LZ,LZ in Eqs. (28) -(32) (see Appendix F).

Single-spin mode truncation
The first approximation required to obtain the ESE theory is the single-spin mode (1SM) approximation, where all couplings between single-spin modes and multi-spin modes are neglected so the theory can be developed within a single-spin Liouville subspace of dimension (3m + 3) for the IP m − I case and (3m + 6) for the ISP m − IS case. Within the SLE framework, the 1SM is not a viable approach, because L α D = 0 in the single-spin subspace. This follows because all single-spin operators have odd SIC parity, but L α D has nonzero matrix elements only between odd and even basis operators.
The single-spin Liouville subspace can be partitioned into a (3 m)-dimensional nonlabile single-spin subspace N 1 and a 3-dimensional (for exchange case IP m − I) or 6dimensional (for exchange case ISP m − IS) labile single-spin subspace L 1 . Two consequences of the rather drastic 1SM approximation can be noted immediately. First, because labile two-spin operators are omitted, only the 2 × 2 single-spin block of the 5 × 5 matrix in Eq. (A.2) remains, so Eq. (12) is obtained without further approximation. Second, because all single-spin operators have odd SIC parity, ∆ α , which is anti-block-diagonal with respect to SIC parity, can be omitted from Eq. (7). The ESE theory, being based on the 1SM approximation, therefore cannot describe the secondary dispersion step at ω 0 ≈ ω D or any other effects of the static dipole couplings.

(I.1)
We need only the L 1 L 1 block of the inverse (Λ α ) −1 , According Eq. (E.4a), R α N 1 N 1 is block-diagonal in spin, with one 3 × 3 block for each of the m nonlabile spins, provided that the basis operators of the N 1 subspace are grouped by spin. Moreover, R α N 1 N 1 only involves self-correlations. If all spins are isochronous, Eq.

S53
(I.2) can therefore be expressed entirely in terms of 3 × 3 matrices as diag(0, 1, −1). Here, we also use the short-hand notation I ≡ L 1 (I) and S ≡ L 1 (S) for the 3-dimensional labile single-spin subspaces. Furthermore, where the sums run over the m nonlabile spins µ. The 3 × 3 relaxation matrices appearing in Eqs. (I.5) and (I.6) can be expressed in terms of the generic auto-spin and cross-spin relaxation matrices in Eqs. (H.8) and (H.9). For the cross-spin matrices, we thus have R α and for the ISP m −IS case, The calculation of the ILRR in the 1SM approximation involves the following steps. Given ω 0 , τ A , and the internuclear geometry, which defines the dipole couplings ω D,X and the Euler angles Ω DX = (γ X , β X , −), we compute the SDF J(n ω 0 ) from Eq. (15) and the angular functions F M M (X) from the expressions in Appendix C. These are then substituted into Eqs. (H.10), (H.11), (H.14), (I.7), and (I.8) to obtain the relevant autospin and cross-spin relaxation matrices. Combining these matrices with Eqs. (I.3) -(I.6), we obtain the 3 × 3 (IP m −I case) or 6 × 6 (ISP m −IS case) matrix [(Λ α ) −1 ] L 1 L 1 , the relevant elements of which are then isotropically averaged, as in Eq. (5), to obtain the matrix elements g np required to calculate the ILRR from Eq. (11) or (12). With the single-spin basis ordered by spin type (first spin I and then spin S for the ISP m −IS case), the matrix elements needed in Eq. (12) Note that the L 1 subspace contains all three components of the labile spin operators, whereas the longitudinal LZ subspace only contains the z component.
For the three-spin case ISP−IS in the MN regime, there are no static dipole couplings and the ILRR only involves single-spin modes. In the MN regime, the 1SM approximation is therefore exact for the ISP −IS case. In contrast, for m ≥ 2, including the IP 2 −I case, the 1SM approximation is not exact even in the MN regime.

Isotropic pre-averaging
The second approximation required to obtain the ESE theory amounts to replacing all local relaxation rates by their isotropic orientational averages, rather than performing this average at the end, as in Eq. (5). Averaging the angular functions as in Eq. (H.4), we find that Eq. (13) in the single-spin subspace (so that X = Y ) reduces to  with the effective local relaxation rater now given bỹ . It is clear, therefore, thatr approaches zero asymptotically, as expected.

S57
The multi-spin ESE theory developed here is equivalent to the theory previously obtained from the extended Solomon equations. 8 For the IP m −I case, the result in Eqs. (I.11) and (I.17) is essentially identical to the result in Eqs. [18] and [23] of the 2006 paper. The only (and trivial) difference is that, in the 2006 paper, we added the bulk relaxation rate and allowed for internal motions by multiplying all dipole couplings with orientational order parameters. For the ISP m −IS case, the 2006 paper suggests that, in Eq. (I.22), This is not the correct result for the model considered in the 2006 paper (and here). However, as seen from Eq. (I.23), it agrees with the correct result in the limit ω 0 τ A → ∞. Equation (I.24) is also exact for the special case of equilateral triangle geometry, where Ω 2 D (0) = (3/2) ω 2 D and Ω 2 D (∞) = 3 ω 2 D . Furthermore, in the USM limit, both the 2006 theory and the present corrected ESE theory yield R 1,I (0) = R 1,IS (0) = P A /τ A . However, for both exchange cases, neither the 2006 theory nor the present corrected ESE theory is correct in the MN regime, since isotropic pre-averaging eliminates cross-mode relaxation.

APPENDIX J: FOUR-SPIN OPERATOR BASIS
As a basis for four-spin Liouville space, we use the 2 8 = 256 irreducible spherical tensor operators (ISTOs) T K Q k I k S K k P {K}k U , constructed by three consecutive couplings of the set of four orthonormal single-spin ISTOs for each spin, e.g., to obtain 4, 9 whereK andK are the ranks of the intermediate tensor operators obtained by first coupling spins I and S and then spin P . In This Appendix, the fourth spin is denoted U .
Here, and in the following, the rank superscript is written in upper case for ISTOs that are normalized in four-spin Liouville space and in lower case for ISTOs that are normalized in single-spin Liouville space. The explicit forms of the 255 basis operators (excluding the identity operator) are listed in Table S3, ordered first by increasing projection index Q and then by increasing rank K. Note that the operators corresponding to the longitudinal magnetization of the I and S spins have sequence numbers 14 and 15 in Table S3, whereas they are assigned sequence numbers 1 and 2 elsewhere. All the operators in Table S3 are normalized in the same four-spin Liouville space, in the sense  049  (2)   123 10 (2K + 1)(2K + 1)(2K + 1)(2K + 1)(2K + 1)(2K + 1) (2k I + 1)(2k I + 1)(2k S + 1)(2k S + 1)(2k P + 1)(2k P + 1)(2k U + 1)(2k U + 1)

APPENDIX L: PROTEIN SPIN SYSTEMS
All SLE and GSRE calculations reported here pertain to protein-derived fragments of proton spins, with m nonlabile protons and one labile side-chain proton (IP m −I case) or one internal water molecule (ISP m −IS case). The nuclear coordinates were extracted either from the BPTI crystal structure 5PTI (A conformer) or from the ubiquitin crystal structure 1UBQ (H atoms added with WHATIF). Proton dipole coupling constants were obtained from the internuclear separation r X as ω D,X (rad For the IP m −I case, the labile proton I was chosen among 20 hydroxyl or carboxyl protons in the serine, threonine, tyrosine, aspartic acid or glutamic acid side-chains of ubiquitin (Table S5). No order parameters were used.
For the ISP m − IS case, the labile protons I and S were identified with the H 2 O protons of five internal water molecules in BPTI (W111, W112, W113 and W122) and ubiquitin (W128). The intramolecular I−S dipole coupling was set to ω D,IS = 2.46 × 10 5 rad s −1 . This value was estimated from the experimental H-H separation in liquid water, r IS = 1.545Å, and a typical orientational order parameter of 0.8 for the I-S vector. 10,11 For each fragment, all protons except the selected labile ones were regarded as nonlabile. Dipole couplings for the four-spin systems used to calibrate the approximate GSRE theory against the exact SLE theory are listed in Tables S4 and S5.

APPENDIX M: ACCURACY OF GSRE THEORY
In this Appendix, we compare the results of GSRE and SLE calculations on the four-spin systems listed in Tables S4 and S5 (Appendix L) to assess the accuracy of the GSRE theory for the IP 3 −I and ISP 2 −IS cases. We also include results from the ESE theory.
1. Exchange case IP 3 −I Figures S2 and S3 show that, for the 20 four-spin systems in Table S5 without SDCs, the simple GSRE result for R 1,I (0), based on Eqs. (39) and (40), remains very close to the exact R 1,I (0) computed from SLE theory over the full τ A range. While the GSRE result is exact in the MN and USM limits, it remains highly accurate also at intermediate τ A values. For example, at τ A = 10 −5 s, the GSRE result is off by less than 2.5 % in 16 out of 20 cases. In contrast, the ESE theory yields a much too large R 1,I (0) near the maximum and is not even correct in the MN limit. Figures S4 -S7 show that, for the 20 four-spin systems in Table S5 without SDCs, a GSRE theory that incorporates the GSDF in Eq. (38) as well as the SND correction in Eq. (42) removes the spurious secondary dispersion (evident in the green profiles, without SND correction) and brings the GSRE theory into good agreement with the exact SLE theory over the full frequency range and for all τ A values. It may be noted that, even though ω D,I τ A ≈ 400 at τ A = 10 −2 s for the labile protons with smallest ω D,I (most Asp and Glu residues), the SLE profiles are not fully in the USM limit, where R 1,I (0) = P A /τ A . In contrast, the GSRE profiles are in this limit, and therefore have a slightly larger R 1,I (0). For τ A = 10 −2 s, the ESE profile nearly coincides with the GSRE profile, although the primary dispersion is slightly upshifted and closer to the SLE profile.
The exact SLE dispersion profiles in Fig. S8, for the four-spin systems of Thr-7 and Asp-39 (Table S5), demonstrate how the effect of the SDCs is diminished as τ A becomes longer.
Figures S9 -S14 show that, for the 20 four-spin systems in Table S5, a GSRE theory that incorporates the GSDF in Eq. (38) as well as the SND correction in Eq. (42) and the SDC renormalization in Eq. (43), provides a good to excellent approximation to the exact SLE R 1,I (ω 0 ) dispersion profile in the full parameter space of the EMOR model. These figures also show the corresponding ESE profiles.
For τ A = 10 −6 s (Figs. S9 and S10), the GSRE theory, unlike the ESE theory, reproduces the hump in the profile, including its fine structure, and accurately describes the primary dispersion step. On the other hand, R 1,I (0) is too large by typically 5 -7 % for those spin systems that include a strong SDC (2.6 × 10 5 rad s −1 , for two protons in a CH 2 or CH 3 group). This discrepancy is a result of insufficient renormalization of these strong SDCs. By using a different functional form in Eq. (43), better agreement between the GSRE and SLE profiles at τ A = 10 −6 s can be achieved. However, the agreement will then deteriorate somewhat for τ A values above the R 1,I (0) maximum, which is the τ A range for most labile protons near neutral pH. Moreover, the renormalization function in Eq. (43) predicts that ω D,µν → 0 when ω D,µν → ∞, thus ensuring that R 1,I (0) = P A /τ A for any m ≥ 3. Without the square in the denominator of Eq. (43), ω D,µν tends to a finite value when ω D,µν → ∞ so the coherent mode transfer matrices X N 1 N 1 and Y N 1 N 1 do not vanish and R 1,I (0) increases with m in the USM limit, contrary to expectation.
For τ A = 10 −4 s (Figs. S11 and S12), the GSRE theory predicts a too small R 1,I (0) for most residues, by as much as 10 -15 % in some cases. However, the ESE theory overestimates R 1,I (0) even more. As for τ A = 10 −6 s, the GSRE theory accurately predicts the primary dispersion frequency. At τ A = 10 −5 s, the GSRE theory substantially overestimates R 1,I (0) for some residues, but it is still much better than the ESE theory.
For τ A = 10 −2 s (Figs. S13 and S14), in the USM regime, all three theories predict that R 1,I (0) = P A /τ A or very nearly so. The SLE profile exhibits a depression with fine structure in the LF regime, features that are not captured by the approximate theories. In most cases, at least the high-frequency part of the primary dispersion is accurately described by the approximate theories. 2. Exchange case ISP 2 −IS Figure S15 shows that, for the five four-spin systems in Table S4 without SDC, the ZF rate R 1,IS (0) computed from the GSRE theory based on the GSDF in Eq. (44) remains very close to the exact R 1,IS (0) computed from SLE theory. This is true over the full range of τ A values, from the MN regime to the USM regime. While the GSRE theory is exact in the MN and USM limits, it remains highly accurate even for τ A = 10 −5 s, erring by at most 2.1 %. Figures S16 and S17 show that, for the five four-spin systems in Table S4 without SDC, a GSRE theory that incorporates the GSDF in Eq. (44) as well as the SND correction in Eq. (46) removes the spurious secondary dispersion (evident in the green profiles, without SND correction) and brings the GSRE theory into good agreement with the exact SLE theory over the full frequency range and for all τ A values.
The exact SLE dispersion profiles in Fig. S18, for the internal water molecules W111 and W113 in BPTI coupled to two nonlabile protons, demonstrate the small effect of the SDC for the ISP 2 −IS case with a dominant I −S coupling. Figure S19 shows that the ZF rate R 1,IS (0) computed with the GSDF in Eq. (44) and the renormalized SDC in Eq. (43) almost coincides with the exact (SLE) result, which is not the case for the ESE theory. For τ A = 10 −6 s, the GSRE theory errs by at most 0.2 % for the five cases. Although not evident from Fig. S19, the GSRE theory is also highly accurate in the USM regime; for τ A = 10 −2 s the GSRE theory errs by at most 0.7 % for the five cases. Figures S20 and S21 show R 1,IS (ω 0 ) dispersion profiles for τ A = 10 −6 , 10 −5 , 10 −4 and 10 −2 s, comparing the predictions of the approximate GSRE and ESE theories with the exact SLE theory. For τ A = 10 −6 s, the GSRE theory is nearly exact and, unlike the ESE theory, it reproduces the tiny hump in the profile for W128. For τ A = 10 −5 s, the GSRE theory remains highly accurate, although the primary dispersion step is slightly stretched to low field. For τ A = 10 −4 s, the picture is much the same as for τ A = 10 −5 s, although the difference between the GSRE and ESE theories is not so large. For τ A = 10 −2 s, in the USM regime, all three theories predict that R 1,I (0) = P A /τ A or very nearly so. The SLE profile exhibits a depression with fine structure in the LF regime, features that are not captured by the approximate theories. In this Appendix, we present dispersion profiles, computed from the multi-spin GSRE theory, for some of the labile protons and internal water molecules listed in Tables S4 and  S5 (Appendix L), but now coupled with up to eight nonlabile protons. We also include results from the ESE theory.  Figure S28 shows the relatively fast convergence with m of the ZF rate R 1,IS (0) for the ISP m −IS case with m = 0 − 12, and Fig. S29 shows the similarly fast convergence of the shape of the R 1,IS (ω 0 ) dispersion profile for m = 0 − 7. In both figures, we use the internal water W122 in BPTI as an example. The other four internal water molecule in Table S4 give rise to similar results.