Sensing of magnetic field effects in radical-pair reactions using a quantum sensor

Magnetic field effects (MFE) in certain chemical reactions have been well established in the last five decades and are attributed to the evolution of transient radical-pairs whose spin dynamics are determined by local and external magnetic fields. The majority of existing experimental techniques used to probe these reactions only provide ensemble averaged reaction parameters and spin chemistry, hindering the observation of the potential presence of quantum coherent phenomena at the single molecule scale. Here, considering a single nitrogen vacancy (NV) centre as quantum sensor, we investigate the prospects and requirements for detection of MFEs on the spin dynamics of radical-pairs at the scale of single and small ensemble of molecules. We employ elaborate and realistic models of radical-pairs, considering its coupling to the local spin environment and the sensor. For two model systems, we derive signals of MFE detectable even in the weak coupling regime between radical-pair and NV quantum sensor, and observe that the dynamics of certain populations, as well as coherence elements, of the density matrix of the radical pair are directly detectable. Our investigations will provide important guidelines for potential detection of spin chemistry of bio-molecules at the single molecule scale, required to witness the hypothesised importance of quantum coherence in biological processes.


I. INTRODUCTION
Investigating the fundamental role of spin interactions in magnetic field effects (MFE) in chemical reactions has a long history [1][2][3].The study of spin-chemical effects provides important insights about structure, kinetics and magnetic properties of transient intermediate chemical species [4,5] and are explored in various interdisciplinary applications including sensitivity enhancement in nuclear magnetic resonance (NMR) [6], quantum computing [7], avian magnetic compass [8], and solar energy conversion kinetics in photosynthetic systems [9].
It is astounding that spin interactions can have decisive effects on the fate of chemical reactions as the energy of spin transitions typically is orders of magnitude smaller than the thermal energy [4].However, spin dependent MFEs in chemical systems, also known as the radical-pair mechanism (RPM) [2,8,10], rely on the creation of transient paramagnetic species in a non-equilibrium state called radicals, which are chemical species with an odd number of electrons.In the RPM, the radical pair (RP) is spatially separated but in a spincorrelated state and the recombination of the radicals back to the molecular precursor state is spin selective.The influence of an external magnetic field then occurs in terms of a modulation of the spin dynamics and consequently an alteration in the yield of products formed from the various spin states [11,12].
The formulation of the RPM started in the late 1960s to explain non-equilibrium magnetic resonance spectra of chemical reactions of organic molecules [13][14][15][16][17][18][19], while the recent interest is fuelled by investigations of the RPM as the most plausible mechanism for MFE in biological reaction kinetics [20,21].These include flavoproteins related to DNA

Radical-pair
Diamond quantum sensor FIG. 1. Simplified illustration for sensing of MFEs of a RP reaction using a NV quantum sensor in diamond.A RP on a host molecule is formed in a spin-correlated state (yellow).Subject to local and external fields, the distinct RP evolution induces a signal that is detectable by a nearby NV centre (red) using appropriate sensing sequences.
photolyases, the involvement of cryptochromes in circadian rhythms [22,23] and their proposed role in animal magnetoreception [24].MFEs have been recorded in cryptochromes [11] and seem to fulfil the structural and dynamical requirements of the RPM [8].Considerable interest exists also in the role of radical pairs in chemical kinetics in photosynthetic reaction centres [9,25,26].Flavoenzymes [27] -Flavin-based enzymes -are responsible for catalytic functions in diverse biological reactions [28] and the involvement of various RPs in these reactions is debated in recent discussions [29,30].
Existing experimental techniques for in vitro probing of MFE in RP dynamics in chemical reactions rely on averaged signals collected from a large ensemble of molecules.These techniques include time-resolved electron paramagnetic resonance (TREPR) spectroscopy [31,32], and optical methods based on absorption [33][34][35][36] and fluorescence detection [37][38][39][40][41]. Studies of RP reactions conducted using these techniques can only provide ensemble averaged information of spin dynamics along with requiring a large quantity (few microliters) of potentially precious biological samples.It is therefore imperative to instead consider single-molecule detection techniques in order to reveal potential quantum coherent spin evolution otherwise hidden by ensemble averaging [42][43][44].Specifically the single negatively charged nitrogen-vacancy (NV) center in diamond [45] has attracted significant interest as potential candidate for the detection of spin-chemical effects of RP reactions at the single molecule scale [43,44].This is achieved due to the excellent bio-compatibility of diamond [46,47], the attainable nanometer scale spatial resolution and remarkable sensitivity to external electromagnetic fields of a single NV center [48][49][50][51][52].
In the present work, we discuss a realistic avenue for the detection of MFE associated with RPs at the scale of single and small ensembles (≤ 100) of molecules using a single NV center in diamond (Fig. 1).We investigate both the weak and the strong coupling regime between the NV center and RPs.Compared to the simplified model employed by Finkler et al. [44] neglecting important spin interactions and considering only the strong coupling regime, we here use an elaborate model governing realistic RP spin dynamics.We include up to three nuclear spins per radical with maximum anisotropic hyperfine coupling along with dipolar and exchange coupling between the radicals.Applying standard sensing protocols, we derive measurable signals received by an NV centre and quantify magnetic field dependent RP dynamics.We show that in the weak coupling regime, the signal generated by an RP on a single-molecule is comfortably within the sensitivity limit of state-of-the-art shallow NV centers and dynamics of certain populations, as well as coherences, of the RP density matrix is directly accessible.Non-trivial features of RP spin dynamics can be observed depending on direction and magnitude of a bias magnetic field.Further, we observe that signal features tend to average out when we consider a small ensemble of RPs, highlighting the importance of single-molecule detection.In the strong coupling regime, we find that although there is an opportunity to probe various RP spin state populations individually, detection becomes challenging especially in larger bio-molecules due to the increased number of unavoidable and spectrally indistinguishable hyperfine interactions within the RP.
The article is organized as follows.In section II, we introduce the theoretical framework for the RP and interaction with the NV center.In section III, we describe the numerical simulations performed with realistic RP systems and present the main results.We conclude in section IV and discuss possible future directions.

