Electron paramagnetic resonance (EPR) spectra of molecular spin centers undergoing reorientational motion are commonly simulated using the stochastic Liouville equation (SLE) with a rigid-body hindered Brownian diffusion model. Current SLE theory applies to specific spin systems such as nitroxides and to high-symmetry orientational potentials. In this work, we extend the SLE theory to arbitrary spin systems with any number of spins and any type of spin Hamiltonian interaction term, as well as to arbitrarily complex orientational potentials. We also examine the limited accuracy of the frequency-to-field conversion used to obtain field-swept EPR spectra and present a more accurate approach. The extensions allow for the simulation of EPR spectra of all types of spin labels (nitroxides, copper^{2+}, and gadolinium^{3+}) attached to proteins in low-symmetry environments.

## I. INTRODUCTION

Continuous-wave electron paramagnetic resonance (CW EPR) spectra reveal important information about the structural and dynamic properties of paramagnetic spin centers. In particular, the spectral line shape can be highly sensitive to the nature and time scale of rotational dynamics of the spin center. An important application in this regard is the study of the dynamics of spin labels attached to soluble or membrane proteins [see Fig. 1(a)].

On the simplest level, the time scale of rotational dynamics of a spin center is characterized by its rotational correlation time, *τ*_{c}, which is related to the rotational diffusion rate constant *R* by *τ*_{c} = 1/6*R*. The shape of the CW EPR spectrum depends on the relation of *τ*_{c} to the width of the rigid-limit spectrum, Δ*ω*. If *τ*_{c}Δ*ω* ≫ 1, then the rotational motion has very little effect on the spin dynamics, the spectrum resembles the one in the immobile limit (*τ*_{c} → *∞*), and the rotational motion can be neglected for spectral simulations. If, on the other hand, *τ*_{c}Δ*ω* ≤ 1, then the motion is fast enough to mostly average out all anisotropies, and the observed spectrum is a sum of individual lines, similar to the one in the isotropic limit (*τ*_{c} = 0). In this fast-motion regime, spectra can be simulated by using an isotropic Hamiltonian and treating the rotational motion as a perturbation. The intermediate regime, where approximately 1 < *τ*_{c}Δ*ω* < 100, is called the slow-motion regime. In this regime, the rotational motion and the spin dynamics are strongly coupled, and the spectrum is sensitive to the details of the rotational motion. For simulating spectra in the slow-motion regime, the spin dynamics and the rotational dynamics have to be treated on an equal footing. For nitroxide radicals at X-band fields (∼0.34 T), Δ*ω*/2*π* ≈ 200 MHz so that the slow-motion regime is around *τ*_{c} ≈ 1–100 ns. This is the range often observed for nitroxides attached to proteins.

Several approaches for the simulation of slow-motion CW EPR spectra have been developed. They are based on motional models that range from full deterministic atomistic molecular dynamics (MD) to simple stochastic reorientation.^{1–5} A simple and, therefore, widely applied model is hindered Brownian diffusion (HBD).^{6} As illustrated in Fig. 1(b), the HBD model represents the tumbling spin center as a single rigid body undergoing Brownian rotational diffusion with an anisotropic rotational diffusion rate *R*. Besides the orientation ** Ω** of the body, no other spatial degrees of freedom are dynamic in this model. All internal degrees of freedom are considered fixed. The nano-environment (such as a protein or lipid environment), also assumed fixed, hinders the rotational diffusion. Its effect is modeled with an orientation-dependent potential energy function

*U*(

**) that imposes an orientation-dependent torque onto the rigid body.**

*Ω*The motional model is combined with spin quantum dynamics to simulate the slow-motion CW EPR spectrum. For the HBD model, the most efficient and widely used approach is a sophisticated frequency-domain method based on the stochastic Liouville equation (SLE), pioneered by Kubo and co-workers.^{7–18} It employs highly efficient numerical methods and has the advantage of significantly lower computational cost compared to other HBD solvers that are based on stochastic diffusive or jump trajectories.

The SLE theory was originally developed to simulate CW EPR spectra for slow-tumbling nitroxide radicals in solution or liquid crystals.^{10,13} The expressions derived were particular to nitroxides, and the orientational potentials employed were of high symmetry. Since then, the theory has been expanded to include two magnetic nuclei^{19,20} and it has been extended and applied to rigid bisnitroxides.^{21–23} The slowly relaxing local structure (SRLS) model was developed to include additional rotational dynamics of the cage encompassing a spin label.^{18,24,25} A nuclear magnetic resonance (NMR)-focused program, Spinach, can solve the SLE for large spin systems in the absence of orientational potentials.^{26} Simulation of slow-motion spectra for high-spin systems (*S* > 1/2) has also been reported. For instance, slow-motion spectra of triplets (*S* = 1) in the absence of magnetic nuclei have been simulated using the SLE^{27} or a discrete-jump approach.^{28–30} The SLE has been employed to calculate CW EPR linewidths for Gd^{3+} centers (*S* = 7/2) undergoing unhindered rotational dynamics by including additional internal degrees of freedom such as vibration^{31} or pseudo-rotation.^{32} Lower-symmetry orientational potentials have been utilized in work on biaxial liquid crystal phases.^{33} Software codes stemming from several of these works are available.

Despite these advances, there exist no comprehensive methodology and software for calculating slow-motion spectra using the SLE approach without constraints on the constitution of the spin center and/or on the complexity of the environment hindering the reorientational motion (represented by the orientational potential). Such an extended SLE theory is needed due to the increased use of spin labels to study protein dynamics and the increased use of spin labels other than nitroxides in recent years.

In this paper, we present expressions that extend SLE theory to (a) orientational potentials of arbitrary complexity and to (b) spin centers of arbitrary composition, i.e., with any number of electron and nuclear spins and spin Hamiltonian interaction terms (Zeeman, hyperfine, zero field, exchange interactions, nuclear quadrupole, etc.). Furthermore, we examine the issue of calculating the field-swept CW EPR spectrum. Traditionally, SLE solvers calculate a frequency spectrum that is then converted to a field-swept spectrum using a first-order approximation. This method is inaccurate for systems with highly anisotropic *g*-tensors, and we present a more accurate way to simulate field-swept spectra.

In the following, Sec. II comprehensively lays out the SLE theory. Section III briefly summarizes implementation details. Section IV illustrates the extended scope of the theory with simulations of high-spin systems, low-symmetry potentials, and multinuclear systems. A few numerical aspects are discussed. Section V contains concluding discussions. All methods presented in this paper are implemented in the open-source software package EasySpin.^{5,34,35}

## II. THEORY

