{'fraction_non_alphanumeric': 0.03882448099218118, 'fraction_numerical': 0.015300620113238069, 'mean_word_length': 4.4587932303164095, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quantum Simulation of the Radical Pair Dynamics of the Avian Compass

The simulation of open quantum dynamics on quantum circuits has attracted wide interests recently with a variety of quantum algorithms developed and demonstrated. Among these, one particular design of a unitary-dilation-based quantum algorithm is capable of simulating general and complex physical systems. In this paper, we apply this quantum algorithm to simulating the dynamics of the radical pair mechanism in the avian compass. This application is demonstrated on the IBM QASM quantum simulator. This work is the first application of any quantum algorithm to simulating the radical pair mechanism in the avian compass, which not only demonstrates the generality of the quantum algorithm, but also opens new opportunities for studying the avian compass with quantum computing devices.


I. INTRODUCTION
An open quantum system is a quantum system that interacts with an external environment or bath. The interaction between the system and the bath is often too complex to be simulated exactly, and thus requires approximations to average out the effects of the bath -this results in the non-unitary dynamics of open quantum systems. Simulating the dynamics of quantum systems has been a main focus of quantum computing research [1][2][3][4][5][6], yet relatively few quantum algorithms have been developed for simulating the dynamics of open quantum systems [7][8][9][10][11][12][13][14]. To this end, we have developed and demonstrated a general quantum algorithm for open quantum dynamics [15][16][17] that is capable of simulating general and complex physical systems. The quantum algorithm leverages the Sz.-Nagy unitary dilation approach to convert non-unitary time evolution operators into corresponding unitary operators, which can then be implemented on a quantum circuit. This quantum algorithm has been applied to a variety of physical systems, including the amplitude damping channel described by the Kraus representation [15], the Jaynes-Cummings model described by the Kraus representation [18], the Fenna-Matthews-Olson (FMO) complex described by the Lindblad master equation [16], and the spin-boson model described by the generalized quantum master equation (GQME) [17].
In this work, we apply the general quantum algorithm to simulating the radical pair mechanism in the avian compass and further demonstrate its generality. The radical pair mechanism (RPM) is a theory proposed to explain the magneto-reception and navigation abilities of certain bird species [19][20][21]. Many animals possess extraordinary abilities to sense the direction by perceiving the geomagnetic field. This is probably the result of natural selection over a very long time of evolution, since In brief, the RPM involves two spatially separated electrons, which are correlated with each other and affected by the external weak magnetic field and internal nuclear spin couplings. The basic scheme of the RPM includes three stages: 1) the photons with certain energies activate a certain type of molecules located in the bird's eyes, enabling an electron transfer reaction and generating a radical pair in the singlet state; 2) the state of the radical pair converts between the singlet state and the triplet state under the influence of the external magnetic field (the geomagnetic field) and the internal magnetic field (the hyperfine coupling effects); 3) the radical pairs in different states will generate different chemical products which can induce a detectable signal for birds to recognize the direction [22]. The RPM is a promising hypothesis that can explain the three unusual properties of the avian compass: 1. the inclination compass: the functional mode of the avian magnetic compass is based on the inclination of the field lines instead of their polarity [23][24][25][26]; 2. the light dependence: light with an energy above a certain threshold is needed for the RPM to work [27][28][29][30][31][32][33]; 3. the narrow range of responsive magnetic field intensities: both higher and lower magnetic fields will disable birds' ability of navigation [34].
To understand the RP mechanism, E. M. Gauger et. al.
proposed a way to model the dynamics of the RPM system with a Lindblad master equation by adding two "shelving states" for the singlet yield and triple yield [35]. In the following, we treat the same Lindblad formulation of the RPM dynamics with our general quantum algorithm for open quantum dynamics and simulate the RPM dynamics on the IBM QASM quantum simulator. To our best knowledge, this is the first ever demonstration of any quantum algorithm applied to simulating the RPM dynamics. This work not only shows the generality of the quantum algorithm, but also opens new potential ways to study the avian compass with quantum computing devices.

