Dr. Ilya Kuprov Prof. Stefan Glaser & Dr. Tim Freegarde Philosophy Doctorate
University of Southampton School of Chemistry Spin Dynamics Group Faculty of
Natural and Environmental Sciences Magnetic Resonance A thesis submitted for
the degree of
Doctor of Philosophy

# Advanced Optimal Control Methods for Spin Systems

(October 2017)

###### Abstract

Work within this thesis advances optimal control algorithms for application to
magnetic resonance systems. Specifically, presenting a quadratically
convergent version of the gradient ascent pulse engineering method. The work
is formulated in a superoperator representation of the Liouville-von Neumann
equation.

A Newton-grape method is developed using efficient calculation of analytical
second directional derivatives. The method is developed to scale with the same
complexity as methods that use only first directional derivatives. Algorithms
to ensure a well-conditioned and positive definite matrix of second
directional derivatives are used so the sufficient conditions of gradient-
based numerical optimisation are met.

A number of applications of optimal control in magnetic resonance are
investigated: solid-state nuclear magnetic resonance, magnetisation-to-singlet
pulses, and electron spin resonance experiments.

###### keywords:

Optimal Control, Magnetic Resonance

####  Declaration of Authorship

I, David L. Goodwin, declare that the thesis entitled _Advanced Optimal
Control Methods for Spin Systems_ and the work presented in the thesis are
both my own, and have been generated by me as the result of my own original
research. I confirm that:

  * •

this work was done wholly or mainly while in candidature for a research degree
at this University;

  * •

where any part of this thesis has previously been submitted for a degree or
any other qualification at this University or any other institution, this has
been clearly stated;

  * •

where I have consulted the published work of others, this is always clearly
attributed;

  * •

where I have quoted from the work of others, the source is always given. With
the exception of such quotations, this thesis is entirely my own work;

  * •

with the oversight of my main supervisor, editorial advice has been sought. No
changes of intellectual content were made as a result of this advice.

  * •

I have acknowledged all main sources of help;

  * •

where the thesis is based on work done by myself jointly with others, I have
made clear exactly what was done by others and what I have contributed myself;

  * •

parts of this work have been published as [GOODWIN15, GOODWIN16, GOODWIN17,
GOODWIN18]

Signed: |  |  |  |
---|---|---|---|---
Date: |  | 1${}^{\text{st}}$ March 2018 |  |
Supervisor: |  | Dr. Ilya Kuprov |  |
Co-supervisor: |  | Dr. Marina Carravetta |  |
External examiner: |  | Prof. Steffen Glaser |  |
Internal examiner: |  | Dr. Tim Freegarde |  |

###### Contents

  1. Declaration of Authorship
  2. 0 Preface
  3. 1 Introduction
    1. 1 Thesis Outline
    2. 2 Early optimal control
    3. 3 Analytical optimal control
    4. 4 Broadband optimal pulses
    5. 5 Gradient based optimal control
      1. 1 Solid-state nuclear magnetic resonance
      2. 2 Nuclear magnetic resonance imaging
      3. 3 Quantum information
      4. 4 Electron spin resonance
    6. 6 Modified methods
      1. 1 Optimal tracking
      2. 2 Cooperative pulses
    7. 7 Other optimal control methods
    8. 8 Numerical computation
    9. 9 Summary
  4. 2 Magnetic Resonance Theory
    1. 1 Angular momentum
      1. 1 Schrödinger’s equation
      2. 2 Euclidean principle of relativity
      3. 3 Orbital angular momentum
      4. 4 The uncertainty principle
      5. 5 Spherical polar coordinates
      6. 6 Intrinsic angular momentum
    2. 2 Irreducible representations of the rotation group
    3. 3 Liouville-von Neumann equation
      1. 1 Superoperator algebra
      2. 2 Composite systems
    4. 4 The Spin Hamiltonian
      1. 1 Zeeman Interactions
      2. 2 Inter-Nuclear Interactions
      3. 3 Microwave and Radiofrequency terms
  5. 3 Optimal Control Theory
    1. 1 Numerical optimisation
      1. 1 The gradient ascent method
      2. 2 Newton’s Method
      3. 3 Quasi-Newton Methods
      4. 4 Line search subproblem
      5. 5 Comparison of Newton-type methods
    2. 2 Gradient Ascent Pulse Engineering
      1. 1 Bilinear systems
      2. 2 Piecewise constant approximation
      3. 3 The Pontryagin maximum principle
      4. 4 Fidelity measures
      5. 5 Fidelity gradient
      6. 6 Algorithm
  6. 4 Hessian Calculations
    1. 1 Fidelity Hessian
      1. 1 Hessian structure
      2. 2 Hessian elements
      3. 3 Control operator commutativity
    2. 2 Hessian Regularisation
      1. 1 Cholesky factorisation
      2. 2 Trust Region Method
      3. 3 Rational Function Optimisation
      4. 4 Hessian Conditioning
    3. 3 Convergence analysis
      1. 1 Scalar coupled three-spin state-to-state population transfer
      2. 2 Quadrupolar spin state-to-state population transfer
      3. 3 Universal rotation
  7. 5 Propagator Directional Derivatives
    1. 1 The Matrix Exponential
      1. 1 Taylor series approximation
      2. 2 Padé approximation
      3. 3 Chebyshev rational approximation
      4. 4 Scaling and squaring
    2. 2 Directional derivatives
      1. 1 Directional derivatives of the matrix exponential
      2. 2 Block triangular auxiliary matrix exponential
    3. 3 Directional derivatives of the Propagator
      1. 1 First propagator directional derivatives
      2. 2 Second directional derivatives
    4. 4 Krylov subspace techniques
    5. 5 Propagator recycling
      1. 1 Fidelity and gradient calculation
      2. 2 Fidelity, gradient and Hessian calculation
      3. 3 Parallelisation
  8. 6 Constrained Optimal Control
    1. 1 Change of variables
      1. 1 Transforms between Cartesian & Polar Coordinates
    2. 2 Power penalties
      1. 1 Norm-square penalty
      2. 2 Spillout-norm-square penalty
    3. 3 Smoothed Controls
      1. 1 Magnetisation-to-singlet transfer
      2. 2 Robust magnetisation-to-singlet transfer
      3. 3 Smooth, robust magnetisation-to-singlet transfer
  9. 7 Solid State Nuclear Magnetic Resonance
    1. 1 Floquet Theory
      1. 1 Powder average
      2. 2 Optimal control of a static powder average
      3. 3 Optimal control under magic angle spinning
    2. 2 The Fokker-Planck formalism
      1. 1 Optimal Control in HMQC-type experiment
  10. 8 Feedback Control
    1. 1 Feedback Control Optimisation
      1. 1 Hardware and Software
      2. 2 Optimisation Method
    2. 2 Experimental Application
      1. 1 Feedback Optimised Hahn Echo
      2. 2 Feedback Optimised OOP-ESEEM
  11. 9 Summary
  12. A Irreducible Spherical Tensor Operators
  13. B Arbitrary Waveform Generator
    1. 1 AWG Schematics
    2. 2 Xepr Python code
    3. 3 Feedback Control Results Tables

###### List of Figures

  1. 1 Finding a minimum of a function using tangents of its gradient plot
  2. 2 Robot manipulator
  3. 1 Spherical coordinates
  4. 2 Stern-Gerlach experiment
  5. 3 Euler angles
  6. 1 Types of stationary points
  7. 2 Comparison of optimisation methods on the Rosenbrock function
  8. 3 Piecewise-constant approximation
  9. 4 A representation of the adjoint state
  10. 1 Block Hessian diagram
  11. 2 Regularised and conditioned Hessian matrices
  12. 3 Control pulses for HCF state transfer
  13. 4 Trajectory analysis of HCF state transfer
  14. 5 BFGS method updating the Hessian and the inverse Hessian
  15. 6 Newton and BFGS methods with Cholesky factorisation
  16. 7 Newton and BFGS methods with regularisation and conditioning
  17. 8 Hessian conditioning iterations
  18. 9 HCF optimisation step lengths
  19. 10 Trajectory analysis of 14N with Newton-step
  20. 11 Convergence of Newton-GRAPE, without a line search.
  21. 12 Convergence of Newton-GRAPE method for a universal rotation.
  22. 1 Parallelism of gradient and Hessian calculations
  23. 1 Norm square penalty and its derivatives
  24. 2 Spillout norm square penalty and its derivatives
  25. 3 Spillout norm cube penalty and its derivatives
  26. 4 The robustness of the SLIC pulse
  27. 5 Pulse shape from an optimised M2S transfer
  28. 6 Smoothed pulse shape from an optimised M2S transfer
  29. 7 The robustness of the optimised M2S transfer
  30. 8 The robustness of the smoothed, optimised M2S transfer
  31. 1 Convergence of Newton-GRAPE method for glycine state transfer
  32. 2 Control pulses for glycine state transfer
  33. 3 Trajectory analysis of glycine state transfer
  34. 4 Trajectory analysis of HMQC-type experiment with hard pulse
  35. 5 Optimal control in an HMQC-type experiment
  36. 6 Phase profile of GRAPE pulses for 2, 4, 6, and 8 rotor periods
  37. 7 Trajectory analysis of GRAPE pulses for 2, 4, 6, and 8 rotor periods
  38. 8 Convergence analysis with 2, 4, 6, and 8 rotor periods
  39. 1 Closed-loop feedback control
  40. 2 Software flow
  41. 3 Simplex algorithms
  42. 4 Echo results with 33 points
  43. 5 Echo results with 21 points
  44. 6 Echo results with 11 points
  45. 7 Transformed 11 point shapes
  46. 8 OOP-ESEEM pulse sequence diagram
  47. 9 OOP-ESEEM from feedback-control
  48. 1 The Bruker SpinJet Schematics.
  49. 2 Convergence of Nelder-Mead and multidirectional search algorithms

@l@>p0.45l@r@ Symbol Description Notes Chapter
Symbol Description Notes Chapter

Magnetic Resonance:
$\langle\,\cdot\,|\,\cdot\,\rangle$ inner product §§ 1, 2, 3, 4, 6, 7
$t$ time interval §§ 2, 3, 5, 7
$i$ imaginary number $=\oldsqrt[\ ]{-1}=\frac{-1}{\oldsqrt[\ ]{-1}}$ §§ 2,
3, 4, 5, 7
$\big{|}\Psi(t)\big{\rangle}$ state of a system at a time $t$ § 2
$\big{|}\Psi_{i}\big{\rangle}$ complete orthonormal basis states (Hilbert
space) § 2
$c_{i}(t)$ time dependent coefficient of state § 2
$\hat{\mathcal{H}}$ Hamiltonian operator § 2
$\hat{R}$ rotation operation § 2
$\mathbf{A}$ (bold typeface) a matrix §§ 2, 3, 4, 5, 6
$\mathbf{A}^{T}$ transpose of a matrix $\mathbf{A}$ $A_{ij}^{T}=A_{ji}$ §§
2, 4
$\mathbf{A}^{\ast}$ complex conjugate of $\mathbf{A}$
$=\mathrm{Re}{(\mathbf{A})}-\mathrm{Im}{(\mathbf{A})}$ §§ 2, 3, 6
$\mathbf{A}^{\dagger}$ Hermitian conjugate of $\mathbf{A}$
$=\big{(}\mathbf{A}^{\ast}\big{)}^{\\!T}$ §§ 2, 3
$\mathbf{A}^{-1}$ matrix inverse of $\mathbf{A}$
$\mathbf{A}=\big{(}\mathbf{A}^{-1}\big{)}^{\\!-1}$ §§ 2, 3, 4, 5, 7
$\big{[}\mathbf{A},\mathbf{B}\big{]}$ commutator of matrices $\mathbf{A}$ and
$\mathbf{B}$ $=\mathbf{A}\mathbf{B}-\mathbf{B}\mathbf{A}$ §§ 2, 4, 5
$\big{|}\Psi\big{\rangle}$ arbitrary state vector §§ 2, 3, 4
$\big{|}\Psi^{\prime}\big{\rangle}$ transformed state vector
$\big{|}\Psi\big{\rangle}$ § 2
$\big{\langle}\Psi\big{|}$ Hermitian conjugate of $\big{|}\Psi\big{\rangle}$
$\big{\langle}\Psi\big{|}=\big{|}\Psi\big{\rangle}^{\\!\dagger}$ §§ 2, 3, 4
$\hat{r}$ displacement operator § 2
$\hat{p}$ momentum operator $=-i\hbar\nabla$ § 2
$\nabla$ gradient operator $=\Big{(}\frac{\partial}{\partial
x},\frac{\partial}{\partial y},\frac{\partial}{\partial z}\Big{)}^{\\!\\!T}$
§§ 2, 3, 4, 5, 6, 7
$\hbar$ reduced Planck constant (Js(rad)-1) $\approx 1.05\times 10^{-34}$ §
2
$\Big{(}x,y,z\Big{)}$ Cartesian coordinates §§ 2, 6
$\Big{(}r,\theta,\varphi\Big{)}$ spherical coordinates §§ 2, 6
$r$ radial distance $=\oldsqrt[\ ]{x^{2}+y^{2}+z^{2}}$ §§ 2, 6
$\theta$ inclination/polar angle $=\arccos{\frac{z}{r}}$ § 2
$\varphi$ azimuthal angle $=\operatorname{atan2}2{\frac{y}{x}}$ §§ 2, 4, 6
$\hat{L}$ orbital angular momentum operator $=\hat{r}\times\hat{p}$ § 2
$\hat{L}_{x}$ $x$-component of angular momentum operator
$=-i\big{[}\hat{L}_{y},\hat{L}_{z}\big{]}$ § 2
$\hat{L}_{y}$ $y$-component of angular momentum operator
$=-i\big{[}\hat{L}_{z},\hat{L}_{x}\big{]}$ § 2
$\hat{L}_{z}$ $z$-component of angular momentum operator
$=-i\big{[}\hat{L}_{x},\hat{L}_{y}\big{]}$ § 2
$\hat{L}^{2}$ Casimir operator
$=\hat{L}_{x}^{2}+\hat{L}_{y}^{2}+\hat{L}_{z}^{2}$ § 2
$\hat{L}_{+}$ raising operator $=\hat{L}_{x}+i\hat{L}_{y}$ § 2
$\hat{L}_{-}$ lowering operator $=\hat{L}_{x}-i\hat{L}_{y}$ § 2
$\Delta x$ change in the variable $x$ § 2
$\big{\langle}A\big{\rangle}$ expected value of $A$ § 2
$Y_{\\!\ell}^{m^{\vphantom{\prime}}}\\!(\theta,\varphi)$ spherical harmonics
§ 2
$\mathbb{N}$ natural numbers $\dots,-2,-1,0,1,2,\dots$ § 2
$\ell$ angular momentum quantum number $\ell\in\mathbb{N}$ § 2
$m$ magnetic quantum number $m=-\ell,\dots,\ell$ § 2
$P_{\\!\ell}(\xi)$ Legendre polynomials
$=\frac{1}{2^{\ell}\ell!}\frac{{\operatorname{d\\!}}^{\>\ell}}{\operatorname{d\\!}\xi^{\ell}}\Big{[}\big{(}\xi^{2}-1\big{)}^{\\!\ell}\Big{]}$
§ 2
$\big{|}\Psi(\vec{r},s)\big{\rangle}$ spinor
$=\\!\sum\limits_{n,k}\\!\big{|}\Psi_{\\!n}(\vec{r}\,)\big{\rangle}\\!\otimes\\!\big{|}\Psi_{\\!k}(s)\big{\rangle}$
§ 2
$s$ spin quantum number $s=0,\frac{1}{2},1,\frac{3}{2},2,\dots$ § 2
$\bm{\mathbbm{1}}$ identity matrix §§ 2, 3, 5, 6
$\mathbf{0}$ vector or matrix of zeros §§ 2, 3
$\hat{S}$ intrinsic angular momentum (spin) § 2
$\hat{S}_{x}$ $x$-component of spin operator
$=-i\big{[}\hat{S}_{y},\hat{S}_{z}\big{]}$ § 2
$\hat{S}_{y}$ $y$-component of spin operator
$=-i\big{[}\hat{S}_{z},\hat{S}_{x}\big{]}$ § 2
$\hat{S}_{z}$ $z$-component of spin operator
$=-i\big{[}\hat{S}_{x},\hat{S}_{y}\big{]}$ § 2
$\hat{S}_{+}$ spin raising operator $=\hat{S}_{x}+i\hat{S}_{y}$ § 2
$\hat{S}_{-}$ spin lowering operator $=\hat{S}_{x}-i\hat{S}_{y}$ § 2
$\big{|}\alpha\big{\rangle}$ $\alpha$ spin state
$=\Big{(}1,0\Big{)}^{\\!\\!T}$ § 2
$\big{|}\beta\big{\rangle}$ $\beta$ spin state $=\Big{(}0,1\Big{)}^{\\!\\!T}$
§ 2
$T$ total time §§ 2, 3
$\hat{\mathcal{H}}(t)$ time-dependent Hamiltonian §§ 2, 3
$\big{|}\psi(t)\big{\rangle}$ time dependent state § 2
$\exp_{(\text{o})}$ time-ordered exponential §§ 2, 3
$\hat{\rho}(t)$ time-dependent density operator
$=\big{|}\psi(t)\big{\rangle}\\!\big{\langle}\psi(t)\big{|}$ §§ 2, 3, 6, 6
$p_{ij}$ probability associated with states $\psi_{i}$ and $\psi_{j}$ § 2
$\hat{\mathcal{U}}(T)$ time propagator over the interval $[0,T]$ §§ 2, 7
$\hat{\rho}_{\textsc{1,2}}$ $2$ non-interacting, uncorrelated subsystems
$=\hat{\rho}_{\textsc{2}}\\!\otimes\\!\hat{\rho}_{\textsc{2}}$ § 2
$\hat{\rho}_{1,2,\dots k}$ $k$ non-interacting, uncorrelated subsystems
$=\hat{\rho}_{\textsc{2}}\\!\otimes\\!\hat{\rho}_{\textsc{2}}\\!\otimes\\!\cdots\\!\otimes\\!\hat{\rho}_{\textsc{k}}$
§ 2
$\hat{S}^{(k)}$ spin operator for the $k^{\text{th}}$ spin § 2
$\hat{\mathcal{H}}_{1,2}$ Hamiltonian of $2$ non-interacting subsystems
$=\hat{\mathcal{H}}_{1}\\!\otimes\\!\bm{\mathbbm{1}}_{2}\\!+\\!\bm{\mathbbm{1}}_{1}\\!\otimes\\!\hat{\mathcal{H}}_{2}$
§ 2
$\hat{T}_{\ell m}^{(12)}$ $2$–particle tensor operators §§ 2, 4
$\mathcal{C}_{m_{1}m_{2}m}^{\ell_{1}\ell_{2}\ell}$ Clebsch-Gordan coefficients
§ 2
$\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.2535pt}}\smash{\hat{\mathcal{H}}}}}{\hat{\vphantom{\rule{1.0pt}{6.2535pt}}\smash{\hat{\mathcal{H}}}}}{\hat{\vphantom{\rule{1.0pt}{6.74754pt}}\smash{\hat{\mathcal{H}}}}}{\hat{\vphantom{\rule{1.0pt}{7.0866pt}}\smash{\hat{\mathcal{H}}}}}$
Hamiltonian superoperator
$=\bm{\mathbbm{1}}\\!\otimes\\!\hat{\mathcal{H}}\\!-\\!\hat{\mathcal{H}}^{\dagger}\\!\otimes\\!\bm{\mathbbm{1}}$
§§ 2, 3, 7
$\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.2535pt}}\smash{\hat{\mathcal{R}}}}}{\hat{\vphantom{\rule{1.0pt}{6.2535pt}}\smash{\hat{\mathcal{R}}}}}{\hat{\vphantom{\rule{1.0pt}{6.74754pt}}\smash{\hat{\mathcal{R}}}}}{\hat{\vphantom{\rule{1.0pt}{7.0866pt}}\smash{\hat{\mathcal{R}}}}}$
relaxation superoperator §§ 2, 3
$\big{|}\rho(t)\big{\rangle}$ time-dependent, vectorised density matrix §§
2, 3, 7
$\mathrm{B_{0}}$ static magnetic field in the $z$-direction §§ 2, 6
$\mathrm{B_{1}}$ pulsed magnetic field §§ 2, 6
$\vec{B}$ external magnetic field vector § 2
$\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{6.66975pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{7.0088pt}}\smash{\vec{S}}}}$
vector of Cartesian spin operators
$=\Big{(}\hat{S}_{x},\hat{S}_{y},\hat{S}_{z}\Big{)}$ §§ 2, 4
$\hat{\mathcal{H}}_{\\!\mathcal{Z}}$ Zeeman interaction Hamiltonian
$=\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{6.66975pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{7.0088pt}}\smash{\vec{S}}}}\\!\cdot\\!\bm{A}\\!\cdot\\!\vec{B}$
§ 2
$\bm{A}$ Zeeman interaction tensor § 2
$\hat{\mathcal{H}}_{\mathcal{CS}}$ chemical shielding Hamiltonian
$=\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{6.66975pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{7.0088pt}}\smash{\vec{S}}}}\\!\cdot\\!\bm{\delta}\\!\cdot\\!\vec{B}$
§ 2
$\bm{\delta}$ chemical shielding tensor § 2
$\hat{\mathcal{H}}_{\mathcal{G}}$ electron Zeeman Hamiltonian
$=\mu_{B}\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.66975pt}}\smash{\vec{L}}}}{\hat{\vphantom{\rule{1.0pt}{7.0088pt}}\smash{\vec{L}}}}\\!\cdot\\!\bm{g}\\!\cdot\\!\vec{B}$
§ 2
$\mu_{B}$ Bohr magneton $=\frac{e\hbar}{2m_{e}}$ § 2
$m_{e}$ electron mass (kg) $\approx 9.11\times 10^{−31}$ § 2
$\bm{g}$ electron g-tensor § 2
$\hat{\mathcal{H}}_{\\!\mathcal{N\\!N}}$ inter-nuclear interaction Hamiltonian
§ 2
$\hat{\mathcal{H}}_{\mathcal{DD}}$ dipole-dipole interaction Hamiltonian § 2
$\gamma$ magnetogyric ratio § 2
$\hat{\mathcal{H}}_{\mathcal{SR}}$ spin-rotation coupling Hamiltonian
$=-\sum\limits_{k}\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.66975pt}}\smash{\vec{L}}}}{\hat{\vphantom{\rule{1.0pt}{7.0088pt}}\smash{\vec{L}}}}\\!\cdot\\!\bm{A}_{k}\\!\cdot\\!\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{6.66975pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{7.0088pt}}\smash{\vec{S}}}}_{k}$
§ 2
$\bm{A}_{k}$ spin-rotation coupling tensors § 2
$\hat{\mathcal{H}}_{\mathcal{Q}}$ quadrupolar coupling Hamiltonian §§ 2, 4
$\hat{\mathcal{H}}_{\mathcal{J}}$ J-coupling Hamiltonian
$a\big{(}\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.66975pt}}\smash{\vec{L}}}}{\hat{\vphantom{\rule{1.0pt}{7.0088pt}}\smash{\vec{L}}}}\\!\cdot\\!\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{6.66975pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{7.0088pt}}\smash{\vec{S}}}}\big{)}$
§§ 2, 4
$\hat{\mathcal{H}}_{\\!\mathcal{N\\!E}}$ electron-nuclear interaction
Hamiltonian § 2
$\hat{\mathcal{H}}_{\mathcal{EE}}$ inter-electron interaction Hamiltonian §
2
$H_{n}$ Hermitian matrices of periodic time-dependent functions § 7
$s(t)$ measured signal as a function of time § 6
$\vec{s}$ measured signal over time § 6
$\vec{r}$ reference signal measured from an experiment with hard pulses § 6
$Q$ characteristic exponents § 7
$\hat{F}_{n}$ ladder operators in Fourier space § 7
$\hat{N}$ number operator § 7