II. DESCRIPTION OF THE MODEL
The general situation considered for the model is illustrated in Fig. 2 (a).Molecules hosting RPs are distributed above the diamond surface and an NV center is situated at a depth d from the surface.We choose to work in the NV frame (shown in the figure) and the total Hamiltonian of the whole system in Reduction (RP creation) " When RP spin dynamics is probed at very low magnetic fields (close to Earth's), a gradient is required for the magnetic field in the z direction to ensure the two level approximation of NV center energy levels is valid.(b) Simplified description of the RPM.The two factors, (i) magnetic field dependent inter-conversion between the singlet state and triplet states, and (ii) spin dependent product yield, make the RPM a plausible mechanism behind MFE in RP reactions.

Recombination
this frame is given as where H NV , H RP , and H C are the Hamiltonians governing the dynamics of the NV center, the RP, and coupling between the NV center and the RP, respectively.We consider the C 3v symmetric, negatively charged NV − center (hereafter called NV) that consists of ground ( 3 A, total spin s = 1), metastable ( 1 A, s = 0) and excited ( 3 E, s = 1) states [45].The ground state triplet {|0 , | ± 1 |} is split at zero magnetic field with D ZFS = 2.87 GHz [48].The | + 1 and | − 1 can be further split by application of an external magnetic field, thus making all the three states accessible by application of appropriate control fields.Upon optical excitation typically with 532nm light, the NV center might relax via a spin-dependent inter-system crossing [53], a process that allows for initialization in |0 and spin-state dependent fluorescence contrast between |0 and the | ± 1 states.These properties combined with an external magnetic field can be exploited to isolate |0 ↔ | + 1 or |0 ↔ | − 1 as a two-level system with effective spin s = 1/2.In addition, the NV center interacts with nuclear spins, predominantly (abundance of 99.6% [54]) the intrinsic 14 N with total spin s = 1, hyperfine coupling tensor A A A N and quadrupolar coupling Q N −5.01 MHz [54], causing further splitting of the | ± 1 states.Due to the axial symmetry of the NV center, A A A N can be expressed in a diagonal form in the principal axis system of the NV axis with diagonal elements [54].Taking the above details into account, the relevant part of the NV center Hamiltonian can be written as where γ e is the gyromagnetic ratio of an electron, J = [J x , J y , J z ] and I N = [I Nx , I Ny , I Nz ] are, respectively, spinoperators of the NV center and 14 N nuclear spin, and B 0 = [B 0x , B 0y , B 0z ] is the vector of the externally applied magnetic field on the NV center.The relatively large zerofield splitting allows to make the secular approximation where we can ignore all terms involving electron-nuclear spin flip-flops (containing A ⊥ N ).Further, we assume B 0z >> B 0x and B 0y and consequently choose |0 ↔ | + 1 (hereafter denoted |1 ) as two-level system.Now the NV center Hamiltonian simplifies to The RP consists of two radicals each containing an unpaired electron.A simplified reaction for the creation and recombination of radicals is shown in Fig. 2(b).The radicals can be generated either by electron transfer or chemical bond breaking from the oxidized state of the molecule.The reduction can be facilitated by various mechanisms, for example, by photoexcitation with light of appropriate wavelength (reduction of flavin based systems [55]) or by chemical means (Haber-Weiss reaction [56]).Depending on the internal molecular dynamics, the unpaired electrons in the radicals are born in either singlet (s = 0) where | ↑ and | ↓ denote the eigenstate of the Pauli z−matrix.After creation, the singlet and triplet states interconvert among each other under the RP Hamiltonian, which for the RP situated at location (r, ζ = {α, β}, see Fig. 2(a)) is given by [57] where S 1 and S 2 are the spin-operators of the unpaired electrons of the 1 st and 2 nd radical, respectively, and B(ζ) denotes the external magnetic field applied on the RP.The radicals couple to each other via exchange interaction (coupling constant J ex (ζ)) and dipolar interaction (coupling tensor D D D(ζ)).
Each unpaired electron is further surrounded by a set of nuclei in the radical and interact with them via hyperfine coupling.The I 1i (A A A 1i (r)) and I 2 j (A A A 2 j (r)) are spin-operators (hyperfine coupling tensors) of the i th nuclei coupled to the unpaired electron in the 1 st and j th nuclei coupled to the unpaired electron in the 2 nd radical, respectively.The magnitude of hyperfine coupling is typically in the range of 2.8 -28 MHz [58] in organic radicals, whereas the strength of the dipolar and exchange interaction is dependent on the separation between the radicals [59].Usually the dipolar and hyperfine coupling tensors are simulated or measured in a coordinate frame associated with the molecule that hosts the RP, hereafter referred to as the RP frame.The NV and RP frames are related by a rotation matrix , where H RP is the RP Hamiltonian in the RP frame.
To gain more insight into the inter-conversion dynamics of spin states of radicals, it is useful to work in the basis spanned by singlet and triplet states for the unpaired electrons.Using this choice of basis, we can divide The first part is diagonal in the singlet-triplet basis for unpaired electrons irrespective of the basis chosen for the nuclei where is the secular part of the dipolar coupling and r RP is the distance between the radicals.The second part H nd RP (ζ) contains terms corresponding to the transverse magnetic field, non-secular dipolar coupling contribution and hyperfine interactions with the surrounding nuclei.The singlet and triplet states are anti-symmetric and symmetric under spin exchange, respectively, and the evolution under H nd RP (ζ) breaks this symmetry which results in interconversion between these states.
Due to interaction with the surrounding environment, after the creation, along with the inter-conversion dynamics, the spin-correlated RP recombine back to the equilibrium state.The rates of recombination depend on whether the RP is in the singlet (rate constant k S ) or the triplet state (rate constant k T ).Products are thus formed with spin-state dependent yield and can be altered by the application of a suitable magnetic field, a process that is called the radical-pair mechanism (RPM) [8].The recombination dynamics can be modeled by treating the RP as an open quantum system with Lindbladian Where N = N 1 + N 2 is total number of nuclei in the RP and {.} denotes the anticommutator.
The electron spin of the NV center and the RP is coupled via dipolar interaction and the corresponding Hamiltonian can be written as where cz represents the effective RP-NV coupling (Fig. 3) and can be determined using techniques developed in recent works for threedimensional molecular localization using the NV center [60][61][62][63].Note that the coupling Hamiltonian can not be approximated to contain only the J z Sz term [44,64] as the strength of cross terms (J z Sx and J z Sy ) may not be smaller compared to various other parameters in H RP at low external magnetic fields.Now, after establishing the general framework, the protocol to study RP dynamics using the NV center is as follows: 1. Prepare the NV center and the RP in a known initial state ρ(0) = ρ NV (0) ⊗ ρ RP (0).The initial state of the unpaired electrons in the RP is either the |S 0 or |T 0 state as described before, while the nuclei start in the maximally mixed state I ⊗N /2 N .The initial state of the NV center can be controlled and depends on the applied sensing scheme.
2. Solve the quantum dynamics in the NV frame using the master equation: Since the NV center has the capability of detecting magnetic signal directly originating from inter-conversion dynamics of the RP instead of relying on the product yield based detection, we can simplify the Lindbla- ) , ρ(t) [57].Further, the rate k may be allowed to include decoherence of the electron spins of NV center as well because its typical time-scale is similar to recombination rates in organic RPs which is in the order of tens of µs.
3. Trace out the RP to calculate the state of the NV center at any time t: ρ NV (t) = Tr RP [ρ(t)] and investigate signatures of the RP evolution on the states of the NV center.
Depending on r, θ and φ, the coupling between an NV center and an RP varies (Fig. 3).The resulting NV-RP dynamics can be divided in two dynamical regimes: Weak coupling : where Γ = 1 π T 2 and T 2 is the dephasing time of the NV center.T 2 depends on factors including applied sensing sequence and magnetic field [65], proximity to surface impurities, presence of paramagnetic defects and nuclear spins around the NV center [45].Typical values of T 2 range from a few microseconds to a few milliseconds depending mainly on the applied sensing sequence, the isotopic purity of the diamond sample used [66] and the surface termination [67], making both weak and strong coupling regimes achievable in practice.

