Simulation of non-radiative energy transfer in photosynthetic systems using a quantum computer

Photosynthesis is an important and complex physical process in nature, whose comprehensive understanding would have many relevant industrial applications, for instance in the field of energy production. In this paper we propose a quantum algorithm for the simulation of the excitonic transport of energy, occurring in the first stage of the process of photosynthesis. The algorithm takes in account the quantum and environmental effects (pure-dephasing), influencing the quantum transport. We performed quantum simulations of such phenomena, for a proof of concept scenario, in an actual quantum computer the IBM Q, of 5 qubits. We validate the results with the Haken-Str\"obl model and discuss the influence of environmental parameters on the efficiency of the energy transport.


Introduction
Photosynthesis is a vital and pervasive complex physical process in nature, where the radiation of the Sun is captured by certain living beings, such as plants and bacteria, and transformed into the necessary carbohydrates needed for their survival [29,35]. From the physics and chemistry perspective, it is a complex process occurring through several stages with several kinds of physical phenomena involved, namely, the light absorption, energy transport, charge separation, photophosphorylation and carbon dioxide fixation [17]. The understanding of such phenomena has greatly progressed in the the past 40 years with the physical characterization of the structure of many photosynthetic complexes [7,12,48]. The comprehension of such processes would allow for many potential huge-impact industrial breakthroughs in the field of energy, from the great efficiency improvement in energy capture of solar panels [32] to the construction of artificial light-harvesting devices and solar fuels [20,21,45,57].
The photosynthesis begins by the absorption of a photon. It occurs via excitation of a pigment molecule, which acts as a light-harvesting antenna connected to the rest of the photosynthetic apparatus by protein molecules. Photosynthetic pigment-protein complexes transfer the absorbed sunlight energy, in the form of molecular electronic excitation, to the reaction center, where charge separation initiates a series of biochemical processes [35]. This work is focused on the first stage of photosynthesis, more precisely, on the transport of the absorbed radiation energy from the antenna to the reaction centre, which proceeds in the form of the so-called Excitonic Energy Transfer (EET), as schematically shown in Fig.1. Figure 1: Schematics of the energy transfer process from light-harvesting antenna (the donor) through a chain of acceptor molecules to the reaction center. The excited states of the participating molecules, denoted m , are broadened and it allows for resonance energy transfer via irreversible Förster-type resonant process of exciton transfer from donor to acceptor even if m = m+1 , which is denoted by the thick arrow labelled FRET. However, if the coupling between the donor and the acceptor molecules is strong enough, the process becomes reversible and the exciton can go to and through many times before it is transferred; this situation is labeled by "reversible EET" and it does not require matching of the energy levels m and m+1 .
This transport is known to be very efficient in photosynthesis, as is the whole process, with the overall quantum efficiency of initiation of charge separation per absorbed photon up to 95% [35]. The absorbed photon creates an exciton on the antenna molecule, which can eventually transfer it to other molecules. In this context, it is called donor, while the others are called acceptors and the EET process can be described by the following reaction equation: The physics of the mechanisms behind Eq. (1) will be discussed in the following section. Here we just notice that EET is a complex process that can be irreversible (i.e. unidirectional) or reversible, i.e. coherent over some period of time as evidenced by experimental observations of long-lived oscillatory features in the dynamical response of several photosynthetic systems [14,23,30]. Moreover, it is strongly influenced by the environment. The donoracceptor pair is not isolated from the rest of the world and is an example of so called open systems [42]. It must be treated as a subsystem of a larger system including a thermal bath. The properties of the latter are crucial because it introduces relaxation and dephasing into the system directly involved in the EET and, therefore, influences the efficiency of the energy transport.
Open quantum systems cannot be described by a wave function because one does not have enough information to specify it, only a (less detailed) description in terms of a density matrix is possible, which represents a statistical mixture of states, or a mixed state (see Supplementary information A.1). The dynamics of such a system can be determined by solving an equation of motion for the density matrix. Such equations of motion are called quantum master equations. Finding exact solutions to the master equations is extremely difficult but there is a wide range of theoretical approaches and techniques available to make their mathematical simplification and numerical simulations. These approaches can be divided into several groups according to the regime under study, characterized by the coupling strength between the bath and the system, and the existence of memory effects in the bath (i.e. whether the system can be considered as Markovian or not). Broadly speaking, for the weak-coupling Markovian regimes, perturbative approaches are applicable, such as the Bloch-Redfield and Lindblad master equations [26,35], which can be extended to medium coupling strengths and non-Markovian regimes by including higher-order system-bath interaction terms [25,44]. For the latter regime, there are also non-perturbative techniques based on the use of path integrals to dissipative systems [16], which can be used to create sets of solvable systems of hierarchical equations, the so-called hierarchical equations of motion (HEOM) [24,50,51]. Open quantum systems that do not have the Markovian property, for example, because of a too small size of the bath effectively coupled to the system, which keeps memory of the past, are much more difficult for theoretical description because the dynamics equation is non-local in time.
However, even within the Markov approximation, the calculations quickly become computationally intractable for realistic photosynthetic systems and environment models. The computational cost of simulating a photosynthetic complex consisting of N molecules with a theoretical tool such as the HEOM grows exponentially with N . A possible computational solution has been arising to bypass this type of problems is the use of quantum simulation, where it is expected to obtain large performance increases in terms of space, as the number of qubits' growth is just polynomial, and in terms of time, where an exponential gain is expected.
The use of quantum mechanics to make calculations about quantum mechanics, promising great computational advantages was firstly proposed by R. Feynman [15,33]. The field of quantum simulation is under a fast paced and intense development, finding already application across all fields of physics and using many different physical implementations [19]. Closer to the present work, there are works on quantum transport [34,47] and on the quantum simulation of dissipative systems [13]. Particularly on the quantum simulation of photosynthesis, we would like to highlight Refs. [37,41,52], using superconducting qubits, and [54] employing a Nuclear Magnetic Resonance (NMR) simulator [58]. The latter is of particular relevance to the work carried out here as it was dedicated to the quantum simulation of the energy transport with environmental actions, where the environment effect is simulated naturally by an appropriate filtering of environmental noise [49], within the NMR system. In this case, the implementation is specific for the EET (i.e. non-universal), and the model Hamiltonian was extracted from spectroscopic data for a photosynthetic system [2]. We are simulating the same Hamiltonian as in Ref. [54] and starting from the same assumptions, however, the simulation algorithm is completely different since we conceived a digital quantum simulation designed to run in a universal quantum computer, the commercially available IBM Q of 5 qubits [11]. Our implementation contains a quantum part, aimed at simulating the unitary part of the system's evolution, and a classical part that simulates the stochastic interaction with the environment, the latter only being able to mimic pure dephasing environmental effects.
The physics of the energy transport in photosynthesis Förster

