A general quantum algorithm for open quantum dynamics demonstrated with the Fenna-Matthews-Olson complex

Using quantum algorithms to simulate complex physical processes and correlations in quantum matter has been a major direction of quantum computing research, towards the promise of a quantum advantage over classical approaches. In this work we develop a generalized quantum algorithm to simulate any dynamical process represented by either the operator sum representation or the Lindblad master equation. We then demonstrate the quantum algorithm by simulating the dynamics of the Fenna-Matthews-Olson (FMO) complex on the IBM QASM quantum simulator. This work represents a first demonstration of a quantum algorithm for open quantum dynamics with a moderately sophisticated dynamical process involving a realistic biological structure. We discuss the complexity of the quantum algorithm relative to the classical method for the same purpose, presenting a decisive query complexity advantage of the quantum approach based on the unique property of quantum measurement.


Introduction
Simulating physical processes with quantum algorithms has been a major focus of quantum computing research [1][2][3][4][5][6]. Open quantum dynamics studies the time evolution of a quantum system interacting with an environment [7]. The complexity of the environment often makes exact treatment impractical, and various approximation approaches have been developed to treat Sabre Kais: kais@purdue.edu the environment as an averaging effect on the system of interest. So far relatively few studies have been done to develop quantum algorithms for open quantum dynamics [8][9][10][11][12][12][13][14][15] despite its importance in modeling realistic physical systems. One main challenge is the time evolution of open quantum systems is non-unitary, while quantum algorithms are realized by unitary quantum gates. To tackle this problem an early study [8] added an auxiliary environment to the system thus making the total evolution unitary. Mimicking the Markovian process, that procedure required a reset of the auxiliary environment at each time step, which could be difficult to achieve if the system and the environment are entangled or the evolution time is long. Other studies [9,10] have proposed novel quantum algorithms to simulate open quantum dynamics, yet the demonstrations are limited to basic 2-level systems. In this work we propose a general quantum algorithm that can be demonstrated with a much more complex and realistic physical model. In a previous study we have developed a quantum algorithm for evolving the dynamics in the operator sum representation [16]. The algorithm uses the Sz.-Nagy dilation approach to convert each Kraus operator into its corresponding unitary dilation matrix, which is then implemented on a quantum circuit. The quantum algorithm was successfully applied to the amplitude damping quantum channel and implemented on the IBM Q quantum simulator and quantum computers [16]. The quantum algorithm was designed with generality in mind and indeed has been adapted through use of an ensemble of Lindbladian trajectories method [17,18] to simulate non-Markovian dynamics as demonstrated on the Jaynes-Cummings model [19]. However, there are two issues that must be solved before our dilation-based quantum algorithm can be applied to more complex dynamical processes. Firstly, the Kraus operators in the amplitude damping model have an explicit dependence on time, which allows each time step to be simulated independently. This explicit dependence however cannot be easily obtained for more complex dynamical models. Secondly, a large class of dynamical models are described by the master equation formulation. To apply our quantum algorithm to such models a connection between the master equation and the operator sum representation must be established. Although there have been studies of converting a widely used type of master equation -the Lindblad equation -into the operator sum representation [20,21], those methods solve the master equation into a time-explicit matrix form and therefore require integration of the superoperators describing the dynamics. Integration of the superoperators is hard for classical algorithms and currently not solvable with any quantum algorithm. In the following we intentionally avoid integrating the superoperators, but instead treat the master equation with the Euler method -this is an essential simplification that allows practical application of the quantum algorithm to more sophisticated dynamical models.
In this work we present a generalized quantum algorithm that can simulate any open quantum dynamics represented by either the operator sum representation or the Lindblad master equation. In this algorithm each Kraus operator can be directly related to a physical process without the integration of superoperators. In addition, the Kraus operators no longer require an explicit dependence on time and the system at a given time is evolved from the initial state by iteratively applying the Kraus operators obtained from the physical processes. To demonstrate the generality of the generalized algorithm, we use it to simulate a moderately sophisticated dynamical process for a realistic biological structure: the Fenna-Matthews-Olson (FMO) complex. The Fenna-Mathews-Olson complex is a trimeric-pigment protein complex found in green sulfur bacteria [22]. In the photosynthetic lightharvesting process, it is responsible for transferring excitonic energy from the antennae complexes to the reaction center with nearly 100% efficiency [23][24][25]. While this is a very well-studied complex [26][27][28][29][30][31][32][33], a more in-depth understanding of this transport process can provide valuable insight for the dynamics of other light harvesting complexes or for the design of artificial photosynthetic systems [34,35].
While some earlier studies exist [36][37][38], the application of our generalized quantum algorithm to the FMO dynamics as demonstrated on the IBM QASM simulator [39] is so far as we know the first successful quantum simulation of a moderately sophisticated dynamical process involving a realistic biological structure. The simulation showcases the generality of the method by mapping the Lindblad master equation of the FMO into the operator sum representation. The same as previously discussed in Ref. [16], in terms of the gate count required to execute the evolution, the generalized quantum algorithm has comparable computational complexity as classical methods. However, here we emphasize that both the previous and the generalized quantum algorithms can have a decisive query complexity advantage over any classical method, when evaluating an observable over the density matrix. Under specific conditions the query complexity advantage can even translate to the total complexity advantage. This complexity advantage is the result of the fundamentally quantum properties of superposition and projection measurement.

Generalized Quantum Algorithm
The basic mechanisms of the generalized quantum algorithm -including the Sz.-Nagy unitary dilation procedures and the observable evaluation procedures -have been introduced in our previous work [16] and reviewed briefly below: We assume the physical composition of the initial density matrix is known and can be expressed by a sum of different pure quantum states weighted by the corresponding probabilities: where each p i is the probability of finding each |φ i in the mixed state of ρ. Now if the dynamical model is given by the operator sum representation: we want to simulate the time evolution of ρ(t) given the initial ρ and the Kraus operators M k 's. To achieve this, we can prepare each |φ i as an input state vector v i in a given basis and then use a quantum circuit to create the quantum state: [16,40], After |φ ik (t) has been obtained for each M k and each v i , the population of each basis state in the current basis can be obtained by calculating the diagonal vector: where diag(|φ ik (t) φ ik (t)|) can be efficiently obtained by applying projection measurements on the first half subspace of U M k (v T i , 0, ..., 0) T . To evaluate the expectation value of an observable A = Tr(Aρ(t)), we consider the operator A = A+I||A|| 2||A|| with the Cholesky decompositioñ A = LL † [41]. We can then evolve: .., 0) T (5) and obtain Ã by: [16]. After this, A can then be obtained from Ã by This procedure was then applied to the amplitude damping channel: where M 0 and M 1 both have an explicit dependence on time t. This concludes the review of the previous algorithm, for more details please see Ref. [16]. Now to introduce the new features that allow the generalized algorithm to simulate more complex dynamics, we first note that the Kraus operators with explicit time dependence as reviewed in Equation (7) are not generally available for more complex dynamics. This is especially true when the dynamics is described by a master equation, where the superoperators have to be integrated first to obtain the Kraus operators with explicit time dependence. To avoid the hard problem of integrating superoperators, a naive and incorrect idea is to directly write the Kraus operators of an arbitrary dynamics in the same basic forms of those used in Equation (7) simulating the amplitude damping channel. However as shown in the Supplementary Information Section S1 for the finite-temperature amplitude damping channel, Kraus operators with naive dependence on time lead to incorrect dynamics for even such a simple deviation from the amplitude damping channel. Indeed, the original formulation of the operator sum representation [42] can be essentially understood as a single physical event that starts with an initial density matrix ρ(0) and ends with a final density matrix ρ(1): , where E 1 is the quantum operation described by the collection of M 1k 's. Now if we use ρ(1) as the initial state and apply another quantum operation, we have ρ (2) 2k . This process can be repeated iteratively such that: where each ρ(s) with s = 1, 2, ..., S can be understood as a discrete time step that samples the dynamics until we reach the final time step. Now an outstanding question is how to determine the Kraus operators M sk for each quantum operation E s . The answer is if the time evolution involves the same physical processes throughout the total time interval being simulated, then we can use the physical processes to determine the collection of M k 's and use the same collection for all the iterations. For example, in the amplitude damping channel, the physical process is a single transition from the excited state to the ground state, and the corresponding Kraus ma- , where p is the probabil-ity of the transition happening over a single time step. The condition k M † k M k = I thus leads Applying the iterative process in Equation (8) to an initial ρ(0) with the same M 0 and M 1 for each iteration we have the excited state population at the final time step: where we assume p is given by a rate constant γ times a small time interval δt and t = Sδt. Equation (9) demonstrates that when the Kraus operators obtained from the physical processes are applied iteratively to the initial density matrix, we obtain the expected exponential decay for the excited population under the condition that the time interval δt between each time step is small compared to the total time t. It is easy to verify that the correct dynamics is also obtained for all other entries of the density matrix. In the Supplementary Information Section S1 the same procedure also produces the correct dynamics for the finite temperature amplitude damping. The relation between the transition probability and the rate constant in Equation (9) also points to a natural relation between the operator sum representation and the Lindblad master equation. Indeed the Lindblad master equation with the form of describes each physical process with the Lindblad operator L k and the rate constant γ k (the {·, ·} is an anticommutator). Now if we associate each L k with an M k , with the details of the physical process described by the form of the operator L k , and the probability of the physical process given by p k = γ k δt, then [7]: as the initial state and incrementing time by δt we get: Equation (12) converges to Equation (10) when δt → 0. At the same time: Therefore the condition of k M † k M k = I is approximately satisfied for the operator sum representation. Equations (10) through (13) show that the Lindblad master equation and the operator sum representation are two different de-scriptions of the same dynamics: when the physical processes are represented by the Lindblad operators L k 's, the Lindblad master equation defines a differential equation on the density operator ρ(t); when the physical processes are repre-sented by the Kraus operators M k 's, the operator sum representation essentially evolves the differential equation with the Euler method. In actual simulation of the dynamics, Equation (11) is especially important as it allows us to easily determine the Kraus operators M k 's with k > 0 for use in the iterative process in Equation (8).
In Equation (13) the condition k M † k M k = I is not exactly satisfied, which may cause problems in our specific quantum algorithm. However this problem is easy to solve as we can enforce the condition by defining: Now with all the Kraus operators defined, we can simulate the dynamics by the iterative process in Equation (8) with the same Kraus operators for each iteration. The first three time steps will look like: For the first time step we can apply the basic unitary dilation procedure described in Equations (1) to (6) or in Ref. [16] to implement each M k term. For the second and later time steps we need to use higher dilations [16,40] to implement each term, e.g. a 2-dilation for each M j M k and a 3-dilation for each M i M j M k . In Equation (15) if the total number of Kraus operators is Kmeaning there are (K −1) physical processes -for the first time step we need to implement K M k terms, for the second step K 2 M j M k terms, for the third step K 3 M i M j M k terms, and in general we need to implement K s terms for the s th time step. This exponential scaling of terms with the number of time steps apparently limits the total number of time steps that can be simulated. In actual simulation however, there are two ways the number of terms can be significantly reduced such that the total number of time steps simulated can be greatly increased (See the Methods section for details). As shown in Equation (12) the time interval δt between time steps must be kept small for the Euler method to work, therefore increasing the total number of time steps is equivalent to increasing the total time for which the dynamics can be evolved. In the following our quantum simulation covers sufficient total time to observe the FMO dynamics. If after all the simplifications the total time is still too short for observing some other physical phenomena, quantum tomography may be applied to the final S th ρ(S) which can be then used as the new initial state. Quantum tomography however will likely take a lot of computational resource and will be beyond the scope of the study.

Simulation of the Fenna-Matthews-Olson Dynamics
In the photosynthetic light-harvesting process, the FMO is a trimer complex which acts as a quantum wire connecting the light-harvesting antennae to the reaction center [43,44]. Each monomer consists of seven bacteriochlorophyll chromophores, where an initial excitation can occur on either chromophore 1 or 6 and is transported to chromophore 3, which is closely coupled to the reaction center [45,46]. This transport is largely driven by environmental interactions including those with neighboring excitations and the protective protein scaffold [47]. Within each monomer of this trimer, there are multiple efficient quantum pathways for the exciton to be transferred to the reaction center. This quantum redundancy has led to a study of the different functional subsystems of the FMO's chromophores, demonstrating that many subsets ex-ist with similar efficiencies to that of the entire monomer [46]. In this work, we consider the three-chromophore subsystem which consists of chromophores 1-3 and has been shown to faithfully reproduce the exciton dynamics on the entire seven-chromophore monomer [46]. A single FMO monomer is shown in Figure 1, where the functional subsystem is depicted in green with the interactions depicted with yellow arrows. Con- sidering the redundancy of the chromophore subunits in the FMO dynamics, we can focus on three chromophore subunits supporting three local states of excitation. Adding a ground state and a sink state, the Hamiltonian of the relevant part of the FMO is: where the state |i with energy ω i is created by the Pauli raising operator σ + i and annihilated by the Pauli lowering operator σ − i , J ij is the coupling strength between |i and |j . The FMO dynamics can be described by the Lindblad master equation: where the [·, ·] is a commutator and the {·, ·} is an anticommutator, and the seven L k 's represent seven physical processes (the rate γ k has been folded into each L k as defined next).  (11) and (14): where the time interval δt is set to be 2000 atomic unit or 48.4 fs. We then apply the procedure in Equation (15) for up to the 6 th time step with a total simulation time of 12000 atomic unit or 290 fs. In principle, with eight Kraus operators, the total number of terms from Equation (15) for the 6 th iteration would be a prohibitively large number 8 6 = 262144. However by considering redundancy of the terms and setting a threshold of the norm of the matrices to be > 0.01 (reduction of terms discussed in the Methods section), we are able to reduce the number of terms to only 679, which is much more manageable. We then construct the quantum circuits from the unitary dilations of the Kraus matrices(for detailed procedures see Equations (1) to (6) or Ref. [16]) and run the circuits on the IBM QASM simulator [39]. An example of the circuits is shown in the Supplementary Section S5. The results are shown in Figure 2 and diagonal elements of the evolved density matrix are obtained by projection measurements into the computational subspace. We first simulate six data points at 2000 (48.4 fs), 4000 (96.8 fs), 6000 (145 fs), 8000 (194 fs), 10000 (242 fs) and 12000 (290 fs) atomic unit respectively. To obtain more data points such that the oscillations can be smoothly represented, we simulate four additional groups with six data points each (details in the Methods section). The results simulated by our generalized quantum algorithm as implemented on the IBM QASM simulator are shown as dots and they agree well with the classically calculated results shown as curves. The results in Figure 2 demonstrate the viability of the generalized quantum algorithm. In particular the iterative procedure in Equation (8) , the formulation of Kraus operators from Lindblad operators in Equations (11) and (14), and the simplification of terms by norm threshold and redundancy, are working together to produce the correct population dynamics with small errors. In Figure 2 we can see the excitation beating between chromophores 1 and 2, as it gradually decays into the sink. This process is driven by a combination of environmental noise and entanglement between the chromophores in the functional subsystem.
In Figure 3, we see the expectation values of the energy observable calculated at different time steps. The evaluation of an observable by projection measurement is a complex process as explained in Equations (5), (6) and Ref. [16]. Again, the quantum simulation results shown as dots agree well with the classical results shown as a curve. Indeed, the way we evaluate the observable in the generalized algorithm is the same as the previous quantum algorithm [16], therefore if the new algorithm can accurately simulate the evolution of the density matrix, then the observable evluation will not introduce new errors not considered before.

The complexity scaling
On the complexity scaling of the quantum algorithm, there are two separate scalings for the evolution of the dynamics: the scaling with the system size n and the scaling with the evolution time step.
For the scaling with the system size n, for an arbitrary n × n Kraus operator M k without any special property, the cost to realize its unitary dilation U M k is O n 2 [16] . In practice however, an M k is often sparse with few non-zero elements. This is particularly true for our current quantum algorithm as we generate the M k 's from the basic physical processes of transition and dephasing. Consequently the complexity scaling of a quantum circuit implementation for each M k can be greatly reduced to O log 2 n [16]. Now we still need to evolve all K number of M k 's, so the total complexity scaling is O Klog 2 n , where K is determined case-by-case by the dynamical model. Note that the different M k 's can be evolved in parallel, thus the scaling in K is a "soft" scaling because it does not contribute to either the depth or the width of each individual quantum circuit.
For the scaling with the evolution time step, without any simplification, Step s requires that K s matrices be evolved, which is an exponential scaling in the time step. Fortunately, the number of terms can be reduced drastically by the sparsity and redundancy of the M k matrices that we discussed in detail under Equation (15) and in the Methods section. In general it is difficult to claim the scaling is polynomial or exponential, because the number of matrix terms depends, in a case-by-case matter, on the actual forms of the M k matrices and the threshold set for the matrix norm. For our particular FMO model we have managed to simplify the number of terms from the prohibitive 8 6 = 262144 to merely 679, which allowed us to implement it on the IBM simulator. Here we remark that the scaling with the evolution time step is the natural result of solving a differential equation with the Euler method. The Euler method is the simplest among the Runge-Kutta methods that are standard tools for the numerical evolution of differential equations. Consequently, the scaling with the evolution time is not an artifact introduced by our method, but an important issue to be investigated if one wishes to adapt classical numerical tools for differential equations to the quantum computing field. With case-specific simplifications to reduce the scaling, our method represents a step in this direction.
So far we have discussed the two complexity scalings for the evolution of the dynamics. In the quantum computing context, the complexity scaling associated with the measurement and retrieval of results -the query complexity -is also important. Here we emphasize that both the previous [16,19] and the current generalized versions of the quantum algorithm can have a decisive polynomial-versus-exponential query complexity advantage over classical methods when used to evaluate an observable over the density matrix. In addition, in specific situations where the evolution complexity is below a threshold, the total complexity is dominated by the query complexity, and therefore the quantum algorithm can have a decisive total complexity advantage over any classical methods. This quantum advantage along with an example system is discussed in greater detail in the Supplementary Information Section S4.

