Commutation simulator for open quantum dynamics

Recent progress in quantum simulation and algorithms has demonstrated a rapid expansion in capabilities. The search continues for new techniques and applications to exploit quantum advantage. Here we propose an innovative method to investigate directly the properties of a time-dependent density operator $\hat{\rho}(t)$. Using generalised quantum commutation simulators, we can directly compute the expectation value of the commutation relation and thus of the rate of change of $\hat{\rho}(t)$. The approach can be utilised as a quantum eigen-vector solver for the von Neumann equation and a decoherence investigator for the Lindblad equation, by using just the statistics of single-qubit measurements. A simple but important example is demonstrated in the single-qubit case and we discuss extension of the method for practical quantum simulation with many qubits, towards investigation of more realistic quantum systems.

A century ago, very early in the development of quantum mechanics, commutation relations emerged in various crucial roles [1,2].Pairs of non-commuting operators (e.g., the position and momentum operators for a particle) describe the complementary nature of their corresponding physical properties, leading to uncertainty relations between these quantities for quantum systems.Commutators also underpin the time evolution of quantum systems, whether this be of general operators in the Heisenberg picture (or the relevant part of the Interaction picture), or the system density operator in the Schrödinger picture, where the time dependence resides in the quantum state or density operator.In this latter picture, the quantum state of a system given by |ψ(t) evolves according to the Schrödinger equation (with units where = 1) [3], given by where Ĥ is the system Hamiltonian.For an initial state defined as |ψ(0) = |ψ 0 at t = 0 and a time-independent Ĥ, the evolution from 0 to t is determined by the unitary operator Û (t), such that |ψ(t) = Û (t)|ψ 0 = e −i Ĥt |ψ 0 .
An equivalent alternative description is via the density operator, so defining this as ρ(t) = Û (t)|ψ 0 ψ 0 |( Û (t)) † the evolution is given by the von Neumann equation [4] as where the commutation relation between x and ŷ is given by [x, ŷ] = xŷ − ŷx.
The form of the von Neumann equation is very interesting because the time-dependence of the system is expressed directly in terms of the commutation relation between the density operator and the Hamiltonian.The density operator approach provides a direct statistical representation because the diagonal parts of d dt ρ(t) give the rate of change of the system probability density.These always correspond to real numbers, which can be measured for the actual physical system, either through repeated measurements on an identically prepared and evolved single pure system, or through measurements on an ensemble of identical systems all equivalently prepared and evolved.We refer to these equivalent approaches as an "ensemble measurement".The density matrix approach can also be used to incorporate classical uncertainty (lack of knowledge), in addition to quantum superposition, via (finite-entropy) mixtures of pure quantum states.In this work we will use the density operator approach ρ(t), both from the perspective of the reversible von Neumann equation ( 2) but also to provide scope for the inclusion of classical uncertainty and irreversible evolution.
In quantum theory, the irreversibility inherent in open systems-those coupled to additional environment degrees of freedom-can be modelled by modification and addition of noise terms to either the Heisenberg equation for system operators or the Schrödinger equation for system states [5].However, the density matrix approach forms a very important method for investigating the dynamics of open quantum systems, beyond just the Schrödinger equation.The Lindblad master equation is a very widely used and applicable example.This commonly describes an open system interacting weakly with its environment, describing the effects of the environment on the system (generally, decoherence mechanisms) using Lindblad operators Lj .These operators modify the von Neumann equation (2) to where the anti-commutation relation between x and ŷ is given by {x, ŷ} = xŷ + ŷx [6,7].In general the Lindblad operators are not Hermitian and act to introduce decoherence to the system, changing its entropy.The particular case of Hermitian Lindblad operators can be used to model quantum measurements, or noisy external source terms in the system Hamiltonian.For Lj = 1 1 (1 1: identity operator), the Lindblad terms disappear and only the unitary term i ρ(t), Ĥ survives, thus returning to the von Neumann equation and the unitary evolution of a closed quantum system in time.
Quantum algorithms for simulating the Schrödinger equation have been developed extensively since the first rigorous idea of quantum simulation in 1996 [8].However, the preparation of pure states would appear to intrinsically prohibit research progress on the general simulation of mixed density matrices in quantum circuits.The von Neumann equation is normally treated as an equivalent equation to the Schrödinger equation, with no decoherence, but containing direct physical interpretation, such as the diagonal elements describing the probability density of the system in the chosen basis.For open quantum systems, several approaches have recently been developed using a vectorized density operator [9] or a trace-out of environmental qubits [10][11][12][13] as the Lindblad-type equations (also named dissipative quantum computation [14,15]).
In this work, we propose a novel method to directly compute, or simulate, matrix elements of d dt ρ(t), by measuring expectation values of the commutation relation in the von Neumann equation ( 2) and the more general Lindblad equation (3).Consider the case where the system of interest comprises L qubits, so the density operator ρ(t) can be represented by a 2 L × 2 L matrix.Our approach provides the (diagonal and off-diagonal) matrix elements d dt ρ n,m (t) = n| d dt ρ(t)|m with n, m = 0, ..., 2 L − 1 ranging over a suitable basis of the system.So, for example, if we seek the expectation of the rate of state change in time, given by Φ| d dt ρ(t)|Φ for some chosen reference state |Φ , we can perform quantum processing to determine this by measuring the expectation value of the commutator in the von Neumann equation ( 2), given by i Φ|[ρ(t), Ĥ]|Φ in the case of closed quantum systems.For the off-diagonal terms, we can compute i Φ|[ρ(t), Ĥ]|Φ by a sum of expectation values given by another controlled-operator gate, with operator Â for |Φ = Â|Φ .For the case of open quantum systems, it is required to perform additional quantum processing to compute the extra Lindblad terms that depend on the Lj .This paper is constructed as follows.In the next section we present the algorithm for simulation of the quantum commutator.Then we present the application of the generalised algorithm to the von Neumann and the Lindblad equations.A specific example is then given, followed by summary and conclusion.

ALGORITHM FOR QUANTUM COMMUTATION SIMULATION
We first provide the protocol describing the algorithm, followed by a detailed explanation of the quantum commutation simulator.The simulation is built upon the following resources, as employed in Fig. 1: The system S, assumed to be of dimension 2 L , or 2 L × 2 L in density matrix form; a separate reference system M of the same size as the system; a separate control qubit C. The state of the total system is denoted by |Ψ .With reference to the full quantum circuit shown in Fig. 1, the protocol for the simulation runs as follows.
( ) 1. Schematic of the generalised quantum circuit to simulate the expectation value of the commutation relation.Three controlled-operators are given by operators N , Â and M and the total number of qubits required in the simulator is 2L + 1 to describe an L-qubit system.The detailed protocol is described in Section .

Protocol
1. Initialise the total system state |Ψ 0 as a product of system state |ψ 0 S , reference |Φ M and control qubit |+ C .
3. Apply controlled-operator gate N between control C and system S as well as two controlled-operator gates Â and M between C and reference M to produce |Ψ 2 .
4. Apply a block controlled-SWAP gate from control C between system S and reference M to produce |Ψ 3 [16].
5. Apply a Hadamard gate Ĥ to control C to produce |Ψ 4 .
6. Measure a single qubit in C in the Pauli-Z gate (the computational) basis, to obtain the expectation value Ẑ .