and Redfield approaches
The molecules of the light-harvesting complexes usually are not electronically coupled to each other and charge transfer via electron tunneling is improbable. Hence, energy transfer can occur between them through electromagnetic interaction, without net charge transport because the whole (neutral) exciton is transferred. Such processes are known to take place between molecules [3] or artificial nanostructures such as quantum dots [46] if appropriate conditions are met, which were first formulated by T. Förster [18]: (i) The distance between the donor and acceptor molecules must be sufficiently small because the transfer probability decreases quickly with the distance between them (R), usually as R −6 ; (ii) There must be a resonance between the excited states of the donor and acceptor molecules; 1 (iii) An increase of the refractive index of the surrounding medium decreases the transfer rate. Förster's approach is based on the second-order perturbation theory (the so called "Fermi's Golden Rule"), where the perturbation operator is the electromagnetic interaction between two transient dipoles corresponding to allowed optical transitions in the donor and acceptor molecule, respectively. It originated the term "Förster resonance energy transfer" (FRET), which applies to an irreversible hopping of an exciton from the donor to the acceptor. The FRET rate (transition probability per unit time) can be expressed by the following relation [35]: where J is the coupling constant, ω is the angular frequency of the electromagnetic field and L D (ω) and I A (ω) denote dimensionless lineshape functions of the donor and acceptor molecules, directly related to the energy spectrum of each molecule. The integral is called the spectral overlap between the molecules. The coupling constant, in the dipole-dipole approximation, is given by [3]: where η is the refractive index of the medium, d D (d A ) is the transient dipole moment of the donor (acceptor) molecule, n = R/R, R is the radius vector between the two molecules and the angular brackets stand for angular average over different orientations of the dipoles. Even though Eq. (2) (and the approach itself) is too simplistic to describe all possible situations in EET, defining this characteristic transfer rate allows for the formulation of the following conditions for FRET to occur: • If the difference between the energy of the excited state energy of the donor ( 0 ) and acceptor ( 1 ) molecules is small, | 0 − 1 | J and they are in resonance, the energy transfer between the molecules can occur with a high probability; , the exciton is trapped in the donor molecule because it has a very low probability of being transferred; in this case it either stays in the molecule and later the donor molecule will decay to the ground state, dissipating the energy, or transfer the energy to a different acceptor nearby.
As pointed out above, the initial idea of Förster was that an exciton is irreversibly transferred from a donor to an acceptor. More recently, it has been shown experimentally that quantum coherent transport, where energy is transported in the form of wave-packets, has a significant role in many important physical effects, including the photosynthesis [9,10,39]. The Förster theory does not apply in this regime, as it simply ignores coherence. Later, in 1957, A. Redfield [43] proposed a transport theory, which applies to the opposite regime of strong coupling between the donor and acceptor [9,31,35] (although originally it appeared in the context of NMR spectroscopy). Within this concept, the exciton forms a coherent state based on the whole donor-acceptor pair and oscillates between the two molecules. This system can be described by the following Hamiltonian: where |m denotes the exciton on the molecule m. The eigenstates of (4) are linear combinations of |0 and |1 . Coherent dynamics corresponds to the presence of non-zero off-diagonal elements in the density matrix describing the evolution of the quantum system. Their oscillation (or quantum beating) is indicative of coherence [35]. For the system with Hamiltonian (4), the description in terms of state vectors is perfectly possible but it will not be the case if interactions with environment are taken into account. Therefore, we may introduce the density matrix description at this point. The evolution of the off-diagonal elements of the system's density matrix, written in the energy basis (where the Hamiltonian is diagonal) and denoted ρ ij , is given by: (see Supplementary Information A.2 for the derivation). These states are perturbed by interactions with the environment (the bath), which destroys their coherence. Mathematically, it is expressed in the form of a master equation, which is known as the Bloch-Redfield equation; its general form can be found e.g. in Ref. [26].
This consideration is extendable to a chain of molecules and can be seen as (partially) coherent transport [9]. The presence of the latter, observable through coherent oscillations of the energy levels of molecules across different sites (the quantum beating), was first conjectured in the 30's [40] and theoretically predicted in more recent works [28,31]. It became possible to observe them more recently, thanks to the advances of optical spectroscopy techniques [14,23,30,53], and it was achieved even at room temperature [39]. In these experiments, it was possible to confirm the substantial impact of such coherent effects on the excitation energy transfer in photosynthetic systems [35]. Moreover, the importance of environmental noise in the quantum transport involving coherence was also discussed more recently [6,42] and it is not fully understood yet.

