Quantum Control of Radical Pair Dynamics beyond Time-Local Optimization

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\'echet 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 ultrafast 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 radiofrequency 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 symmetrybreaking interactions-in particular, the hyperfine interactions with surrounding magnetic nuclei.Quite recently, quantum beats reflecting the coherent singlet-triplet in-terconversion 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 radiofrequency 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 exciplexforming 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 lowfield 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 singlettriplet 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 timeindependent 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 fluorescencebased organic light emitting diodes with large magnetoelectroluminescence [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.

A. 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.

FIG. 1. Reaction Scheme.
The corresponding spin dynamics is described in terms of the time-dependent spin density matrix ρ(t) by with Ĥ(t) denoting the spin Hamiltonian and accounting for radical pair recombination.Here, k S and k T denote the reaction rate constants in the singlet and triplet state, respectively, and PS , PT are the projection operators onto the singlet and triplet states.We introduce an effective Hamiltonian given by Â allowing us to reformulate eq. ( 1) as where The formal solution to eq. ( 4) is given by where the time evolution operator in terms of Â(t) is Here, 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 where p S (t) = T r[ PS ρ(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 Ĥ0 = ĤA + ĤB + ĤAB (8) where Ĥi , i ∈ {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.Ĥ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 Îik and Ŝi by and where ω i = −γ i B, with γ 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 where ) is the total number of nuclear spin states in radical i and is the singlet projection operator.In the context of spinbiological systems, we also consider directional magnetic field effects of systems with anisotropic hyperfine interactions.In this case, Ĥi is of the form where n(ϑ, φ) is an unit vector with polar angle ϑ and azimuth φ defining the magnetic field direction and the A ik s are the hyperfine tensors.

B. 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 where Ĥ0 is the drift Hamiltonian describing the timeindependent 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 Vl 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 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, with constants u l,n and the indicator function χ n (t) defined such that χ n (t) = 1 for t n−1 ≤ t < t n 0 otherwise. .
The corresponding propagators Û (t n , 0) are discretized on the temporal grid such that where with ∆t i = t i − t i−1 and Âi = Â(t i−1 ).The gradients of G depend on the propagator derivatives, evaluated as where L(X, E) is the Fréchet derivative of the matrix exponent defined, for every E, by More explicitly, we have [35] L(X, E) 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.
which also allows us to demonstrate that Combining the above equations, the gradient of the survived singlet probability, for k ≤ n, is thus where Ŵi,l = i Û −1 i L i,l and Ûm .
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 Ŵi,l = Vl .
For our approach, we use the full Fréchet derivative.Finally, this allows us to approximate the singlet yield by with w n denoting suitably chose quadrature weights.The gradient of the singlet yield is then obtained as

III. RESULTS
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 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.

A. 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 singletprobabilities 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 FIG. 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 l,n are simultaneously updated, whereby B disjoint blocks of 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 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 mT.The RF control field of 0.1 mT amplitude comprised 1024 control pulses of 1ns 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 µs − 1.The Py/DCB system is discussed in more detail in the Results section.
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 gradientbased 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 timesteps 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].

B. Protocol applied to a prototypical spin system
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 highfield scenario using time-local optimization.Specifically, we assume a radical pair comprising 7 spin-1 2 particles, coupled through isotropic hyperfine (0.2, 0.5, 1 mT and 0.2 & 0.3 mT for the two radicals, respectively) and exchange interactions (j ex /(2π) = 1 MHz), undergoing spin-selective recombination in the singlet state with rate constant k b = 1 µs −1 or spin independent escape with rate constant k f = 1 µs −1 , such that k S = k b + k f and k 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 through a piece-wise constant control field applied perpendicular to a static biasing field, i.e.L = 1 and , where ω 1 is the maximal Rabi frequency of the control field; control amplitudes are sub- We chose a discretization time step of 1 ns and controlled the first 5µ 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 → ∞ ∼ t max limit with t max = 14µ 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 mT and 10 mT.
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 supe-FIG.4. Controlled recombination dynamics of a 7-spin radical pair system under 0.5 mT and 10 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. pS(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 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.FIG. 5.For the 7 and 5-spin model systems, a maximal control field amplitude of 0.25 mT, with biasing fields of 0.5 and 10 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 ±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.
rior.This is, for example, seen in Fig. S4 for a biasing field of 1 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 mT, is more sensitive to the sampling interval, S, than for 10 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 Ma-terial 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).
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 ns and sampling the recombination yield from every 25th time point.The maximal control amplitude was 0.25 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 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 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 mT and FIG. 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 mT control field, and show yield minimization against increasing B0 (blue) in comparison to the un-controlled reference (orange).In (b) we show yield minimization against increasing control field amplitude, B1, for a biasing field of B0 = 10 mT and 50 µT.In (c) we show the dependence of the singlet recombination yield at B0 = 1 mT on the speed limit u−1 max .For all optimizations, N = 5000, B = 5 and S = 25.
50 µT, again realized by optimizing 5 blocks of 1000 time steps, sampling every 25 th .Appreciable minimization for B 0 = 10 mT requires control amplitudes of 0.1 mT; controls with field strength below ≈ 50 µT do not elicit appreciable effects.This finding is in line with the wellstudied 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, τ , if -γB 1 /(2π)τ ⪆ 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 −1 f = 1 µs and k −1 b = 1 µs, we can expect effects for B 1 exceeding 36 µT to 72 µ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 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 µT) slightly larger control fields exceeding 0.12 mT are required.
Above we have assumed that controls are applied with a time-resolution of ∆t = 1 ns.While AWGs are available with sufficient sample rates and analog bandwidths (e.g.50 GS/s at 15 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 | < umax ∆t.For B 0 = 1 mT, Fig. 6c summarizes the dependence of the singlet yield on this speed limit.We find that transition times of u−1 max ⪅ 8 ns are required to elicit significant control; fast controls only deliver marginal additional gains.The minimal speed limit corresponds to a bandwidth of approximately umax 2 = 63 MHz, which is of the order of, but smaller than, the magnitude of the spread of energy eigenvalues of the radical pair (93 MHz).

C. Protocol applied to a realistic exciplex-forming radical pair
We now present results for coherently controlling the symmetric recombination of a large spin system consisting of 18 spins.We consider the pyrene/paradicyanobenzene (PY/DCB) exciplex-forming donoracceptor systems in a solvent of moderate polarity [16,37,38].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 S = k 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 α,β (t)R (2) where 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, Û accounts for the coherent evolution only, i.e.Â(T ) = Ĥ(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].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 [9,16,18,37,38] by building on the time-resolved MFE measurements and studies of radio-frequency effects that have already been realized.

D. 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.
The flavin semiquinone (FADH • )/superoxide (O •− 2 ) 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 , 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, FIG. 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 µs −1 and spin-selective recombination in the singlet state with rate constant k b = 1 µs −1 .The first 5 µs after radical pair generation were controlled using a step-wise constant control field with 1 ns time resolution and a maximal amplitude 0.25 mT, applied perpendicular to the static magnetic field.The control was optimized using the L-BFGS-B algorithm applied to 5 blocks of 1 µs, sampling the yield every 25 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 mT biasing field in absence and presence of the optimized control field, scaled to remove the spin-independent decay, i.e. pS(t) exp(k f t).The associated singlet recombination yields amount to 0.218 and 0.115 with and without controls, respectively.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 [41,48] 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 2 ) 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.
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 n max and n min , i.e. ∆ S = Y S (n max )−Y S (n 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, ∆ 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 ∆ 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 pa-FIG.9. Change in singlet yield contrast ∆S for FAD •− /W •+ and FADH • /Z • radical pairs subject to controls optimized to maximise or minimise the singlet recombination yield contrast ∆S in the geomagnetic field for one of the pairs (as indicated in the axis labels).The geomagnetic field was assumed as 50 µT and the radical pair lifetime 1 µs, i.e. kS = kT = 1 µs −1 .The control was applied parallel to the static field with an amplitude of 12.5 µT.The first 5 µs after radical pair generation in the singlet state were controlled with a time-step of 1 ns.Yields were calculated by integrating to tmax = 14 µs.The optimization was realized through 40 steps of steepest descent/ascent.The unperturbed singlet yield contrasts are ∆S = 0.017 and 0.13 for FAD •− /W •+ and FADH • /Z • , respectively.In b), we show the accumulated reaction yield differences for each, kS t 0 pS(nmax, t) − pS(nmin, 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 ∆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 MHz, demonstrating that designed perturbations are significantly more disruptive than unspecific perturbations.
rameters [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 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 radiofrequency 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 [4,57,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 [17,60] qubits highly sensitive to applied magnetic fields.This could translate to further applications [55,[61][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 [7,61,[63][64][65][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 ∈ [0, ∞), 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 timelocal approach from [8,24,25]), 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 % 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 timelocal 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 7spin 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 [65,70,71].In this scenario, quantum control could allow us to derive con-trol 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 GRAPEinspired 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 %, 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.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 matrix multiplications, i.e. for every tn, (n − 1) matrix multiplications are required to evaluate Û , (n − 1) to evaluate Û+ Ŵn,l Û− (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 Û and Û+ Ŵn,l Û− with a singular matrix multiplication (instead of n − 1), thus requiring 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 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 ≪ N , to realize a substantial speedup, whilst delivering adequate practical optimization of the reaction yield.Furthermore, for S ≫ 3, the complexity of the algorithm approaches independence of S.

B. 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 e X = e 2 −s X 2 s (A4) by s ∈ N repeated squaring starting from e 2 −s X calculated as being a convenient splitting of X = −idt Âi (cf.eq. ( 19).Practically, Â′ and Â′′ represent the drift and control Hamiltonians, respectively.If the latter is based on coupling applied magnetic fields via the Zeeman interaction, the propagator under Â′′ can be expressed as a direct product of evolution operators for the individual electron spins based on expression where ω ∈ 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], 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 ∆tA ′′ .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 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.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.
FIG. 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 mT using 1024 control pulses of 1 ns width.

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.
These additional fields comprise identically and independently normally distributed field components of mean 0 and variance σ 2 i = 1 3 j a 2 i,j 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 mT, for which we increased σ i for one or both radicals up to 150 rad/µs.The maximal hyperfine field exceeds the hyperfine field of the dominant interactions, which amounts to 100 and 32 rad/µ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 with Li ∈ Ŝ1x , Ŝ1y , Ŝ1z , Ŝ2x , Ŝ2y , Ŝ2z , impact the control success when added to the equation of motion in the main text.Assuming a common relaxation rate for all collapse operators Li (isotropic random field noise), Fig. S6 shows the effect of increasing the relaxation rate γ 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 µs −1 , which exceeds the radical pair lifetime, and consequently strongly reduces the intrinsic magnetic ), 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.
TABLE II.Recombination yields of the 7-spin radical pair for B0 = 0 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.First, we applied the local optimization approach for t f in the range from 1 to 5 µs, choosing the amplitude parameter A (cf. eq.(S6) such that a comparable maximal control amplitude of max |u(t)| ≈ 1 resulted, and as A = −2.5, which yielded larger control amplitudes for most target times t f .Furthermore, for t f = 1 µ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 µ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 µ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 µ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 µ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 µ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)| ≈ 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 III.Recombination yields of the 7-spin radical pair for B0 = 0 mT subject to the optimal controls as derived using the time-local approach targeting the singlet probability for consecutive time blocks.The first 5 µ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
FIG.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.

FIG. 7 .
FIG.7.Recombination yield minimization of the 18 spin PY/DCB system in a bias field of 1 mT.(a) shows the singlet probability, scaled to remove the spin-independent decay due to radical recombination, i.e. pS(t) exp(kt), with and without control.Control amplitudes were optimized using 32 samples and 1024 time-steps of 1 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.

Fig. 7
shows the minimization of the PY/DCB singlet recombination yield for k −1 = 200 ns, a control amplitude of 0.2 mT, and a static basing field of 10 mT.The control was extended to the first 1024 time steps of 1 ns each.

ACKNOWLEDGMENTS
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.
FIG. 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 mT, with biasing fields of 50 µT, 10 mT and 1 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.

FIG. S5 .
FIG.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 mT (top), 2 mT (middle row) and 10 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.
FIG.S7.Time-local optimization of the 7-spin model system studied in the main text.The optimization was realized for B0 = 0 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 µs.The singlet recombination yield was reduced from 0.287 to 0.242 using a maximal control amplitude max |u(t)| ≈ 2.3 (the GRAPE-inspired approach yielded superior minimization with YS = 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 mT), and the optimization criterion y(t), which evolves monotonously in time.

TABLE I .
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 mT and the maximal control field amplitude to 0.25 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) ≤ f tol, 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.