Fokker-Planck formalism in magnetic resonance simulations

This paper presents an overview of the Fokker-Planck formalism for non-biological magnetic resonance simulations, describes its existing applications and proposes some novel ones. The most attractive feature of Fokker-Planck theory compared to the commonly used Liouville - von Neumann equation is that, for all relevant types of spatial dynamics (spinning, diffusion, flow, etc.), the corresponding Fokker-Planck Hamiltonian is time-independent. Many difficult NMR, EPR and MRI simulation problems (multiple rotation NMR, ultrafast NMR, gradient-based zero-quantum filters, diffusion and flow NMR, off-resonance soft microwave pulses in EPR, spin-spin coupling effects in MRI, etc.) are simplified significantly in Fokker-Planck space. The paper also summarises the author's experiences with writing and using the corresponding modules of the Spinach library - the methods described below have enabled a large variety of simulations previously considered too complicated for routine practical use.


Introduction
Good theory papers have two essential features: they are readable and computable.That is the reason why the Liouville -von Neumann equation (on the coherent side) and Bloch-Redfield-Wangsness / Lipari-Szabo theories (on the relaxation side) dominate magnetic resonance -there are papers and books that describe them with the eloquence and elegance of a well written detective story [1][2][3][4].
In this respect, Fokker-Planck theory is a hard sell.Its applications to magnetic resonance feature impractically large matrices [5], highly technical derivations [6], and very flexible dynamical models that require a level of knowledge few people would ever have about their systems [7].Yet the prize is tempting -Fokker-Planck theory captures, in one general equation, almost everything there is to model in NMR, EPR and MRI.Most importantly, it treats spatial dynamics at the same conceptual level as spin dynamics, and includes relaxation processes in a natural way that is free of perturbative assumptions [8].
The problem with the status quo is that spatial degrees of freedom are too often an afterthought in the Liouville -von Neumann formalism.For condensed phase spin systems it is commonly assumed that the spin state has no influence on spatial dynamics, but that spatial dynamics has an effect on the spin Hamiltonian [1,2].In other words, it is assumed that "something happens" to spatial coordinates that makes the spin Hamiltonian time-dependent.The resulting equation of motion is where   t ρ is the density operator and   t H is the spin Hamiltonian commutation superoperator.When the Hamiltonian contains stochastic terms (for example, from rotational diffusion in liquids), the effect of those terms is represented by a relaxation superoperator R [1,3,9] that can be modified to drive the solution to some thermal equilibrium state [10,11]: where the overbar indicates ensemble averaging and K is the kinetics superoperator that accounts for the possible presence of chemical processes in the system [2].This equation is currently the central pillar of most magnetic resonance simulation frameworks [12][13][14][15][16][17].It is deterministic, and many methods exist (Floquet theory [18], time slicing [15], COMPUTE [19], etc.) for calculating its exact solution analytically or numerically.The biggest source of complications here is that the "invisible hand" of spatial dynamics makes even the ensemble-averaged spin Hamiltonian time-dependent in non-trivial ways.This can lead to spectacularly complex analytical solutions -the excellent Equation 38 in the recent paper by Scholz, Meier and Ernst [20] is a good indication that perhaps the time has come to take a closer look at a numerical formalism that promises to avoid it.
The most attractive feature of Fokker-Planck theory compared to the Liouville -von Neumann formalism is that, for all common types of spatial dynamics (spinning, diffusion, flow, etc.), the corresponding Fokker-Planck evolution generator is time-independent.It is demonstrated below that many difficult magnetic resonance simulation problems (multiple rotation NMR, gradient chirps, ultrafast NMR, soft off-resonance microwave pulses in EPR, overtone NMR, etc.) are simplified significantly in Fokker-Planck space.It also produces major generalisations and simplifications across the simulation code -recent versions of Spinach [12] owe their versatility to the Fokker-Planck formalism.

Fokker-Planck equation
The state of an individual system in classical physics can be described by a vector of state variables that we will denote x .The associated equation of motion involves time derivatives and other operators acting on this vector.This formalism goes back to Newton [21], but describing an ensemble of systems in this way is technically difficult because the number of state variables becomes very large.A less detailed, but more convenient description of a classical ensemble may be formulated in terms of the probability density   , p t x of finding a system in a state x at time t .As the ensemble evolves, this probability flows around the state space in a way that depends on the interactions present.The resulting probability flux   ,t j x is determined by the local velocity   and the velocity depends on the equation of motion.
The Fokker-Planck equation, proposed independently by Adriaan Fokker [22] and Max Planck [23], can be formulated as the continuity equation for the probability flux: It essentially means that probability cannot be destroyed or created.It can only be moved around: the local rate of decrease in the probability density is equal to the divergence of its flux.