Numerical Optimisation:
$\vphantom{+}{}^{\vphantom{|}}_{\vphantom{|}}J(\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}\bm{c})$
scalar objective function §§ 3, 8
$\vphantom{+}{}^{\vphantom{|}}_{\vphantom{|}}\bm{c}$ objective variable
$=\Big{(}c_{1},c_{2},\dots,c_{n}\Big{)}^{\\!\\!T}$ §§ 3, 8
$\mathbb{R}^{n}$ $n$D space of real numbers §§ 3, 4, 5
$\mathbb{C}^{n}$ $n$D space of complex numbers § 5
$\min\limits_{\vphantom{+}{}^{\vphantom{|}}_{\vphantom{|}}\bm{c}}\Big{\\{}\\!\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}J(\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}\bm{c})\Big{\\}}$
minimum of a function
$\vphantom{+}{}^{\vphantom{|}}_{\vphantom{|}}J(\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}\bm{c})$
§§ 3, 4, 5
$\max\limits_{\vphantom{+}{}^{\vphantom{|}}_{\vphantom{|}}\bm{c}}\Big{\\{}\\!\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}J(\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}\bm{c})\Big{\\}}$
maximum of a function
$\vphantom{+}{}^{\vphantom{|}}_{\vphantom{|}}J(\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}\bm{c})$
§§ 3, 5, 6
$\bm{h}$ small displacement vector § 3
$\vphantom{+}{}^{\vphantom{|}}_{\vphantom{|}}\nabla\\!\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}J(\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}\bm{c})$
gradient vector of objective §§ 3, 4
$\vphantom{+}{}^{\vphantom{|}}_{\vphantom{|}}\nabla^{2}\\!\vphantom{+}{}^{\vphantom{|}}_{\vphantom{|}}J(\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}\bm{c})$
Hessian matrix of objective §§ 3, 4
$\nabla$ del operator $=\Big{(}\frac{\partial}{\partial
c_{1}},\frac{\partial}{\partial c_{2}},\dots,\frac{\partial}{\partial
c_{n}}\Big{)}^{\\!\\!T}$ §§ 3, 4, 5, 6
$\nabla^{2}$ Hessian operator
$=\big{|}\nabla\big{\rangle}\\!\big{\langle}\nabla\big{|}$ §§ 3, 4, 5, 6
$\nabla^{r}$ $r^{\text{th}}$ order derivative §§ 3, 5
$\jmath$ optimisation iteration §§ 3, 4
$\bm{c}^{\vphantom{|}}_{\jmath}$ objective function at an iteration $\jmath$
§ 3
$\bm{c}^{\vphantom{|}}_{\ast}$ minimiser/maximiser of objective function § 3
$\big{|}\bm{c}^{\vphantom{|}}_{\jmath}\big{\rangle}$ direction at a point of
surface, $\bm{c}^{\vphantom{|}}_{\jmath}$ § 3
$\big{\|}\cdot\big{\|}_{n}$ $n$-norm §§ 4, 5, 6
$\big{\|}\cdot\big{\|}_{\infty}$ infinity norm §§ 3, 4
$\big{\|}\cdot\big{\|}_{\text{F}}$ Frobenius norm § 4
$\lambda$ eigenvalues §§ 3, 4
$\vphantom{+}{}^{\vphantom{|}}_{\vphantom{|}}\bm{d}$ line search direction §
3
$\bm{d}^{\vphantom{|}}_{\jmath}$ line search direction at iteration $\jmath$
§ 3
$\varepsilon$ line search step length §§ 3, 4
$\vphantom{+}{}^{\vphantom{|}}_{\vphantom{|}}\bm{\delta}$ Newton step § 3
$\mathcal{O}(\cdot)$ error term § 3
eps double-precision floating point accuracy $\approx 2.22\times 10^{-16}$
§§ 3, 4
$\vphantom{+}{}^{\vphantom{|}}_{\vphantom{|}}\Delta\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}\bm{c}$
change in objective variable
$=\bm{c}^{\vphantom{|}}_{\jmath+1}-\bm{c}^{\vphantom{|}}_{\jmath}$ § 3
$\vphantom{+}{}^{\vphantom{|}}_{\vphantom{|}}\Delta\mathbf{g}$ change in
gradient
$=\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}\nabla\\!\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}J(\bm{c}^{\vphantom{|}}_{\jmath+1})\\!-\\!\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}\nabla\\!\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}J(\bm{c}^{\vphantom{|}}_{\jmath})$
§ 3
$\mathbf{H}_{\vphantom{|}}$ approximation to the Hessian matrix § 3
$\mathbf{I}_{\vphantom{|}}$ approximation to the inverse of the Hessian matrix
§ 3
$\vphantom{+}{}^{\vphantom{|}}_{\vphantom{|}}\Delta\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}\bm{C}_{\\!\jmath}$
store of most recent
$\vphantom{+}{}^{\vphantom{|}}_{\vphantom{|}}\Delta\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}\bm{c}$
§ 3
$\vphantom{+}{}^{\vphantom{|}}_{\vphantom{|}}\Delta\mathbf{G}_{\\!\jmath}$
store of most recent
$\vphantom{+}{}^{\vphantom{|}}_{\vphantom{|}}\Delta\mathbf{g}$ § 3
$\mathbf{L}$ lower triangular matrix § 4
$\mathbf{\Lambda}$ diagonal matrix of eigenvalues § 4
$\mathbf{Q}$ matrix with columns of eigenvectors § 4
$\mathbf{S}$ scaling matrix § 4
$D_{\mathbf{B}}(f(\mathbf{A}))$ derivative of the function $f$ at $\mathbf{A}$
in the direction $\mathbf{B}$ § 5
$D_{\mathbf{B}}^{r}(f(\mathbf{A}))$ $r^{\text{th}}$-derivative of the function
$f$ at $\mathbf{A}$ in the direction $\mathbf{B}$ § 5
$\beta$ simplex expansion factor § 8
$\gamma$ simplex contraction factor § 8
$\bm{c}^{\vphantom{|}}_{j}$ simplex vertex § 8
$\bm{c}^{\vphantom{|}}_{0}$ initial simplex vertex § 8
$\bm{c}^{\vphantom{|}}_{c}$ simplex centroid § 8
$\mathbf{X}$ simplex § 8
$\mathbf{R}$ reflected simplex § 8
$\mathbf{E}$ expanded simplex § 8
$\mathbf{C}$ contracted simplex § 8
$\hat{\mathbf{\Delta}}$ Finite difference matrix/differentiation approximation
§ 6

Optimal Control:
$\hat{\mathcal{H}}_{k}$ Control operator §§ 2, 3
$\hat{\mathcal{H}}_{x}^{(k)},\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.2535pt}}\smash{\hat{\mathcal{H}}}}}{\hat{\vphantom{\rule{1.0pt}{6.2535pt}}\smash{\hat{\mathcal{H}}}}}{\hat{\vphantom{\rule{1.0pt}{6.74754pt}}\smash{\hat{\mathcal{H}}}}}{\hat{\vphantom{\rule{1.0pt}{7.0866pt}}\smash{\hat{\mathcal{H}}}}}_{x}^{(k)}$
$k^{\text{th}}$ control operator in the $x$-direction §§ 2, 3, 4, 6
$\hat{\mathcal{H}}_{y}^{(k)},\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.2535pt}}\smash{\hat{\mathcal{H}}}}}{\hat{\vphantom{\rule{1.0pt}{6.2535pt}}\smash{\hat{\mathcal{H}}}}}{\hat{\vphantom{\rule{1.0pt}{6.74754pt}}\smash{\hat{\mathcal{H}}}}}{\hat{\vphantom{\rule{1.0pt}{7.0866pt}}\smash{\hat{\mathcal{H}}}}}_{y}^{(k)}$
$k^{\text{th}}$ control operator in the $y$-direction §§ 2, 3, 4, 6
${c}_{\scriptscriptstyle k}^{\vphantom{|}}(t)$ time-dependent control
amplitudes §§ 3, 6
$K$ number of control operators §§ 3, 4, 5
$\big{|}\rho_{0}\big{\rangle}$ initial state of system § 3
$\big{|}\sigma\big{\rangle}$ desired state of system §§ 3, 6
$\Delta t$ small time slice §§ 3, 5
$N$ number of time slices §§ 3, 4, 5, 6
$t_{n}$ $n^{\text{th}}$ time interval § 3
$\bm{c}$ 2D vector of $K\times N$ control amplitudes §§ 3, 4, 6
$\bm{c}_{k}^{\vphantom{|}}$ row vector of $N$ control amplitudes §§ 3, 6
$\bm{c}_{n}^{\vphantom{|}}$ column vector of $K$ control amplitudes § 3
${c}_{\scriptscriptstyle k,n}^{\vphantom{|}}$ individual element of $\bm{c}$
§§ 3, 4, 5
$\dot{\rho}(t)$ time derivative of $\rho$ § 3
$\mathcal{L}\big{(}{c}(t)\big{)}$ Lagrangian as a function of control ${c}(t)$
§ 3
$\mathcal{J}{}\big{(}\rho(T)\big{)}$ fidelity: terminal cost of Lagrangian
§§ 3, 4, 5, 6
$\mathcal{K}$ penalty: running cost of Lagrangian §§ 3, 6
$\big{|}\chi(t)\big{\rangle}$ costate/adjoint state of system
$\big{|}\chi(T)\big{\rangle}=\big{|}\sigma\big{\rangle}$ § 3
$\big{|}\rho_{n}\big{\rangle}$ state at a time interval $n$ §§ 3, 4, 5
$\big{|}\chi_{n}\big{\rangle}$ adjoint state at a time interval $n$ §§ 3, 4,
5
$\mathcal{H}_{\textsc{p}}$ pseudo-Hamiltonian § 3
$\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.2535pt}}\smash{\hat{\mathcal{P}}}}}{\hat{\vphantom{\rule{1.0pt}{6.2535pt}}\smash{\hat{\mathcal{P}}}}}{\hat{\vphantom{\rule{1.0pt}{6.74754pt}}\smash{\hat{\mathcal{P}}}}}{\hat{\vphantom{\rule{1.0pt}{7.0866pt}}\smash{\hat{\mathcal{P}}}}}\\!_{n}(\bm{c}_{n}^{\vphantom{|}})$
time propagator of controls §§ 3, 4
$\tilde{\bm{c}}$ normalised control amplitude vector §§ 3, 6
$p_{c}$ nominal power level of control pulses § 3
$\bm{c}^{(\jmath)}$ control vector at an optimisation iteration $\jmath$ § 4
$\mathbf{D}_{m,n}$ Hessian matrix of controls at a single time point, $K\times
K$ § 5
$\mathbf{P}\\!_{n}$ time propagator matrix § 5
$\mathbf{H}$ Hamiltonian matrix § 5

## Preface

Do not think that what is hard for you to master is humanly impossible; but if
a thing is humanly possible, consider it to be within your reach.

– Marcus Aurelius, Meditations IV

### Resources

The main resource used within this thesis is the Spinach software
library111Available at www.spindynamics.org [HOGBEN11], developed by Ilya
Kuprov, Hannah Hogben, Luke Edwards, Matthew Krzystyniak, Peter J. Hore,
Gareth T.P. Charnock, Dmitry Savostyanov, Sergey Dolgov, Liza Sutarina, and
Zenawi Welderufael.

The majority of numerical simulation in this work was completed with Matlab
2016-2017a/b. All random numbers were seeded with rng(1) from within a Matlab
session, unless otherwise stated.

Esr spectroscopy measurements were performed in the Centre for Advanced esr
(CAESR) in the Department of Chemistry, University of Oxford, using a Bruker
BioSpin EleXSys II E580 spectrometer operating with a SpinJet arbitrary
waveform generator (awg) that is based on a SPDevices SDR14 PCI board, with a
0.625 ns time base. Samples were held at 85 K in a Oxford Instruments CF935O
cryostat under a flow of cold N2 gas, controlled by an Oxford Instruments
Mercury instrument temperature controller. At X-band, the Bruker Biospin
ER4118-MD5-W1 resonator was used, which is a sapphire dielectric resonator of
dimensions 5 mm ID, 10 mm OD and 13 mm height, that was overcoupled for pulsed
measurements to a Q-value of about 200. At Q-band the resonator was a ER5107D2
of typically 80-100 MHz bandwidths.

