#### Abstract

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, IBM *Q*, of 5 qubits. We validate the results with the Haken-Ströbl model and discuss the influence of environmental parameters on the efficiency of the energy transport.

#### 1. 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 [1, 2]. 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 [3]. The understanding of such phenomena has greatly progressed in the past 40 years with the physical characterization of the structure of many photosynthetic complexes [4–6]. 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 [7] to the construction of artificial light-harvesting devices and solar fuels [8–11].

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 [2]. 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 center, which proceeds in the form of the so-called Excitonic Energy Transfer (EET), as schematically shown in Figure 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% [2]. 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 equation (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, that is, coherent over some period of time, as evidenced by experimental observations of long-lived oscillatory features in the dynamical response of several photosynthetic systems [12–14]. Moreover, it is strongly influenced by the environment. The donor-acceptor pair is not isolated from the rest of the world and is an example of so-called open systems [15]. 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 [2, 16], which can be extended to medium coupling strengths and non-Markovian regimes by including higher-order system-bath interaction terms [17, 18]. For the latter regime, there are also nonperturbative techniques based on the use of path integrals to dissipative systems [19], which can be used to create sets of solvable systems of hierarchical equations, the so-called hierarchical equations of motion (HEOM) [20–22]. 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 nonlocal 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 molecules with a theoretical tool such as the HEOM grows exponentially with . A possible computational solution that 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 Feynman [23, 24]. 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 [25]. Closer to the present work, there are works on quantum transport [26, 27] and on the quantum simulation of dissipative systems [28]. Particularly on the quantum simulation of photosynthesis, we would like to highlight the papers in [29–31], using superconducting qubits and [32] employing a Nuclear Magnetic Resonance (NMR) simulator [33]. 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 [34], within the NMR system. In this case, the implementation is specific for the EET (i.e., nonuniversal), and the model Hamiltonian was extracted from spectroscopic data for a photosynthetic system [35]. We are simulating the same Hamiltonian as in [32] 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 [36]. 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.

#### 2. Physics of the Energy Transport in Photosynthesis

##### 2.1. 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 [37] or artificial nanostructures such as quantum dots [38] if appropriate conditions are met, which were first formulated by Förster [39]:(i)The distance between the donor and acceptor molecules must be sufficiently small because the transfer probability decreases quickly with the distance between them , usually as (ii)There must be a resonance between the excited states of the donor and acceptor molecules (“Resonance” here means that the energy spectra of the two molecules, broadened because of a number of natural reasons, overlap-see Figure 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 molecules, 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 [2]:where is the coupling constant, is the angular frequency of the electromagnetic field, and and 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 [37]where is the refractive index of the medium, is the transient dipole moment of the donor (acceptor) molecule, , is the radius vector between the two molecules, and the angular brackets stand for angular average over different orientations of the dipoles.

Even though equation (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:(i)If the difference between the energies of the excited state energy of the donor and acceptor molecules is small, , and these states are in resonance, the energy transfer between the molecules can occur with a high probability(ii)If (off-resonance), 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 transfers 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 [40–42]. The Förster theory does not apply in this regime, as it simply ignores coherence. Later, in 1957, Redfield [43] proposed a transport theory, which applies to the opposite regime of strong coupling between the donor and acceptor [2, 40, 44] (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 denotes the exciton on the molecule . The eigenstates of (4) are linear combinations of and . Coherent dynamics correspond to the presence of nonzero off-diagonal elements in the density matrix describing the evolution of the quantum system. Their oscillation (or quantum beating) is indicative of coherence [2]. 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 as , 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, for example, in [16]. This consideration is extendable to a chain of molecules and can be seen as (partially) coherent transport [40]. 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 30s [45] and theoretically predicted in more recent works [44, 46]. It became possible to observe them more recently, thanks to the advances of optical spectroscopy techniques [12–14, 47], and it was achieved even at room temperature [42]. In these experiments, it was possible to confirm the substantial impact of such coherent effects on the excitation energy transfer in photosynthetic systems [2]. Moreover, the importance of environmental noise in the quantum transport involving coherence was also discussed more recently [15, 48] and it is not fully understood yet.

##### 2.2. Decoherence

Processes caused by the molecules’ environment may destroy coherence and thus influence this type of energy transport [2, 15, 49]; 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 in 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 [50].

###### 2.2.1. 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) [51]. For a two-level system (e.g., a qubit), this type of decoherence contracts the Bloch sphere along the -axis (see Supplementary Information A.1).

###### 2.2.2. 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; that is, 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 plane.

###### 2.2.3. Depolarization

This type of decoherence changes system’s state, which initially is pure, to a mixed state, with a probability of another pure state and the probability 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 ). 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 suppression). It eventually turns the diagonal matrix elements (populations) into (noncorrelated) 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 [2]. First, it yields random fluctuations in the energy spectrum of each molecule, which can bridge the energy gap between the molecules, momentarily turning a nonresonant system into a resonant one. Second, dephasing can also help in 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 [2]. 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.

#### 3. 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 following equation:for a time-independent Hamiltonian. Here and 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, , 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 , for times, of the computational Hamiltonian, such that . The process is finished by the observation of the desired properties, that is, a set of measurements, in the appropriate basis, on the end state.

Concerning the particular qubit encoding chosen, a chain of molecules is encoded by a set of qubits, where corresponds to the excitation (exciton) on the -th molecule; for example, for a two-molecule chain, state represents the exciton on the first molecule and represents the excitation on the second one, and a possible successful transport of energy would correspond to the transition of the state to the state . 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 . Also, it is convenient to measure the energies/frequencies in , as this is common in spectroscopy.

##### 3.1. No-Decoherence Hamiltonian

Considering a small chain of molecules, the system’s Hamiltonian in the site basis reads as follows:where is the first excited state energy of the molecule and is the electronic coupling between the molecules and . The Hamiltonian (7) for just two molecules (1 qubit), identical to equation (4), in the matrix form, reads

Its evolution operator is given by

Although the Hamiltonian (8) possesses nondiagonal elements, finding a good approximation in terms of quantum circuits is relatively straightforward. A possible strategy for this is by finding a diagonalizing transformation, , of the Hamiltonian, such thatwhere is the diagonal Hamiltonian. Therefore, the evolution operator can be rewritten as follows:

The problem now reduces to the approximation of the operator (and its adjoint) and the Hamiltonian , 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 as

The and matrices can be implemented by simple rotations, and , for a two-molecule system. However, for a higher number of molecules, a rotational decomposition algorithm together with the Gray code [51], 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 molecules is [51]. On the other hand, the diagonalized evolution operator,translates into trivial phase rotations over each of the energy eigenstates of the system with the respective energy eigenvalues . This operator can be constructed as a sequence of gates applied to an ancilla qubit (initialized at ), where the angle is given by , . The 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 molecules, is .

For the whole circuit, resulting from the sequencing of , the number of qubits required to simulate a molecular chain of elements is and the gate count scales with single qubit and CNOT gates. The transformations and , 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.

##### 3.2. 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 [15, 44, 50]. The actual algorithm to be used is the one of [52], 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 readsand it consists of the system Hamiltonian, , of the previous section and the perturbation of a bistable fluctuator environment, . The latter simply shifts the energy by a constant value for each molecule, , as illustrated in Figure 3. Explicitly,where is the projection operator and, considering one fluctuator interacting with each molecule ,

The function switches the fluctuator between the positive and negative values (appearing randomly) at a given fixed rate and is the fluctuation strength (or the coupling strength to a molecule ). Physically, the action of the fluctuations is typically stronger for the excited states [44, 53] and can be larger than the donor-acceptor coupling .

The implementation of such random bivalued function can be done in a straightforward way by a classical pseudorandom numbers generator with a probability of of the values and . 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 [54]. Under this approximation, the unitary evolution operator of the Hamiltonian, for a time , where is the number of iterations and is the iteration time-step, becomeswhere denote the eigenvalues of the system Hamiltonian.

Note that the projection operator is not present in the evolution operator (17) because the latter is used in its eigenbasis, that is, the site basis. The fluctuator interaction evolution operator is a selective rotational gate over a molecule , which can be implemented by a set of *X* gates and a controlled gate with angle , applied over an ancilla qubit initialized at . The whole circuit is presented in Figure 4 for one iteration. The fluctuator waiting time (interval of time between switches), that is, , can only be equal to or higher than the iteration time-step, . The switching in the fluctuator-molecule coupling strength is performed at every iteration, where .

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 [51, 52]; however, in this case, the open system is simulated in a closed form so, similar 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 , using an iteration time-step , 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: single qubit and CNOT gates, where is the number of molecules and is the number of fluctuators interacting with each one.

In the implementation of the system with decoherence, the algorithm gate resources complexity is for a single run. This simulation, yet again, possesses a very high circuit depth, which makes its application unfeasible in quantum computers. The number of necessary qubits is the same as that in the no-decoherence simulation .

It also requires random numbers to be classically generated, where is the number of runs of the algorithm and is the switching rate of the fluctuator interacting with the molecule. The number of required simulation runs to average the results and obtain an error is predicted to scale as . This complexity is calculated based on the possible nondegenerate energy state outcomes of the entire chain in the simulation for a time . These outcomes are caused by the bistable random fluctuations; therefore, the possible nondegenerate energy state outcomes for each molecule obey a discrete Gaussian probability distribution.

#### 4. 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, IBM *Q* of 5 qubits, while the pure dephasing scenario was simulated on the QASM quantum simulator in both the near-resonant and nonresonant 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 [15, 55]. The simulations and circuits involved, encoded in the Qiskit platform [36], can be performed in the following url: https://github.com/jakumin/Photosynthesis-quantum-simulation.

##### 4.1. Coherent Regime

The scenario for this regime was simulated with a simple chain of two molecules. As discussed in Materials and Methods and using the parameters proposed in [32], we define the system’s Hamiltonian as follows:

Near-resonant regime:

Nonresonant regime:The results for both regimes were obtained using an actual quantum device (IBMQ London of 5 qubits) and can be seen in Figures 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 and , respectively.Taking the fluctuator’s switching rate to be or the fluctuator-molecule coupling strength to be , one has the coherent regime. These simulations show the limiting case of the Redfield regime, that is, the very weak system-environment coupling, . 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 [56].

In the performed simulations, the near-resonant and nonresonant regimes have a maximum probability of and , respectively, of the energy being transferred to the acceptor molecule. Using the quantum Liouville equation [2] (see Supplementary Information A.2), the period of the quantum beating is for the near-resonant regime and for the nonresonant regime. These periods are in the femtosecond timescale of the experimentally observable quantum beatings [12, 41, 42]. The simulation results show a similar behaviour to that of those predicted by the Schrödinger and quantum Liouville equations, where the off-curve points are predominantly originated by errors in the quantum hardware.

##### 4.2. 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 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 Lindblad master equation [15, 55]. The Lindblad equations were solved in a classical computer using QuTiP [57], a quantum open-systems software framework. The set of Lindblad equations, correspondent to the model in this setting, had one free parameter regarding the environment, the dephasing rate, . The Lindblad equation in the Haken-Ströbl model readswhere are the Lindblad operators, responsible for the system-environment interaction. The system Hamiltonian, , is given by the matrix in (18) for the near-resonant system and the matrix in (19) for the nonresonant 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 evolutions 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 . As mentioned above, the dephasing rate, , for the Lindblad equation is adjusted to the behaviour of the system under the action of a fluctuation strength . For the fluctuation strength in the algorithm varying within the interval , and the corresponding dephasing rate of the Haken-Ströbl model lies in the range. Due to the existence of random fluctuations, a large number of samples had to be generated. The algorithm was implemented with 250 runs, where 5000 shots were performed for each time . Figures 7 and 8 present the simulation results for different values of the fluctuation strength, along with the theoretical evolution dynamics, for the near-resonant and nonresonant systems, respectively.

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(d)**

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 . Here is the Boltzmann constant, is the temperature of the bath, and is a normalization constant [50]. Taking the limit of very high temperatures, the population terms approach the Boltzmann distribution , which is compatible with the results obtained. The relaxation cannot be fully observed in Figures 7(a), 8(a), and 8(b) 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 times larger than the transfer rate, (i.e., the fluctuator waiting time must be shorter than ). 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 (Figures 5 and 6).

The time that coherence lasts in the system is essentially defined by the fluctuation strength, : in Figures 7(a), 7(b), 8(a), and 8(b) (lower ) the coherence is maintained for some time, while in Figures 7(c), 7(d), 8(c), and 8(d) (higher ) 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 nonoscillating behaviour) is about in Figure 7(b) (near-resonant system) and in Figure 8(b) (nonresonant system), with a fluctuation strength . At a longer time, it has been experimentally observed to persist ( [42]), a timescale which could be modeled in the present simulation by changing the environment parameters, that is, lowering the fluctuation strength , as can be observed in Figures 7(a) and 8(a).