Decoherence
Processes caused by the molecules' environment may destroy coherence and thus influence this type of energy transport [35,36,42], moreover, they can foster it. Indeed, completely coherent oscillations (called Rabi floppings in atomic physics) between different molecular sites do not correspond to an energy flux. Breaking the oscillatory evolution at some moment may help transferring the exciton along the molecular chain.
If interactions exist between a system and its environment, they affect the (pure) states of the system, introducing "errors" and making these states mixed. It means the so-called phenomenon of decoherence, which, by the way, has been the main obstacle to the success of quantum computation. Decoherence processes can be divided into three categories: (i) amplitude damping, (ii) dephasing, and (iii) depolarization, which are briefly described below. [5] Amplitude damping. Environment interactions with the system may cause a loss of the amplitude of one or more system's states. The spontaneous emission of a photon from the system (i.e. from one of the molecules) to the environment is an example of this kind of process, so that the system returns to its ground state (without exciton) [38]. For a two level system (e.g. a qubit), this type of decoherence contracts the Bloch sphere along the z axis (see Supplementary Information A.1).
Phase damping or dephasing. Such interactions conserve the energy of the system, contrary to the amplitude damping. A phase damping channel removes the superposition of the system state, i.e. the off-diagonal terms of the system's density matrix decay over time down to zero. It is a process of removing the coherence of the system, causing a classical probability distribution of states and, therefore, imposing some classical behaviour in a quantum system. A simple way to look at this type of decoherence is also to think of the system interacting with the environment where the relative phases of the system's states become randomized by the environment. This randomness comes from a distribution of energy eigenvalues of the environment. As a result, the evolution of the quantum system's Rabi cycle ceases but the time-average populations of the states may not change and this is the case of the pure dephasing. For a two-level system with a pure dephasing interaction, the Bloch sphere contracts in the x − y plane.
Depolarization. This type of decoherence changes system's state, which initially is pure, to a mixed state, with a probability P of another pure state and the probability (1 − P ) of the initial state of the system. It is equivalent to saying that, for a single qubit, an initial pure state represented on the Bloch sphere has suffered a contraction over all dimensions of the sphere (with the contraction degree that depends on the probability P ). It can be thought of as a combination of the other two types of decoherence.
The amplitude damping is certainly detrimental for EET since the energy is simply dissipated into the environment. The action of dephasing processes progressively eliminates the coherence (off-diagonal) elements in the system's density matrix, causing the oscillation amplitude to decay (beating supression). It eventually turns the diagonal matrix elements (populations) into (non-correlated) classical probabilities, a process known as thermal relaxation, for which the existence of coherence in a system is time limited. On the other hand, it has also been shown that dephasing processes can have a positive role in the coherent transport of energy [35]. First, it yields random fluctuations in the energy spectrum of each molecule, which can bridge the energy gap between the molecules, momentarily turning a non resonant system into a resonant one. Secondly, dephasing can also help avoiding the existence of the so called coherence traps in a molecular chain, a kind of deadlocks in energy transport where the exciton can be confined [35]. Thus, the result of action of a decoherence source on an EET system is not obvious a priori. Below we shall consider a simple model of pure dephasing consisting in a telegraph-type classical noise affecting the donor-acceptor pair.

