Anisotropic and Coherent Control of Radical Pairs by Optimized RF Fields

Abstract

Radical pair kinetics is determined by the coherent and incoherent spin dynamics of spin pair and spin-selective chemical reactions. In a previous paper, reaction control and nuclear spin state selection by designed radiofrequency (RF) magnetic resonance was proposed. Here, we present two novel types of reaction control calculated by the local optimization method. One is anisotropic reaction control and the other is coherent path control. In both cases, the weighting parameters for the target states play an important role in the optimizing of the RF field. In the anisotropic control of radical pairs, the weighting parameters play an important role in the selection of the sub-ensemble. In coherent control, one can set the parameters for the intermediate states, and it is possible to specify the path to reach a final state by adjusting the weighting parameters. The global optimization of the weighting parameters for coherent control has been studied. These manifest calculations show the possibility of controlling the chemical reactions of radical pair intermediates in different ways.

1. Introduction

A radical pair (RP) is composed of two radicals, each of which has an unpaired electron. When RPs are produced photochemically, the initial electron spin state of the RPs preserves the spin manifold of the precursor excited state. Although this fact has been known for a long time, it has been re-evaluated and has attracted much attention because this state is the so-called “entanglement state” [1,2,3,4,5]. On the other hand, the formation of RPs is also known in photoreceptor protein systems, such as photosynthetic reaction centers [6,7,8,9] and cryptochromes (CRYs) [10,11,12,13,14]. In particular, CRYs are candidate molecules for magnetoreception in animals, and the relationship between magnetism and biological systems has been discussed [15,16,17,18,19]. The influence of not only magnetic fields but also that of electromagnetic waves has been discussed [20]. In the geomagnetic perception of migratory birds, the influence of electromagnetic noise is still under discussion [21].

The most probable effect of electromagnetic waves on RPs is magnetic resonance [22]. It is called reaction yield detected magnetic resonance (RYDMR) [23,24,25,26,27,28,29] and is measured at different external magnetic fields and resonance frequencies. Figure 1 shows a schematic diagram of the spin state and chemical reactivity of RPs. Here, we consider the typical case where a singlet-born radical pair reacts spin-selectively from the singlet state under the magnetic field. At the high magnetic field, the microwave frequency matching with the Zeeman splitting energy induces the resonance phenomena, producing ${\mathrm{T}}_{+}$ and ${\mathrm{T}}_{-}$ states, which do not have singlet properties and are retained for a long time or become free radicals by radical diffusion. Therefore, recombination reactions from the singlet state are significantly suppressed by the magnetic resonance phenomena. On the other hand, in a small magnetic field at the geomagnetic level, Zeeman splitting is smaller than hyperfine coupling, so the resonance of the radiofrequency to hyperfine coupling is important. Such magnetic resonance phenomena of RPs are complex, but the resulting change in the singlet character of RPs affects the reaction yield. In contrast to Figure 1, the spin states of general RPs are more complex. The resonance frequency of RPs is the same as the electron spin resonance frequency of the two constituent radicals when the interaction between the radicals is negligibly small. However, they are complicatedly split by hyperfine coupling and its anisotropy. If high-field approximation holds and mixing between the nuclear spin states is small, these resonance lines can be considered as sub-ensembles. In such a case, the selection of these sub-ensembles can be conducted by selecting the RF frequency.

The spin states of a radical pair can be regarded as a quantum mechanical system. Therefore, the control of chemical reactions associated with the spin manipulation of RPs is said to be a quantum state manipulation. In the field of magnetic resonance, an optimization technique called GRAPE (Gradient Ascent Pulse Engineering) [30] is widely used. GRAPE optimizes the waveform globally on the time axis, but it tends to be trapped in local optimal states due to the combination of many segments of the RF field. In this sense, it is not a perfect global optimization for multiple variables. Therefore, in many cases, it is necessary to set many waveforms as initial conditions to find the global optimal waveform, which inevitably increases the computational complexity. Therefore, the application of this method is limited to a small number of short pulse designs and is difficult to optimize the RF field for the chemical reaction control to occur over long periods of time.

On the other hand, a local optimization method with respect to the time axis has been used to design an infrared laser field for vibronic state control [31]. This method works very well for relatively simple quantum states and the optimal wave can be obtained in a very short computation time. In contrast to the GRAPE method, this method requires essentially no iterative computation. Thus, the optimized waveform and the resulting time evolution of the system are obtained simultaneously in a single calculation. The method is based on the concept of “living for the day” and monotonically increases fidelity, i.e., the overlap of the system with the target state. In fact, although it is a very crude method, it shows excellent results for quantum systems achieved with relatively simple paths.

We have applied this local optimization method to control the chemical reaction of RPs. This method has the following features:

(1)

Smooth waveforms can be obtained because waveforms can be computed with very high time resolution. In contrast, the waveforms computed in GRAPE are optimized for relatively coarse time steps.

(2)