The error analysis
There are two main sources of errors in the current study. The first source is the error introduced by evolving the master equation with the Euler method and this error is well known to be O δt 2 . As a future direction, this error could be potentially reduced by introducing high-order methods such as the Runge-Kutta methods. Note that this source of error is not unique to our quantum algorithm but is also present in many classical algorithms used to numerically solve a differential equation. The second source is the measurement error when the IBM QASM simulator simulates a quantum projection measurement. This error is notably only scaling inversely with the number of measurements taken, and not scaling with the system size, as we have discussed in Section S3 of the Supplementary Information. This error is also not unique to our quantum algorithm but always present when one wants to extract physical information from the evolved density matrix. Overall, these two errors did not prevent us from obtaining reasonably good results compared to the classical benchmarks. A third source of error -the gate error -would become important when one tries to run the quantum algorithm on an actual quantum computer. Unfortunately the quantum computers available to us at the writing of the manuscript have non-trivial gate errors such that an experimental demonstration of our quantum algorithm has not been achieved. Nonetheless, the demonstration on the IBM QASM simulator shows that our quantum algorithm is fully programmable (without blind spots or "oracles") and cost-efficient (with casespecific simplifications), such that it can be implemented on quantum computers with smaller gate error or better error-mitigation capabilities in the near future.

Conclusions and Outlook
In this work we have developed a generalized quantum algorithm for open quantum dynamics to simulate more complex dynamical models. In particular we formulate the Kraus operators from basic physical processes that constitute the dynamics, and then realize the time evolution by an iterative process. We then relate the Lindblad operators to the Kraus operators in a simple manner with the Lindblad master equation being connected to the operator sum representation through the correspondence between a differential equation and its integration by the Euler method. The generalized quantum algorithm works with any dynamical model represented by the operator sum representation or the Lindblad master equation, and thus is much more general than the previous algorithm. We demonstrate the generalized algorithm on the FMO dynamical model using the IBM QASM simulator, which is so far as we know the first successful quantum simulation of a moderately sophisticated dynamical process involving a realistic biological structure. Finally we analyze the query complexity of the algorithm, and constructs an example where the quantum algorithm showcases decisive complexity advantage over any classical method. (15) To reduce the number of terms in Equation (15) and increase the total simulation time, the first way is to notice that the M k 's with k > 0 typically correspond to basic physical processes such as state transition or dephasing, and thus an M k typically contains only one or few non-zero elements. For example, a transition from the state |1 to |0 corresponds to the operator |0 1| whose matrix M has only one non-zero element at the M [0, 1] location regardless of the total dimension of the matrix. Matrices with few non-zero elements multiplying each other will often produce the zero matrix or matrices with negligible norms. Hence by setting a threshold for the norm of the matrix, we can significantly reduce the number of terms to implement as in Equation (15) while maintaining a reasonable accuracy. The second way for term reduction is to notice that many of the K s terms for the s th time step are different by only a constant. This is again due to the fact that the M k 's typically contain few non-zero elements, such that repeated multiplications of certain M k 's tend to form closed groups. Hence by organizing the terms in Equation (15) into types within which the matrices only differ by a constant, we can also significantly reduce the number of terms to implement. The actual effects of both ways of term reduction are case specific, and we have shown the numbers in the FMO simulation section.