III. SIMULATIONS
In this section we present details and results of the simulations.We consider two RP systems, respectively with isotropic and anisotropic hyperfine coupling tensors.We investigate MFEs in RP systems in the weak as well as strong coupling regimes and discuss appropriate detection strategies.In the weak coupling regime, we study the MFE as a function of strength and direction of magnetic field due to a single as well as a small ensemble of RP.For the strong coupling regime, we point out potential challenges in the regime of single molecule detection that arise for larger bio-molecules.

A. RP systems
The two RP systems we investigate are: 1. Flavin adenine dinucleotide -Tryptophan (FAD •− -TrPH •+ ) with anisotropic hyperfine couplings, which has been subjected to various MFE studies using absorption spectroscopy in recent years [11].This system is considered to be the RP behind avian magnetoreception [8].

B. Weak coupling regime
In the weak coupling regime, no observable splitting of the NV |1 level occurs upon interaction with the RP and the effect of spin-dynamics of the RP on the NV center dynamics then appears only as a time-dependent classical magnetic field.The coupling Hamiltonian simplifies to where The detectable magnetic signal generated from the RP situated at a location (r, ζ) can now be calculated as (in units of Tesla): where i = x, y, z, ξ is the number density of RPs, and (r 2 − r 1 ) is the sensing radius of the NV center (see Fig. 2(a)).The maximum signal is generated when R ζ = I , where I is the identity matrix, i.e. when the RP frame is aligned with the NV frame: The possibility of achieving such a condition is discussed in Appendix VI A.
A simple protocol to detect the above signal is shown in Fig. 2(d).First, we initialize the NV center in |0 and the RP in either the |S 0 or |T 0 state.An appropriate sensing sequence is then applied on the NV.The choice depends on the Fourier spectrum of the detectable signal, denoted by X i (ω, ζ) here.Finally, a projective measurement is made on the NV to yield information about the accumulated phase, which is proportional to X i (t, ζ), acquired during the sensing sequence.
As evident from Eq. ( 11), only specific combinations of RP density matrix elements corresponding to {S ji , j = 1, 2 and i = x, y, z} contribute to X(t) and are directly measurable by the NV center.Explicitly, in the singlet-triplet basis, the directly measurable elements are: which measures the magnetization corresponding to the population difference of outer triplet states and which measure the magnetization corresponding coherences between the outer and central triplet states.Thus the detection of magnetization corresponding to elements S 1x + S 2x and S 1y + S 2y can reveal involvement of quantum coherent phenomenon in RP dynamics in singlet-triplet basis.The full state tomography of density matrix of RP, however, requires conversion of various elements onto directly detectable ones (Eq.14 and 15) and is beyond the scope of this work.
We consider RPs to be statistically distributed above the surface of the diamond sample.The coupling strength of an RP outside the sensing radius of three-four times the depth d of the NV centre (cf.Fig. 2(a)) decreases by an order of magnitude compared to RPs on the surface.The number of RPs falling in the sensing volume depends on the size of the host protein.For the case of proteins with many amino acids and large molecular weights (>500kDa, the average radius > 5nm [70]), only one of them might fall in the sensing volume with high probability.On the other hand, for the case of smaller proteins (<50kDa, the average radius being 2-3 nm), tens of RPs can be accommodated in the sensing volume.The following two subsections discuss how the generated signal behaves in these two cases.is of the order of the observed spin relaxation time in behavioural studies in migratory birds [71].To calculate the total signal received (or phase accumulated by the NV center) at a given strength and direction of applied magnetic field, we define the time integrated signal per second as In Fig. 5, we show the strength of the magnetic signal received by the NV from an RP on a single molecule at a distance r = 10nm.For such shallow NV centres, the achievable sensitivity is of the order of a few nT/ √ Hz to a slowly alternating signal [64,72].For the case of φ = 0, as we assume here, the detectable signal is only generated by X I x and X I z , which, respectively, corresponds to the dynamics of RP density matrix elements S 1x + S 2x and S 1z + S 2z .
To determine the suitable NV sensing sequence to probe the generated signal, Fig. 5(a) plots the temporal evolution of X z (t) for a bias magnetic field that maximizes the time integrated signal X I z (see Fig. 5(c)).The unpaired electrons in the RP are initialized in the singlet-state whose population decays with time as it is converted into other density matrix elements under the evolution of H RP , thus resulting in the build-up of X I z .The generated signal (hundreds of nT) is within the sensitivity achieved for < 10nm deep NV centres.The corresponding spectrum is shown in Fig. 5(b).Because we are interested in low frequency dynamics (kHz to tens of MHz) resulting from the dipolar and hyperfine coupling terms (typically 1-10mT) in H RP , the maximum component of the signal in the frequency domain is limited to 500 MHz.The cut-off frequency of these oscillations observed is in the order of 10MHz.Sequences including Ramsey, spin-echo, and dynamical decoupling are thus suitable for probing such RP dynamics [65].
Variation with strength of the magnetic field − To study the dependence of the signal on the strength of the applied magnetic field, we set θ = 0, i.e., B 0x = 0 and B 0y = 0 so that there is no spin mixing of the NV center states keeping the two-level approximation valid and B z is varied from 0 to 50mT.The dependence of the time-integrated signal X I z X I x = X I y = 0 in this case as θ = 0 is plotted as a function of B z in Fig. 5(c).
The features of the time-integrated signal X I z can be qualitatively explained by the dynamics of the RP under H nd RP in singlet-triplet basis as shown in Fig. 5(d).At B 0z ≈ 0, all three triplet states are almost degenerate and the energy gap between |S 0 and |T 0 states is proportional to J ex and D s .With the unpaired electrons in the RP initialized in the singlet state, singlet-triplet oscillations occur due to evolution of the RP under H nd RP containing both a non-secular dipolar coupling contribution and hyperfine interactions with the surrounding nuclei (Eq.4).However, for small B 0z , the RP only generates a very low signal X I z as outer triplet states (|T ± ) remain close to degeneracy.As B 0z increases, the outer triplet states move apart and give rise to additional manifolds for singlet-triplet mixing [73,74].This opens the possibility of information (population and coherence) transfer among |S 0 and |T 0 states to only one of the outer triplet states, creating a population imbalance between outer triplet states and giving rise to an increase in X I z .At a certain B 0z value, there is a maximum transfer between singlet-triplet oscillations to population difference of outer triplet states, giving rise to the peak-like features in X I z (for example the peak at ∼ 1.16 mT and ∼ 0.25 mT, respectively, for FAD •− -TrPH •+ RP and PY •− -DMA •+ RP).This phenomenon, which is closely related to low field effects (LFE) in the context of singlet yield [75], occurs due to level crossing between |S 0 and either of the outer triplet [76] states.A further increase in B 0z results in energetic isolation of the outer triplet states which causes steady decrease in the probability of information transfer between them, and |S 0 and |T 0 states, eventually vanishing at very high magnetic fields Variation with direction of the magnetic field − In Fig. 5(e), we investigate the signal dependence on the direction of the applied magnetic field.Here we restrict the magnetic field to the XZ plane by assuming φ = 0 as the following discussion holds for the general case.Although B 0x is nonzero for certain values of θ, its magnitude should be less than that required for mixing of the |1 and | − 1 levels of the NV center.We fix the magnetic field strength to ∼ 1.16 mT and ∼ 0.25 mT, respectively, for FAD •− -TrPH •+ RP and PY •− -DMA •+ RP where the maximum MFE is expected based on the results shown in Fig. 5(c).The dependence of the time-integrated signal on θ is plotted in Fig. 5(e).The shape of the signal received X I x and X I z (X I y = 0 since φ = 0) is dominated by the shape of D cx and D cz , respectively.To reveal the true dynamics of the RP with θ, we calculate the normalized signals instead in Fig. 5(f).
Once again, using the level dynamics under H nd RP in the singlet-triplet basis, a qualitative explanation of various features of the generated signals is possible using the following arguments: (i) for a given magnetic field, θ close to 90 • essentially means presence of an extra channel for singlet-triplet mixing due to nonzero B 0x and B 0y and (ii) at θ close to 90 • , the outer triplet states are almost degenerate with maximum splitting occurring at θ close to 0 • or 180 • .The shape of X I x is a direct consequence of argument (i) as nonzero transverse fields close to 90 • increase the probability of singlet-triplet as well as triplet-triplet mixing.The signal X I z is a consequence of argument (ii) as higher splitting of outer triplet states opens new manifolds for information transfer between them and |S 0 and |T 0 states as described earlier.
For parameter values chosen here and described above, we do not observe the recently discovered spike feature [58] in the generated signals X I z and X I x presumably due to the high value of magnetic field used, and the non-zero exchange and dipolar coupling in the H RP [77].However, the spike feature can be observed when we consider a simpler RP system with only one nuclear spin in each RP in earth's magnetic field as summarized in Appendix VI B.

