#### § 3 Smooth, robust magnetisation-to-singlet transfer

Figure 8: A plot of fidelity indicating the robustness of the smoothed
optimised m2s pulses from Fig. 6. The robustness is measured as a function of
the nominal power level of the pulse amplitude in Hz on the $y$-axis, and the
$\mathrm{B_{0}}$ resonance offsets in Hz on the $x$-axis.

With the assumption that non-smooth pulses will suffer amplifier non-
linearity, smoothed pulse shapes should perform better experimentally.

As was set out in §§ 4, 5, optimisation with a Newton-grape method converges
to a local maximum. The optimisation strategy is to first gain a quick
solution to the state-transfer problem with a $\ell$-bfgs-grape method
[FOUQUIERES11], robust to $\mathrm{B_{0}}$ inhomogeneity, then to feed this
solution into a quadratically convergent Newton-grape method with a quadratic
smoothing penalty of Eq. 21.

The pulses produced by this two-step optimisation are shown in Fig. 6 and are
evidently smooth, much more smooth than in Fig. 5. As in Fig. 5, the pulse
duration is twice that of the slic pulse for the same reason.

The robustness of the pulse shape from Fig. 6 is shown in Fig. 8, showing
very similar robustness to that of Fig. 7. It is expected that these pulses
perform much better in experiments than both the slic pulse and the non-
smoothed robust optimal pulses of Fig. 7.

The point of using the Newton-grape method is that similar smoothed pulses are
not converged to in the $\ell$-bfgs-grape method – the penalty competes with
the fidelity measure during a line search and a pulse is found that is neither
robust, as designed, nor smooth: finding a solution that is a trade-off
between the two. The Newton-grape method seems to be able to cope with this
trade-off, identifying a place where the maximum fidelity and the minimum
penalty coexist. One reason for this is that the Newton method uses
regularisation and conditioning of § 4, and is able to avoid saddle-points.

## Solid StateNuclear Magnetic Resonance

“I must not fear. Fear is the mind-killer. Fear is the little-death that
brings total obliteration. I will face my fear. I will permit it to pass over
me and through me. And when it has gone past I will turn the inner eye to see
its path. Where the fear has gone there will be nothing. Only I will remain.”

– Frank Herbert, Dune

Solid-state-nmr can be thought of as an ensemble of systems at different
orientations, rotating relative to a laboratory frame. The rotation gives a
periodically time-dependent Hamiltonian. A Lebedev quadrature [LEBEDEV75,
LEBEDEV76] can be used to characterise a surface integral over a 3D sphere,
with weighted values at particular points on the surface (§ 1). This surface
integral can be used to characterise an ensemble in the rotating frame in a
Floquet space, [FLOQUET1883, SHIRLEY65].

This chapter is divided into two sections. The first section is a comparison
of the convergence of the Newton-grape and $\ell$-bfgs-grape methods,
simulating the optimal control problem over a Floquet space describing the
crystalline orientations of a solid-state-nmr sample (§ 1). The success of
the Newton-grape method comes with a large expense and is set out in § 2.
The second outlines the description of the solid-state-nmr sample by a Fokker-
Planck formalism from [EDWARDS13, KUPROV16] (§ 2), for a magic-angle-
spinning experiment [ODELL11b] § 3 – using optimal control to excite the
overtone transition in § 1.

### § 1 Floquet Theory

The different orientation rotating relative to a laboratory frame gives
periodically time-dependent coefficients of the Schrödinger equation [ERNST87,
LEVANTE95]:

$\frac{d}{dt}\hat{\mathcal{U}}(t)=-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}}}}}(t)\hat{\mathcal{U}}(t)$
(1)

suggesting the time evolution operator as a Fourier series expansion of the
Hamiltonian [LESKES10]:

$\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)=\sum_{n}H_{n}e^{in\omega
t}$ (2)

where $H_{n}$ are Hermitian matrices of periodic time-dependent functions. The
solution to Eq. 1 is given in the following form [FLOQUET1883, SHIRLEY65]:

$\displaystyle\hat{\mathcal{U}}(t)=\hat{\mathcal{P}}(t)e^{-iQt}\hat{\mathcal{P}}^{-1}(t),$
$\displaystyle\hat{\mathcal{P}}(t)=\sum_{n}\hat{\mathcal{P}}_{n}e^{in\omega
t}$ (3)

This formulation allows a similar representation to be used when expanding the
density matrix [VINOGRADOV05] from Eq. 54

$\rho(t)=\hat{\mathcal{U}}(t)\rho(0)\hat{\mathcal{U}}^{-1}(t)=\sum_{n}\rho_{n}(t)e^{in\omega
t}$ (4)

It is easy to see that
$\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+\tau)=\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)e^{in\omega\tau}$,
defining the periodic nature of the system. The diagonal elements of the
diagonal matrix, $Q$, are called characteristic exponents [SHIRLEY65].
Subtracting multiples of $\omega$ from the diagonal of $Q$ can be compensated
by multiplying $\hat{\mathcal{P}}(t)$ by $e^{ik\omega t}$.

The direct product of the basis states $\big{\\{}|\rho_{i}\rangle\big{\\}}$
with an infinite set of Fourier states $\big{\\{}|n\rangle\big{\\}}$,
resulting from the Fourier expansion in Eq. 2, give the full basis set of
states:

$\Big{\\{}\big{|}\rho_{i},n\big{\rangle}\Big{\\}}=\Big{\\{}\big{|}\rho_{i}\big{\rangle}\Big{\\}}\otimes\Big{\\{}\big{|}n\big{\rangle}\Big{\\}}$
(5)

In Fourier space, the Hamiltonian and density matrix can be written as
[VINOGRADOV05]

$\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}}}}}_{\textsc{r}}(t)=\sum_{n}H_{n}\hat{F}_{n}e^{in\omega
t}=e^{i\omega\hat{N}t}\bigg{\\{}\sum_{n}H_{n}\hat{F}_{n}\bigg{\\}}e^{-i\omega\hat{N}t}$
(6a)
$\displaystyle\rho_{\textsc{r}}(t)=\sum_{n}\rho_{n}(t)\hat{F}_{n}e^{in\omega
t}=e^{i\omega\hat{N}t}\bigg{\\{}\sum_{n}\rho_{n}(t)\hat{F}_{n}\bigg{\\}}e^{-i\omega\hat{N}t}$
(6b)

where $\hat{F}_{n}$ are ladder operators in the Fourier space; $\langle
v+n|\hat{F}_{n}|v\rangle=\bm{\mathbbm{1}}$, and $\hat{N}$ is the number
operator; $\langle
v|\hat{N}|v^{\prime}\rangle=v\delta_{vv^{\prime}}\bm{\mathbbm{1}}$, and their
commutation is $\big{[}\hat{N},\hat{F}_{n}\big{]}=n\hat{F}_{n}$. The
Liouville-von Neumann equation Eq. 58 in Fourier space is

$\frac{d}{dt}\rho_{\textsc{r}}(t)=-i\Big{[}\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}}}}}_{\textsc{r}}(t),\rho_{\textsc{r}}(t)\Big{]}$
(7)

Using this Fourier representation of the Liouville-von Neumann equation, the
interaction Hamiltonian becomes time independent [VINOGRADOV05]. This time
independent Hamiltonian in Floquet space,
$\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}}}}}_{\textsc{f}}$
is

$\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}}}}}_{\textsc{q}}=\sum_{n}H_{n}\hat{F}_{n}+\omega\hat{N}$
(8)

giving the Floquet representation of the Liouville-von Neumann equation. This
defines the evolution of the density operator;

$\frac{d}{dt}\rho_{\textsc{q}}(t)=-i\Big{[}\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}}}}}_{\textsc{q}},\rho_{\textsc{q}}(t)\Big{]}$
(9)

The Floquet density matrix can be evaluated with propagation over a time $t$
using

$\displaystyle\rho_{\textsc{q}}(t)=U_{\textsc{q}}(t)\rho_{\textsc{q}}(0)U_{\textsc{q}}^{-1}(t),$
$\displaystyle U_{\textsc{q}}(t)=e^{-i\hat{H}_{\textsc{q}}t}$ (10)

Floquet theory has been successfully applied to describe solid-state-nmr
[SCHMIDT92, BAIN01, SCHOLZ10]. It should be noted that Kuprov pointed out in
[GOODWIN15] that the auxiliary matrix exponential formalism in Eq. 21 can be
used to find solutions to rotating frame transformations without the usual
method of nested commutators.

#### § 1 Powder average

Various methods exist to discretise a sphere emulating the crystalline
orientations of a powder average [BAK97, PONTI99, EDEN03]; a grid uniformly
distributed over a unit sphere used in solid-state-nmr simulations:

  1. 1

A spiral circling the $z$-axis of a unit sphere [MOMBOURQUETTE92].

  2. 2

Finite triangular elements partitioning a unit sphere [WANG95] (sophe).

  3. 3

Simulating many subspectra, then averaged to give the powder pattern
[VARNER96] (sums).

  4. 4

Iteratively moving particles on the surface of a sphere under the influence of
a repulsive force from other particles on the surface of a sphere [BAK97]
(repulsion).

  5. 5

The use of Gaussian spherical quadrature [LEBEDEV75, LEBEDEV76] to weight the
distribution over the surface of the sphere [EDEN98]

The following work uses the Lebedev quadrature to weight a distribution over
the spherical grid [EDEN98, EDEN03]; the Euler angles and weights are shown in
Tab. 1.

$\alpha$ | $\beta$ | $w$
---|---|---
0.000000 | 1.570796 | 0.066667
3.141593 | 1.570796 | 0.066667
1.570796 | 1.570796 | 0.066667
-1.570796 | 1.570796 | 0.066667
1.570796 | 0.000000 | 0.066667
1.570796 | 3.141593 | 0.066667
0.785398 | 0.955317 | 0.075000
0.785398 | 2.186276 | 0.075000
-0.785398 | 0.955317 | 0.075000
-0.785398 | 2.186276 | 0.075000
2.356194 | 0.955317 | 0.075000
2.356194 | 2.186276 | 0.075000
-2.356194 | 0.955317 | 0.075000
-2.356194 | 2.186276 | 0.075000
Table 1: Orientations, $\alpha$ and $\beta$ (radians), and weights $w$ of the
Lebedev grid quadrature.

#### § 2 Optimal control of astatic powder average

With the Newton-grape developed in §§ 4, 5, it seems a good place to start
to test the method on a large system. As set out above, an optimisation
problem involving an average over large grid will take many Hessian matrix
evaluations, magnifying computational efficiency.

The system set out to test this is similar to the one in § 2, except that
the system is averaged over crystalline orientations. Optimal control of the
powder average has a fidelity, gradients, and Hessians weighted over the
crystalline orientations – rank 17 Lebedev grid with 110 points.

Figure 1: Comparison of optimisation methods for a 14N quadrupolar overtone
signal in static glycine powder. Two methods are compared: Newton-grape method
using a bracketing and sectioning line search; $\ell$-bfgs-grape using a
bracketing and sectioning line search.

The optimal control problem is to excite 14N from a state $\hat{T}_{1,0}$ to
$\hat{T}_{2,2}$ (§ A) through the nuclear quadrupolar interaction with
appropriate control operators. The optimisation if performed over 400 time
points for total pulse duration of 40 $\mu$s at a power level of 40 kHz.

At a particularly short total pulse duration, the optimal solution is not at a
maximum fidelity $=1$. The comparison of results for the convergence of
$\ell$-bfgs-grape and Newton-grape is shown in Fig. 1111One convergence
profile is shown – but convergence over the averaged grid gives very similar
convergence from random starting pulses. It should be noted that the averaging
has dictated a smooth convergence profile.. Newton-grape is quicker in the
initial stages of the optimisation, as is expected compared to a quasi-Newton
method, and has a similar convergence characteristic to the universal rotation
pulses over a power miscalibration grid in § 3. The terminal fidelity is
also slightly better than that of the quasi-Newton method, again similar to
the results of § 3.

However, Newton-grape took much longer (wall clock time) to reach a fidelity
only a few percent better than the quasi-Newton approach – a couple of hours
for $\ell$-bfgs-grape to reach 100 iterations, and half a day for Newton-grape
to reach 100 iterations222Both methods used a 32 core computer parallelised
calculations over the time steps..

#### § 3 Optimal control undermagic angle spinning

Using Floquet space to describe the states of the solid-state-nmr system with
a powder average, then to spin the system at the magic angle333the magic angle
is $\theta_{\text{magic}}=\arccos{\tfrac{1}{\oldsqrt[\ ]{3}}}=54.7^{\circ}$,
is the test of optimal control set out in this section.

The powder average should be modelled by an appropriate Lebedev grid
describing the crystalline orientations in the power. The local fidelity and
the local gradient of the fidelity at each orientation of the grid to produce
a global value of the fidelity and its gradient – the global measure is a
weighted average over the crystalline grid. The local optimal control problem
is performed with the $\ell$-bfgs-grape method.

Using the parameters set out in [JARVIS13] for the overtone excitation in a
glycine molecule, the optimal control problem here is a cross-polarisation
experiment between ${}^{13})C$ and the quadrupolar 14N. A spinning rate of
frequency $f=25\times 10^{3}$ Hz is used to simulate the magic angle spinning
experiment of [JARVIS13], and the Fourier rank at which the Floquet matrix is
truncated is 5 and the crystalline grid is averaged over a Lebedev grid rank 5
$(\alpha,\beta,\gamma)$ relative to the rotor frame (as in Tab. 1).

The initial state and final states of the system are the direct products with
the identity matrix of size $(2\times\text{max-rank})+1$, where $\text{max-
rank}=5$ in this case. This gives the optimal control boundaries as

$\displaystyle\rho(0)=\left(\bm{\mathbbm{1}}\otimes L_{z}^{\text{(C)}}\right)$
(11a) $\displaystyle\sigma{}=\left(\bm{\mathbbm{1}}\otimes
L_{z}^{\text{(N)}}\right)$ (11b)

where both states should be normalised. Control operators are in the Floquet
formalism irradiating the resonant frequencies of 13C and 14N, again as a
direct product with the identity matrix:

$\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}}}}}_{1}=\bm{\mathbbm{1}}\otimes\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}^{\text{(C)}},$
$\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}}}}}_{2}=\bm{\mathbbm{1}}\otimes\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}^{\text{(C)}},$
$\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}}}}}_{3}=\bm{\mathbbm{1}}\otimes\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}^{\text{(N)}},$
$\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}}}}}_{4}=\bm{\mathbbm{1}}\otimes\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}^{\text{(N)}}$
(12)

with the nominal power of 10 kHz. Time steps are integer multiples of the
characteristic time from the spinning rate, 5000 in this case, and the total
pulse duration is 20 ms.

The $\ell$-bfgs-grape method is used; following from the long wall clock time
of the Newton-grape, the limited memory quasi-Newton method is more
appropriate considering the large number of optimisation variables, $4\times
5000$. The optimisation was run for 256 iterations, although the optimisation
had not fully converged at this point, a good set of results were obtained
without waiting for super-linear convergence to find small gains.

Figure 2: The four channel control pulse sequences producing the trajectories
in Fig. 3. The pulses are a result of the $\ell$-bfgs-grape optimisation,
which reached high fidelity, state-to-state transfer from 13C to 14N of
glycine in a solid-state-nmr simulation.

The control pulses produced at the end of the 256 iterations are shown in
Fig. 2 – looking like noise, apart from the 14N pulses generally needing less
power.

Figure 3: Trajectories analysis for a state-to-state transfer from 13C to 14N
of glycine in a solid-state-nmr simulation. Averaged over the crystalline
orientations using a Lebedev quadrature. 10 kHz power level of control pulses
and 14.1 MHz magnet.  a Local spin population.  b Spin correlation order
amplitude.

The trajectory analysis [KUPROV13] of the pulses shows more useful
information. As intended, the local population of the magnetisation is
transferred from 13C to 14N in Fig. 3 a, with a smooth transfer except for a
small amplitude, high frequency oscillation on the 14N trajectory.

The correlation orders in Fig. 3 b clearly show a 1-spin order giving way to
the 2-spin order, to mediate the 13C and 14N coupling, then after about half
of the pulse duration, the 2-spin order is transferred back to 1-spin order.
There is an initial ramp in the 2-spin order and then a linear decay to zero
at time $=T$. Again, the correlation orders develop a small amplitude, high
frequency oscillation in the trajectories, this time after a maximum 2-spin
order has been reached.

### § 2 The Fokker-Planck formalism

Kuprov [KUPROV16] set out a review of the Fokker-Planck [FOKKER14, PLANCK17]
formalism used in magnetic resonance simulations, finding the very general
formulation that should be preferred to the Liouville-von Neumann equation.
For this reason, the optimal control methods set out in this thesis should be
constructed to be compatible with the Fokker-Planck equation of motion. As
pointed out in [KUPROV16], there is a close relation between the Fokker-Planck
equation and Floquet theory of the previous section [LEVANTE95] – the latter
is a special case of the former.

#### § 1 Optimal Control in HMQC-type experiment

Figure 4: Trajectory analysis an hmqc-type experiment to indirectly excite the
overtone transition in glycine with hard pulses. Coherence orders are
represented as the effect of a constant amplitude and zero phase. Magnified
section of the plot shows higher coherence orders populated near the end of
the control pulse trajectory (upper)

The testing set out in this section is an hmqc-type experiment (Hetronuclear
Multiple Quantum Coherence) to indirectly excite the overtone transition in
glycine [HAIES15] (results from hard pulses alone shown in Fig. 4). Grape
calculations of optimal control in the Fokker-Planck formalism are used to
find efficient pulses for an excitation of the $+2$ coherence order of 14N in
the hmqc-type experiment. Optimal pulses are designed to be resilient to
radiofrequency pulse miscalibration of $\pm 2.5\%$. The total pulse length of
the optimal set is investigated with respect to numerical convergence, maximum
fidelity after 200 iterations, and coherence order trajectories.

The aim is to find optimal control pulses that spy on 14N through echoes on a
connected 1H by exciting high coherence orders of 14N. The optimal control
pulses form part of a larger pulse sequence, shown in Fig. 5, where the
second section of the optimal control pulses is the time reversed version of
the first section.

Figure 5: Pulse set including optimal control pulse. Pulses vary in phase
only, and the amplitude should be considered constant. The second section of
optimal pulses is the time-reversed version of the first section.

The spinning rate of the system is 35 kHz at the magic angle, in a 14.1 tesla
magnet. Relaxation is included in the form of a diagonal approximation with a
rate of 300 Hz, and the pulse duration is varied with different integer values
of the rotor period. Optimal control of the initial state to a target state is
performed as a sum of states picking out the $+2$ coherence excitation on the
14N. The initial to target states are set out as

$\displaystyle L_{+}^{\text{(H)}}$
$\displaystyle\xrightarrow{\hphantom{aaa}\textsc{ grape }\hphantom{aaa}}$
$\displaystyle
L_{+}^{\text{(H)}}I_{+}^{\text{(N)}}+L_{-}^{\text{(H)}}I_{+}^{\text{(N)}}$

Control operators are defined in the Cartesian representation, allowing pulses
on both 14N and 1H to work in the rotor frame at each of the weighted
crystalline orientations in the rank-5 powder average:

$\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}}}}}_{x}^{\text{(N)}}=\bm{\mathbbm{1}}\otimes\left[\cos(\theta)\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}}}}}_{z}^{\text{(N)}}+\sin(\theta)\left(\frac{\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}}}}}_{+}^{\text{(N)}}+\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}}}}}_{-}^{\text{(N)}}}{2}\right)\right]$
(13a)
$\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}}}}}_{y}^{\text{(N)}}=\bm{\mathbbm{1}}\otimes\left[\cos(\theta)\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}}}}}_{z}^{\text{(N)}}+\sin(\theta)\left(\frac{\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}}}}}_{+}^{\text{(N)}}-\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}}}}}_{-}^{\text{(N)}}}{2i}\right)\right]$
(13b)
$\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}}}}}_{x}^{\text{(H)}}=\bm{\mathbbm{1}}\otimes\left[\cos(\theta)\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}}}}}_{z}^{\text{(H)}}+\sin(\theta)\left(\frac{\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}}}}}_{+}^{\text{(H)}}+\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}}}}}_{-}^{\text{(H)}}}{2}\right)\right]$
(13c)
$\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}}}}}_{y}^{\text{(H)}}=\bm{\mathbbm{1}}\otimes\left[\cos(\theta)\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}}}}}_{z}^{\text{(H)}}+\sin(\theta)\left(\frac{\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}}}}}_{+}^{\text{(H)}}-\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}}}}}_{-}^{\text{(H)}}}{2i}\right)\right]$
(13d)

where $\theta$ is the magic angle of the rotor spinning axis and
$\bm{\mathbbm{1}}$ is the identity operator of size $d$, being the size of the
spatial dimension of the Fokker-Planck formalism.

Amplitude (polar) | total time | pulses |  $\nabla$ calc. | $\%\text{ state transfer }$
---|---|---|---|---
(kHz) |  ($\mu$s) | (#) | (s) |  initial state $=\left\\{L_{+}^{\text{(H)}}\right\\}$
spin rate = 35 kHz |  |  |  targets $=\left\\{L_{+}^{\text{(H)}}L_{+}^{\text{(N)}},L_{-}^{\text{(H)}}L_{+}^{\text{(N)}}\right\\}$
$\oldsqrt[\ ]{2}\times$ 50 | 57.143 | 200 | 47 | 7.89 %
$\oldsqrt[\ ]{2}\times$ 50 | 114.286 | 400 | 92 | 15.86 %
$\oldsqrt[\ ]{2}\times$ 50 | 171.429 | 600 | 140 | 18.67 %
$\oldsqrt[\ ]{2}\times$ 50 | 228.571 | 800 | 189 | 21.57 %
Table 2: Simulation performance of each of four test optimisations: 2, 4, 6,
and 8 rotor periods, each with 100 time slices per rotor period. Note: results
are not converged and a greater fidelity would be achieved if the optimisation
were allowed to run for more than $200$ iterations, particularly for the
longer pulse sets (see Fig. 8)  Figure 6: Phase profile resulting from the
$\ell$-bfgs-grape method with 2, 4, 6, and 8 rotor period for total pulse
duration (left to right). The amplitude profile is constant in all cases,
$A=\oldsqrt[\ ]{2}\cdot 50$ kHz, and the phase pulses on 14N and 1H (upper and
middle) should be tolerant to a $\pm 2.5\%$ miscalibration on the pulse
amplitudes (when translated to a Cartesian representation). The total density
of magnetisation at each spin is shown as an effect of the pulses (lower).
Figure 7: Trajectory analysis of an l-bfgs optimisation with 2, 4, 6, and 8
rotor periods for the total pulse duration (left to right). Coherence orders
are represented as the effect of phase profiles shown in Fig. 6. Magnified
section of the plot shows higher coherence orders populated near the end of
the control pulse trajectory (upper). NOTE: Trajectories are the last
trajectory calculated i.e. that at $+2.5\%$ pulse miscalibration and for the
target state $L_{-}^{\text{(H)}}L_{+}^{\text{(N)}}$.

Normalised initial states and target states are similarly defined via the
Kronecker product:

$\displaystyle\rho(0)=\frac{1}{\oldsqrt[\ ]{d}}\left(\bm{\mathbbm{1}}\otimes
L_{+}^{\text{(H)}}\right)$ (14a) $\displaystyle\rho_{1}(T)=\frac{1}{\oldsqrt[\
]{d}}\left(\bm{\mathbbm{1}}\otimes
L_{+}^{\text{(H)}}L_{+}^{\text{(N)}}\right)$ (14b)
$\displaystyle\rho_{2}(T)=\frac{1}{\oldsqrt[\
]{d}}\left(\bm{\mathbbm{1}}\otimes
L_{-}^{\text{(H)}}L_{+}^{\text{(N)}}\right)$ (14c)

where $\bm{\mathbbm{1}}$ a unity vector of length $d$, $\rho(0)$ is the
initial state at time $t=0$, $\rho_{1}(T)$ and $\rho_{2}(T)$ are the two
target states at the time $t=T$. The $\ell$-bfgs-grape method is used, and
performs two fidelity and gradient calculations from the initial state to the
two final states, then to average the two fidelities and two gradients. The
full experimental state transfer problem can be envisaged as the pulse set
being part of the whole experimental pulse sequence in Fig. 5.

Four optimisation tests were performed with the only difference being the
total pulse time allowed for the control pulses. The time step remained equal
in all cases, having $100$ time steps per rotor period ($35$ kHz spinning
rate). The resulting fidelity from each test is shown in the Tab. 2 (final
column).

The polar amplitude was kept to a constant, $\oldsqrt[\ ]{2}\cdot 50$ kHz
444corresponding to $\pm 50$ kHz maximum/minimum pulse amplitude in the
Cartesian representation, and only the phase profiles of pulses on each spin
were varied. The phase profiles resulting from the optimisation are shown in
Fig. 6. In addition, the coherence pathways for these pulses are shown in
Fig. 7 – this is the main indication of the success of the simulations in
comparison to the hard pulses of Fig. 4.

The four tests were run to $200$ iterations, rather than a converged set of
solutions. Fig. 8 shows how far each of the test optimisations converged. In
addition to showing the final fidelity reached at $200$ iterations, the norm
of the gradient, $\|\nabla\|_{2}$, is shown at each iteration – the norm of
the gradient is a measure of how optimal the solution is 555a maximum or
minimum proper would have $\nabla=0$ and hence $\|\nabla\|_{2}=0$ – the main
point is the smaller this number, the more optimal the solution.

Figure 8: Convergence characteristics. Note that the optimisations have not
yet achieved convergence to a maximum (apart from the optimisation with 2
rotor periods).

## Feedback Control

“Most men will not swim before they are able to. Is that not witty? Naturally,
they won’t swim! They are born for the solid earth, not for the water. And
naturally they won’t think. They are made for life, not for thought. Yes, and
he who thinks, what’s more, he who makes thought his business, he may go far
in it, but he has bartered the solid earth for the water all the same, and one
day he will drown.”

– Hermann Hesse, Steppenwolf

A significant current problem in high-field electron spin resonance (esr)
spectroscopy is the difficulty of achieving uniform and quantitative signal
excitation using microwave pulses [SCHWEIGER01, BRUSTOLON09]. The greatest
instrumentally available microwave power, in terms of the electron spin
nutation frequency, is about 50 MHz; the shortest realistic
${}^{\pi}\\!/_{\\!2}$ pulse is therefore $5$ ns long [CRUICKSHANK09],
corresponding to the excitation bandwidth of $200$ MHz – enough to affect a
significant portion of many solid state esr signals, but insufficient to
excite such signals uniformly and quantitatively. The consequences of partial
excitation of the spectrum include useful orientation selection effects
[SCHIEMANN09, POLYHACH07, BODE08], but also reduced sensitivity and diminished
modulation depth in two-electron dipolar spectroscopy [THURNAUR80, SALIKOV92,
TAN94].

The time resolution of the best available microwave pulse shaping equipment is
of the order of $30$ ps 111Keysight M8196A $92$ GSa/s Arbitrary Waveform
Generator. This work uses Bruker SpinJet awg (Arbitrary Waveform Generator)
with $0.625$ ns time resolution § 1 – in combination with non-resonant
cavities it enables generation of shaped pulses with the bandwidth of about 1
GHz [CRUICKSHANK09] and allows many broadband excitation schemes originally
developed for nuclear magnetic resonance (nmr) spectroscopy [WIMPERIS94] to be
used with only minor modifications [SPINDLER13, DOLL13, DOLL14, SCHOPS15,
JESCHKE15]. Numerically designed “optimal control” microwave pulses
[KHANEJA05, FOUQUIERES11, GOODWIN16] are also possible [SPINDLER12,
KAUFMANN13], but a complication specific to esr is that the waveforms received
by the sample are very different from those sent by the awg – the response
function of the esr instrument cannot be ignored [SPINDLER12].

One way around this is to construct a transfer matrix or a response function
that connects, under the linear response approximation, the ideal pulse
emitted by the computer to the real pulse seen by the sample. The transfer
matrix 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].

The transfer matrix approach is not perfect – different samples alter the
dynamical properties of the resonator in different ways, as does the antenna –
but the linear and quasi-linear models work well in practice. The standard
procedure is to take the desired pulse, reverse-distort it through the
transfer matrices, send the result out of the waveform generator and hope that
a good rendering of the intended pulse shape arrives at the sample point. It
usually does [SPINDLER12, KAUFMANN13], but the logistical overhead of
measuring the transfer matrix is significant. Accurate measurement of the
instrument response function in the ways described above is time-consuming
(hours), and esr resonators, particularly at high frequencies, tend to have
strongly sample-dependent response functions.

In this chapter, a different microwave pulse shape refinement strategy is
explored, which does not require explicit knowledge of the transfer matrix §
1. It relies instead on the possibility of repeating an esr experiment
hundreds of times per second, and recognises the fact, discussed in detail
below, that microwave pulses in esr need very few discretisation points due to
the significant width of the instrument response function. Further details of
this work are in § B, including software code and experimental results.

Numerically optimised microwave pulses are used to increase excitation
efficiency and modulation depth in electron spin resonance experiments
performed on a spectrometer equipped with an arbitrary waveform generator (§
1). The optimisation procedure is sample-specific and reminiscent of the
magnet shimming process used in the early days of nuclear magnetic resonance –
an objective function (for example, echo integral in a spin echo experiment)
is defined and optimised numerically as a function of the pulse waveform
vector using noise-resilient gradient-free methods (§ 2). The resulting
shaped microwave pulses achieve higher excitation bandwidth and better echo
modulation depth than hard pulses. Although the method is theoretically less
sophisticated than quantum optimal control techniques, the rapid electron
relaxation means that the optimisation takes only a few seconds, and the
knowledge of either the instrument response function or spin system ensemble
parameters is not required. This makes the procedure fast, convenient, and
easy to use. Improvements in broadband excitation efficiency, spin echo
intensity [HAHN50] in § 1 and modulation depth are demonstrated with optimal
conditions of oop-eseem [TIMMEL98] in § 2, at the instrument time cost not
exceeding the time it used to take to auto-shim an nmr magnet.

The method is known as “feedback control” [SKOGESTAD07, FRANKLIN94, DOYLE13];
in its electron spin resonance adaptation it is similar in principle to the
well known (in the nmr circles) task of maximising a deuterium lock signal
during the magnet shimming process [DEROME13] – a target variable is chosen
and maximised, using a noise-resilient optimisation algorithm, with respect to
the variables of interest. In the esr case, these variables are amplitudes of
the microwave field at each waveform discretisation point.

### § 1 Feedback Control Optimisation

Figure 1: Block diagram of awg closed-loop feedback control.  a An initial,
random, pulse set sent to the awg.  b Pulse-set is processed and interpolated,
then sent to the sample (resonator).  c The pulse-set arriving at the sample
experiences shape distortions from system noise and hardware limitations.  d
The effect of pulse-set induced excitations on the sample is measured as a
signal.  e The optimisation algorithm calculates a new, modified set of pulses
that are the sent back to the awg.

One simple way to iteratively improve the result of experimental measurements
is to use feedback control [SKOGESTAD07, FRANKLIN94, DOYLE13]. The idea is to
modify inputs to the experiment according to the result of an experimental
measurement: controlling the experiment “on the fly”. The Bruker SpinJet
arbitrary waveform generator (awg) has the ability to shape pulses of the esr
experiment – giving a flexible control method. However, distortions to the
pulse shape are expected in two forms – constant distortions that are the same
for repeated measurements, such as shape pulse interpolation, and those
associated with non-linear effects such as noise from electronic components
and sample specific miscalibration. Fig. 1 sets out a block diagram of a
closed-loop feedback control using an awg.

#### § 1 Hardware and Software

The Bruker SpinJet awg is a two channel arbitrary waveform generator, allowing
the input of in-phase and quadrature components of shaped pulses to a
spectrometer (Fig. 1). The Bruker SpinJet awg used in this work has a time
resolution of $0.625$ ns, $14$-bit amplitude resolution, $1.6$ GS/s sampling
rate, and $\pm 400$ MHz bandwidth around the carrier frequency. The Bruker
EleXSys E580 esr spectrometer has a $2$ ns time base and, in combination with
the awg, resolves the time resolution mismatch by downsampling the pulse
waveform onto a $1$ ns increment time grid.

The software used in this work was written in-house, and has a flow of
communication between Spinach [HOGBEN11] and Xepr Python libraries, shown in
Fig. 2. The master process runs in Matlab and calls Xepr Python functions as
necessary to control the instrument. Experimental data is written by Xepr into
plain text files and subsequently parsed by Matlab. Optimisation restart
capability is implemented using an md5 hash table of the previously submitted
experimental settings and outcomes [GOODWIN16] – an interrupted optimisation
can therefore retrace its steps quickly without re-running previously executed
experiments.

Figure 2: Software flow diagram.  a Set initial pulses. A good guess of an
ideal pulse set based on prior knowledge of the sample.  b 1) Send this pulse-
set to the awg using PulseSPEL: Load, show, and compile shape file. 2) Modify
PulseSPEL variables and compile. 3) Validate and compile PulseSPEL program.  c
Run experiment and wait for experiment completion, then read signal data.  d
Interpret data in optimiser.  e If the optimisation algorithm hasn’t
converged, calculate a new, modified set of pulses that are the sent back to
the awg. Python code is detailed in § 2.

The simplex optimiser allows both Nelder-Mead [NELDER65, LAGARIAS98] and
multidirectional search [TORCZON89, DENNIS91] methods, with the option of
evaluating independent sets of objective variables in parallel [HIGHAM93].

Many shapes can be stored within the awg in the form of a plain-text shape-
file, limited to $262,144$ bytes of physical disk size. A limit of $115$
consecutive shape-file compile events must be enforced to avoid overload of
the awg/console shared memory. Considering this last point, and that there may
be many independent shapes involved in performing an optimisation iteration –
it will become advantageous to perform a 2D experiment with these pre-
calculated, shapes.

Two primary sources of waveform distortions may be identified: the static
ones, introduced by the instrument electronics (pulse shape interpolation,
hardware response function, etc.), and those associated with transient
effects, e.g. electrical noise and sample-specific magnetic susceptibility
effects. The presence and the transient nature of these distortions makes
common esr objective functions (signal integral, modulation depth, etc.)
impossible to differentiate numerically and necessitates the use of gradient-
free noise-resilient optimisation methods.

#### § 2 Optimisation Method

Figure 3: Different simplex steps for the Nelder-Mead and multidirectional
search algorithms, each simplex has $n+1$ vertices where $n$ is the number of
variables. Faces of the initial simplex are scaled to
$\big{|}\bm{c}^{\vphantom{|}}_{0}\big{|}_{\infty}$, and form an equilateral
simplex. Expansion and contraction factors set to $\beta=2$ and $\gamma=0.5$
respectively. The best/worst vertex and the centroid
$\bm{c}^{\vphantom{|}}_{c}=\sum_{0}^{n}\bm{c}^{\vphantom{|}}_{j}$ of the
simplex are used when reflecting, expanding, and contracting.

An optimisation problem usually attempts to find extrema of an objective
function [GILL81, FLETCHER87, NOCEDAL06] – in financial terms these extrema
could be the minimum of a cost function, or a maximum of a profit function.
This function should map a set of variables to a scalar number designed as a
convenient metric characterising the “performance” of a solution formed from
those variables.

The intensity of an echo makes a good echo experiment [SCHWEIGER01,
BRUSTOLON09] and the integral of the in-phase part of the echo signal makes a
convenient metric to maximise. This is defined as the 1-norm of the signal
plot vector [DOYLE13]:

$\|\vec{s}\,\|_{1}:=\int\limits_{0}^{T}\big{|}s(t)\big{|}\mathrm{d}t$ (1)

Any optimisation should be properly scaled [FLETCHER87], and scaling for these
experiments simply involves normalising the objective metric with the signal
produced by a reference experiment; in this case, an experiment with
conventional hard pulses. This gives the metric of the quality of an echo as

$\vphantom{+}{}^{\vphantom{|}}_{\vphantom{|}}J(\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}\bm{c})=\frac{\|\vec{s}\,\|_{1}}{\|\vec{r}\,\|_{1}}$
(2)

where
$\vphantom{+}{}^{\vphantom{|}}_{\vphantom{|}}J(\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}\bm{c})$
is the objective measure to maximise, $\vec{s}$ is the signal measured within
an objective function and $\vec{r}$ is the reference signal measured from an
experiment with hard pulses.
$\vphantom{+}{}^{\vphantom{|}}_{\vphantom{|}}J(\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}\bm{c})>1$
is defined as a good echo, better than that from hard pulses alone, and
$0<\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}J(\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}\bm{c})<1$
is defined as a bad echo, worse than that from hard pulses.

There is no obvious gradient with respect to the control parameters when
considering an experimental measurement. Gradient-free optimisation strategies
[SWANN72, GILL81, FLETCHER87, NOCEDAL06, BRENT13] are the only choice for this
real-time optimisation. A simple Nelder-Mead algorithm [NELDER65, LAGARIAS98]
(a member of a family of direct search methods termed simplex
methods[SWANN72], polytope methods [GILL81] and ad-hoc methods [FLETCHER87])
has the benefit of a relatively small number of experimental evaluations for a
single iteration compared to other popular algorithms e.g. genetic algorithms
[JUDSON92, BRIF10] or simulated annealing [METROPOLIS53]. However, convergence
is not guaranteed and is linear at best [DENNIS87]. A desirable benefit of the
Nelder-Mead algorithm is its tolerance to random noise (on a smooth function)
[GILMORE95, BORTZ98]. In addition; a restart step has been shown to alleviate
stagnation of simplex algorithms [KELLEY99].

The shape file characterising the optimisation variables has the constraint
that it can only contain values in the range $\big{[}-1,+1\big{]}$ and these
constraints can complicate the formulation considered in the optimisation
problem. However, gradients are not considered in for the optimisation problem
stated in this work, and a transform can easily classify the bounded pulse-
shape variables [FLETCHER87]: a repeating, mirrored transform can allow all
real numbers to be evaluated by the objective function in a continuous manner.
This waveform constraint transform can be classified as in Eq. 2:

$\tilde{\bm{c}}=\Bigg{|}\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}\bm{c}-4\bigg{\lfloor}\frac{\vphantom{+}{}^{\vphantom{|}}_{\vphantom{|}}\bm{c}-1}{4}\bigg{\rfloor}-3\Bigg{|}-1$

where $\vphantom{+}{}^{\vphantom{|}}_{\vphantom{|}}\bm{c}$ is the unbounded
waveform and
$\vphantom{+}{}^{\vphantom{|}}_{\vphantom{|}}\bm{c}-n\lfloor\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}\bm{c}/n\rfloor$
is the remainder from dividing
$\vphantom{+}{}^{\vphantom{|}}_{\vphantom{|}}\bm{c}$ by $n$.

A modified Nelder-Mead simplex algorithm, the multidirectional search simplex
method [TORCZON89, DENNIS91, HIGHAM93], has an efficient characteristic in
that it can be used as a partially parallel simplex method, essentially
creating a search grid of simplices. The added computational cost is that many
more function evaluations must be made per simplex, depicted in Fig. 3. This
multidirectional search algorithm has a comparable performance to finite-
difference gradient algorithms and is designed to be tolerant to problems with
many local minima [GILMORE95]. A further useful property of the
multidirectional search algorithm is that is has guaranteed convergence
[TORCZON91], similar to gradient following algorithms. The common algorithmic
steps of the Nelder-Mead and multidirectional search algorithms are shown in
Fig. 3.

To improve stability of the optimisation process, the pulse shape variables
are stored in a database using an md5 hash function during the optimisation
iterations [GOODWIN16] – in the event of the optimisation algorithm halting
mid-algorithm due to unexpected hardware incidents, the optimisation can
effectively be resumed from the same stage after hardware problems have been
resolved.

### § 2 Experimental Application

#### § 1 Feedback Optimised Hahn Echo

Figure 4:  a In-phase,  b quadrature, and  c phase parts of the pulse shape at
each iteration. The spin echo experiment performed on a nitroxide radical and
the optimisation variables are $33$ discrete shape points of a $32$ ns
$\pi$-pulse.  d Real,  e imaginary, and  f magnitude of the measured signal at
each iteration.

A two-pulse Hahn echo [HAHN50, ROWAN65] is a good test for this closed-loop
feedback control, using the 1-norm of the echo signal as a figure of merit, as
in Eq. 1. An optimised $32$ ns $\pi$-pulse is the objective of the test,
using discretised in-phase and quadrature components of a shaped waveform. In
emulating the time resolution of the spectrometer, a shape with $33$ discrete
points should be used as the finest grid; one point at the start of every time
point plus one point at the end of the last time point.

The optimisation starts at a random guess within the allowed bounds, and
proceeds to search for shapes that produce better echoes. This first test of
the feedback control methods uses the Nelder-Mead simplex algorithm, with
termination conditions of a simplex size less than $0.05$, or more than $100$
consecutive contraction steps without finding a better solution (indicative of
stagnation of the algorithm). The echoes produced at each of the optimisation
iterations are shown in Fig. 4.

The results of the feedback control optimisation are not encouraging for these
pulse shapes § 3. The most obvious failure is that the integral of the in-
phase part of the echo does not improve on that from hard pulses. The
optimisation is considered converged to a local maximum, evident in Tabs. 2,
1, the latter part of the optimisation makes many functional evaluations with
very little gain. This is typical of a simplex optimisation which tries many
contraction steps before an improved solution can be found (see Fig. 3).
There are two other characteristics of the optimisation:

  1. 1

The quadrature part of the spin echo indicates an echo without the correct
phase calibration.

  2. 2

The pulse shape does not change by an appreciable amount during the
optimisation.

The first of these characteristics is not so much of a concern: a phase
calibration can be performed after the optimisation and can only give a larger
in-phase part of the echo integral than that shown in Fig. 4 a. The second
characteristic may have two possible reasons: either the variables for the
optimisation are not normalised for an effective optimisation, making too
small steps at each iteration, or there are just too many variables for the
optimiser to effectively handle. The Nelder-Mead algorithm is known to be
effective only for low dimensional problems [LAGARIAS98] – there is no
gradient to direct the algorithm and it can only be a sophisticated ad-hoc
search [FLETCHER87]. Directed by the limits of the Nelder-Mead algorithm,
further tests of lower dimension are required.

Two further tests are performed, one with approximately ${}^{2}\\!/_{\\!3}$ of
the number of time points and one with approximately ${}^{1}\\!/_{\\!3}$ of
the time points. In both cases the total pulse duration remains the same. The
first of these is a $21$ point shape and the results at each iteration of the
Nelder-Mead algorithm are shown in Fig. 5.

Figure 5:  a In-phase,  b quadrature, and  c phase parts of the pulse shape at
each iteration. The spin echo experiment performed on a nitroxide radical and
the optimisation variables are $21$ discrete shape points of a $32$ ns $\pi$
pulse.  d Real,  e imaginary, and  f magnitude of the measured signal at each
iteration.

Results are more encouraging for a pulse shape using $21$ discrete points §
3 – achieving an echo integral 27% better than that of hard pulses (see Tab.
2). However, this optimisation makes many function evaluations giving a total
optimisation time of more than $20$ minutes (much longer than an experienced
spectroscopist would take to find a good echo signal). The majority of the
function evaluations are from inner contraction steps of the simplex algorithm
– making smaller and smaller simplices with little gain. In fact, this
optimisation terminates because the relative simplex size was less than the
defined numerical tolerance. This optimisation can be considered stagnated
within the noise level of the signal.

Figure 6:  a In-phase,  b quadrature, and  c phase parts of the pulse shape at
each iteration. The spin echo experiment performed on a nitroxide radical and
the optimisation variables are $11$ discrete shape points of a $32$ ns $\pi$
pulse.  d Real,  e imaginary, and  f magnitude of the measured signal at each
iteration.

The final feedback controlled spin echo experiment shows a further insight
into a hidden sophistication of the optimisation algorithm, when used within
its limits. Again reducing the amount of dimensions of the optimisation
problem, to approximately $1/3$ of the time points, a shape with $11$ discrete
points is used as the variables to the Nelder-Mead algorithm.

These results, shown in Fig. 6 and Tab. 2, produce a good echo integral in a
reasonable amount of time ($6$-$7$ minutes). They are considered converged for
the same reason as the $22$ point pulse shape – indicating that this number a
variables is low enough for the Nelder-Mead algorithm to function effectively.

A final note, which leads indirectly to the next section of results, is that
the variation of the quadrature part of the echo in Fig. 6 covers all phase
from positive to negative offset, finishing at a zero phase offset.
Investigating further – it is evident from the large variation in the
quadrature channel of the pulse shape that the feedback control optimisation
is trying to compensate for the ineffective initial phase calibration
(calibrated to that for hard pulses without a quadrature element). This is
also evident to a lesser degree from the $22$ point pulse shape in Fig. 5.

In effect, these feedback control optimisations are giving better echo signals
by being able to utilise the quadrature pulse channel at the same time as the
in-phase channel, and the only real optimisation is to find the correct signal
phase offset to give the largest echo integral. An important conclusion is
that the development of the waveform shaping equipment in esr spectroscopy
must at this point focus on improving the instrument response function – there
is little to gain by seeking faster arbitrary waveform generators until the
width of the convolution kernel is reduced. The convolution of the $11$ point
optimum pulse from Fig. 6 with the experimentally measured response function
for our EleXSys II E580 spectrometer is shown in Fig. 7.

Figure 7: In-phase and quadrature parts of pulse-shapes from Fig. 6 a- b,
transformed using a $2$ ns resolution transfer matrix in  c- d, followed by
de-noising with a regularised least squares filter in  e- f.

#### § 2 Feedback Optimised OOP-ESEEM

Figure 8: Out-of-phase ESEEM pulse sequence diagram, showing feedback control
pulse (adapted from [BRUSTOLON09]) Figure 9:  a Amplitude, and  b phase parts
of the wurst shape at each iteration.  c Real, and  d imaginary parts of the
measured oop-eseem signal at each iteration, normalised to the peak intensity
of the real part of the reference signal from hard pulses. Opaque lines show
the raw data and bold lines show data smoothed with a multi-pass Savitzky-
Golay filter [SAVITZY64] ($1^{\text{st}}$ pass: $25$-point cubic;
$2^{\text{nd}}$ pass: $75$-point quintic).

As an application of optimised spin echoes the oop-eseem experiment, shown in
Fig. 8, would need to optimise two echoes, one being a maximum echo at a
delay time in the expected maximum modulation, and the other a minimum at the
first minimum of the modulated signal. The objective metric to optimise would
reflect that of a full oop-eseem experiment, but only at these two delays:

$\vphantom{+}{}^{\vphantom{|}}_{\vphantom{|}}J(\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}\bm{c})=\frac{\|\vec{s}_{d_{1}}\\\
|_{1}^{(g)}}{\|\vec{r}_{d_{1}}\|_{1}^{(g)}}-\frac{\|\vec{s}_{d_{2}}\|_{1}^{(g)}}{\|\vec{r}_{d_{2}}\|_{1}^{(g)}}$
(3)

where $\vec{s}_{d_{1}}$ is the echo with a delay $d_{1}$, corresponding to the
first modulation maximum, $\vec{s}_{d_{2}}$ is the echo with a delay $d_{2}$,
corresponding to the first modulation minimum, and similarly for the hard
pulse reference $\vec{r}$. The first maximum is placed at $16$ ns and the
second at $144$ ns (approximate values are found from an oop-eseem experiment
with hard pulses, shown in Fig. 9). The superscript $(g)$ on the 1-norm
indicates that the integral is only performed over the integrator gate width
(as would be done for a full oop-eseem experiment, set to $32$ ns for these
experiments). The time stepping of the oop-eseem experiments is set to a
coarse $6$ ns, and only $22$ shots-per-point/transient averages are used. In
these feedback controlled optimisations, the initial
${}^{\pi}\\!/_{\\!2}$-pulse is optimised and the $\pi$-pulse is left as a hard
pulse. The reason for investigating the initial pulse character is that it has
previously been investigated [TIMMEL98] and found to be important for oop-
eseem experiments.

Considering the observation that results from the spin echo optimisation are
finding an appropriate signal phase within their quadrature pulse channel, and
that the in-phase pulse channel does not change appreciably from its initial
shape; this should be defined explicitly as part of the optimisation
variables.

Chirped pulses [DOLL13, DOLL14, KHANEJA17] and wurst pulse [SPINDLER12,
SPINDLER13] have an easily controllable phase and will be used as the basis of
pulse shapes in this section of work. A wurst-type chirped pulse can be
characterised by its amplitude and phase as

