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.

Introduction

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 environment[1]. 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 Feynman[2] and has received massive attention in recent years[3-17] for its promise of outperforming the best available classical algorithms. However relatively few studies[18-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 study[18] 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 theorem[24]. 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:where ρ(t) is the system density matrix at time t, ρ is the initial system density matrix, and M’s are the Kraus operators that satisfy:

It is also common that the time evolution of ρ(t) is described by a master equation: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:where p is the probability of finding the state |ϕ〉 in the mixture with the condition . Note here the different |ϕ〉’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 tomography[28]. 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 Qiskit[29] quantum simulator and the IBM Q 5 Tenerife[30] 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 in Eq. (4), which represents a knowledge of the initial physical composition of the system. Here we prepare each |ϕ〉 as an input state v in a given basis and want to evolve:

First note that each Kraus operator M 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 . By Eq. (2) we have , which implies M is a contraction (see the supplementary information (SI) for the proof). By the Sz.-Nagy dilation theorem[31,32], any contraction A of a Hilbert space H has a corresponding unitary operator UA in a larger Hilbert space K such that:where P 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 U 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 with N = 1 is:where is called the defect operator of A. We can easily verify that UA is unitary and . However if we want to apply , or where UB is a unitary dilation of B, then we need minimal unitary dilations with N = 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-dilation[31]. 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 with unitary gates we construct a 1-dilation of the form in Eq. (7) with and evolve:

If has the dimension by , then the 1-dilation is by . Now 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 setup[33-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 lower-triangular part of the gate[3,33]. For each of the form in Eq. (7) with , we have non-zero elements from , non-zero elements from , non-zero elements from , and the total gate count is for each and . The classical complexity of evaluating is incidentally also for each and , 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 , the maximum number[3] of Kraus operators in Eq. (1) is , and if the number of pure states in the physical composition Eq. (4) is comparable to , we have totally two-level by gates. As a comparison, the conventional Stinespring dilation converts the whole quantum operation into a unitary operation by considering the tensor product space of the initial and an auxiliary environment[3]. By results from ref. [3], the dimension of is given by where is the total number of Kraus operators in the sum . In the most general case , is by which requires two-level by unitary gates for the gate decomposition[3]. Therefore 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 by two-level unitary is significantly smaller than that for realizing a by two-level unitary on conventional quantum computers using 1-qubit and 2-qubit elementary gates, due to the complication involved in multi-qubit controlled gates[3]. Finally, the greatest advantage of our method is the ability to break the whole quantum operation into processes (each costing 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 and measured by quantum tomography[28]. 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 without actually determining the state of .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 is:

Equation (10) means we can obtain the diagonal element of by applying a projection measurement on the entry in the first -dimensional subspace of . Using an optical setup such as in ref. [36] the probability of measuring each entry in can be efficiently obtained by recording the photon distribution at the output of the optical modes. Adding the results through and gives the diagonal elements of :which gives the populations in the current basis. Although the off-diagonal elements of cannot be directly obtained without quantum tomography, they are nonetheless carried by 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 before measuring the diagonal elements:

The additional unitary applied to requires no dilation and increases the quantum gate count by (non-zero elements in the lower-triangular part of ) to a total of for each and . The classical complexity of evaluating is doubled from to for each and . Note here the quantum algorithm outperforms the classical one by taking advantage of the unitarity of .

Next to evaluate the expectation value of an observable , we recognize 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: , such that we define:where Õ is always a contraction and positive-semidefinite (see the SI for a detailed proof) so we can apply the Cholesky decomposition[37]: , where is a lower triangular matrix. Note we could have defined with in place of in Eq. (13), but requires fewer arithmetic steps to calculate. Now we have:where each is positive-semidefinite. is obviously a contraction because Õ is a contraction. Now evolve with two 2-dilations of the form in Eq. (8) with and :and we can evaluate 〈Õ〉 by:where the trace of 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 -dimensional space of . This can potentially save measurement costs by reducing the number of inquiries needed on the output vector. Finally 〈O〉 is calculated with:where we have successfully obtained for the original observable.

The gate ( is upper triangular requiring reduced number of two-level unitaries for the decomposition) plus the additional level of dilation for increases the quantum gate count to for each and . A classical overhead (independent from the and counts) cost of for and for the Cholesky decomposition[37] should also be counted towards the total cost of the quantum algorithm. On the other hand the classical complexity of evaluating is for each and (taking into account that is upper triangular), plus an overhead of .

Complexity analysis

So far in evaluating the complexity of our method, each Kraus operator is assumed to be completely general, and thus the number of 2-level unitary gates required to realize is polynomial in the system dimension n (O(n2)). If using a conventional qubit-based hardware, qubits can represent a -dimensional system, and the number of 2-level unitaries is exponential in the qubit number (O(4)). We note that this exponential scaling in the qubit number is necessary if is a general matrix without any special properties, because approximating arbitrary unitary gates is generically exponential[3]. That said, in the case when is not a full matrix with arbitrary entries ( is sparse), the algorithm can become polynomial in . This happens in many realistic settings (like the amplitude damping channel shown below) when each only mediates the transition process between quantum states, where is a constant of and . This effectively reduces the dimension of the system from to for each and the number of 2-level unitaries becomes O(l2), without a dependence on and . Note that we still need to run all the number of ’s through all the initial composite states ’s, but this can be done in parallel and the complexity of the circuit for each 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 is constant in and . This desirable result is possible if using an optical beamsplitter setup[33-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 (O(log2 n)) and polynomial in (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.

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:where is the spontaneous emission rate, is the Pauli raising operator that mediates the transition from the excited state to the ground state, and . In the operator sum representation:

To calculate the populations in the initial basis we can construct and of the 1-dilation form in Eq. (7) with and :

Using an arbitrary initial and assuming the physical composition is known, the input states are:where m = 2 matching the size of the vectors with the dilation . We set (typical nanosecond lifetime), numerically calculate from 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 . The results are shown as the smooth lines in Fig. 1.

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).

Figure 4

Showing the quantum circuit for as used in the IBM simulator and quantum device. For a full list of the quantum circuits please see the SI.

To calculate the populations in another basis where , we need the transformation matrix . We numerically calculate and obtain the populations of the states shown as the smooth lines in Fig. 2.

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).

Now we evaluate the expectation value of an observable for as an example. With , we define and find through Cholesky decomposition . Next we construct and of the 2-dilation form in Eq. (8) and apply them to the initial state in the form of Eq. (21) with . Numerically calculating the output vector will give us by Eqs. (16) and (17). The results are shown as the smooth line in Fig. 3:

Figure 3

Showing the expectation values . 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.

The overall simplicity of the unitary dilation gates used in our algorithm allows us to further demonstrate it using the IBM Qiskit[29] quantum simulator and the IBM Q 5 Tenerife[30] quantum device. So far we have been using the 2-level unitaries obtained from the and gates as elementary gates in our complexity evaluation. This is possible if we use an optical platform for quantum computation involving beamsplitters[33,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 evolution and 30 elementary gates for the basis transformation . An example of the circuit for is shown in Fig. 4, and a full list of the circuits can be found in 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 we construct a circuit of 3 qubits with an average of 182 elementary gates for the operation (see the SI). The increase in the number of gates is due to the increased size of the and 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 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 and extracting physical information from the output. The method takes each physical composition of as the input and evolves it through minimal Sz.-Nagy dilations of the Kraus operators . The input vector has the base length (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 and , which is a significant improvement over the O(n6) scaling of the conventional Stinespring dilation. In cases when involves a reduced number of quantum states in the system ( 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 ’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 ’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 baths[38,39] – the latter helps to understand natural light harvesting complexes and exploit quantum coherence effect to improve light harvesting efficiency in artificial photocells[38,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.

Supplementary information.