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.

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.arXiv: 1904.00910[quant-ph] 14 May 2019

I. 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][4][5][6][7][8][9][10][11][12][13][14][15][16][17] for its promise of outperforming the best available classical algorithms.However relatively few studies [18][19][20][21][22] 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 [23].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 k 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 L 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) [24,25], 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. [26].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 i p is the probability of finding the state i  in the mixture with the condition Note here the different i  's are not necessarily orthogonal to each other, and Eq. ( 4) should be understood as a knowledge of the initial physical composition of the system that is easily available from system preparation, not as a diagonal decomposition of  that requires considerable 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 [27].However, we remark that the most important physical  .We then estimate the gate complexity of our method to be significantly lower than the conventional Stinespring dilation and comparable to the classical method.Finally we demonstrate the application of the quantum algorithm to an amplitude damping channel with implementation on the IBM Qiskit [28] quantum simulator and the IBM Q 5 Tenerife [29] quantum device.

II. 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 i  as an input state i v in a given basis and want to evolve: First note that each Kraus operator k 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 sup contraction (see the supplementary information (SI) for the proof).By the Sz.-Nagy dilation theorem [30,31], any contraction A of a Hilbert space H has a corresponding unitary operator A U in a larger Hilbert space K such that: , where H 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 A 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 A with 1 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 [30].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 k i M v with unitary gates we construct a 1-dilation A U of the form in Eq. ( 7) with k  A M and evolve: U can be further decomposed into sequences of two-level unitary gates with a procedure illustrated in Ref. [3,32].The two-level unitary gates can be used as elementary gates in an optical beamsplitter setup [32][33][34][35] 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 twolevel 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,32].For each k M U of the form in Eq. ( 7) with M , and the total gate count is m n  we need 2 n circuits to evolve all the k M 's, but this can be done in parallel, thus the complexity of each single circuit is greatly reduced.Now to put the pieces back together, the full evolved density matrix can be assembled by and measured by quantum tomography.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 from the output  not by quantum tomography, but by projection measurements instead, thus saving considerable computational costs.To get the populations of states in the current basis, note the diagonal vector of Eq. ( 10) means we can obtain the diagonal element by applying a projection measurement on the th j entry in the first n -dimensional subspace of Using an optical setup such as in Ref. [35] the probability of measuring each entry in v 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   which gives the populations in the current basis.
The additional unitary T applied to where each  L with two 2-dilations of the form in Eq. ( 8) with and we can evaluate O  by: where the trace of

III. 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, 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: To calculate the populations in the initial basis   0 , 1 we can construct 0 M U and 1

M
U of the 1- dilation form in Eq. ( 7) with Using an arbitrary initial 1 1 and assuming the physical composition where 2 m  matching the size of the vectors with the dilation k M U .We set    8) and apply them to the initial state i v in the form of Eq. ( 21) with 4 m  .Numerically calculating the output vector will give us O by Eq. ( 16) and ( 17).The results are shown as the smooth line in Figure 3: The overall simplicity of the unitary dilation gates used in our algorithm allows us to further demonstrate it using the IBM Qiskit [28] quantum simulator and the IBM Q 5 Tenerife [29] quantum device.So far we have been using the 2-level unitaries obtained from the 0 M U and 1 M U gates as elementary gates in our complexity evaluation.This is possible if we use an optical platform for quantum computation involving beamsplitters [32,34].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  We implement these circuits on both the simulator and the quantum devices and show the results as the crosses (simulator) and dots (device) in Figure 1  matrices.The decomposition of 2-level unitaries into 1-qubit and 2-qubit elementary gates grows rapidly with the overall size of the 2-level unitaries because of the complexity to implement multiqubit control gates.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 Figure 3. On the other hand the 2-level unitaries can be more easily implemented with a multiport photonic device as demonstrated in Refs.[34,35].We will seek to demonstrate our methods on such a preferred device in a future study.
In Figure 1, Figure 2 and Figure 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.

IV. Discussion and conclusion
In this work we have presented a quantum algorithm for evolving the density matrix   M M and extracting physical information from the output.The method takes each physical composition i  of  '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 the method needs to be modified.As details shown in the SI, the modified method uses the flattened vector of the initial density matrix as the input and requires   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 [37,38] -the latter helps to understand natural light harvesting complexes and exploit quantum coherence effect to improve light harvesting efficiency in artificial photocells [37,39].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. V.
This supplementary document supports the discussion in the main text by providing technical details.Section 1 provides the proof that each Kraus operator k M is a contraction.Section 2 proves that O  is a contraction and positive-semidefinite.Section 3 presents the modified algorithm for the initial state given in a general matrix form.Section 4 lists the quantum circuits used on the IBM Qiskit simulator and the IBM Q 5 Tenerife device.

Proof that each Kraus operator k M is a contraction
As defined in the main text, an operator A is a contraction if it shrinks or preserves the norm of any vector such that the operator norm sup 1