Materials and Methods
We aim at exploring the energy transport underlying the photosynthesis, throughout time, under two regimes: (i) in an isolated system and (ii) under an action of the environment causing decoherence. In the "no decoherence" case (i), one can study the evolution of system's state vector, which obeys the equation for a time-independent Hamiltonian. Here t f and t i are the upper and lower limits of the time interval to study. In order to be able to do the calculation of the system's evolution on a quantum computer, it is necessary to provide a suitable qubit encoding for the possible states and an approximation for the Hamiltonian evolution operator,Û , in terms of quantum gates and circuits (computational Hamiltonian). The time, a continuous entity in equation (6), has to be discretized onto a set of intervals, ∆t, where the Hamiltonian of interest can be approximated as constant.
The actual computational process is given by the repeated application of the evolution operator on the prepared state |Ψ i , for s times, of the computational Hamiltonian, such that t i + s · ∆t = t f . The process is finished by the observation of the desired properties, i.e. a set of measurements, in the appropriate basis, on the end state.
Concerning the particular qubit encoding chosen, a chain of N = 2 q molecules is encoded by a set of q qubits, where |m corresponds to the excitation (exciton) on the m-th molecule, e.g. for a two-molecule chain, state |0 represents the exciton on the first molecule and |1 on the second one, and a possible successful transport of energy would correspond to the transition of the state |0 to the state |1 . We denote this as the site basis. The computational Hamiltonians under this encoding for the cases under study are discussed in the following sections. From now on, we shall set = 1. Also, it is convenient to measure the energies/frequencies in cm −1 , as it is common in spectroscopy.

No-decoherence Hamiltonian
Considering a small chain of N molecules, the system's Hamiltonian in the site basis reads as follows, where m is the first excited state energy of the molecule m and J nm is the electronic coupling between the molecules n and m. The Hamiltonian (7) for just two molecules (1 qubit), identical to Eq.(4), in the 2 × 2 matrix form, reads:Ĥ Its evolution operator is given by Although the Hamiltonian (44) possesses non-diagonal elements, finding a good approximation in terms of quantum circuits is relatively straightforward. A possible strategy for this is by finding a diagonalizing transformation, T , of the Hamiltonian, such that, whereĤ S−diag is the diagonal Hamiltonian. Therefore, the evolution operator can be rewritten as follows: The problem now reduces to the approximation of the T operator (and its adjoint) and the HamiltonianĤ S−diag , which can all be efficiently approximated in quantum circuits. The latter operator is diagonal in the site basis, thus the unitary evolution operator can be expressed aŝ The T and T † matrices can be implemented by simple rotations, R y (θ) and R y (−θ), for a two-molecule system. However, for a higher number of molecules, a rotational decomposition algorithm together with the Gray code [38], which decomposes a matrix in the multiplication of a single qubit and CNOT gates, has to be used. Using this particular algorithm the gate complexity for N molecules is O(N 2 log 2 [N ]) [38]. On the other hand, the diagonalized evolution operator,Û translates into trivial phase rotations over each of the energy eigenstates |E i of the system with the respective energy eigenvalues E i . This operator can be constructed as a sequence of CR Z (φ i ) gates applied to an ancilla qubit (initialized at |1 ), where the angle is given by φ i = −2E i t, i = 1, 2. The X gates are used to "select" the eigenvector to which the controlled rotation is to be applied. The circuit implementation of the operator defined in (13) is illustrated in Figure 2. The gate complexity of this operator, in terms of single qubit and CNOT gates for N molecules, is O(N log[N ]). Figure 2: Implementation of the system's evolution operator. |q system is the state vector of the system's qubit in the energy eigenbasis.
For the whole circuit, resulting from the sequencing of T †Ĥ S−diag T , the number of qubits required to simulate a molecular chain of N elements is 2 log 2 N and the gate count scales with O(N 2 log 2 2 N ) single qubit and CNOT gates. The transformations T and T † , in the general case, possess a high circuit depth, which makes the system hard to simulate accurately, with low error rate, in the current available quantum computers.

Introducing decoherence into the system
We shall implement artificial decoherence as pure-dephasing by adding Markovian fluctuations to the Hamiltonian. This approach is considered a good approximation in the high-temperature regime for the bath [5,31,42]. The actual algorithm to be used is the one of [55], which is used to simulate open quantum systems, with pure dephasing, modeling the action of the decoherence as classical random fluctuations (a telegraph-type classical noise affecting the system). The actual Hamiltonian for this system reads aŝ and it consists of the system Hamiltonian,Ĥ S , of the previous section and the perturbation of a bi-stable fluctuator environment,Ĥ F . The latter simply shifts the energy by a constant value for each molecule, ±g m /2, as illustrated in Fig. 3. Explicitly,Ĥ whereÂ m |m m| is the projection operator and considering one fluctuator interacting with each molecule m, The function ξ m (t) switches the fluctuator between the positive and negative values (appearing randomly) at a given fixed rate γ and g m is the fluctuation strength (or the coupling strength to a molecule m). Physically, the action of the fluctuations is typically stronger for the excited states [1,31] and g can be larger than the donor-acceptor coupling J. Figure 3: Uncorrelated random fluctuations applied to donor and aceptor's excited state energies, 0 and 1 . Each molecule is affected by one fluctuator, which generates a telegraph-type classical noise. The fluctuators switch randomly between the positive and negative value at a given fixed rate, so that the periods of time when the molecule energy is constant, m + g m /2 or m − g m /2, are random. J is the coupling strength between the molecules that can be seen as the rate of hoppings between these fluctuating energy levels.
The implementation of such random bi-valued function ξ m (t), can be done in a straightforward way by a classical pseudo-random numbers generator with a probability of 50% of the values −1/2 and 1/2. For circuit generation purposes, the values resulting from the random sampling have to be provided in advance of the quantum simulation. The fluctuator interaction Hamiltonian and the system Hamiltonian do not commute, so, in order to generate an appropriate quantum circuit, one needs to apply an approximation technique such as the Trotter product formula [56]. Under this approximation, the unitary evolution operator of the Hamiltonian, for a time t = N i ∆t, where N i is the number of iterations and ∆t is the iteration time-step, becomes where E m denote the eigenvalues of the system Hamiltonian.
Note that the projection operatorÂ m is not present in the evolution operator (17) because the latter is used in its eigenbasis, i.e. the site basis. The fluctuator interaction evolution operator e ±i gm 2 ∆t is a selective rotational gate over a molecule m (|m ), which can be implemented by a set of X gates and a controlled gate CR Z (φ m ) with angle φ m = ±g m ∆t, applied over an ancilla qubit initialized at |1 . The whole circuit is presented in Fig. 4 for one iteration. The fluctuator waiting time (interval of time between switches), i.e. 1 γ , can only be equal or higher than the iteration time-step, ∆t. The switching in the fluctuator-molecule coupling strength is performed at every 1 γ∆t iterations, where a∆t = 1 γ , a ∈ N.
Usually in the study of open quantum systems with a dilated system's Hilbert space (as is the case here), different measurement techniques are required [38,55], however in this case, the open system is simulated in a closed form so, similarly to the no decoherence case, the measurement over the site basis suffices. The full algorithm (random values generator plus the actual simulation) must be performed several times, so that the results of all runs are averaged.
Let us consider the simulation for a time t, using an iteration time step ∆t, and assuming that the environment can have more than one fluctuator interacting with each molecule as well as the chain can have more than just two elements. Then the fluctuator interaction evolution operator requires the following gate resource complexity for a single run: O( t ∆t [N (log 2 N + F )]) single qubit and CNOT gates, where N is the number of molecules and F is the number of fluctuators interacting with each one. Figure 4: Implementation of one iteration of the system with decoherence algorithm. Here |q system represents the system's qubit state vector in the site basis.
In the implementation of the system with decoherence, the algorithm gate resources complexity is O( t ∆t [N 2 × log 2 2 N + N F ]) for a single run. This simulation, yet again, possess a very high circuit depth which makes its application unfeasible in quantum computers. The number of necessary qubits is the same as in the no decoherence simulation (2 log 2 N ).
It also requires O(N R F j=0 tγ j ) random numbers to be classically generated, where R is the number of runs of the algorithm and γ j is the switching rate of the fluctuator j interacting with the molecule. The number of required simulation runs to average the results and obtain an error > 0, is predicted to scale as O [F t ∆t ] 2 / 2 . This complexity is calculated based on the possible non-degenerate energy state outcomes of the entire chain in the simulation for a time t. These outcomes are caused by the bi-stable random fluctuations, therefore, the possible non-degenerate energy state outcomes for each molecule obey a discrete Gaussian probability distribution.

Results
We conducted simulation experiments for the quantum transport in a molecular chain using the algorithm described in the previous section. We executed the simulation for the coherent system on a real quantum computer, the IBM Q of 5 qubits, while the pure dephasing scenario was simulated on the QASM quantum simulator, both in the nearresonant and non-resonant regimes. For the validation purposes, we compared the results for the coherent system with the theoretical predictions obtained by solving the Schrödinger equation (see Supplementary Information A.2).
As for the decoherent regime, we used a classical computation of the stochastic Haken-Ströbl model [22,42]. The simulations and circuits involved, encoded in the Qiskit platform [11], can be performed in the following url: https://github.com/jakumin/Photosynthesis-quantum-simulation.

Coherent regime
The scenario for this regime was simulated with a simple chain of two molecules. As discussed in section Materials and Methods and using the parameters as proposed in [54], we define the system's Hamiltonian as follows: (Near-resonant regime) The results for both regimes were obtained using an actual quantum device (the IBMQ london of 5 qubits) and can be seen in Figs. 5 and 6, respectively. Due to the stochastic nature of quantum computers, the experiments were conducted with 2048 shots for each time value. The specific optimized quantum circuits used in this experiment are presented in Supplementary information A.3. In the following results, the probability of the donor and acceptor molecules being excited is denoted by P (0) = 0| ρ S (t) |0 and P (1) = 1| ρ S (t) |1 , respectively.  Taking the fluctuator's switching rate to be γ = 0 or the fluctuator-molecule coupling strength to be g = 0, one has the coherent regime. These simulations show the limiting case of the Redfield regime, i.e. the very weak system-environment coupling, g J. The quantum beatings, observed in the simulation results, can be thought of as a reversible transfer of energy between the molecules, where the excitation goes back and forth across the molecules [8].
In the performed simulations, the near-resonant and non-resonant regimes have a maximum probability of ∼ 90% and ∼ 20%, respectively, of the energy being transferred to the acceptor molecule. Using the quantum Liouville equation [35] (see Supplementary information A.2), the period of the quantum beating is T near−res ≈ 123 f s for the near-resonant regime and T non−res ≈ 51 f s for the non-resonant regime. These periods are in the femtosecond timescale of the experimentally observable quantum beatings [10,14,39]. The simulation results show a similar behaviour as those predicted by the Schrödinger and quantum Liouville equations, where the off curve points are predominantly originated by errors in the quantum hardware.

Decoherent regime
The scenario for the regime with decoherence introduced is, in some respect, similar to the one presented for the coherent regime for a chain of two molecules. No further changes are made to the Hamiltonian discussed in the section Introduction of decoherence in the system. The quantum simulation results are compared with a theoretical evolution based on the stochastic Haken-Ströbl model, in the form of the Lindbland master equation [22,42]. The Lindbland equations were solved in a classical computer using Qutip [27], a quantum open systems software framework. The set of Lindbland equations, correspondent to the model in this setting, had one free parameter regarding the environment, the dephasing rate, γ deph . The Lindbland equation in the Haken-Ströbl model reads: where L m = |m m| are the Lindbland operators, responsible for the system-environment interaction. The system Hamiltonian, H S , is given by the matrix (18) for the near-resonant system and the matrix (19) for the non-resonant system.
For each quantum simulation performed, a fitting process has been employed by adjusting the dephasing rate of the Haken-Ströbl model, so that the system's evolution in both classical and quantum algorithms have similar behaviours. This enables one to perform a direct comparison between both theories and to find the actual dephasing rate of the modeled environment over the various regimes considered in this work.
The environment contains only one fluctuator interacting with each molecule with switching rate γ = 125 THz. As mentioned above, the dephasing rate, γ deph , for the Lindbland equation is adjusted to the behaviour of the system under the action of a fluctuation strength g.   It is seen in Figures 7 and 8 that oscillation amplitudes decay over time, as expected, due to the loss of relative phase coherence between the excited states of the two molecules, evidenced by the disappearance of the quantum beatings. This is associated with the irreversible evolution when the system loses its capacity of performing coherent transport. Additionally, it is clear that the system is led to a classical distribution of the populations in the site eigenbasis.
In the regime under the study, where the environment is assumed to be at thermal equilibrium, the final probability distribution is calculated in the limit of the classical Boltzmann distribution m| ρ S (t → ∞) |m = const × e − m k B T . Here k B is the Boltzmann constant, T is the temperature of the bath and const is a normalization constant [5]. Taking the limit of very high temperatures, the population terms approach the Boltzmann distribution 0| ρ S (t → ∞) |0 ≈ 1| ρ S (t → ∞) |1 ≈ 1 2 , which is compatible with the results obtained. The relaxation can not be fully observed in Figs. 7a, 8a and 8b because a very large number of iterations would be required for this.
The switching rate must be high enough to observe the dephasing effects. Here we used a value ≈ 33 times larger than the transfer rate, J (that is, the fluctuator waiting time must be shorter than J −1 ). As observed in the simulations, it is a suitable value for observing the relevant effects of random fluctuations in the system. At very low rates, it leads the system's evolution to a behaviour similar to the previously observed in the no-decoherence regime, Figs. 5 and 6.
The time that coherence lasts in the system is essentially defined by the fluctuation strength, g: in Figs. 7a, 7b, 8a and 8b (lower g) the coherence is maintained for some time, while in Figures 7c, 7d, 8c and 8d (higher g) it is quickly suppressed. In the latter regime, an approximated diffusive motion drives the system's evolution, where quantum beating is practically absent. The time that the quantum beating lasts in these simulations (until it reaches an approximate non-oscillating behaviour), is about 350 f s in Figure 7b (near-resonant system) and 200 f s in Figure 8b (non-resonant system), with a fluctuation strength g = 300 cm −1 . At a longer time, it has been experimentally observed to persist (t > 660 f s [39]), a timescale which could be modeled in the present simulation by changing the environment parameters, i.e. lowering the fluctuation strength g as can be observed in Figures 7a and 8a.

Discussion and Conclusions
Two main conclusions can be drawn from the presented results: • There is a very good agreement between the solution of the Schrödinger equation and the coherent quantum algorithm results in the reproduction of the purely oscillatory evolution of the isolated quantum system.
• There also is a good agreement between the results obtained by the Haken-Ströbl model and the quantum algorithm. The increase of the dephasing rate imply an increase in the fluctuation strengths thus, a faster suppression of the quantum beatings can be observed, as predicted theoretically [42].
Therefore, the correctness of the results obtained in the quantum simulations is verified.
The results obtained in the present work were not directly compared with Ref. [54], due to the different timescales used. The major difference lies in the physical implementation, NMR vs universal quantum computer, where there might be an advantage for the former from the point view of the scalability and reliability, at the current state of quantum technology. However, there is a clear advantage of the quantum computer, from the point of easiness of implementation, as it is also possible to implement circuits of arbitrary precision, harder to do with the NMR simulator, which is dependent on a Hamiltonian mapping process. The computational advantage verified for the NMR simulator still holds after the present work, as the number of executions for the algorithm of this work is polynomial on the precision required, although the circuit generation maybe problematic, as a matrix diagonalization operation is necessary (complexity estimated in O(N 3 )).
To conclude, we proposed a quantum algorithm to simulate the energy transfer phenomenon present in general photosynthesis, under the presence of quantum coherence between the molecules and the decoherence effects caused by environmental interference. Using this algorithm we also performed simulations in the commercially available quantum computer of IBM, the IBM Q of 5 qubits for the coherent scenario and in the quantum simulator (QASM) for the decoherent scenario. For validation purposes, we also computed the evolution of analogous systems using well-established (classical) methods in literature, obtaining quite similar results between the methods. The results obtained were also in agreement with the predictions that can be found in literature, for the role of the quantum coherent and dephasing effects in the energy transport of photosynthesis: for the high temperature environment here defined, it was clear that dephasing, modelled as energy fluctuations in the site energies, limited the time quantum coherence lasts in light-harvesting antenna. Moreover, it was also verified that the fluctuation strength and the switching rate of the Markovian fluctuator environment are directly related with the energy transfer efficiency, allowing the simulation of different transport regimes by setting them appropriately.
Similar to Ref. [54], this setting revealed itself as an interesting platform for the study the quantum and environmental effects in a small photosynthetic system, and therefore we consider, that the use of quantum simulations may be a feasible alternative in systems with medium-strong coupling and non-Markovian systems, in the future. However, the algorithm obtained, due to the high requirements of gates and qubits, is not scalable to real world photosynthetic systems, with the current state of quantum technology. Hence, this simulation should be seen as a proof of concept, since a realistic quantum simulation of a photosynthetic system would have to involve hundreds of light-harvesting molecules, which is beyond the current quantum technology. Furthermore, the algorithm only effectively simulates pure-dephasing baths. For future work, we aim at extending it to new types of bath, e.g. those allowing for higher exciton recombination rates and non-Markovian effects as well as to new geometries of photosynthetic systems, in particular, to the Fenna-Matthews-Olson complex [35].

A.1 -Brief introduction to quantum systems
There are several formalisms in which quantum mechanics can be expressed, including the one making use of Hilbert spaces and density matrices [38], which we shall denote as Hilbert for clarity. In the Hilbert formalism, the theory is expressed in terms of states, |Ψ , and transitions between them. The states may be expanded on an eigenbasis of the Hilbert space, where i |λ i | 2 = 1. The eigenbasis states are often taken as the eigenstates of a Hamiltonian and known in quantum mechanics as the stationary states, while |Ψ is an example of the so-called superposition states.
Furthermore, quantum systems can also be composed and form new systems through the Kronecker product of the Hilbert spaces, corresponding to two sub-systems (in this case 1 and 2): The analysis of this operator quickly leads to another very relevant property of quantum mechanics, which is the existence of non-separable states, i.e. states that cannot be written as a Cartesian separable product of the subsystems 1 and 2, as for instance happens in the so-called Bell states, which are a possible state of the ⊗ on two qubits. This phenomenon is the so-called entanglement, a vital component of many quantum technologies such as quantum information. An extension of pure states represented by vectors in the Hilbert space leads to the density matrices, which can also describe the so called mixed states. A mixed state can be considered as an mixture of several pure states, with certain statistical weights, therefore it cannot be described by a vector in the Hilbert space. A (Hermitian) density matrix (or operator) can be written aŝ where P α are real numbers between 0 and 1 and |Φ α Φ α | is known as the projection operator of a state Φ α . For a pure state all Ψ α = 0 except for one, which is equal to unity. By expanding Φ α on a basis of stationary states, one can obtain the matrix representation ofρ in that basis. The density matrix of a system composed of two parts, A and B, allows for a simple criterion for checking whether these parts are entangled or not in a particular state,ρ A+B (see, e.g. [4]). To distinguish between entangled and separable states, one needs to calculate the density matrix of subsystem A (or B) as the partial trace ofρ A+B and check the trace ofρ 2 A . If it is equal to unity, the subsystem A is in a pure state and the whole system is separable.
While the Hamiltonian formalism of quantum mechanics permits properly describe closed systems (and, in principle, universe as a whole), in practice, for instance in quantum technology, one deals with systems that are not closed, i.e. there is a free uncontrollable element, not necessarily conservative, interacting the system under study. The theory that studies these kind of systems is a branch of quantum mechanics known as open quantum systems [5]. In fact, an open system is nothing more than a closed system composed of two subsystems, the system and the environment, in which the system plus environment Hilbert space is expressed as H S+E . The overall effect of the environment in the system can be obtained by a tracing operation over the environment Hilbert space: A.2 -Solution of the Schrödinger and Liouville equations for the two-molecule system Making use of the Schrödinger equation ( = 1): and writing the state |Ψ in the site basis as: where |α| 2 + |β| 2 = 1. One gets by applying the Schrödinger equation to the state |Ψ and multiplying it by 0| and 1|: iα = 0 α + Jβ (26) By making the substitions α = ae −i 0t and β = be −i 1t so that: one has, taking into account the substitutions and Eqs. (26) and (27): where w = 1 − 0 . Then the variable c is introduced as c = ae iwt . One can observe through Eq. (30) thaṫ a = (ċ − iwc)e −iwt so that i(ċ − iwc) = Jb. Differentiating with respect to time and using Eq. (31), one has: By seeking solutions as c = e iλt : − λ 2 + wλ + J 2 = 0 one gets: where Ω = √ w 2 + 4J 2 . The solution of Eq. (32) is: then, a(t) = Ae i Ω−w 2 t + Be −i Ω+w 2 t (36) and b(t) is found by using equation (31): The constants A and B depend on the initial conditions. If the donor molecule is the only one excited (|0 ) at the initial time t = 0, then α(0) = 1 and β(0) = 0. Thus, a(0) = 1 and b(0) = 0 and one gets A = Ω+w 2Ω and B = Ω−w 2Ω . The solutions are: and the system's density matrix becomes: The time dependence of the non-diagonal elements of the density matrix in the energy basis can be obtained much easier. Consider the quantum master Liouville equation applied to a closed system ρ subjected to a Hamiltonian H, defined as: Writing the density matrix in the energy eigenbasis (where the Hamiltonian is diagonal) yields: where E i is the energy eigenvalue of the ith energy eigenstate. The populations (matrix elements ρ ij , where i = j) remain constant and the coherence terms (matrix elements ρ ij , where i = j) evolve oscillating in time as:

A.3 -Circuits for the coherent regime
This annex contains the optimized circuits used to build the quantum simulations of section No decoherence regime for the near resonance (figure 10) and off resonance (figure 11) regimes. Figure 10: Optimized quantum circuit for the near-resonance system simulation. Figure 11: Optimized quantum circuit for the non-resonance system simulation.