A. The General Quantum Algorithm for Open Quantum Dynamics
The general quantum algorithm for open quantum dynamics has different versions that can evolve the Kraus representation [15], the Lindblad master equation [16], and the generalized quantum master equation [17]. In this work we use the version for the Lindblad master equation [16].
We first review how the generalized quantum algorithm converts the non-unitary matrix that encodes the open quantum dynamics into a unitary evaluation based on the Sz.-Nagy unitary dilation procedure [15,16]. We assume the initial density matrix that describes the physical system is composed of a set of unique pure quantum states |φ i that are weighted by their corresponding probabilities p i : we want to simulate the time evolution of ρ(t) given the initial ρ and the Kraus operators M k 's. This task can be achieved by preparing each input state |φ i in a vector form v i in a given basis and then building a quantum circuit that generates the quantum state: The U M k is generated via the 1-dilation of M k : where . After obtaining each |φ ik (t) , we can calculate the population of each basis state in the current basis from the diagonal vector: where diag(|φ ik (t) φ ik (t)|) can be efficiently obtained by applying projection measurements on the first half sub- For the dynamics of an open quantum system, the time evolution of the density matrix can be represented as: where ρ(s) is the density matrix at time step s, and δs is considered as the discrete time step, during which the Kraus operators M sk are assumed to be constant. The formula eq.(4) can be used iteratively until reaching the time of interest. Explicitly, the dynamics of the density matrix is described as: Here without losing any generality, the Kraus operators M sk are indexed by the time step s, which allows each time step to have a different set of Kraus operators. However, as will be discussed in Section II C, the RPM dynamical model used in this work is a Markovian process described by the Lindblad master equation, therefore all the time steps have the same set of Kraus operators M k where the time step index s has been removed.

B. The Radical Pair Mechanism Theory and Dynamics
The basic scheme of the RPM is shown in Fig. 1. To simplify the fields, we assume that only the electron near the donor interacts with the nucleus, and the electron away from the donor is not affected by the anisotropic hyperfine coupling [35]. Therefore, the Hamiltonian of the system is where A is the anisotropic hyperfine tensor coupling the nucleus and one of two spatially separated electrons, and B 0 (cos ϕ sin θ, sin ϕ sin θ, cos θ), and B 0 is the magnitude of the geomagnetic field; ϕ is the angle between the xaxis of the radical pair and the external magnetic field; θ is the angle between the z-axis of the radical pair and the external magnetic field; γ = 1 2 µ 0 g, and µ 0 is the Bohr magneton and g = 2 is the electron-spin g-factor.
To model the dynamics of the system with a quantum master equation formulation, two "shelving states" were added to the 8-dimensional Hilbert space of the three spins (two electron spins and one nuclear spin) [35]. We employ operators as shown in Eq.(9) to represent the spin-selective relaxation into the singlet shelf |S from the electron singlet state, or the triplet shelf |T from the electron triplet state. The final populations of |S and |T give the singlet and triplet yields.
With the usual definition of singlet |s and triplet states |t i in the electronic subspace, while |↑ and |↓ describing the states of the nuclear spin, we define the following decay operators: This gives a standard Lindblad master equation: where k d is the decay rate of the singlet and triplet states. Note the decay rate k d is independent of the radical pair states, so we have assigned the same decay rate k d to all eight projectors.

C. Calculation and Simulation
Now to apply the general quantum algorithm to the RPM dynamics, we first consider the non-unitary part on the right side of Eq. (10) which can be rewritten as: Given a very small δt, Eq. (11) becomes: + O(δt 2 ). Now assuming M 0 = I − 1 2 k d δt 8 k=1 P † k P k and M k = √ k d δtP k for k > 0, and ignoring the second order of δt as δt → 0, Eq. ((12)) can be rewritten as: Eq.(13) is in the same form of Eq.(4), thus knowing the initial state ρ(0) we can evolve the density matrix to a certain time with the iterative procedure described in Eq.(5) to Eq. (7).

Parameter Details Symbol
Description Values Ax Anisotropic hyperfine tensor 1 × 10 −4 T B0 Magnitude of the geomagnetic field 5 × 10 −5 T γ Half of the product of the Bohr magneton and electron-spin g-factor 9.27 × 10 −24 J/T Reduced Planck constant 1.05457 × 10 −32 J · s ϕ Angle between x-axis of the radical pair and the magnetic field 0 k d Decay Rate of the singlet and triplet states 1 × 10 4 s −1 In the meantime, When δt → 0, according to the above equation, we have: With Eq. (15), we can formulate the Kraus operators to satisfy the condition k M † k M k = I. Thus, we have defined all the 9 Kraus operators M k required to describe the RPM dynamics.
There is one additional term containing the Hamiltonian, − i [H, ρ], in Eq. (10) as compared with Eq. (11). This "oscillating part" of the dynamics is unitary and thus can be easily realized by multiplying each Kraus operator by a unitary matrix obtained through the diagonalization of the Hamiltonian [16].
With the parameters in Table I, we simulated the RPM dynamics by iteratively applying Eq. (13) on the IBM QASM quantum simulator and then used the output results to calculate the singlet and triplet yields. Also, we assume the initial state of the two electron spins is 1 √ 2 (|↑↓ − |↓↑ ), and the initial state of the nuclear spin is 1 √ 2 (|↑ − |↓ ). In our simulation, the time interval δt is set as 0.5/k d = 5 × 10 −5 s. We then apply the procedure shown in Eq.(5), Eq.(6) and Eq.(7). As mentioned above, the populations of | S and | T is the singlet and triplet yields respectively, after the system reaches the steady state. The populations are calculated by the procedure explained in Eq. (3), where the diagonal elements of the density matrix are obtained by projection measurements into the computational subspace.
The results are compared with those obtained from classical methods in Fig. 2, where the quantum algorithm results are highly consistent with the classical method results. Fig. 3 shows the dynamical evolution of the singlet and triplet yields (when θ = π 2 ) as simulated by the general quantum algorithm on the IBM quantum simulator. After 2 × 10 −4 s, the yields almost reach steadystate, which is consistent with the chosen decay rate of k = 1 × 10 4 s −1 .

D. Complexity analysis
One factor that contributes to the complexity of the quantum algorithm is the system's size. For a density matrix of the size n × n, the cost to realize the unitary dilation of a most general n × n Kraus operator M k is O(n 2 ) [15,16]. However, in our calculation of the dynamics of the RPM, the Kraus operators each represents a single elementary physical process and thus the M k matrices are often sparse with few non-zero elements. This means the practical complexity scaling of implementing each M k matrix on a quantum circuit can be greatly reduced to O(log 2 n). Taking into account the total K number of M k matrices to be simulated on the quantum circuit, the total complexity scaling is O(K log 2 n) for our given system. It is worth noting that the K is determined on a case-by-case basis by the dynamical model and different M k matrices can be evolved in parallel, therefore the scaling in K is a "soft" scaling that does not contribute to either the depth or the width of each individual quantum circuit [15,16]. Another contributing factor to the quantum algorithm's complexity is the number of time steps. In the most general case, as can be seen from Eq.(5) to Eq. (7), taking s steps requires K s matrices to be evolved, which is an exponential scaling in the number of time steps. However, fortunately the actual number of matrix terms to be simulated can be greatly reduced once again due to the sparsity of the M k matrices. As mentioned above, the Kraus operators represent elementary physical processes and thus the M k matrices are often sparse with very few non-zero elements: this means that most matrix product terms in e.g. Eq.(7) are zero matrices or matrices with negligible norms. The actual number of matrix products we need to evaluate is determined on a case-by-case basis. In the current simulation, in theory the total number of terms in n-th iteration will be 9 n with nine Kraus operators {M k |k = 0, 1, 2, ..., 8}. However, since the product of each pair of the the decay operators {P k |k = 1, 2, 3, ..., 8} is 0, the product of each pair of the Kraus operator {M k |k = 1, 2, 3, ..., 8} is 0. Therefore, there will be only 8 more terms when adding one more iteration. Thus, there will be 8 × n + 1 terms in n-th iteration, reducing the terms significantly. More simulation details are in the Supplementary Information.

III. CONCLUSION
Based on the calculations and results, we can conclude that the general quantum algorithm based on the Sz.-Nagy dilation can accurately simulate the RPM dynamics described by the Lindblad master equation. As discussed above, the RPM is an important theory that can explain the magneto-reception process of the avian compass. The RPM also acts as an ideal model to help explaining the isotope effects in xenon anaesthesia and lithium treatment of hyperactivity, magnetic field effects on the circadian clock, as well as hypomagnetic field effects on neurogenesis and microtubule assembly [37,38].
Our simulation of the RPM model with the quantum algorithm not only demonstrates the generality of the algorithm but also helps bridging the gap between applying tools of quantum-information science to the investigation of new areas of quantum biology.