Sub-ensemble systems can be easily selected by resonance frequency in inhomogeneous RP systems.

(3)

The calculation can incorporate strategies for the time evolution of the quantum system with respect to its path.

(4)

Calculations in a rotational coordinate system are easy because there is no computational burden due to the multidimensionality of radio waves (microwaves).

In particular, sub-ensemble selection (2) and path control in coherent control (3) can be achieved by using the same sets of parameters, called weighting parameters {$q_i$}, which will be discussed in the next section. In a previous paper [32], we demonstrated (2) nuclear spin state selections under (3) high-field approximation. In the present paper, we present anisotropic reaction control and coherent control in an extremely low magnetic field, focusing on the weighting parameters. Furthermore, we have improved the local-optimization-based coherent control of the spin system of RPs via the global optimization of a few weighting parameters. Of course, this method is not perfect global optimization, but it is a significant improvement. The model system used here is rather simple compared to realistic systems. However, once the methodology is developed and a strategy for reaction control is established, it should be easy to apply them to a realistic system because the computational burden is significantly lower.

2. Theory

We have previously applied the local optimization theory, which was used for the control of vibronic systems using an infrared laser [31], to the design of RF magnetic field for the control of RP chemical reactions [32]. In order to extend computational flexibility and to take into account spin-selective chemical reactions and spin relaxations, the basic theory of the local optimization theory has been reformulated in Liouville space. In this chapter, we provide an overview of the theory and describe the important weighting parameters in control. Then, anisotropic control and coherent pathway control are formulated as an example of control using weighting parameters.

2.1. Local Optimization Theory

We begin with an overview of the local optimization theory for our reaction control [31,32]. In the local optimization theory, the degree of achievement of a desired quantum state is set as a performance index, and control is realized by designing an external oscillating field such that the performance index increases monotonically. Defining ${|P\rangle }_{ }$ as the projection operator of the target quantum state (where ${|P\rangle }_{ }$ considers eigenstates of system Hamiltonian but can be extended to non-eigenstates), the performance index $y\left(t\right)$ is described as follows.

$$y\left(t\right)=\langle P\left|\rho \right(t)\rangle$$

(1)

$\left|\rho \right(t)\rangle$ is the density matrix of the system, and its time evolution follows the Liouville von Neuman equation in Liouville space.

$$\frac{d{|\rho \left(t\right)\rangle }_{ }}{dt}=-i\left[{\widehat{\widehat{L}}}_0+u\left(t\right){\widehat{\widehat{L}}}_V\right]{|\rho \left(t\right)\rangle }_{ }$$

(2)

where ${\widehat{\widehat{L}}}_0+u\left(t\right){\widehat{\widehat{L}}}_V$ is composed of the static term, ${\widehat{\widehat{L}}}_0$, which drives the time evolution by the static Hamiltonian $\widehat{H}$, and time dependent term, $u\left(t\right){\widehat{\widehat{L}}}_V$, which represents the interaction with the optimized oscillating RF fields $u\left(t\right)$:

$${\widehat{\widehat{L}}}_0^{ }=\widehat{H}-\tilde{H},\hspace{1em}{\widehat{\widehat{L}}}_V=\frac{g_A{\mu }_B}{\hslash }\left({\widehat{S}}_x^A-{\tilde{S}}_x^A\right)+\frac{g_B{\mu }_B}{\hslash }\left({\widehat{S}}_x^B-{\tilde{S}}_x^B\right)$$

(3)

where hat and tilde denote the operation to ket and bra in Hilbert spaces, respectively.

From ${\widehat{\widehat{L}}}_0^{ }{|P\rangle }_{ }=0$, the time derivative of the performance index $y\left(t\right)$ is written by:

$$\frac{dy\left(t\right)}{dt}=-iu\left(t\right)\langle P|{\widehat{\widehat{L}}}_V{|\rho \left(t\right)\rangle }_{ }$$

(4)

Now, we introduce an amplitude parameter $A$, which determines the degree of optimization as well as the weighting parameter, which will be discussed later, and if $u\left(t\right)$ is given by

$$u\left(t\right)=iA\langle P|{\widehat{\widehat{L}}}_V|\rho \left(t\right)\rangle \hspace{1em}\left(A>0\right)$$

(5)

a monotonical increase in $y\left(t\right)$ is guaranteed.

$$\frac{dy\left(t\right)}{dt}=A{\left[\langle P|{\widehat{\widehat{L}}}_V|\rho \left(t\right)\rangle \right]}^2\ge 0$$

(6)

Substituting Equation (5) into Equation (2), we obtain the nonlinear differential equation shown in Equation (7).

$$\frac{d{|\rho \left(t\right)\rangle }_{ }}{dt}=-i\left[{\widehat{\widehat{L}}}_0+iA{\widehat{\widehat{L}}}_V\langle P|{\widehat{\widehat{L}}}_V|\rho \left(t\right)\rangle \right]|\rho \left(t\right)\rangle$$

(7)

Once we solve the time evolution of $|\rho \left(t\right)\rangle$, we can obtain $u\left(t\right)$ by Equation (5).

2.2. Weighting Parameters

Equation (1) was an evaluation function of control targeting a single quantum state, but to target multiple quantum states, we extend Equation (1) as follows.

$$y\left(t\right)={\displaystyle \sum }_i{q}_i\langle P_i|\rho \left(t\right)\rangle$$

(8)

Here, we introduce the weighting parameters $q_i$, which determine the priority to the target states ${|P_i\rangle }_{ }$. They are useful for the following two types of control.

One is sub-ensemble control, in which each sub-ensemble, such as different nuclear spin configurations under a high field approximation, isotope selection, and any inhomogeneous spectral broadening of the spin systems. Here, we focus on the orientation of an anisotropic system with respect to laboratory coordinates. By setting weighting parameters for each sub-ensemble, it is possible to induce ($q_i>0$) or suppress ($q_i<0$) transitions only in certain sub-ensembles.

The other is path control. By setting the weighting parameters $q_1<q_2<\cdots <q_N$ for the quantum states $|P_1\rangle ,|P_2\rangle ,\cdots ,|P_N\rangle$, transitions can occur in the order $|P_1\rangle \to |P_2\rangle \to \cdots \to |P_N\rangle$. In other words, these parameters can be thought of as guiding the transition path to produce a particular spin state at a particular time. From the next chapter, we will formulate these controls in a radical pair system.

2.3. Anisotropic Sub-Ensemble Control

Here, we consider RPs with an anisotropic hyperfine interaction. The sub-ensemble can be determined by the orientation with respect to the laboratory coordinates. The Hamiltonian of the anisotropic hyperfine interaction is described as follows.

$${\widehat{H}}_{HFC}=\widehat{S}·{\widehat{A}}_{arb}·\widehat{I}$$

(9)

The operators $\widehat{S}$ and $\widehat{I}$ correspond to the electron and nucleus spin operators, respectively, and ${\widehat{A}}_{arb}$ is a hyperfine tensor. ${\widehat{A}}_{arb}$ is diagonalized by multiplying the appropriate rotation matrix:

$$\widehat{A}=\left(\begin{array}{ccc}{A}_{x_M{x}_M}& 0& 0\\ 0& A_{y_M{y}_M}& 0\\ 0& 0& A_{z_M{z}_M}\end{array}\right)$$

(10)

The subscripts $x_M, y_M, z_M$ correspond to $x, y, z$ in the molecular coordinate system. To describe the anisotropic response of a molecule to the direction of a static magnetic field $B_0$, the molecular coordinate system needs to be transformed into a laboratory frame. The orientation of these two coordinates systems is defined by the following Euler angles as shown in Figure 2. The conversion of the coordinates can be described by a rotational matrix $R\left(\alpha , \beta , \gamma \right)$:

$$R\left(\alpha ,\beta ,\gamma \right)=\left(\begin{array}{ccc}\cos \gamma \cos \beta \cos \alpha -\sin \gamma \sin \alpha & \cos \gamma \cos \beta \sin \alpha +\sin \gamma \cos \alpha & -\cos \gamma \sin \beta \\ -\sin \gamma \cos \beta \cos \alpha -\cos \gamma \sin \alpha & -\sin \gamma \cos \beta \sin \alpha +\cos \gamma \cos \alpha & \sin \gamma \sin \beta \\ \sin \beta \cos \alpha & \sin \beta \sin \alpha & \cos \beta \end{array}\right)$$

(11)

By using $R\left(\alpha , \beta , \gamma \right)$, the hyperfine tensor in laboratory frame ${\widehat{A}}_L$ can be calculated by:

$${\widehat{A}}_L\left(\alpha ,\beta ,\gamma \right)=R\left(\alpha ,\beta ,\gamma \right) \widehat{A} R^{-1}\left(\alpha ,\beta ,\gamma \right)$$

(12)

Here, we introduce the new index $k=\left\{\alpha , \beta , \gamma \right\}$ for the sub-ensembles of RPs. The Hamiltonian of a radical pair with one nuclear spin ($I=1/2$) for a sub-ensemble $k$ is given by:

$${\widehat{H}}_L^k=\frac{g_A{\mu }_B{B}_0}{\hslash }{\widehat{S}}_z^A+\frac{g_B{\mu }_B{B}_0}{\hslash }{\widehat{S}}_z^B+\frac{g_N{\mu }_N{B}_0}{\hslash }{\widehat{I}}_z^A+{\widehat{S}}^A·{\widehat{A}}_L\left(k\right)·{\widehat{I}}^A-J\left(\frac{1}{2}+2{\widehat{S}}^A·{\widehat{S}}^B\right)$$

(13)

The first and second terms correspond to the electron spin Zeeman interaction of radical A and radical B, the third term indicates the nuclear Zeeman interaction, the fourth term denotes the hyperfine interaction, and the fifth term represents the exchange interaction of the RPs.

Since the index $k$ represents the sub-ensemble of orientations relative to the magnetic field direction of the system, the total density vector ${|\rho \rangle }_{ }$ can be divided into the vector ${|{\rho }_k\rangle }_{ }$ for each sub-ensemble $k$. Thus, the performance index $y\left(t\right)$ can be defined by:

$$\begin{array}{ll}y\left(t\right)& ={{\displaystyle \int }}_k^{ }{q}_k\langle P{}_k|{\rho }_k{\left(t\right)\rangle }_{ }\left|\sin \beta \right|\mathrm{d}k\\ & ={\displaystyle \sum }_k{q}_k^j\langle P{}_k|{\rho }_k{\left(t\right)\rangle }_{ }\end{array}$$

(14)

$P_k$ is the projection operator for the target states, and $q_k$ is the weighting factor for $P_k$. $q_k^j$ is the product of the weighting factor $q_k$ and the Jacobian $\left|\sin \beta \right|$ ($q_k^j=q_k\times \left|\sin \beta \right|$). The $P_k$ can be constructed as follows:

• The eigenstates of ${\widehat{H}}_L^k$ were calculated.
• The two eigenstates with the largest ${\mathrm{T}}_{+\mathsf{\alpha }}$ and ${\mathrm{T}}_{-\mathsf{\beta }}$ characteristics were chosen, where the state ${\mathrm{T}}_{+\mathsf{\alpha }}$ and ${\mathrm{T}}_{-\mathsf{\beta }}$ indicate ${\mathrm{T}}_{+}$ electron spin state with α nuclear spin and ${\mathrm{T}}_{-}$ with β nuclear spin, respectively.
• The sum of the projection operators onto these two eigenstates was determined to be the target, $P_k$.

The magnetic resonance transitions from the singlet rich eigenstates to the target steady states $P_k$ increase $y\left(t\right)$ and provide long-lived RP sub-ensembles. By an appropriate setting of the weighting parameters, $q_k$ one can select the specific sub-ensembles and transfer the population to the stable eigenstates that are free from the recombination reaction.

The kinetics of the RPs reaction of each sub-ensemble $k$ can be calculated by the Haberkorn super operator, ${\widehat{\widehat{L}}}_K,$ using the projection operator ${\widehat{P}}_S=\left|S\rangle \langle S\right|$ [33,34,35].

$${\widehat{\widehat{L}}}_K=\frac{-i{k}_S}{2}\left({\widehat{P}}_S+{\tilde{P}}_S\right)$$

(15)

Since $\left({\widehat{\widehat{L}}}_0^k+{\widehat{\widehat{L}}}_K\right){|P\rangle }_L\ne 0$, occasionally $y\left(t\right)$ decreases by the recombination kinetics. Therefore, this calculation is limited to the condition where $k_S$ is much smaller than the frequency of spin dynamics. In this model calculation, we have set $k_S=1\times {10}^6 {\mathrm{s}}^{-1}$, such that the introduction of recombination reactions has a minimal effect on the optimized waveforms and does not break down the calculation.

The calculation process of anisotropic reaction control is performed in the same way as the previous calculation, i.e., solving the Liouville von Neuman equation, Equation (17) using optimized $u\left(t\right)$, which guarantees a monotonous increase in $y\left(t\right)$:

$$u\left(t\right)=iA{\displaystyle \sum }_k{q}_k^j\langle P_k|{\widehat{\widehat{L}}}_V|{\rho }_k{\left(t\right)\rangle }_{ }\hspace{1em}\left(A>0\right)$$

(16)

$$\frac{d{|{\rho }_k\left(t\right)\rangle }_{ }}{dt}=-i\left[{\widehat{\widehat{L}}}_0^k+{\widehat{\widehat{L}}}_K+iA{\widehat{\widehat{L}}}_V{\displaystyle \sum }_k{q}_k^j\langle P_k|{\widehat{\widehat{L}}}_V|{\rho }_k\left(t\right)\rangle \right]{|{\rho }_k\left(t\right)\rangle }_{ }$$

(17)

2.4. Coherent Pathway Control

In coherent control, the non-eigenstate of the Hamiltonian is set as the target state. Since the non-eigenstate evolves automatically in time, the target state ${|f\rangle }_{ }$ should be set at a specific time $t_f$, and reverse time evolution should be introduced.

$${|{\rho }_f\rangle }_{ }=exp\left\{i{\widehat{\widehat{L}}}_0\left(t_f-t\right)\right\}{|f\rangle }_{ }$$

(18)