Photoexcitation of the oop-eseem sample was accomplished with a Continuum
Shurlite Nd:YAG emitting a 7 ns laser pulse at the 3rd harmonic of 1064 nm,
355 nm, at a rate of 10 Hz, attenuated with a $\lambda$/2 plate and finally
depolarized to a 1 mJ pulse energy. The beam was found to match the 5 mm
cryostat window and was not further focused on the sample. Synchronisation of
the laser and oop-eseem measurement involved using the E580 console PatternJet
board user-defined channel as an external trigger to the Stanford Research
DG645 delay generator.

The author acknowledges the use of the iridis High Performance Computing
Facility, and associated support services at the University of Southampton, in
the completion of some parts of this work.

This report was typeset by the author using LaTeX and the TeX mathematical
typesetting system of Donald Knuth222D.E.Knuth, tex and metafont, New
Directions in Typesetting, American Mathematical Society and Digital Press,
Bedford, Massachusetts (1979), and vector graphics were generated with

### Acknowledgements

A number of people have influenced this work, particularly my supervisor Prof.
Ilya Kuprov; without his advice, I would not have had confidence to undertake
this project.

Academics working in the area of optimal control that have had a direct
involvement in this work through collaboration and conversation are Mads Sloth
Vinding, Thomas Schulte-Herbrüggen, Tony Reinsperger, Robert Zeier, David
Tannor, Reuven Eitan, Sean Wright, Jack Saywell. Less direct involvement
through discussion, although sometimes more substantial; Thomas Heydenreich,
Wolfgang Kallies, Sophie Schirmer, Frank Langbein, Martin Koos, Burkhard Luy,
Daniel Wisniewski, Shai Machines, and Tim Freegarde.

Magnetisation-to-singlet state work was an outcome of interesting conversation
with Stuart Elliott, Chiristian Bengs, Gabriele Stevanato, Giuseppe Pileio,
and Malcolm Levitt.

Work on solid state nmr was a collaboration with James Jarvis, Phil
Williamson, Marina Carravetta, and Maria Concistre.

In developing the work using the esr awg, I am most grateful to Will Myers for
his extensive knowledge of all hardware under his command, but mostly for his
patience – I’m sure he saw emptiness behind my eyes while he tried to explain
how to configure the spectrometer to perform basic Hahn echo to me. I don’t
know if he knows that this was the only physical experiment I have done since
my undergraduate days with physics experiments in 2001–2004. In enabling me to
do this work I am also grateful to Christiane Timmel, Arzhang Ardavan, Alice
Bowen, Thomas Prisner, Anton Tcholakov, Gavin Morley, and Gunnar Jeschke for
insightful discussions, and technical support of Bruker through Peter Höfer
and Patrick Carl.

Academics to influence me and allow me to believe in myself are Phil Dawson,
Tobias Brandes, Mike Godfrey, Clive Saunders, Klaus Gernoth (whose
undergraduate course in advanced quantum mechanics and path integrals seduced
me into enjoying difficult, and possibly insoluble, problems). Also, Jason
Ralph, Neil Oxtoby, and Charles Hill for allowing me to work in self imposed
solitude; a mistake I needed to make at that time.

Finally, personal acknowledgement is given to my wife Katarzyna – who helped
me in an uncountable number of ways – my father Martin, and my mother Jane.
These people have become a bedrock for supporting me and accepting my
eccentricities without question. Further acknowledgement is also given to
Martin Goodwin and Emmaleen Goodwin, who helped in proofreading this document.

This thesis is dedicated to my son Ruadhán, born three months before
submission of this Ph.D. thesis, and who I hope will surpass all of my
achievements.
Southampton David L. GoodwinOctober 2017

## Introduction

In fact, the mere act of opening the box will determine the state of the cat,
although in this case there were three determinate states the cat could be in:
these being Alive, Dead, and Bloody Furious.

– Terry Pratchett, Lords and Ladies

The subject of this thesis is optimal control and its aim is advancing
existing methods of optimal control for spin systems of magnetic resonance. A
novel method with improved convergence properties will be presented with a
focus on computational efficiency and avoiding numerical problems. The
computationally expensive task of calculating the exponential of a matrix will
be reviewed and concluded with particularly useful numerical efficiencies
which decrease simulation time. The thesis will be completed by investigating
realistic applications of optimal control solutions, including methods to
overcome experimental limitations.

Many descriptions of optimisation use the analogy of hiking in a mountainous
terrain to highlight the character of numerical optimisation. The goal of the
hiking experience could be to get to the highest peak within a range of
mountains, or descend to the lowest altitude. The first of these goals is
finding the maximum altitude within the landscape, and the second is finding
the minimum; here we optimise the altitude as a function of the coordinates.
Although this analogy is easy to relate to, it is not a good example of
optimisation: In navigating the mountain we have a map and possibly a compass,
more importantly we have a clear view of the mountain and can pick our path to
a visible peak; we know the easiest path to the maximum on the outset.
Furthermore, the mountain range has many peaks – an initial path may only find
a low-lying peak which becomes apparent from the summit.

##### Example of numerical optimisation

A better example is to follow that of Isaac Newton, finding the roots to a
simple polynomial function of a single variable, $f(x)$. Assuming the function
is well behaved, smooth with no discontinuities, the maximum and minimum
should reside in places where the gradient of the function is equal to zero.
As a first step in finding these solutions is to use an algorithm to find
iteratively better estimates to the roots of a function, the places where
$f(x)=0$: This algorithm is known as the Newton-Raphson method (also known as
the Newton method) and is attributed to Isaac Newton [NEWTON1671] and Joseph
Raphson [RAPHSON1690]:

$x_{n+1}=x_{n}-\frac{f(x_{n})}{f^{\prime}(x_{n})}$

where $f(x_{n})$ if the function of a real variable $x_{n}$ at the iteration
$n$.
$f^{\prime}(x_{n})=\tfrac{\operatorname{d\\!}f}{\operatorname{d\\!}x}\big{|}_{x_{n}}$
is the tangent to the function $f(x)$ at $x_{n}$. It is instructive to show
the algorithm with an example111The function used by Newton in his original
example [NEWTON1671] was $f(x)=x^{3}-2x+5$.. For the function
$f(x)=\tfrac{1}{3}x^{3}-x+1$, the tangent to the function can be found
analytically by differentiation.

Figure 1: Finding a minimum of a function using tangents of the gradient from
an initial trial point.  a the function $f(x)=\tfrac{1}{3}x^{3}-x+1$ and its
gradient $x^{2}-1$ with two extremal solutions at $x=\pm 1$.  b First
iteration: An initial guess to the extrema at $x=2$. The tangent to the
gradient shows where to look for a better point, at $x=1.25$.  c Second
iteration: tangent of the gradient at $x=1.25$, finding a better point at
$x=1.025$.  d Iterations shown as the tangents on the functional, converging
closer to a tangential solution with zero gradient.

However, the aim is not to find the roots of the function, the aim is to find
the roots of the derivative of the function: the place where the gradient is
$=0$222Pre-dating Isaac Newton’s and Gottfried Leibniz’ work on calculus, in
the middle of the seventeenth century, de Fermat’s work in analytical geometry
used tangents to show that at the extreme points of various curves, the
gradient becomes zero.. Substituting this into the previous equation we find
an iterative algorithm to find the extrema of a function $f(x)$:

$x_{n+1}=x_{n}-\frac{f^{\prime}(x_{n})}{f^{\prime\prime}(x_{n})}$

where the second derivative, the derivative of the gradient is
$f^{\prime\prime}(x_{n})=\tfrac{\operatorname{d\\!}}{\operatorname{d\\!}x}\tfrac{\operatorname{d\\!}f}{\operatorname{d\\!}x}\big{|}_{x_{n}}=\tfrac{\operatorname{d\\!}^{2}f}{\operatorname{d\\!}x^{2}}\big{|}_{x_{n}}$
evaluated at $x_{n}$. The Newton-Raphson method can be used to find these
roots to the gradient of the function, and the method is shown for the
function $f(x)=\tfrac{1}{3}x^{3}-x+1$ in Fig. 1.

The algorithm converges to an extremum at $x=+1$; within an acceptable level
of accuracy at $4$ iterations and further than the machine precision333this is
floating-point arithmetic, also referred to as eps, is about 16 decimal places
on a 64-bit computer after $5$ iterations. Notice that starting from a guess
$x_{0}<0$ converges to a minimum at $-1$, a starting guess of $x_{0}=0$ will
be undefined at the first iteration, and a starting guess $x_{0}>0$ converges
to a maximum at $+1$. The reason is clearly seen with a graphical
representation of the algorithm in Fig. 1.

##### Example of control theory

Optimal control theory is more subtle. Although optimisation is a part of
optimal control theory, the functional of the optimisation is more than
navigating an $n$-dimensional space. Control theory separates the controllable
and uncontrollable physics of the system, allowing the controllable physics to
become part of the optimisation functional, and the uncontrollable physics of
the system become inherent “drift”.

Figure 2: Diagram of a robot manipulator control task. Upper row is the
control task of taking the manipulator from 1 an initial state to 2 a defined
target state. Central row shows a control strategy, moving each control motor
in turn, 3, 4, then 5. Lower row shows an efficient control strategy, allowing
all control motors to be moved simultaneously where possible, 6, then a final
single control 7. Assuming all motors move at the same speed, the control
strategy of the lower low takes 60% of the time it takes the central row
control task.

A classical example of a control theory is shown in Fig. 2, where the task of
taking a robot manipulator from an initial position to a target position can
be more efficient if control theory is used. The robot manipulator in Fig. 2
has three degrees of freedom, three controllable motors to position the
grippers. The easiest way to control the robot manipulator is to move each
pivot in turn, controlling each motor to the desired final position in turn.
However, if motors are allowed to be controlled simultaneously, the task is
performed quicker. This can be viewed as an application of control theory,
although simplistic, reducing the amount of time the control task takes.

A proper use of control theory would take into account the weight of each
section of the manipulator, the weight of any potential load at the end of the
manipulator, the surrounding environment, and ensuring movement of the
manipulator is smooth with motors started slowly and finishing slowly. This
information, the physics, is easily programmed into a robot manipulator which
performs repetitive tasks.

The task of finding the most efficient, the most optimal, controls to move
from the initial state to the target state is the task of optimal control.

### § 1 Thesis Outline

Chapter 1

An introduction to optimal control, broadly reviewing the literature of the
last fifteen years, where optimal control methods have been applied to the
area of magnetic resonance. The review focuses on the optimal control method
know as grape, but includes a number of other methods also used for optimal
control, some well established and some newer.

Chapter 2

This chapter will give an overview of magnetic resonance theory. It will form
a particularly theoretical chapter that will be referred to in later chapters,
and also serve the purpose of setting out mathematical formalisms that are
used in the software toolbox Spinach – the software used in all numerical
simulation within this thesis.

Chapter 3

A presentation of the grape method in the context of magnetic resonance. The
chapter starts from the mathematical area of numerical optimisation,
progressing through control theory to optimal control theory. The grape method
in this thesis uses a superoperator algebra representation of the Liouville-
von Neumann equation and the mathematical derivations of the chapter are all
in the context of this formalism.

Chapter 4

Work developed by the author for this thesis and published in [GOODWIN16]:
advancing the grape method from the previous best superlinearly convergent
method [FOUQUIERES11] to quadratic convergence. The method presented is shown
to scale linearly with a similar computational complexity to previous methods.
The Newton-grape method, unique work for this thesis, will be presented along
with well established numerical procedures to avoid predicted singularities in
the second-order optimal control method. Convergence analysis will be shown
for simulations of magnetic resonance pulse experiments using this optimal
control method.

Chapter 5

The matrix exponential is far from a trivial mathematical operator. This
chapter reviews a number series methods from the “nineteen dubious ways to
compute the exponential of a matrix” [MOLER78, MOLER03], giving recommendation
on which method to use for simulation in magnetic resonance systems – the
method used in the software toolbox Spinach. The chapter then goes on to
present work unique to this thesis, published in [GOODWIN15], of analytical
second directional derivative for the grape method. §§ 4, 5 set out the
focus of this thesis – the Newton-grape method of optimal control.

Chapter 6

An overview of popular quadratic penalty methods is presented, optimal control
of phase in an amplitude-phase representation of control pulses, and smoothed
and shaped pulse solutions with application to finding smooth and robust
pulses to generate the singlet state.

Chapter 7

This chapter sets out the application of optimal control to solid-state nmr,
simulating the optimisation algorithm over the crystalline orientations in a
powder sample, and also for magic angle spinning experiments. The chapter
builds to finding optimal control pulses to excite the overtone transition in
an hmqc-type experiment.

Chapter 8

Work on optimal control in the area of electron spin resonance is far from
trivial – the pulses sent through the spectrometer are not the same shape as
those seen by the sample. This is not a problem for shapes that do not vary
rapidly – but optimal control pulses tend to be very specific, with high
frequency components that are very specific to the outcome of the experiment.
To investigate this, a real-time feedback control loop is set up and
experiments are optimised in a brute-force way. Experimental results are shown
for Hahn echo experiments with shaped pulses produced from this type of real-
time feedback loop.

Chapter 9

Summary of main findings of the thesis.

The remainder of this introductory chapter will outline a survey of published
literature in the area of optimal control. Although not extensive in its
applications, the story of the development of optimal control theory in the
context of finding optimal pulses for magnetic resonance systems is fully
considered. The reader is directed to a more detailed perspective of modern
optimal control theory in the report Training Schrödinger’s cat: Quantum
optimal control: Strategic report on current status, visions and goals for
research in Europe [GLASER15].

### § 2 Early optimal control

Optimal control theory in magnetic resonance started at a similar time to
hardware implementation of arbitrarily shaped pulses [SILVER84, TEMPS84].
Although not optimal control theory in itself, numerical optimisation was used
to design highly selective excitation pulses [LURIE85] and population
inversion pulses [WARREN84]. A simple approach was to minimise the square of
the difference between an ideal signal and a measured signal,

$\min_{\text{pulse shape}}\Big{(}\textsc{Ideal Signal}-\textsc{Measured
Signal}\Big{)}^{2}$

using a Nelder-Mead simplex method [NELDER65] varying parameters describing a
shaped pulse [LURIE85]. This approach to optimisation will be revisited in §
8.

The development of appropriate hardware, and the early successes of
experiments using non-rectangular pulses, led researchers to explore the use
of optimal control theory to design “better” pulses for magnetic resonance
experiments. Selective excitation problems using numerically designed
${}^{\pi}\\!/_{\\!2}$ and $\pi$ pulses used optimal control theory to solve
Bloch equations444Bloch equations describe magnetisation of isolated spins.
with a piecewise constant approximation of an arbitrary pulse shape
[CONOLLY86]. The optimisation metric used was defined as a minimum-distance
measure between the current magnetisation vector and a desired magnetisation
vector, then solving a Lagrange equation to find an optimum. Later, a
modification to the minimum-distance measure to include a measure of
adiabaticity gave adiabatic and band selective pulse shape solutions
[ROSENFELD96].

In magnetic resonance systems an ensemble of spins is manipulated to achieve a
desired ensemble state by the application of unitary operations [ERNST87].
Investigation into the controllability of quantum systems showed that, given a
defined initial state and enough pulsing time, any eigenstate of the system
can be populated when
$\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.2535pt}}\smash{\hat{\mathcal{H}}}}}{\hat{\vphantom{\rule{1.0pt}{6.2535pt}}\smash{\hat{\mathcal{H}}}}}{\hat{\vphantom{\rule{1.0pt}{6.74754pt}}\smash{\hat{\mathcal{H}}}}}{\hat{\vphantom{\rule{1.0pt}{7.0866pt}}\smash{\hat{\mathcal{H}}}}}_{x}$
and
$\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.2535pt}}\smash{\hat{\mathcal{H}}}}}{\hat{\vphantom{\rule{1.0pt}{6.2535pt}}\smash{\hat{\mathcal{H}}}}}{\hat{\vphantom{\rule{1.0pt}{6.74754pt}}\smash{\hat{\mathcal{H}}}}}{\hat{\vphantom{\rule{1.0pt}{7.0866pt}}\smash{\hat{\mathcal{H}}}}}_{y}$
controls are available, and their commutators with the drift span the active
space [JUDSON90]555Unitary operations,
$\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.14584pt}}\smash{\hat{\mathcal{H}}}}}{\hat{\vphantom{\rule{1.0pt}{6.14584pt}}\smash{\hat{\mathcal{H}}}}}{\hat{\vphantom{\rule{1.0pt}{6.6948pt}}\smash{\hat{\mathcal{H}}}}}{\hat{\vphantom{\rule{1.0pt}{7.07153pt}}\smash{\hat{\mathcal{H}}}}}_{x}$
controls, and
$\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.14584pt}}\smash{\hat{\mathcal{H}}}}}{\hat{\vphantom{\rule{1.0pt}{6.14584pt}}\smash{\hat{\mathcal{H}}}}}{\hat{\vphantom{\rule{1.0pt}{6.6948pt}}\smash{\hat{\mathcal{H}}}}}{\hat{\vphantom{\rule{1.0pt}{7.07153pt}}\smash{\hat{\mathcal{H}}}}}_{y}$
controls will be introduced in § 2.

### § 3 Analytical optimal control

One of the first successful applications of optimal control theory in nmr used
an existing deuterium decoupling pulse sequence [SCHENKER87] as an initial
guess to a gradient-based numerical optimisation algorithm [LEVANTE96]. The
optimisation metric was defined by the eigenvalues and eigenvectors of the
propagator from an average Hamiltonian. The analytical derivatives of the
linear operators [AIZU63] used for the gradient-based optimisation were
derived for a piecewise constant approximation of an arbitrary pulse sequence.

Using a method developed for laser spectroscopy [JUDSON92, WARREN93],
analytical design of optimal unitary operations for use in nmr proved to be a
difficult task for complicated, coherence-order selective experiments
[UNTIDT98, UNTIDT99]. However, the maximal unitary bounds on general nmr
polarisation transfer experiments were shown to to be much larger than the
apparent limits of nmr polarisation transfer using state-of-the-art
experiments of that time [GLASER98, UNTIDT99]. The search for solutions nearer
to these bounds, finding pulsed solutions to give better experimental results,
and the search for solutions at these bounds, would become the work of the
following decade.

Given that a system is controllable and can be steered from one defined state
to another [JUDSON90], optimal solutions must exist. With the aim of designing
optimally short pulses, to minimise the effects of relaxation, unitary
transforms are found analytically by separating out a Lie subgroup
corresponding to controllable operations, then deducing pulses from the
unitary operator [KHANEJA01]. Pulses with the shortest possible time were
found for coherence transfer experiments in heteronuclear two-spin systems
[KHANEJA01] and linear three-spin systems [KHANEJA02].

### § 4 Broadband optimal pulses

Optimal control pulses can be designed to be broadband, robust to resonance
offsets, by simulating the optimal control algorithm over a range of offsets
[SKINNER03]. These pulses are called Broadband Excitation By Optimised Pulses
(bebop). Further to this, pulses designed by bebop were designed to have a
shorted duration and less power [SKINNER04].

A similar method was used to create Broadband Inversion By Optimised Pulse
(bibop) [KOBZAR04]. This optimal control problem is simulated over a range of
control pulse power levels – creating pulses robust to miscalibration.
Allowing only the phase of the pulses to vary, keeping a constant amplitude,
allows pulses designed with bebop to be calibration-free [SKINNER05,
SKINNER06].

### § 5 Gradient based optimal control

The original paper setting out the theory for Gradient Ascent Pulse
Engineering666also termed Gradient Assisted Pulse Engineering (grape) followed
the work on analytical optimisation [KHANEJA01] \- setting out an elegant base
for Newton-type optimisation methods applied to quantum systems [KHANEJA05].
The essence of the method involved discretising the system into slices of
time, each with a constant Hamiltonian [CONOLLY86], and splitting the system
into an uncontrollable drift Hamiltonian, and control operators with
associated amplitude constants. The validity of the approximation defines the
validity of grape. Along with similar methods, the measure to maximise – to
make the highest possible value – is the overlap between two states termed the
fidelity of the two states [GLASER98]:

$\max_{\text{pulse shape}}\Big{\langle}\textsc{Target State }\Big{|}\textsc{
Current State}\Big{\rangle}$

The angled brackets, $\langle\,\cdot\,|\,\cdot\,\rangle$, are a formal
definition of the overlap between two vectors, denoting the inner
product777The inner product, $\langle\,\cdot\,|\,\cdot\,\rangle$, is a
generalisation of the scalar product between two vectors. In this context,
$|\,\cdot\,\rangle$ is a row vector and $\langle\,\cdot\,|$ is a column
vector, using the Dirac bra-ket notation [DIRAC30]. of the desired Target
state vector and the simulated Current state of the system.

Following its publication in 2005, grape immediately found applications in
creating broadband pulses designed with bebop, as any optimal control method
should, to systems requiring robust excitation pulses – pulses that are
tolerant to miscalibration of the nominal pulse power. The grape method of
optimal control can easily be modified to find broadband pulse excitation with
tolerance to radio-frequency ($\mathrm{B_{1}}$) miscalibration/inhomogeneity
[KOBZAR05]. Although the system and optimal control task are simple,
transferring longitudinal magnetisation of an uncoupled spin to transverse
magnetisation, the example shows the possible flexibility of grape, where
rectangular pulses alone would fail to be adequately robust.

Transferring a single initial to final magnetisation component, point-to-point
state transfer, may not entirely describe an adequate desired final state of
the system [KOBZAR04, KOBZAR08]. Simple point-to-point888also termed state-to-
state problems; these define an initial state and a target state each as
single, or linear combination of, members of the basis set optimal control
solutions can be used within the framework of finding universal rotation
pulses, required for mixing or refocusing pulses [LUY05]. Universal rotation
pulses are solutions to problems where each of the three magnetisation vectors
of the Bloch equation must have a defined rotation [MACHNES11,
SCHULTE_HERBRUEGGEN11]. When these universal rotations are the target state of
a grape simulation, the resulting pulse shapes are found to exhibit symmetry
around the centre of the pulse [KOBZAR12]. Designing universal rotation pulses
to be robust is called Broadband Universal Rotation By Optimised Pulses
(burbop) [JANICH11, JANICH12, SKINNER12a].

The cpmg sequence (Carr-Purcell-Meiboom-Gill [CARR54, MEIBOOM58]), useful in
investigation of relaxation rates999Relaxation is called decoherence in some
literature, was redesigned with grape resulting in universal rotations. This
application extended the cpmg robustness to a larger range of static
$\mathrm{B_{0}}$ resonance offsets and $\mathrm{B_{1}}$ pulse power
miscalibration [BORNEMAN10]. Robust shaped pulses, intended to be direct
replacements for hard ${}^{\pi}\\!/_{\\!2}$ and $\pi$ pulses, were designed
using the grape method [NIMBALKAR13]. Furthermore, broadband optimal solutions
for heteronuclear spin systems can be designed to compensate for coupling
evolution over the duration of a pulse shape, with the strategy of using a
point-to-point optimisation with its time reversal for finding universal
rotation pulses [EHNI13].

Spin systems simulated with grape speed up calculation of coherence transfer
in larger spin and pseudospin chains [SCHULTE_HERBRUEGGEN05]. Further to this,
a method of using grape with partitioned subsystems used to describe a large
spin system proves effective in finding optimal control solutions for very
large spin systems [RYAN08]. The long-range state transfer can be improved
with the temporary occupation of higher order states for larger systems
[KHANI09].

Pulses produced by grape may be difficult to implement on real hardware.
Hardware constraints may interpolate a shape to a smooth realisation of
pulses, rather than the square pulses produced by the grape method. Smooth,
shaped pulses can be produced as analytic functions with a modification to the
fidelity functional [SKINNER10]. Alternatively, this hardware constraint can
be accounted for in simulations by increasing the discretisation of time
evolution above that of the discrete control pulses [MOTZOI11], using
experimentally obtained transfer matrices within grape [BORNEMAN12], or by
employing a feedback control loop after a designed pulse sequence is
implemented [EGGER14].

#### § 1 Solid-state nuclear magnetic resonance

When nmr is performed on a solid powder, the anisotropic spin interactions
result in large line broadening of the spectrum. Spinning the sample at the
magic angle101010the magic angle is
$\theta_{\text{magic}}=\arccos{\tfrac{1}{\oldsqrt[\ ]{3}}}=54.7^{\circ}$
relative to the static $\mathrm{B_{0}}$ axis averages out many of these
interactions, resulting in sharper spectral line shapes [ERNST87]. This method
is termed mas (Magic Angle Spinning), and can be useful for investigating
dipolar and quadrupolar systems [FRYDMAN95].

In the same year as the publication of the grape method, it found application
to solid-state mas nmr, accounting for anisotropic components of spin
interactions, sample spinning and specifically addressing problems of
$\mathrm{B_{1}}$ inhomogeneity [VOSEGAARD05] within simulation. Design of
pulses robust to hardware miscalibration of $\mathrm{B_{1}}$ averages the
grape method over a range of $\mathrm{B_{1}}$ power levels, in effect,
performing a grape simulation at each power level, then averaging the
fidelities and the gradients [VOSEGAARD05]. A range of chemical shift offsets
can be treated in a similar way to $\mathrm{B_{1}}$ miscalibration
[SKINNER06]. Chemical shift offsets, $\mathrm{B_{1}}$ inhomogeneity, and
dispersion from powder averaging, can also be simulated with the grape
framework [TOSNER06], in this case using effective Hamiltonians. A combination
of optimal control pulses and average Hamiltonian theory can be used to
analyse and then remove unwanted second and third order coupling terms
produced by optimal control solutions [BJERRING13].

A method of solid-state nmr can take advantage of the Overhauser effect
[OVERHAUSER53] to transfer polarisation from unpaired electrons to nearby
nuclei, in a method termed dnp (Dynamic Nuclear Polarisation). The framework
to control dnp systems [GRIESINGER12] is set out in [KHANEJA07b, MAXIMOV08],
and simple dnp systems can be handled with grape to give optimal solutions
[HODGES08] in the presence of relaxation [POMPLUN08]. Optimal control proves
useful in combining incoherent and coherent transfer schemes, investigated in
a system of one electron and two nuclear spins [POMPLUN10].

Again in the context of solid-state nmr, optimal pulses found use in
experiments with labelled proteins [KEHLET07, HANSEN07]. Biomolecular samples,
where high power control pulses can overheat the sample, dictate the need for
low power pulse solutions. These can be designed with grape by penalising high
power pulse solutions [KEHLET05, NIELSEN09], and can be extended to
$\mathrm{{}^{2}H-^{13}C}$ cross-polarisation mas experiments [WEI11]. The
quadratic power penalty is also used in calculating optimal pulses for
conversion of parahydrogen induced longitudinal two-spin order to single spin
order [BRETSCHNEIDER12], and in designing robust pulses for transfer of
magnetisation to singlet order [LAUSTSEN14]. Further radio frequency power and
amplitude restriction are investigated for broadband excitation pulses
[KOOS15].

Methods of simulating solid-state mas to find optimal control pulses will be
given in § 7 and an overview of useful penalty methods used in conjunction
with grape will be given § 6.

#### § 2 Nuclear magnetic resonance imaging

A useful application of nmr in medicine is non-invasive, macroscopic imaging
of the human body [ERNST87]. Termed mri (Magnetic Resonance Imaging), the
method uses magnetic field gradients to record a projection of proton density
such that resonant frequencies are a linear functions of spatial coordinates
[LAUTERBUR73]. Extending grape into an average optimisation over an ensemble
of Hamiltonians can compensate for the dispersion in the system dynamics in
mri simulations [LI06, KHANI12].

When radio-frequency pulse durations are similar to relaxation rates of
quadrupolar nuclei, as in the case of mri applications, previous work on
optimal control with mas [VOSEGAARD05] fails to be relevant. This was
investigated, resulting in a novel method of using grape [LEE08]. The idea is
to find optimal solutions with grape for long total pulse duration, then
decreasing the total pulse duration but using the same shape. When the total
pulse duration was reduced to less than the inverse of the quadrupolar
frequency, the resulting pulse shapes were simple, resembling two discrete
pulses separated by a delay. Results from that work showed an increase in
signal intensity of a $\mathrm{spin\\!-\\!{}^{3}\\!/_{\\!2}}$ central
transition while suppressing its satellite transitions. In the context of mri,
this is useful in separating the signal of $\mathrm{{}^{23}Na}$ in cartilage
from that of free $\mathrm{{}^{23}Na}$ within the image. Although found by
optimal control, the solutions to the problem indicated simple pulse shapes
that do not need optimal control [LEE09]. Further investigation proved optimal
control theory found its use in designing robust pulses [LEE10], where simple
pulse shapes could not. Robust pulses were designed with grape to excite the
central transition of arbitrary $\mathrm{spin\\!-\\!{}^{3}\\!/_{\\!2}}$
systems with static powder patterns [ODELL10]. Dft (Density Functional
Theory)111111a popular computational method to investigate electronic
structure. calculations can be used to predict the necessary interaction
parameters for these grape simulations [ODELL11], giving a very general
approach to numerically finding optimal control pulses.

A practicable realisation of an efficient optimal control algorithm for mri
purposes should be fast, working in a small window of time, to be able to take
patient specific distortions into account. The optimisation should be
performed as the patient waits in the instrument, and timing is critical –
should the patient move during the optimal control simulation, the distortions
also change and the optimal control solution becomes ineffective [VINDING12].
Furthermore, advanced arrangements of control coils used in modern
mri121212the arrangement is called pTX (parallel transmission), and can have
an arrangement of up to 6 control coils, rather the the conventional 2 in mri.
need to be taken into account when designing optimal pulses [MASSIRE13]. In
comparison to conventional adiabatic pulses, optimal control pulses can
achieve similar results with reduced sar131313Patient safety guidelines have
strict limitations on sar (specific-absorption-rate) – local hot-spots of
temperature, due to the distribution of sar when using pTX, will be different
for every patient and every scan. [LEE16] and within local and global sar
limits [VINDING17a, VINDING17b].

#### § 3 Quantum information

The research area of quantum information processing is based using two-level
quantum systems, and can be thought of as extending the language of classical
computation to these quantum bits [NIELSEN00, WISEMAN10]. Just as the language
of classical computation is based on circuitry and logic gates, so too is the
area of quantum information processing. The difference between classical bits
and quantum bits is that quantum bits can exist in a superposition of the
classical binary states.

Although optimal control was applied to the area of quantum information before
the advent of grape, and there have been many applications of optimal control
to this area, this review is intended to show applications of grape to
magnetic resonance systems – specifically to ensembles of spins141414Magnetic
resonance on single spins is an emerging area, but beyond the scope of this
work.. As such, only a modest number of applications of grape to this area
will be mentioned; those the author judges general enough to have use in
optimal control of an ensemble system in magnetic resonance.

Compared with the previous best methods for control of superconducting two-
level systems, grape can achieve lower information leakage, higher fidelity,
and faster gates [SPORL07]. A strategy for reducing information leakage
entails using an extra control channel, proportional to the derivative of the
main control, within grape [MOTZOI09]. A chain of 3 coupled superconducting
two-level systems can use grape to reduce unwanted “crosstalk” between the two
outer two-level systems, so mediating transfer in a serial manner
[GROSZKOWSKI11]. High fidelity spin entanglement can be achieved with
universal rotations designed with grape, demonstrated experimentally for
nitrogen-vacancy centres in diamond [DOLDE14].

The robust nature of grape has been shown to be tolerant to errors in pulse
length and off-resonance pulses for nitrogen-vacancy centres in diamond
[SAID09, DOLDE14], atomic vapour ensembles [SCHONFELDT09], and design of fault
tolerant quantum gates using a modified grape optimisation update rule
[NIGMATULLIN09, SCHIRMER09a]. Time optimal solutions [KHANEJA01] giving the
minimal pulse length to achieve high fidelity transfer of Ising coupled
cluster states are verified with grape simulated at different pulse lengths
[FISHER09], with exponential falloff below a threshold pulse length [KHANI09].
Realistic pulse implementation can be accounted for in the form of power and
rise-time limits [FISHER10]. Optimal grape solutions for these realistic
pulses are realised by averaging optimisation over individual detunings of
resonance frequencies (with longitudinal magnetisation control), with optimal
pulses using only a single $x$-magnetisation control.

Although the grape method is designed for closed systems, in the case of
coupling to an environment the method was advanced to open non-Markovian
[REBENTROST09], and open Markovian systems [SCHULTE_HERBRUEGGEN11]. Control
pulses that maximise transfer starting with a high-temperature thermal
ensemble are automatically optimal for lower-temperature initial ensembles
[WANG10]. Further, the control landscape of an open system is similar to a
large closed system [WU12] and there is evidence that noise can be suppressed
with control [ARENZ14].

#### § 4 Electron spin resonance

The particularly challenging application of optimal control to esr (Electron
Spin Resonance)151515also termed epr (Electron Paramagnetic Resonance). Esr is
used in this thesis to avoid conflict with the quantum mechanics term “epr-
pair” derived from the Einstein-Podolsky-Rosen thought experiment [EINSTEIN35]
includes the problem of the pulse response function of the microwave setup –
the pulses seen by the sample differ from those sent to the pulse shaping
hardware. The transfer matrix (response function) may be measured either by
adding an antenna to the resonator [SPINDLER12], or by using a sample with a
narrow esr line to pick up the intensity of each spectral component
[KAUFMANN13]. Quasi-linear responses, such as the phase variation across the
excitation bandwidth in nutation frequency experiments, can be described with
additional transfer matrices [DOLL13]. Once measured, the transfer matrix can
be used within an optimal control algorithm to transform pulse shapes before
trajectory time evolution is calculated, then be used to design broadband esr
excitation pulses [SPINDLER12].

This relatively new area of applying optimal control pulses to esr systems
will be investigated in § 8, finding a complimentary approach to those
mentioned in the previous paragraph.

### § 6 Modified methods

Effective Hamiltonian methods [TOSNER06] and the grape algorithm [KHANEJA05]
have showed great success in finding optimal control solutions. However, more
powerful algorithms have the potential to take full advantage of optimal
control. Further modifications of the grape method spawned two useful optimal
control methods, optimal tracking and cooperative pulses.

#### § 1 Optimal tracking

Optimal tracking is a modified grape method that takes an input control field
that will approximate a desired output trajectory rather than a state
[NEVES09]. Heteronuclear decoupling is an example of a limitation of the
original grape algorithm. The control problem involves taking a system to a
desired intermediate state then back to the original state. In this context,
defining the cost functional as the overlap between the initial and final
states is not enough, since the initial and desired final states are the same
– the system may simply do nothing and produce a maximal overlap while not
penetrating the intermediate state. As an extension of the grape algorithm,
optimal tracking [NEVES09] propagates the system to an intermediate desired
state while simultaneously propagating backward from the final state to this
desired intermediate state – so ensuring that the initial and final states are
always the same. Optimal tracking used for heteronuclear decoupling sequences,
with linear offset-dependent scaling factors, can encode chemical shift
correlations from an indirect dimension [SCHILLING12, ZHANG13].

#### § 2 Cooperative pulses

Cooperative pulses take into account multi-pulse experiments, generalising the
experimental method of phase cycling and cooperatively compensating for each
pulse imperfection [BRAUN10]. By generalising the grape method, it can be
applied to broadband and band-selective experiments [BRAUN10], quantum filter
type experiments [BRAUN14], and pseudo-pure state preparation in the context
of quantum information processing [WEI14].

### § 7 Other optimal control methods