Quantum commutation simulator
Using qubit terminology, we first explain the operation of the generalised quantum commutation simulator, shown in Fig. 1.Generally, by default, qubits are assumed to be initialised in |0 for a quantum circuit, but here we assume some additional preparation.The control qubit C is prepared in the state |+ C = Ĥ|0 C using a Hadamard gate Ĥ.As shown in Fig. 1, there are two L-qubit states, for the system S and the reference M .For the system S, the initial state |ψ 0 is assumed to be created by a suitable prior quantum circuit, specified by the chosen initial conditions of the target problem to be simulated at t = 0.For the reference M , we can simply utilise one of the computational basis states (e.g., |0 ⊗L ), or any other interesting reference state |Φ to be evolved dynamically.Note that only the control qubit C is measured at the end of the process and thus the outcome provides us with the expectation value of a quantum operator for the other degrees of freedom, effectively given by a 2 L × 2 L matrix for each of system S and reference M .
For the simulation, a total input quantum state is therefore prepared as In the first step, we apply a single-qubit gate R(χ) on qubit C.This is a rotation by χ around the qubit Z-axis with an additional phase, R(χ) = e −iχ/2 RZ (χ).In parallel, the unitary operator Û (t) = e −i Ĥt is applied to the system |ψ 0 , giving the total state as Next, we apply three controlled-operator gates in order to produce |Ψ 2 .These are controlled by qubit C and applied to the system qubits S and to the reference qubits M .The operators are represented by N , Â and M and these controlled operators are each 2 L × 2 L matrices operating on the L reference qubits in either S or M .We discuss the roles of these operators and their implementation in more detail in the next section of the paper.
For the next step of the simulation, to produce |Ψ 3 , we apply a block controlled-SWAP gate between the system S and the reference M , as shown in Fig. 1 [16].Following this, a Hadamard gate Ĥ is applied on control qubit C, to produce the final (entangled) state of the total system, before measurement Ensemble measurement of the control qubit C in the computational basis will generate the probabilities of the outcomes |0 C and |1 C , defined respectively as P 0 and P 1 .The expectation value of Pauli operator Ẑ is equal to Since the results are given as a difference of scalar values (details in Appendix ), we can interchange these to reformulate the expectation value of Ẑ for qubit C as defining a new quantum operator as Note that this new operator contains actions of the controlled gates Â and M , in a manner that depends on the chosen rotation angle χ.
As an example, for an identity N = Â = 1 1 and M being a Hermitian operator M † = M , the value of χ can then determine whether the result delivers the expectation value of the commutation or anti-commutation relation between the time-dependent density matrix ρ(t) and the operator M .These follow from the statistics of single-qubit measurements through For N = 1 1 and Â = We now apply the simulation to the case of the von Neumann equation (2).We therefore focus on the case of M = Ĥ as the system Hamiltonian operator.If the Hamiltonian is decomposed in terms of its eigenvalues λ l and eigenvectors |λ l , as then the unitary evolution generated by Ĥ is given by The von Neumann equation ( 2) can then be written in the form where α k = λ k |ψ 0 and α * l = ψ 0 |λ l for eigenvectors |λ k and |λ l .In (11), the matrix of d dt ρ(t) comprises all off-diagonal terms in the eigenvector representation of the von Neumann equation.
Note that, independent of the reference state |Φ , the expression i Φ|[ρ(t), Ĥ]|Φ = 0 if either α k or α l is zero in each term of the summation.This follows if the initial state |ψ 0 is chosen to equal one of the eigenvectors of Ĥ in the quantum commutation simulator, rather than being a superposition of two or more eigenvectors with different eigenvalues.Therefore, we can search for eigenvectors of Ĥ by tuning the parameters in |ψ 0 .In general, if |ψ 0 is tuned to include amplitudes of both eigenvectors |λ k and |λ l , the expectation value i Φ|[ρ(t), Ĥ]|Φ oscillates in time with the frequency λ l − λ k and the amplitude ( Clearly this is conditional on |Φ being of a suitable form so that the amplitude is non-zero. Interestingly, if |Φ is selected as an eigenvector |λ m of Ĥ, the expectation value of the right hand side of equation ( 11) is always zero, regardless of the time t and the form of |ψ 0 , due to the orthogonality of eigenvectors, giving This result implies that |Φ = |λ m is a stationary state in the unitary dynamics and one of the eigenvectors of the closed system described by Ĥ.
In addition, the off-diagonal elements of the commutator in the eigenstate basis can be computed.The coherence between the two eigenvectors |λ k and |λ l is given by λ k ρ(t), Ĥ λ l = (λ l − λ k ) α k α * l e i(λ l −λ k )t for λ l |λ k = 0 with l = k.Note that again the amplitude is proportional to the difference of the two eigenvalues λ l − λ k , as well as this difference setting the frequency of the time-dependent part.
Consider now the specific case where the reference state |Φ is chosen as a computational basis state |n (n = 0, ..., 2 L − 1 for L qubits), instead of a general computational basis state.Then the expectation value defines the associated probability density rate d dt ρ n,n (t).Based on the quantum commutation simulator in Fig. 1, taking N = Â = 1 1 and M = Ĥ and with χ = π/2 in (8), we can calculate the diagonal elements of d dt ρ(t) as Thus, if we are able to implement the Hamiltonian operator Ĥ in the contolled-operator gate of Fig. 1, the value of d dt ρ n,n (t) follows from the statistics of single-qubit measurements of C, through n| Ẑπ/2 1 1 |n .Moreover, for N = 1 1 and M = Ĥ with χ = π/2, if Â is not set as the identity operator but chosen e.g. as the translation operator Â [16,17], such that Â |n = |n + 1 , we can also investigate any off-diagonal matrix elements, such as d dt ρ n,m (t) for n = m.Specifically, if we replace Â → ( Â) p , as p actions of the translation operator Â for the L-qubit system (p = 1, ..., L − 1), we can compute all the off-diagonal elements of (11), such as d dt ρ n,n+p (t).In any actual implementation of the quantum commutation simulator, one challenge is the construction of the (controlled) operator Ĥ in the quantum circuit, because the Hamiltonian operator is in general not a unitary operator.To resolve this issue, we use the fact that the Hermitian operator Ĥ can be split into a sum of unitary gates (e.g., Pauli matrices for qubits).Thus, we need to build a decomposition of the Hamiltonian Ĥ = k c k Ûk , with coefficients c k where each Hamiltonian term Ûk is represented by a unitary gate.Then, the additivity of commutation relations can be used in equation (13) to give ρ(t), Ĥ = k c k ρ(t), Ûk .Through this decomposition it is then feasible to simulate the required commutator, via the implementation of controlled-Ûk gates in appropriate quantum circuits.

Dynamics of an open quantum system: decoherence investigator
The Lindblad equation ( 3) can describe a decoherence mechanism for an open quantum system weakly coupled with a large environment [6,7].As shown in Appendix A, the derivation of the Lindblad master equation can be described by a series of short time periods δt.In more detail, the initial density matrix is given in ρ0 = |ψ 0 ψ 0 | and coherently evolves to ρ(δt) in the period δt, with the dynamics of the quantum system δ ρ(δt) δt given by the von Neumann simulation.During the period between δt and 2δt, the time evolution description of the quantum system is split into two terms: a unitary term with the commutation relation (the same as the von Neumann case) and a sum of the Lindblad terms including the anti-commutation relation [18].With these steps, we are able to rewrite the Lindblad equation at t = 2δt in the form The key point is that we can compute all the expectation values of the right hand side terms of equation ( 14) in the quantum commutation simulator, using both commutation and anti-commutation relations.
In order to obtain the first part of the Lindblad term, let us choose the controlled-operator set { N , M , Â} = { Lj , L † j , 1 1} in equation ( 5).In the case of Â = 1 1, the expectation value of the diagonal elements is given by which becomes n| Lj ρ(δt) L † j |n for |Φ = |n as a computational basis state.For the off-diagonal elements of the first Lindblad terms, we replace Â with ( Â) p .This enables us to generate the expectation value of any off-diagonal part of the first Lindblad term, which is given by a combination of two different expectation values in the form where Φ|Φ = 0 due to |Φ = ( Â) p |Φ .For |Φ = |n in the computational basis, n| Lj ρ(δt) L † j |n + p is given by the combination of n| Ẑχ A |n for χ = 0, π/2.Lastly, we calculate the anti-commutator component of the second part of the Lindblad term.For the diagonal elements, if { N , M , Â} = { L † j Lj , 1 1, 1 1}, the expectation value of the anti-commutation relation is given by For the off-diagonal elements, if we select { N , M , Â} = { L † j Lj , 1 1, ( Â) p } for |Φ = Â|Φ , a combination of two expectation values gives the result Thus, adding the results of equations ( 18) and ( 19), we obtain the expectation value of the off-diagonal part in the second Lindblad term.Therefore, the expectation values can be computed for all the elements of δ δt ρ (2δt) in equation ( 14).

1-QUBIT EXAMPLE Theory: amplitude damping for spontaneous emission
Let us examine the algorithm for the quantum commutation simulator in the case of a single-qubit system under amplitude damping.We consider a general initial system state given by |ψ 0 = cos θ 2 |0 +e iφ sin θ 2 |1 , with Hamiltonian Ĥ = − ω 2 Ẑ.This gives the unitary operator Û (t) = exp(−i Ĥt) = exp(i ω 2 Ẑt).Writing the evolved density matrix in the form ρ(t) = Û (t)|ψ 0 ψ 0 |( Û (t)) † , the right side of the von Neumann equation at time δt is given by For the Lindblad equation with a single Lindblad operator L, we utilse equation ( 14) at t = 2δt and obtain a 2 × 2 matrix such as

Quantum commutation simulation for the von Neumann equation
For the unitary time evolution of a single-qubit system in the von Neumann equation, we first choose N = Â = 1 1, M = Ẑ = − 2 ω Ĥ and Û (t) = exp(i ω 2 Ẑt) in Fig. 1.In equation ( 13), the expectation value Φ| Ẑπ/2 1 1 |Φ provides the rate of change of the probability density for |Φ , so that for the computational basis states |Φ = |0 and |1 .Therefore, regardless of the initial state parameters θ and φ in |ψ 0 , these expectation values are always zero.This implies (as expected) that the two reference states (|0 and |1 ) are stationary states in time and the eigenvectors of Ĥ.
In addition, by choosing Â = X, we can compute the off-diagonal elements, given by the combination of four expectation values (each one for either a real or imaginary part in n|[ρ(t), M ]|m ), such that For this case of qubit amplitude damping, Fig. 2 demonstrates examples of relevant non-zero matrix elements of the Lindblad equation Eq. ( 21), for specific parameter choices and as functions of t and θ.In Fig. 2(a) we observe that the real part of Eq. ( 24) oscillates in time for θ = 0 and π.The static cases (θ = 0 , π) identify the eigenvectors |0 and |1 of H, whereas ω follows from the period of the oscillation in sin(ωt) (with ω = −2 and φ = 0 in this example).Thus, the eigenvalue information ω can be also extracted from the landscape of measurement outcomes.
Therefore, for the single-qubit open quantum system under spontaneous emission, the full Lindblad equation (21) at t = 2δt is given by Let us make various comments at this stage.
• The result of equation ( 35) is exactly that which we would expect for amplitude damping of a qubit.
• We note that the rates of change of the diagonal density matrix elements are no longer zero at t = 2δt, due to the decoherence from the amplitude damping in the open quantum system.
• The most important point to stress is that we have calculated the terms that contribute to equation (35) in a manner that is amenable to quantum simulations, through various specifically chosen cases of Fig. 1.This provides a route for simulation of Lindblad evolution.
• We chose an initial (t = 0) pure state, but an initial mixed state could be simulated as a mixture over a suitable decomposition of pure states.
Finally, to fully justify our first comment above, we construct the actual resultant density matrix ρ(t).In effect, although we used the terminology of a t = 0 initial state, what equation (35) represents is the relationship of the rates of change of the density matrix elements at some given time to the actual matrix elements at that time.Therefore, equation ( 35) provides (in this case) the four anticipated first order differential equations for the density matrix, that can be integrated to give the expected result for a qubit with amplitude damping Analogous to the controlled-Hamiltonian gate in Section , the Lindblad operators L and L † are used here.The non-Hermitian operators describing the spontaneous emission with damping parameter κ are given by L = Then, the implementation of these operators is always feasible in the commutation simulator, given by the decomposition of the operators ij c jk Σj ρ(δt) Σk and Σ denotes the basis for 2 × 2 Hermitian matrices, given by the three Pauli operators plus the identity 1 1. system, via the Lindblad equation.For a large quantum system, the von Neumann method is beneficial for computing a transition rate between two specific quantum states and this result can be also used for the study of its open quantum system.For example, although the sizes of the system qubits state |ψ 0 and the reference qubits state |Φ is each given by a 2 L × 1 column vector for L qubits, the probability transition rate between two specific states is always given by the expectation value of the off-diagonal elements in a 2 × 2 matrix form.Correspondingly, the changes of each state probability follow from the diagonal elements.
Clearly in general Lindblad evolution is irreversible, with a change in mixture (entropy) of the density matrix.In the simulation approach this changing mixture is introduced because the outputs of different simulations have to be combined, and each of these simulations involve measurements, to compute the expectation values that comprise the various matrix elements of the (rate of change of the) density matrix.
Extensions of the Lindblad equation simulation method have the potential to generate innovative approaches for general purpose master equations.This is because the simulation approach preserves a probabilistic interpretation for the system even for open systems, generating the evolution of the system probabilities from the diagonal elements of a density.The investigation of possible quantum advantage in such open system applications will be an interesting topic for future study, for example to investigate the quantum speed limit of simulating an open quantum system [19,20].