where ${\widehat{\widehat{L}}}_0$ is a Liouvillian of field-free Hamiltonian ${\widehat{H}}_0$. We call this time-dependent target state ${|{\rho }_f\rangle }_{ }$ the moving target. If the density vector $|\rho \rangle$ of the system can be reached to the moving target by applying the external field at time $t$, the desired state ${|f\rangle }_{ }$ can be obtained automatically at time $t_f$. In a magnetic field, the singlet state becomes an unsteady state due to S-T₀ mixing. Therefore, the singlet state is a candidate for the target of coherent control. If the singlet state can be maximized at $t_f$, a spin-selective reaction, such as a recombination reaction, can be promoted. However, depending on the target state, it may be difficult to directly reach the moving target. Even in such a case, the moving target can be obtained by introducing the intermediate state $|i\rangle$ and the weighting factor $q$ [36].

$${|{\rho }_{i\to f}\rangle }_{ }=exp\left\{i{\widehat{\widehat{L}}}_0\left(t_f-t\right)\right\}\left[q_i{|i\rangle }_{ }+q_f|f\rangle \right]$$

(19)

Here, the roles of $q_i, q_f$ are not the same as $q_k$ in Section 2.3, but are rather parameters for specifying the transition path. By setting $q_f>q_i$, it becomes possible to specify a transition path that transitions to the state ${|f\rangle }_{ }$ at time $t_f$ via the state $|i\rangle$. In such a way, one can input strategies into the spin dynamics of RPs.

To evaluate the degree of achievement of $|{\rho }_{i\to f}\rangle$, the performance index $y\left(t\right)$ is defined as:

$$y\left(t\right)=\langle {\rho }_{i\to f}|\rho \left(t\right)\rangle$$

(20)

The time evolution of ${|\rho \left(t\right)\rangle }_{ }$ is described by the Liouville von Neuman equation, so the time derivative of $y\left(t\right)$ is given by:

$$\begin{array}{ll}\frac{dy\left(t\right)}{dt}& =\left[q_i\langle i|+q_f\langle f|\right]i{\widehat{\widehat{L}}}_{0}exp\left\{i{\widehat{\widehat{L}}}_0\left(t-t_f\right)\right\}{|\rho \left(t\right)\rangle }_{ }+\langle {\rho }_{i\to f}|\frac{d{|\rho \left(t\right)\rangle }_{ }}{dt}\\ & =i\langle {\rho }_{i\to f}|{\widehat{\widehat{L}}}_0{|\rho \left(t\right)\rangle }_{ }-i\langle {\rho }_{i\to f}|\left[{\widehat{\widehat{L}}}_0+{\widehat{\widehat{L}}}_{V}u\left(t\right)\right]{|\rho \left(t\right)\rangle }_{ }\\ & =-iu\left(t\right)\langle {\rho }_{i\to f}|{\widehat{\widehat{L}}}_V{|\rho \left(t\right)\rangle }_{ }\end{array}$$

(21)

As in previous discussions, by defining $u\left(t\right)$ as Equation (22), a monotonous increase in Equation (20) is guaranteed (Equation (23)).

$$u\left(t\right)=iA\langle {\rho }_{i\to f}|{\widehat{\widehat{L}}}_V{|\rho \left(t\right)\rangle }_{ }\hspace{1em}\left(A>0\right)$$

(22)

$$\frac{dy\left(t\right)}{dt}=A{\left[\langle {\rho }_{i\to f}|{\widehat{\widehat{L}}}_V|\rho \left(t\right)\rangle \right]}^2\ge 0$$

(23)

As a result, the equation that the coherent controlled system follows is following a nonlinear differential equation.

$$\frac{d{|\rho \left(t\right)\rangle }_{ }}{dt}=-i\left[{\widehat{\widehat{L}}}_0+iA{\widehat{\widehat{L}}}_V\langle {\rho }_{i\to f}|{\widehat{\widehat{L}}}_V|\rho \left(t\right)\rangle \right]|\rho \left(t\right)\rangle$$

(24)

In the coherent control calculations, we assume a single nuclear spin as in the anisotropy calculations. However, we consider only isotropic hyperfine coupling. The spin Hamiltonian is given by:

$${\widehat{H}}_0=\frac{g_A{\mu }_B{B}_0}{\hslash }{\widehat{S}}_z^A+\frac{g_B{\mu }_B{B}_0}{\hslash }{\widehat{S}}_z^B+\frac{g_N{\mu }_N{B}_0}{\hslash }{\widehat{I}}_z^A+a_A {\widehat{S}}^A·{\widehat{I}}^A-J\left(\frac{1}{2}+2{\widehat{S}}^A·{\widehat{S}}^B\right)$$

(25)

In addition, for coherent control, we did not consider the recombination reaction.

2.5. Parameter Global Optimization

The parameters used for control, such as the weighting parameter $q_i$. and the amplitude parameter $A$, have a strong influence on the fidelity of the desired quantum state $|f\rangle$ depending on their values. To obtain the desired state with high fidelity, it is necessary to use an appropriate combination of parameters. Therefore, we define the fidelity of the target state obtained by the one control as a function of $q_i$ and $A$, Equation (26), and use global optimization (ga function in MATLAB Global Optimization Toolbox [37,38,39]) to find the combination of parameters that maximizes the fidelity of $|f\rangle$.