In this section, we present the key expressions of SLE theory,^{13,16,17,37} including our extensions. We will use a series of space-fixed and body-fixed frames, which are shown in Fig. 2(a). The laboratory frame (L) is a space-fixed frame with its *z* axis along the external static magnetic field and its *y* axis along the oscillatory magnetic-field component of the microwave radiation. The potential frame (U) is a frame that is attached to the nano-environment hindering the motion of the spin center, e.g., a protein, a membrane bilayer, a polymer matrix, or a liquid-crystal phase. We limit ourselves to situations where proteins and membranes are immobile, i.e., where the potential frame is space-fixed. Often, the frame U is also referred to as the director frame. The potential frame will be our space-fixed reference frame. The mobile spin center is associated with a series of body-fixed frames that have time-dependent orientation: There is a molecular frame (M) associated with the molecular geometry. Each of the interaction tensors (*g*-tensors, hyperfine tensors, etc.) has a body-fixed eigenframe. The diffusion tensor (*vide infra*) is body-fixed as well, and the associated eigenframe is called the diffusion frame (D). It will serve as our body-fixed reference frame. To specify the relative orientation between any two frames, we use a triplet of Euler angles, indicated by ** Ω** = (

*α*,

*β*,

*γ*) and defined in Fig. 2(b). Some of the relative orientations used in the following are indicated in Fig. 2(a).

### A. Spin and rotational dynamics

The spin state of a tumbling spin center is described by the quantum spin density operator *ρ*(*t*). Its orientational state is described classically by the orientation ** Ω**(

*t*).

**indicates the orientation of a body-fixed reference frame relative to a space-fixed reference frame, as described above.**

*Ω*The evolution in time of *ρ*(*t*) is described by the Liouville–von Neumann equation

Here, *∂*_{t} indicates the time derivative, and [*H*, *ρ*] is the commutator *Hρ* − *ρH*. The spin dynamics is coupled to the rotational dynamics via the orientational trajectory of the spin centers, ** Ω**(

*t*).

*H*is the spin Hamiltonian operator (in angular-frequency units) summing all EPR-relevant interactions within the spin center and between the spin center and the external static magnetic field.

*H*is implicitly time-dependent through the time dependence of

**. During irradiation with the microwave field,**

*Ω**H*is explicitly time-dependent. However, for our purpose of simulating CW EPR spectra, we do not need to consider the interaction of the spin system with the microwave irradiation. It is possible to rewrite Eq. (1) as

where *H*^{×} is called the Hamiltonian commutation superoperator and is defined by its operation on *ρ*, *H*^{×}*ρ* ≡ [*H*, *ρ*]. The space spanned by all possible *ρ* is called the Liouville space. *H*^{×} is an operator on this space.

Instead of dealing with explicit orientational trajectories ** Ω**(

*t*) and the associated dynamic equation, the rotational rigid-body dynamics is modeled using the orientational distribution of the spin centers in the sample at time

*t*,

*P*(

**,**

*Ω**t*). The time evolution of this distribution is described by the Fokker–Planck-type differential equation

^{6}

where *Γ*(** Ω**) is the diffusion operator representing the rotational dynamics.

*Γ*is assumed to be independent of the spin degrees of freedom, reflecting the reasonable and accurate assumption that the rotational dynamics of a spin center is independent of its spin state.

The explicit form of *Γ* depends on the model that is used to describe the rotational dynamics.^{13} Here, we focus on hindered Brownian rotational diffusion (HBD) in the presence of an orienting potential *U*(** Ω**). Brownian motion assumes the absence of inertial motion, i.e., it assumes that the diffusion process is Markovian (memoryless).

*U*(

**) describes an orientation-dependent effective potential energy for the spin center which is a result of the interaction of the spin center with its immediate nanoenvironment. This potential encodes that different orientations of the spin center have different energies. It leads to a systematic torque on the spin center (in the downhill direction on the potential-energy surface). The associated**

*Ω**Γ*is

^{6}

Here, *J*_{x}, *J*_{y}, and *J*_{z} are differential angular-momentum operators around the axes *x*, *y*, and *z* of a body-fixed reference frame. *R*_{ij} are the real-valued elements of the body-fixed diffusion tensor, which is symmetric (*R*_{ij} = *R*_{ji}) and is assumed to be time-independent.

During a CW EPR experiment, the orientational distribution is at thermal equilibrium at all times. The thermal-equilibrium orientational distribution *P*_{eq} is stationary (*∂*_{t}*P*_{eq} = 0) and is given by

with the partition function

where *k*_{B} is the Boltzmann constant and *T* is the temperature.

The dynamical equations for the spin and the orientational degrees of freedom [Eqs. (2) and (3)] can be combined into a single equation, the stochastic Liouville equation (SLE)^{8,10,16}

Here, *ρ*(** Ω**,

*t*) is the total spin density operator for all spin centers with orientation

**at time**

*Ω**t*, no matter which orientation

*Ω*_{0}they had initially,

*P*(** Ω**′, 0) is the initial orientational distribution at time zero,

*P*(

**,**

*Ω**t*|

**′, 0) is the distribution at time**

*Ω**t*given that the orientation was

**′ initially, and**

*Ω**ρ*(

*t*|

**′) is the spin density matrix at time**

*Ω**t*given that the orientation was

**′ initially. Note that**

*Ω**ρ*(

**,**

*Ω**t*) is different from

*ρ*(

*t*) in Eq. (1).

In the presence of a potential, *Γ* of Eq. (4) is not Hermitian. It is advantageous to transform it to Hermitian form.^{6} This can be achieved by the transformation $\Gamma \u0303=Peq\u22121/2\Gamma Peq1/2$. The diffusion equation then becomes $\u2202tP\u0303=\u2212\Gamma \u0303P\u0303$, with the scaled distribution function $P\u0303=Peq\u22121/2P$. Defining in addition the scaled density $\rho \u0303=Peq\u22121/2\rho $, the SLE with $\Gamma \u0303$ reads

where now both *H*^{×} and $\Gamma \u0303$ are Hermitian. Since *H*^{×} and $\Gamma \u0303$ are time-independent, the integral of this equation is

(We omit from now on the explicit indication of the dependence on ** Ω**.) An explicit form of the Hermitianized diffusion superoperator is

^{6}

In the body-fixed diffusion frame, the diffusion tensor is diagonal, i.e., *R*_{ij} = *R*_{i}*δ*_{ij}. We use this as the body-fixed reference frame without loss of generality but with significant simplifications in the expressions for $\Gamma \u0303$. In the diffusion frame, $\Gamma \u0303$ from Eq. (11) simplifies to^{37}

where we have separated the potential-independent and potential-dependent terms.

The potential function *U* is best expanded in a complete orthogonal basis of Wigner functions $DM,KL$ because $DM,KL$ have simple transformation properties under rotation—they are simultaneous eigenfunctions of *J*^{2} and *J*_{z}. The expansion is

where the − sign is a matter of convention. Here, *Ω*_{U→D} is the triplet of Euler angles (*α*, *β*, *γ*) that describes the orientation of the body-fixed diffusion frame (D) relative to the space-fixed potential frame (U) (see Fig. 2). The Wigner functions $DM,KL$ are defined as^{38,39}