Krotov type methods are algorithms of optimal control [KROTOV96] that were
developed in parallel to the grape method reviewed above. The methods were
successfully developed for application to laser spectroscopy [TANNOR92] and
general two-level quantum systems [ZHU98]. Krotov methods recursively find
pulses necessary to perform a state-transfer problem by solving the Lagrange
multiplier problem [PEIRCE88]. The algorithm was recently formulated in the
context of magnetic resonance [MAXIMOV08] and is an alternative to grape. It
was further developed to a hybrid type method, including a quasi-Newton
character to improve convergence [EITAN11]. The total pulse time of this
method can be critical in finding high fidelity solutions to optimal control
problems, and a characteristic threshold must be passed to find these
solutions using Krotov type methods [CANEVA09]. Krotov based methods can
appear unattractive in that they need a seemingly arbitrary constant in the
form of a Lagrange multiplier – and although the constant is the same for
similar systems, any new system needs investigation to find its own Lagrange
multiplier.

The crab (Chopped Random Basis) algorithm shows promise in finding optimal
control solutions to many-body problems [DORIA11, CANEVA11]. Specifically, the
method lends itself to tensor-network descriptions of the system, also known
as time-dependent density matrix renormalisation group [SCHOLLWOCK05] or the
tensor train formalism [SAVOSTYANOV14] – an almost approximation-free
description. An extension of the algorithm avoids local extrema by
reformulating pulse bounds [RACH15]. Pulse shapes are described by Fourier
harmonics and these become the variables of optimisation. Calculating the
gradient elements in this formalism was expected to be prohibitively
expensive, therefore the method uses a gradient-free optimisation method such
as the Nelder-Mead simplex method [NELDER65]. A further modification includes
gradient calculations in the goat method (Gradient Optimisation of Analytical
controls) [MACHNES15] giving improved convergence.

A method for efficiently calculating gradients will be presented in § 5.

### § 8 Numerical computation

A number of public software packages include functionality to perform optimal
control calculations for magnetic resonance problems. A selection is listed
below:

Dynamo

[MACHNES11, SCHULTE_HERBRUEGGEN12]
Using Matlab, includes grape, Krotov, and hybrid methods for quantum optimal
control simulations of open and closed systems.

QuTiP

[JOHANSSON12, JOHANSSON13]
Using Python, includes the grape and crab algorithms using l-bfgs and Nelder-
Mead algorithms, respectively.

Simpson

[BAK00]
The release of the magnetic resonance simulation package simpson progressed to
include optimal control in the form of the grape method [TOSNER09, TOSNER14].
Includes the ability to simulate liquid- and solid-state nmr, mri, and quantum
computation. Simpson is coded in C with the Tcl scripting interface.

Spinach

[HOGBEN11]
In the context of optimal control, Spinach solves the Liouville-von Neumann
equation using state-space restriction techniques in an irreducible spherical
tensor basis using Matlab. Grape-l-bfgs [FOUQUIERES11] and Krotov methods were
included before the author became involved in the project161616The Spinach
software package is used and developed by the author of this document..

The first systematic look at the algorithmic complexity of the grape algorithm
[KHANEJA05] was published in [GRADL06]. Efficient implementation would take
advantage of inherent parallelisation of the grape algorithm. The bottleneck
in computation was identified as performing a matrix exponential – a number of
approaches to numerically compute this were compared: scaling and squaring
based on the Padé approximation and series expansion with either the Taylor
expansion or Chebyshev polynomials [SCHULTE_HERBRUEGGEN09, AUCKENTHALER10].

The obvious way to create a gradient for a gradient ascent algorithm is to
define it by a central finite-difference approximation – however, for anything
but small systems, this gives a computationally expensive calculation171717A
two point finite difference requires two functional evaluations, a four point
requires four etc. whereas an analytic gradient would require one functional
evaluation by definition. – making optimal control of larger systems
impractical. The need for analytic gradients is paramount for any gradient
ascent method to become useful in larger systems [KUPROV09]. Based on an
auxiliary matrix formalism, analytic gradient elements can be calculated by
taking the exponential of a $2\times 2$ triangular block matrix, with the
drift Hamiltonian on the diagonal and the control operator on the off-diagonal
[FLOETHER12].

A Quasi-Newton optimisation of grape, using the $\ell$-bfgs algorithm, was
developed to make full use of the gradient history to give second-order
derivative information and improve convergence to superlinear [FOUQUIERES11].
A Newton-Raphson root finding approach, using the conditioned Jacobian matrix
of first-order partial derivatives as a least-squares optimisation problem
[FOUQUIERES12], improved convergence to quadratic for a chain of five coupled
two-level systems using two controls. The Jacobian matrix is also used in
Krylov-Newton methods for closed spin systems, showing impressive convergence
without calculation of a full Hessian [CIARAMELLA15].

Many examples in the literature above, when simulating high dimensional spin
systems, produce optimal control solutions which can be described to appear
like “noise”. Trajectory analysis of the “noisy” pulse sequence evolution
shows an underlying order to the spin dynamics [KUPROV13], and a time-
frequency representation of the pulse shapes reveals underlying patterns which
are not apparent in a conventional Cartesian representation of pulses
[KOCHER14].

### § 9 Summary

Optimal control has become a useful tool in exploring the limits of, among
other technologies, magnetic resonance spectroscopy and imaging. The design of
pulsed magnetic resonance experiments with either composite pulses, shaped
pulses, or adiabatic pulses has shown improvement when designed with optimal
control methods.

The grape method of optimal control has become one of the leaders in its field
for finding computationally efficient solutions, able to find state-of-the-art
pulse sequences for large and complicated spin systems. Furthermore, pulse
sequences produced can be made robust to experimental variability.

In addition to becoming an indispensable tool of chemical spectroscopy, the
emerging area of optimal control of mri shows promise in reducing medical
costs at the same time as increasing diagnostic abilities. And it is this last
point that becomes the motivation of the work presented in this thesis –
finding more accurate optimal control solutions with decreased computational
complexity – with the vision of large, highly parallel computational
infrastructure installed at medical facilities, real-time, bespoke optimal
control solutions will become affordable. With reduced experiment time on the
expensive apparatus, and increased throughput: the dream of scientific
facilitators working in business and finance will become more real.

## Magnetic Resonance Theory

A youth who had begun to read geometry with Euclid, when he had learnt the
first proposition, inquired, “What do I get by learning these things?” So
Euclid called a slave and said “Give him three pence, since he must make a
gain out of what he learns.”

– Joannes Stobaeus, Physical and Moral Extracts

Optimal control can be thought of as the study of controlling a system to find
an optimally desired solution. A formal definition of optimal control will be
given in § 3, culminating in the grape method. § 1 was a review of
published literature on the subject of optimal control applied to magnetic
resonance systems. Normally, when studying magnetic resonance systems, an
ensemble of particle spins are manipulated to give molecular information of
how these spins interact with one another.

Magnetic resonance spectroscopy classically uses rectangular pulses of
electromagnetic radiation to manipulate spins located within a static magnetic
field. Non-rectangular shaped pulses are also used in magnetic resonance with
a review given in [FREEMAN98]. The idea of applying optimal control to
magnetic resonance is to optimise desired experimental outcomes with respect
to the pulse shape.

This chapter will introduce features of magnetic resonance theory that will be
used within this thesis. Much of this chapter is taken from standard textbooks
on quantum mechanics [LANDAU58, MERZBACHER61, ARFKEN99, TANNOR07], nuclear
magnetic resonance (nmr) [ABRAGAM61, ERNST87, HORE95] and electron spin
resonance (esr) [SCHWEIGER01, BRUSTOLON09].

### § 1 Angular momentum

Angular momentum in quantum mechanics can be derived from many starting
points. Here, the starting point is the Schrödinger equation and Euclidean
geometry.

#### § 1 Schrödinger’s equation

In 1926 Erwin Schrödinger formulated the famous Schrödinger equation
[SCHRODINGER26] in an attempt to find solutions to a Hamiltonian equation
using Louis de Broglie’s wave-particle duality [DEBROGLIE25]. The Hamiltonian,
$\hat{\mathcal{H}}$, is an energy operator which dictates dynamics in the
equation of motion of that system in the form of a Hamiltonian equation. If
$|\Psi(t)\rangle$ describes the state of a system at a particular time $t$,
then the time-dependent Schrödinger equation is

$\frac{\partial}{\partial
t}\big{|}\Psi(t)\big{\rangle}=-i\hat{\mathcal{H}}\big{|}\Psi(t)\big{\rangle}$
(1)

It is common to measure energy eigenvalues in angular frequency units within
the area of magnetic resonance, and the usual factor of $\hbar$ is therefore
dropped from Eq. 1. The state $|\Psi(t)\rangle$ can be expanded to a complete
orthonormal basis

$\displaystyle\big{|}\Psi(t)\big{\rangle}=\sum\limits_{i=1}^{n}c_{i}(t)\big{|}\Psi_{i}\big{\rangle},$
$\displaystyle i=1,2,\dots,n$ (2)

where the time dependence is now described with the coefficients $c_{i}(t)$ of
the complete orthonormal basis states $\Psi_{i}$. This basis spans a Hilbert
space of dimension $n$.

#### § 2 Euclidean principle of relativity

Quantum mechanics is built on the principle that space is subject to the laws
of Euclidean geometry – and that space itself is isotropic and homogeneous.
This assumption is called the Euclidean principle of relativity. A physical
interpretation is – in the absence of external fields and perturbations the
effect of moving a physical system, such as experimental apparatus, should not
result in measurements that depend on the location or orientation of the
experiment. Or more precisely; a translation or rotation of the coordinate
system should not change the Hamiltonian describing the energy of particles
within the system. The Hamiltonian is invariant under rotation and
translation.

For a rotation operation denoted by $\hat{R}$ the mapping of a state vector
$\Psi$ to $\Psi^{\prime}$ is

$\displaystyle\big{|}\Psi\big{\rangle}\mapsto\hat{R}\big{|}\Psi\big{\rangle}=\big{|}\Psi^{\prime}\big{\rangle}$
$\displaystyle\Rightarrow$
$\displaystyle\big{\langle}\Psi\big{|}\hat{R}^{\dagger}\hat{\mathcal{H}}\hat{R}\big{|}\Psi\big{\rangle}=\big{\langle}\Psi\big{|}\hat{\mathcal{H}}\big{|}\Psi\big{\rangle}$
(3)

This relation must hold for any state vectors, and therefore $\hat{R}$ must
commute with the system Hamiltonian $\hat{\mathcal{H}}$

$\displaystyle\hat{R}^{\dagger}\hat{\mathcal{H}}\hat{R}=\hat{\mathcal{H}}$
$\displaystyle\Rightarrow$
$\displaystyle\hat{R}^{-1}\\!\hat{\mathcal{H}}\hat{R}=\hat{\mathcal{H}}$
$\displaystyle\Rightarrow$
$\displaystyle\big{[}\\!\hat{\mathcal{H}},\hat{R}\big{]}=0$ (4)

Using the Schrödinger Eq. 1, and the measurement of an observable111An
Hermitian operator that possesses a complete set of eigenfunctions is an
observable. property of the system being its expected value222The expected
value of an hermitian operator is
$\big{\langle}\Psi\big{|}\hat{A}\big{|}\Psi\big{\rangle}$ – the following
relation can be derived

$-i\Big{\langle}\Psi\Big{|}\big{[}\hat{\mathcal{H}},\hat{R}\big{]}\Big{|}\Psi\Big{\rangle}=0$
(5)

showing the conservation of the observable property associated with the
rotation $\hat{R}$.

#### § 3 Orbital angular momentum

The property conserved for the above rotation, from the assumption of
isotropic space, is angular momentum . The quantum mechanical angular momentum
operator, $\hat{L}$, is analogous to its classical vector counterpart
[LANDAU60]

$\displaystyle\hat{L}=\hat{r}\times\hat{p}=\hat{r}\times\frac{\hbar}{i}\nabla,$
$\displaystyle\nabla\triangleq\begin{pmatrix}\frac{\partial}{\partial x}\\\
\frac{\partial}{\partial y}\\\ \frac{\partial}{\partial z}\\\ \end{pmatrix}$
(6)

The Cartesian components of the angular momentum operator can be derived by
considering infinitesimal rotations giving

$\displaystyle i\hat{L}_{x}=\bigg{(}y\frac{\partial}{\partial
z}-z\frac{\partial}{\partial y}\bigg{)},$ $\displaystyle
i\hat{L}_{y}=\bigg{(}z\frac{\partial}{\partial x}-x\frac{\partial}{\partial
z}\bigg{)},$ $\displaystyle i\hat{L}_{z}=\bigg{(}x\frac{\partial}{\partial
y}-y\frac{\partial}{\partial x}\bigg{)}$ (7)

where the constant unit of $\hbar$ is dropped, as is normal in the subject of
magnetic resonance. The solution to a differential equation for the
infinitesimal rotation $\varphi$ is

$\displaystyle\frac{\operatorname{d\\!}}{\operatorname{d\\!}\varphi}\hat{R}(\varphi)=-i\hat{L}_{z}\hat{R}$
$\displaystyle\Rightarrow$
$\displaystyle\hat{R}(\varphi)=\operatorname{e}^{-i\hat{L}_{z}\varphi}$ (8)

This can be extended to show all rotation operators are exponentials of the
angular momentum operators, and the angular momentum operators are the basis
of the Lie algebra corresponding to the Lie group of rotations. By inspection
of Eq. 7, the following commutation relations can be formed:

$\displaystyle i\hat{L}_{x}=\big{[}\hat{L}_{y},\hat{L}_{z}\big{]},$
$\displaystyle i\hat{L}_{y}=\big{[}\hat{L}_{z},\hat{L}_{x}\big{]},$
$\displaystyle i\hat{L}_{z}=\big{[}\hat{L}_{x},\hat{L}_{y}\big{]}$ (9)

where the commutator of two elements $a$ and $b$ defines a Lie algebra, with
$[a,b]=ab-ba$. The square of the angular momentum, the modulus, can be
constructed from the angular momentum operator components:

$\hat{L}^{2}=\hat{L}_{x}^{2}+\hat{L}_{y}^{2}+\hat{L}_{z}^{2}$ (10)

and this operator, called the Casimir operator, does commute with each of the
angular momentum operators

$\displaystyle\big{[}\hat{L}^{2},\hat{L}_{x}\big{]}=0,$
$\displaystyle\big{[}\hat{L}^{2},\hat{L}_{y}\big{]}=0,$
$\displaystyle\big{[}\hat{L}^{2},\hat{L}_{z}\big{]}=0$ (11)

A convenient set of non-Hermitian operators, named the raising and lowering
operators, can be defined in terms of $\hat{L}_{x}$ and $\hat{L}_{y}$

$\displaystyle\hat{L}_{+}=\hat{L}_{x}+i\hat{L}_{y},$
$\displaystyle\hat{L}_{-}=\hat{L}_{x}-i\hat{L}_{y}$ (12)

with the following commutation relations involving $\hat{L}_{z}$

$\displaystyle\big{[}\hat{L}_{+},\hat{L}_{-}\big{]}=2\hat{L}_{z},$
$\displaystyle\big{[}\hat{L}_{z},\hat{L}_{+}\big{]}=\hat{L}_{+},$
$\displaystyle\big{[}\hat{L}_{z},\hat{L}_{-}\big{]}=-\hat{L}_{-}$ (13)

#### § 4 The uncertainty principle

Given an Hermitian operator $\hat{L}$, representing an observable property of
the system in a state $\big{|}\Psi\big{\rangle}$; the uncertainty in a
measurement of that observable property can be represented as a standard
deviation [WEYL27]:

$\Delta\hat{L}=\oldsqrt[\
]{\big{\langle}\hat{L}^{2}\big{\rangle}-\big{\langle}\hat{L}\big{\rangle}^{2}}=\oldsqrt[\
]{\big{\langle}\psi\big{|}\hat{L}^{2}\big{|}\psi\big{\rangle}-\big{\langle}\psi\big{|}\hat{L}\big{|}\psi\big{\rangle}^{2}}$
(14)

Heisenberg’s uncertainty principle [HEISENBERG27] states that for two
hermitian operators:

$\displaystyle\Delta\hat{A}\Delta\hat{B}\geqslant\frac{\hbar}{2}\bigg{|}\Big{\langle}\big{[}\hat{A},\hat{B}\big{]}\Big{\rangle}\bigg{|}$
$\displaystyle\big{[}\hat{A},\hat{B}\big{]}=ic$ (15)

where $c$ is a constant. This can be seen with angular momentum operators
applied to the Heisenberg uncertainty relation in atomic units, where
$\hbar=1$:

$\displaystyle\Delta\hat{L}_{x}=$ $\displaystyle\oldsqrt[\
]{\big{\langle}\hat{L}_{x}^{2}\big{\rangle}-\big{\langle}\hat{L}_{x}\big{\rangle}^{2}}$
$\displaystyle\Delta\hat{L}_{y}=$ $\displaystyle\oldsqrt[\
]{\big{\langle}\hat{L}_{y}^{2}\big{\rangle}-\big{\langle}\hat{L}_{y}\big{\rangle}^{2}}$
$\displaystyle\Rightarrow\quad\Delta\hat{L}_{x}\Delta\hat{L}_{y}\geqslant$
$\displaystyle\frac{1}{2}\bigg{|}\Big{\langle}\big{[}\hat{L}_{x},\hat{L}_{y}\big{]}\Big{\rangle}\bigg{|}=\frac{1}{2}\Big{|}\big{\langle}\hat{L}_{z}\big{\rangle}\Big{|}$
(16)

The angular momentum operator components do no commute with each other. This
indicates that all components of angular momentum cannot be simultaneously
defined, unless all equal zero – which sets it apart from linear
momentum333The result obtained by considering infinitesimal translations and
the homogeneity of space, would give relations for the linear momentum
operator.. However, the modulus of the angular momentum, the Casimir operator,
can have a definite value at the same time as each component of the angular
momentum and is therefore simultaneously observable.

#### § 5 Spherical polar coordinates

Figure 1: Diagram of a spherical coordinate system in relation to a Cartesian
coordinate system.