$\displaystyle\bm{c}^{\vphantom{|}}_{r}=2\pi\oldsqrt[\
]{\frac{\omega_{\text{bw}}}{T}}\left(1-\left|\sin^{a}{\left(\frac{\pi
t}{T}\right)}\right|\right)$
$\displaystyle\bm{c}^{\vphantom{|}}_{\varphi}=\pi\frac{\omega_{\text{bw}}}{T}t^{2}$
(4)

where $T$ is the total pulse duration, $\omega_{\text{bw}}$ is the pulse
excitation bandwidth, $a$ is the sine power used for the amplitude envelope
(set to 25% in these experiments), and $t$ is a time grid from
$-{}^{T}\\!/_{\\!2}$ to $+{}^{T}\\!/_{\\!2}$. The variables of these shapes
used for the objective function are the pulse excitation bandwidth,
$\omega_{\text{bw}}$, the total pulse duration, $T$, the phase offset of the
pulse shape, ${\bm{c}^{\vphantom{|}}_{\varphi}}_{0}$, and the nominal pulse
amplitude of the waveform shape, ${\bm{c}^{\vphantom{|}}_{r}}_{\max}$ (defined
as a percentage of the maximum pulse amplitude in the spectrometer
configuration). These shapes are characterised by few variables, which should
enable simplex methods to perform effectively after appropriate scaling.

A first test of this optimisation using the Nelder-Mead algorithm proved
unstable – the discontinuity of the objective variables forced the optimiser
into local maxima, either contracting onto the pulse shape, the pulse length
or the phase offset. This should be expected with discontinuous variables in
the Nelder-Mead algorithm [FLETCHER87].

Now that the dimensionality of the problem has reduced substantially by using
chirped pulses, the multidirectional search algorithm can now become less
expensive. Furthermore, the pattern-search type nature of the algorithm allows
it to deal effectively with the discontinuous variables of this optimisation.
The results of using this optimisation method are shown in Fig. 9.

The use of a phase offset with chirped-wurst pulses as variables of a
multidirectional search optimisation can be considered quite successful in
increasing the modulation depth of the oop-eseem experiment compared to a
similar experiment with hard pulses alone. The number of iterations is very
low and utilising a partially parallel version of the simplex method reduces
the total time taken to less than $2$ minutes § 3.

The optimisation can be considered converged because the simplex size reduces
to less than the defined tolerance, which is a feature of the multidirectional
search convergence [TORCZON91]. These solutions should not be considered
unique: the optimisation was started from a chirped pulse with a very small
bandwidth, almost a hard pulse. When the optimisation is started from other
random starts, it will find other solutions that give a good modulation depth
– these optimisation methods are not global searches.

It should be noted that the results shown in Fig. 9 are explicit oop-eseem
experiments performed at each improved solution. This is for display in this
thesis rather than a need of the optimisation – this only needs two of the
$256$ echo experiments needed for an oop-eseem experiment. This extra time to
perform the oop-eseem experiment is not included in the timing of Tab. 1.
Furthermore, only $22$ transient averages were used in the echo signal -
because this method takes a very short amount of time, a much larger number of
transient averages could be used, removing much of the noise in the raw data
of Fig. 9 at very little cost compared to increasing the shots-per-point by a
similar amount in the oop-eseem experiment.

Using a simple feedback control with a gradient-free optimisation algorithm of
a shaped pulse can improve the signal intensity of a spin echo esr experiment.
The simplicity, almost naïvety, of this feedback controlled optimisation is
that experiment and sample specific pulse distortions need not be classified,
which could be time-consuming and require rare expertise. Furthermore, this
blind method of finding broadband pulses does not need to be time-consuming
and can become a useful tool, similar to the convenience of an auto-shim.

This work may be considered an extra tool in the esr optimisation toolbox.
Code to produce examples within this manuscript are a combination of Matlab
functions of the Spinach libraries [HOGBEN11] to parse data, optimise an
objective function, and invoke Python scripts to communicate with the awg
using the XeprAPI libraries. The function of Matlab in this context is a
triviality, and code is not shown here. However, using Python scripts to
communicate with the awg in ways set out in this manuscript is not well
documented. These communications scripts are bundled with Spinach. A user may
implement these methods with relative ease and very little specialist
knowledge of the hardware electronics or optimisation methods. Furthermore,
this feedback control can be used in conjunction with transfer function
methods as a fine-tuning of the numerical method [HELLGREN13, EGGER14].

## Summary

Finally, from so little sleeping and so much reading, his brain dried up and
he went completely out of his mind.

– Miguel de Cervantes Saavedra, Don Quixote

A Newton-grape method was developed using efficient calculation of analytical
second directional derivatives. The method was 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 Irreducible SphericalTensor Operators

| $\\!d_{m,2}^{(2)}\\!$ | $d_{m,1}^{(2)}$ | $d_{m,0}^{(2)}$ | $d_{m,-1}^{(2)}$ | $d_{m,-2}^{(2)}$
---|---|---|---|---|---
$\\!d_{2,n}^{(2)}\\!$ | $\frac{(1+\cos{\beta})^{2}}{4}$ | $-\frac{(1+\cos{\beta})\sin{\beta}}{2}$ | $\oldsqrt[\ ]{\frac{3}{8}}\sin^{2}{\beta}$ | $-\frac{(1-\cos{\beta})\sin{\beta}}{2}$ | $\frac{(1-\cos{\beta})^{2}}{4}$
$\\!d_{1,n}^{(2)}\\!$ | $\frac{(1+\cos{\beta})\sin{\beta}}{2}$ | $\frac{\cos{\beta}-1}{2}+\cos^{2}{\beta}$ | $-\oldsqrt[\ ]{\frac{3}{8}}\sin{2\beta}$ | $\frac{\cos{\beta}+1}{2}-\cos^{2}{\beta}$ | $-\frac{(1-\cos{\beta})\sin{\beta}}{2}$
$\\!d_{0,n}^{(2)}\\!$ | $\oldsqrt[\ ]{\frac{3}{8}}\sin^{2}{\beta}$ | $\oldsqrt[\ ]{\frac{3}{8}}\sin{2\beta}$ | $\frac{3\cos^{2}{\beta}-1}{2}$ | $-\oldsqrt[\ ]{\frac{3}{8}}\sin{2\beta}$ | $\oldsqrt[\ ]{\frac{3}{8}}\sin^{2}{\beta}$
$\\!d_{-1,n}^{(2)}\\!$ | $\frac{(1-\cos{\beta})\sin{\beta}}{2}$ | $\frac{\cos{\beta}+1}{2}-\cos^{2}{\beta}$ | $\oldsqrt[\ ]{\frac{3}{8}}\sin{2\beta}$ | $\frac{\cos{\beta}-1}{2}+\cos^{2}{\beta}$ | $-\frac{(1+\cos{\beta})\sin{\beta}}{2}$
$\\!d_{-2,n}^{(2)}\\!$ | $\frac{(1-\cos{\beta})^{2}}{4}$ | $\frac{(1-\cos{\beta})\sin{\beta}}{2}$ | $\oldsqrt[\ ]{\frac{3}{8}}\sin^{2}{\beta}$ | $\frac{(1+\cos{\beta})\sin{\beta}}{2}$ | $\frac{(1+\cos{\beta})^{2}}{4}$
Table 1: Second-rank Wigner functions are defined in terms of reduced functions as $\mathfrak{D}_{m,n}^{(2)}(\alpha,\beta,\gamma)=\operatorname{e}^{-im\alpha}d_{m,n}^{(2)}(\beta)\operatorname{e}^{-in\gamma}$ $\ell$ | $m$ | $Y_{\\!\ell}^{m^{\vphantom{\prime}}}\\!(\theta,\varphi)$ | $Y_{\\!\ell}^{m^{\vphantom{\prime}}}\\!(x,y,z)$ | $\hat{T}_{\ell m}(\hat{L})$
---|---|---|---|---
0 | 0 | $\oldsqrt[\ ]{\frac{1}{4\pi}}$ | $\oldsqrt[\ ]{\frac{1}{4\pi}}$ | $\bm{\mathbbm{1}}$
1 | 1 | $-\oldsqrt[\ ]{\frac{3}{8\pi}}\operatorname{e}^{i\varphi}\sin{\theta}$ | $-\oldsqrt[\ ]{\frac{3}{8\pi}}\frac{x+iy}{r}$ | $-\oldsqrt[\ ]{\frac{1}{2}}\hat{L}_{+}$
1 | 0 | $\oldsqrt[\ ]{\frac{3}{4\pi}}\cos{\theta}$ | $\oldsqrt[\ ]{\frac{3}{4\pi}}\frac{z}{r}$ | $\hat{L}_{z}$
1 | -1 | $\oldsqrt[\ ]{\frac{3}{8\pi}}\operatorname{e}^{-i\varphi}\sin{\theta}$ | $\oldsqrt[\ ]{\frac{3}{8\pi}}\frac{x-iy}{r}$ | $\oldsqrt[\ ]{\frac{1}{2}}\hat{L}_{-}$
2 | 2 | $\oldsqrt[\ ]{\frac{15}{32\pi}}\operatorname{e}^{2i\varphi}\sin^{2}{\theta}$ | $\oldsqrt[\ ]{\frac{15}{32\pi}}\frac{(x+iy)^{2}}{r^{2}}$ | $\frac{1}{2}\hat{L}_{+}^{2}$
2 | 1 | $-\oldsqrt[\ ]{\frac{15}{8\pi}}\operatorname{e}^{i\varphi}\sin{\theta}\cos{\theta}$ | $-\oldsqrt[\ ]{\frac{15}{8\pi}}\frac{(x+iy)z}{r^{2}}$ | $-\frac{1}{2}\big{(}\hat{L}_{z}\hat{L}_{+}+\hat{L}_{+}\hat{L}_{z}\big{)}$
2 | 0 | $\oldsqrt[\ ]{\frac{5}{16\pi}}(3\cos^{2}{\theta}-1)$ | $\oldsqrt[\ ]{\frac{5}{16\pi}}\frac{2z^{2}-x^{2}-y^{2}}{r^{2}}$ | $\oldsqrt[\ ]{\frac{2}{3}}\\!\Big{(}\hat{L}_{z}^{2}-\frac{1}{4}\big{(}\hat{L}_{+}\hat{L}_{-}+\hat{L}_{-}\hat{L}_{+}\big{)}\\!\Big{)}$
2 | -1 | $\oldsqrt[\ ]{\frac{15}{8\pi}}\operatorname{e}^{-i\varphi}\sin{\theta}\cos{\theta}$ | $\oldsqrt[\ ]{\frac{15}{8\pi}}\frac{(x-iy)z}{r^{2}}$ | $\frac{1}{2}\big{(}\hat{L}_{z}\hat{L}_{-}+\hat{L}_{-}\hat{L}_{z}\big{)}$
2 | -2 | $\oldsqrt[\ ]{\frac{15}{32\pi}}\operatorname{e}^{-2i\varphi}\sin^{2}{\theta}$ | $\oldsqrt[\ ]{\frac{15}{32\pi}}\frac{(x-iy)^{2}}{r^{2}}$ | $\frac{1}{2}\hat{L}_{-}^{2}$
Table 2: Spherical harmonics may be written in terms of Cartesian coordinates, then to be replaced by corresponding Cartesian operators, to give rank-1 irreducible spherical tensor operators, for a single spin. $\ell$ | $m$ | $a_{m}^{(\ell)}$ | $\hat{T}_{\ell m}(\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}}}},\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}}}})$
---|---|---|---
0 | 0 | $-\frac{1}{\oldsqrt[\ ]{3}}\big{(}a_{xx}+a_{yy}+a_{zz}\big{)}$ | $-\frac{1}{\oldsqrt[\ ]{3}}\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\\!{\begin{pmatrix}1&0&0\\\ 0&1&0\\\ 0&0&1\end{pmatrix}}\\!\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}}}}$
1 | 1 | $-\frac{1}{2}\big{(}a_{zx}-a_{xz}-i(a_{zy}-a_{yz}\,)\\!\big{)}$ | $-\frac{1}{2}\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\\!{\begin{pmatrix}0&0&-1\\\ 0&0&-i\\\ 1&i&0\end{pmatrix}}\\!\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}}}}$
1 | 0 | $-\frac{i}{\oldsqrt[\ ]{2}}\big{(}a_{xy}-a_{yx}\big{)}$ | $-\frac{1}{\oldsqrt[\ ]{2}}\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\\!{\begin{pmatrix}0&-i&0\\\ i&0&0\\\ 0&0&0\end{pmatrix}}\\!\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}}}}$
1 | -1 | $-\frac{1}{2}\big{(}a_{zx}-a_{xz}+i(a_{zy}-a_{yz}\,)\\!\big{)}$ | $-\frac{1}{2}\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\\!{\begin{pmatrix}0&0&-1\\\ 0&0&i\\\ 1&-i&0\end{pmatrix}}\\!\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 | 2 | $+\frac{1}{2}\big{(}a_{xx}-a_{yy}-i(a_{xy}+a_{yx}\,)\\!\big{)}$ | $+\frac{1}{2}\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\\!{\begin{pmatrix}1&i&0\\\ i&-1&0\\\ 0&0&0\end{pmatrix}}\\!\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 | 1 | $-\frac{1}{2}\big{(}a_{xz}+a_{zx}-i(a_{yz}+a_{zy}\,)\\!\big{)}$ | $-\frac{1}{2}\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\\!{\begin{pmatrix}0&0&1\\\ 0&0&i\\\ 1&i&0\end{pmatrix}}\\!\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 | 0 | $+\frac{i}{\oldsqrt[\ ]{6}}\big{(}2a_{zz}-(a_{xx}+a_{yy}\,)\\!\big{)}$ | $+\frac{1}{\oldsqrt[\ ]{6}}\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\\!{\begin{pmatrix}-1&0&0\\\ 0&-1&0\\\ 0&0&2\end{pmatrix}}\\!\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 | -1 | $+\frac{1}{2}\big{(}a_{xz}+a_{zx}+i(a_{yz}+a_{zy}\,)\\!\big{)}$ | $+\frac{1}{2}\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\\!{\begin{pmatrix}0&0&1\\\ 0&0&-i\\\ 1&-i&0\end{pmatrix}}\\!\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 | -2 | $+\frac{1}{2}\big{(}a_{xx}-a_{yy}+i(a_{xy}+a_{yx}\,)\\!\big{)}$ | $+\frac{1}{2}\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\\!{\begin{pmatrix}1&-i&0\\\ -i&-1&0\\\ 0&0&0\end{pmatrix}}\\!\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}}}}$
Table 3: The coefficients required for the transformation between the standard
$3\times 3$ matrix representation and rank-2 irreducible spherical tensors
operators, for 2 spins.

##  B The Arbitrary WaveformGenerator in Optimisation

### § 1 AWG Schematics

Figure 1: The Bruker SpinJet Schematics.

A schematic block diagram of an awg is shown in Fig. 1 and has the following
function:

  1. 1

A set of pulses, in the form of a specified code, is sent to the awg through
an Ethernet port, processed by on-board computation hardware.

  2. 2

Pulses are formed on an in-phase channel (I channel) and an $90^{\circ}$ out-
of-phase, quadrature channel (Q channel).

$\displaystyle\textsc{awg}_{I}=\cos{\omega_{a}t}$
$\displaystyle\textsc{awg}_{Q}=\sin{\omega_{a}t}$

  3. 3

Two awg channels are connected to the I and Q ports of an IQ-mixer.

  4. 4

The IQ-mixer takes the amplified signal from a microwave source, then
splitting this into an in-phase channel and its quadrature channel

$\displaystyle LO_{I}=\cos{\omega_{o}t}$ $\displaystyle
LO_{Q}=\sin{\omega_{o}t}$

where $\omega_{o}$ is the angular frequency of the carrier, microwave signal
(local oscillator).

  5. 5

Separately for the I and Q channels, the awg channel and corresponding
microwave channel are mixed.

$\displaystyle
I=\cos{(\omega_{o}-\omega_{a})t}+\cos{(\omega_{o}+\omega_{a})t},$
$\displaystyle
Q=\sin{(\omega_{o}-\omega_{a})t}-\sin{(\omega_{o}+\omega_{a})t}$

  6. 6

These two mixed signals are then recombined at the final stage of the IQ-
mixer.

### § 2 Xepr Python code

To programme a feedback loop on to the combination of awg and spectrometer,
Python interpreted XeprAPI scripting can be used for operation of the Xepr
application by Bruker BioSpin. XeprAPI can be easily combined with the
language ProDeL111PROcedure DEscription Language embedded within the Xepr
software used in PulseSPEL files as a back-end to the Xepr console software.

##### Load current experiment

A header script must be invoked at the start of every Xepr Python script, to
import libraries needed and to locate and load the Xepr API module. After the
header, XeprAPI.Xepr() creates the Xepr object and connects to the Xepr
software, the API has to be enabled within the Xepr GUI menu222Processing
$\to$ XeprAPI $\to$ Enable Xepr API. Also common to Xepr scripts is the
command to create/access an experiment object of the Xepr application, invoked
with cur_exp = Xepr.XeprExperiment(), to access the current Xepr experiment or
config = Xepr.XeprExperiment(‘‘AcqHidden’’) to access Xepr’s hidden experiment
object.

The current experiment object must be imported to the current workspace,
assuming there is a current experiment on the spectrometer. The hidden
experiment object will become useful in modifying some parts of the
experiment.

⬇

1# load the Xepr API module into Xepr structure

2Xepr = XeprAPI.Xepr()

3

4# define the current and hidden experiments

5currentExp = Xepr.XeprExperiment()

6hiddenExp = Xepr.XeprExperiment("AcqHidden")

##### Acquire signal data

Measure the 1D or 2D signal as an array of real numbers, corresponding to
relative intensity of the measured signal. The current experiment object
should already be imported to the workspace. The signal data should be stored
in a plain text file to be used later, within Matlab.

Having established a connection to Xepr through the XeprAPI module, and
accessed a current experiment within the Xepr GUI333An initial experiment
should be set up within the Xepr GUI before these scripts are invoked; the
intention is to modify values of variables already set within the Xepr GUI.,
the next task is to modify values of variables that will become the variables
of the objective function.

⬇

1# run compiled experiment, wait until finished

2currentExp.aqExpRunAndWait()

3

4# acquire current signal dataset

5dset = Xepr.XeprDataset()

6

7# separate data into x, real(y), and imag(y)

8ordinate = dset.O

9x_signal = dset.X

10rY_signal = ordinate.real

11iy_signal = ordinate.imag

##### Modify PulseSPEL variables

PulseSPEL definitions from the definition file can be modified by initially
parsing all PulseSPEL parameters contained in the experiment and then matching
them with inputs to the Python script, the matched parameters can then be
modified as instructed. The first input to the Python script is the variable
name and the second is the new value.

⬇

1# get current PulseSPEL variables

2Defs = currentExp.getParam("PlsSPELGlbTxt").value

3Defs = Defs.split("\n")

4

5# replace matched PulseSPEL variables

6for value in Defs:

7 if str(sys.argv[1]) in value:

8 new = str(sys.argv[1])+"␣=␣"+str(sys.argv[2])

9 currentExp["ftEPR.PlsSPELSetVar"].value = new

The script will read comma-separated variable names stored in the first line
of the text file data_varsfile.txt, which can be created from within Matlab,
then match them to similarly ordered comma-separated variable values in the
second line of the text file. Although it is not ideal to pass variables from
one script to another in this way, it is a simple way which works given a
waiting time between commands. The script also contains an error handling for
variables defined in data_varsfile.txt which have not been previously defined
in the Xepr experiment. The command cur_exp[‘‘ftesr.PlsSPELSetVar’’].value =
cmdStr sets the current experiment’s PulseSPEL variables to those scraped from
the externally created data_varsfile.txt file to the string cmdStr.

##### Load shape file

Assuming a plain text shape file has already been created, the code to send
this shape to the experiment involves loading this shape file to the
workspace, then compiling it. The argument of the Python script is the name of
the shape file. When a maximum number of uploads is reached, 115, the current
experiment should be reset (not shown here).

⬇

1# Load shapefile into PulseSPEL

2Xepr.XeprCmds.aqPgShpLoad(str(sys.argv[1]))

3

4# Show shapefile in PulseSPEL window

5Xepr.XeprCmds.aqPgShowShp()

6

7# Compile shapefile into experiment

8Xepr.XeprCmds.aqPgCompile()

The final step is to run the current experiment with the modified variables
with the command cur_exp.aqExpRunAndWait() and wait for the experiment to
finish.

Within an objective function, a Matlab optimisation should modify the values
of the variables in the file data_varsfile.txt at each iteration, then to call
the Python script to send these modified values to the spectrometer. The
result of the modified experiment should pass a metric to be
maximised/minimised back to Matlab’s optimisation; the objective of the
objective function. A simple metric to use would be that which gives the best
echo signal i.e. maximising the echo.

The echo signal data must be read from the Xepr software, again, this is
achieved with the XeprAPI. Script that acquires the Xepr data with the command
dset = Xepr.XeprDataset(), then to sort and store relevant data in three text
files to be read by Matlab’s objective function an interpreted as an echo. A
simple echo signal maximisation would maximise the integral of the echo shape.

The optimal control strategy set out so far is no more than an automatic
calibration, which a moderately experienced experimenter can achieve in very
little time through the Xepr GUI. In fact, with error handling and waiting
times included in the Python scripts outlined so far, each experimental
modification will take a number of seconds: a simplex method would take a long
time to search through this variable space to its optimal experimental
settings. However, this is just a simple example to text the two-way
communication between an optimisation algorithm and its “hardware objective
function”.

### § 3 Feedback Control Results Tables

Experiment: | oop-eseem
---|---
Pulse sequence: | $h\nu\to\frac{\pi}{2}\to\pi\to s(t)$
Optimised Pulse: |  $\frac{\pi}{2}$-pulse
$\frac{\pi}{2}$ pulse length: |  2ns$\leq$variable$\leq$160ns
$\pi$ pulse length: | 32ns
Phase cycling: | 16-step
Transient averages: | 22
Algorithm: | NM | MDS
Pulse shape: | wurst | wurst
Phase offset: | nome |  $-\frac{\pi}{2}<$variable$\leq+\frac{\pi}{2}$
Iteration | $\mathbf{X}$ |  #$f(\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}\bm{c})$ | $\vphantom{+}{}^{\vphantom{|}}_{\vphantom{|}}J(\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}\bm{c})$ | $\mathbf{X}$ |  #$f(\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}\bm{c})$ | $\vphantom{+}{}^{\vphantom{|}}_{\vphantom{|}}J(\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}\bm{c})$
1: | $\mathbf{I}$ | 4 | 0.20 | $\mathbf{I}$ | 5 | 1.25
2: | $\mathbf{E}$ | 2 | 0.39 | $\mathbf{R}$ | 8 | 1.29
3: | $\mathbf{C}$ | 2 | 0.56 | $\mathbf{C}$ | 16 | 1.46
4: | $\mathbf{R}$ | 2 | 0.79 | $\mathbf{C}$ | 8 | 1.63
5: | $\mathbf{C}$ | 2 | 1.03 |  |  |
6: | $\mathbf{S}$ | 7 | 1.20 |  |  |
7: | $\mathbf{S}$ | 7 | 1.20 |  |  |
8: | $\mathbf{S}$ | 7 | 1.20 |  |  |
9: | $\mathbf{C}$ | 2 | 1.20 |  |  |
Total experiments: |  $2\times$42 |  $2\times$37
Total time: | 218s | 99s
Average f(x) time: | 5.2s | 2.7s
Table 1: Results of the feedback control optimisation of oop-eseem experiments
using Nelder-Mead (NM) and multidirectional search (MDS) algorithms.
$\vphantom{+}{}^{\vphantom{|}}_{\vphantom{|}}J(\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}\bm{c})$
is the quality of the echo signal.
#$f(\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}\bm{c})$ is the number of
functional evaluations (2 experiments for oop-eseem). The step is the type of
simplex $\mathbf{X}=$ at that iteration – $\mathbf{I}=$ initial, $\mathbf{R}=$
reflect, $\mathbf{E}=$ expand, $\mathbf{C}=$ contract, $\mathbf{S}=$ shrink.
The termination condition of the algorithms was a simplex size $<0.05$ or no
progress was made after $150$ experiment evaluations were made during the
contraction steps. In each case the objective vale was normalised to the same
measure from the same experiments using conventional hard pulses.

Tabs. 2, 1 outline the results obtained in § 2 for the feedback control of
2-pulse Hahn echo experiments and oop-eseem experiments. The variables of
optimisation are the piecewise triangular shape increments of a $\pi$-pulse in
the echo experiments, and the wurst pulse shape parameters (pulse length,
amplitude, phase offset, and pulse excitation bandwidth from Eq. 4) in the
oop-eseem experiments.

Experiment: | 2-pulse echo
---|---
Pulse sequence: | $\frac{\pi}{2}\to\pi\to s(t)$
Optimised Pulse: |  $\pi$-pulse
$\frac{\pi}{2}$ pulse length: | 16ns
$\pi$ pulse length: | 32ns
Phase cycling: | none
Transient averages: | 200
Algorithm: | NM | NM | NM
Pulse shape: | 33 point | 21 point | 11 point
Phase offset: | none | none | none
Iteration | $\mathbf{X}$ |  #$f(\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}\bm{c})$ | $\vphantom{+}{}^{\vphantom{|}}_{\vphantom{|}}J(\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}\bm{c})$ | $\mathbf{X}$ |  #$f(\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}\bm{c})$ | $\vphantom{+}{}^{\vphantom{|}}_{\vphantom{|}}J(\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}\bm{c})$ | $\mathbf{X}$ |  #$f(\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}\bm{c})$ | $\vphantom{+}{}^{\vphantom{|}}_{\vphantom{|}}J(\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}\bm{c})$
1: | $\mathbf{I}$ | 67 | 0.43 | $\mathbf{I}$ | 43 | 0.75 | $\mathbf{I}$ | 23 | 0.71
2: | $\mathbf{R}$ | 2 | 0.44 | $\mathbf{E}$ | 24 | 0.83 | $\mathbf{E}$ | 4 | 0.71
3: | $\mathbf{E}$ | 20 | 0.47 | $\mathbf{E}$ | 16 | 0.89 | $\mathbf{E}$ | 2 | 0.80
4: | $\mathbf{E}$ | 20 | 0.52 | $\mathbf{E}$ | 18 | 1.09 | $\mathbf{E}$ | 6 | 0.84
5: | $\mathbf{E}$ | 16 | 0.59 | $\mathbf{E}$ | 40 | 1.09 | $\mathbf{E}$ | 2 | 1.02
6: | $\mathbf{E}$ | 26 | 0.72 | $\mathbf{E}$ | 58 | 1.16 | $\mathbf{E}$ | 10 | 1.03
7: | $\mathbf{E}$ | 60 | 0.93 | $\mathbf{E}$ | 58 | 1.21 | $\mathbf{E}$ | 2 | 1.23
8: | $\mathbf{E}$ | 104 | 0.98 | $\mathbf{E}$ | 20 | 1.21 | $\mathbf{E}$ | 28 | 1.27
9: | $\mathbf{R}$ | 24 | 0.98 | $\mathbf{R}$ | 2 | 1.22 | $\mathbf{R}$ | 14 | 1.28
10: | $\mathbf{R}$ | 2 | 0.98 | $\mathbf{E}$ | 36 | 1.24 | $\mathbf{C}$ | 26 | 1.28
11: |  |  |  | $\mathbf{R}$ | 14 | 1.24 |  |  |
12: |  |  |  | $\mathbf{E}$ | 16 | 1.24 |  |  |
13: |  |  |  | $\mathbf{R}$ | 18 | 1.25 |  |  |
14: |  |  |  | $\mathbf{R}$ | 16 | 1.27 |  |  |
Total experiments: | 341 | 389 | 119
Total time: | 1030s | 1273s | 380s
Average f(x) time: | 3.0s | 3.3s | 3.2s
Table 2: Results of the feedback control optimisation of 2-pulse echo
experiments using the Nelder-Mead (NM) algorithm. A number of contraction
steps are omitted from the 2-pulse echo experiments, where no objective
improvement was made.
$\vphantom{+}{}^{\vphantom{|}}_{\vphantom{|}}J(\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}\bm{c})$
is the quality of the echo signal.
#$f(\vphantom{+}^{\vphantom{|}}_{\vphantom{|}}\bm{c})$ is the number of
functional evaluations (1 experiment for the 2-pulse echoes). The step is the
type of simplex $\mathbf{X}=$ at that iteration – $\mathbf{I}=$ initial,
$\mathbf{R}=$ reflect, $\mathbf{E}=$ expand, $\mathbf{C}=$ contract,
$\mathbf{S}=$ shrink. The termination condition of the algorithms was a
simplex size $<0.05$ or no progress was made after $150$ experiment
evaluations were made during the contraction steps. In each case the objective
vale was normalised to the same measure from the same experiments using
conventional hard pulses. Figure 2: Convergence of Nelder-Mead (NM) and
multidirectional search (MDS) algorithms.  a shows the quality of the echo
signal and  c is its change.  b shows the number of functional evaluations and
d is a measure of simplex size.

Echo experiments all start from a random $\pi$-pulse shape and oop-eseem
experiments start with random parameters of the wurst pulse at 32 ns, seeded
from the same random number. All wurst pulse shapes have a amplitude envelope
of 25% rise/fall time. Convergence characteristics of these feedback control
experiments are shown in Fig. 2.