A Quantum Algorithm for Evolving

A quantum algorithm for evolving open quantum dynamics on quantum computing devices Designing quantum algorithms for simulating quantum systems has seen enormous progress, yet few studies have been done to develop quantum algorithms for open quantum dynamics despite its importance in modeling the system-environment interaction found in most realistic physical models. in this work we propose and demonstrate a general quantum algorithm to evolve open quantum dynamics on quantum computing devices. the Kraus operators governing the time evolution can be converted into unitary matrices with minimal dilation guaranteed by the Sz.-nagy theorem. this allows the evolution of the initial state through unitary quantum gates, while using significantly less resource than required by the conventional Stinespring dilation. We demonstrate the algorithm on an amplitude damping channel using the IBM Qiskit quantum simulator and the IBM Q 5 Tenerife quantum device. The proposed algorithm does not require particular models of dynamics or decomposition of the quantum channel, and thus can be easily generalized to other open quantum dynamical models. The time evolution of quantum systems is a century old subject that has been extensively studied for both fundamental and practical purposes. Open quantum dynamics is an important subfield of quantum physics that studies the time evolution of a system interacting with an environment1. Because the environment is usually too large to be treated exactly, in open quantum dynamics we often make approximations by averaging out the environment’s effect on the system. The resultant time evolution of the system density matrix is non-unitary and often governed by a master equation. The idea of simulating quantum systems with quantum algorithms was first proposed by Feynman2 and has received massive attention in recent years3–17 for its promise of outperforming the best available classical algorithms. However relatively few studies18–23 have been done to develop quantum algorithms for open quantum dynamics despite its importance. A key difficulty is the evolution of an open quantum system is often non-unitary, while quantum algorithms are mostly realized by unitary quantum gates. An early study18 tackled this problem by including the environment in the quantum simulation process therefore making the evolution unitary. Focused on Markovian dynamics, their procedure required a reset on the environment for every time step, which can become expensive if the system becomes entangled with the environment or the evolution time is large. It is known that any non-unitary quantum operation can be made into a unitary one by the Stinespring dilation theorem24. However, due to the large increase of the dimension of the Hilbert space, the computational cost required by naïve application of Stinespring dilation to the quantum operation can be prohibitive for actual implementation on a quantum computing device. Perhaps for this reason, most algorithms developed (see e.g. ref. 19,20) so far rely on the knowledge of the decomposition of the quantum channel which may not be generally available for a large system without costly computational efforts. So far as we know a general quantum algorithm to simulate an arbitrary quantum channel for a general density matrix with minimal resource has not been proposed and demonstrated. In this study we propose and demonstrate such a quantum algorithm utilizing the Sz.-Nagy dilation theorem – being a variation of the Stinespring dilation theorem the Sz.-Nagy dilation requires much smaller dimension increase and can save significant computational resources. Without assuming any specific property of the quantum channel, we work with the most general form of the time evolution for a density matrix – the operator sum representation: ρ(t ) = ∑ρk(t ) = ∑MkρM†k k k (1) 1 where ρ(t) is the system density matrix at time t, ρ is the initial system density matrix, and Mk’s are the Kraus operators that satisfy: ∑M†kMk = I (2) k It is also common that the time evolution of ρ(t) is described by a master equation: ρ(t ) = ρ(t ) (3) where the time derivative of ρ(t) is given by the superoperator  applying to ρ(t) itself. It is known that any master equation of the form of Eq. (3) can be converted into the form of Eq. (1) 25,26, affirming the generality of Eq. (1). An example of converting a widely used type of master equation – the Lindblad equation – into the operator sum representation has been provided in ref. 27. In this work we focus on the more general operator sum representation. Our method starts with an initial ρ expressed by a sum of different pure quantum states weighted by the associated probabilities: ρ= ∑pi φi i φi (4) where pi is the probability of finding the state |φi〉 in the mixture with the condition ∑ i pi = 1. Note here the different |φi〉’s are not necessarily orthogonal to each other, and Eq. (4) can be understood either as a knowledge of the initial physical composition of the system that is easily available from system preparation, or as a diagonal decomposition of ρ that requires additional resource to compute. We show in the following that the quantum algorithm proposed can simulate the time evolution of ρ(t). The outputs of the algorithm carry the full information of ρ(t) that can be extracted by quantum tomography28. However, we remark that the most important physical information carried by ρ(t)–e.g. the populations of different quantum states and the expectation value of an observable 〈O〉 = Tr(Oρ(t))–can be obtained without quantum tomography but by projection measurements instead, which greatly reduces the resources needed. This is realized by exploiting the positive-semidefiniteness of ρ(t). We then estimate the gate complexity of our method to be significantly lower than the conventional Stinespring dilation. Finally we demonstrate the application of the quantum algorithm to an amplitude damping channel with implementation on the IBM Qiskit29 quantum simulator and the IBM Q 5 Tenerife30 quantum device. Results Theory for the algorithm. In this section we present the quantum algorithm that evolves ρ(t) with the initial ρ given in the form ρ = ∑ i pi φi φi in Eq. (4), which represents a knowledge of the initial physical composition of the system. Here we prepare each |φi〉 as an input state vi in a given basis and want to evolve: φik(t ) = Mkvi = (cik1, cik2, … , cikn)T (5) First note that each Kraus operator Mk is a contraction such that it can be dilated into a unitary matrix. An operator A is a contraction if it shrinks or preserves the norm of any vector such that the operator norm Av † A = sup v ≤ 1. By Eq. (2) we have ∑kMkMk = I, which implies Mk is a contraction (see the supplementary information (SI) for the proof). By the Sz.-Nagy dilation theorem31,32, any contraction A of a Hilbert space H has a corresponding unitary operator UA in a larger Hilbert space K such that: An = PH U nAPH , n≤N (6) where PH is the projection operator into space H. The physical meaning of Eq. (6) is that the effect of a contraction A applied up to N times on a smaller space H can be replicated by a unitary UA applied up to N times on a larger space K, given the input vector lies entirely in H and the output vector is projected into H. For the purpose of creating a quantum circuit the Sz.-Nagy dilation theorem allows us to simulate the effect of any non-unitary matrix by a unitary quantum gate, because every operator on a finite dimensional space is bounded and therefore can be made into a contraction which has a unitary dilation. The Sz.-Nagy theorem also guarantees the existence of a minimal dilation in the sense that the space K has the smallest dimension to achieve Eq. (6). An example of a minimal unitary dilation of A with N = 1 is:  A DA†   UA =  †  DA −A  (7) † where DA = I − A A is called the defect operator of A. We can easily verify that UA is unitary and A = PH UAPH . However if we want to apply A2 = PH UA2 PH , or BA = PH UBUAPH where UB is a unitary dilation of B, then we need minimal unitary dilations with N = 2:  A 0 DA†    UA =  D 0 −A†  , UB = A    0 I 0   B 0 DB†     D 0 −B†   B    0 I 0  (8) 2 We see that the number N in Eq. (6) is an important parameter that defines both the form and the applicability of the minimal unitary dilation. In the following we will refer to a minimal unitary dilation with a given N an N-dilation31. A general rule is that to simulate the effect of N contractions multiplying successively, we need to convert all of them to N-dilations. Going back to our original goal of simulating Mkvi with unitary gates we construct a 1-dilation UA of the form in Eq. (7) with A = Mk and evolve: unitary dilation φik(t ) = Mkvi    → UMk(vTi , 0, . . . , 0)T (9) If Mk has the dimension n by n, then the 1-dilation UMk is 2n by 2n. Now UMk can be further decomposed into sequences of two-level unitary gates with a procedure illustrated in ref. 3,33. The two-level unitary gates can be used as elementary gates in an optical beamsplitter setup33–36 and in the following we use the number of two-level unitary gates in the decomposition of a unitary gate to represent the gate complexity. Generally the number of two-level unitary gates needed to decompose a unitary gate is equal to the number of non-zero elements in the 2 lower-triangular part of the gate3,33. For each UMk of the form in Eq. (7) with A = Mk, we have n − n non-zero 2 2 elements from Mk, n2 non-zero elements from DMk, n − n non-zero elements from −M†k, and the total gate count 2 is 2n2 − n for each k and i. The classical complexity of evaluating φik(t ) = Mkvi is incidentally also 2n2 − n for each k and i, counting all the multiplications and additions needed for a matrix multiplying a vector. Now we have the gate count for a single Kraus operator on a single pure state Mkvi , the maximum number3 of Kraus operators in Eq. (1) is m = n2, and if the number of pure states in the physical composition Eq. (4) is comparable to n, we have totally 2n5 − n4 two-level 2n by 2n gates. As a comparison, the conventional Stinespring dilation converts the whole quantum operation (ρ ) = ∑kMkρ M†k into a unitary operation U(ρ ⊗ e e )U † by considering the tensor product space of the initial ρ and an auxiliary environment3. By results from ref. 3, the dimension of U is given by n ⋅ m where m is the total number of Kraus operators Mk in the sum ∑kMkρ M†k. In the most general case 6 3 m = n2, U is n3 by n3 which requires n − n two-level n3 by n3 unitary gates for the gate decomposition3. Therefore 2 comparing the total cost for the entire quantum operation our method has a moderate advantage (O(n5) vs. O(n6)) over the conventional Stinespring method. In addition, the complexity for realizing a 2n by 2n two-level unitary is significantly smaller than that for realizing a n3 by n3 two-level unitary on conventional quantum computers using 1-qubit and 2-qubit elementary gates, due to the complication involved in multi-qubit controlled gates3. Finally, the greatest advantage of our method is the ability to break the whole quantum operation into n3 processes (each Mkvi costing 2n2 − n two-level unitaries) that can be simulated in parallel. Because all current quantum computing platforms have a limit on the maximum circuit depth due to gate errors, our method can enable the simulation of more complicated open quantum systems by replacing a single O(n6) process with n3 parallel O(n2) processes. Now to put the pieces back together, the full evolved density matrix can be assembled by ρ(t ) = ∑ ikpi ⋅ φik(t ) φik(t ) and measured by quantum tomography28. However, we remark that the most important physical information such as the populations of different states and the expectation value of an observable can be extracted by performing projection measurements on the output φik(t ) without actually determining the state of φik(t ) .This means that we do not need the phases of the coefficients associated with each basis state, but only the absolute values, thus saving considerable computational costs. To get the populations of states in the current basis, note the diagonal vector of φik(t ) φik(t ) is: diag ( φik(t ) φik(t ) ) = ( cik1 2 , cik2 2 , … , cikn 2 )T (10) 2 Equation (10) means we can obtain the diagonal element cikj of φik(t ) φik(t ) by applying a projection measurement on the jth entry in the first n-dimensional subspace of UMk(vTi , 0, … , 0)T . Using an optical setup such as in ref. 36 the probability of measuring each entry in UMk(vTi , 0, … , 0)T can be efficiently obtained by recording the photon distribution at the output of the optical modes. Adding the results through k and i gives the diagonal elements of ρ(t ): diag (ρ(t )) = ∑pi ik ⋅ diag ( φik(t ) φik(t ) ) (11) which gives the populations in the current basis. Although the off-diagonal elements of ρ(t ) cannot be directly obtained without quantum tomography, they are nonetheless carried by φik(t ) and can become physically important. For example, if we want to get the populations in another basis, a unitary basis transformation can be applied to each φik(t ) before measuring the diagonal elements: ( ) unitary dilation T φik(t ) = TMkvi    → T 0 UMk(vTi , 0, . . . , 0)T 0 I diag (Tρ(t )T †) = ∑pi ⋅ diag (T φik(t ) φik(t ) T †) ik (12) 2 The additional unitary T applied to φik(t ) requires no dilation and increases the quantum gate count by n − n 2 2 (non-zero elements in the lower-triangular part of T) to a total of 5n − 3n for each k and i. The classical complexity 2 of evaluating TMkvi is doubled from 2n2 − n to 4n2 − 2n for each k and i. Note here the quantum algorithm outperforms the classical one by taking advantage of the unitarity of T. Next to evaluate the expectation value of an observable O = Tr(Oρ(t )), we recognize Oρ(t ) is not always positive-semidefinite and hence additional processing must be done before the trace can be obtained from projection measurements on the output quantum state. The idea here is to see the operator norm is bounded by the Hilbert-Schmidt norm: O ≤ O HS, such that we define: O + I O HS ∼ O= 2 O HS (13) where Õ is always a contraction and positive-semidefinite (see the SI for a detailed proof) so we can apply the ∼ ∼ Cholesky decomposition37: O = LL†, where L is a lower triangular matrix. Note we could have defined O with O in place of O HS in Eq. (13), but O HS requires fewer arithmetic steps to calculate. Now we have: ∼ ∼ O = Tr(Oρ(t )) = Tr(LL†ρ(t )) = Tr(L†ρ(t )L) = ∑Tr(L†ρk(t )L) (14) k † where each L†ρk(t )L is positive-semidefinite. L is obviously a contraction because Õ is a contraction. Now evolve L† φik(t ) with two 2-dilations of the form in Eq. (8) with A = Mk and B = L†: unitary dilation L† φik(t ) = L†Mkvi    → UL†UMk(vTi , 0, . . . , 0)T (15) and we can evaluate 〈Õ〉 by: ∼ ∼ O = Tr(Oρ(t )) = ∑Tr(pi i ,k ⋅ L† φik(t ) φik(t ) L) (16) where the trace of L† φik(t ) φik(t ) L can be obtained by projection measurements in a way similar to what has been explained for Eq. (10). The difference here is that we do not need to measure individual diagonal elements and sum them together for the trace, but instead can measure the total probability of projecting into the first n T -dimensional space of UL†UMk(vTi , 0, … , 0) . This can potentially save measurement costs by reducing the number of inquiries needed on the output vector. Finally 〈O〉 is calculated with: O = Tr(Oρ(t )) = Tr((2 O ∼ HS O −I O HS )ρ(t )) =2 O ∼ HS Tr (Oρ(t )) − O HS (17) where we have successfully obtained O for the original observable. The L† gate (L† is upper triangular requiring reduced number of two-level unitaries for the decomposition) plus 2 the additional level of dilation for Mk increases the quantum gate count to 5n + 3n for each k and i. A classical 23 overhead (independent from the k and i counts) cost of 2n2 − 1 for O HS and n for the Cholesky decomposi3 tion37 should also be counted towards the total cost of the quantum algorithm. On the other hand the classical 2 complexity of evaluating L†Mkvi is 3n − n for each k and i (taking into account that L† is upper triangular), plus 3 an overhead of n + 2n2 − 1. 3 Complexity analysis. So far in evaluating the complexity of our method, each Kraus operator Mk is assumed to be completely general, and thus the number of 2-level unitary gates required to realize UMk is polynomial in the system dimension n (O(n2)). If using a conventional qubit-based hardware, q qubits can represent a n = 2q -dimensional system, and the number of 2-level unitaries is exponential in the qubit number q (O(4q)). We note that this exponential scaling in the qubit number is necessary if Mk is a general matrix without any special properties, because approximating arbitrary unitary gates is generically exponential3. That said, in the case when Mk is not a full matrix with arbitrary entries (Mk is sparse), the algorithm can become polynomial in q. This happens in many realistic settings (like the amplitude damping channel shown below) when each Mk only mediates the transition process between l quantum states, where l < n is a constant of n and q. This effectively reduces the dimension of the system from n to l for each Mk and the number of 2-level unitaries becomes O(l2), without a dependence on n and q. Note that we still need to run all the m number of Mk’s through all the initial composite states φi ’s, but this can be done in parallel and the complexity of the circuit for each Mk is greatly reduced as compared to the conventional Stinespring dilation. If we stop the analysis on the level of 2-level unitary gates the complexity for a sparse Mk is constant in n and q. This desirable result is possible if using an optical beamsplitter setup33–36 to implement the 2-level unitaries as elementary gates. If a conventional qubit-based quantum computer is used, the 2-level unitries need to be further decomposed into 1-qubit and 2-qubit elementary gates, in which case a Gray code sequence and a multi-control gate sequence are constructed with the complexity logarithmic in n (O(log2 n)) and polynomial in q (O(q2))3. When implemented on a qubit-based quantum computer, our method also requires ancillary qubits for the Sz.-Nagy dilation process that enlarges the Hilbert space. For the simple evolution described in Eqs. (9) and (12) we need to increase the dimension by two-fold and one ancillary qubit is needed. For the observable evaluation described in Eq. (15) we need to increase the dimension by three-fold and two ancillary qubits are needed. Note that our use of the ancillary qubits is very simple as they do not scale with the system size or the complexity of the quantum operation. 4 Figure 1. Showing the populations of the ground and excited states for the amplitude damping model. The smooth lines are obtained by classical numerical calculations of the output vectors. These agree exactly with analytic results and are used as benchmarks. The crosses are obtained by the IBM Qiskit simulator. The dots are obtained by the IBM Q 5 Tenerife device. The quantum circuits include 2 qubits and on average 13 elementary gates (see Fig. 4 for an example and the SI for all the circuits). Application to the amplitude damping channel. In this section we use a quantum channel model to demonstrate the application of the proposed method. The amplitude-damping channel models the spontaneous emission of a 2-level atom. The Lindblad master equation for this model is:  1  ρ(t ) = γσ +ρ(t )σ− − {σ−σ +, ρ(t )}   2 (18) where γ is the spontaneous emission rate, σ + = 0 1 is the Pauli raising operator that mediates the transition from the excited state to the ground state, and σ− = (σ +)†. In the operator sum representation: ρ(t ) = M 0ρ M †0 + M1ρ M1† M0 = M1 = e−γt e−γt 0   e−γt  1 σz =  0 2 2  −γt  1 − e−γt σ + = 0 1 − e  0  0 1+ I+ 1− (19) To calculate the populations in the initial basis { 0 , 1 } we can construct UM 0 and UM1 of the 1-dilation form in Eq. (7) with A = Mk and DMk = I − M†kMk : 0 DM 0 =   0 Using an arbitrary initial ρ = known, the input states are: 1 4 0    , D = 1 M1  0  −γt   1−e 0   e−γt  (11 13) and assuming the physical composition ρ = (20) 1 (|1〉〈1|+ 2 T | +〉〈+|) is T m m        1  φ1 = 1 → v1 = 0, 1, 0, . . . , 0  φ2 = + → v2 = 1, 1, 0, . . . , 0    2   (21) where m = 2 matching the size of the vectors with the dilation UMk. We set γ = 1.52 × 109 s−1 (typical nanosecond lifetime), numerically calculate UMkvi from t = 0 to t = 1000 ps with a time step of 10 picosecond, and obtain the populations of the ground and excited states from the first two entries of UMkvi . The results are shown as the smooth lines in Fig. 1. To calculate the populations in another basis { + , − } where |± = 1 ( 0 ± 1 ), we need the transforma2 tion matrix T = 1 1 1 . We numerically calculate T 0 UM vi and obtain the populations of the ± states k 2 1 −1 0 I shown as the smooth lines in Fig. 2. Now we evaluate the expectation value of an observable O for O = −2 0.5 as an example. With O HS = 0. 5 1 O + I‖O‖HS 22 ~ 0.0740 0.107 and find L ≈ 0.271 0 , we define O through Cholesky ≈ 2 . 35 = ≈ 2 2‖O‖HS 0.107 0.713 0.393 0.748 ∼ † decomposition O = LL . Next we construct UL† and UMk of the 2-dilation form in Eq. (8) and apply them to the ( ( ) ) ( ) ( ( ) ) 5 Figure 2. Showing the populations of the + and − states for the amplitude damping model. The smooth lines are obtained by classical numerical calculations of the output vectors. These agree exactly with analytic results and are used as benchmarks. The crosses are obtained by the IBM Qiskit simulator. The dots are obtained by the IBM Q 5 Tenerife device. The quantum circuits include 2 qubits and on average 30 elementary gates (see the SI for the circuits). Figure 3. Showing the expectation values O . The smooth line is obtained by classical numerical calculations of the output vectors. This agrees exactly with analytic results and is used as a benchmark. The crosses are obtained by the IBM Qiskit simulator. The quantum circuits include 3 qubits and on average 184 elementary gates (see the SI for the circuits). Due to the large number of gates required the quantum device is not used for these results. initial state vi in the form of Eq. (21) with m = 4. Numerically calculating the output vector will give us O by Eqs. (16) and (17). The results are shown as the smooth line in Fig. 3: The overall simplicity of the unitary dilation gates used in our algorithm allows us to further demonstrate it using the IBM Qiskit29 quantum simulator and the IBM Q 5 Tenerife30 quantum device. So far we have been using the 2-level unitaries obtained from the UM 0 and UM1 gates as elementary gates in our complexity evaluation. This is possible if we use an optical platform for quantum computation involving beamsplitters33,35. To implement the algorithm on a conventional quantum platform such as the IBM simulator and devices, we further decompose the 2-level unitaries into elementary 1-qubit and 2-qubit gates. The constructed circuits require 2 qubits and on average 13 elementary gates for the basic UMkvi evolution and 30 elementary gates for the basis transformation T 0 U v . An example of the circuit for U v is shown in Fig. 4, and a full list of the circuits can be found in M0 1 0 I Mk i the SI. We implement these circuits on both the simulator and the quantum devices and show the results as the ∼ crosses (simulator) and dots (device) in Figs. 1 and 2. To obtain O we construct a circuit of 3 qubits with an average of 182 elementary gates for the UL†UMkvi operation (see the SI). The increase in the number of gates is due to the increased size of the UL† and UMk matrices and the associated 2-level unitaries. As mentioned earlier the decomposition of 2-level unitaries into 1-qubit and 2-qubit elementary gates scales as O(log2 n). The large number of gates for this circuit prevents us from running it on the quantum device so that only the simulator results are shown as the crosses in Fig. 3. On the other hand the 2-level unitaries can be more easily implemented with a ( ) 6 Figure 4. Showing the quantum circuit for UM 0v1 as used in the IBM simulator and quantum device. For a full list of the quantum circuits please see the SI. multiport photonic device as demonstrated in refs. 35,36. We will seek to demonstrate our methods on such a preferred device in a future study. In Figs. 1, 2 and 3, the numerical results (smooth lines) agree exactly with analytic results and are used as benchmarks for the simulator and device results. The simulator results (crosses) fit the benchmark well while including the probabilistic error from the projection measurements. The quantum device results (dots) fit the benchmark reasonably well considering all experimental errors from gate fault, qubit decoherence, and measurement. These results demonstrate the ability of the proposed algorithm to simulate the time evolution of a density matrix and evaluate physical observables from the outputs. Discussion In this work we have presented a quantum algorithm for evolving the density matrix ρ(t ) = ∑kMkρ M†k and extracting physical information from the output. The method takes each physical composition φi of as the input and evolves it through minimal Sz.-Nagy dilations of the Kraus operators Mk. The ρ = ∑ i pi φi φi input vector has the base length n (plus additional zeros matching the dimension of the evolution matrices) and various realizations of the time evolution have the quantum gate count of O(n2) for each k and i, which is a significant improvement over the O(n6) scaling of the conventional Stinespring dilation. In cases when Mk involves a reduced number of quantum states in the system (Mk is sparse) the complexity scaling can become O(log2 n), or polynomial in the number of qubits. The requirement of knowing the initial physical composition as a probability-weighted mixture of not necessarily orthogonal φi ’s should be easy to satisfy from system preparation. If indeed the initial physical composition is unknown and we have to work with the most general matrix form of the density matrix, then we can either perform the eigen-decomposition on ρ to obtain the form in Eq. (4) or use a modified quantum algorithm. As details shown in the SI, the modified method uses the flattened vector of the initial density matrix as the input and requires O(n3) gates per Kraus operator for various realizations of the time evolution. Both methods can be easily generalized to other open quantum dynamical models because the procedures involved are essentially the same – only the Kraus operators Mk’s for the operator sum representation are different for different models. The generality of our methods–in the sense of not requiring particular dynamical models or costly decompositions of the density matrix and the quantum channel–opens up the possibility of simulating more interesting systems such as decohering qubits or excitonic structures interacting with multiple baths38,39 – the latter helps to understand natural light harvesting complexes and exploit quantum coherence effect to improve light harvesting efficiency in artificial photocells38,40. Finally we have demonstrated the implementation of the algorithm on the IBM Qiskit simulator and the IBM Q 5 Tenerife device. Although the gate complexity is larger than calculated with the preferred optical setup, the results show reasonable agreements with the analytic benchmarks considering gate fault, qubit decoherence, and measurement error. In future studies we will seek to demonstrate our quantum algorithms on the preferred photonic devices that can implement 2-level unitaries as elementary gates. Methods Theoretical derivations and numerical details essential to the study are presented in the Results section. More details can be found in the Supplementary Information. For the parameters and configurations of the IBM simulator and quantum machine, please go to IBM Quantum Experience at http://www.research.ibm.com/quantum. References 1. Breuer, H.-P. & Petruccione, F. The Theory of Open Quantum Systems. (Oxford University Press, 2007). 2. Feynman, R. P. Simulating physics with computers. International Journal of Theoretical Physics 21, 467–488, https://doi.org/10.1007/ BF02650179 (1982). 3. Nielsen, M. A. & Chuang, I. L. Quantum Computation and Quantum Information: 10th Anniversary Edition. (Cambridge University Press, 2011). 4. Lloyd, S. Universal Quantum Simulators. Science 273, 1073, https://doi.org/10.1126/science.273.5278.1073 (1996). 5. Kais, S. In Quantum Information and Computation for Chemistry (John Wiley & Sons, 2014). 6. Kassal, I., Whitfield, J. D., Perdomo-Ortiz, A., Yung, M.-H. & Aspuru-Guzik, A. Simulating Chemistry Using Quantum Computers. Annual Review of Physical Chemistry 62, 185–207, https://doi.org/10.1146/annurev-physchem-032210-103512 (2011). 7. Georgescu, I. M., Ashhab, S. & Nori, F. Quantum simulation. Reviews of Modern Physics 86, 153–185, https://doi.org/10.1103/ RevModPhys.86.153 (2014). 8. Abrams, D. S. & Lloyd, S. Simulation of Many-Body Fermi Systems on a Universal Quantum Computer. Physical Review Letters 79, 2586–2589, https://doi.org/10.1103/PhysRevLett.79.2586 (1997). 7 9. King, A. D. et al. Observation of topological phenomena in a programmable lattice of 1,800 qubits. Nature 560, 456–460, https://doi. org/10.1038/s41586-018-0410-x (2018). 10. Bravyi, S., DiVincenzo, D. P., Loss, D. & Terhal, B. M. Quantum Simulation of Many-Body Hamiltonians Using Perturbation Theory with Bounded-Strength Interactions. Physical Review Letters 101, 070503, https://doi.org/10.1103/PhysRevLett.101.070503 (2008). 11. Wu, L. A., Byrd, M. S. & Lidar, D. A. Polynomial-Time Simulation of Pairing Models on a Quantum Computer. Physical Review Letters 89, 057904, https://doi.org/10.1103/PhysRevLett.89.057904 (2002). 12. Babbush, R., Love, P. J. & Aspuru-Guzik, A. Adiabatic Quantum Simulation of Quantum Chemistry. Scientific Reports 4, 6603, https://doi.org/10.1038/srep06603 (2014). 13. O’Malley, P. J. J. et al. Scalable Quantum Simulation of Molecular Energies. Physical Review X 6, 031007, https://doi.org/10.1103/ PhysRevX.6.031007 (2016). 14. Young, K. C., Sarovar, M., Aytac, J., Herdman, C. M. & Whaley, K. B. Finite temperature quantum simulation of stabilizer Hamiltonians. Journal of Physics B: Atomic, Molecular and Optical Physics 45, 154012, https://doi.org/10.1088/09534075/45/15/154012 (2012). 15. Xia, R. & Kais, S. Quantum machine learning for electronic structure calculations. Nature Communications 9, 4195, https://doi. org/10.1038/s41467-018-06598-z (2018). 16. Xia, R., Bian, T. & Kais, S. Electronic Structure Calculations and the Ising Hamiltonian. The Journal of Physical Chemistry B 122, 3384–3395, https://doi.org/10.1021/acs.jpcb.7b10371 (2018). 17. Karra, M., Sharma, K., Friedrich, B., Kais, S. & Herschbach, D. Prospects for quantum computing with an array of ultracold polar paramagnetic molecules. The Journal of Chemical Physics 144, 094301, https://doi.org/10.1063/1.4942928 (2016). 18. Wang, H., Ashhab, S. & Nori, F. Quantum algorithm for simulating the dynamics of an open quantum system. Physical Review A 83, 062317 (2011). 19. Wang, D.-S., Berry, D. W., de Oliveira, M. C. & Sanders, B. C. Solovay-Kitaev Decomposition Strategy for Single-Qubit Channels. Physical Review Letters 111, 130504, https://doi.org/10.1103/PhysRevLett.111.130504 (2013). 20. Wei, S.-J., Ruan, D. & Long, G.-L. Duality quantum algorithm efficiently simulates open quantum systems. Scientific Reports 6, 30727, https://doi.org/10.1038/srep30727 (2016). 21. Di Candia, R., Pedernales, J. S., del Campo, A., Solano, E. & Casanova, J. Quantum Simulation of Dissipative Processes without Reservoir Engineering. Scientific Reports 5, 9981, https://doi.org/10.1038/srep09981 https://www.nature.com/articles/ srep09981#supplementary-information (2015). 22. Sweke, R., Sinayskiy, I., Bernard, D. & Petruccione, F. Universal simulation of Markovian open quantum systems. Physical Review A 91, 062308, https://doi.org/10.1103/PhysRevA.91.062308 (2015). 23. Wei, S.-J., Xin, T. & Long, G.-L. Efficient universal quantum channel simulation in IBM’s cloud quantum computer. Science China Physics, Mechanics & Astronomy 61, 70311, https://doi.org/10.1007/s11433-017-9181-9 (2018). 24. Buscemi, F., D’Ariano, G. M. & Sacchi, M. F. Physical realizations of quantum operations. Physical Review A 68, 042113, https://doi. org/10.1103/PhysRevA.68.042113 (2003). 25. Tong, D. M., Kwek, L. C., Oh, C. H., Chen, J.-L. & Ma, L. Operator-sum representation of time-dependent density operators and its applications. Physical Review A 69, 054102, https://doi.org/10.1103/PhysRevA.69.054102 (2004). 26. Lidar, D. A., Bihary, Z. & Whaley, K. B. From completely positive maps to the quantum Markovian semigroup master equation. Chemical Physics 268, 35–53, https://doi.org/10.1016/S0301-0104(01)00330-5 (2001). 27. Nakazato, H. et al. Solution of the Lindblad equation in the Kraus representation. Physical Review A 74, 062113, https://doi. org/10.1103/PhysRevA.74.062113 (2006). 28. Lvovsky, A. I. & Raymer, M. G. Continuous-variable optical quantum-state tomography. Reviews of Modern Physics 81, 299–332, https://doi.org/10.1103/RevModPhys.81.299 (2009). 29. Aleksandrowicz, G. et al. Qiskit: An Open-source Framework for Quantum Computing. IBM Qiskit (2019). 30. IBM Quantum Experience http://www.research.ibm.com/quantum. 31. Levy, E. & Shalit, O. Dilation theory in finite dimensions: The possible, the impossible and the unknown. Rocky Mountain Journal of Mathematics 44, 203–221 (2014). 32. Langer, H., Nagy-Sz, B. & Foias, C. Harmonic Analysis of Operators on Hilbert Space. VIII+387S. Budapest/Amsterdam/London 1970. Akadémiai Kiadó/North-Holland Publishing Company. ZAMM - Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik 52, 501–501, https://doi.org/10.1002/zamm.19720520821 (1972). 33. Reck, M., Zeilinger, A., Bernstein, H. J. & Bertani, P. Experimental realization of any discrete unitary operator. Physical Review Letters 73, 58–61, https://doi.org/10.1103/PhysRevLett.73.58 (1994). 34. Clements, W. R., Humphreys, P. C., Metcalf, B. J., Kolthammer, W. S. & Walmsley, I. A. Optimal design for universal multiport interferometers. Optica 3, 1460–1465, https://doi.org/10.1364/OPTICA.3.001460 (2016). 35. Carolan, J. et al. Universal linear optics. Science 349, 711, https://doi.org/10.1126/science.aab3642 (2015). 36. Sparrow, C. et al. Simulating the vibrational quantum dynamics of molecules using photonics. Nature 557, 660–667, https://doi. org/10.1038/s41586-018-0152-9 (2018). 37. Krishnamoorthy, A. & Menon, D. In Signal Processing: Algorithms, Architectures, Arrangements, and Applications (SPA). 70–72 (2013). 38. Hu, Z., Engel, G. S. & Kais, S. Double-excitation manifold’s effect on exciton transfer dynamics and the efficiency of coherent light harvestin. g. Physical Chemistry Chemical Physics 20, 30032–30040, https://doi.org/10.1039/C8CP05535A (2018). 39. Yeh, S.-H., Hoehn, R. D., Allodi, M. A., Engel, G. S. & Kais, S. Elucidation of near-resonance vibronic coherence lifetimes by nonadiabatic electronic-vibrational state character mixing. Proceedings of the National Academy of Sciences, 201701390, https://doi. org/10.1073/pnas.1701390115 (2018). 40. Hu, Z., Engel, G. S., Alharbi, F. H. & Kais, S. Dark states and delocalization: Competing effects of quantum coherence on the efficiency of light harvesting systems. The Journal of Chemical Physics 148, 064304, https://doi.org/10.1063/1.5009903 (2018). competing interests The authors declare no competing interests.