contraction and positive-semidefinite
To prove 2 Where we have used the triangle inequality of the operator norm and the fact that Note we could have made 2 Therefore O  is indeed positive-semidefinite.

Modified algorithm for the initial state in a general matrix form
In this section we present the quantum algorithm that evolves   t  with the initial  given in a general matrix form: We first flatten  into a vector form: , ... , , , ... , , ... ... , , ... , for which the norm of  v is given by: where  HS is the Hilbert-Schmidt norm of the density matrix  .Eq. ( 26) connects the norm of  v to the purity     M M which are the populations of the final system state in the basis currently in use.
Although the off-diagonal elements of where each   † k t  L L is positive-semidefinite and its diagonal elements can be obtained by projection measurements.† L is obviously a contraction because O  is a contraction.To realize we need all four matrices to be in the 4-dilation form: † † 0 0 0 where Now after With an initial 1 1 1 1 3 4 HS   , the input state is:   0 1 1 , 0,..., 0 1,1,1,3, 0,..., 0 2 3 where 8 m  here is the number of zeros after T  v .We use the same parameters as in the main text: information carried by   t  --e.g. the populations of different quantum states and the expectation value of an observable 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

M
If k has the dimension n by n , then the 1-dilation k M U is 2n by 2n .Now k M

U
is incidentally also2  2n n  for each k and i , counting all the multiplications and additions needed for a matrix multiplying a vector.As a comparison, the conventional Stinespring dilation converts the whole quantum operation   † by considering the tensor product space of the initial  and an auxiliary environment[3].The dimension of U is given by n m  where m is the total number of Kraus operators k M in the sum † k k k   M M .In the most general case 2 unitary gates for the gate decomposition.We see that our procedure involving the Sz.-Nagy dilation offers a major simplification over the conventional is achieved by breaking the quantum operation   † them separately.Of course for 2 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 the additional level of dilation for k M increases the quantum gate count to k and i .A classical overhead (independent from the k and i counts) Cholesky decomposition[36] should also be counted towards the total cost of the quantum algorithm.On the other hand the classical complexity of evaluating † k i L M v is 2 3n n  for each k and i (taking into account that † L is upper triangular) picosecond, and obtain the populations of the ground and excited states from the first two entries of k i M U v .The results are shown as the smooth lines in Figure 1.

Figure 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 Figure4for an example and the SI for all the circuits).

Figure 2 .U
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 .
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.

Figure 4 ,
Figure4, and a full list of the circuits can be found in the SI.

Figure 4 .
Figure 4. Showing the quantum circuit for

and Figure 2 .
To obtain O  we construct a circuit of 3 qubits with an average of 182 elementary gates for the † k i M L U U v operation (see the SI).The drastic increase in the number of gates is due to the increased size of the †

2 O 6 O
and evolves it through minimal Sz.-Nagy dilations of the Kraus operators k M .The 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   n for each k and i , which is a significant improvement over the   n scaling of the conventional Stinespring dilation.The requirement of knowing the initial physical composition as a probability-weighted mixture of not necessarily orthogonal i

3 O
n 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 k M 's for the operator sum representation are different for different models.
successfully obtained O for the original observable.The   †  L I and   †  I L gates ( † L is upper triangular requiring reduced number of two-level unitaries for the decomposition) plus the additional two levels of dilation for M k and N k increase the quantum gate should be added as an overhead cost (independent from the k count) when evaluating the total cost of the quantum algorithm.In the meanwhile the classical complexity of evaluating k .Next we demonstrate the method proposed above on the same amplitude damping model used in the main text.To calculate the populations in the current basis   0 , 1 we can construct k U M and k U N of the 2-dilation form in Eq. (8) using k k

.,U 4 .
with a time step of 10 picosecond, and obtain the populations of the ground and excited states from the first and fourth entries of 0 The results are the same as the smooth lines in Figure1in the main text.To calculate the populations in another basis   obtain the populations of the  states.The results are the same as the smooth lines in Figure2in the main text.Now we evaluate the expectation value of an observable O for 2 N of the form in Eq. (31) and apply them to the initial state 0 v in the form of Eq. (36) with 16 m  .Numerically calculating the output vector will give us O by Eq.(33).The results are the same as the smooth line in Figure3in the main text.Quantum circuits used for the IBM Qiskit and Q 5 Tenerife device.To implement Method 2 on the IBM simulator and quantum device, we further decompose the unitary gates into 1-qubit and 2-qubit elementary gates and construct the quantum circuits.All the gates used below are standard.For each circuit, only the  parameter of the the time evolution.The circuits below show the 3 U at the last time step at 991ps.First for evolution in the original basis: the evolution with O  evaluation, due to the large number of gates involved we only show the Although the off-diagonal elements of   (10)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 -dimensional space of can potentially save measurement costs by reducing the number of inquiries needed on the output vector.Finally O is calculated with: