SSupplementary References

# Quantum Control of Radical Pair Dynamics beyond Time-Local Optimization

Farhan T. Chowdhury  0000-0001-8229-2374 Matt C. J. Denton
0009-0009-0246-3387 Daniel C. Bonser  0009-0002-0749-9273 Daniel R. Kattnig
0000-0003-4236-2627 Department of Physics and Astronomy, University of Exeter,
Stocker Road, Exeter EX4 4QL, UK

###### Abstract

We realize arbitrary waveform-based control of spin-selective recombination
reactions of radical pairs in the low magnetic field regime. To this end, we
extend the Gradient Ascent Pulse Engineering (GRAPE) paradigm to allow for
optimizing reaction yields. This overcomes drawbacks of previously suggested
time-local optimization approaches for the reaction control of radical pairs,
which were limited to high biasing fields. We demonstrate how efficient time-
global optimization of the recombination yields can be realized by gradient
based methods augmented by time-blocking, sparse sampling of the yield, and
evaluation of the central single time-step propagators and their Fréchet
derivatives using iterated Trotter-Suzuki splittings. Results are shown for
both a toy model, previously used to demonstrate coherent control of radical
pair reactions in the simpler high-field scenario, and furthermore for a
realistic exciplex-forming donor-acceptor system comprising 16 nuclear spins.
This raises prospects for the spin-control of actual radical pair systems in
ambient magnetic fields, by suppressing or boosting radical reaction yields
using purpose-specific radio-frequency waveforms, paving the way for reaction-
yield-dependent quantum magnetometry and potentially applications of quantum
control to biochemical radical pair reactions. We demonstrate the latter
aspect for two radical pairs implicated in quantum biology.

##  I Introduction

The rapidly evolving field of optimal quantum control (OQC) seeks to
manipulate phenomena at the quantum scale by devising and implementing
perturbations, typically in the form of electromagnetic pulses, to steer a
given quantum system to a desired target. Fuelled by demonstrable successes in
nuclear magnetic resonance (NMR) pulse engineering [1, 2, 3] and the control
of ultra-fast excited-state reactions by laser fields [4], pioneering
developments in “coherent control” were first conceived in chemistry. Today
OQC has developed into a mature discipline that is central to modern quantum
technologies associated with the second quantum revolution [5], with numerous
applications across quantum information processing and quantum metrology [1].
Following on from an earlier proposal for controlling the dynamics of radical
pair reactions through modulating the exchange interaction via an optically
switchable bridging group [6, 7], a recent suggestion extends coherent control
to radical pair reactions mediated through radio-frequency magnetic fields
[8], thereby putting a renewed spotlight on the quantum control of spin-
chemical reactions. In the present work we go on to demonstrate the
implementation of a Gradient Ascent Pulse Engineering (GRAPE)-inspired
approach, made computationally feasible via novel algorithmic modifications,
towards realizing open-loop control of singlet recombination yields for
radical pair dynamics in weak magnetic fields.

Reactions involving the recombination of radical pair intermediates are well-
known to depend on spin degrees of freedom and their intrinsic quantum
dynamics. Specifically, the electronic singlet and triplet states associated
with the unpaired electrons for each radical can undergo coherent inter-
conversion as a consequence of symmetry-breaking interactions—in particular,
the hyperfine interactions with surrounding magnetic nuclei. Quite recently,
quantum beats reflecting the coherent singlet-triplet interconversion in
radical pairs have been directly revealed through pump-push spectroscopy [9].
In the case of singlet and triplet states exhibiting differential chemical
reactivity, which is usually the case since radical pair recombination
preserves spin-multiplicity in the absence of strong spin-orbit coupling, the
spin dynamics is reflected in the reaction yields realized via the singlet and
triplet channels. By coupling to magnetic fields via the Zeeman interaction,
radical pairs furthermore acquire sensitivity to static and oscillatory
applied magnetic fields. The former is central to various magnetic field
effects (MFEs) associated with radical pair reactions, while the latter are
pertinent to the prospect of radical pair reaction control through radio-
frequency magnetic fields. Studies of effects of oscillatory magnetic fields
on radical pair reactions have been realized for chemical model systems and
play a discriminatory role in identifying magnetosensitive radical pair
reactions in spin-biological systems [10, 11].

While a majority of model studies on oscillatory magnetic field effects have
applied monochromatic radio-frequency magnetic fields in the presence of a
strong biasing field [12], often referred to as reaction yield detected
magnetic resonance (RYDMR), some studies have been realized in a weak static
field, in particular for exciplex-forming systems [13, 14, 15, 16, 17]. In
such systems, consisting of an electronically excited donor-acceptor complex
characterized by partial charge transfer and photo-emissivity [18], the
radical pair is in equilibrium with or can populate an exciplex state. In this
way, radical pair dynamics become accessible/measurable via the exciplex
emission, i.e. the recombination yield is proportional to the fluorescence
emission under steady-state conditions [19]. Radical pair mediated mechanisms,
or variations thereof [20], are also hypothesized to underpin a compass sense
in various animals, including migratory songbirds [21]. In this case, the
magnetosensitive radical pair is thought to originate from a photo-induced
electron transfer reaction in the flavo-protein cryptochrome. Even though the
exact nature of the underlying radical pair has remained unclear, the fact
that this compass sense is interruptible by weak monochromatic [22] or broad-
band radio-frequency [23] magnetic fields strongly supports an underlying spin
mediated mechanism. These observations support the prospect of controlling
chemical reactions involving radical pairs, possibly in the near future, in
the spin-biological context, and the challenging low-field regime, provided
that an efficient approach to pulse engineering can be developed, as we set
out to do here.

Recent work [8, 24] suggested the control of radical pair reactions based on
model calculations for a prototypical spin system. They focus on the
traditional RYDMR scenario, for which the applied electromagnetic field is in
resonance with the electronic Zeeman splitting produced by a strong static
field. This perturbs the singlet-triplet interconversion of the radical pair
and thus affects the yield of geminate recombination. Their control approach,
building upon the theoretical approach developed by Sugawara [25], propagates
the spin density operator whereby at every moment the control amplitude is
chosen such that an optimization criterion is bound to not decrease. Since the
algorithm optimizes the controls in a time-local fashion, it is dubbed as
_local optimization_. By construction, the employed algorithm allows
optimizing the trajectory of an observable commuting with the time-independent
dynamics generator, i.e. the drift Hamiltonian, or alternatively, an arbitrary
observable at _one_ chosen moment in time (via back-propagation). While the
combined population of the singlet ($S$) and one of the triplet states
($T_{0}$) is hence controllable over time in the high-field limit, the more
relevant singlet population falls in the second category. Through our present
work, we overcome ostensible drawbacks of the time-local optimization approach
by permitting time-global optimization that simultaneously optimizes the
controls for all times and by focusing on the singlet recombination yield,
i.e. the actual experimental observable (rather than $S$ and
$T_{0}$-populations over time or the singlet probability at a moment of time).
We achieve this through using piece-wise constant control amplitudes and
gradient information, by building upon Gradient Ascent Pulse Engineering
(GRAPE) [26]. To realize quantum control of radical pair dynamics in a weak
static magnetic field, these developments are critical.

Our approach, even in the framework of controls derived for prototypical
models not explicitly considering weakly coupled nuclear spins and spin
relaxation processes, promises robustness in applications to more realistic
system descriptions. This opens avenues for an OQC based model-tomography
approach to distinguish between candidate radical pairs in systems where
magnetosensitivity is phenomenologically observable, but the underpinning
mechanisms remain opaque. Enabling the discriminatory shaping of
magnetosensitive responses could prove to be particularly fruitful for pinning
down key facilitating mechanisms, with applications ranging from designing
organic spintronic materials [27], to enhancing the imaging of magnetic
nanostructures [28]. Further potential avenues for exploration are magnetic
field measurements, e.g. via delayed fluorescence-based organic light emitting
diodes with large magneto-electroluminescence [29], or for medical and
biological applications [30], where selected pathways could potentially be
controlled, possibly to endpoints not reachable via static magnetic field
exposure alone.

We deviate from traditional GRAPE by optimizing reaction yields (rather than a
fidelity measure defined at a given moment in time), accounting for asymmetric
radical pair recombination, which gives rise to non-unitary propagators, using
exact gradients rather than approximate ones, and, following de Fouquieres et
al. [31], by using curvature information of the loss function (rather than
steepest decent/ascent). To overcome a major objection to GRAPE-based reaction
control of radical pairs, namely high computational demands, as also voiced by
Masuzawa et al. [8], we suggest a number of practical optimizations such as a
block-optimization scheme and sparse sampling of the reaction yield to improve
the efficiency of the approach while retaining adequate control fidelity. With
these improvements, realistically large radical pair systems can be yield-
controlled in weak magnetic fields. We demonstrate the success of the approach
for a large spin system in an exciplex-forming complex and further provide
examples motivated by quantum biology.

##  II Theory

###  II.1 Radical Pair Spin Dynamics

We consider a system comprising radicals A∙- and B∙+ subject to singlet-
triplet interconversion, undergoing spin-selective recombination reactions as
per Fig. 1.

Figure 1: Reaction Scheme.

The corresponding spin dynamics is described in terms of the time-dependent
spin density matrix $\hat{\rho}(t)$ by

$\displaystyle\frac{d\hat{\rho}(t)}{dt}=-i[\hat{H}(t),{\hat{\rho}(t)}]-\\{\hat{K},\hat{\rho}(t)\\},$
(1)

with $\hat{H}(t)$ denoting the spin Hamiltonian and

$\displaystyle\hat{K}=\frac{k_{\rm S}}{2}\hat{P}_{\rm S}+\frac{k_{\rm
T}}{2}\hat{P}_{\rm T}$ (2)

accounting for radical pair recombination. Here, $k_{\rm S}$ and $k_{\rm T}$
denote the reaction rate constants in the singlet and triplet state,
respectively, and $\hat{P}_{\rm S},\hat{P}_{\rm T}$ are the projection
operators onto the singlet and triplet states. We introduce an effective
Hamiltonian given by

$\displaystyle\hat{A}(t)=\hat{H}(t)-i\hat{K},$ (3)

allowing us to reformulate eq. (1) as

$\displaystyle\frac{\partial\hat{\rho}(t)}{\partial
t}=-i\llbracket\hat{A}(t),\hat{\rho}(t)\rrbracket$
$\displaystyle=-i(\hat{A}(t)\hat{\rho}(t)-\hat{\rho}(t)\hat{A}(T)^{\dagger}),$
(4)

where $\llbracket A,B\rrbracket=AB-(AB)^{\dagger}=AB-BA^{\dagger}$ $\forall$
$B=B^{\dagger}$. The formal solution to eq. (4) is given by

$\displaystyle\hat{\rho}(t_{n})=\hat{U}(t_{n},0)\hat{\rho}(0)\hat{U}^{\dagger}(t_{n},0)$
(5)

where the time evolution operator in terms of $\hat{A}(t)$ is

$\displaystyle\hat{U}(t,0)=\mathcal{T}\mathrm{exp}\left[-i\int_{0}^{t}\hat{A}(\tau)\mathrm{d}\tau\right].$
(6)

Here, $\mathcal{T}$ denotes the time ordering operator, that reorders products
of time-dependent operators such that their time arguments decrease going from
left to right.

The singlet channel recombination yield is given by

$\displaystyle Y_{\rm S}=k_{\rm S}\int_{0}^{\infty}p_{\rm S}(t)dt$ (7)

where $p_{\rm S}(t)=Tr[\hat{P}_{\rm S}\hat{\rho}(t)]$ is the survived singlet
probability of the radical pair. Note that eq. (4) is of truncated Lindblad
form (by introducing shelving states, an equivalent formulation in traditional
Lindblad form is also possible [32]). Thus, the essentially open quantum
system description of radical pair spin dynamics in the presence of spin-
selective recombination eq. (1) can be treated in terms of the dynamics of a
closed system, albeit with a non-Hermitian dynamics generator, non-unitary
evolution operators and non-conserved trace.

Here, we consider radical pairs for which the static part of the Hamiltonian,
$H_{0}$, is of the form

$\displaystyle\hat{H}_{0}=\hat{H}_{A}+\hat{H}_{B}+\hat{H}_{AB}$ (8)

where $\hat{H}_{i},i\in\\{A,B\\}$ is local to radical $i$ and comprises the
Zeeman interaction with the static magnetic field and isotropic hyperfine
couplings with surrounding nuclear spins. $\hat{H}_{AB}$ accounts for inter-
radical couplings in the form of the exchange interaction. Specifically, each
radical is described in terms of corresponding nuclear and electron spin
angular momentum operators $\hat{\bf I}_{ik}$ and $\hat{\bf S}_{i}$ by

$\displaystyle\hat{H}_{i}={\omega}_{i}\hat{S}_{iz}+\sum_{k=1}^{n_{i}}a_{ik}\hat{\bf
I}_{ik}\cdot\hat{\bf S}_{i}.$ (9)

and

$\displaystyle\hat{H}_{AB}=-j_{ex}\left(2\hat{\bf S}_{A}\cdot\hat{\bf
S}_{B}+\frac{1}{2}\right)$ (10)

where ${\omega}_{i}=-\gamma_{i}B$, with $\gamma_{i}$ denoting the gyromagnetic
ratio of the electron in radical $i$ and $B$ the applied magnetic field,
$a_{ik}$ is the isotropic hyperfine coupling constant between electron spin
$i$ and the $k$-th nuclear spin (out of total $n_{i}$), and the exchange
coupling constant is denoted by $j_{ex}$. The initial density operator is
assumed as

$\displaystyle\hat{\rho}(0)=\frac{1}{Z_{1}Z_{2}}\hat{P}_{\rm S},$ (11)

where $Z_{i}=\prod_{k=1}^{n_{i}}(2I_{ik}+1)$ is the total number of nuclear
spin states in radical $i$ and

$\displaystyle\hat{P}_{\rm S}=\frac{1}{4}-\hat{\bf S}_{1}\cdot\hat{\bf
S}_{2},$ (12)

is the singlet projection operator. In the context of spin-biological systems,
we also consider directional magnetic field effects of systems with
anisotropic hyperfine interactions. In this case, $\hat{H}_{i}$ is of the form

$\displaystyle\hat{H}_{i}={\omega}_{i}{\bf n}(\vartheta,\varphi)\cdot\hat{\bf
S}_{i}+\sum_{k=1}^{n_{i}}\hat{\bf I}_{ik}\cdot{\bf A}_{ik}\cdot\hat{\bf
S}_{i}.$ (13)

where ${\bf n}(\vartheta,\varphi)$ is an unit vector with polar angle
$\vartheta$ and azimuth $\varphi$ defining the magnetic field direction and
the ${\bf A}_{ik}$s are the hyperfine tensors.

###  II.2 Gradient-based coherent control

We here outline pertinent features and the required extensions for our GRAPE-
inspired approach. Borrowing its basic set up from GRAPE, our approach entails
discretizing controls into time slices that are improved in a concurrent-
update scheme, i.e. involving controls for all time points (thus “time-
global”), with each update determined only by information from the previous
iteration, so they can be evaluated at each point independently. Specifically
in GRAPE, the gradient of the fidelity (usually defined at a predefined final
time) with respect to all control variables is calculated (or rather
approximated) in each iteration and a step is made in the variable space in
the direction of steepest ascent (or descent). For quantum systems
controllable via convex optimization, first-order gradient information can
sufficiently guarantee that the global maximum at the top of the control
landscape can be located from any given direction, with the extremes of a
given control landscape corresponding to either perfect control or no-control
[33].

Without loss of generality, we can describe a given quantum system undergoing
spin dynamics under coherent control with

$\displaystyle\hat{H}(t)=\hat{H}_{0}+\sum_{l=1}^{M}{u}_{l}(t)\hat{V}_{l},$
(14)

where $\hat{H}_{0}$ is the drift Hamiltonian describing the time-independent
part of the system, i.e. the intrinsic radical pair dynamics in the static
bias field, as introduced above, ${u}_{l}(t)$ are the control amplitudes and
$\hat{V}_{l}$ is the set of control Hamiltonians that couple the controls to
the system, for example via the Zeeman interaction. We control parameters of
only the coherent evolution of the open system dynamics of the radical pair
undergoing spin selective recombination reactions. In principle, the
recombination could also be directly controlled. However, direct control of
recombination is less amenable to experimental manipulation and therefore not
pursued here.

The goal of OQC is to optimize some control objective functional
$G[\\{{u}_{l}(t)\\}]$, with a control target subject to physical constraints
and other criteria. For our purposes, we focus on minimizing or maximizing the
singlet recombination yield of the radical pair, i.e. $G=Y_{\rm S}$. However,
the approach is easily generalized to focus on other reaction outcomes, e.g.
accumulated nuclear polarization. Significant prior work has been concentrated
on optimizing $G$ subject to the assumption that the values taken by the
controls over time may be parameterized by piece-wise constant control
amplitudes in the time domain [34]. This paradigm, which is central to GRAPE
as the first widely applied quantum control algorithm [26], is well aligned
with the practical aspect of controlling radical pairs, as standard arbitrary
waveform generator (AWG) based control schemes represent pulse inputs as
piece-wise constant functions. We thus discretize the time axis,
$t_{0}=0<t_{1}<t_{2}<\cdots<t_{N}=t_{max}$ and we set

$\displaystyle\hat{H}(t)\approx\hat{H}_{0}+\sum_{l=1}^{M}\sum_{n=1}^{N}u_{l,n}\chi_{n}(t)\hat{V_{l}}$
(15)

with constants ${u}_{l,n}$ and the indicator function $\chi_{n}(t)$ defined
such that $\chi_{n}(t)=\begin{cases}1&\text{for $t_{n-1}\leq t<t_{n}$}\\\
0&\text{otherwise}.\end{cases}.$

Guess controls ${u}_{l}(t)$ Time evolve with given controls eq. (17) Evaluate
cost eq. (27) Precondition Gradients via L-BFGS-B Efficient propagator
evaluation via Trotter-Suzuki Splitting eq. (A4) Compute Gradients eq. (23)
Output convergence

Figure 2: A simplified flowchart for optimizing radical pair reaction yields.
Initially, the control variables are guessed, and the system is time evolved.
The fidelity is then computed in terms of an objective function, and we get an
output where convergence is achieved. Otherwise, the gradient with respect to
control variables is computed and utilized, together with iteratively
accumulated curvature information, to update the controls, using which the
system is again time evolved and the process repeats.

The corresponding propagators $\hat{U}(t_{n},0)$ are discretized on the
temporal grid such that

$\displaystyle\hat{U}(t_{n},0)=\mathcal{T}\prod_{i=1}^{n}\hat{U}_{i}=\hat{U}_{n}\cdots\hat{U}_{1}$
(16)

where

$\displaystyle\hat{U}_{i}=e^{-i\hat{A}_{i}\Delta t_{i}}$ (17)

with $\Delta t_{i}=t_{i}-t_{i-1}$ and $\hat{A}_{i}=\hat{A}(t_{i-1})$. The
gradients of $G$ depend on the propagator derivatives, evaluated as

$\displaystyle\frac{\partial\hat{U}_{i}}{\partial{u}_{l,j}}=\delta_{i,j}L_{i,l}=\delta_{i,j}L(i\hat{A}_{i}\Delta
t_{i},i\hat{V_{l}}\Delta t_{i})$ (18)

where $L(X,E)$ is the Fréchet derivative of the matrix exponent defined, for
every $E$, by

$\displaystyle e^{(X+E)}-e^{X}=L(X,E)+o(\lVert E\rVert).$ (19)

More explicitly, we have [35]

$\displaystyle L(X,E)=\int_{0}^{1}{e^{X(1-s)}}Ee^{Xs}ds$
$\displaystyle=\sum_{n=1}^{\infty}\frac{1}{n!}\sum_{j=1}^{n}X^{j-1}EX^{n-j},$
(20)

which also allows us to demonstrate that

$\displaystyle\frac{\partial\hat{U}_{i}^{\dagger}}{\partial{u}_{l,j}}=\delta_{i,j}L(+i\hat{A}_{i}^{\dagger}\Delta
t_{i};+i\hat{V_{l}}\Delta t_{i})$
$\displaystyle=\delta_{i,j}L^{\dagger}(-i\hat{A}_{i}\Delta
t_{i};-i\hat{V_{l}}\Delta t_{i})$
$\displaystyle=\delta_{i,j}L_{i,l}^{\dagger}$ (21)

Combining the above equations, the gradient of the survived singlet
probability,

$\displaystyle p_{\rm S}(t_{n})=Tr[\hat{P}_{\rm
S}\hat{U}(t_{n},0)\hat{\rho}(0)\hat{U}^{\dagger}(t_{n},0)],$ (22)

for $k\leq n$, is thus

$\displaystyle\frac{\partial p_{\rm
S}}{\partial{u}_{l,k}}=-i\textrm{Tr}(\hat{P}_{\rm
S}\hat{U}_{+}\llbracket\hat{W}_{k,l},\hat{U}_{-}\hat{\rho}_{0}\hat{U}_{-}^{\dagger}\rrbracket\hat{U}_{+}^{\dagger})$
$\displaystyle=-i\textrm{Tr}(\hat{U}_{+}^{\dagger}\hat{P}_{\rm
S}\hat{U}_{+}\llbracket\hat{W}_{k,l},\hat{U}_{-}\hat{\rho}_{0}\hat{U}^{\dagger}_{-}\rrbracket)$
$\displaystyle=2\Im[\textrm{Tr}(\hat{U}_{+}^{\dagger}\hat{P}_{\rm
S}\hat{U}_{+}\hat{W}_{k,l}\hat{U}_{-}\hat{\rho}_{0}\hat{U}^{\dagger}_{-})]$
$\displaystyle=2\Im[\textrm{Tr}(\hat{P}_{\rm
S}\hat{U}_{+}\hat{W}_{k,l}\hat{U}_{-}\hat{\rho}_{0}\hat{U}^{\dagger})]$ (23)

where $\hat{W}_{i,l}=i\hat{U}_{i}^{-1}L_{i,l}$ and

$\displaystyle\hat{U}_{+}=\mathcal{T}\prod_{m=k}^{n}\hat{U}_{m}$ (24)
$\displaystyle\hat{U}_{-}=\mathcal{T}\prod_{m=1}^{k-1}\hat{U}_{m}.$ (25)

For $k>n$, the gradient obviously vanishes. In the original formulation of
GRAPE, a first-order approximation of $L_{k,l}$ was used [26], which is
obtained with $\hat{W}_{i,l}=\hat{V_{l}}$. For our approach, we use the full
Fréchet derivative. Finally, this allows us to approximate the singlet yield
by

$\displaystyle Y_{\rm S}=k_{\rm S}\int_{0}^{\infty}p_{\rm S}(t)dt$
$\displaystyle\approx k_{\rm S}\sum_{n=0}^{N}{w}_{n}p_{\rm S}(t_{n})$ (26)

with ${w}_{n}$ denoting suitably chose quadrature weights. The gradient of the
singlet yield is then obtained as

$\displaystyle\frac{\partial Y_{\rm S}}{\partial{u}_{l,k}}\approx k_{\rm
S}\sum_{n=k}^{N}{w}_{n}\frac{\partial p_{\rm S}}{\partial{u}_{l,k}}.$ (27)

##  III Results

Figure 3: Schematic illustration of the time-blocking and sparse-sampling
approach used to reduce the computational expense of gradient-based
optimization of the radical pair recombination yield. The control amplitudes
${u}_{\text{l,n}}$ are simultaneously updated, whereby $B$ disjoint blocks of
$\frac{N}{B}$ time-steps each are considered in succession, with the initial
state taken to be that of the system subject to the preceding optimized
controls. A sparsely sampled approximation of the reaction yields based on
$\frac{N}{BS}$ samples per block suffices to effectively optimize the reaction
yield, in particular for larger spin systems for which the singlet probability
is barely oscillatory. The bottom panel illustrates the unperturbed (orange
solid line) and minimized (blue) singlet probability for Py/DCB, an 18 spin
exciplex system in a static biasing field of $1\mskip 3.0mu$mT. The RF control
field of $0.1\mskip 3.0mu$mT amplitude comprised $1024$ control pulses of
$\mskip 3.0mu1$ns width, the amplitude of which was optimized such that the
singlet yield approximately calculated from $256$ samples of the singlet
probability was minimized assuming a recombination rate of $k=1\mskip
3.0mu\mu$s${}^{-}1$. The Py/DCB system is discussed in more detail in the
Results section.

The control of radical pair dynamics using time-global optimization methods is
expected to be computationally costly, which is why the possibility of
adapting GRAPE was dismissed in favour of time-local optimization in earlier
attempts [8, 24] at realizing radical pair-reaction control. The computational
demand is particularly large if, as we do here (and unlike [8, 24]), the goal
is to optimize yields, which are time-integral quantities (cf. eq. (7)),
unlike the typical fidelity measures in quantum control applications
pertaining to the implementation of a state transfer or a unitary operation
for a fixed time $t_{\rm max}$. Nevertheless, by leveraging various
algorithmic modifications, we succeed in demonstrating the feasibility of a
time-global approach for controlling the dynamics of radical pair spin systems
of realistic complexity.

###  III.1 Algorithmic efficiency improvements

We require computing $N$ single-time-step propagators alongside their Fréchet
derivatives with respect to $L$ control Hamiltonians, and scores of repeated
matrix multiplications to evaluate the time-dependent singlet-probabilities
and associated gradients, which sum to the recombination yield and its
associated gradients. In the Appendix we evaluate the number of required
matrix multiplications (scaling with $N^{2}$), leading to the key insight that
for the considered scenario and typical $N$s and problem sizes the cost of
gradient/expectation value evaluation will exceed that of evaluating the
elementary propagators and their derivatives (scaling with $N$). To allow for
efficiently optimizing the reaction yield of radical pairs, we thus suggest
the adaptation of the gradient-based optimization procedure as illustrated in
Fig. 3.

First, we argue that the optimization can be realized in practice in terms of
$B$ disjoint blocks of $N/B$ time-steps each (time blocking). Second, we
suggest that a sparsely sampled, i.e. crude, approximation of the reaction
yield based on $N/(BS)$ samples per block, i.e. sampling every $S$-th time-
step, is sufficient to adequately optimize the reaction yield. Here, $B$ and
$S$ are integers that divide $N$ and $N/B$, respectively. While this approach
formally does not yield the optimum of the yield as defined in eq. (7),
henceforth referred to as the complete optimum, we argue that, for suitably
chosen $B$ and $S$, yields optimized through sparsely sampled time blocking,
realized at a fraction of computational effort, are sufficiently close to the
true optimum. Third, to efficiently calculate the exponential propagators and
the associated Fréchet derivatives, we use a scaling and squaring approach in
combination with a Trotter-Suzuki splitting of the elementary propagator, as
detailed in the Appendix. Lastly, to accelerate the optimization we utilise
L-BFGS-B, a limited memory variant of the Broyden–Fletcher–Goldfarb–Shanno
algorithm, which preconditions gradients with curvature information in lieu of
steepest descent [31].

###  III.2 Protocol applied to a prototypical spin system

Figure 4: Controlled recombination dynamics of a 7-spin radical pair system
under $0.5\mskip 3.0mu$mT and $10\mskip 3.0mu$mT static biasing fields in (a)
and (b), respectively. The plots depict the survived singlet probability (cf.
eq. 23) with (blue) and without (orange) coherent control, scaled to remove
the spin-independent decay, i.e. $p_{\rm S}(t)\exp(k_{f}t)$. The recombination
yield was minimized over $N=5000$ time-steps split into $B=2$ blocks of $2500$
time-steps each, sampled sparsely at $S=5$ and $S=10$ for (a) and (b),
respectively. Accompanying spectrograms display the control field magnitudes,
computed via discrete Fourier transforms of 64 points, overlapped by 60
points, and windowed using a Tukey function with a shape parameter of 0.25.
Using a control amplitude of $0.25\mskip 3.0mu$mT, the singlet recombination
yield could be reduced to $0.156$ (from $0.269$ without control) and $0.097$
(from $0.346$) for the two biasing fields.

We consider the minimization of singlet recombination yield for a
generalization of the prototypical radical pair studied earlier in the work by
Masuzawa et al. [8], where coherent control of a fidelity measure, rather than
the reaction yield, was achieved in the simpler high-field scenario using
time-local optimization. Specifically, we assume a radical pair comprising 7
spin-$\frac{1}{2}$ particles, coupled through isotropic hyperfine ($0.2$,
$0.5$, $1\mskip 3.0mu$mT and $0.2$ & $0.3\mskip 3.0mu$mT for the two radicals,
respectively) and exchange interactions ($j_{ex}/(2\pi)=1\mskip 3.0mu$MHz),
undergoing spin-selective recombination in the singlet state with rate
constant $k_{b}=1\mskip 3.0mu\mu$s-1 or spin independent escape with rate
constant $k_{f}=1\mskip 3.0mu\mu$s-1, such that $k_{\rm S}=k_{b}+k_{f}$ and
$k_{\rm T}=k_{f}$. Given the high computational demands of completely
optimizing the 7-spin system, we have additionally studied the 5-spin system
resulting from the 7-spin system by retaining the three largest hyperfine
coupling constants. Unlike Masuzawa et al. [8], we seek to directly minimise
the singlet recombination yield

$\displaystyle Y_{b}=\frac{k_{b}}{k_{b}+k_{f}}Y_{\rm S}$ (28)

through a piece-wise constant control field applied perpendicular to a static
biasing field, i.e. $L=1$ and $V\equiv
V_{1}=\omega_{1}(\hat{S}_{1,x}+\hat{S}_{2,X})$, where $\omega_{1}$ is the
maximal Rabi frequency of the control field; control amplitudes are subject to
$|{u}_{1,n}|\equiv|{u}_{n}|\leq 1$. We chose a discretization time step of
$1\mskip 3.0mu$ns and controlled the first $5\mu\mskip 3.0mu$s ($N=5000$)
after initialization of the spin system in the singlet electronic state (thus
accounting for more than 99.3% of radical decay; final yields still evaluated
in the $t\to\infty\sim t_{max}$ limit with $t_{max}=14\mu\mskip 3.0mu$s). Fig.
4 illustrates exemplary time dependent singlet probabilities, controlled and
unperturbed, associated with this scenario for the biasing field intensities
$B_{0}=0.5\mskip 3.0mu$mT and $10\mskip 3.0mu$mT.

Figure 5: For the 7 and 5-spin model systems, a maximal control field
amplitude of $0.25\mskip 3.0mu$mT, with biasing fields of $0.5$ and $10\mskip
3.0mu$mT for (a) and (b), respectively, we show singlet recombination yield as
a function of $S$ for yield minimizations using $B$ blocks of $N/B$ steps
each, the latter value reported in the captions ($N=5000$). The plots present
the best minimization result from $6$ independent replications with randomly
chosen initial control amplitudes (chosen from a normal random distribution
truncated at $\pm 1$); with dots indicate the 75-percentile range. The best
optimization results allowed a reduction of the singlet recombination yield
from $0.269$ to $0.156$ and from $0.348$ to $0.071$ for (a) and (b),
respectively.

Efficient evaluation of single time-step propagators, and their Fréchet
derivatives, was made possible via iterated Trotter-Suzuki splittings, and
time blocking and sparse sampling of reaction yields, as introduced above. In
view of Fig. 5 and comparable datasets for other biasing fields provided in
the Supplementary Material (cf. Fig. S4), it is apparent that excellent
control results can be achieved for all $N/B$ block sizes exceeding $500$.
Neither $B=1$ and $N=B$ are strict necessities for practical control of the
reaction outcome, substantially smaller blocks yield comparable optimization
results at a fractional cost. On the other hand, block sizes of $50$ and $100$
are unsatisfactory, indicating that a time-local optimization is insufficient
to yield optimal control results in the low to intermediate field region. For
these inadequate block sizes, the results of the utilised L-BFGS-B
optimization is furthermore strongly dependent on the initial condition. Note
that intermediate block sizes $N/B$ are occasionally preferred over larger
optimization blocks, as the latter are less prone to get stuck in local
minima, such that the overall optimization outcome can be superior. This is,
for example, seen in Fig. S4 for a biasing field of $1\mskip 3.0mu$mT, for
which the best result of 6 replicated optimizations with $B=2$ and $B=5$
happened to surpass those with $B=1$. Note furthermore that the optimization
results for smaller biasing fields, e.g. $0.5\mskip 3.0mu$mT, is more
sensitive to the sampling interval, $S$, than for $10\mskip 3.0mu$mT. This is
not surprising in view of the more varied singlet probability as a function of
time observed for the lower field. This is encouraging insofar as more complex
spin systems with diverse hyperfine interactions are expected to be less
oscillatory. In any case, it appears that sampling intervals $S$ of $25$ and
less are adequate here.

For the 7-spin system we found it impractical to realize a complete set of
optimizations for $S=1$ and $B=1$, i.e. complete, unblocked optimization, as
optimization runs initiated from random configurations required more than a
week to finish on the hardware utilized. However, a complete optimization
starting from the best blocked optimization finished swiftly for the used
tolerances without yielding a significant additional improvement. Furthermore,
a single complete optimization yielded results comparable to those of the good
blocked and sparsely sampled optimizations. For the smaller 5-spin system
system, it is apparent that the sparsely sampled, blocked optimizations
smoothly converge to the completely sampled optimum and that, again, $B>1$ and
$S>1$ allow excellent practical optimization outcomes relative to this
(expensive) limit. See Fig. 5b and the Supplemental Material for these and
further optimization results achieved for the 5 and 7-spin systems in weak
magnetic fields, including the geomagnetic field (cf. Fig. S1 and Fig. S2).

Figure 6: Dependence on the optimization outcome on the field strength of the
biasing field and driving field, and the control speed for the 7-spin model
system. In (a), we consider a $0.25\mskip 3.0mu$mT control field, and show
yield minimization against increasing $B_{0}$ (blue) in comparison to the un-
controlled reference (orange). In (b) we show yield minimization against
increasing control field amplitude, $B_{1}$, for a biasing field of
$B_{0}=10\mskip 3.0mu$mT and $50\mskip 3.0mu\mu$T. In (c) we show the
dependence of the singlet recombination yield at $B_{0}=1\mskip 3.0mu$mT on
the speed limit $\dot{u}_{max}^{-1}$. For all optimizations, $N=5000$, $B=5$
and $S=25$.

With the approach established, we set out to explore the prospects offered by
controlling radical pair reaction yields in the weak to intermediate field
case. Fig. 6a illustrates the results of minimizing the singlet recombination
yield as a function of the applied static magnetic field relative to the
unperturbed scenario. These results were realized using $5$ blocks of $1000$
time steps of $1\mskip 3.0mu$ns and sampling the recombination yield from
every $25$th time point. The maximal control amplitude was $0.25\mskip
3.0mu$mT. The plotted whiskers indicate 75-percentiles from 7 optimizations
starting from random initial controls. Evidently, the recombination yield can
be markedly reduced for all applied fields, the lowest yield realizable being
$0.0966$, realized for a biasing field of $9\mskip 3.0mu$mT. This amounts to a
reduction by more than a factor of $3.5$ compared to the uncontrolled yield at
the same field ($0.346$). Generally, it is remarkable that the yields
controlled for minimal recombination are all markedly lower than the minimal
recombination yield of $0.252$ realizable by applying a static magnetic field
of $0.2\mskip 3.0mu$mT, i.e. the field matching the dip in the singlet yield
attributed to the low-field effect.

We have further analysed the dependence of the optimization outcome on the
strength of the control field. Fig. 6b shows $Y_{b}$ as a function of the
maximal control amplitude for the case of a static biasing field of $10\mskip
3.0mu$mT and $50\mskip 3.0mu\mu$T, again realized by optimizing $5$ blocks of
$1000$ time steps, sampling every $25^{\mathrm{th}}$. Appreciable minimization
for $B_{0}=10\mskip 3.0mu$mT requires control amplitudes of $0.1\mskip
3.0mu$mT; controls with field strength below $\approx 50\mskip 3.0mu\mu$T do
not elicit appreciable effects. This finding is in line with the well-studied
effects of radio-frequency magnetic fields on radical pairs. In this scenario
too, a resonant perturbation can only appreciably impact the spin evolution
during the radical pair lifetime, $\tau$, if -$\gamma
B_{1}/(2\pi)\tau\gtrapprox 1$, i.e. if the lifetime is sufficiently large to
permit at least one Larmor precession in the field associated with the
perturbation. As here $k_{f}^{-1}=1\mskip 3.0mu\mu$s and $k_{b}^{-1}=1\mskip
3.0mu\mu$s, we can expect effects for $B_{1}$ exceeding 36$\mskip 3.0mu\mu$T
to 72$\mskip 3.0mu\mu$T to become effective, which is in tentative agreement
with the observed $B_{1}$-dependence. On the other hand, driving fields larger
than $0.1\mskip 3.0mu$mT allow very significant reductions of the
recombination yield. It is interesting to note that for controlling reaction
outcomes in the geomagnetic field ($B_{0}=50\mskip 3.0mu\mu$T) slightly larger
control fields exceeding $0.12\mskip 3.0mu$mT are required.

Above we have assumed that controls are applied with a time-resolution of
$\Delta t=1\mskip 3.0mu$ns. While AWGs are available with sufficient sample
rates and analog bandwidths (e.g. $50\mskip 3.0mu$GS/s at $15\mskip 3.0mu$GHz
bandwidth) and the delivery of broadband magnetic fields at the required
intensity has been realized (e.g. for high-field ENDOR spectroscopy [36]), the
question of the minimal bandwidth to control radical pair reactions is a
pertinent one. To this end, we have studied the dependence of the control
result on the “speed limit” enforced on the control signal. Specifically, we
stipulated that the controls must obey
$|{u}_{n}-{u}_{n-1}|<\dot{u}_{max}\Delta t$. For $B_{0}=1\mskip 3.0mu$mT, Fig.
6c summarizes the dependence of the singlet yield on this speed limit. We find
that transition times of $\dot{u}_{max}^{-1}\lessapprox 8\mskip 3.0mu$ns are
required to elicit significant control; fast controls only deliver marginal
additional gains. The minimal speed limit corresponds to a bandwidth of
approximately $\frac{\dot{u}_{max}}{2}=63\mskip 3.0mu$MHz, which is of the
order of, but smaller than, the magnitude of the spread of energy eigenvalues
of the radical pair ($93\mskip 3.0mu$MHz).

###  III.3 Protocol applied to a realistic exciplex-forming radical pair

Figure 7: Recombination yield minimization of the 18 spin PY/DCB system in a
bias field of $1\mskip 3.0mu$mT. (a) shows the singlet probability, scaled to
remove the spin-independent decay due to radical recombination, i.e. $p_{\rm
S}(t)\exp(kt)$, with and without control. Control amplitudes were optimized
using $32$ samples and $1024$ time-steps of $1\mskip 3.0mu$ns; a spectrogram
of the control field magnitude is shown as insert. (b) shows the final yield
as realized using L-BFGS-B complemented gradient optimization for increasing
sampling intervals and various block sizes, $N/B$, as indicated in the legend.

We now present results for coherently controlling the symmetric recombination
of a large spin system consisting of 18 spins. We consider the pyrene/para-
dicyanobenzene (PY/DCB) exciplex-forming donor-acceptor systems in a solvent
of moderate polarity [37, 38, 16]. The magnetosensitive, delayed fluorescence
emission of this system has been widely studied. To make this large scale
simulation tractable, we limit ourselves to the scenario of homogeneous decay,
i.e. $k_{\rm S}=k_{\rm T}=k$, and neglect inter-radical interactions. While
this may appear crude at first glance, models typical of this have seen
considerable success in the interpretation of experimental data [39]. The
approximations are adequate since for the majority of its lifetime, the
radical pair is in fact diffusively well separated, and the moderate reaction
asymmetry of the real system (resulting from different electron transfer rate
constants in the singlet and triplet state; the latter to the locally excited
triplet state) has all but a small effect on the magnetic field sensitivity.
Thus, in this limit, the singlet probability can be expressed as

$\displaystyle p_{\rm
S}(t)=e^{-kt}\left(\frac{1}{4}+\sum_{\alpha\beta\in\\{x,y,z\\}}R^{(1)}_{\alpha,\beta}(t)R^{(2)}_{\alpha\beta}(t)\right),$
(29)

where

$\displaystyle R^{(i)}_{\alpha\beta}(t)=\frac{1}{Z_{i}}{\rm
Tr}\left[\hat{S}_{i\alpha}(0)\hat{S}_{i\beta}(t)\right]$
$\displaystyle=\frac{1}{Z_{i}}{\rm
Tr}\left[\hat{S}_{i\alpha}\hat{U}^{\dagger}(t,0)\hat{S}_{i\beta}\hat{U}(t,0)\right]$
(30)

are the electron spin correlation tensors for each spin in radical $i$, which
are defined or calculable in the Hilbert space of the individual radicals.
Here, $\hat{U}$ accounts for the coherent evolution only, i.e.
$\hat{A}(T)=\hat{H}(T)$ (cf. eq. (6)). Thus, control of the radical pair
reaction yield can be realized by separately calculating the gradients for
each of the spin correlation tensors of two recombining radicals, and
assembling the gradient of the reaction yield using the chain rule. Again, we
follow the ideas of blocking the time axis and sparse sampling to increase the
efficiency. In addition to the Trotter-Suzuki splitting discussed above, we
use the fact that the Hamiltonian is only a function of the total nuclear spin
associated with groups of completely equivalent spins, which permits an
efficient calculation in the “coupled representation” [40]. Fig. 7 shows the
minimization of the PY/DCB singlet recombination yield for $k^{-1}=200\mskip
3.0mu$ns, a control amplitude of $0.2\mskip 3.0mu$mT, and a static basing
field of $10\mskip 3.0mu$mT. The control was extended to the first $1024$ time
steps of $1\mskip 3.0mu$ns each. Varying block and sample sizes, again
confirms practically sufficient fidelity from non-complete optimization, as is
demonstrated in Fig. 7b (also see Fig. S3). In our view, exciplex-systems of
the type discussed here would be ideally suited to realize first proof-of-
principle experimental realizations of OQC of radical pair reaction yields
[16, 37, 38, 18, 9] by building on the time-resolved MFE measurements and
studies of radio-frequency effects that have already been realized.

###  III.4 Applications to Quantum Biology

Many observations of biochemical magnetosensitivity appear to be attributed to
radical pair intermediates in complex biophysical environments [21, 30]. We
anticipate that the control of radical pair reactions will develop into a
versatile applied field, particularly in the spin-biological context. Static
and oscillatory magnetic fields can and have been used to perturb such
systems, but these naïve controls are neither selective for specific pathways
nor can reaction yields be changed drastically. But our OQC approach promises
distinguishability of reaction pathways involving different radical pairs and
enlarged effect sizes in weak biasing fields, such as the geomagnetic field.
We demonstrate this potential here in terms of two selected examples with a
quantum biology underpinning.

Figure 8: Controlled and unperturbed singlet recombination yield of a flavin
semiquinone/superoxide radical pair, abstracted as FADH∙/Z∙, initialized in
the triplet state, as a function of the applied static biasing field. The
radical pair was assumed to undergo a spin-independent decay/escape reaction
with rate constant $k_{f}=1\mskip 3.0mu\mu$s-1 and spin-selective
recombination in the singlet state with rate constant $k_{b}=1\mskip
3.0mu\mu$s-1. The first $5\mskip 3.0mu\mu$s after radical pair generation were
controlled using a step-wise constant control field with $1\mskip 3.0mu$ns
time resolution and a maximal amplitude $0.25\mskip 3.0mu$mT, applied
perpendicular to the static magnetic field. The control was optimized using
the L-BFGS-B algorithm applied to $5$ blocks of $1\mskip 3.0mu\mu$s, sampling
the yield every $25\mskip 3.0mu$ns. Optimizations were repeated $5$ to $6$
times starting from random initial controls; with the error bars indicating
the $75$-percentile range of best results (the variability between runs is
noticeably small compared to the change in recombination yield produced
through control). In b), we show the time-dependence of the singlet
probability of the radical pair in a $4\mskip 3.0mu$mT biasing field in
absence and presence of the optimized control field, scaled to remove the
spin-independent decay, i.e. $p_{\rm S}(t)\exp(k_{f}t)$. The associated
singlet recombination yields amount to $0.218$ and $0.115$ with and without
controls, respectively.

The flavin semiquinone (FADH∙)/superoxide (O${}_{2}^{\bullet-}$) radical pair
has been suggested to underpin various biomagnetic processes, as implicated by
MFEs in neurogenesis [41] and cellular energetics [42]. Being of so-called
reference-probe type, this radical pair is renowned for its sensitivity to
weak magnetic fields [43], but questions concerning its potentially fast spin
relaxation remain [44]. We consider the limit of slow spin relaxation,
conventionally referred to as the FADH∙/Z∙-model. Here, Z∙ denotes a radical
that is devoid of hyperfine interactions and slowly relaxing, i.e. an
idealization of O${}_{2}^{\bullet-}$, as it might be approached in the real
system through immobilization. Fig. 8 shows the ability of our approach to
maximize the recombination yield of the triplet-born radical pair. This model,
with parameters taken from [45], exhibits a large low-field effect in response
to static magnetic field (orange line). Controlling the recombination (blue
line) allows to boost the recombination efficiency, realizing recombination
yields that are not achievable by the application of static magnetic fields of
any strength. In a biophysical context, this would lead to a decrease of
superoxide release from flavin re-oxidation reactions, thereby reducing the
concentration of free superoxide. Superoxide is a central reactive oxygen
species that assumes important roles as a signalling molecule, thus linking
radical pair reactions to central biological control pathways [46]. Since
magnetosensitivity in a physiological setting is currently not well understood
beyond circumstantial evidence, the latter is a more distant aim, albeit an
auspicious one, because it could allow targeting specific redox pathways in
isolation for various applications, e.g. wound healing [47] or the suppression
of the adverse effects of hypomagnetic field [48, 41] exposure for space
travel. Other future applications might be the control of chemical reactions
for synthetic purposes, e.g. in the context of flavin photocatalysis [49].

In sensory biology, a certain consensus has been reached attributing an
inclination compass in birds and various other species to the proposed radical
pair mechanism involving the flavo-enzyme cryptochrome [21]. However, the
identity of the actual magnetosensitive radical pair is still subject to
debate. _In vitro_ experiments on cryptochromes have demonstrated
magnetosensitivity of a photo-induced radical pair comprising the flavin anion
radical (FAD∙-) and a surface exposed tryptophan cation radical (W∙+),
implicating FAD∙-/W∙+ as the sensory species [50]. However, magnetoreception
appears to also be possible in the dark, suggesting the potential involvement
of a re-oxidation pathway, implicating once more the flavin
semiquinone/superoxide radical pair (FADH∙/O${}_{2}^{\bullet-}$) in
magnetoreception. Furthermore, various alternative radical pairs have been
suggested in theoretical studies demonstrating superior magnetic field
sensitivity, e.g. for ascorbyl and Z∙ containing radical pairs [43]. Given the
multitude of suggestions and hypothesis, the distinguishability of various
radical pairs in spin biology is a timely question. This question cannot be
resolved by direct spectroscopic means as the radical pair is a transient
intermediate in a complex biophysical environment in a living organism. We
here demonstrate that OQC could in principle provide the sought insights.

Figure 9: Change in singlet yield contrast $\Delta_{S}$ for FAD∙-/W∙+ and
FADH∙/Z∙ radical pairs subject to controls optimized to maximise or minimise
the singlet recombination yield contrast $\Delta_{S}$ in the geomagnetic field
for one of the pairs (as indicated in the axis labels). The geomagnetic field
was assumed as $50\mskip 3.0mu\mu$T and the radical pair lifetime $1\mskip
3.0mu\mu$s, i.e. $k_{S}=k_{T}=1\mskip 3.0mu\mu$s-1. The control was applied
parallel to the static field with an amplitude of $12.5\mskip 3.0mu\mu$T. The
first $5\mskip 3.0mu\mu$s after radical pair generation in the singlet state
were controlled with a time-step of $1\mskip 3.0mu$ns. Yields were calculated
by integrating to $t_{\text{max}}=14\mskip 3.0mu\mu$s. The optimization was
realized through $40$ steps of steepest descent/ascent. The unperturbed
singlet yield contrasts are $\Delta_{S}=0.017$ and $0.13$ for FAD∙-/W∙+ and
FADH∙/Z∙, respectively. In b), we show the accumulated reaction yield
differences for each,
$k_{S}\int_{0}^{t}p_{S}(n_{\text{max}},t)-p_{S}(n_{\text{min}},t)\,dt$, i.e.
the accumulated yield difference up to time $t$ as a function of time for the
FAD∙-/W∙+ (top panel) and FADH∙/Z∙ (bottom panel) radical pair. The red lines
for each radical pair correspond to the case of the radical pair dynamics
being perturbed by the OQC-derived perturbation targeting the minimization of
$\Delta_{S}$ for FAD∙-/W∙+. The FADH∙/Z∙ pair is practically unresponsive to
this perturbation (bottom panel), while FAD∙-/W∙+’s directional sensitivity is
mostly attenuated. For FAD∙-/W∙+, the grey line applies to the scenario of
applying a harmonic perturbation with frequency $1.4\mskip 3.0mu$MHz,
demonstrating that designed perturbations are significantly more disruptive
than unspecific perturbations.

We considered models of the FAD∙-/W∙+ and FADH∙/Z∙ radical pairs, including
the anisotropic hyperfine interactions of the dominant nuclear spins in a
relative orientation as found in a cryptochrome (parameters taken from [51]).
The recombination yield of this system depends on the orientation of the
geomagnetic field relative to the protein frame. The difference of the maximal
and minimal yield, realized for field directions ${\bf n}_{\text{max}}$ and
${\bf n}_{\text{min}}$, i.e. $\Delta_{S}=Y_{S}({\bf
n}_{\text{max}})-Y_{S}({\bf n}_{\text{min}})$ is considered a measure of the
sensitivity of the radical pair-based compass. In Fig. 9, we demonstrate the
ability of OQC to control, i.e. maximize or minimize, this performance
parameter for the two radical pairs considered. The control fields were
applied parallel to the geomagnetic field. Excitation profiles were optimized
for one radical pair and its performance subsequently evaluated in both. We
find that OQC can significantly alter, enhance or attenuate, $\Delta_{S}$ for
both radical pairs. Importantly, sequences optimized for one pair turn out to
only have a minor effect of the spin dynamics of the other. Thus, the control
fields are selective, permitting to distinguish the two models insofar as they
induce a large effect on the singlet recombination yield contrast in the
targeted pair and a negligible effect in the competitor. In particular, the
controls derived for FAD∙-/W∙+ appear as an apt tool to realize such an
experiment, because the contrast $\Delta_{S}$ is strongly altered, the system
specificity high, the time of radical pair generation is controllable using
pulsed photo-excitation, and the parameters pertaining to FAD∙-/W∙+ are best
known, as they are available from _in vitro_ experiments and derivable from
crystal structures and computational chemistry methods. Clearly, here we limit
ourselves to only demonstrate principal feasibility; more detailed models will
have to be studied before instigating actual experiments.

The robustness of our derived control pulses to small changes is central for
future applications. While control engineering generally aims to capture the
intricacies of the control target as closely as possible, in practice our
knowledge of the actual spin Hamiltonian parameters [52] might be limited, or
the available computational resources might restrict the optimization to a
simplified spin system comprising only the strongest hyperfine-coupled nuclear
spins. In this context, controls derived from simplified systems, in
particular those not accounting for weakly coupled nuclear spins or spin
relaxation processes, are crucial. As evidenced in the Supplemental Material,
controls derived using our approach are indeed robust when applied to
simulations accounting for the relevant inhomogeneous hyperfine interactions
via the Schulten-Wolynes framework [53] (cf. Fig. S5). The controls are also
found to be robust against noise in the form of random field fluctuations,
which induce decoherence and spin relaxation (cf. Fig. S6). We particularly
find that controls derived for idealized systems retain their unique property
of permitting recombination yields that cannot be realized using static
magnetic field of any flux density even in the presence of spin relaxation and
unresolved hyperfine interactions.

##  IV Discussion

In electronic and nuclear spins of organic radical pairs, coherence lifetimes
on the order of $100\mskip 3.0mu$ns to microseconds are typical, even in
complex biophysical environments [54], opening a window of opportunity to
realize control by applying time-dependent magnetic fields to manipulate
reaction outcomes. Yet, the approach is still in its infancy. Experimental
approaches have been limited to perturbing radical pair spin systems by radio-
frequency magnetic fields [17, 23]. While positioned far from optimal control,
these experiments are crucial as they demonstrate the principal possibility of
influencing radical pair dynamics by RF-magnetic fields in weak magnetic
fields. Examples of RF-magnetic field effects are furthermore well-established
for spin-biological radical pairs _in vivo_ [17], such as those supposedly
underpinning the avian compass sense. The latter are particularly noteworthy,
because surprisingly weak RF-magnetic fields have been found to interrupt
magnetoreception, suggesting an exquisitely sensitive underlying system of
remarkable coherence time.

Optimal quantum control protocols could be used to drive chemical and spin-
biological systems in a controlled and selective way, whereby reaction
outcomes not attainable by static magnetic field perturbations alone become
accessible. For quantum information processing in near term quantum hardware,
among the many hurdles faced in counteracting environment-induced decoherence
is the need to realize precise timing of quantum gate operations [55, 56] in
nanosecond timescales. To evade this bottleneck, precise or adaptive [57, 4,
58] quantum control protocols may be harnessed to change the frequency of an
applied control field dynamically throughout its pulse duration. But this is
often hard to realize in practice, since formulating the appropriate optimal
quantum control protocols is often highly non-trivial. In nature, radical
pairs seem to evade the effects of decoherence to retain exceptionally long
coherence times, resulting in exquisite and unparalleled magnetometric
sensitivity even in extremely noisy biological settings [59]. Thus, optimal
quantum control protocols, formulated with the primary aim of reproducing this
sensitivity in spin-chemical reactions in solution, could also inform methods
for similarly realizing artificial or bio-molecular [60, 17] qubits highly
sensitive to applied magnetic fields. This could translate to further
applications [61, 55, 62, 63] in both quantum metrology and quantum
information processing. In particular, potential engineering principles
derived from gleaning deeper physical insights gained from artificially
driving radical pair reactions could conceivably inspire a new generation of
quantum enhanced sensing technologies for magnetometry boasting heightened
sensitivity and robustness against noise [61, 64, 7, 65, 63, 66, 67]. But to
realize this goal, and future applications of optimal quantum control in a
potentially biophysical context, efficient approaches for control design are
vital.

The control of radical pair reactions is distinct from the majority of OQC
approaches in the sense that the optimization target is not defined at a
particular instant of time, but encoded in the time-evolution from radical-
pair creation to its eventual decay, $t\in[0,\infty)$, i.e. the reaction
yield. Here, we have generalized the GRAPE-paradigm of using piece-wise
constant control amplitudes to realize reaction control. The suggested
algorithm aims to overcome limitations of previous approaches by targeting
actual reaction yields (rather than continuously updating the control to
optimize a fidelity function involving observables that commute with the
Hamiltonian or the singlet probability at one instance of time) and permitting
control in weak magnetic fields. A time-global update scheme of control
parameters (as opposed to the time-local approach from [25, 8, 24]), central
to GRAPE, has here been found essential to realize the latter efficiently. The
computational cost of controlling radical pair reaction yields by time-global
optimization techniques such as GRAPE has been viewed as prohibitive [8]. And
indeed, significant computational costs are associated with the evaluation of
the numerous single time-step propagators and their Fréchet derivatives, and,
for the optimization of reaction yields in particular, the calculation of the
time-integral recombination yield and its gradient through iterated matrix
multiplications involving the former. We have introduced a block-optimization
scheme and sparse sampling of the reaction yield to improve the efficiency of
this step while retaining adequate control fidelity. To facilitate the
efficient evaluation of propagators and their derivatives, we have utilized an
iterated Trotter-Suzuki splitting. We further leverage curvature information
of the loss function in lieu of standard steepest decent/ascent approaches to
reduce the number of optimization steps required, as previously suggested
[31].

The additional algorithmic modifications we incorporate into our GRAPE-
inspired approach are vital for addressing the high computational costs, which
had impeded earlier attempts to realize GRAPE-based reaction control of
radical pair dynamics. We note that while sparse sampling of the reaction
yield demonstrably provides a useful and computationally feasible proxy to the
actual reaction yield, sampling of the singlet probability at a single moment
in time is oftentimes misleading the optimization in weak magnetic fields as
the increased/decreased probability at the chosen time does not necessarily
correlate with increased/decreased reaction yield. On the other hand, as we
have shown, sparse sampling at a moderate sampling density reduces the
computational cost while yielding close to optimal optimization outcomes and a
smooth progression for $B>1$. We have further demonstrated the applicability
of the approach as implemented here to realistically complex radical pair
systems, such as the widely studied PY/DCB exciplex system subject to
symmetric recombination. Our suggested approach can serve as a blueprint for
future developments and optimizations. Although not explored here, further
algorithmic modification could be attempted, such as bottom-up strategies for
which the sampling density is increased as the optimization progresses.
Additionally, the optimization blocks could be overlapped. Nonetheless, as our
approach yielded control parameters that upon complete optimization did not
change significantly subject to the utilized termination criteria of the
optimizer, these additional improvements did not appear imminently pertinent.
Although, they could further speed up the rate at which the target is reached.
Additional computational efficiency could furthermore be realized by
reformulating the singlet yield evaluation in terms of the propagation of
wave-packets, in particular in combination with using an incomplete basis for
the nuclear spin functions and trace sampling [68]. We are currently
developing such approaches with the aim of extending control to the complete
open systems [69] description of radical pair systems, including via the
incorporation of direct control of the dissipation processes through
incoherent control. As for the calculation of propagators and derivatives,
higher order operator splitting approaches could be used and the order of the
approach, i.e. accuracy, to the optimization process increased as the target
function converges, as has been successfully demonstrated in [58] for a state
transfer problem. However, since the generation of propagators and derivatives
was not identified as the most time-consuming step for the considered problem
sizes, which is related to the evaluation of reaction yields and gradients
instead, here we have again refrained from such additional optimizations. From
the reaction-control point of view, it is noteworthy that control has allowed
us to suppress the singlet recombination to levels below those achievable by
application of static magnetic fields of any strength. Furthermore, the
achievable recombination yields are well below the $25\mskip 3.0mu$% yield
that is expected to ensue for complete randomization of the spin states. This
shows that the control fields here do not merely enhance spin relaxation, as
one might expect, but drive the system to optimal reaction outcomes.

Finally, we comment on the alternative approach of time-local optimization as
introduced in [8]. Even with the described optimization tweaks, the presented
GRAPE-based approach to radical pair reaction control remains computationally
more demanding than the time-local approach. However, it captures the time-
integrated nature of the reaction yield, the observable quantity, rather than
merely the singlet probability at one moment in time, which is not necessarily
a well-defined proxy of the overall recombination yield. Another limitation of
the time-local approach is conceptual and related to the property of
permitting only controls that at every moment of time increase the
optimization target function. This property permits the straight-forward
derivation of the optimal control amplitudes but constrains the optimizations
by excluding controls that would permit a better overall optimization at the
cost of a momentary reduction in optimality. We argue that this constraint is
critical in the optimization of weak biasing fields. In the Supporting
Material we provide a generalized version of the time-local algorithm from [8]
that permits the optimization of the singlet probability, still for one
instance in time, in the presence of spin selective recombination, as assumed
here. Comparing this generalized algorithm to the GRAPE-based approach, we
find that it is inferior for optimizing the recombination yield in weak
biasing fields, generally requiring larger control amplitudes while still
falling short in terms of optimization fidelities. Please refer to the
Supporting Material for the detailed quantitative comparison undertaken for
the 7-spin model system. As shown there, superiority of the GRAPE-based
approach is realized not only when optimizing the singlet probability at a
single time, but also when applying block-wise approach with successively
redefined “moving target”. While this time-blocking approach appears
auspicious for the time-local approach too, the GRAPE-based approach delivers
markedly better yield optimization in weak biasing fields.

The quantum control of radical pair reactions is still in its infancy. Our
present work focuses on the control of the coherent evolution via Zeeman
interactions, as this is a modality in principle addressable with currently
available experimental tools without complex changes to the chemistry. We
anticipate that optimal control of radical pair reactions will find
applications in the context of spin-biological radical pair reactions, either
as an analytical tool or to drive the systems towards desired outcomes. For
example, for radical pair reactions implicated in magnetoreception, a key
problem to be addressed concerns distinguishing between different incarnations
of radical reactions to better identify the underlying mechanism in a complex
biological environment [70, 65, 71]. In this scenario, quantum control could
allow us to derive control parameters that suppress magnetic field effects in
one mechanism while leaving those of competing alternative models unaltered.
Ultimately, controls could allow us to turn MFEs into a tool to potentially
control reaction processes of physiological relevance for applications in
biology and medicine, by facilitating selective stimulation or suppression of
cellular functions linked to radical pair processes. We hope that this work,
by building on established techniques from quantum control like GRAPE to go
beyond existing local optimization techniques, will lay the groundwork for
realizing this.

##  V Conclusions

We have suggested and implemented a GRAPE-inspired approach to control
reaction yields of radical pair recombination reactions through synchronized
application of radio-frequency magnetic fields in a direction perpendicular to
a static biasing field. We deviate from previous approaches by optimizing
reaction yields, i.e. the actual experimental observables, instead of fidelity
measures defined at a particular time instants. Computational efficiency has
been realized by block optimization, sparse sampling, a split operator
approach for calculating propagators and by utilizing curvature information in
the non-linear optimization problem. The approach thus realized is applicable
to radical pair systems of realistic complexity. We have demonstrated that
control permits the reduction of the singlet recombination yield of a singlet-
born model system to $10\mskip 3.0mu$%, a yield unattainable by static
magnetic field application alone, when using moderate control amplitudes and
bandwidth (cf. Fig. 6). We hope that our approach helps accelerate the
adoption of control approaches for radical pair recombination reactions in low
external fields. We anticipate that these future applications will utilize
control beyond mere reaction control, for example to realize
distinguishability of radical pairs in complex environments and for
facilitating model selection for competing hypotheses based on experiments
with optimized contrast power.

###### Acknowledgements.

This work was supported by the Office of Naval Research (ONR Award Number
N62909-21-1-2018), the Leverhulme Trust (RPG-2020-261), and the Engineering
and Physical Sciences Research Council (EPSRC grants EP/V047175/1 and
EP/X027376/1). We acknowledge use of University of Exeter’s HPC facility.

## Appendix

###  V.1 Sparse Blocking

The control of reaction yields via numerical approaches like GRAPE is expected
to be costly. The computational demand is particularly large if, as we do here
(and unlike [8, 24]), the goal is to optimize recombination yields of radical
pair reactions, which are time-integral quantities (cf. eq. (7); unlike the
typical fidelity measures in quantum control applications pertaining to the
implementation of a state transfer or a unitary operation for a fixed time
$t$). The approach requires iteratively repeated matrix multiplications for
evaluating the time-dependent singlet-probabilities that sum to the
recombination yield (and its associated gradients).

Formally, eq. (26) implies

$\displaystyle\sum_{n=1}^{N}[(n-1)+3]+nL[(n-1)+3]=O(N^{3})$ (A1)

matrix multiplications, i.e. for every $t_{n}$, $(n-1)$ matrix multiplications
are required to evaluate $\hat{U}$, $(n-1)$ to evaluate
$\hat{U}_{+}\hat{W}_{n,l}\hat{U}_{-}$ (via $L_{n,l}$) for every $l$ and $n$;
and an additional three each to evaluate the expectation value/gradients. In
practical implementations, the estimate in eq. (A1) can be improved upon by
preserving terms from previous $n$, i.e. $n-1$, allowing the update of
$\hat{U}$ and $\hat{U}_{+}\hat{W}_{n,l}\hat{U}_{-}$ with a singular matrix
multiplication (instead of $n-1$), thus requiring

$\displaystyle\sum_{n=1}^{N}4(1+nL)=2LN^{2}+2(L+2)N$ (A2)

matrix multiplications. Nevertheless, due to the quadratic scaling of the
evaluation of eq. (27) with $N$, the gradient/expectation value evaluation
will often exceed the cost of evaluating the elementary propagators and their
derivatives (scaling as $N$) for large $N$ and thus be the limiting factor.

To efficiently optimize the reaction yield of complex radical pairs, we
suggest the adaptation of a gradient-based optimization procedure. We posit
that the optimization can be realized in terms of $B$ disjoint blocks of $N/B$
time-steps each (time blocking), with the initial state taken as that of the
system subject to the preceding optimized controls, sparsely sampled using
$N/(BS)$ samples per block, as illustrated in Fig. 3. Here, $B$ and $S$ are
integers that divide $N$ and $B$, respectively. The number of matrix
multiplications in this approach is reduced to

$\displaystyle
B\left(\sum_{n=1}^{N/B}(1+nL)[1+3(\delta_{n\;\mathrm{mod}\;s,0})]\right)$
$\displaystyle=B\left(\sum_{n=1}^{N/B}(1+nL)+3\sum_{m=1}^{N/BS}(1+LS)\right)$
$\displaystyle=\frac{N^{2}L(S+3)}{2BS}+N\left(2L+1+\frac{3}{S}\right).$ (A3)

The number of matrix multiplications necessary in general still scales
quadratically with $N$, i.e. as $O(N{{}^{2}}/B)$. Linear scaling is only
obtained in the limit that $B=N$, i.e. time-local optimization. But in
practice, the block size can be chosen much smaller than $N$, i.e. $N/B\ll N$,
to realize a substantial speedup, whilst delivering adequate practical
optimization of the reaction yield. Furthermore, for $S\gg 3$, the complexity
of the algorithm approaches independence of $S$.

###  V.2 Trotter-Suzuki Splitting

To efficiently calculate the propagators and the associates Fréchet
derivatives, we use a scaling and squaring approach in combination with a
Trotter-Suzuki splitting of the elementary propagator. Specifically, we
compute the matrix exponent

$\displaystyle e^{X}=\left(e^{{2^{-s}X}}\right)^{2^{s}}$ (A4)

by $s\in\mathbb{N}$ repeated squaring starting from $e^{{2^{-s}X}}$ calculated
as

$\displaystyle e^{{2^{-s}X}}=e^{{-i2^{-s-1}}\Delta
t\hat{A}^{\prime}}e^{{-i2^{-s}}\Delta
t\hat{A}^{\prime\prime}}e^{{-i2^{-s-1}}\Delta
t\hat{A}^{\prime}}+O\left(\frac{\Delta t^{2}}{2^{2s}}\right),$

with

$\displaystyle X=-i\hat{A}^{\prime}dt-i\hat{A}^{\prime\prime}dt$ (A5)

being a convenient splitting of $X=-idt\hat{A}_{i}$ (cf. eq. (19).
Practically, $\hat{A}^{\prime}$ and $\hat{A}^{\prime\prime}$ represent the
drift and control Hamiltonians, respectively. If the latter is based on
coupling applied magnetic fields via the Zeeman interaction, the propagator
under $\hat{A}^{\prime\prime}$ can be expressed as a direct product of
evolution operators for the individual electron spins based on expression

$\displaystyle e^{-i\omega\hat{\rm n}\cdot
S}=\cos\left(\frac{\omega}{2}\right)I-2i\sin\left(\frac{\omega}{2}\right)\hat{\rm
n}\cdot S$ (A6)

where $\omega\in\mathbb{N}$ is a parameter related to the strength of the
control field and $n$ in a unit vector specifying the orientation of the
perturbation in the laboratory frame.

Likewise, the Fréchet derivative of the exponential $e^{X}$ can be evaluated
using an $s$-fold recurrence based on [72],

$\displaystyle L(X,E)=e^{X/2}L(X/2,E/2)+L(X/2,E/2)e^{X/2},$ (A7)

obtained by differentiating $e^{X}=({e^{X/2})}^{2}$ using the chain rule. In
combination with eq. (A4), the only derivatives that are explicitly required
in this process are those of $e^{{-i2^{-s}}\Delta tA^{\prime\prime}}$. For
Zeeman coupling with a single spin propagator as given by eq. (A6), the
necessary Fréchet derivatives can be obtained in analytic form, in combination
with the chain rule, from

$\displaystyle L(\omega\hat{\rm n}\cdot S,-i\hat{\rm n}\cdot
S)=-i\cos\left(\frac{\omega}{2}\right)\hat{\rm n}\cdot
S-\frac{1}{2}\sin\left(\frac{\omega}{2}\right)\mathbb{1}.$ (A8)

We optimize using the limited memory variant of the
Broyden–Fletcher–Goldfarb–Shanno (L-BFGS-B) algorithm, that preconditions
gradients using curvature information, and use the SciPy [73] implementation
with a stopping criterion of $10^{-6}$ and $10^{-7}$ for the change of the
sparsely sampled yield and the projected gradient, respectively. The maximum
number of iterations was set to 200. We further ran simulations in parallel
for varying number of blocks $B$ and the sampling interval $S$, as seen in
Figs. 5(a) and 5(b), to assess the performance and scalability of our
approach.

## References

  * Koch et al. [2022] Koch et al., Quantum Optimal Control in Quantum Technologies. Strategic report on current status, visions and goals for research in Europe, EPJ Quantum Tech. 9, 19 (2022).
  * Tannor  _et al._ [1986] D. J. Tannor, R. Kosloff, and S. A. Rice, Coherent pulse sequence induced control of selectivity of reactions: Exact quantum mechanical calculations, J. Chem. Phys. 85, 5805 (1986).
  * Tannor and Rice [1985] D. J. Tannor and S. A. Rice, Control of selectivity of chemical reaction via control of wave packet evolution, J. Chem. Phys. 83, 5013 (1985).
  * Brixner and Gerber [2003] T. Brixner and G. Gerber, Quantum control of gas-phase and liquid-phase femtochemistry, ChemPhysChem 4, 418 (2003).
  * Dowling and Milburn [2003] J. P. Dowling and G. J. Milburn, Quantum Technology: The Second Quantum Revolution, Philos. Trans. Royal Soc. A 361, 1655 (2003).
  * Guerreschi  _et al._ [2013] G. G. Guerreschi, M. Tiersch, U. E. Steiner, and H. J. Briegel, Optical switching of radical pair conformation enhances magnetic sensitivity, Chem. Phys. Lett. 572, 106 (2013).
  * Vitalis and Kominis [2017] K. M. Vitalis and I. K. Kominis, Quantum-limited biochemical magnetometers designed using the fisher information and quantum reaction control, Phys. Rev. A 95, 032129 (2017).
  * Masuzawa  _et al._ [2020] K. Masuzawa, M. Sato, M. Sugawara, and K. Maeda, Quantum Control of Radical Pair Reactions by Local Optimization Theory, J. Chem. Phys. 152, 014301 (2020).
  * Mims  _et al._ [2021] D. Mims, J. Herpich, N. N. Lukzen, U. E. Steiner, and C. Lambert, Readout of spin quantum beats in a charge-separated radical pair by pump-push spectroscopy., Science 374, 1470 (2021).
  * Pophof et al. [2023a] Pophof et al., Biological effects of electric, magnetic, and electromagnetic fields from 0 to 100 mhz on fauna and flora: Workshop report, Health Phys. 124, 39 (2023a).
  * Pophof et al. [2023b] Pophof et al., Biological effects of radiofrequency electromagnetic fields above 100 mhz on fauna and flora: Workshop report, Health Phys. 124, 31 (2023b).
  * Wasielewski  _et al._ [1983] M. R. Wasielewski, C. H. Bock, M. K. Bowman, and J. R. Norris, Controlling the duration of photosynthetic charge separation with microwave radiation, Nature 303, 520 (1983).
  * Henbest  _et al._ [2004] K. B. Henbest, P. Kukura, C. T. Rodgers, P. J. Hore, and C. R. Timmel, Radio frequency magnetic field effects on a radical recombination reaction: A diagnostic test for the radical pair mechanism, J. Amer. Chem. Soc. 126, 8102 (2004), pMID: 15225036.
  * Rodgers et al. [2009] Rodgers et al., Radiofrequency polarization effects in zero-field electron paramagnetic resonance, Phys. Chem. Chem. Phys. 11, 6569 (2009).
  * Wedge  _et al._ [2009] C. Wedge, C. Rodgers, S. Norman, N. Baker, K. Maeda, K. Henbest, C. Timmel, and P. Hore, Radiofrequency polarization effects in low-field electron paramagnetic resonance, Phys. Chem. Chem. Phys. 11, 6573 (2009).
  * Grampp  _et al._ [2002] G. Grampp, M. Justinek, and S. Landgraf, Magnetic field effects on the pyrene—dicyanobenzene system: determination of electron self-exchange rates by mary spectroscopy, Mol. Phys. 100, 1063 (2002).
  * Mani [2022] T. Mani, Molecular qubits based on photogenerated spin-correlated radical pairs for quantum sensing, Chem. Phys. Rev. 3, 021301 (2022).
  * Rugg  _et al._ [2019] B. K. Rugg, M. D. Krzyaniak, B. T. Phelan, M. A. Ratner, R. M. Young, and M. R. Wasielewski, Photodriven quantum teleportation of an electron spin state in a covalent donor-acceptor-radical system, Nat. Chem. 11, 981 (2019).
  * Wedge  _et al._ [2013] C. J. Wedge, J. C. S. Lau, K.-A. Ferguson, S. A. Norman, P. J. Hore, and C. R. Timmel, Spin-locking in low-frequency reaction yield detected magnetic resonance, Phys. Chem. Chem. Phys. 15, 16043 (2013).
  * Kattnig [2021] D. R. Kattnig, F-cluster: Reaction-induced spin correlation in multi-radical systems, J. Chem. Phys. 154, 204105 (2021).
  * Hore and Mouritsen [2016] P. J. Hore and H. Mouritsen, The radical-pair mechanism of magnetoreception., Annual Rev. Biophys. 45, 299 (2016).
  * Ritz et al. [2009] Ritz et al., Magnetic compass of birds is based on a molecule with optimal directional sensitivity., Biophys. J. 96, 3451 (2009).
  * Engels et al. [2014] Engels et al., Anthropogenic electromagnetic noise disrupts magnetic compass orientation in a migratory bird., Nature 509, 353 (2014).
  * Tateno  _et al._ [2023] A. Tateno, K. Masuzawa, H. Nagashima, and K. Maeda, Anisotropic and coherent control of radical pairs by optimized rf fields, Int. J. Mol. Sci. 24, 10.3390/ijms24119700 (2023).
  * Sugawara [2003] M. Sugawara, General formulation of locally designed coherent control theory for quantum system, J. Chem. Phys. 118, 6784 (2003).
  * Khaneja  _et al._ [2005] N. Khaneja, T. Reiss, C. Kehlet, T. Schulte-Herbrüggen, and S. J. Glaser, Optimal control of coupled spin dynamics: design of NMR pulse sequences by gradient ascent algorithms, J. Magn. Reson. 172, 296 (2005).
  * Geng  _et al._ [2016] R. Geng, T. T. Daugherty, K. Do, H. M. Luong, and T. D. Nguyen, A review on organic spintronic materials and devices: I. magnetic field effect on organic light emitting diodes, J. Sci. Adv. Mater. Devices 1, 128 (2016), special Issue in Memory of Dr. P.E. BROMMER.
  * Lee  _et al._ [2011] H. Lee, N. Yang, and A. E. Cohen, Mapping nanomagnetic fields using a radical pair reaction, Nano Lett. 11, 5367 (2011).
  * Engmann  _et al._ [2023] S. Engmann, E. G. Bittle, and D. J. Gundlach, Magnetic field sensor based on a oled/organic photodetector stack, ACS Appl. Electron. Mater. 5, 4595 (2023).
  * Zadeh-Haghighi and Simon [2022] H. Zadeh-Haghighi and C. Simon, Magnetic field effects in biology from the perspective of the radical pair mechanism, J. R. Soc. Interface 19, 10.1098/rsif.2022.0325 (2022).
  * de Fouquieres  _et al._ [2011] P. de Fouquieres, S. Schirmer, S. Glaser, and I. Kuprov, Second order gradient ascent pulse engineering, J. Magn. Res. 212, 412 (2011).
  * Gauger  _et al._ [2011] E. M. Gauger, E. Rieper, J. J. L. Morton, S. C. Benjamin, and V. Vedral, Sustained quantum coherence and entanglement in the avian compass, Phys. Rev. Lett. 106, 040503 (2011).
  * Rabitz  _et al._ [2004] H. A. Rabitz, M. M. Hsieh, and C. M. Rosenthal, Quantum optimally controlled transition landscapes, Science 303, 1998 (2004).
  * Machnes et al. [2011] Machnes et al., Comparing, optimizing, and benchmarking quantum-control algorithms in a unifying programming framework, Phys. Rev. A 84, 022305 (2011).
  * Wilcox [1967] R. M. Wilcox, Exponential Operators and Parameter Differentiation in Quantum Physics, J. Math. Phys. 8, 962 (1967).
  * Hertel  _et al._ [2005] M. M. Hertel, V. P. Denysenkov, M. Bennati, and T. F. Prisner, Pulsed 180-ghz epr/endor/peldor spectroscopy, Magn. Res. Chem. 43, S248 (2005).
  * Kattnig  _et al._ [2011] D. R. Kattnig, A. Rosspeintner, and G. Grampp, Magnetic field effects on exciplex-forming systems: the effect on the locally excited fluorophore and its dependence on free energy, Phys. Chem. Chem. Phys. 13, 3446 (2011).
  * Hoang  _et al._ [2014] H. M. Hoang, T. B. V. Pham, G. Grampp, and D. R. Kattnig, Exciplexes versus loose ion pairs: How does the driving force impact the initial product ratio of photoinduced charge separation reactions?, J. Phys. Chem. Lett. 5, 3188 (2014).
  * Hayashi [2004] H. Hayashi,  _Introduction to Dynamic Spin Chemistry : Magnetic Field Effects on Chemical and Biochemical Reactions._ , World Scientific Lecture and Course Notes in Chemistry No. vol. 8 (World Scientific Publishing Company, 2004).
  * Atkins  _et al._ [2019] C. Atkins, K. Bajpai, J. Rumball, and D. R. Kattnig, On the optimal relative orientation of radicals in the cryptochrome magnetic compass, J. Chem. Phys. 151, 065103 (2019).
  * Ramsay and Kattnig [2022] J. Ramsay and D. R. Kattnig, Radical triads, not pairs, may explain effects of hypomagnetic fields on neurogenesis, PLoS Comput. Biol. 18, 1 (2022).
  * Usselman  _et al._ [2016] R. J. Usselman, C. Chavarriaga, P. R. Castello, M. Procopio, T. Ritz, E. A. Dratz, D. J. Singel, and C. F. Martino, The quantum biology of reactive oxygen species partitioning impacts cellular bioenergetics, Scientific Reports 6, 10.1038/srep38543 (2016).
  * Lee  _et al._ [2014] A. A. Lee, J. C. S. Lau, H. J. Hogben, T. Biskup, D. R. Kattnig, and P. J. Hore, Alternative radical pairs for cryptochrome-based magnetoreception, J. R. Soc. Interface 11, 20131063 (2014).
  * Hogben  _et al._ [2009] H. J. Hogben, O. Efimova, N. Wagner-Rundell, C. R. Timmel, and P. Hore, Possible involvement of superoxide and dioxygen with cryptochrome in avian magnetoreception: Origin of zeeman resonances observed by in vivo epr spectroscopy, Chem. Phys. Lett. 480, 118 (2009).
  * Deviers  _et al._ [2022] J. Deviers, F. Cailliez, B. Z. Gutiérrez, D. R. Kattnig, and A. de la Lande, Ab initio derivation of flavin hyperfine interactions for the protein magnetosensor cryptochrome, Phys. Chem. Chem. Phys. 24, 16784 (2022).
  * Winterbourn [2020] C. C. Winterbourn, Biological chemistry of superoxide radicals, ChemTexts 6, 10.1007/s40828-019-0101-8 (2020).
  * Huizen  _et al._ [2019] A. V. V. Huizen, J. M. Morton, L. J. Kinsey, D. G. V. Kannon, M. A. Saad, T. R. Birkholz, J. M. Czajka, J. Cyrus, F. S. Barnes, and W. S. Beane, Weak magnetic fields alter stem cell–mediated growth, Sci. Adv. 5, eaau7201 (2019), https://www.science.org/doi/pdf/10.1126/sciadv.aau7201 .
  * Zhang  _et al._ [2021] B. Zhang, L. Wang, A. Zhan, M. Wang, L. Tian, W. Guo, and Y. Pan, Long-term exposure to a hypomagnetic field attenuates adult hippocampal neurogenesis and cognition, Nat. Commun. 12, 1174 (2021).
  * König  _et al._ [2018] B. König, S. Kümmel, E. Svobodová, and R. Cibulka, Flavin photocatalysis, Phys. Sci. Rev. 3, 20170168 (2018).
  * Xu  _et al._ [2012] C. Xu, X. Yin, Y. Lv, C. Wu, Y. Zhang, and T. Song, A near-null magnetic field affects cryptochrome-related hypocotyl growth and flowering in arabidopsis, Adv. Space Res. 49, 834 (2012).
  * Babcock and Kattnig [2021] N. S. Babcock and D. R. Kattnig, Radical scavenging could answer the challenge posed by electron-electron dipolar interactions in the cryptochrome compass model, JACS Au 1, 2033 (2021).
  * Martino  _et al._ [2023] C. F. Martino, P. Jimenez, M. Goldfarb, and U. G. Abdulla, Optimization of parameters in coherent spin dynamics of radical pairs in quantum biology, PLOS ONE 18, 1 (2023).
  * Schulten and Wolynes [1978] K. Schulten and P. G. Wolynes, Semiclassical description of electron spin motion in radicals including the effect of electron hopping, J. Chem. Phys. 68, 3292 (1978), https://pubs.aip.org/aip/jcp/article-pdf/68/7/3292/11350283/3292_1_online.pdf .
  * Kattnig  _et al._ [2016] D. R. Kattnig, I. A. Solov’yov, and P. J. Hore, Electron spin relaxation in cryptochrome-based magnetoreception, Phys. Chem. Chem. Phys. 18, 12443 (2016).
  * Mao  _et al._ [2023] H. Mao, G. J. Pažėra, R. M. Young, M. D. Krzyaniak, and M. R. Wasielewski, Quantum gate operations on a spectrally addressable photogenerated molecular electron spin-qubit pair, J. Am. Chem. Soc. 145, 6585 (2023).
  * Xuereb  _et al._ [2023] J. Xuereb, P. Erker, F. Meier, M. T. Mitchison, and M. Huber, Impact of imperfect timekeeping on quantum control, Phys. Rev. Lett. 131, 160204 (2023).
  * Brixner and Gerber [2004] T. Brixner and G. Gerber, Adaptive quantum control of femtochemistry, Phys. Scr. 2004, 101 (2004).
  * Goodwin  _et al._ [2022] D. L. Goodwin, P. Singh, and M. Foroozandeh, Adaptive optimal control of entangled qubits, Sci. Adv. 8, eabq4244 (2022).
  * Smith  _et al._ [2024] L. D. Smith, J. Glatthard, F. T. Chowdhury, and D. R. Kattnig, On the optimality of the radical-pair quantum compass, arXiv:2401.02923 (2024).
  * Wasielewski et al. [2020] Wasielewski et al., Exploiting chemistry and molecular systems for quantum information science, Nature Rev. Chem. 4, 490 (2020).
  * Cai  _et al._ [2010] J. Cai, G. G. Guerreschi, and H. J. Briegel, Quantum control and entanglement in a chemical compass, Phys. Rev. Lett. 104, 220502 (2010).
  * Wu  _et al._ [2022] J.-Y. Wu, X.-Y. Hu, H.-Y. Zhu, R.-Q. Deng, and Q. Ai, A bionic compass based on multiradicals, J. Phys. Chem. B 126, 10327 (2022).
  * Li  _et al._ [2022] Z. Li, H. Yu, T. Shen, Y. Li, and W. Zhang, Bionic Magnetic Compass Algorithm Based on Radical Pair Theory, IEEE Sensors J. 22, 23812 (2022).
  * Aiello  _et al._ [2013] C. D. Aiello, M. Hirose, and P. Cappellaro, Composite-pulse magnetometry with a solid-state quantum sensor, Nat. Commun. 4, 1419 (2013).
  * Finkler and Dasari [2021] A. Finkler and D. Dasari, Quantum sensing and control of spin-state dynamics in the radical-pair mechanism, Phys. Rev. Applied 15, 034066 (2021).
  * Smith  _et al._ [2022] L. D. Smith, F. T. Chowdhury, I. Peasgood, N. Dawkins, and D. R. Kattnig, Driven radical motion enhances cryptochrome magnetoreception: Toward live quantum sensing, J. Phys. Chem. Lett. 13, 10500 (2022).
  * Winkler et al. [2023] Winkler et al., A review of the current state of magnetic force microscopy to unravel the magnetic properties of nanomaterials applied in biological systems and future directions for quantum technologies, Nanomaterials 13, 10.3390/nano13182585 (2023).
  * Fay  _et al._ [2021] T. P. Fay, L. P. Lindoy, and D. E. Manolopoulos, Spin relaxation in radical pairs from the stochastic schrödinger equation., J. Chem. Phys. 154, 084121 (2021).
  * Kallush  _et al._ [2022] S. Kallush, R. Dann, and R. Kosloff, Controlling the uncontrollable: Quantum control of open-system dynamics, Sci. Adv. 8, eadd0828 (2022).
  * Hiscock  _et al._ [2019] H. G. Hiscock, T. W. Hiscock, D. R. Kattnig, T. Scrivener, A. M. Lewis, D. E. Manolopoulos, and P. J. Hore, Navigating at night: fundamental limits on the sensitivity of radical pair magnetoreception under dim light, Q. Rev. Biophys. 52, e9 (2019).
  * Xu et al. [2021] Xu et al., Magnetic Sensitivity of Cryptochrome 4 From a Migratory Songbird, Nature 594, 535 (2021).
  * Al-Mohy and Higham [2009] A. H. Al-Mohy and N. J. Higham, Computing the fréchet derivative of the matrix exponential, with an application to condition number estimation, SIAM J. Matrix Anal. Appl. 30, 1639 (2009).
  * Virtanen et al. [2020] Virtanen et al., SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python, Nat. Methods 17, 261 (2020).

Supplementary Material

Quantum Control of Radical Pair Dynamics beyond Time-local Optimization

Farhan T. Chowdhury, Matt C. J. Denton, Daniel C. Bonser, and Daniel R.
Kattnig∗

Department of Physics and Living Systems Institute

University of Exeter, Stocker Road, Exeter, Devon, EX4 4QD, United Kingdom

E-mail<EMAIL_ADDRESS>
## Additional Results

We provide further results for the minimization of the recombination yield of
the 5-spin system in the low-magnetic field regime, including the geomagnetic
field. This is followed by additional results for the 5-spin system controlled
under $100\mskip 3.0mu\mu$T, with and without a speed limit imposed. An
additional figure shows the frequency profile and recombination yields
realized using L-BFGS-B complemented GRAPE for increasing sample intervals for
the 18 spin PY/DCB system driven by an RF control field of $0.1\mskip
3.0mu$mT.

Figure S1: Using GRAPE to control the 5-proton spin system in the geomagnetic field ($50\mskip 3.0mu\mu$T), with the singlet probability optimized over 5000 time-steps of $1\mskip 3.0mu$ns using 2 cycles of 2500 time-steps and 250 sampling points. The control amplitude amounted to $0.25\mskip 3.0mu$mT. Figure S2: Using GRAPE to control the 5-proton spin system under 100$\mu$T magnetic field, subject to a speed limit of $\dot{u}_{max}=0.2\mskip 3.0mu$ns-1 in (a), and without in (b), with the singlet probability optimized using 2 cycles of 2500 time-steps and 250 sampling points. Figure S3: We show the frequency profile and recombination yield minimization using L-BFGS-B complemented GRAPE for increasing sample intervals for the 18 spin PY/DCB system, driven by an RF control field of 0.1 mT in a static biasing field of $1\mskip 3.0mu$mT using 1024 control pulses of $1\mskip 3.0mu$ns width. Table 1: Dependence of the optimized recombination yield and runtime of the L-BFGS-B algorithm, as implemented in scipy.optimize.minimize, on the tolerancy ftol. 15 attempts to minimize the recombination yield of the 7-spin model system starting from random initial states uniformly distributed in $[-1,1]$ were recorded. The biasing field amounted to $1\mskip 3.0mu$mT and the maximal control field amplitude to $0.25\mskip 3.0mu$mT. The minimal yield realized was $0.1336$. The iterative algorithm stops when $(f^{k}-f^{k+1})/max(|f^{k}|,|f^{k+1}|,1)\leq ftol$, where $f^{k}$ is the sparse approximation of the yield at iteration $k$. For the application considered here, the denominator equals 1. $N=5000$, $B=5$, $S=200$. ftol |  Recombination Yield (minimal) |  Recombination Yield (average) |  runtime (s)
---|---|---|---
4 |  0.2665 |  0.2738 ± 0.0039 |  75.6 ± 0.8
4.5 |  0.1354 |  0.1496 ± 0.0367 |  1522 ± 480
5 |  0.1345 |  0.1376 ± 0.0019 |  2132 ± 298
6 |  0.1343 |  0.1374 ± 0.0022 |  6483 ± 543
7 |  0.1338 |  0.1375 ± 0.0022 |  8886 ± 447
8 |  0.1336 |  0.1372 ± 0.0020 |  9959 ± 113

We explore the question of robustness for our derived controls for two
pertinent scenarios, firstly in the case of “inhomogeneous broadening” due to
many (weaker) additional hyperfine interactions, and secondly in the case of
spin relaxation as induced by random field fluctuations. As demonstrated in
the next section, controls derived for the idealized systems are not
invalidated, obviously provided that their effect does not eradicate the
magnetic field effects in the absence of control. To further support the case
for global optimization over time-local methods, which roughly correspond to
what is achievable with GRAPE for insufficient block sizes, we further present
plots for 5-spin systems showing singlet recombination yields as a function of
$S$ for yield minimizations using $B$ blocks of $N/B$ steps each. The
superiority of our approach is contingent upon its ability to capture the
time-integrated nature of the reaction yield, the observable quantity, rather
than merely the singlet probability at one moment in time. We further provide
a quantitative comparison of a generalized version of the time-local
algorithm, with required adaptations of the approach described in an
additional section, which permits the optimization of the singlet probability,
for one instance in time, in the presence of spin selective recombination.

Figure S4: For 5-proton and 3-proton model systems in (a)-(b) and (c),
respectively, with a a maximal control field amplitude of $0.25\mskip
3.0mu$mT, with biasing fields of 50$\mskip 3.0mu\mu$T, $10\mskip 3.0mu$mT and
$1\mskip 3.0mu$mT for (a), (b) and (c), respectively, we show singlet
recombination yield as a function of $S$ for yield minimizations using $B$
blocks of $N/B$ steps each.

## Robustness to Noise and Unresolved Hyperfine Interactions

We consider the effects of additional hyperfine interactions, specifically
hyperfine interactions that are small compared to the principal hyperfine
interactions (which are used when deriving the controls) but numerous, which
is typical for the radicals of interest in organic radical recombination
reactions. We follow the approach of Schulten and Wolynes [2], applicable for
the assumed large number of nuclear spins, according to which the effect of
the additional hyperfine interactions can be approximated in terms of
additional static random magnetic fields, applied to the two electron spins.

Figure S5: Effect of unresolved hyperfine interactions on the singlet
recombination yield in the absence (blue) and presence (green) of the control
fields for the 7-spin model system studied in the main text. The controls were
derived for the idealized systems devoid of inhomogeneous hyperfine
interactions, which, if present, were accounted for in the Schulten-Wolynes
framework [2]. The biasing field amounted to $0\mskip 3.0mu$mT (top), $2\mskip
3.0mu$mT (middle row) and $10\mskip 3.0mu$mT (bottom row); the additional
hyperfine fields were added for both radicals (left), or radical A (with 3
dominant hyperfine interactions; central column) or radical B alone (2
dominant hyperfine interactions; right column). The effect of the additional
hyperfine interactions was evaluated through Monte Carlo sampling of 100
hyperfine field realizations; the shaded regions correspond to 95$\%$
confidence intervals of the mean recombination yield. The largest added
hyperfine fields exceed the intrinsic hyperfine fields of the radicals,
thereby demonstrating remarkable robustness of the approach to unresolved
hyperfine interactions. Figure S6: Dependence of the singlet recombination
yield of the 7-spin model spin system studied in Fig. 4 of the main text on
the rate of random-field relaxation, $\gamma$, in the absence (green) and
presence (blue) of controls. Controls were derived for the idealized system
without spin relaxation. $\gamma$ increases with increasing lightness of the
curves, comprising the values $0$, $0.25$, $0.5$, $1$, $2\mskip 3.0mu\mu$s-1.
The latter exceeds the radical pair lifetime ($k_{S}=k_{f}=1\mskip
3.0mu\mu$s-1), thus strongly attenuating the magnetic field sensitivity. Yet,
the controls remain functional and lead to reaction yields that cannot be
realized using static magnetic field of any flux density.

These additional fields comprise identically and independently normally
distributed field components of mean 0 and variance
$\sigma_{i}^{2}=\frac{1}{3}\sum_{j}a_{i,j}^{2}I_{i,j}(I_{i,j}+1)$, where the
sum runs over all additional hyperfine interactions of radical $i$. Fig. S5
illustrate the effect of such unresolved hyperfine on the singlet
recombination yield for the 7-spin model system, based on Monte Carlo sampling
hyperfine fields according to the above paradigm. The controls have been
derived based on the idealized system devoid of the additional hyperfine
fields and their resilience tested in their presence. We investigated three
different biasing field strengths, $0$, $2$ and $10\mskip 3.0mu$mT, for which
we increased $\sigma_{i}$ for one or both radicals up to $150\mskip
3.0mu$rad/$\mu$s. The maximal hyperfine field exceeds the hyperfine field of
the dominant interactions, which amounts to $100$ and $32\mskip
3.0mu$rad/$\mu$s (but realized in terms of 3 or 2 hyperfine interactions
only). As is apparent from Fig. S5, the controls remain robust even in the
presence of strong inhomogeneous broadening, realizing recombination yields
that cannot be attained through application of static fields alone. To assess
the effects of the presence of spin relaxation, we study how noise, modeled
through the Lindblad form

$\hat{\hat{R}}\hat{\rho}(t)=\sum_{i}\gamma_{i}\left[\hat{L}_{i}\hat{\rho}(t)\hat{L}_{i}^{\dagger}-\frac{1}{2}\left\\{\hat{L}_{i}^{\dagger}\hat{L}_{i},\hat{\rho}(t)\right\\}\right]$
(S1)

with
$\hat{L}_{i}\in\left\\{\hat{S}_{1x},\hat{S}_{1y},\hat{S}_{1z},\hat{S}_{2x},\hat{S}_{2y},\hat{S}_{2z}\right\\}$,
impact the control success when added to the equation of motion in the main
text. Assuming a common relaxation rate for all collapse operators
$\hat{L}_{i}$ (isotropic random field noise), Fig. S6 shows the effect of
increasing the relaxation rate $\gamma$ on the singlet recombination yield for
the controlled and unperturbed 7-spin system (cf. Fig. 4 in the main text). We
chose relaxation rates ranging up to $2\mskip 3.0mu\mu$s-1, which exceeds the
radical pair lifetime, and consequently strongly reduces the intrinsic
magnetic

field sensitivity, i.e. the curves describing the recombination yield as a
function of the biasing fields are flattened and the low-field effect
attenuated. Nonetheless, controls derived for the idealized non-relaxing
system remain effective in reducing the recombination yield beyond the yields
realizable through applying static magnetic fields, i.e. the controls are
robust.

##  VI Comparisons to Time-Local Optimization

We generalize the approach of [3] to account for spin-selective recombination
and observables that do not commute with the effective Hamiltonian. We limit
the exposition to the essential facts; the reader is referred to the original
publication for details. The spin density operator of a radical pair subject
to coherent evolution under the drift Hamiltonian $H_{0}$ and a single control
field with Hamiltonian $H_{1}$ and amplitude $u(t)$ obeys the following
equation of motion in the Liouville representation

$\frac{d}{dt}|\rho(t)\rangle=-i\left(\hat{L}_{0}+u(t)\hat{L}_{1}-i\hat{K}\right)|\rho(t)\rangle=-i\left(\hat{L}+u(t)\hat{L}_{1}\right)|\rho(t)\rangle,$
(S2)

where $L_{0}$ and $L_{1}$ are commutation superoperators, i.e.
$\hat{L}_{i}|\rho\rangle=|[H_{i},\rho]\rangle$, $\hat{K}$ accounts for
recombination, i.e. $\hat{K}|\rho\rangle=|\\{K,\rho\\}\rangle$, and
$\hat{L}=\hat{L}_{0}-i\hat{K}$ was introduced in the final equality. We use
single hats for super-operators and omit them for usual Hilbert space
operators, which induce the ket-states of the Lioville space, denoted by the
usual Dirac kets. To optimize the observable $\Omega$ at time $t_{f}$, we
introduce a moving target

$|\Omega(t)\rangle=\exp(i\hat{L}^{\dagger}(t_{f}-t))|\Omega\rangle,$ (S3)

where the adjoint operation must be applied as $L$ is not anti-Hermitian as a
consequence of the presence of $\hat{K}$ accounting for the radical
recombination. With the target state redefined as above, we follow [3] and
introduce the optimization quantity

$y=\langle\Omega(t)|\rho(t)\rangle,$ (S4)

which for $t=t_{f}$ equals the expectation value of observable $\Omega$. The
time evolution of $y$ is found as

$\frac{d}{dt}y(t)=-iu(t)\langle\Omega(t)|\hat{L}_{1}|\rho(t)\rangle.$ (S5)

Hence, choosing $u(t)$ to obey

$u(t)=iA\langle\Omega(t)|\hat{L}_{1}|\rho(t)\rangle^{*}\in\mathbb{R}$ (S6)

for any real constant $A>0$, guarantees that

$\frac{d}{dt}y(t)=A|\langle\Omega(t)|\hat{L}_{1}|\rho(t)\rangle|^{2}\geq 0.$
(S7)

Consequently, $y(t)$ will monotonously increase over time, thereby enabling a
(constraint) maximization of $y(t_{f})$. The optimization is constrained in
the sense that solutions for which $y(t)$ is non-monotonic are excluded. Note
furthermore that for $\rho(t=0)\propto P_{S}$ and $H_{1}\propto S_{x}$,
$u(t=0)=0$. For $A<0$, a minimization will be realized. Apart from the
redefinition of the target state, the approach is identical to [3], to which
we refer for additional details.

Table 2: Recombination yields of the 7-spin radical pair for $B_{0}=0\mskip 3.0mu$mT subject to the optimal controls as derived using the generalized time-local approach targeting the singlet probability at $t_{f}$. $A$ is the amplitude parameter (cf. eq. (S6) and $S$ is the relative success expressed as the change of yield induced by the controls relative to the change of yields realized via the GRAPE-derived controls from the main text. $t_{f}$ ($\mu$s) |  A |  $Y_{S}$ |  S(%) |  $\max|u(t)|$
---|---|---|---|---
1 |  -1.48 |  0.283 |  3.3 |  1.0
1 |  -1.5 |  0.279 |  6.7 |  1.5
1 |  -2 |  0.258 |  25 |  2.2
1 |  -2.5 |  0.247 |  35 |  3.0
1 |  -3 |  0.249 |  33 |  5.6
2 |  -1.5 |  0.265 |  19 |  1.2
2 |  -2.5 |  0.242 |  39 |  2.3
3 |  -2.5 |  0.279 |  6.5 |  0.9
3 |  -2.7 |  0.279 |  6.8 |  1.0
4 |  -2 |  0.282 |  4.5 |  1.1
4 |  -2.5 |  0.277 |  8.4 |  1.4
5 |  -2.5 |  0.287 |  0.0 |  0.0
5 |  -100 |  0.286 |  0.4 |  2.1

We implement the time-local optimization approach using the zvode ODE-solver
in SciPy [1], which employs the Adams/BDF-method of order 12. We consider the
problem of minimizing the singlet recombination yield of the 7-spin model
system, proxied through the singlet probability at time $t_{f}$. For no
biasing field ($B_{0}=0$), our GRAPE-inspired approach yielded a reduction of
the recombination yield from $0.287$ to $0.172$ using controls of maximal
amplitude $0.25\mskip 3.0mu$mT (corresponding to $u(t)=\pm 1$) when
controlling the first $5\mskip 3.0mu\mu$s, with controls established through
optimization of 5 blocks of 1000 control amplitudes each, and sparse sampling
of the yield.

Figure S7: Time-local optimization of the 7-spin model system studied in the
main text. The optimization was realized for $B_{0}=0\mskip 3.0mu$mT through a
generalization of the approach from [3], detailed above, with $A=-2.5$. The
moving target (cf. eq. (S3) optimized the singlet probability at time
$t_{f}=2\mskip 3.0mu\mu$s. The singlet recombination yield was reduced from
$0.287$ to $0.242$ using a maximal control amplitude $\max|u(t)|\approx 2.3$
(the GRAPE-inspired approach yielded superior minimization with
$Y_{S}=0.172$). The panels give the singlet probability, rescaled to remove
the effect of the spin-independent reaction, the relative control amplitude
(multiplying the nominal amplitude of $0.25\mskip 3.0mu$mT), and the
optimization criterion $y(t)$, which evolves monotonously in time.

First, we applied the local optimization approach for $t_{f}$ in the range
from $1$ to $5\mskip 3.0mu\mu$s, choosing the amplitude parameter $A$ (cf. eq.
(S6) such that a comparable maximal control amplitude of $\max|u(t)|\approx 1$
resulted, and as $A=-2.5$, which yielded larger control amplitudes for most
target times $t_{f}$. Furthermore, for $t_{f}=1\mskip 3.0mu\mu$s, additional
$A$-parameters up to $|A|=3$ were tried. Inspection of the control results,
summarized in Table. (II), leads to the quantitative insight that the time-
local control approach is inferior for all $t_{f}$. When limited to the
approximately same maximal control amplitude, the “best” reduction of the
recombination yields realized did not permit $Y_{S}$ below 0.265. Better
outcomes could only be achieved using (markedly) larger control amplitudes.
For example, for $t_{f}=1\mskip 3.0mu\mu$s, for $A=-2.5$ we obtained
$Y_{S}=0.247$, realized using a control amplitude that exceeded those used in
our approach by nearly a factor of 3. For $t_{f}=2\mskip 3.0mu\mu$s, $A=-2.5$
yielded $Y_{S}=0.242$ at more than twice the intended control amplitude.

We note that further increasing $|A|$ did not allow approaching the GRAPE-
based results, as is apparent from the data for $t_{f}=1\mskip 3.0mu\mu$s in
Table. (II). Even for the 5.6-fold increased maximal amplitude realized for
$A=-3$, the yield is limited to $0.249$. Fig. (S7) illustrates the
optimization approach for $t_{f}=2\mskip 3.0mu\mu$s and $A=-2.5$. Note that
the optimized $u(t)$ is approximately zero for the initial time period, i.e. a
phase of unperturbed evolution precedes the control part, which is a
characteristic feature observed for all time-local optimizations realized here
and also seen in the test cases from [3]. Given that controls start out at
zero and observing that overall better results are obtained for the shorter
$t_{f}$ tried, we wondered if our blocking approach successfully applied for
the new algorithm, as laid out in the main text, could likewise enhance the
optimization of the time-local approach, as discussed here. To this end, we
focused on the first $5\mskip 3.0mu\mu$s of the evolution, as for the GRAPE-
inspired approach, and subdivided this interval in non-overlapping blocks. The
optimization was carried out by re-defining the moving target for every block
such that it yielded the singlet probability at the respective block end time.
Indeed, in this way, we were able to realize an improvement over the results
for a single $t_{f}$. However, again, larger amplitudes were required, and the
results still fell short of the GRAPE-based approach. Table. (III) summarizes
the optimization results for $A=-2.5$, which gives rise to amplitudes markedly
exceeding $|u|=1$, for the smallest block size by a factor of 23. The best
result was obtained for 80 blocks, giving $Y_{S}=0.208$ with
$\max|u(t)|\approx 8.9$. While this yield still significantly overshoots the
GRAPE-based $Y_{S}=0.172$, we think that time-blocking might be an interesting
option for the time-local optimization approach.

Table 3: Recombination yields of the 7-spin radical pair for $B_{0}=0\mskip 3.0mu$mT subject to the optimal controls as derived using the time-local approach targeting the singlet probability for consecutive time blocks. The first $5\mskip 3.0mu\mu$s of the time evolution where controlled with controls optimized for blocks of length $t_{b}$. $A$ is the amplitude parameter (cf. eq. (S6) and $S$ is the relative success expressed as the change of yield induced by the controls relative to the change of yields realized via the GRAPE-derived controls from the main text. No. of blocks |  $t_{b}$ ($\mu$s) |  A |  $Y_{S}$ |  S (%) |  $\max|u(t)|$
---|---|---|---|---|---
5 |  1 |  -1.48 |  0.272 |  13 |  1.9
5 |  1 |  -2 |  0.258 |  25 |  2.2
5 |  1 |  -2.5 |  0.248 |  33 |  3.0
5 |  1 |  -3 |  0.248 |  34 |  5.6
10 |  0.5 |  -2.5 |  0.239 |  42 |  3.3
20 |  0.25 |  -2.5 |  0.214 |  64 |  8.4
40 |  0.125 |  -2.5 |  0.228 |  51 |  8.2
80 |  0.0625 |  -2.5 |  0.208 |  69 |  8.9
160 |  0.03125 |  -2.5 |  0.208 |  69 |  14.7
320 |  0.015625 |  -2.5 |  0.215 |  62 |  21.1
Figure S8: Time-blocked time-local optimization of the 7-spin model system
studied in the main text. The optimization was realized for $B_{0}=0\mskip
3.0mu$mT through a generalization of the approach from Ref. [3], detailed in
the appendix, with $A=-2.5$, and application of a time-blocking approach like
that used here for the GRAPE-inspired optimization. The first $5\mskip
3.0mu\mu$s of evolution were subdivided in 10 blocks, to which the time-local
optimization with moving target optimizing the singlet probability at the
blocks’ ends was applied. The redefinitions of the targets give rise to jumps
in $y(t)$, which otherwise decreases monotonously, as enforced by the time-
local approach. The singlet recombination yield was reduced from $0.287$ to
$0.239$ using a maximal control amplitude (the GRAPE-inspired approach yielded
superior minimization with $Y_{S}=0.172$). Using 20 to 160 blocks gave rise to
better optimization results (cf. Table. III), but the characteristic features
of $y(t)$ are less apparent.

## References

  * [1] Virtanen et al.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  Nat. Methods, 17:261–272, 2020.  DOI: https://doi.org/10.1038/s41592-019-0686-2.  URL https://rdcu.be/b08Wh.
  * [2] Klaus Schulten and Peter G. Wolynes.  Semiclassical description of electron spin motion in radicals including the effect of electron hopping.  The Journal of Chemical Physics, 68(7):3292–3297, 08 1978.  DOI: 10.1063/1.436135.  URL https://doi.org/10.1063/1.436135.
  * [3] Kenta Masuzawa, Masaya Sato, Michihiko Sugawara, and Kiminori Maeda.  Quantum Control of Radical Pair Reactions by Local Optimization Theory.  J. Chem. Phys., 152(1):014301, Jan 2020.  DOI: 10.1063/1.5131557.  URL https://doi.org/10.1063/1.5131557.