where $dM,KL$ are real-valued functions consisting of sums of products of cos(*β*/2) and sin(*β*/2). The sum in Eq. (13) runs over all possible combinations of integer rank *L* and projections *M* and *K* (*L* ≥ 0; −*L* ≤ *M* ≤ *L*; −*L* ≤ *K* ≤ *L*). $\lambda M,KL$ are dimensionless coefficients that may be complex-valued. Since *U* is real-valued and $d\u2212M,\u2212KL=(\u22121)M\u2212KdM,KL$, the coefficients satisfy $\lambda \u2212M,\u2212KL=(\u22121)M\u2212K\lambda M,KL*$.

Applying the *J*_{i} operators to the Wigner expansion of *U*, the potential-dependent part of $\Gamma \u0303$ in Eq. (12) reduces to a linear combination of Wigner functions. $\Gamma \u0303$ becomes

The scalar expansion coefficients $X\u0303M,KL$ depend on *R*_{i} and on $\lambda M,KL$ and are given in Appendix C. They have the symmetry $X\u0303\u2212M,\u2212KL=(\u22121)M\u2212K(X\u0303M,KL)*$. If all $\lambda M,KL$ are real-valued, then all $X\u0303M,KL$ are real-valued as well. The largest *L* for which $X\u0303M,KL$ is non-zero is twice the largest *L* for which $\lambda M,KL$ is non-zero. Note that the potential-dependent part of $\Gamma \u0303$ is a purely multiplicative operator, whereas the potential-independent part contains the differential operators *J*_{i}.

The use of the complete expansion [Eq. (13)] and the ensuing expressions for $X\u0303M,KL$ constitute extensions of the existing theory, where the expansion is limited to terms with *M* = 0, with low even values of *L* and *K* (*L* = 2, 4; *K* = 0, 2, 4) and with real-valued coefficients $\lambda M,KL$. These specializations stem from the fact that the theory was initially developed for spin probes in uniaxial liquid crystals and model membranes with nanoenvironments of cylindrical symmetry (D_{∞h}) and required only high-symmetry potentials.^{40} However, the nanoenvironment hindering the rotational motion of spin labels on proteins is generally of much lower symmetry (C_{1}), and the orientational potential expansions, consequently, needs more terms. Therefore, the general expansion of Eq. (13) is necessary for being able to broadly apply the HBD model to protein-attached spin labels.

### B. The spin Hamiltonian

A general EPR spin Hamiltonian for a system of coupled electrons and magnetic nuclei is

The first two terms are the electron and nuclear Zeeman interactions, with the Bohr magneton *μ*_{B}, the external magnetic field vector ** B**, the

*i*th electron

*g*-tensor

*g*

^{(i)}, the

*i*th electron spin vector

*S*^{(i)}, the nuclear magneton

*μ*

_{N}, the

*k*th nuclear g factor $gn(k)$, and the

*k*th nuclear spin vector

*I*^{(k)}. The third term contains all zero-field interactions for electron spins

*S*

^{(i)}≥ 1, with the zero-field interaction tensors

*D*

^{(i)}. The fourth term collects all hyperfine interactions between electrons and nuclei, with the hyperfine tensors

*A*

^{(i,k)}. The fifth term sums over all electron–electron interactions, with the associated coupling tensors

*J*

^{(i,j)}. The last term represents the electric quadrupole interactions for nuclei with

*I*

^{(k)}≥ 1, with the quadrupole tensors

*P*

^{(k)}. We assume that all vectors and tensors are represented in the same space-fixed frame. For a tumbling molecule,

**is time-independent in this frame, and all the interaction tensors (**

*B**g*,

*D*,

*A*,

*J*,

*P*) are time-dependent.

The spin Hamiltonian of Eq. (16) is a sum over bilinear forms,

where *a*_{μ} and *b*_{μ} are the spin vector operators or the magnetic field vector, *F*_{μ} are the time-dependent Cartesian tensors, and the axes *x*, *y*, and *z* refer to a space-fixed frame. The index *μ* runs over all interaction terms. Table I lists the correspondence between Eqs. (17) and (16).

Interaction type . | a_{μ}
. | F_{μ}
. | b_{μ}
. |
---|---|---|---|

Electron Zeeman | B | +μ_{B}g^{(i)} | S^{(i)} |

Nuclear Zeeman | B | $\u2212\mu Ngn(k)$ | I^{(k)} |

Zero-field splitting | S^{(i)} | D^{(i)} | S^{(i)} |

Hyperfine | S^{(i)} | A^{(i,k)} | I^{(k)} |

Electron–electron coupling | S^{(i)} | J^{(i,j)} | S^{(j)} |

Nuclear quadrupole | I^{(k)} | P^{(k)} | I^{(k)} |

Interaction type . | a_{μ}
. | F_{μ}
. | b_{μ}
. |
---|---|---|---|

Electron Zeeman | B | +μ_{B}g^{(i)} | S^{(i)} |

Nuclear Zeeman | B | $\u2212\mu Ngn(k)$ | I^{(k)} |

Zero-field splitting | S^{(i)} | D^{(i)} | S^{(i)} |

Hyperfine | S^{(i)} | A^{(i,k)} | I^{(k)} |

Electron–electron coupling | S^{(i)} | J^{(i,j)} | S^{(j)} |

Nuclear quadrupole | I^{(k)} | P^{(k)} | I^{(k)} |

For modeling rotational dynamics, it is advantageous to reformulate *H* from a sum over terms with distinct physical origin to a sum over terms with distinct rotational properties. For this, we rewrite *H* as a sum over scalar products of irreducible spherical tensors and tensor operators,^{13,18,41}

Here, $F\mu (l)$ are spherical tensors constructed from the matrix elements of *F*_{μ}, and *T*^{(l)}(*a*_{μ}, *b*_{μ}) are spherical tensor operators constructed from the Cartesian components of *a*_{μ} and *b*_{μ}. *l* is the rank of the spherical tensor. Each Cartesian tensor is decomposed into three spherical tensors with ranks *l* = 0, 1, and 2. These three tensors have distinct rotational properties. Each spherical tensor has 2*l* + 1 components (*m* = −*l*, …, *l*), indicated by *F*^{(l,m)} and *T*^{(l,m)}. The scalars $F\mu (l,m)$ and the operators *T*^{(l,m)}(*a*_{μ}, *b*_{μ}) can be constructed in a straightforward fashion from *F*_{μ}, *a*_{μ}, and *b*_{μ}. All the required expressions are listed in Appendix A.

As in Eqs. (16) and (17), also in Eq. (18), all vectors and interaction matrices are represented in a space-fixed frame. Choosing the laboratory frame (L), indicating it by an additional subscript, and using (−1)^{m}*F*^{(l,−m)} = (−1)^{l}*F*^{(l,m)*} (see Appendix A), we write the Hamiltonian as

The laboratory frame components $F\mu ,L(l,m)$ of $F\mu (l)$ are expressed in terms of the time-independent diffusion frame components $F\mu ,D(l,m)$ using

where $Dm,m\u2033l$ are again Wigner functions as defined in Eq. (14) and *Ω*_{L→D} represents the Euler angles that parameterize the transformation from the laboratory frame L to the body-fixed diffusion frame D (see Fig. 2). We decompose this transformation into two consecutive transformations via the intermediate potential frame (U),

Combining this with the previous two equations and defining the operators

and

gives the compact expression

Since both L and U are space-fixed stationary frames, the only time dependence is in the stochastically varying orientation *Ω*_{U→D}. The rotational time dependence is now fully isolated in the Wigner function prefactors. The time-independent operators $Qm\u2032,m\u2033l$, which we call rotational basis operators (RBOs), contain all the internal specifics of the spin system as well as the orientation of the potential frame relative to the lab frame. Note that this expression is similar to the potential-dependent term of the diffusion operator in Eq. (15).

Equation (24) is general in the sense that the RBOs can be constructed using the same procedure regardless of the number of spins in the spin system, the number of interaction terms, or the nature or relative size of those terms. No matter how large or complex the tumbling spin system is, all of the information about the system that is needed to calculate the EPR spectral response is collected into 35 RBOs (one rank-0, 9 rank-1, and 25 rank-2). If all the interaction matrices are symmetric, as is commonly the case, all rank-1 RBOs vanish. If the symmetry is high (for example, axial and collinear interaction tensors), then the number of non-zero RBOs reduces further. Equation (24) covers two separate situations: (i) a single potential-frame orientation $\Omega U\u2192D$ (e.g., an oriented membrane or a protein crystal) and (ii) an orientational distribution of potential-frame orientations (such as a solution of essentially immobile proteins or liposomes). The latter model with a disordered static distribution of proteins, but mobile spin labels, has been termed the MOMD (microscopic order macroscopic disorder) model.^{14} For calculation in such cases, the $Pm,m\u2033l$ operators are precomputed, and the $Qm\u2032,m\u2033l$ operators are efficiently computed from $Pm,m\u2033l$ for each potential-frame (protein) orientation, $\Omega U\u2192D$, without the need of recomputing any matrix elements.

The expression for the Hamiltonian presented above is quite general: it includes situations with tensor eigenframes tilted arbitrarily with respect to the diffusion frame; it includes rank-1 terms; and it permits any number of spins and types/strengths of interactions. In addition, the pre-calculation of RBOs renders simulations efficient, particularly for disordered samples.

### C. The CW EPR spectrum

The SLE, together with the expressions for the diffusion operator and the spin Hamiltonian, can be used to calculate the signal from any type of EPR experiment (pulse or CW). Here, we focus on CW EPR. The frequency-swept CW EPR spectrum is proportional to the Fourier–Laplace transform of the free-induction decay (FID) following a non-selective microwave pulse. (Alternatively, it can be calculated using linear-response theory.^{16,42}) We start with the equilibrium density distribution *ρ*_{eq} ≈ *S*_{z}*P*_{eq} (corresponding to $\rho \u0303eq\u2248SzPeq1/2$), where *z* is the laboratory *z* axis. The detectable part of the state immediately after a microwave pulse with oscillating magnetic-field component along the laboratory *y* axis is *S*_{x}*P*_{eq}, where *x* is the laboratory *x* axis. The FID signal is the expectation value of *S*_{x} (or of $\u2211iSx(i)$ if there are multiple electron spins) for this state,

where *H*^{×} is the static spin Hamiltonian superoperator in the absence of the microwave field. The notation ⟨*u*|$v$⟩ indicates integration of *u*^{†}$v$ over all orientations and trace over all spin space.

The real part of the Fourier–Laplace transform of the FID gives the spectrum

Here, *ω* is the microwave angular frequency and *B* is the magnitude of the applied static magnetic field.

### D. Basis expansion

In order to evaluate the expression in Eq. (26), the quantities involved (*H*^{×}, $\Gamma \u0303$, and $SxPeq1/2$) are expanded in an appropriate basis that encompasses both the quantum spin states and the classical orientational degrees of freedom. This transforms the integral expression in Eq. (26) into a linear-algebra expression (vector × matrix inverse × vector) that can be solved efficiently using numerical methods.

We use a basis where each basis function $\sigma \xi $ is a direct product of a spin basis function $\sigma $ and an orientational basis function $\xi $,

The spin basis functions, $\sigma $, are given by the complete set of single-transition operators in the Zeeman (high-field) basis.^{43} Each single-transition operator is parameterized by a pair of projection numbers for each spin in the system

where the spin projection quantum numbers (along the *z* axis of the potential frame) $mi\u2032$ and $mi\u2032\u2032$ for spin *i* run from −*S*_{i} to +*S*_{i} or −*I*_{i} to *I*_{i} for nuclear spins. The basis functions in Eq. (28) are orthonormal, and there are a total of $\u220fi(2Si+1)2\u220fk(2Ik+1)2$ of them. Instead of using $mi\u2032$ and $mi\u2032\u2032$, the single-transition operators can also be indexed by *p*_{i} = $mi\u2032$ − $mi\u2032\u2032$ and *q*_{i} = $mi\u2032$ + $mi\u2032\u2032$. This facilitates the application of spin-space truncation schemes, such as limiting the basis to single-quantum transitions (|*p*_{i}| ≤ 1) and application of the high-field approximation (only *p*_{i} = +1).^{13}

The orientational basis functions are normalized Wigner functions [see Eq. (14)] of *Ω*_{U→D},

with the full range of *L*, *M*, and *K*. These functions are orthonormal, i.e.,

They are a convenient choice because they have simple properties under rotation and are also the coefficients for the transformation of the spherical tensors and tensor operators between coordinate frames. The complete set of Wigner functions is infinite, and in practice, the basis must be truncated. A simple truncation scheme is to impose separate upper limits on even *L*, odd *L*, |*K*|, and |*M*|. A more sophisticated pruning procedure has been developed.^{15,17} In general, slower rotational diffusion and larger interaction tensor anisotropies require larger orientational basis sets to produce accurate and converged simulated spectra. A crucial practical point is the confirmation that a simulated spectrum is sufficiently converged as a function of orientational basis size.

with Δ*M* = *M*_{1} − *M*_{2} and Δ*K* = *K*_{1} − *K*_{2}. The sum runs over all *l* = 0, 1, 2 that also satisfy |*L*_{1} − *L*_{2}| ≤ *l* ≤ *L*_{1} + *L*_{2} and *l* ≥ |Δ*M*| and *l* ≥ |Δ*K*|. Since *l* ≤ 2, the matrix is at most pentadiagonal in each of *L*, *M*, and *K*. Any zero RBO $Q\xd7l\Delta M,\Delta K$ will further thin this nonzero bandedness along *M* and *K*.

The expressions for the Wigner function matrix elements needed in Eq. (31) are given in Appendix B. The matrix representations of $Q\xd7l\Delta M,\Delta K$ are constructed as follows: (1) The spin matrices $Sx(i)$, $Sy(i)$, $Sz(i)$, $Ix(k)$, $Iy(k)$, and $Iz(k)$ (where *x*, *y*, and *z* refer to the laboratory frame) in the standard Hilbert-space $mi$ representation are constructed for each spin, and spherical tensor matrices are constructed using the expressions from Table III; (2) the $Qm\u2032,m\u2033l$ matrices are constructed using Eqs. (22) and (23); (3) the corresponding Liouville-space matrices are constructed using the Kronecker product ⊗ according to

in which I is the identity matrix of the same size as the *Q* matrix; and (4) any spin-space truncations are applied.

The matrix elements of $\Gamma \u0303$ in the chosen basis are

with the required matrix elements of *Γ*_{0} given in Appendix C. $\Gamma \u0303$ is diagonal in the spin quantum numbers. The sum runs over all *L* that satisfy |*L*_{1} − *L*_{2}| ≤ *L* ≤ *L*_{1} + *L*_{2} and *L* ≥ |Δ*M*| and *L* ≥ |Δ*K*|.

In general, the $\Gamma \u0303$ matrix is Hermitian. In our choice of frame where *R* is diagonal (the diffusion frame) and in the absence of a potential, *Γ* is diagonal if the diffusion tensor is isotropic or axial; otherwise, it is tridiagonal. For a non-zero potential, if the potential contains only terms with *M* = *K* = 0, or if all *λ* coefficients are real-valued, then $\Gamma \u0303$ is real-valued and symmetric; otherwise, it is complex-valued and Hermitian.

Finally, the vector elements of $SxPeq1/2$ in the chosen basis are

$\sigma |Sx$ is obtained by vectorizing the $mi$ matrix representation of *S*_{x} in column-major order, i.e., by stacking the columns of the matrix. The vector elements of $Peq1/2$ are given by

Analytical evaluation of the triple integral in this expression is only possible for very simple potentials.^{17} In general, the integral needs to be evaluated numerically. If *U* does not depend on *α* and/or *γ* (representing the cylindrical symmetry of the environment and/or the spin label), then the dimensionality of the integral reduces.

The norm-squared of the $Peq1/2$ vector equals 1 in the infinite-basis limit.^{38} If the basis is excessively truncated, this will be significantly less than 1. This can serve as a diagnostic of basis overtruncation.

### E. Numerical evaluation

With matrix representations of all relevant quantities in hand, we now describe how to evaluate Eq. (26). Denoting with ** A** the matrix representation of $iH\xd7+\Gamma \u0303$ in the chosen truncated basis, with

**the vector representation of $SxPeq1/2$ in the same basis, and with**

*b***I**the identity matrix of the same size as

**, the spectral line shape function of Eq. (26) is**

*A*Evaluation of this as a function of *ω* for a fixed external magnetic field (i.e., fixed ** A**) gives the frequency-swept EPR spectrum.

A variety of approaches are available to numerically evaluate the right-hand side of Eq. (36). The straightforward diagonalization method^{16,44} computes eigenvectors and eigenvalues of ** A** (such that

**=**

*A*

*UΛU*^{−1}with the diagonal matrix of eigenvalues

**and the matrix of eigenvectors**

*Λ***) and transforms the expression into the eigenbasis of**

*U***, resulting in a sum over Lorentzian lines**

*A*which can be easily evaluated as a function of *ω*. In the above equation, both *S*_{k} and *Λ*_{kk} are in general complex-valued. The real and imaginary parts of *Λ*_{kk} give the linewidth and the line position, respectively. *S*_{k} determines the amplitude and phase of the line. Another method uses a linear-equation solver to calculate the solution ** x**(

*ω*) of (

**− i**

*A**ω*

**I**)

**(**

*x**ω*) =

**for each**

*b**ω*and then obtains the spectral values via

*I*(

*ω*) = Re(

*b*^{†}

**(**

*x**ω*)). This method does not explicitly break the signal down into its line shape components. Both methods are general and work for any

**without symmetry restrictions. A variety of well-established numerical algorithms (conjugate gradient, Lanczos, etc.) can be used for diagonalization and as linear solvers.**

*A*A much faster and very elegant method is applicable if ** A** is complex symmetric. In this case, the complex symmetric Lanczos algorithm

^{45}is used to tridiagonalize

**with**

*A***as the starting vector.**

*b*^{11,12,16,17}The spectral function is then represented as a continued-fraction expansion containing the matrix elements of the tridiagonal matrix, and it can be evaluated very efficiently left-to-right using the modified Lentz method.

^{46}This method is fast because it converges rapidly as a function of iteration count and because the spectral function can be cheaply evaluated for the entire desired range of frequencies with only one tridiagonalization. The continued-fraction expansion is possible only if

**is complex symmetric, i.e., if the matrix representations of both**

*A**H*

^{×}and $\Gamma \u0303$ are real-valued. In the

*σξ*basis, this is only the case for a very small class of high-symmetry situations. However, under certain conditions, it is possible to convert

**to complex symmetric form by transforming it to a new**

*A**K*-symmetrized basis with functions $\sigma LMK\xafjK$ that are parameterized by

*L*,

*M*, $K\xaf$, and

*j*

^{K}, where $K\xaf=|K|$, and

*j*

^{K}= ±1 for

*K*≠ 0 and

*j*

^{K}= (−1)

^{L}for

*K*= 0.

^{13,17}The transformation is achieved using $A\u2032=TKATK\u2020$, with the elements of the unitary transformation matrix

*T*_{K}given by

In a basis ordering where the functions with identical *L*, *M*, and $K\xaf$ are adjacent, the matrix *T*_{K} is block diagonal, consisting of a sequence of 2 × 2 blocks (for $K\xaf\u22600$) and 1 × 1 blocks (for $K\xaf=0$) along the diagonal.

Prior work derived explicit expressions for the matrix elements of *H*^{×} and $\Gamma \u0303$ in this new basis for some special cases.^{13,17,37} Here, we perform this transformation numerically at the matrix level. Since all matrices are sparse (the number of non-zero elements in *T*_{K} only grows linearly with basis size), this operation is very efficient.

In the *K*-symmetrized basis, the matrix of any Hamiltonian of the form given in Eq. (17) is real-valued, irrespective of the number, symmetry, and relative orientation of interaction tensors—as long as the tilt between the potential frame and the lab frame only involves a *β* angle (*α* = *γ* = 0). The situation is more complicated for $\Gamma \u0303$. In the *LMK* basis, $\Gamma \u0303$ is real-valued for arbitrary complicated potentials. However, in the *K*-symmetrized basis, it is real-valued only if all potential coefficients are real-valued and the potential contains only terms with either all *M* = 0 or all *K* = 0 (or both). Otherwise, the transformation will generate a complex-valued Hermitian $\Gamma \u0303$. In these cases, the transformation is not useful, and the spectral line shape function must be calculated using the diagonalization or linear-solver methods instead of the Lanczos method.

### F. Field sweep

The line shape of Eq. (36) is expressed as a function of microwave frequency *ω*, in the presence of a constant applied magnetic field. In practice, CW EPR experiments are performed by varying the external field strength *B* over a range from *B*_{min} to *B*_{max} while irradiating the sample at fixed frequency *ω*_{mw}.

For the field-swept spectrum, a separate Hamiltonian at each field point *B* is required, as evident from Eq. (26). Since the Hamiltonian is linear in the magnetic field, Hamiltonians for different field values can be constructed very efficiently. For this, each RBO $Qm\u2032,\u2032m\u2033l$ is separated into a field-dependent term and a field-independent term,

The field-independent RBOs $Qm\u2032,m\u20330,l$ and $Qm\u2032,m\u20331,l$ are pre-calculated, and then Eq. (39) is used to assemble $Qm\u2032,m\u2033l$ for each field point, thus providing minimal overhead in Hamiltonian re-calculation. To minimize the number of field points for which the spectral function needs to be evaluated, an adaptive iterative bisection interpolation method analogous to one used for constructing Zeeman energy diagrams can be used.^{47}

If the complex symmetric Lanczos method is applicable (i.e., if the matrix representation of $iH\xd7+\Gamma \u0303$ is complex symmetric), then the calculation of a frequency-swept spectrum is orders of magnitude more efficient than the calculation of a field-swept spectrum. In this case, the field-swept spectrum can be approximated by a frequency-swept spectrum,

Here, *B*_{0} = (*B*_{min} + *B*_{max})/2 is the center of the desired field range, and the auxiliary frequency *ω* is obtained from *B* in one of two ways,

*g*_{iso} is the isotropic g-value of the spin system. Approach (i) is exact for systems with g anisotropy and no other anisotropy (hyperfine, zero-field). For other systems, it leads to systematic errors in the resonance field positions. Approach (ii) is accurate if g is isotropic and the Zeeman interaction dominates; otherwise, there are systematic errors in the resonance field positions. At X-band, the errors in approach (ii) are very small for nitroxides, but they can be significant for copper. Approach (i) and approach (ii) give identical frequencies for *B* = *B*_{0} and for *B* = *ℏω*_{mw}/*μ*_{B}*g*_{iso}.

Finally, the field-modulated first-harmonic spectrum is obtained from the field-swept absorption spectrum via pseudo-field modulation^{48} or, in the limit of zero modulation amplitude, via numerical differentiation.

## III. METHODS

The theory outlined above is implemented in version 6 of the open-source MATLAB-based EPR simulation package EasySpin^{5,34,35} in the function chili.

All matrices are sparse and are handled in MATLAB’s compressed sparse column (CSC) format. The spectral line shape can be calculated using either the diagonalization, the linear-solver, the Lanczos method (only if applicable, given the symmetry of $\Gamma \u0303$), or the biconjugate gradient stabilized method as implemented in the MATLAB function bicgstab.

The integrals necessary for the vector elements in Eq. (34) are evaluated numerically using the 15-point Gauss–Kronrod quadrature as implemented in the MATLAB function integral, with absolute and relative error tolerances of 10^{−6}. Even after using selection rules to reduce the number and dimensionality of necessary numerical integrations, the evaluation of the root-equilibrium-distribution vector [Eq. (34)] constitutes the major performance bottleneck in the overall calculation in the presence of low-symmetry potentials.

Wigner 3-j symbols for Eqs. (B4) and (C2) are evaluated using specialized expressions for small angular momenta [min(*L*_{1}, *L*, *L*_{2}) ≤ 2]^{38,39,49} and using a general expression otherwise.^{50} For large angular momenta [max(*L*_{1}, *L*, *L*_{2}) > 20], arbitrary-precision integer arithmetic (using the Java class BigInteger) is used for intermediate results. This results in a significant computational bottleneck for large angular momenta. Values of 3-j symbols calculated via these methods in EasySpin have been extensively tested for angular-momentum values up to 5000 against arbitrary-precision results obtained using Wolfram Mathematica 11.

## IV. APPLICATION EXAMPLES

In this section, we illustrate the scope of the extended theory by demonstrating examples of simulated CW EPR spectra for a high-spin system, for lower-symmetry orienting potentials, and for a spin system with multiple nuclei. We also discuss some numerical aspects.

### A. High-spin system: Gadolinium(III)

As an example of a high-spin system, Fig. 3 shows a series of simulated spectra of a prototypical Gd(III) complex over a range of rotational time scales from the fast-motion regime to the quasi-rigid limit. Gd(III) is a high-spin 4f^{7} ion with an *S* = 7/2 ground state. Its magnetic properties are described by the spin Hamiltonian *ℏH* = *μ*_{B}*g*** BS** +

*S**D*

**with isotropic g-value and the diagonal zero-field tensor**

*S**D*= diag(−

*D*+

*E*, −

*D*−

*E*, 2

*D*).

The simulations in Fig. 3 are plotted as absorption spectra. With a rotational correlation time *τ*_{c} = 0.1 ns (*R* = 1/6*τ*_{c} = 1.7 rad^{2} ns^{−1}), the spectrum consists of a single peak because rotational motion is so rapid that the zero-field interaction is almost averaged out. At *τ*_{c} = 100 ns (*R* = 1.7 rad^{2} *μ*s^{−1}), the motional EPR spectrum is indistinguishable from the rigid-limit simulation.

### B. Orientational potentials

Figure 4 illustrates the effect of different potential energy functions on the X-band spectrum of a nitroxide radical in the quasi-rigid limit with *τ*_{c} = 10 ns. Each simulation has a different orientational potential defined by the pair of coefficients $\lambda M,KL=(\u22121)M\u2212K\lambda M,KL$ = 1, with various combinations of *L*, *M*, and *K*. The gray spectrum is calculated in the absence of an orientational potential. The spectral line shape depends strongly on the particular term. The size of the deviation from the unhindered potential-free spectrum is dependent on (*L*, *M*, *K*) as well. For a given *λ*, this deviation is generally smaller for odd values of *M* and *K* (and *L*) than for even values.

### C. Multiple nuclei

As an example of slow-motion simulations for spin systems with multiple magnetic nuclei, Fig. 5 shows simulations of a model copper(II) complex in which two of the ligand atoms are nitrogen-14 nuclei and the other ligand atoms are non-magnetic. The copper-63 nucleus has nuclear spin *I* = 3/2, and each of the nitrogen-14 nuclei has nuclear spin *I* = 1. An orientational potential with the coefficients $\lambda 2,22=\lambda \u22122,\u221222=1$ was used.

The simulations show how the spectrum changes with rotational correlation times from the fast-motion regime (small *τ*_{c}) to the quasi-rigid limit (large *τ*_{c}). At *τ*_{c} = 0.1 ns, the spectrum is close to the fast-motion regime. Most of the structure in the spectrum is due to averaged nitrogen hyperfine interactions, with some weak residual features in the low-field region from copper hyperfine splittings. As *τ*_{c} increases, the copper peaks become more distinct and reveal the non-isotropic orientational potential.

A major challenge with the simulation of multi-spin spectra using the SLE is the rapid factorial growth of the spin basis size with the number of spins. Three approaches at truncating this space can be applied individually or simultaneously: (a) limit the coherence orders |*p*_{k}| for each nucleus *k*; (b) limit the total maximum nuclear coherence order |*∑*_{k}*p*_{k}|; and (c) utilize the high-field approximation where for the electron spin *S* only the subspace with *p*_{S} = +1 is utilized. Truncation approaches based on subpartitioning based on the spin coupling topology, well established in NMR,^{26} are only marginally beneficial for EPR systems since the latter typically have a star topology (all nuclei coupled to a central electron spin) that cannot be partitioned into subgraphs.

A perturbation-based approach at reducing the spin-space dimension is applicable if the hyperfine couplings are such that a subset of nuclei have hyperfine couplings *A* small enough to fall into the fast-motion regime (*A* ≪ 1/*τ*_{c}). In this approach, the SLE-simulated spectrum of the spin system containing only the nuclei that fall into the slow-motion regime is convolved with the perturbation-based simulated spectrum of the remaining fast-motion-regime nuclei. This post-convolution approach factors the spin space and has been implemented for isotropic diffusion in the absence of a potential.^{51,52} It is not applicable if the hyperfine tensors in the system are of comparable magnitude.

The application of these truncations and approximations is crucial for making many-spin simulations feasible. Numerical experiments must be performed to carefully ascertain that any truncation does not degrade the convergence of the simulated spectrum. Practically, simulations of systems with more than about six spins are very memory- and time-consuming. Luckily, very few slow-motion CW EPR spectra resolve splittings or other features from more than a few electron or nuclear spins.

### D. Numerical aspects

Lanczos and conjugate-gradient methods are numerically unstable under certain conditions since they do not produce a fully orthogonal basis. In EPR simulations, it is known that this instability can manifest, for instance, in the presence of anisotropic diffusion tensors and very slow motion, leading to a failure to obtain a converged spectrum. In case numerical instabilities are encountered, a slight change in the orientational basis (by adding a few rotational basis functions) often mitigates the issue. As a safe fallback, other linear solvers or diagonalization can be used as described above. However, this comes at a significant computational cost. Our implementation shows reasonable stability of the Lanczos-based approach in difficult regimes. For example, Fig. 6 illustrates simulations in the presence of a strongly anisotropic diffusion tensor, with correlation times up to 1000 ns. With the chosen basis, which is not larger than necessary to obtain smooth spectra, convergence is not an issue.

Another practical numerical aspect is the time cost of simulations. It strongly depends on whether the continued-fraction expansion of the spectral function can be computed (via, e.g., the complex symmetric Lanczos method). If yes, only one tridiagonalization is necessary for simulating the entire spectrum. If not, diagonalization or linear-solver methods need to be used at every field/frequency point and result in significant runtime slowdowns. For the spectra in Fig. 3, the slowdown is >100× (0.92 s vs 123.5 s). For the simulations in Fig. 4, the slowdowns range between 30× and 7000×.

## V. DISCUSSION

In this paper, we have laid out expressions to apply the SLE methodology for HBD models to spin systems with any number or nature of interaction terms. From the familiarly defined lab-frame spin vector operators and molecular-frame Cartesian interaction tensors, ISTOs and RBOs are constructed. From the RBOs and the choice of orientational basis, the elements of the Hamiltonian superoperator and of the rotational diffusion operator are calculated. Once the calculation of the elements of *H*^{×} and *Γ* is complete, the spectral line shape function is evaluated.

We have also presented an expression for the diffusion operator in the presence of orienting potentials of any form, even of very low symmetry. This is useful for two scenarios: (1) for fitting an effective potential to an orientational histogram obtained from a molecular-dynamics trajectory and (2) for fitting a potential expansion to an experimental CW EPR spectrum. For the first case, many terms in the potential expansion are needed. For the second case, one has to consider the limited information content and the given signal-to-noise ratio of the CW EPR spectrum and limit modeling to potentials with relatively few terms.

Care must be taken to assure the orientational basis and the spin basis are large enough to provide sufficient convergence of the spectrum. The allowable truncation level depends on the spin Hamiltonian parameters, the rotational diffusion rates, and the orientational potential.

For systems with highly anisotropic *g*-tensors, care must be taken to simulate explicit field sweeps when fitting experimental spectra. Using an approximate field sweep method, it is possible that the fit parameters extracted from the best-fit simulation will be inaccurate due to the error in the peak positions and spectral width of the frequency-to-field converted simulation. Explicit field sweep simulations are more time consuming, but the construction of field-independent and field-dependent RBOs allows for reasonable simulation times in small to medium-sized spin systems.

Several important numerical challenges remain. (1) One is the efficient evaluation of the integrals for the basis representation of the root equilibrium distribution $Peq1/2$, as given in Eq. (34). Numerical integration for this, as done in this work, constitutes a significant computational bottleneck, particularly for large orientational basis sets where the basis contains highly oscillatory large-*L* basis functions. (2) Another challenge is the identification of additional symmetrizing basis transformations, similar to the *K*-symmetrization in Eq. (38), and the *M*-symmetrization/truncation, which is applicable only in the high-field regime^{13} for *M* = 0 potentials. The goal of the symmetrizations is to allow for efficient basis truncation and for the application of the complex symmetric Lanczos method. (3) The factorial growth of the spin space with the number of spins represents a challenge for large spin systems since the simulation times grow factorially, despite the spin-space truncations outlined above. (4) The Lanczos algorithm can be numerically unstable and fail to converge occasionally. Future advances addressing these challenges will improve the performance and robustness of the simulations.

It is important to reiterate the assumptions underlying the HBD model: (1) The spin dynamics does not affect the rotational dynamics. This assumption is accurate in all known cases. (2) The spin center is internally rigid, i.e., the principal values and relative orientations of the interaction tensors are time-independent in a body-fixed frame. This assumption is violated in some Gd(III) complexes, where there can be significant internal conformational dynamics.^{31,32} In such cases, the HBD model needs to be expanded by additional degrees of freedom, with associated potential energy functions and interaction tensor dependencies. (3) The rotational motion is Brownian, i.e., there are no inertial effects. This assumption is reasonably well satisfied for spin labels in solution and for spin labels at mobile solvent-exposed sites in proteins but is less appropriate for more restricted sites, where a multi-site Markov jump model might be more successful.^{4,5,53} (4) The motion can be described with a single diffusion tensor (i.e., there is a single time scale). Since the rotational motion of a spin label is a consequence of simultaneous torsional motions around several different bonds, rotational motion may occur over multiple time scales. To what degree an assumption about a single time scale is adequate remains to be explored. (5) The external potential is time-independent, at least on the time scale of the rotational diffusion. This assumption disregards the fact that different parts of the nano-environment might move at different time scales and that some of these time scales might be on the order of the spin label diffusion time scale, which would complicate the situation.

Our extensions of the SLE theory now allow assumptions (4) and (5) to be tested against atomistic molecular-dynamics simulations and against experimental data. The HBD model can then be quantitatively compared to other approaches, such as Markov state models.

## ACKNOWLEDGMENTS

This work was supported by the National Science Foundation (Grant No. CHE-1452967, S.S.), the National Institutes of Health (Grant No. GM125753, S.S.), and the Research Corporation for Science Advancement (Award No. 23447, S.S.). We thank David E. Budil (Northeastern University) and Keith A. Earle (University at Albany) for many helpful discussions as well as EasySpin users for bug reports.

### APPENDIX A: SPHERICAL TENSORS

Tables II and III list irreducible spherical tensor components.^{41} The tables give expressions for the conversion of vectors and tensors from their Cartesian form to irreducible spherical tensor form. Table III can be derived from Table II using the Cartesian tensor *C* with *C*_{ij} = *b*_{i}*a*_{j} (not *a*_{i}*b*_{j}).

(l, m)
. | F^{(l,m)}
. |
---|---|

(0, 0) | $\u221213(Fxx+Fyy+Fzz)$ |

(1, 0) | $\u2212i2(Fxy\u2212Fyx)$ |

(1, ±1) | $\u221212[(Fzx\u2212Fxz)\xb1i(Fzy\u2212Fyz)]$ |

(2, 0) | $+23[Fzz\u221212(Fxx+Fyy)]$ |

(2, ±1) | $\u221312[(Fxz+Fzx)\xb1i(Fyz+Fzy)]$ |

(2, ±2) | $+12[(Fxx\u2212Fyy)\xb1i(Fxy+Fyx)]$ |

(l, m)
. | F^{(l,m)}
. |
---|---|

(0, 0) | $\u221213(Fxx+Fyy+Fzz)$ |

(1, 0) | $\u2212i2(Fxy\u2212Fyx)$ |

(1, ±1) | $\u221212[(Fzx\u2212Fxz)\xb1i(Fzy\u2212Fyz)]$ |

(2, 0) | $+23[Fzz\u221212(Fxx+Fyy)]$ |

(2, ±1) | $\u221312[(Fxz+Fzx)\xb1i(Fyz+Fzy)]$ |

(2, ±2) | $+12[(Fxx\u2212Fyy)\xb1i(Fxy+Fyx)]$ |

(l, m)
. | T^{(l,m)} (, a), from Cartesian
. b | T^{(l,m)} (, a), from polar
. b |
---|---|---|

(0, 0) | $\u221213(axbx+ayby+azbz)$ | $\u221213azbz+12(a+b\u2212+a\u2212b+)$ |

(1, 0) | $+i2(axby\u2212aybx)$ | $\u2212122(a+b\u2212\u2212a\u2212b+)$ |

(1, ±1) | $+12[(azbx\u2212axbz)\xb1i(azby\u2212aybz)]$ | $\u221212(a\xb1bz\u2212azb\xb1)$ |

(2, 0) | $+23azbz\u221212(axbx+ayby)$ | $+23azbz\u221214(a+b\u2212+a\u2212b+)$ |

(2, ±1) | $\u221312[(axbz+azbx)\xb1i(aybz+azby)]$ | $\u221312(a\xb1bz+azb\xb1)$ |

(2, ±2) | $+12[(axbx\u2212ayby)\xb1i(axby+aybx)]$ | $+12a\xb1b\xb1$ |

(l, m)
. | T^{(l,m)} (, a), from Cartesian
. b | T^{(l,m)} (, a), from polar
. b |
---|---|---|

(0, 0) | $\u221213(axbx+ayby+azbz)$ | $\u221213azbz+12(a+b\u2212+a\u2212b+)$ |

(1, 0) | $+i2(axby\u2212aybx)$ | $\u2212122(a+b\u2212\u2212a\u2212b+)$ |

(1, ±1) | $+12[(azbx\u2212axbz)\xb1i(azby\u2212aybz)]$ | $\u221212(a\xb1bz\u2212azb\xb1)$ |

(2, 0) | $+23azbz\u221212(axbx+ayby)$ | $+23azbz\u221214(a+b\u2212+a\u2212b+)$ |

(2, ±1) | $\u221312[(axbz+azbx)\xb1i(aybz+azby)]$ | $\u221312(a\xb1bz+azb\xb1)$ |

(2, ±2) | $+12[(axbx\u2212ayby)\xb1i(axby+aybx)]$ | $+12a\xb1b\xb1$ |

*F*^{(l,m)} and *T*^{(l,m)} satisfy the symmetry relations

All *T*^{(l,m)} are real-valued in the conventional Zeeman basis.

The frame transformation of the spherical tensors is as follows:

where we have used *D*^{l}(*Ω*_{AB}) = *D*^{l†}(*Ω*_{BA}) with *Ω*_{BA} = (*α*, *β*, *γ*) as in Fig. 2 and *Ω*_{AB} = (−*γ*, −*β*, −*α*).

The Wigner matrix for the composition of two transformations is given in terms of the Wigner matrices of the component transformations by

### APPENDIX B: OPERATOR MATRIX ELEMENTS IN THE WIGNER BASIS

Using the shorthand $\xi i=Li,Mi,Ki$ [see Eq. (29)], the matrix elements for the angular-momentum operators are^{38,39}

where *J*_{±} = *J*_{x} ± i*J*_{y}. *x*, *y*, and *z* refer to a body-fixed frame and $\delta LM=\delta L1,L2\delta M1,M2$. Note the negative sign in the equation involving *J*_{z}.

The matrix elements of $DM,KL$ are

where the expressions in parentheses are Wigner 3-j symbols.^{38,39} Due to the selection rules for the 3-j symbols, the matrix elements can be non-zero only if |*L*_{1} − *L*_{2}| ≤ *L* ≤ *L*_{1} + *L*_{2} and *M* = *M*_{1} − *M*_{2} and *K* = *K*_{1} − *K*_{2}. They are all real-valued and possess the symmetry

Additionally, the following matrix elements are useful:

### APPENDIX C: DETAILS ABOUT THE DIFFUSION OPERATOR

The matrix elements of the isotropic part of the diffusion operator in the *LMK* basis are

with *R*_{d} = (*R*_{x} − *R*_{y})/4, *R*_{⊥} = (*R*_{x} + *R*_{y})/2, and $cL,K\xb1=L(L+1)\u2212K(K\xb11)$.

The expansion coefficients $X\u0303M,KL$ in Eq. (15) are given by

## REFERENCES

_{2}-Aib-OMe in acetonitrile

_{10}-helical peptides with TOAC nitroxide spin labels

^{1}H NMRD profiles of DOTA-Gd derivatives by means of the slow motion theory