The three Cartesian components of the angular momentum operators can be
expressed in a spherical polar coordinate system, more natural than a
Cartesian coordinate system when considering rotations from § 2. Depicted in
Fig. 1, their formulation is:

$\displaystyle\hat{L}_{x}=-i\bigg{(}-\sin{\varphi}\frac{\partial}{\partial\theta}-\cos{\varphi}\cot{\theta}\frac{\partial}{\partial\varphi}\bigg{)}$
(17)
$\displaystyle\hat{L}_{y}=-i\bigg{(}\cos{\varphi}\frac{\partial}{\partial\theta}-\sin{\varphi}\cot{\theta}\frac{\partial}{\partial\varphi}\bigg{)}$
(18) $\displaystyle\hat{L}_{z}=-i\frac{\partial}{\partial\varphi}$ (19)

Similarly, the transform of the Casimir operator to spherical polar
coordinates gives

$\hat{L}^{2}=-\frac{1}{\sin^{2}{\theta}}\frac{\partial^{2}}{\partial\varphi^{2}}-\frac{1}{\sin{\theta}}\frac{\partial}{\partial\theta}\bigg{(}\sin{\theta}\frac{\partial}{\partial\theta}\bigg{)}$
(20)

The simultaneous eigenvalues of $\hat{L}^{2}$ and $\hat{L}_{z}$ are

$\displaystyle\begin{rcases}\hat{L}^{2}\big{|}Y_{\\!\ell}^{m^{\vphantom{\prime}}}\\!(\theta,\varphi)\big{\rangle}=\ell(\ell+1)\big{|}Y_{\\!\ell}^{m^{\vphantom{\prime}}}\\!(\theta,\varphi)\big{\rangle}\\\
\hat{L}_{z}\big{|}Y_{\\!\ell}^{m^{\vphantom{\prime}}}\\!(\theta,\varphi)\big{\rangle}=m\big{|}Y_{\\!\ell}^{m^{\vphantom{\prime}}}\\!(\theta,\varphi)\big{\rangle}\end{rcases}$
$\displaystyle\begin{array}[]{l}\ell\in\mathbb{N}\\\
m=-\ell,-\ell+1,\dotsc,\ell\end{array}$ (23)

where the normalised, simultaneous eigenfunctions of the Casimir operator,
$\hat{L}^{2}$, and the $z-$projection operator, $\hat{L}_{z}$, are the
spherical harmonics

$\displaystyle
Y_{\\!\ell}^{m^{\vphantom{\prime}}}\\!(\theta,\varphi)=(-1)^{\\!\frac{m+|m|}{2}}\\!\oldsqrt[\
]{\frac{2\ell+1}{4\pi}\frac{\big{(}\ell-\big{|}m\big{|}\big{)}!}{\big{(}\ell+\big{|}m\big{|}\big{)}!}}P_{\\!\ell}^{m^{\vphantom{\prime}}}\\!\big{(}\cos{\theta}\big{)}\operatorname{e}^{im\varphi}$
(24) $\displaystyle
P_{\\!\ell}^{m^{\vphantom{\prime}}}\\!(\xi)=\big{(}1-\xi^{2}\big{)}^{\\!\frac{|m|}{2}}\frac{\operatorname{d\\!}^{\>|m|}}{\operatorname{d\\!}\xi^{|m|}}P_{\\!\ell}(\xi)$
(25)

where $P_{\\!\ell}(\xi)$ are the Legendre polynomials which can be defined by
the formula

$P_{\\!\ell}(\xi)=\frac{1}{2^{\ell}\ell!}\frac{{\operatorname{d\\!}}^{\>\ell}}{\operatorname{d\\!}\xi^{\ell}}\Big{[}\big{(}\xi^{2}-1\big{)}^{\\!\ell}\Big{]}$
(26)

$\ell$ is the total momentum quantum number, $m$ is the projection quantum
number, $\ell\in\mathbb{N},m=-\ell,-\ell\\!+\\!1,\dotsc,\ell$.

Following Eqs. 17, 18 the transform of raising and lowering operators to
spherical polar coordinates are

$\displaystyle\hat{L}_{+}=$
$\displaystyle\operatorname{e}^{i\varphi}\big{(}\frac{\partial}{\partial\theta}+i\cot{\theta}\frac{\partial}{\partial\varphi}\big{)}$
$\displaystyle\Rightarrow\hat{L}_{+}\big{|}Y_{\\!\ell}^{m^{\vphantom{\prime}}}\big{\rangle}=$
$\displaystyle\oldsqrt[\
]{\ell(\ell+1)-m(m+1)}\big{|}Y_{\\!\ell}^{m+1^{\vphantom{\prime}}}\big{\rangle}$
(27) $\displaystyle\hat{L}_{-}=$
$\displaystyle-\operatorname{e}^{i\varphi}\big{(}\frac{\partial}{\partial\theta}-i\cot{\theta}\frac{\partial}{\partial\varphi}\big{)}$
$\displaystyle\Rightarrow\hat{L}_{-}\big{|}Y_{\\!\ell}^{m^{\vphantom{\prime}}}\big{\rangle}=$
$\displaystyle\oldsqrt[\
]{\ell(\ell+1)-m(m-1)}\big{|}Y_{\\!\ell}^{m-1^{\vphantom{\prime}}}\big{\rangle}$
(28)

The raising and lowering operators perform operations that increment and
decrement the magnetic quantum number $m$.

#### § 6 Intrinsic angular momentum

In dealing with the subject of magnetic resonance the idea of intrinsic
angular momentum (as opposed to orbital angular momentum discussed above),
termed spin, must first be introduced [UHLENBECK25, PAULI27]. The most direct
evidence for this spin is the Stern-Gerlach experiment, proposed by Otto Stern
[STERN21]: Silver atoms are sent in a focused beam through an non-uniform
magnetic field, $B_{0}$ [GERLACH22]. Classically, the atoms should give a
continuous trace, between $+\mu$ and $-\mu$, when they emerge from the
magnetic field. However, observation of the emerging atoms shows a number of
distinct traces, rather than the expected continuous trace, spaced equally
between $+\mu$ and $-\mu$. Conventionally, $\mu$ is considered to be the
magnetic moment of the particle. A diagram of the experiment is shown in Fig.
2.

Figure 2: Schematic diagram of the Stern-Gerlach experiment showing a beam of
silver atoms split by an inhomogeneous magnetic field. The classical picture
of this experiment would show a continuous, rather than discrete, trace on the
photographic plate.

This is attributed to the two spin states of atoms within the incident beam –
spin up and spin down. The “up” and “down” could also be termed $+1$ and $-1$,
or $\alpha$ and $\beta$. The nomenclature is only to indicate that there is a
quantised property of the silver atoms, having two states which become
physically separated in the presence of an external magnetic field.

This spatial splitting of the wave functions describing the silver atoms must
be included in addition to, and independent of (commuting with), the spatial
coordinate description of the wave function. The discrete spin variable
describing the splitting in the Stern-Gerlach can take two values444for a
spin-${}^{1}\\!/_{\\!2}$ particle, which is the case for silver atoms, $+1$ or
$-1$.

Spin degrees of freedom can be split from the spatial degrees of freedom,
forming a spinor:

$\big{|}\Psi(\vec{r},s)\big{\rangle}=\sum\limits_{n,k}\big{|}\Psi_{\\!n}(\vec{r}\,)\big{\rangle}\otimes\big{|}\Psi_{\\!k}(s)\big{\rangle}$
(29)

The spin part of the this construct can be described by a pair of orthogonal
vectors projecting onto the two discrete values of the spin,

$\chi=\begin{pmatrix}c_{{}^{+1}\\!/_{\\!2}}\\\
c_{{}^{-1}\\!/_{\\!2}}\end{pmatrix}=c_{{}^{+1}\\!/_{\\!2}}\begin{pmatrix}1\\\
0\end{pmatrix}+c_{{}^{-1}\\!/_{\\!2}}\begin{pmatrix}0\\\
1\end{pmatrix}=c_{{}^{+1}\\!/_{\\!2}}\big{|}\alpha\big{\rangle}+c_{{}^{-1}\\!/_{\\!2}}\big{|}\beta\big{\rangle}$
(30)

which must be normalised, requiring
$\chi^{\dagger}\chi=|c_{{}^{+1}\\!/_{\\!2}}|^{2}+|c_{{}^{-1}\\!/_{\\!2}}|^{2}=1$.
In the case of a spin-${}^{1}\\!/_{\\!2}$ particle, this basis can be used for
describing the two spin eigenstates of the Stern-Gerlach experiment. Without
derivation, we make the assertion that the spin quantum number can have $2s+1$
values, and for a spin-${}^{1}\\!/_{\\!2}$ particle these values are
$s=\pm{}^{1}\\!/_{\\!2}$. This imposes the condition that
$s=0,\frac{1}{2},1,\frac{3}{2},2,\dots$.

In analogy to orbital angular momentum, for intrinsic angular momentum a
matrix representation can be formulated with raising and lowering operators to
move between the $|\alpha\rangle$ and $|\beta\rangle$ states

$\displaystyle\hat{S}_{+}\big{|}\alpha\big{\rangle}=$
$\displaystyle\big{|}\bm{0}\big{\rangle},$
$\displaystyle\hat{S}_{+}\big{|}\beta\big{\rangle}=\big{|}\alpha\big{\rangle},$
$\displaystyle\Rightarrow$ $\displaystyle\hat{S}_{+}=\begin{bmatrix}0&1\\\
0&0\end{bmatrix}$ (31) $\displaystyle\hat{S}_{-}\big{|}\alpha\big{\rangle}=$
$\displaystyle\big{|}\beta\big{\rangle},$
$\displaystyle\hat{S}_{-}\big{|}\beta\big{\rangle}=\big{|}\bm{0}\big{\rangle},$
$\displaystyle\Rightarrow$ $\displaystyle\hat{S}_{-}=\begin{bmatrix}0&0\\\
1&0\end{bmatrix}$ (32)

and with analogy to orbital angular momentum, the raising and lowering
operators can be formulated in terms of $\hat{S}_{x}$ and $\hat{S}_{y}$:

$\displaystyle\hat{S}_{+}=$ $\displaystyle\hat{S}_{x}+i\hat{S}_{y},$
$\displaystyle\hat{S}_{-}=\hat{S}_{x}-i\hat{S}_{y}$ (33)

The Cartesian operators can be constructed from the raising and lowering
operators, and by requiring the eigenvalues of the $z$-component to be the
allowed values of spin $\pm{}^{1}\\!/_{\\!2}$, giving

$\displaystyle\hat{S}_{x}=$ $\displaystyle\frac{\hat{S}_{+}+\hat{S}_{-}}{2},$
$\displaystyle\hat{S}_{y}=\frac{\hat{S}_{+}-\hat{S}_{-}}{2i},$
$\displaystyle\begin{matrix}\hat{S}_{z}\big{|}\alpha\big{\rangle}=+\frac{1}{2}\big{|}\alpha\big{\rangle}\\\
\hat{S}_{z}\big{|}\beta\big{\rangle}=-\frac{1}{2}\big{|}\beta\big{\rangle}\end{matrix}$
$\displaystyle\hat{S}_{x}=$ $\displaystyle\frac{1}{2}\begin{bmatrix}0&1\\\
1&0\end{bmatrix},$ $\displaystyle\hat{S}_{y}=\frac{1}{2}\begin{bmatrix}0&-i\\\
i&0\end{bmatrix},$ $\displaystyle\hat{S}_{z}=\frac{1}{2}\begin{bmatrix}1&0\\\
0&-1\end{bmatrix}$ (34)

These are called the Pauli matrices of a spin-${}^{1}\\!/_{\\!2}$ particle. It
is easy to see that the square of each of these operators is the same and
proportional to the unit matrix, $\bm{\mathbbm{1}}$. Unlike orbital angular
momentum, the square of spin is also proportional to the identity

$\hat{S}^{2}=\hat{S}_{x}^{2}+\hat{S}_{y}^{2}+\hat{S}_{z}^{2}=\frac{1}{4}\bm{\mathbbm{1}}+\frac{1}{4}\bm{\mathbbm{1}}+\frac{1}{4}\bm{\mathbbm{1}}=\frac{3}{4}\bm{\mathbbm{1}}$
(35)

For spin greater than ${}^{1}\\!/_{\\!2}$, the matrix dimension becomes
larger, and similar operators can be constructed e.g. for spin-$1$ particles:

$\displaystyle\hat{S}_{x}=$ $\displaystyle\frac{1}{\oldsqrt[\
]{2}}\begin{bmatrix}0&1&0\\\ 1&0&1\\\ 0&1&0\end{bmatrix},$
$\displaystyle\hat{S}_{y}=\frac{1}{\oldsqrt[\ ]{2}}\begin{bmatrix}0&-i&0\\\
i&0&-i\\\ 0&i&0\end{bmatrix},$
$\displaystyle\hat{S}_{z}=\frac{1}{2}\begin{bmatrix}1&0&0\\\ 0&0&0\\\
0&0&-1\end{bmatrix}$ (36)

The commutation relations follow by inspection:

$\displaystyle\Big{[}\hat{S}_{x},\hat{S}_{y}\Big{]}=i\hat{S}_{z}$
$\displaystyle\Big{[}\hat{S}_{z},\hat{S}_{x}\Big{]}=i\hat{S}_{y}$
$\displaystyle\Big{[}\hat{S}_{y},\hat{S}_{z}\Big{]}=i\hat{S}_{x}$
$\displaystyle\Big{[}\hat{S}^{2},\hat{S}_{\\{x,y,z\\}}\Big{]}=0$ (37)

### § 2 Irreducible representationsof the rotation group

Returning to the Euclidean principle of relativity, space is required to be
isotropic and homogeneous. The consequence is that when a state vector
$\big{|}\Psi\big{\rangle}$ undergoes a rotation operation in space, $\hat{R}$,
all probabilities must be invariant. Mathematically, this means that all inner
products of two rotated states remain invariant

$\displaystyle\big{|}\Psi\big{\rangle}\longmapsto\hat{R}\big{|}\Psi\big{\rangle}=\big{|}\Psi^{\prime}\big{\rangle}$
(38)
$\displaystyle\Big{|}\big{\langle}\Psi^{\prime}\big{|}\Phi^{\prime}\big{\rangle}\Big{|}^{2}=\Big{|}\big{\langle}\Psi\big{|}\Phi\big{\rangle}\Big{|}^{2}$
(39)

for every pair of state vectors, $\Psi$ and $\Phi$. This mapping is called
isometry and must be reversible. A second rotation results only in a phase
change [MERZBACHER61]:

$\big{|}\Psi^{\prime\prime}\big{\rangle}=\operatorname{e}^{i\varphi}\big{|}\Psi^{\prime}\big{\rangle}$
(40)

The transform of a sum of two vectors is equal to the sum of the transforms of
two vectors. It follows that [MERZBACHER61]

$\displaystyle\big{\langle}\Psi_{a}^{\prime\prime}\big{|}\Psi_{b}^{\prime\prime}\big{\rangle}+\Big{(}\big{\langle}\Psi_{a}^{\prime\prime}\big{|}\Psi_{b}^{\prime\prime}\big{\rangle}\Big{)}^{\ast}$
$\displaystyle=\big{\langle}\Psi_{a}\big{|}\Psi_{b}\big{\rangle}+\Big{(}\big{\langle}\Psi_{a}\big{|}\Psi_{b}\big{\rangle}\Big{)}^{\ast}$
(41)
$\displaystyle\begin{split}\mathrm{Re}\big{\langle}\Psi_{a}^{\prime\prime}\big{|}\Psi_{b}^{\prime\prime}\big{\rangle}&=\mathrm{Re}\big{\langle}\Psi_{a}\big{|}\Psi_{b}\big{\rangle}\\\
\mathrm{Im}\big{\langle}\Psi_{a}^{\prime\prime}\big{|}\Psi_{b}^{\prime\prime}\big{\rangle}&=\pm\mathrm{Im}\big{\langle}\Psi_{a}\big{|}\Psi_{b}\big{\rangle}\end{split}$
(42)
$\displaystyle\Rightarrow\Big{(}\lambda\big{|}\Psi\big{\rangle}\Big{)}^{\prime\prime}$
$\displaystyle=\lambda\big{|}\Psi^{\prime\prime}\big{\rangle}$ (43)
$\displaystyle\Big{(}\lambda\big{|}\Psi\big{\rangle}\Big{)}^{\prime\prime}$
$\displaystyle=\lambda^{\ast}\big{|}\Psi^{\prime\prime}\big{\rangle}$ (44)

where ∗ has been used here to denote the complex conjugate. Eqs. 41, 43 are
the properties of a linear operator, and Eq. 43 implies the transform is
unitary. A unitary transformation operator $\hat{\mathcal{U}}$ satisfies
$\hat{\mathcal{U}}^{\dagger}\hat{\mathcal{U}}=\hat{\mathcal{U}}\hat{\mathcal{U}}^{\dagger}=\bm{\mathbbm{1}}$
and gives a symmetry operation555A transformation that leaves the physics
unaltered is called a symmetry operation, and a symmetry operation applied to
all eigenstates of a basis must leave the Schrödinger equation invariant. This
condition is satisfied if the system Hamiltonian, $\hat{\mathcal{H}}$,
commutes with the transform: $[\hat{\mathcal{H}},\hat{\mathcal{U}}_{R}]=0$, an
important property of a rotation. Eq. 44 characterises an antilinear
operator, not considered in generating rotation operations, which is important
in describing a system under time reversal.

Spin does not rotate when a particle is rotated, but anisotropic interactions
(see § 4) do rotate. To represent these rotations, the following complex
Hermitian matrices will become generators of the rotations

$\displaystyle\hat{J}_{1}\triangleq\begin{bmatrix}0&0&0\\\ 0&0&i\\\
0&-i&0\end{bmatrix},$
$\displaystyle\hat{J}_{2}\triangleq\begin{bmatrix}0&0&i\\\ 0&0&0\\\
-i&0&0\end{bmatrix},$
$\displaystyle\hat{J}_{3}\triangleq\begin{bmatrix}0&i&0\\\ -i&0&0\\\
0&0&0\end{bmatrix}$ (45)