Simplified procedure for implementing products of Kraus operators
To implement each term in Equation (15) such as the M i M j M k for the 3 rd time step, we need to use a 3-dilation unitary having four times the dimension of M k , and this increases the gate count of the quantum circuit. Here expecting a large gate count that may exceed the capability of the current noisy intermediate-scaled quantum (NISQ) devices [48], we calculate the product of the matrices such as M i M j M k classically and only implement the product as a single 1-dilation unitary. Note that the sparsity of the M k matrices ensures the classical calculation of the product is simple, and with more powerful quantum devices in the future we can easily switch back to an all quantum implementation of the matrix products with higher dilations.

Procedure to obtain more data points
To obtain more data points such that the oscillations can be smoothly represented, we simulate four additional groups with six data points each: e.g. the first group starting with 400 atomic unit (9. where p 1 = γ 1 δt and p 2 = γ 2 δt are the probabilities of the corresponding transition happening over a small time interval δt. Now the population of the ground state after S iterations is: Recognizing that δt = t S , substituting in p 1 = γ 1 δt and p 2 = γ 2 δt, we have: which is the correct solution as in Equation S4. Therefore we see that the iterative procedure described in Equation (8) in the main text will give the correct solution if δt is small compared to the total time t.
S2 Parameters used for the FMO dynamics simulation in the main text.
The matrix form of the Hamiltonian of the FMO in the unit of eV:

S3
Number of projection measurements required to achieve a certain error limit.
A quantum state |φ = (c 1 , c 2 , ..., c n ) defines a probability distribution of the basis states (e.g. (1, 0, 0, .., 0) T , (0, 1, 0, .., 0) T , etc.) with the probability of getting the i th basis state equal to |c i | 2 . Applying projection measurements on such a quantum state is essentially sampling this probability distribution. It is a basic result that if we want to deduce the original probability distribution with some finite number of sampling measurements, then the error decreases with increasing number of measurements. In particular, if the error of the sampling can be represented by the standard error of the mean σ mean , and the original probability distribution defines a standard deviation σ, then: where P is the number of measurements. We then have the result mentioned in the following Section S4 that P does not scale with the dimension n of the quantum state, but only depends on the error σ mean we can tolerate.

S4 Quantum advantage in query complexity
Compared to the previous quantum algorithm, the generalized quantum algorithm generalizes to more complex dynamical models by utilizing an iterative process as in Equation (8) in the main text and formulating the Kraus operators with physical processes as in Equation (11)in the main text. Therefore, at the core of the generalized algorithm is still implementing matrices multiplying vectors by unitary dilations and an execution complexity comparable to classical methods is maintained at O(n 2 ) (n being the dimension of the vectors used to represent the initial states). However, here we emphasize that both versions of the quantum algorithm can have a significant query complexity advantage over classical methods, when used to evaluate an observable over the density matrix. As shown in Equations (5) and (6) in the main text, after the observable A has been converted intoÃ = L † L, we can evaluate A by calculating the trace Tr p i · L † |φ ik (t) φ ik (t)|L and then summing over i and k. The core idea here is classically taking the trace of the matrix of L † |φ ik (t) φ ik (t)|L involves summing over the modulus squares of all the coefficients of the vector L † |φ ik (t) , while quantum mechanically we can simply measure the probability of the final state U L † U M k (v T i , 0, ..., 0) T projecting into the first n-dimensional space. In other words, with the quantum method we do not have to look into individual coefficients of L † |φ ik (t) , but can treat the first n-dimensional space of U L † U M k (v T i , 0, ..., 0) T as a whole -this is an exploitation of quantum superposition and quantum measurement. Note to obtain the probability of projecting into a subspace, a number of quantum measurements are required. This number of measurements P however, does not scale with the dimension n but only depends on the error that we can tolerate, as explained in Section S3. Consequently the query complexity of the quantum method to determine Tr L † |φ ik (t) φ ik (t)|L is P that is a constant determined by the error tolerated, while the query complexity of any classical method cannot be reduced below O(n) (which is in general O(2 q ), exponential in the number of qubits q) for calculating and summing over the modulus square of each individual coefficient of L † |φ ik (t) . This means that the quantum method has a decisive exponential query complexity advantage over classical methods. Now if we define the total complexity of the algorithm to include the execution complexity and the query complexity, and then compare the quantum algorithm with classical methods, we see that applying P measurements requires P parallel implementations of the circuit for U L † U M k (v T i , 0, ..., 0) T (because multiple projection measurements cannot be applied to the same quantum state). This means that the total complexity of the quantum algorithm is the query complexity P multiplied by the execution complexity Q, while the total complexity of a classical algorithm with the same execution complexity is Q plus O(n). We see that if Q is O(n) or above, the quantum algorithm does not have a total complexity advantage over classical methods. However, when Q is O(log(n)) or below, which happens when the implementations of both U M k and U L † are very simple, the quantum algorithm will have a decisive advantage over classical methods. An illustrative example of this scenario can be constructed in the following. Suppose the initial density matrix is given by Equation (1) in the main text, where the physical composition of a few pure states is known. In fact, since the number of pure states in the mixture must be O(log(n)) or below for ρ to be representable with a cost that scales polynomially with the qubit number, we may assume there is only one pure state |φ = (c 1 , c 2 , ..., c n ) in the density matrix for simplicity. Now also suppose there is no dynamical evolution and we are interested in the expectation value of an observable on φ itself, then the execution cost for U M k is removed. Given that the system can be represented by q qubits such that 2 q = n, consider the observable: Now defineÃ = A+I||A|| 2||A|| , then the Cholesky decomposition gives: and Ã = Tr(L † |φ φ|L) = n 2 i=1 |c i | 2 (S12) where c i 's are the coefficients of the initial state |φ . The sum n 2 i=1 |c i | 2 in Equation (S12) has a clear physical meaning of the probability of projecting into the first half space where the first qubit has the value |0 , and thus can be measured by a projection measurement efficiently. Note that here the observable A is particularly simple that no actual implementation of L † is required. The total complexity of the quantum algorithm evaluating A over |φ is then only the query complexity of P .
On the other hand, for any classical method, evaluating the sum n 2 i=1 |c i | 2 costs O(n) steps and the total complexity is O(n) = O(2 q ). Consequently the quantum algorithm has a desicive total complexity advantage over any classical method for the example we have constructed. In addition to this special example, as discussed in the section on the generalized quantum algorithm in the main text, the Kraus operators are formulated with the basic physical processes that constitute the dynamical model, and thus typically contain only one or few non-zero elements. This means the execution cost for each M k is typically O(log(n)) or below, therefore if the number of M k 's is O(log(n)) or below, and implementing the observable A with L † also costs O(log(n)) or below, we will be able to demonstrate the quantum advantage on a dynamical evolution. To find such an observable A with a realistic physical meaning would be an interesting subject of a future study.
S5 An example of the quantum circuits used to evolve the Kraus operators Below we give an example of the quantum gate sequence used to build the quantum circuit for evolving M 1 = √ αδt |1 1| at the first time step at 400 atomic unit. After being multiplied by the unitary matrix accounting for the coherent part as in Equation (17)  (S14) Using the dilation procedure described around Equation 3 in the main text this becomes a 10 × 10 unitary matrix U M 1 . Using 4 qubits to cover the 10 dimensions and 2 ancilla qubits for the gate decomposition, we obtain a gate sequence of 899 gates as shown below. Note this sequence is not unique and may be further optimized, but it is good enough for the purpose of producing the results as shown in the main text when run on the IBM QASM simulator. The details of the decomposition of all the quantum circuits used are available from the corresponding author on reasonable request. Figure S1: An example of the quantum gate sequence used to evolve a Kraus operator