#### 5. Discussion and Conclusions

Two main conclusions can be drawn from the presented results:(i)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.(ii)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 implies an increase in the fluctuation strengths; thus, a faster suppression of the quantum beatings can be observed, as predicted theoretically [15].

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 [32] due to the different timescales used. The major difference lies in the physical implementation, NMR versus universal quantum computer, where there might be an advantage for the former from the viewpoint of the scalability and reliability, at the current state of quantum technology. However, there is a clear advantage of the quantum computer, from the viewpoint 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 may be problematic, as a matrix diagonalization operation is necessary (complexity estimated in ).

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, 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 defined here, it was clear that dephasing, modeled 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 to the energy transfer efficiency, allowing the simulation of different transport regimes by setting them appropriately.

Similar to [32], this setting revealed itself as an interesting platform for the study of 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, that is, 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 [2].

#### Data Availability

The data, codes, and other tools used in this work are available upon request.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest.

#### Acknowledgments

Carlos Tavares was funded by Fundação para a Ciência e a Tecnologia (FCT) under Grant no. SFRH/BD/116367/2016, funded under the POCH programme and MCTES national funds. This work was financed by the European Regional Development Fund (ERDF) through the Operational Programme for Competitiveness and Internationalisation (COMPETE 2020 Programme) and by national funds through the Portuguese funding agency, Fundação para a Ciência e a Tecnologia (FCT), within project KLEE (POCI-01-0145-FEDER-030947) and the Strategic Funding UIDB/04650/2020 of the Centre of Physics.

#### Supplementary Materials

A.1: brief introduction to quantum systems. A.2: solution of the Schrödinger and Liouville equations for the two-molecule system. A.3: circuits for the coherent regime.* (Supplementary Materials)*