A rotation from a state $|\Psi\rangle$ to $|\Psi^{\prime}\rangle$ can be
expressed as a unitary operator,
$|\Psi\rangle=\hat{\mathcal{U}}_{R}|\Psi^{\prime}\rangle$. From the generators
in Eq. 45 and the rotation in Eq. 40, the unitary rotation operators
representing rotations through angles $\\{\alpha,\beta,\gamma\\}$ are

$\displaystyle\operatorname{e}^{i\hat{J}_{1}\alpha}=\begin{bmatrix}1&0&0\\\
0&\cos{\alpha}&-\sin{\alpha}\\\ 0&\sin{\alpha}&\cos{\alpha}\end{bmatrix}$
$\displaystyle\operatorname{e}^{i\hat{J}_{2}\beta}=\begin{bmatrix}\cos{\beta}&0&-\sin{\beta}\\\
0&1&0\\\ \sin{\beta}&0&\cos{\beta}\end{bmatrix}$ (46)
$\displaystyle\operatorname{e}^{i\hat{J}_{3}\gamma}=\begin{bmatrix}\cos{\gamma}&-\sin{\gamma}&0\\\
\sin{\gamma}&\cos{\gamma}&0\\\ 0&0&1\end{bmatrix}$

$\hat{J}$ must satisfy the general commutation relations to give a well-
defined transform:

$\displaystyle\big{[}\hat{J}_{i},\hat{J}_{j}\big{]}=i\varepsilon_{ijk}\hat{J}_{k},$
$\displaystyle\varepsilon_{ijk}=\begin{cases}1&ijk\in\big{\\{}123,231,312\big{\\}}\\\
-1&ijk\in\big{\\{}321,213,132\big{\\}}\\\ 0&\text{all else}\end{cases}$ (47)

where $\hat{J}$ is a rotation operator and $i,j,k$ take the values $1,2,3$ and
represent the Cartesian components $x,y,z$. Further to the unitary operator
representing a rotation in Eq. 46, we also require the same rotation at
further $2\pi$ angles

$\hat{\mathcal{U}}_{R}=\exp{\big{[}-i\hat{J}\big{(}\varphi+2\pi
k\big{)}\big{]}}=\exp{\big{[}-i\hat{J}\varphi\big{]}}\exp{\big{[}-i\hat{J}2\pi
k\big{]}}$ (48)

where $k$ is an arbitrary integer. The additional factor in Eq. 48 has the
effect of multiplying an eigenstate $Y_{\ell}^{m}$ by $(-1)^{2km}$, where $m$
is the magnetic quantum number and the eigenvalue of an operator
$\hat{J}_{z}$.

Any rotation $\hat{R}$ operation in 3D can be described by three Euler
angles666other rotation conventions exist, such are quaternions or directional
cosine matrices, but Euler angles are shown here for convenience., shown in
Fig. 3, which represent three consecutive rotations by angles
$\alpha,\beta,\gamma$ [ARFKEN99],

$\hat{R}(\alpha,\beta,\gamma)=\begin{bmatrix}\cos{\gamma}&-\sin{\gamma}&0\\\
\sin{\gamma}&\cos{\gamma}&0\\\
0&0&1\end{bmatrix}\begin{bmatrix}\cos{\beta}&0&-\sin{\beta}\\\ 0&1&0\\\
\sin{\beta}&0&\cos{\beta}\end{bmatrix}\begin{bmatrix}\cos{\alpha}&-\sin{\alpha}&0\\\
\sin{\alpha}&\cos{\alpha}&0\\\ 0&0&1\end{bmatrix}$ (49)

Figure 3: Euler angles: This rotation convention corresponds to 1 rotation of
angle $\alpha$ around the $z$-axis, 2 rotation of angle $\beta$ around the
$y^{\prime}$-axis, 3 rotation of angle $\gamma$ around the
$z^{\prime\prime}$-axis.

### § 3 Liouville-von Neumann equation

Magnetic resonance is concerned with ensemble averages of spin particles.
Simply describing the basis in a Hilbert space of § 1 is not enough, and a
more descriptive, larger space must be constructed.

Given the initial state of the system, $|\psi(0)\rangle$, the solution to this
Eq. 1 at a future time $T>0$ is

$\big{|}\psi(T)\big{\rangle}=\exp_{(\text{o})}{\\!\\!\Bigg{[}\\!-i\\!\\!\int\limits_{0}^{T}\\!\\!\hat{\mathcal{H}}(t)\operatorname{d\\!}t\Bigg{]}}\big{|}\psi(0)\big{\rangle}$
(50)

where $\exp_{(\text{o})}{\\![\dots]}$ indicates a time ordered exponential
[DYSON49, DYSON49b, ERNST87, TANNOR07].

The density operator is a description of the state of a system ensemble

$\displaystyle\hat{\rho}(t)$
$\displaystyle=\sum\limits_{i,j}p_{ij}\big{|}\psi_{i}(t)\big{\rangle}\\!\big{\langle}\psi_{\\!j}(t)\big{|}$
(51) $\displaystyle=\big{|}\psi(t)\big{\rangle}\\!\big{\langle}\psi(t)\big{|}$
(52)

where $p_{ii}$ are the associated probabilities for each state $\psi_{i}(t)$,
normalised so that $\sum\limits_{i}p_{ii}=1$. The diagonal elements of
$\hat{\rho}(t)$ are the probabilities of being in an eigenstate of the
Hamiltonian, and the off-diagonal elements of $\hat{\rho}(t)$ indicate a
superposition of states [ERNST87].

The equation of motion for the density operator can be found from the time-
dependent Schrödinger equation:

$\displaystyle\frac{\partial}{\partial t}\hat{\rho}(t)=$
$\displaystyle\frac{\partial}{\partial
t}\Big{(}\big{|}\psi(t)\big{\rangle}\\!\big{\langle}\psi(t)\big{|}\Big{)}=\bigg{(}\frac{\partial}{\partial
t}\big{|}\psi(t)\big{\rangle}\bigg{)}\\!\big{|}\psi(t)\big{\rangle}\\!\big{\langle}\psi(t)\big{|}\\!\bigg{(}\frac{\partial}{\partial
t}\big{\langle}\psi(t)\big{|}\bigg{)}$ $\displaystyle=$
$\displaystyle-i\hat{\mathcal{H}}(t)\big{|}\psi(t)\big{\rangle}\\!\big{\langle}\psi(t)\big{|}+i\big{|}\psi(t)\big{\rangle}\\!\big{\langle}\psi(t)\big{|}\hat{\mathcal{H}}(t)$
$\displaystyle=$
$\displaystyle-i\Big{(}\hat{\mathcal{H}}(t)\hat{\rho}(t)-\hat{\rho}(t)\hat{\mathcal{H}}(t)\Big{)}$

Identifying the commutator, the equation of motion can be written as

$\frac{\partial}{\partial
t}\hat{\rho}(t)=-i\Big{[}\hat{\mathcal{H}}(t),\hat{\rho}(t)\Big{]}$ (53)

This is known as the Liouville-von Neumann equation. In a similar manner to
finding the solutions of Eq. 1 with Eq. 50, the solutions to Eq. 53 are

$\hat{\rho}(T)=\hat{\mathcal{U}}(T)\hat{\rho}(0)\hat{\mathcal{U}}^{-1}(T)$
(54)

where $\hat{\mathcal{U}}(T)$ is the time propagator over a time interval
$[0,T]$ defined by

$\hat{\mathcal{U}}(T)=\exp_{(\text{o})}\\!\\!\Bigg{[}\\!-i\\!\\!\int\limits_{0}^{T}\\!\\!\hat{\mathcal{H}}(t)\operatorname{d\\!}t\Bigg{]}$
(55)

#### § 1 Superoperator algebra

The transform in Eq. 54 can be represented as the multiplication of three
$n\times n$ matrices in the following schematic:

$\scalebox{0.75}{\mbox{$\displaystyle\begin{bmatrix}\rho_{11}&\rho_{12}&\cdots&\rho_{1n}\\\
\rho_{21}&\rho_{22}&\cdots&\rho_{2n}\\\
\vdots&\vdots&\ddots\vphantom{\Bigg{|}}&\vdots\\\
\rho_{n1}&\rho_{n2}&\cdots&\rho_{nn}\end{bmatrix}$}}=\scalebox{0.75}{\mbox{$\displaystyle\begin{bmatrix}&\hphantom{\text{\hskip
45.00008pt}}&\\\ &&\\\ &\hat{\mathcal{U}}&\\\ &\textsc{size: }(n\times n)&\\\
&&\end{bmatrix}$}}\scalebox{0.75}{\mbox{$\displaystyle\begin{bmatrix}\sigma_{11}&\sigma_{12}&\cdots&\sigma_{1n}\\\
\sigma_{21}&\sigma_{22}&\cdots&\sigma_{2n}\\\
\vdots&\vdots&\ddots\vphantom{\Bigg{|}}&\vdots\\\
\sigma_{n1}&\sigma_{n2}&\cdots&\sigma_{nn}\end{bmatrix}$}}\scalebox{0.75}{\mbox{$\displaystyle\begin{bmatrix}&\hphantom{\text{\hskip
45.00008pt}}&\\\ &&\\\ &\hat{\mathcal{U}}^{-1}&\\\ &\textsc{size: }(n\times
n)&\\\ &&\end{bmatrix}$}}$ (56)

The same transform in Eq. 54 can also be written as a matrix-vector
multiplication, numerically more appealing than the product of three matrices,
by using the linear transform of density matrix vectorisation:

$\scalebox{0.75}{\mbox{$\displaystyle\begin{pmatrix}\rho_{11}\\\ \rho_{21}\\\
\vdots\\\ \rho_{n1}\\\ \hdashline[1pt/1.5pt]\rho_{12}\\\ \rho_{22}\\\
\vdots\\\ \rho_{n2}\\\ \hdashline[1pt/1.5pt]\\\ \vdots\\\ \\\
\hdashline[1pt/1.5pt]\rho_{1n}\\\ \rho_{2n}\\\ \vdots\\\ \rho_{nn}\\\
\end{pmatrix}$}}=\scalebox{0.75}{\mbox{$\displaystyle\begin{bmatrix}\hphantom{\text{\hskip
225.00038pt}}\\\ \\\ \\\ \\\ \\\ \\\ \\\
\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.2535pt}}\smash{\hat{\mathcal{U}}}}}{\hat{\vphantom{\rule{1.0pt}{6.2535pt}}\smash{\hat{\mathcal{U}}}}}{\hat{\vphantom{\rule{1.0pt}{6.74754pt}}\smash{\hat{\mathcal{U}}}}}{\hat{\vphantom{\rule{1.0pt}{7.0866pt}}\smash{\hat{\mathcal{U}}}}}\\\
\vspace{-0.5em}\\\ \textsc{size: }(n^{2}\times n^{2})\\\ \\\ \\\ \\\ \\\ \\\
\\\
\end{bmatrix}$}}\scalebox{0.75}{\mbox{$\displaystyle\begin{pmatrix}\sigma_{11}\\\
\sigma_{21}\\\ \vdots\\\ \sigma_{n1}\\\ \hdashline[1pt/1.5pt]\sigma_{12}\\\
\sigma_{22}\\\ \vdots\\\ \sigma_{n2}\\\ \hdashline[1pt/1.5pt]\\\ \vdots\\\ \\\
\hdashline[1pt/1.5pt]\sigma_{1n}\\\ \sigma_{2n}\\\ \vdots\\\ \sigma_{nn}\\\
\end{pmatrix}$}}$ (57)

The commutation relation in Eq. 53 can be represented as a matrix in itself,
called a superoperator, and transforms the Liouville-von Neumann equation into
a matrix-vector multiplication

$\frac{\partial}{\partial
t}\hat{\rho}=-i\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.2535pt}}\smash{\hat{\mathcal{H}}}}}{\hat{\vphantom{\rule{1.0pt}{6.2535pt}}\smash{\hat{\mathcal{H}}}}}{\hat{\vphantom{\rule{1.0pt}{6.74754pt}}\smash{\hat{\mathcal{H}}}}}{\hat{\vphantom{\rule{1.0pt}{7.0866pt}}\smash{\hat{\mathcal{H}}}}}\hat{\rho}$
(58)

where the “double-hat” superoperator notation has been introduced, with the
Liouville superoperator acting on another operator defined as a commutation
$\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.2535pt}}\smash{\hat{\mathcal{H}}}}}{\hat{\vphantom{\rule{1.0pt}{6.2535pt}}\smash{\hat{\mathcal{H}}}}}{\hat{\vphantom{\rule{1.0pt}{6.74754pt}}\smash{\hat{\mathcal{H}}}}}{\hat{\vphantom{\rule{1.0pt}{7.0866pt}}\smash{\hat{\mathcal{H}}}}}\hat{\rho}=[\hat{\mathcal{H}},\hat{\rho}]$,
and is represented by a matrix of size $n^{2}\times n^{2}$. Using this
superoperator representation Eq. 54 becomes

$\displaystyle\big{|}{\rho}(T)\big{\rangle}=\exp_{(\text{o})}\\!\\!\Bigg{[}\\!-i\\!\\!\int\limits_{0}^{T}\\!\\!\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.2535pt}}\smash{\hat{\mathcal{H}}}}}{\hat{\vphantom{\rule{1.0pt}{6.2535pt}}\smash{\hat{\mathcal{H}}}}}{\hat{\vphantom{\rule{1.0pt}{6.74754pt}}\smash{\hat{\mathcal{H}}}}}{\hat{\vphantom{\rule{1.0pt}{7.0866pt}}\smash{\hat{\mathcal{H}}}}}(t)\operatorname{d\\!}t\Bigg{]}\big{|}{\rho}(0)\big{\rangle}$
$\displaystyle\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.2535pt}}\smash{\hat{\mathcal{H}}}}}{\hat{\vphantom{\rule{1.0pt}{6.2535pt}}\smash{\hat{\mathcal{H}}}}}{\hat{\vphantom{\rule{1.0pt}{6.74754pt}}\smash{\hat{\mathcal{H}}}}}{\hat{\vphantom{\rule{1.0pt}{7.0866pt}}\smash{\hat{\mathcal{H}}}}}=\bm{\mathbbm{1}}\otimes\hat{\mathcal{H}}-\hat{\mathcal{H}}^{\dagger}\otimes\bm{\mathbbm{1}}+i\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.2535pt}}\smash{\hat{\mathcal{R}}}}}{\hat{\vphantom{\rule{1.0pt}{6.2535pt}}\smash{\hat{\mathcal{R}}}}}{\hat{\vphantom{\rule{1.0pt}{6.74754pt}}\smash{\hat{\mathcal{R}}}}}{\hat{\vphantom{\rule{1.0pt}{7.0866pt}}\smash{\hat{\mathcal{R}}}}}$
(59)

where
$\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.2535pt}}\smash{\hat{\mathcal{R}}}}}{\hat{\vphantom{\rule{1.0pt}{6.2535pt}}\smash{\hat{\mathcal{R}}}}}{\hat{\vphantom{\rule{1.0pt}{6.74754pt}}\smash{\hat{\mathcal{R}}}}}{\hat{\vphantom{\rule{1.0pt}{7.0866pt}}\smash{\hat{\mathcal{R}}}}}$
is the relaxation superoperator, describing damping and dissipations effects,
and $\bm{\mathbbm{1}}$ is the identity matrix777The identity matrix,
$\bm{\mathbbm{1}}$, is a matrix with ones on every diagonal element and zeros
everywhere else.. The Liouville superoperator defined in this way, as a direct
product, allows interacting composite systems described in a numerically
efficient way.

Continuing from the spherical harmonics of § 5, the vector space of Eqs.
23, 27, 28 can be decomposed to $(2\ell+1)$-dimensional subspaces, which are
invariant under rotation. Under a rotation, a state in the spherical polar
coordinate basis transforms as

$\hat{R}\big{|}Y_{\ell}^{m^{\vphantom{\prime}}}\big{\rangle}=\sum\limits_{m^{\prime}=-\ell}^{\ell}\big{|}Y_{\ell}^{m^{\prime}}\big{\rangle}\mathfrak{D}_{m^{\prime}m}^{(\ell)}(\alpha,\beta,\gamma)$
(60)

where the matrix $\mathfrak{D}_{m^{\prime}m}^{(\ell)}(\alpha,\beta,\gamma)$ is
the Wigner rotation matrix (Tab. 1) with the elements described by

$\mathfrak{D}_{m^{\prime}m}^{(\ell)}(\alpha,\beta,\gamma)=\big{\langle}Y_{\ell}^{m^{\prime}}\big{|}\hat{R}\big{|}Y_{\ell}^{m}\big{\rangle}$
(61)

shown explicitly in § A. A rotation in one of these subspaces can be
expressed by an irreducible spherical tensor operator $\hat{T}_{\ell m}$ in a
matrix representation, which is defined as [SANCTUARY76]

$\hat{R}\hat{T}_{\ell
m}\hat{R}^{-1}=\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.2535pt}}\smash{\hat{R}}}}{\hat{\vphantom{\rule{1.0pt}{6.2535pt}}\smash{\hat{R}}}}{\hat{\vphantom{\rule{1.0pt}{6.74754pt}}\smash{\hat{R}}}}{\hat{\vphantom{\rule{1.0pt}{7.0866pt}}\smash{\hat{R}}}}\hat{T}_{\ell
m}=\sum\limits_{m^{\prime}=-\ell}^{\ell}\hat{T}_{\ell
m^{\prime}}\mathfrak{D}_{m^{\prime}m}^{(\ell)}(\alpha,\beta,\gamma)$ (62)

where $\alpha,\beta,\gamma$ are Euler angles in Eq. 49. For the infinitesimal
rotations used in §§ 3, 2, $\hat{J}$, Eq. 62 simplifies to

$\big{[}\hat{J},\hat{T}_{\ell
m}\big{]}=\sum\limits_{m^{\prime}=-\ell}^{\ell}\hat{T}_{\ell
m^{\prime}}\big{\langle}Y_{\ell}^{m^{\prime}}\big{|}\hat{J}\big{|}Y_{\ell}^{m}\big{\rangle}$
(63)