$$f\left(q_i, A\right)=\langle f|\rho \left(t_f\right)\rangle$$

(26)

While global optimization requires iterative computations and is computationally expensive, it is more likely to obtain the desired parameters without being trapped by a local minimum due to its holistic search. In contrast, computations based on the local optimization theory do not require iterative computations, so the one-time computational cost is very low and practical control results can be obtained, but whether the obtained results are the best results cannot be determined in a single computation. Therefore, a global search of the parameters necessary for local optimization is a good method that combines the advantages of both local and global optimization since it can calculate optimal values while reducing the computational cost of each iteration of the search. Using this method, the optimal parameters for coherent control were obtained ($q_S, q_{T_{\pm }},q_{T_{+\alpha }, T_{-\beta }}, A$ described below, which was obtained by this method).

3. Results and Discussion

3.1. Anisotropic Reaction Control

Model RPs g-factors are set as $g_A=2.000, g_B=2.002$, and external magnetic field $B_0$ was fixed to 3.6 mT. The following anisotropic hyperfine tensor was used.

$$\widehat{A}\left(/\mathrm{mT}\right)=\left(\begin{array}{ccc}1.2& 0& 0\\ 0& 1.2& 0\\ 0& 0& 1.8 \end{array}\right)$$

(27)

The exchange interaction $J$ was fixed to 0.1 MHz. The initial state of RP is assumed to be singlet state.

First, we designed a waveform that provides a transfer of the radical pair from S-T₀ mixed states to $P_k$ states homogeneously. This control corresponds with keeping the RP population isotropic in the magnetic field direction by setting $q_k=1$ for all $k$. Figure 3 shows the results of the optimized RF field and decay of RPs with and without RF.

Figure 3 shows that 2 $\mathsf{\mu }$s after radical pair formation, the population of RPs, which is 38% in the case of no RF irradiation, is successfully maintained at 94% by RF irradiation. This is a result of the RF-induced transition of the radical pair from the singlet to states with high ${\mathrm{T}}_{+\mathsf{\alpha }}$ and ${\mathrm{T}}_{-\mathsf{\beta }}$ character, which are hardly involved in S-T₀ mixing, and the suppression of the recombination reaction.

We now examine the anisotropy of the RP population using the coordinates shown in the top panel of Figure 4. Figure 4 shows the RP population at $t$ = 2 $\mathsf{\mu }$s during RF irradiation and non-irradiation, indicating that the RP population is maintained isotropically during RF irradiation.

Figure 5 shows the results of a Fourier transform of the irradiated RF signals, and it can be confirmed that the frequency components corresponding to the transitions of radical B are dominant. This is because the transition to radical B, which has isotropic interactions, is advantageous for maintaining an isotropic population.

From now on, we consider the control of the RP population anisotropic with respect to the external static magnetic field direction ($z_L$). At first, we tried to keep the RPs with nearly parallel orientation with respect to the external field, i.e., $z_L\parallel z_{\mathit{M}}$. We set $q_k=5$ for $\beta$ = 0~50°, 130~180°, and $q_k=-3$ for $\beta$ = 60~120°. The negative value directly means a penalty for the undesired transition.

Figure 6 shows the results of the anisotropic control. The RP population decays by recombination for sub-ensembles $\beta =$ 60~120° and the component, for $\beta$ = 0~50°, 130~180° stayed longer as RPs. This is because the RF-induced transitions to the eigenstates with relatively high ${\mathrm{T}}_{+\mathsf{\alpha }}$ and ${\mathrm{T}}_{-\mathsf{\beta }}$ character for the sub-ensemble $\beta$ = 0~50°, 130~180° remarkably suppress the recombination efficiency, while the recombination was not suppressed for the sub-ensemble of $\beta$ = 60~120° because the transitions were not induced.

The Fourier transform of the optimized waveform is shown in Figure 7. In contrast to the isotropic control case, the algorithm for the optimized waveform has switched the target radical transition from radical B to radical A, because radical A has the frequency distribution of the transition due to the anisotropic hyperfine interaction.

Next, we tried to keep the RPs of the $\beta$ = 60~120° sub-ensemble longer. For this, the time weighting parameters were set to $q_k=-2.5$ for $\beta$ = 0~50°, 130~180° sub-ensemble and $q_k=5$ for the $\beta$ = 60~120° sub-ensemble.

In contrast to the previous result, RPs whose hyperfine tensors are relatively oriented perpendicular to the magnetic field direction are retained, while RPs whose hyperfine tensors are oriented parallel to the magnetic field direction are reduced, as shown in Figure 8. Figure 9 shows the optimized waveform and the result of the Fourier transform. The frequency components corresponding to the transitions of radical A tend to dominate, but the components corresponding to the transitions of radical B are increased compared to the previous control. In order to confirm the timing of the frequency components used in the control, we performed short-time Fourier transform, as shown in Figure 10.