Ensemble of RPs
Assuming the minimum size of the protein to be 3-4 nm (small protein with molecular weight < 50kDa [70]), ≈ 100 RPs (ξ = 5 × 10 −2 /nm 3 ) might fall in the NV sensing volume.In Fig. 6 we plot the mean and variance of the generated signal as a function of strength and direction of applied magnetic field for an ensemble of RPs.We consider two types of ensembles: (i) All the RPs are oriented along the direction of the NV axis, i.e., R ζ = I (solid lines in Fig .6), and (ii) all RPs are randomly oriented with respect to the NV axis (dashed lines in Fig. 6 , respectively, are rotation matrices about the x, y and z axes with randomly chosen Euler angles α, β, and γ.For both cases, the distance between RPs and the NV were sampled from a uniform distribution in the range 5 − 20nm. For the uniformly oriented ensemble (i), the mean of the generated signal is obviously amplified in contrast to the randomly oriented ensemble in (ii) where the features due to the The variance of the generated signal as function of strength of the magnetic field is within the sensitivity range of a shallow NV center for both types of ensembles and the LFE behaviour akin to the single-molecule case (Fig. 5 (c)) is still observable after averaging.On the contrary, although there is an observable signal when θ is varied, the shape of the variation is erratic with no similarity to the single molecule case (Fig. 5 (f)).Therefore, although the ensemble (ii) is simply realized by drop casting a molecule sample on the surface of the diamond, detecting a MFE in this arrangement is more challenging.The situation might be compared to absorption or EPR studies where the size of the probed ensemble is very large (micro litres of sample) and hence MFE detection relies on statistical signals generated from a very small number of molecules.This limits the detection of MFE to a few percent [35] along with causing wastage of precious biological samples.Hence, it becomes important to appropriately posi-tion the proteins in the desired orientation with respect to the NV center to maximize the probability of successful detection of an MFE (Appendix VI A).

C. Strong coupling regime
In the strong coupling regime, the energy levels of the NV and RP mix, giving rise to a level structure as illustrated in Fig. 7. Due to the s = 1 structure of the NV center, the RP components of eigenstates in |0 and |1 are quantized along different axes determined by the strength of the coupling tensor.The total Hamiltonian in this regime can be reduced to The magnitude of the spin resonance peaks is proportional to the population difference of |ψ n and |ψ n states.Taking the contribution from all RPs into account, it can be calculated as where P ψ i = |ψ i ψ i | is the projector corresponding to the state |ψ i .There are a couple of important details about the detection of RP dynamics in this regime which should be pointed out: (ii) As the size of the bio-molecule increases, the number of unavoidable hyperfine interactions also grows.As a consequence, the spectral features will start to overlap, especially when the applied magnetic fields are comparable to the hyperfine couplings in H RP (ζ) (see Fig. 7 (c) for FAD •− -TrPH •+ ).In this case it becomes challenging to access the various electronic transitions of the RP separately.One could consider decoupling the nuclear spin bath [78,79] during the readout time by driving the RP at a frequency larger than the strongest hyperfine coupling, however the short readout duration (typically ≈ 300ns for NV ceneter [45,80]) will render the averaging inefficient.
In the light of the second argument, although the detection of RP dynamics in the strong coupling regime may be challenging, it offers an opportunity to track various populations of the RP density matrix individually and the possibility of quantum control of the RP using an NV center [44].

IV. CONCLUSION AND OUTLOOK
In summary, we analysed the realistic prospects of MFE detection at the single molecule scale in RP reactions in biomolecules using a single NV center in diamond.To realistically describe the bio-molecular spin dynamics, our RP model includes three nuclear spins per radical that have the largest hyperfine couplings.Including dipolar and exchange coupling between the two electrons forming the radical pair, the RP spin dynamics becomes correlated and can not be simulated by treating each of the radicals independently [43].Depending on distance and orientation between NV and RP, two coupling regimes of dynamics can be considered requiring accordingly tailored NV sensing approaches.We find that in the weak coupling regime, RP dynamics can be seen as a classical magnetic field by the NV.This regime is most suitable for large bio-molecules with unavoidable and spectrally indistinguishable spin interactions.The signal generated even from a single RP is well within the sensitivity achievable by state of the art NV centers and distinct features of RP spin dynamics are thus observable at the single molecule and a small ensemble (≤ 100 molecules) scale.
The simulations performed in our work can be further improved by incorporating more nuclear spins and using the recently developed coherent state sampling method [81,82] for efficient computation of spin dynamics.Adding to the recent theoretical works [42][43][44], we expect our analysis to pave the way towards experimental detection of MFE in bio-molecules at the single molecule level and unravel the importance of quantum coherent effects in biochemical processes [20,21].

V. ACKNOWLEDGMENTS
We acknowledge financial support from the Novo Nordisk Foundation through the projects bioQ, bio-mag, and Quant-BioEng, the Danish National Research Foundation (DNRF) through the center for Macroscopic Quantum States (bigQ, Grant No. DNRF0142), the Villum Foundation through the project bioCompass, and the EMPIR program co-financed by the Participating States and from the European Union's Horizon 2020 research and innovation programme via the QADeT project (Grant No. 20IND05).

VI. APPENDIX A. Orientation of RP with respect to NV
As described in the main text, the orientation of the molecule hosting the RP is crucial for obtaining maximum signal.Recent studies involving investigation of structure and motion of single proteins using NV centers have employed either statistical placement of proteins via diamond surface functionalization [64,83,84] or controlled alignment via attaching the protein to a solid host and then using a nanopositioning system, for example a tuning fork scanning probe of an atomic force microscope [85].The latter method is preferred because, in the former case, there are two main disadvantages: (i) The RP might end up in an orientation where one or more coupling parameters (D ci , i = x, y, z) are zero, thereby reducing the coupling to the NV center (ii) The functionalization of the diamond surface may not be possible for all types of proteins.FIG. 8.The generated signal X I x and X I z for various types of hyperfine coupling tensors for the model system.

B. Spike features in the generated signal
Here we consider a simple RP model system with only one nuclear spin in each radical, with hyperfine tensors A A A 1 1 1 and A A A 2 2 2 , to analyze which parameter ranges give rise to spike-like features of the generated signal as a function of the direction of the applied magnetic field.This feature has been suggested to be behind the precision of an avian magnetic compass [58].
In Fig 8, we plot the generated signal X I x and X I z , respectively, in (a) and (b) for various types of hyperfine coupling tensors chosen by varying principal axis components set (A xx , A yy , A zz ).For the anisotropic hyperfine coupling case, the principal axis system is that of the N 5 in FAD •− [58].A qualitative description of the seen behaviour can be provided by looking at the symmetry of the RP Hamiltonian.In the absence of dipolar coupling, there are three interactions in the RP Hamiltonian, namely, Zeeman, exchange and hyperfine.Zeeman (same magnetic field on all the spins) and exchange interactions are symmetrical with respect to electronic spin exchange while hyperfine interaction is symmetry breaking in general.When the hyperfine interaction is also assumed to be isotropic, the full Hamiltonian is symmetric, and singlet and triplet are eigenstates.As a result, there is no singlet-triplet FIG. 9.The generated signal X I x and X I z , and yield for various values of exchange coupling for the model system.Here A A 1 1 1 A A A 2 2 2 and chosen according to the case 'axial 3 ' (see caption of Fig. 8).
oscillations possible for isotropic hyperfine coupling as the sub-spaces of different symmetry remain unconnected and as consequence, X I x and X I z remains zero.As the anisotropy is included in the Hamiltonian, a non-zero signal is generated along with appearance of spike-features with the most prominent around θ = 90 • .Again, the amplitude of this spike is within the sensitivity achieved for < 10nm deep NV centres.Now we analyze the behaviour of the generated signal and spike as a of various parameters of the nian.
The amplitude of the spike increases and it becomes narrower as the principal components A xx and A yy are increased, similar to the behaviour observed in the singlet yield in Ref. [58].This observation hints that the origin of the spike might also be similar, i.e, an avoided crossing of energy levels as function of direction of the magnetic field.However, further investigations are required to draw firm conclusions.
In and triplet states, and reduction (increase) in the signal (singlet yield) is observed.The spike feature in the generated signal at θ = 90 • is most prominent when J ex is comparable to the A xx and A yy principal component of the hyperfine coupling.This behaviour hints that the spike in the generated signal depends on exchange interaction (and also dipolar in general) as well as anisotropic hyperfine interaction.This is in contrast to the singlet yield where the spike depends strongly on hyperfine anisotropy as it disappears even for small J ex , comparable to A xx and A yy parameters [77].
In Fig 10, we show the generated signal X I x and X I z along with the singlet yield [86] (c) with increasing value of RP lifetime τ.A contrasting behaviour is observed in the generated signal and singlet yield in addition to the former being relatively less sensitive to the lifetime of the radical.This observation seems to suggest that unpaired electrons prefer to stay in the singlet state compared to other populations and coherences of the RP.However, again, further studies are required in this direction to draw concrete conclusions which are beyond purview of current work.

FIG. 2 .
FIG. 2. (a)General setting for detection of RP spin dynamics using an NV center situated at depth d under diamond surface.The light blue hemisphere shows the sensing volume of the NV center and a RP situated in this volume contribute to the detectable signal by the NV center.Without loss of generality, we choose to work in the axis system where normal to the diamond surface is aligned along the NV axis ([1 1 1] crystal axis), named 'NV frame' here.The center of the coordinate system associated with the molecule that hosts a RP is situated at location ( r, ζ = {α, β}) with respect to the NV center, named 'RP frame'.The two axis systems are related by the rotation matrix R ζ which is different for each RP.The magnetic field is applied in direction (θ, φ).When RP spin dynamics is probed at very low magnetic fields (close to Earth's), a gradient is required for the magnetic field in the z direction to ensure the two level approximation of NV center energy levels is valid.(b) Simplified description of the RPM.The two factors, (i) magnetic field dependent inter-conversion between the singlet state and triplet states, and (ii) spin dependent product yield, make the RPM a plausible mechanism behind MFE in RP reactions.
where D D D c c c is the dipolar coupling tensor in the NV frame which depends on the direction of the applied magnetic field (θ, φ) and the NV centre distance to the RP.Under the secular approximation |D ZFS + γ e B z | >> |D D D c c c | and assumption |B z | >> {|B x |, |B y |}, the coupling Hamiltonian simplifies to

FIG. 3 .
FIG. 3. The effective coupling strength g eff /2π = D r D 2 cx + D 2 cz /π for different relative distances r between the RP and NV sensor and angles θ of the magnetic field (φ = 0 is assumed).

1 .FIG. 4 .
FIG. 4. Energy level structure of an NV center in the weak coupling regime.The RP causes a relative phases between |0 and |1 states of NV center that is proportional to H C (ζ,t).The sequences to detect the RP dynamics in each regime are shown in (b).The RP is born in |S 0 or |T 0 state and the NV is initialized in |0 state.Then an appropriate sensing sequence is applied which is chosen based on the frequency profile of X i (ω, ζ).During sensing, the NV center acquires a phase which contains the information about the RP dynamics which is then read out by a projective measurement.

FIG. 5 .
FIG. 5. Signal received by the NV center from a single RP molecule of model systems FAD •− -TrPH •+ (solid lines) and PY •− -DMA •+ RP) (dashed lines).Time (a) and frequency (b) profile of the signal at magnetic field B 0z = 1.16mT for FAD •− -TrPH •+ RP, where maximum of time integrated signal X I z and X I x is obtained as plotted in (c).θ = 0 • is assumed in (a), (b), and (c).(d) Level dynamics of the RP in singlet-triplet basis as function of the strength of the applied magnetic field.Dependence of the integrated signal on the angle of the applied magnetic field with (e) and without (f) θ correction at magnetic field |B 0 | = 1.16mT and |B 0 | = 0.25mT, respectively, for FAD •− -TrPH •+ and PY •− -DMA •+ RP.

FIG. 6 .
FIG. 6. Dynamics of the statistical signal generated by an ensemble of RPs in the sensing volume of the NV center.To reduce computation time, two nuclear spins were included for each radical: N 5 , N 10 for FAD •− , and N 1 , H 1 for TrPH •+ ).MFE against variation in (a and c) strength of magnetic field (b and d) direction of the magnetic field (θ corrected).The solid line indicate the average MFE when the coordinate frames of all the RPs are aligned ( R ζ = I) with the NV frame, while the dashed lines show the average MFE from a randomly oriented ensemble (mean and variance of 50 random realizations of R ζ ).

FIG. 7 .
FIG. 7. (a) Energy level structure in the strong coupling regime.The RP components of an eigenstate in |0 manifold (|ψ i ) are determined by H RP alone whereas in the eigenstate |1 manifold, (|ψ i ) are determined by both H RP and H c together.Here {(∆ i − ∆ i )} depend on the strength of the coupling tensor D D D. The sequence to detect RP dynamics in this regime is shown in (b).(c) The number (each circle represent one peak) of peaks in spin resonance spectrum of NV center resulting from interaction of the NV center with an RP as a function of applied magnetic field for a single RP.Depending on the magnitude of various couplings in H RP (ζ), the |0 and |1 states of the NV center are further split into a maximum of 2 N+2 levels, as shown in Fig. 7 (a).The (i) Since we wish to study RP dynamics exclusively under an applied magnetic field, a pulsed sensing sequence (Fig 7 (b)) is preferred as it ensures that the NV center stays in the |0 state during the time of evolution, implying no backaction of the NV on the RP.In a Ramsay type of sequence where the NV center is prepared in (|0 + |1 )/ √ 2, the RP experiences an additional induced field from the |1 state component [44] which is different for each RP, complicating the probing of MFEs.

Fig 9 ,FIG. 10 .
FIG. 10.The generated signal X Ix and X I z , and singlet yield for increasing values of RP lifetime.Here B 0 = 50µT and A A A 1 1 1 = A A A 2 2 2 and chosen according to the case 'axial 3 ' (see caption of Fig.8).