giving an alternate form of Eqs. 23, 27, 28 as commutation relations for the
spherical tensor operator:

$\displaystyle\big{[}\hat{J}_{z},\hat{T}_{\ell m}\big{]}=m\hat{T}_{\ell m}$
(64) $\displaystyle\big{[}\hat{J}_{+},\hat{T}_{\ell m}\big{]}=\oldsqrt[\
]{\ell(\ell+1)-m(m+1)}\,\hat{T}_{\ell m+1}$ (65)
$\displaystyle\big{[}\hat{J}_{-},\hat{T}_{\ell m}\big{]}=\oldsqrt[\
]{\ell(\ell+1)-m(m-1)}\,\hat{T}_{\ell m-1}$ (66)

The single particle spherical tensor operators of rank-1 are

$\displaystyle\hat{T}_{00}=\bm{\mathbbm{1}}$ (67)
$\displaystyle\hat{T}_{10}=\hat{L}_{z}$ (68) $\displaystyle\hat{T}_{1\pm
1}=\mp\frac{1}{\oldsqrt[\ ]{2}}\hat{L}_{\pm}$ (69)

The full list of irreducible spherical tensor operators is given in § A
Tab. 2. Within the software toolbox Spinach [HOGBEN11], all rotations are
active, meaning the system is rotated and not the coordinates.

#### § 2 Composite systems

The total wavefunction of two non-interacting uncorrelated subsystems is a
direct product of the subsystems888Also called the Kronecker product.:

$\big{|}\psi_{\textsc{1,2}}\big{\rangle}=\big{|}\psi_{\textsc{1}}\big{\rangle}\otimes\big{|}\psi_{\textsc{2}}\big{\rangle}$
(70)

and the density matrices are

$\hat{\rho}_{\textsc{1,2}}=\big{|}\psi_{\textsc{1,2}}\big{\rangle}\\!\big{\langle}\psi_{\textsc{1,2}}\big{|}=\big{(}\big{|}\psi_{\textsc{1}}\big{\rangle}\\!\big{\langle}\psi_{\textsc{1}}\big{|}\big{)}\\!\otimes\\!\big{(}\big{|}\psi_{\textsc{2}}\big{\rangle}\\!\big{\langle}\psi_{\textsc{2}}\big{|}\big{)}=\hat{\rho}_{\textsc{2}}\otimes\hat{\rho}_{\textsc{2}}$
(71)

For a system of $k$ non-interacting uncorrelated subsystems, the system is
described by a chain of $k$ direct products:

$\hat{\rho}_{1,2,\dots
k}=\hat{\rho}_{\textsc{2}}\otimes\hat{\rho}_{\textsc{2}}\otimes\cdots\otimes\hat{\rho}_{\textsc{k}}$
(72)

Similarly, when the angular momentum of a composite system is represented by a
direct product space, so to are the normalised simultaneous eigenvectors:

$\big{|}Y_{\ell_{1},\ell_{2}}^{m_{1},m_{2}}\big{\rangle}=\big{|}Y_{\ell_{1}}^{m_{1}}\big{\rangle}\otimes\big{|}Y_{\ell_{2}}^{m_{2}}\big{\rangle}$
(73)

and the state for $n$ particles can be expressed with $n$ direct products

$\big{|}Y_{\ell_{1},\ell_{2},\dots,\ell_{n}}^{m_{1},m_{2},\dots,m_{n}}\big{\rangle}=\big{|}Y_{\ell_{1}}^{m_{1}}\big{\rangle}\otimes\big{|}Y_{\ell_{2}}^{m_{2}}\big{\rangle}\otimes\cdots\otimes\big{|}Y_{\ell_{n}}^{m_{n}}\big{\rangle}$
(74)

which form the basis of this direct product space.

The operators for a multi-particle system are direct products, where the
operator of the $k^{\text{th}}$ subsystem occurs at the $k^{\text{th}}$
position in the direct product chain e.g.

$\displaystyle\hat{S}_{x}^{(k)}=\bm{\mathbbm{1}}\\!\otimes\\!\bm{\mathbbm{1}}\\!\otimes\\!\dotsm\\!\otimes\\!\hat{S}_{x}\\!\otimes\\!\dotsm\\!\otimes\\!\bm{\mathbbm{1}}$
$\displaystyle\hat{S}_{y}^{(k)}=\bm{\mathbbm{1}}\\!\otimes\\!\bm{\mathbbm{1}}\\!\otimes\\!\dotsm\\!\otimes\\!\hat{S}_{y}\\!\otimes\\!\dotsm\\!\otimes\\!\bm{\mathbbm{1}}$
(75)
$\displaystyle\hat{S}_{z}^{(k)}=\bm{\mathbbm{1}}\\!\otimes\\!\bm{\mathbbm{1}}\\!\otimes\\!\dotsm\\!\otimes\\!\hat{S}_{z}\\!\otimes\\!\dotsm\\!\otimes\\!\bm{\mathbbm{1}}$

where $\hat{S}_{\\{x,y,z\\}}^{(k)}$ is the operator representing an $x-$,$y-$,
or $z-$ projection of spin on the $k^{\text{th}}$ particle in the system of
particles.

The Hamiltonians of two non-interacting uncorrelated subsystems must not allow
the subsystems to affect each other, giving a composite Hamiltonian of two
subsystems as

$\hat{\mathcal{H}}_{1,2}=\hat{\mathcal{H}}_{1}\\!\otimes\\!\bm{\mathbbm{1}}_{2}+\bm{\mathbbm{1}}_{1}\\!\otimes\\!\hat{\mathcal{H}}_{2}$
(76)

where $\bm{\mathbbm{1}}$ is the identity operator.

In extending this description of spherical tensor operators to coupled
particles, Clebsch-Gordan coefficients can be used to couple two tensors of a
given rank to another rank [ARFKEN99]. The Clebsch-Gordan coefficients can be
interpreted as a unitary transformation from a linear combination of products
of single particle operators to two particle operators [ERNST87]. Two particle
tensor operators, $\hat{T}_{\ell m}^{(12)}$, from the single particle
operators $\hat{T}_{\ell_{1}m_{1}}^{(1)}$ and $\hat{T}_{\ell_{2}m_{2}}^{(2)}$
defined in this way are

$\hat{T}_{\ell
m}^{(12)}=\sum\limits_{m_{1}}\mathcal{C}_{m_{1}m_{2}m}^{\ell_{1}\ell_{2}\ell}\hat{T}_{\ell_{1}m_{1}}^{(1)}\hat{T}_{\ell_{2}m-m_{1}}^{(2)}$
(77)

where $\mathcal{C}_{m_{1}m_{2}m}^{\ell_{1}\ell_{2}\ell}$ are Clebsch-Gordan
coefficients. The two particle spherical tensor operators of rank-2 are

$\displaystyle\hat{T}_{2\pm 2}^{(2)}=\frac{1}{2}\hat{L}_{\pm}^{2}$ (78)
$\displaystyle\hat{T}_{2\pm
1}^{(2)}=\mp\frac{1}{2}\big{(}\hat{L}_{z}\hat{L}_{\pm}+\hat{L}_{\pm}\hat{L}_{z}\big{)}$
(79) $\displaystyle\hat{T}_{20}^{(2)}=\oldsqrt[\
]{\frac{2}{3}}\big{(}\hat{L}_{z}^{2}-\frac{1}{4}\big{(}\hat{L}_{+}\hat{L}_{-}+\hat{L}_{-}\hat{L}_{+}\big{)}\big{)}$
(80)

The full list of irreducible spherical tensor operators is given in § A
Tab. 3.

### § 4 The Spin Hamiltonian

Spin Hamiltonians are subject to three types of interactions, and can only
have three generic algebraic forms [ERNST87]:

Linear in the spin operators

$\longrightarrow$ Interactions with the static magnetic field,
$\mathrm{B_{0}}$, and the irradiated radio/microwave frequency
$\mathrm{B_{1}}(t)$ field.

Bilinear in the spin operators

$\longrightarrow$ Interactions between two spins.

Quadratic in the spin operators

$\longrightarrow$ Indirect interactions e.g. with field gradient, interpreted
as the interaction of a spin and itself.

In generating the total Hamiltonian of a arbitrary spin system, it is useful
to partition the system into general parts corresponding to different physical
processes

$\hat{\mathcal{H}}=\hat{\mathcal{H}}_{\\!\mathcal{Z}}+\hat{\mathcal{H}}_{\\!\mathcal{N\\!N}}+\hat{\mathcal{H}}_{\\!\mathcal{N\\!E}}+\hat{\mathcal{H}}_{\mathcal{EE}}+\hat{\mathcal{H}}_{k}$
(81)

Each of the split Hamiltonian parts represents the following interactions

$\hat{\mathcal{H}}_{\\!\mathcal{Z}}=$ Zeeman interactions

$\longrightarrow$ Chemical shielding tensors and g-tensors.

$\hat{\mathcal{H}}_{\\!\mathcal{N\\!N}}=$ Inter-Nuclear interactions

$\longrightarrow$ J-couplings, quadrupolar interactions, dipolar interactions.

$\hat{\mathcal{H}}_{\\!\mathcal{N\\!E}}=$ Electron-Nuclear interactions

$\longrightarrow$ Isotropic and anisotropic hyperfine couplings.

$\hat{\mathcal{H}}_{\mathcal{EE}}=$ Inter-Electron interactions

$\longrightarrow$ Exchange interaction, zero field splitting, dipolar
interactions.

$\hat{\mathcal{H}}_{k}=$ Microwave/Radiofrequency irradiation

$\longrightarrow$ Irradiated pulses with amplitude and phase, or control
pulses with associated amplitude coefficients.

In describing these interactions it is useful to introduce the notation for a
vector of spin operators,
$\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{6.66975pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{7.0088pt}}\smash{\vec{S}}}}$,
which contains the three Cartesian operators as elements, ordered
$\big{(}\hat{S}_{x},\hat{S}_{y},\hat{S}_{z}\big{)}$.

#### § 1 Zeeman Interactions

An applied external magnetic field, $\vec{B}$ induces an energy difference
$\Delta E$ between the quantised magnetic moment values, with the energy is
given by

$E=-\vec{\mu}\cdot\vec{B}$ (82)

where $\vec{B}$ is the external magnetic field. For nucleons, this is the
nuclear Zeeman interaction, and for electrons, this is called the electron
Zeeman interaction. Its Hamiltonian is

$\hat{\mathcal{H}}_{\mathcal{Z}}=\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{6.66975pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{7.0088pt}}\smash{\vec{S}}}}\\!\cdot\\!\bm{A}\\!\cdot\\!\vec{B}=\begin{pmatrix}\hat{S}_{x}&\hat{S}_{y}&\hat{S}_{z}\end{pmatrix}\begin{pmatrix}a_{xx}&a_{xy}&a_{xz}\\\
a_{yx}&a_{yy}&a_{yz}\\\
a_{zx}&a_{zy}&a_{zz}\end{pmatrix}\begin{pmatrix}B_{x}\\\ B_{y}\\\
B_{z}\end{pmatrix}$ (83)

where the tensor $\bm{A}$ is the Zeeman interaction tensor. This is the
coupling of the spin dipole moment to the external magnetic field. Without any
further interaction, the Zeeman tensor reduces to the constant
$\frac{g\mu}{\hbar}$, where $g$ is the $g-$factor and $\mu$ is the nuclear
magneton for nuclei and the Bohr magneton for electrons.

##### Chemical shielding tensor

The magnetic field also perturbs the electronic structure of the nucleus, with
the induced field $\bm{\delta}\vec{B}_{0}$. Using perturbation theory
[MERZBACHER61] gives the spin interaction from the contribution of
$\bm{\delta}\vec{B}_{0}$ as

$\hat{\mathcal{H}}_{\mathcal{CS}}=\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{6.66975pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{7.0088pt}}\smash{\vec{S}}}}\\!\cdot\\!\bm{\delta}\\!\cdot\\!\vec{B}=\begin{pmatrix}\hat{S}_{x}&\hat{S}_{y}&\hat{S}_{z}\end{pmatrix}\begin{pmatrix}\delta_{xx}&\delta_{xy}&\delta_{xz}\\\
\delta_{yx}&\delta_{yy}&\delta_{yz}\\\
\delta_{zx}&\delta_{zy}&\delta_{zz}\end{pmatrix}\begin{pmatrix}B_{x}\\\
B_{y}\\\ B_{z}\end{pmatrix}$ (84)

where $\bm{\delta}$ is known as the chemical shielding tensor.

##### Electron g-tensor

A similar interaction to the perturbation from chemical shielding can be
derived for electrons. The energy operator for the electron dipole in an
external magnetic field is

$\hat{\mathcal{H}}_{\mathcal{G}}=\mu_{B}\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.66975pt}}\smash{\vec{L}}}}{\hat{\vphantom{\rule{1.0pt}{7.0088pt}}\smash{\vec{L}}}}\\!\cdot\\!\bm{g}\\!\cdot\\!\vec{B}=\begin{pmatrix}\hat{L}_{x}&\hat{L}_{y}&\hat{L}_{z}\end{pmatrix}\begin{pmatrix}g_{xx}&g_{xy}&g_{xz}\\\
g_{yx}&g_{yy}&g_{yz}\\\
g_{zx}&g_{zy}&g_{zz}\end{pmatrix}\begin{pmatrix}B_{x}\\\ B_{y}\\\
B_{z}\end{pmatrix}$ (85)

where $\bm{g}$ is known as the $g-$tensor and the constant $\mu_{B}$ is the
Bohr magneton.

#### § 2 Inter-Nuclear Interactions

##### Inter-nuclear dipolar interactions

Treating nuclear spins as point magnetic dipoles, $\vec{\mu}$, the energy of
the interaction between two nuclear dipoles, $\vec{\mu}_{1}$ and
$\vec{\mu}_{2}$ is

$E=-\vec{\mu}\\!\cdot\\!\vec{B}=-\frac{\mu_{0}}{4\pi}\frac{3\big{(}\vec{\mu}_{1}\\!\cdot\\!\vec{r}\big{)}\big{(}\vec{r}\\!\cdot\\!\vec{\mu}_{2}\big{)}-r^{3}\big{(}\vec{\mu}_{1}\\!\cdot\\!\vec{\mu}_{2}\big{)}}{r^{5}}$
(86)

where $\vec{r}$ is the vector that connects the dipoles, and the constant
$\mu_{0}$ is the permeability of free space. Hamiltonian is for this
interaction between two spins,
$\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{6.66975pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{7.0088pt}}\smash{\vec{S}}}}^{(1)}$
and
$\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{6.66975pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{7.0088pt}}\smash{\vec{S}}}}^{(2)}$

$\hat{\mathcal{H}}_{\mathcal{DD}}=-\frac{\mu_{0}}{4\pi}\frac{\gamma_{1}\gamma_{2}\hbar}{r^{5}}\Big{(}3\big{(}\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{6.66975pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{7.0088pt}}\smash{\vec{S}}}}^{(1)}\\!\cdot\\!\vec{r}\big{)}\big{(}\vec{r}\\!\cdot\\!\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{6.66975pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{7.0088pt}}\smash{\vec{S}}}}^{(2)}\big{)}-r^{2}\big{(}\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{6.66975pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{7.0088pt}}\smash{\vec{S}}}}^{(1)}\\!\cdot\\!\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{6.66975pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{7.0088pt}}\smash{\vec{S}}}}^{(2)}\big{)}\Big{)}$
(87)

where $\gamma_{1}$ and $\gamma_{2}$ are the magnetogyric ratios of the two
spins.

##### Spin-rotation coupling

Collections of particles have angular momentum, just as single particles do. A
magnetic induction is generated by a collection charges moving at velocity.
Using Eq. 83 and considering the motion of the $n^{\text{th}}$ charge with
mass $m_{n}$ and angular momentum operator
$\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.66975pt}}\smash{\vec{L}}}}{\hat{\vphantom{\rule{1.0pt}{7.0088pt}}\smash{\vec{L}}}}_{n}$
– the spin rotation coupling interaction is

$\hat{\mathcal{H}}_{\mathcal{SR}}=\sum\limits_{k}\gamma_{k}\vec{B}\\!\cdot\\!\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{6.66975pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{7.0088pt}}\smash{\vec{S}}}}_{k}=-\frac{\mu_{0}}{4\pi}\sum\limits_{n,k}\frac{\gamma_{k}q_{n}}{m_{n}r_{n}^{3}}\Big{(}\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.66975pt}}\smash{\vec{L}}}}{\hat{\vphantom{\rule{1.0pt}{7.0088pt}}\smash{\vec{L}}}}_{n}\\!\cdot\\!\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{6.66975pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{7.0088pt}}\smash{\vec{S}}}}_{k}\Big{)}$
(88)

where $\gamma_{k}$ is the magnetogyric ratio of the $k^{\text{th}}$ nucleus
and
$\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{6.66975pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{7.0088pt}}\smash{\vec{S}}}}_{k}$
are their spin operators. Evaluating the sum over the $n$ moving charges, the
spin-rotation Hamiltonian is

$\hat{\mathcal{H}}_{\mathcal{SR}}=-\sum\limits_{k}\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.66975pt}}\smash{\vec{L}}}}{\hat{\vphantom{\rule{1.0pt}{7.0088pt}}\smash{\vec{L}}}}\\!\cdot\\!\bm{A}_{k}\\!\cdot\\!\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{6.1757pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{6.66975pt}}\smash{\vec{S}}}}{\hat{\vphantom{\rule{1.0pt}{7.0088pt}}\smash{\vec{S}}}}_{k}$
(89)

where $\bm{A}_{k}$ are spin-rotation coupling tensors.

##### Quadrupolar coupling

Nuclei with spin$>\mathrm{spin\\!-\\!{}^{1}\\!/_{\\!2}}$ have non-spherical
charge distribution, and the direction of this nuclear quadrupole moment is
fixed, it interacts with the spin itself.