Spin degrees of freedom
In situations where spin degrees of freedom are present in the system, Equation (4) may be extended to include the quantum mechanical density matrix ρ as a state variable: The velocity in the spin space is given by Equation (1); the velocity in the lab space is determined by whatever is happening there.The average density matrix   ,t ρ x at every point x in the lab space is calculated by taking an integral over the probability density with respect to the spin degrees of freedom: We need the equation of motion for this quantity; it is obtained by differentiating Equation ( 6) with respect to time, then using the expression for the derivative from Equation ( 5), and then going through a few rounds of simplifications.That is a surprisingly convoluted procedure -Appendix A in [5] contains a detailed discussion.Here we would simply present the final result: R K is the Liouvillian that is responsible for spin dynamics and  

,t M x
is the spatial dynamics generator that controls diffusion, flow, sample spinning, and other types of classical mechanics in the laboratory space.This is the equation of motion that will be referred to as the "Fokker-Planck equation" for the rest of this paper.In the special case when

 
,t M x is the diffusion operator, Equation ( 7) is known as the stochastic Liouville equation [24].

Spatial dynamics generators
Expressions for   ,t M x come directly from classical mechanics [25].A few relevant equations of motion for the probability density are (the spatial dynamics generator that we need is in square brackets): 1. Flow along a coordinate x and circular motion with a phase   where v is linear velocity and  is angular velocity.
2. Circular motion around a specific axis where   , , J J J are angular momentum operators and  is the angular velocity.

Translational diffusion in a potential  
U r in an isotropic medium where  is the amplitude of the fluctuating force that generates the diffusion and  is the friction constant of the medium.

Isotropic rotational diffusion in three dimensions
where  is a shorthand for Euler angles or other orientation parameters,   X Y Z , , J J J are angular momentum operators and R D is the rotational diffusion coefficient.
Probability density evolution equations for other types of motion may be found in the specialist literature [26,27].For our purposes here, these equations yield the spatial dynamics generator that is to be inserted into Equation (7).For example, the stochastic Liouville equation used in electron spin resonance [24] is obtained by using the rotational diffusion term: Although spatial dynamics generators can depend on time, they are usually static in practice.One important thing to note is that this includes not only diffusion, but deterministic motion as well.For example, magic angle spinning and phase rotation of an external radiofrequency field both have time-independent generators within the Fokker-Planck formalism.
Quantum mechanical processes involving spatial degrees of freedom, such as quantum rotor [28] and quantum oscillator [29] dynamics, may be accounted for by using the Schrödinger Hamiltonian as the generator.In that case, the Fokker-Planck equation yields back the Liouville -von Neumann formalism.

Matrix representations of spatial dynamics generators
Analytical solutions to Equation ( 7) can be complicated and uninviting; that is the primary reason for the frosty reception that the Fokker-Planck formalism has seen so far in the magnetic resonance community.
In this brief overview we would argue that numerical solutions, where spatial coordinates are discretised on finite grids [30], are easier to understand, implement and use because all differential operators of the kind shown in the previous section become matrices.
Of course, selecting a basis of continuous functions would also lead to a matrix representation -for example, Wigner D-functions are a popular choice for rotational dynamics problems [6,8,24].Finite grids are advocated here because they are easier to automate -there are standard libraries that generate differential operators on any given grid.For uniform grids with periodic boundary conditions, Fourier differentiation matrices [31] are recommended: Operators for higher derivatives are obtained by taking powers of this matrix.For non-periodic coordinates or non-uniform grids, finite difference matrices are convenient -a general algorithm for building them has been published by Fornberg [32].Implementations of Equation ( 13) and of the Fornberg method are available in Spinach [12].Fourier differentiation matrices are an exact representation of the derivative operator on a given uniform grid -using them for rotational dynamics is important because magic angle spinning NMR simulations can go through thousands of rotor cycles.The number of discretisation points required for each grid is dictated by the accuracy required.Formal error bounds are available [31,32], but in practice the grid size is simply increased until convergence is achieved on the output.

Algebraic structure of the composite problem
From the algebraic perspective, the Fokker-Planck equation operates in a direct product of spatial and spin coordinates.For logistical reasons (to do with data layouts in Matlab), it is preferable to have spatial coordinates as the first term in the Kronecker product and spin coordinates as the second term.Thus, the general form of a purely spatial operator is  A 1, the general form of a pure spin operator is  1 B and the form of an operator that correlates spatial and spin degrees of freedom (for example, a pulsed field gradient) is  A B. The spin subspace has the usual direct product structure [1,2,15] and for the spatial coordinates it is reasonable to put all periodic coordinates ahead of all non-periodic ones.All of this results in the following product algebra for the evolution generators where the first term in the direct product accounts for the (possibly multiple) radiofrequency irradiations and uniaxial sample rotations, all generated by  

 so
Lie algebras, the second term is responsible for translational and diffusional dynamics generated by the general linear algebra   k  gl acting on a k - dimensional real space, and the last term is the usual special unitary algebra   n  su generating the dynamics of an n -state quantum system [33].In human language:

Fokker-Planck space
Liouville space  spatial degrees of freedom where the Liouville space [34][35][36] may be formally defined as the vectorisation of the density operator space acted upon by the adjoint Lie algebra of   n  su .Once the matrix representations for all operators are built, the Fokker-Plank equation of motion has the following general form: The process of recovering the average spin density matrix from the Fokker-Planck state vector is particularly straightforward in Matlab: the matrix is reshaped to have spatial degrees of freedom along one dimension and spin degrees of freedom along the other, and summed over the spatial degrees of freedom.
The direct product structure in Equation ( 14) can produce matrices of large dimension.As an example, using ten points in the RF phase grid, 100 slices in the sample position grid and a five-spin system (1024 × 1024 matrices in Liouville space) would produce a Fokker-Planck matrix dimension of 10 × 100 × 1024 ≈ 10 6 .Because finite difference matrices and spin operators are very sparse (there are hard guarantees for spin Hamiltonians in the Pauli basis [37]), this is not a problem -as the examples given below demonstrate, modern state space reduction [38][39][40][41], sparse matrix manipulation [42] and diagonalization-free simulation [43,44] techniques are sufficient.

Case study 1: magic angle spinning NMR
Fokker-Planck theory is the primary working formalism for MAS NMR simulations in Spinach [12].Two different formulations are possible [5].In Cartesian coordinates in three dimensions, the spatial dynamics generator comes from Equation ( 9): where  are Euler angles or any other directional parameters, and    H is the spin Hamiltonian commutation superoperator.To make the rotation generator J turn the spin sys- tem around the magic angle, the axis n must be directed appropriately, for example   Alternatively, we could observe that, once the spinning axis direction is chosen, the motion is actually uniparametric -the only spatial variable that changes with time is the rotor phase  .This leads to a more elegant formulation with a simpler spatial dynamics generator coming from Equation ( 8): The physical meaning of    is that of a phase increment generator -it is responsible for moving the rotor phase forward at the rate MAS Selecting a uniform periodic phase grid   j  makes this operator a constant matrix given by Equation (13) and results in the following expression for the spin Liouvillians at the phase grid points: , , , , , , where the 25 irreducible spherical components km Q of the Hamiltonian are defined in our earlier paper dealing with a large-scale implementation of Redfield theory [44].The best practical way of obtaining them is to use the Hamiltonian generation function of Spinach [12] that returns them as a cell array.The three sets of Wigner D-matrices accomplish the following roles.

 
lab2rot n D then rotates the result into the orientation that corresponds to the rotor direction vector n .Finally, the matrix representation for the rotor turning generator    is given by Equation ( 13).When the system dynamics permits,  angle averaging may be performed efficiently by averaging over the rotor phase in the same way as happens in Floquet theory [45].Only a two-angle spherical grid is then required for powder averaging.
A schematic illustration of the events taking place under (A) Equation (16), where the spatial dynamics operator generates a population flux through the entire orientation set and every point effectively hosts a Liouville -von Neumann equation for its own angles, and (B) Equation (18), where the rotor turning operator cycles the populations through the rotor orientations at every point on a two-angle spherical grid.Both approaches may be shown to be algebraically equivalent to, and about as fast computationally as, Floquet theory [5,18], with the notable difference that perturbative corrections to the rotating frame approximations are easier to implement within the Fokker-Planck formalism -Spinach [12] permits arbitrary-order corrections.
An important feature of the operators in square brackets in Equations ( 16) and ( 18) is that they are timeindependent and work as a part of the background Hamiltonian -in the Fokker-Planck picture, magic angle spinning effectively looks like a static interaction.All rotor cycle time slicing concerns intrinsic to the Liouville -von Neumann approach to MAS NMR simulation simply disappear.The practical convenience of this fact is hard to overstate -a good real-life example is given in our recent paper [46].  1N nucleus, assuming proton decoupling.Full simulation details are given in our recent paper [46].The ellipsoid plot illustrates the orientation of the 14 N NQI tensor relative to the molecular geometry in the glycine zwitterion.
Equation ( 16) has a further advantage -if the entire three-angle manifold is considered and the spatial basis set is chosen (following Freed's treatment of the mathematically related rotational diffusion process [8]) to be Wigner D-functions of those angles, Equation ( 16) does not require a powder averaging grid: the spherical average is obtained by taking the coefficient in front of   0 0,0 D .Given the titanic amount of work that has gone into optimising spherical grids in NMR [47][48][49][50][51][52][53][54], this is remarkable.Further details on this "grid-free" formalism are given in our recent paper [5] -because matrix dimensions go up significantly, it may be viewed as striking a different balance between memory utilisation and CPU time.
At the numerical implementation level, the evolution generator is structured as a sparse block-diagonal matrix with each block corresponding to a point on the rotor phase grid.The job of the rotor turning generator is to move block populations around during evolution as illustrated in Figure 1B.The matrix form of Equation ( 18) at each point of the spherical grid is: where N is the number of points in the rotor phase grid, the expression for the rotor turning generator matrix      is given in Equation ( 13) and 1 is a unit matrix of the same dimension as the spin Hamiltonian superoperator.External time-dependent events, such as shaped pulses, should be accounted for by adding the corresponding Hamiltonian to each block of the Liouvillian.Any rotating frame transformations or other Hamiltonian modifications should likewise be applied to each block.An extensively annotated numerical implementation of Equation ( 21) is available in the singlerot.mcontext function of Spinach [12], which performs all required procedures automatically.
It could be argued that a time-independent rotor turning generator is also a feature of Floquet theory [45].The two approaches do not actually compete -Floquet theory is obtained from Equation (18) when the basis set is chosen to be complex exponentials of the rotor phase (a formal proof is given in [5]), meaning that Floquet theory is a special case of the Fokker-Planck formalism.The advantage of the finite grid approach is that Equation ( 21) does not become any more complicated when second-and higherorder corrections to the rotating frame transformation are required -the corrections simply need to be applied to every   j  L block individually.
A good applications example is given by quadrupolar overtone simulations, where high-order perturbative corrections to rotating frame transformations are inevitable and a rather complicated effective Hamiltonian treatment must be performed to make it possible to simulate the spectrum in reasonable time [46].
Figure 2 illustrates the non-trivial dependence of the quadrupolar overtone MAS NMR spectrum on the spinning direction of the MAS rotor [55,56] -the sideband intensity pattern is reflected when the spinning direction is changed from clockwise to counter-clockwise.Further simulation details are given in [46], the extensively annotated Matlab code is a part of the standard example set supplied with Spinach [12].

Case study 2: multiple angle spinning NMR
The equation of motion for spin dynamics in a multiple angle spinning NMR experiment [57,58] is a simple extension of Equation ( 18).In the case of double rotation, we have two rotor turning generators: Ln n ρ (22) where O n is the direction vector specifying the orientation of the outer rotor relative to the laboratory frame of reference, I n is the direction vector specifying the orientation of the inner rotor relative to the outer rotor, O,I  are outer and inner rotor phases, O,I  are outer and inner rotor angular frequencies, ρ is the state vector, and L is the spin dynamics Liouvillian comprising the superoperators for free evolu- tion, relaxation theory and chemical kinetics: , , , , , , , The spin Hamiltonian commutation superoperator is: where the 25 irreducible spherical components km Q of the Hamiltonian [44] are the same as those in Equation (20).The five sets of Wigner D-matrices in Equation ( 24) accomplish the following transformations.are given by Equation ( 13).Here it is also possible to perform inexpensive  angle and relative rotor phase averaging efficiently by averaging over the phase grid points of both rotors, meaning that only a two-angle spherical grid is required in many cases.Setting Equation ( 24) up numerically is a significant undertaking and the reader is therefore advised to inspect and use the extensively annotated double- rot.m context function of Spinach [12] that performs this procedure automatically.

 
The matrix form of Equation ( 22) looks similar to Equation ( 21), but with two rotor turning operators: , where O N is the number of phase grid points for the outer rotor, I N is the number of phase grid points   14 N overtone NMR spectrum of glycine, reproduced from the paper by Carravetta and co-authors [59].The simulation (A) was carried out using ideal pulses, but due to the low transition moment across the forbidden overtone transition, it is instrumentally impossible to excite more than one line at a time.Experimental overtone spectra (B,C) were therefore recorded selectively at the most intense spinning sidebands.The ellipsoid plot illustrates the orientation of the 14 N NQI tensor relative to the molecular geometry in the glycine zwitterion.
Much as in the case of MAS described in the previous section, the increase in the matrix dimension is compensated for by the convenience of having a time-independent rotor turning generator -the 25) is just another constant matrix.From the programming point of view, setting up a DOR pulse sequence is at this point no harder than setting it up for liquid state NMR.The fact that a complicated spinning process is present no longer generates any housekeeping difficulties because it is a part of what looks like a static background Hamiltonian.
A sophisticated DOR simulation example from our recent paper [59] is shown in Figure 3. Quadrupolar overtone simulations are slow [55,56], and the most efficient way of accelerating them (Equations 6 and 7 in [46]) requires the evolution generator to be time-independent, which it would not be in the Liouville -von Neumann formalism.However, Equation ( 22) fits that requirement perfectly.

Case study 3: spatio-temporal NMR experiments
A number of advanced magnetic resonance experiments (Thrippleton-Keeler zero-quantum filter [60], ultrafast NMR sequences [61], pure shift NMR sequences [62], slice selection in MRI [63], etc.) use frequency-modulated pulses in the presence of magnetic field gradients.A significant source of magnetization loss in such experiments is the inevitable presence of diffusion and hydrodynamics that reduce the efficiency of gradient refocusing steps.Good theoretical models of this process must simultaneously account for coordinate-dependent quantum mechanical spin evolution, diffusion and hydrodynamics.That is not a simple task within the Liouville -von Neumann formalism, where formally correct but computationally expensive Monte-Carlo sums over stochastic trajectories have so far been used [64].
In the simplest case, the Fokker-Planck equation for spin evolution under a radiofrequency pulse in the presence of a magnetic field gradient, translational diffusion and uniform flow in one dimension would borrow spatial dynamics generators from Equations ( 8) and ( 10): where T D is the translational diffusion coefficient, v is the flow velocity, z is the sample coordinate,  is the phase of the radiofrequency pulse and   that contains the time-independent Hamiltonian 0 H (centre point chemical shifts, J-couplings and other static interactions), the radiofrequency term with initial phase 0  and a time-dependent amplitude   a t , a gradient term with a time-dependent amplitude   Z g t , the chemical kinetics superoperator K and the relaxation superoperator R .Here, the advantage of Equation ( 26) over the equivalent Liouville -von Neumann equation formulation is not in the elimination of time dependence, but in the fact that diffusion and flow are constant operators working in the background.
For numerical calculations using Equation ( 26) to become possible, matrix representations must be obtained for all of its constituent operators.Spin operators are matrices by definition.The differential operators acquire matrix representations when a grid of points is chosen for z and  .Finite difference ma- trices [32] are sufficient for the spatial coordinate and Equation ( 13) is recommended for the radiofrequency phase derivative operator.The resulting matrices must be projected into the Fokker-Planck space by taking Kronecker products with unit matrices on all coordinates that are unaffected by the operator.
The procedure results in the following matrix representation for Equation (26): in which M is the number of points in the RF phase grid, N is the number of points in the coordinate grid and K is the dimension of the spin state space.The block structure of the spin Liouvillian is: The responsibility for using sufficient number of grid points along both coordinates rests with the userin practice the number of grid points is increased until the answer stops changing.At the end of the calculation, a partial trace over the RF phase coordinate and the sample position coordinate (i.e. a sum of all

 
, k m z  blocks in ρ) produces the average spin state that is visible to spectrometer coils.Equation ( 26) has the following advantages over other equivalent formulations: 1.The radiofrequency pulse appears in frequency-amplitude, rather than Cartesian, parameterisation.The difference between the two is illustrated in Figure 4.It is clear that frequency-swept pulses are easier to discretise in frequency-amplitude coordinates.
2. Diffusion and flow generators are constant matrices that may simply be viewed as a part of the background Hamiltonian.For all practical purposes, the calculation looks like a standard simulation of a shaped radiofrequency pulse in NMR spectroscopy.Expensive Monte-Carlo averages over stochastic trajectories are unnecessary.
3. Diffusion and flow are non-periodic motions -while it could be argued in the sections above that Floquet theory provides an established alternative for periodic processes in the spin Hamiltonian, that is not the case here.A practical example of Equation ( 26) enabling the simulation of a sophisticated spatio-temporal NMR experiment is ultrafast NMR spectroscopy [61,65].The processes that make ultrafast NMR possible take place in the direct product  of the one-dimensional real space  and the spin state space : Once the spatial grid is chosen, the Fokker-Planck equation of motion for an ensemble of spin systems diffusing and flowing in one spatial dimension in the presence of a magnetic field gradient acquires the following matrix form: where H is the spin Hamiltonian commutation superoperator,   Z g t is the gradient amplitude, Z is the spatial coordinate operator (a diagonal matrix with the coordinates of grid points on the diagonal), k  are magnetogyric ratios of the spins in the system,   Z k S are the corresponding longitudinal spin commutation superoperators, T D is the diffusion coefficient, v is the flow velocity, Z D is a matrix represen- tation of the finite difference first derivative operator on our grid [32], R is the spin relaxation superop- erator, K is the chemical kinetics superoperator, 1  is a unit spin superoperator and 1  is the unit op- erator on the spatial grid.This equation may be solved using standard time propagation methods [66].
Figure 5 demonstrates the effect of diffusion and flow on the ultrafast COSY spectrum of a simple threespin system.The J-modulation effects -a quantum mechanical phenomenon that essentially requires full density matrix treatment in the spin subspace -are clearly visible [65].Another interesting conclusion is that fast diffusion and fast flow can cancel each other's deleterious effects to some extent.This is likely because they compensate each other for some small fraction of the sample.
A three-dimensional version of Equation ( 26) would include anisotropic diffusion and hydrodynamics [67].
The diffusion term is just T T

   D
, where T D is the translational diffusion tensor.However, hydrody- namics equations are in general non-linear and therefore undesirable within the linear structure of Equation (26).Their solution is a difficult task that requires specialised numerical methods [68]; the prospect of having to integrate a Navier-Stokes solver into a magnetic resonance package is dire.Fortunately, an elegant workaround is available because flows in magnetic resonance systems are usually stationary [69,70].A hydrodynamics solver may be run to obtain the stationary velocity field   v r which may then be combined with Equation ( 5) to yield the following extra term in the evolution generator: The operator in square brackets does not depend on time and therefore fits naturally into the background evolution operator of Equation (26).Its matrix representation is easy to obtain using finite difference matrices on a suitable grid.In the case of one-dimensional flow of constant velocity, the flow term from Equation ( 26) is produced.

Case study 4: orientation selection in double-electron resonance
Another relevant example where the Fokker-Planck approach is advantageous is the off-resonance irradiation and orientation selection that take place during DEER experiments [71].Because large zero-field splittings and soft microwave pulses at two distinct frequencies may be involved [72], Liouville -von Neumann formalism simulations must proceed by explicitly discretising both waveforms and running a torturously slow simulation with a time-dependent Hamiltonian in the laboratory frame.However, the Fokker-Planck evolution generator, coming from Equation (18), is again time-independent: where the Liouvillian now includes the background Hamiltonian and the microwave irradiation with the phase  that is incremented at the frequency MW  : where 0 H is the spin Hamiltonian containing all interactions intrinsic to the spin system,   X Y , S S are electron spin superoperators, a is the amplitude of the microwave pulse and 0  is its initial phase.A matrix representation of this problem is obtained by discretising    using Equation (13): in which M is the number of points in the microwave phase grid (fewer than ten are required in practice) and K is the dimension of the spin state space.The block structure of the spin Liouvillian is: Equation ( 33) should be invoked for each of the three microwave pulses; the usual Liouville -von Neumann equation is sufficient for the free evolution period.A partial trace over the microwave phase grid at the end of each pulse produces the resulting spin density matrix.The DEER trace has a modulation depth of approximately 10% of the total echo intensity.The associated Matlab source code is available in the example set supplied with versions 1.8 and later of Spinach library [12].
Equation ( 33) is currently the default formalism for soft-pulse DEER simulation in Spinach [12] -the corresponding functions are extensively annotated and should be consulted for the details of the numerical implementation.One of the examples included with the package is detailed in Figure 6.It is clear that the description of the soft pulses is very accurate -even the sinc wiggles are visible.

Case study 5: overtone cross-polarisation under MAS
Another good example of the Fokker-Planck formalism dealing successfully with a combination of spatial motion, quantum mechanics and radiofrequency irradiation is quadrupolar overtone cross-polarisation.
Even basic overtone NMR simulations are difficult because they involve a radiofrequency that cannot be assumed to be close to the Zeeman frequency of the nucleus in question [46,55,56,59].This impedes the first step normally taken in magnetic resonance simulations -the rotating frame transformation -and makes Liouville -von Neumann equation simulations slow.By that standard, cross-polarisation overtone MAS simulations are the stuff of nightmares -there are three non-commensurate and very different frequencies (MAS, 1 H, 14 N overtone) in the time-dependent part of the Hamiltonian.The system must be propagated for several milliseconds in the laboratory frame and the entire calculation averaged over a three-angle spherical grid.The situation may be simplified somewhat by going into the partial rotating frame with respect to the proton Zeeman frequency and taking care to retain the pseudosecular terms in the various couplings involving protons, but the two remaining frequencies are still problematic.
The advantage of the Fokker-Planck formulation in this context is threefold.Firstly, it makes time dependence disappear completely from the evolution generator because phases of the RF and the rotor are now incremented by time-independent derivative terms introduced in Equations ( 26) and (18).Secondly, only a two-angle spherical grid is required because rotor phase averaging is built-in.Thirdly, the lack of explicit time dependence in the evolution generator permits frequency domain detection which, for overtone peaks, is much faster than the time domain approach [46,59].The equation of motion is: with the following spin Liouvillian (assuming, to match the available experimental literature, 14 N to be the quadrupolar nucleus and 1 H to be the polarisation source): , , , cos cos sin sin cos sin where , ρ is the density matrix,  is the magic angle, MAS  is the MAS rotor phase, N RF  is the nitrogen radiofrequency phase, MAS  and N RF  are the corresponding frequencies, 0

H
is the orientation-independent part of the Hamiltonian commutation superoperator (J-couplings, isotropic chemical shifts, etc.), km Q are the 25 irreducible spherical components [44] of the anisotropic part (dipolar interactions, quadrupolar interactions, chemical shift anisotropies, etc.), The matrix representation for Equation (37) is built in a similar way to the previous sections, by discretising the RF phase and the rotor phase on finite grids and using Equation ( 13) to obtain matrix representations for the phase derivative operators.In practical calculations on 1 H-14 N two-spin systems in glycine and Nacetylvaline (the experimental paper will be published in due course), 15 points were found to be sufficient for the rotor phase grid and 5 points for the RF phase grid, meaning that the total spatial dynamics subspace dimension is 75.Given the Liouville space dimension of 36 for the spin subspace, the total dimension for the matrix representation of Equation ( 37) is 2700, which is easy to handle using modern sparse matrix manipulation methods, and likely tensor train methods in due course [75,76].
Because the time evolution generator in Equation ( 37) is time-independent, time propagation for the duration of the cross-polarisation period may be carried out with a single matrix exponential: where F is the matrix representation of the Fokker-Planck evolution generator.Because the complexity of the scaling and squaring procedure [77,78] for matrix exponentiation is asymptotically logarithmic in the duration of the time interval CP t , the fact that is not in practice a problem [43].
Once the cross-polarisation pulse is over, the equation of motion becomes an instance of Equation ( 18): with a Liouvillian that no longer contains any radiofrequency terms: , , , where now . For each point on the spherical averaging grid, the matrix representation of this equation is identical to Equation ( 21).If we denote the Fokker-Planck evolution generator F , the following equation gives the Fourier transform of the free induction decay on 14 N: where 1 is a unit matrix of the same dimension as F .The matrix dimension at this stage is smaller than at the cross-polarisation stage because a partial trace over the RF degrees of freedom is taken at the end of the CP period and the Liouvillian in Equation ( 42) no longer contains RF terms.
A significant advantage of Equation ( 42) over time-domain detection is that very few points are required (a few hundred in the frequency domain compared to billions in the time domain) and that the sparse matrix-inverse-times-vector operation is very fast when modern iterative solvers, such as ILU preconditioned GMRES [79], are used.Because a time-independent evolution generator is a requirement, this is only possible within the Fokker-Planck formalism.

Potential further applications
This section contains a speculative overview of other types of magnetic resonance simulations and applications that could be enabled, or made simpler, by the Fokker-Planck formalism.One such area is optimal control, where the problem of the interpretation of numerically optimised microwave and radiofrequency pulses is significant because their immediate appearance is obscure [44].All current techniques rely either on the time-frequency representation of the pulse itself [80], or on the analysis of the spin system trajectory under the pulse [44].In both cases the analysis step is performed after the optimal control solution has been obtained.The Fokker-Planck formalism is potentially useful here because it makes it possible to formulate the commonly used gradient ascent pulse engineering (GRAPE) algorithm [81,82] directly in the frequency-amplitude representation.
The standard setting for the quantum optimal control problem [81] involves a system with the Hamiltonian partitioned into the "drift" part 0 H that we cannot influence, and the "control" part, in which the spectrometer can vary the coefficients   k c t in front of some operators k H [81]: Analytical or numerical optimization of "control sequences"   k c t to achieve some pre-defined experi- mental objectives is the subject of optimal control theory [83].In most magnetic resonance cases, the drift Hamiltonian contains the magnet Zeeman terms and the spin-spin coupling terms.The control channels correspond to radiofrequency and microwave irradiation: where the summation indices run over all spins in the system, k Z are Zeeman tensors, nk A are coupling tensors and

S S 
are spin operators.The current wisdom, illustrated well in [80], and also in Figure 4, is that Cartesian representation of the waveform Pulse amplitudes and frequencies enter this equation linearly and may therefore be optimised by the GRAPE procedure [82], as well as its recent enhancements [81,84,85].The first group of control operators consists of terms discretised on a finite phase grid: The second group of control operators contains spectral finite difference representations of the derivative operators obtained using Equation (13).With those matrix representations in place, Equation (45) acquires the form that is directly usable by the existing optimal control software [12,17,86].
In our testing so far, this frequency-amplitude representation for control sequences did not significantly outperform either the Cartesian or the phase-amplitude control version of GRAPE.It should, however, certainly be mentioned here and kept on the books until such time as a system presents itself for which amplitude-frequency coordinates would be in some sense natural.
Another interesting, if exotic, hypothetical case of a spatial variable acting as a control channel in magnetic resonance experiments is the rotor phase.As per Equation ( 18), the spinning operator enters the Fokker-Plank equation of motion for a MAS experiment in the same algebraic way as the radiofrequency phase enters Equation (45).The matrix representation of the rotor phase derivative may therefore be used as a control operator and the associated spinning rate as a control coefficient.Another control coefficient is provided by the directional cosine of the spinning axis.
While variable spinning angle experiments are firmly established in NMR [87,88] and some applications of variable spinning rate experiments do exist [89], it is at this point unclear whether the use of optimal control theory here would result in any improvement, or whether the hardware of sufficient agility to control both the angle and the rate in real time could be built at a reasonable financial and time cost.We have not explored this matter any further, but it bears notice that the possibility exists.

Conclusions and outlook
Our experience of designing and coding Spinach [12] indicates that a high level of abstraction in the fundamental equations of motion, exotic and unwieldy though they may at first appear, is worth it in the long run because the resulting framework is general, flexible, extensible and maintainable.Throughout the last ten years we had to resist significant temptation to introduce case-specific tweaks that could have made some simulations faster at the cost of fragmenting the code into incompatible chunks.The result of this holistic approach is a software package that simulates everything there is in magnetic resonance.One factor that made this possible is Matlab.The other is the Fokker-Planck equation.
A formalism that simultaneously supports orientation-selective DEER [72], ultrafast NMR [61,65], singlet state flow imaging [64] and gadolinium cross-effect DNP (see the example set supplied with Spinach [12]), with an option of including the most general spin relaxation theory in existence [8], is very valuable.It constitutes a sweeping generalisation that also supports multi-site exchange problems [90], spin isomers [91], quantised spatial degrees of freedom [28], contains Floquet theory [45] and stochastic Liouville equation [8] as special cases, and accommodates the whole body of theoretical methods for NMR spectroscopy of diffusion and transport in complex media [69,70].The basic principles have been known since the seminal papers and books by Fokker [22], Planck [23], Callaghan [69,70] and particularly Freed [8,24], but it was not until recently that the computers became powerful enough, matrix dimension reduction tools arrived [38][39][40][41] and sparse matrix libraries became sufficiently convenient.In our opinion, the Fokker-Planck equation should now be preferred to the Liouville -von Neumann formalism as the fundamental equation of motion for numerical simulations in magnetic resonance.
frame, in which the interactions are defined, into each of the individual crystallite orientations specified by the spherical averaging grid.result, around the Z axis, into the orientation that corresponds to the specified rotor phase.