Figure 10 shows that the frequency components of radical A are used for a long time in both the controls shown in Figure 6 and Figure 8, while the frequency components of radical B are used only in the early time (see around 100 MHz). In the early stage, the transition to radical B contributes to the total RPs being maintained for a longer time, thus contributing to the increase in the performance index $y\left(t\right)$. On the other hand, in the later part, isotropic transitions to radical B are avoided in order to avoid an increase in the sub-ensemble with penalized weighting parameters.

Our calculation can increase the performance index by selecting the RF frequency precisely over a long period of time, ~2 μs. In practice, however, most of the control is achieved relatively early. In the three examples shown in Figure 3, Figure 7 and Figure 9, the performance index y(t) at t = 0.25 μs is 97%, 79%, and 92% of y(t) at t = 2 μs, respectively.

3.2. Coherent Pathway Control

Now, we discuss coherent control. For coherent control, the hyperfine coupling constant $a_A$ was set to 15 MHz ($\approx$0.54 mT), the RP’s g-factors were set as $g_A, g_B=2.000$, and the exchange interaction $J$ was set to 0.1 MHz. Under this condition, only the hyperfine interaction contributes to S-T₀ mixing. The calculation was performed under external magnetic fields of 20 mT and 0.05 mT. As in the previous section, the initial state of the RP was set to the singlet state, and the calculation was performed so that the RP returned to the singlet state again after 1 $\mathsf{\mu }$s via the T_± state. To handle this control, the moving target was set as Equation (28).

$${|{\rho }_{T_{\pm } \to S}\rangle }_{ }=exp\left\{i{\widehat{\widehat{L}}}_0\left(t_f-t\right)\right\}[q_{T_{\pm } }{|T_{\pm }\rangle }_{ }+q_S|S\rangle ]$$

(28)

Figure 11 shows the result of coherent control under $B_0=20 \mathrm{mT}$.

In the high-field condition, the time evolution does not produce the probability of T_±. Therefore, the strategy of coherent control is as follows.

• S-T₀ mixing increases the population of T₀ states;
• The RF provides a transition from the T₀ state to the T_± state;
• The transition from T_± to T₀ is again triggered by the RF, timed to maximize the singlet at the target time.

The population of the singlet state after 1$\mathsf{\mu }$s without RF irradiation is less than 10% due to S-T₀ mixing. Optimized RF induces T₀ – T_± transitions around 0.2$\mathsf{\mu }$s, and the subsequent T_± − T₀ transition shifts the phase of S-T₀ mixing, eventually returning about 98% of the RP to the singlet state at 1 $\mathsf{\mu }$s.

Figure 12a,b show the power spectrum of the optimized RF and short-time Fourier transform, respectively. The optimized RF is mainly composed of radical B, which has no hyperfine interaction. The peak of this frequency component rises from 0.2 $\mathsf{\mu }$s and decays around 0.6 $\mathsf{\mu }$s, which is in good agreement with the time evolution of the spin state. From these results, we can confirm that the local control theory has attained the desired control by causing transitions to the radical B step by step in the control at a high magnetic field.

In the extremely low magnetic field condition ($B_0$ = 0.05 mT), the energy separation Zeeman splitting becomes small, and hyperfine interactions dominate. As a result, not only S-T₀ mixing but also other state mixings occur, as shown in Figure 13. The mixing between spin states looks more complicated than at 20 mT. In the model system, singlet states mix with ${\mathrm{T}}_{+\mathsf{\beta },-\mathsf{\alpha } }$, but ${\mathrm{T}}_{+\mathsf{\alpha },-\mathsf{\beta } }$ states remain isolated even in a low magnetic field. Only the RF field can make transitions to the isolated ${\mathrm{T}}_{+\mathsf{\alpha },-\mathsf{\beta } }$. Therefore, these states are useful as the intermediate quantum states for the coherent control of the RPs. To achieve more efficient population control, the moving target with weighting factors for ${\mathrm{T}}_{+\mathsf{\alpha },-\mathsf{\beta } }$ was employed.

$${|{\rho }_{T \to S}\rangle }_{ }=exp\left\{i{\widehat{\widehat{L}}}_0\left(t_f-t\right)\right\}\left[q_{T_{+\alpha }, T_{-\beta }}|T_{+\alpha ,-\beta }\rangle +q_S{|S\rangle }_{ }\right]$$

(29)

Figure 13 shows the result of coherent control under $B_0=0.05 \mathrm{mT}$. By optimized RF, RP’s population were stored to ${\mathrm{T}}_{+\mathsf{\alpha },-\mathsf{\beta } }$ states for 0.4 $\mathsf{\mu }$s to 0.9 $\mathsf{\mu }$s. Due to the transitions, the phase of S-T₀ mixing was changed around 0.5 $\mathsf{\mu }$s, and singlet fidelity was improved from 39% to 96%.

Figure 14 shows the power spectrum and the short-time Fourier transform of optimized RF shown Figure 13. In contrast to the high field case, the frequency component corresponding to the transition of radical A mainly constitutes optimized RF. This is considered because the Zeeman splitting under the geomagnetic field is so small that the RF prioritizes the control of transitions between states split by nuclear spins. Thus, under very low magnetic fields such as the geomagnetic field, the broad frequency component due to the hyperfine interaction of the system can affect the reaction rather than the Zeeman splitting of the external magnetic field [21], and the present control results suggest this.

3.3. Parameter Optimization

Figure 15 shows the evolution of the weighting parameters $\left\{q_i\right\}$ and amplitude parameter A and Singlet fidelity at the target time (1 $\mathsf{\mu }$s) during global optimization.

At generation zero, we set arbitrary weighting parameters and amplitude parameters. The genetic algorithm (ga function) of the MATLAB global optimization toolbox generated fifty sets of parameters {$q_S, q_{T_{\pm }}, A \mathrm{or} q_S, q_{T_{+\alpha }, T_{-\beta }}, A$} and calculated the locally optimized RF field. In fifty sets of the calculations, the maximum and minimum fidelities, $f_{max}$ and $f_{min}$, were plotted as shown in Figure 15. This process was repeated fifty times and we can confirm the convergence of the fidelities with respect to the parameter sets. As shown in Figure 15, the number of generations required for the convergence of the parameter sets is greater for the low-field control (Figure 15b) than for the high-field control (Figure 15a). The high field control is a simple problem of how to make a phase shift of S-T₀ mixing oscillating by a single frequency. In contrast, the control at a low magnetic field is more difficult because the coherent oscillations contain both the hyperfine coupling and the external magnetic field, and the value of the control parameter has a large effect on the fidelity of the singlet. Therefore, the high-field results converged early, while the low-field results required more generations to converge the parameters. Here, we define $\left[P_{max}\right]$ as the parameter combination $\left\{q_S, q_{T_{+\alpha }, T_{-\beta }}, A\right\}=\left\{0.93, 0.89, 1.1\right\}$, which provides the maximum value of singlet fidelity and $\left[P_{min}\right]$as the parameter combination $\left\{q_S, q_{T_{+\alpha }, T_{-\beta }}, A\right\}=\left\{0.91, 0.80, 0.46\right\}$, which provides the minimum value of singlet fidelity in low-magnetic-field control.

Figure 16 shows the time evolution of the spin system when $\left[P_{max}\right]$ and $\left[P_{min}\right]$ are used as control parameters. The transition to the ${\mathrm{T}}_{+\mathsf{\alpha },-\mathsf{\beta } }$ state is slower in the control using $\left[P_{min}\right]$ compared to the control using $\left[P_{max}\right]$, and the target time is reached before the RP populations return to singlet state again. This reflects the importance of the optimal combination of parameters, especially when coherent mixing is significant. Therefore, the global search for local optimization parameters used in this study is a method that can compensate for the weakness of the local optimization method and obtain better control results.

4. Concluding Remarks

In this paper, we focused on the weighting parameters used to control local optimization and demonstrated sub-ensemble and coherent control for a model system of RPs. By using appropriate weighting parameters, we achieved control over the selection of anisotropic sub-ensembles and the design of transition paths, albeit in a simple system.

The anisotropic response of RPs to magnetic fields is of interest from the perspective of avian magnetic compasses, and the possibility of oscillating magnetic fields affecting their dynamics is related to test experiments on animal behavior. The anisotropic induction of RP by RF irradiation presented in this paper, together with the fact that broadband RF irradiation disrupts the orientation of migratory birds [15,20,21], implies that RF irradiation of a certain waveform may be able to induce bird orientation. However, in practice, it would be difficult to apply this method to animal experiments in geomagnetic fields.

Coherent control reveals the RF waveform that maximizes the singlet state at a specific time. When singlet recombination is slow, this coherent control has little effect on the actual reaction yield. However, it has recently been reported that the coherent mixing of RPs has been directly observed by pump-push spectroscopy [40]. In this technique, the spin-selective chemical reaction is enhanced by a push laser pulse instantaneously. Therefore, we could reflect the coherent spin control into the reaction kinetics in an efficient way. The combination of this method with the coherent control by AWG-RYDMR could open up a new methodology for chemical reaction control. This calculation and the experimental verification are in progress.

Although our method uses an almost completely known spin Hamiltonian to control the chemistry of radical pairs, an inverse problem is the possibility of searching for spin Hamiltonians from this AWG-based control of chemical reactions. The intention is to optimize the AWG signal by experimental feedback and, inversely, to obtain information about the radical pairs. However, it is challenging to obtain information in this way that cannot be obtained by conventional EPR spectroscopy. Rather, such a method may be useful for obtaining kinetic information, such as recombination reactivity, or when several radical pairs with different reactivities are mixed.