Figure 2 .
Figure 2.An illustration of the dependence of the 14 N quadrupolar overtone MAS NMR spectrum on the rotor spinning direction, simulated using Equation (18) with C Q = 1.18 MHz,  Q = 0.53 and a chemical shift of 32.4 ppm for the14 N nucleus, assuming proton decoupling.Full simulation details are given in our recent paper[46].The ellipsoid plot illustrates the orientation of the14 N NQI tensor relative to the molecular geometry in the glycine zwitterion.

D
rotates the molecule from the reference frame, in which the interactions are defined, into each of the individual crystallite orientations specified by the spherical averaging grid.result, around the Z axis, into the orientation that corresponds to the specified inner rotor phase.result into the orientation that corresponds to the inner rotor direction vector I n relative to the outer rotor.
result, around the Z axis, into the orientation that corresponds to the specified outer rotor phase.Finally, rotates the result into the orientation that corresponds to the outer rotor direction vector O n relative to the laboratory frame of reference.Matrix representations for the rotor turning generators for the inner rotor, S N is the dimension of the spin Hamiltonian commutation superoperator, 1 are identity superoperators of indicated dimensions,      are Fourier spectral differentiation matrices of indicated dimensions obtained using Equation (13), state vectors in the order of increasing index of the inner rotor phase grid points, followed by the increasing index of the outer rotor phase grid points, and the block-diagonal matrix in the middle is obtained by concatenating individual phase grid point Liouvillians order as the state vectors.External timedependent events, such as shaped pulses, should be accounted for by adding the corresponding Hamiltonian superoperator to each matrix block.Any rotating frame transformations or other Hamiltonian modifications should likewise simply be applied to each block.

Figure 3 .
Figure 3. Double rotation (1425 Hz, 6950 Hz)14 N overtone NMR spectrum of glycine, reproduced from the paper by Carravetta and co-authors[59].The simulation (A) was carried out using ideal pulses, but due to the low transition moment across the forbidden overtone transition, it is instrumentally impossible to excite more than one line at a time.Experimental overtone spectra (B,C) were therefore recorded selectively at the most intense spinning sidebands.The ellipsoid plot illustrates the orientation of the14 N NQI tensor relative to the molecular geometry in the glycine zwitterion.

Figure 4 .
Figure 4. Cartesian (left) vs. frequency-amplitude (right) representation of a frequency-swept pulse with the frequency range of 2 MHz around the centre and a duration of 1 ms.

Figure 5 .
Figure 5.Ultrafast COSY[61] spectrum of a three-spin system (δ 1 = 3.70 ppm, δ 2 = 3.92 ppm, δ 1 = 4.50 ppm, J 1,2 = 10 Hz, J 2,3 = 12 Hz, J 1,3 = 4 Hz) as a function of the diffusion coefficient and the linear flow velocity.The spectrum was simulated for a sample distributed in one spatial dimension with a 500-point grid with absorptive boundary conditions used to discretise a 15 mm interval representing the active volume of a typical NMR sample.The ultrafast COSY sequence was recorded using 128 gradient readout loops, 512 acquisition points each, with a dwell time of 1 µs.Acquisition gradient amplitude was set to 0.10 T/m, encoding gradient amplitude to 0.01 T/m, and the coherence selection gradient amplitude to 0.47 T/m.The duration of the coherence selection gradient pulses was 1 ms.Spatial encoding was achieved with a WURST pulse with a smoothing factor of 40 and a bandwidth of 10 kHz acting over a 15 ms time interval and discretised over a 1000-point time grid.Technical details about ultrafast COSY simulations in Spinach are available in[65].

Figure 6 .
Figure 6.A three-pulse DEER trace and associated diagnostic information for a two-electron system, simulated as described in the main text, at a magnetic induction of 0.34518 Tesla with the following parameters: first electron g-tensor eigenvalues[2.284,2.123, 2.075], Euler angles [45°, 90°, 135°], second electron g-tensor eigenvalues [2.035, 2.013, 1.975], Euler angles [120°, 60°, 30°], inter-electron distance of 20 Angstrom along the X axis, pulse durations 20 ns, 40 ns and 40 ns, pulse amplitude 8MHz for all three pulses, pulse frequencies 9.623 GHz, 9.092 GHz and 9.623 GHz, 1 µs gap between the first and the third pulse, 100 steps in the position of the second pulse, 100 points in the echo sampling window.The DEER trace has a modulation depth of approximately 10% of the total echo intensity.The associated Matlab source code is available in the example set supplied with versions 1.8 and later of Spinach library[12].
D-functions of the three Euler angles specifying the crystallite orientation within the powder grid[73], Wigner D-functions defining rotor orientation and phase using angle-axis parameterisation[74] with MAS 2 are nitrogen Zeeman commutation superoperators, H RF a is proton RF amplitude and N RF a is nitrogen RF amplitude.
fact, this last one appears to be the most convenient[80].It follows from Equation (26) that the amplitude-frequency representation appears naturally within the Fokker-Planck formalism: