Digital simulation of convex mixtures of Markovian and non-Markovian single qubit Pauli channels on NISQ devices

Quantum algorithms for simulating quantum systems provide a clear and provable advantage over classical algorithms in fault-tolerant settings. There is also interest in quantum algorithms and their implementation in Noisy Intermediate Scale Quantum (NISQ) settings. In these settings, various noise sources and errors must be accounted for when executing any experiments. Recently, NISQ devices have been verified as versatile testbeds for simulating open quantum systems and have been used to simulate simple quantum channels. Our goal is to solve the more complicated problem of simulating convex mixtures of single qubit Pauli channels on NISQ devices. We consider two specific cases: mixtures of Markovian channels that result in a non-Markovian channel (M+M=nM) and mixtures of non-Markovian channels that result in a Markovian channel (nM+nM=M). For the first case, we consider mixtures of Markovian single-qubit Pauli channels; for the second case, we consider mixtures of Non-Markovian single-qubit depolarising channels, which is a special case of the single-qubit Pauli channel. We show that efficient circuits, which account for the topology of currently available devices and current levels of decoherence, can be constructed by heuristic approaches that reduce the number of CNOT gates used in our circuit. We also present a strategy for regularising the process matrix so that the process tomography yields a completely positive and trace-preserving (CPTP) channel.


Introduction
Simulating large, complex quantum systems with classical computers is a computationally challenging problem.The computational resources required to perform these simulations would scale exponentially with the number of quantum particles, and these simulations would quickly become classically intractable.For this reason, Quantum Simulation, the simulation of quantum systems with quantum computers, has been proposed [1,2].The computational resources required to use quantum computers would scale only polynomially with the number of quantum particles.Since this discovery, quantum simulation has become one of the main motivations for developing quantum computers [3].
Researchers have developed many algorithms to simulate quantum systems using quantum computers [4][5][6][7][8][9][10][11][12][13].However, these algorithms are best suited for fault-tolerant settings.In these settings, quantum computers would provide a clear advantage over classical computers in the simulation of quantum systems.However, despite the advancements made in quantum hardware, much more work needs to be done to reach fault-tolerant settings.This, and the growing accessibility of Noisy Intermediate Scale Quantum (NISQ) computers through the cloud from platforms such as the IBM Quantum Experience (IB-MQE), has inspired interest in the use of NISQ devices for simulating quantum systems [14][15][16][17].
In these settings, the simulations are restricted by the size of the systems, various error rates and noise sources.Despite this, even currently available quantum computers provide versatile test beds for various theories in quantum physics [18][19][20].
A majority of recent work has been devoted to the simulation of closed quantum systems in the NISQ era [14][15][16][17].However, there has been less work done in the simulation of open quantum systems in the NISQ era [20][21][22].Open quantum systems [23,24], are systems that are allowed to interact with their environment and to simulate them we need to be able to simulate their evolution on a quantum computer.
A master equation describes the dynamics of an open quantum system.The solution of this master equation is a dynamical map, also known as a quantum channel, which describes the evolution of an open quantum system.Under certain assumptions, such as the Born-Markov approximation, the Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) form of the master equation can be derived [25,26].The GKSL form of the master equation describes Markovian dynamics, where all memory effects are neglected.Non-Markovian dynamics differs from Markovian dynamics in that it allows information to flow back into the system from the environment and does not neglect memory effects [27].While everyone agrees that the time-independent GKSL generator with non-negative damping constants describes the quantum Markov process, there are several approaches to defining the non-Markovianity in quantum physics [28][29][30][31].This makes the study of non-Markovianity in quantum physics a highly non-trivial problem.
There has been much interest in the study of the non-Markovian dynamics of an open quantum system [28][29][30][32][33][34][35][36][37].There are many different descriptions of non-Markovianity [28,29,37] however, in this work, we make use of CP divisibility [29] to characterise a channel as Markovian or non-Markovian.We say that a channel is Markovian if it is CP-divisible and non-Markovian if it is CP-indivisible.
Recently, simulations of open quantum systems have been used to understand non-Markovian dynamics further [20,38].There is a need for more algorithms for simulating open quantum systems with NISQ devices and more experiments to test existing algo-rithms [20].This would demonstrate that NISQ devices can be used for more practical problems, taking us further toward using quantum computers as practical computing devices.
Recently, there has been much interest in studying the mixtures of quantum channels and how mixing channels with specific properties leads to a new channel with counterintuitive and surprising properties [39][40][41].In [41], the author studies mixtures of commutative, unital and Markovian quantum channels, and they show that Markovianity can be recovered by mixing a certain number of phase covariant channels.The work done in [39] calculates the number of negative decay rates that qudit channels can have while still being physically legitimate quantum channels.The author of [40] provides connections between the non-Markovianity degree of general phase damping qubit maps and their legitimate mixtures.Necessary and sufficient conditions for Pauli maps to satisfy divisibility criteria are formulated.In [39], mixtures of Pauli channels are studied, but a prescription for how to design these mixtures of channels for simulations and experimental demonstrations is not provided.In this work, we will provide a prescription for designing mixtures of depolarising channels to mix two non-Markovian depolarising channels and yield a Markovian channel.We also show that it is impossible to mix two Markovian depolarising channels and obtain a non-Markovian channel.Of interest to us is the work done by [38], where they study mixtures of Markovian channels that lead to Non-Markovian channels and vice versa.
In this work, we use NISQ devices, provided by the IBMQE through the cloud, to simulate the mixtures of Markovian quantum channels that yield a Non-Markovian channel (M + M = nM) and the mixtures of Non-Markovian channels that yield a Markovian channel (nM + nM = M).We will also use the Python package Qiskit to interact with the NISQ devices [42].For the case of M + M = nM, we consider mixtures of Pauli channels because of their relevance to bit-flip and phase-flip error encountered on quantum devices, while for the case of nM + nM = M, we consider mixtures of depolarising channels because of the common occurrence of depolarising noise on quantum devices.
Each mixture was then simulated on a NISQ device.We show that efficient circuits that account for the topology of currently available devices and current levels of decoherence can be constructed by carefully considering device properties when applying the unitary dilation, such as Stinespring dilation [43].These NISQ device motivated circuit constructions provide an improvement in the quality of the results obtained in simulations.
The channel that resulted from the simulation was then reconstructed.Previous attempts at simulating convex mixtures employed Maximum Likelihood Estimation (MLE) [38] to reconstruct the channel.This has since been shown, experimentally, to be unideal.In this work, we reconstruct the channel by solving a convex optimisation problem with problem-specific constraints [44].With this approach, a CPTP channel is obtained.The simulation was then verified by characterising the reconstructed channel as either Markovian or Non-Markovian using the CP-divisibility criteria [29].We also show that the reconstructed channels have high fidelity to the theoretical channel and have Choi matrices with the required trace norm of one.This indicates that the channels simulated are physical.The rest of this paper is outlined as follows.In Sect.2, we introduce some background information and theory related to quantum channels and their generators, and we also present the theoretical design of the quantum channels to be simulated for both cases M + M = nM and nM + nM = M.In Sect.3, we discuss the simulation of the channels on the NISQ device.We discuss the quantum circuits designed to simulate the quantum channels in Sect. 3 A. In Sect. 3 B, we outline the method for process tomography and the convex optimisation technique we use to reconstruct a CPTP map.Section 4 will describe the method to characterise the simulated channels as Markovian or non-Markovian using the CP-divisibility criteria.The results from the simulation will be presented in detail in Sect. 5. Lastly, in Sect.6, we make some concluding remarks and discuss possible extensions of this work.

Preliminaries
The evolution of an isolated (closed) quantum system is described by the von-Neumann equation, where H S is the Hamiltonian of the system and ρ S (t) is the density matrix of the system at some time t ≥ 0. Formally, the solution to equation ( 1) is a unitary transformation, where U(t) = exp(-itH S ).If the system is not isolated and is allowed to interact with some environment then the system and environment together undergo unitary evolution as an isolated system.To describe the dynamics of the system we need to "trace out" the environmental degrees of freedom obtaining the evolution of the open quantum system.
Without loss of generality one can assume that the total state ρ tot (t), at t = 0, of the system and environment is a product state i.e.
where ρ E is the density matrix of the environment and ρ S is the density matrix of the system.The total unitary evolution U tot (t) of the system and environment is given by, U tot (t) = exp(-itH tot ) where H tot is the Hamiltonian of the system and environment and is, where H S is the system Hamiltonian, H E is the Hamiltonian of the environment and H I is the Hamiltonian describing the interaction between the system and environment.The total system and environment undergo unitary evolution i.e. ρ tot (t) = U tot (t)ρ tot (0)U † tot (t), to obtain the dynamics of just the system we must trace out the environment using the partial trace so that, Equation ( 5) allows us to describe the effective dynamics of the system using a dynamical map t where, The dynamical map t where t ≥ 0 and 0 = 1 are a family of single parameter completely positive and trace preserving (CPTP) maps.From this point on since we are only interested in the dynamics of the system we can drop the subscript and represent the state of the system at some time t ≥ 0, with just ρ(t) so that, if ρ(0) is the initial state of the system then ρ(t) = t ρ(0) [23].
A dynamical map is also referred to as a quantum channel.These shall be used interchangeably throughout this work.One can assume that the map t , in most practical cases, satisfies the time-local master equation, where L(t) is the time-local generator and has the known form, where Ĥ(t) is the time-dependent Hamiltonian, γ k (t) and Vk (t) are the time-dependent decay rates and noise operators respectively.The form of the generator in ( 2) is very general and gives both Markovian and non-Markovian dynamics.It is widely accepted that the famous Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) form of the generator corresponds to Markovian dynamics [25,26].The GKSL form of the generator is, Since the GKSL form of the generator corresponds to Markovian dynamics, all memory effects are neglected.A crucial part of this work will involve classifying the simulated dynamical maps t as either Markovian or non-Markovian.As stated above, we shall use the CP divisibility criteria to classify a dynamical map as Markovian or non-Markovian [29], which shall be outlined in more detail in the next section.Another approach in studying the Markovianity of a dynamical map is to model the total system and environment dynamics, this work does not take this approach, but one could look at [30,31] for more information.The aim of this work can be stated formally as follows.For two channels (1)   t and (2)  t we are interested in classifying the total channel: where η ∈ [0, 1].We consider the following two cases: (i) (1)  t , (2)  t are Markovian and (T) t is non-Markovian, which we shall refer to as (M + M = nM).
(ii) (1)  t , (2)  t are non-Markovian and (T) t is Markovian, which shall be referred to as (nM + nM = M).These two cases above arise from the non-convex geometry of the set of Markovian and non-Markovian channels [32].Past works that have studied this non-convex geometry have provided examples of mixtures of Markovian channels leading to a non-Markovian channel, and vice versa [34,[45][46][47][48][49].

Design of single qubit channels to simulate
Since we are only focused on single qubit channels, our channel of interest in this work will be the single qubit Pauli channel.These are the simplest non-trivial channels that we can use for our simulations.The single qubit Pauli channel is of the form, where σ 0 = 1 and p α (t) is a time dependent probability distribution such that p 0 (0) = 1 and p i (0) = 0 for i ∈ {1, 2, 3} and 3 i=0 p i (t) = 1 for all t ≥ 0. This channel was studied extensively in [36].It should also be noted that this channel also falls into a larger class of channels called random unitary channels [35].In this work, we consider a subset of single qubit channels by focusing on single qubit Pauli channels; the reason we consider these channels is that the Pauli channel in equation ( 11) is the most general unital stochastic single qubit quantum channel [36,50].This means that it is the most general single qubit channel that preserves the identity operator and the most general single qubit channel with a Kraus operator, which is proportional to the identity; this is what it means for a channel to be a stochastic channel [51].Also, since the space of stochastic channels is manifestly convex [51], taking convex mixtures of single qubit Pauli channels will again be a stochastic channel.Therefore, simulating efficiently on a NISQ device, a single qubit Pauli channel provides insight into a broad class of single qubit channels, which is why we only consider single qubit Pauli channels.
To make a statement about whether the Pauli channel is Markovian or not, we need to analyze the decay rates of the time local generator L(t) of the Pauli channel in equation (11).From the time-local master equation (7), we see that t , this tells us that we need to compute -1  t to calculate the generator [33].Let us note that: where the time-dependent eigenvalues are, with λ 0 (t) = 1 and H αβ being the Hadamard matrix defined as: We should note that λ 0 (t) = 1 and |λ k (t)| ≤ 1 for k = 1, 2, 3. Now it is clear that: where μ α (t) = λα (t) λ α (t) and in particular μ 0 (t) = 0 since λ 0 (t) = 1.Now by introducing the decay rate γ α (t) we get the local generator as: where the decay rates γ α (t) are defined as, By observing that 3 α=0 γ α = 0, we arrive at the standard form for the generator of the Pauli channel [33] as, Hence we get the expression for γ k (t) as: for k = 1, 2, 3. We use the following result from [33,35,36].The Pauli channel (11) is Markovian (i.e t is divisible) if and only if γ k (t) ≥ 0 for all t ≥ 0 and for k = 1, 2, 3.This tells us that if the decay rates of the generator of a channel are non-negative, then that channel is Markovian, and any deviation from this leads to the channel becoming non-Markovian.We know that from the time local generator L(t) we can get the channel t from the relation, where T is the chronological time ordering operator.Now we note that a linear combinations of Markovian generators L 1 (t) and L 2 (t) i.e. α 1 L 1 (t) + α 2 L 2 (t) with α 1 , α 2 ≥ 0, is also a Markovian generator.Hence Markovian generators form a convex set in the space of admissible generators.The same is not true on the channel level, since for two Markovian channels (1)  t and (2)  t it is not always the case that their convex linear combination will be Markovian, i.e. η (1)  t + (1η) (2)  t , η ∈ [0, 1], is not necessarily Markovian [34].We can now consider the two cases of convex mixing of channels outlined above, i.e.M + M = nM and nM + nM = M.We shall find Pauli channels that satisfy these two cases.

Markovian channel addition (M + M = nM)
The goal is to find two Markovian channels that, when convexly combined, yield a non-Markovian channel.We shall use the channels from [38] that demonstrate this.We start by defining the following two channels, (1) (2) where the probability p(t) = 1+e -t 2 .Using equations ( 19) and ( 18) above, we can find the generators of these channels.We find that for channels (1)  t = e tL 1 and (2)  t = e tL 2 the generators are given by: Since the decay rates in both these generators are non-negative for all t ≥ 0, the channels in equations ( 21) and ( 22) are Markovian.We can now consider the convex combination of these channels in the following way: (T) (1) (2) where (T) t is the total channel.Now for this total channel, the generator is: Since the decay rate γ 3 (t) < 0 in equation ( 25), the total channel is non-Markovian more so this channel is eternally non-Markovian [33,34].

Non-Markovian channel addition (nM + nM = M)
To design channels for the non-Markovian channel addition, we need to consider a special case of the Pauli channel in equation (11), i.e. the depolarizing channel [20].We can obtain the depolarizing channel from equation ( 11) by parameterizing the probabilities p α (t) as follows: where 0 ≤ p(t) ≤ 1 for all times t ≥ 0. Now from the fact that 0 = 1 we see that p(0) = 0.The decay rates for the depolarizing channel can now be calculated using equation (19), for k = 1, 2, 3. From equation ( 27), we can tighten the bounds on p(t), i.e. 0 ≤ p(t) < 1.Now the form of the decay rate in equation ( 27) tells us that for the channel to be Markovian, the function p(t) should satisfy ṗ(t) ≥ 0 for all t ≥ 0 and for the channel to be non-Markovian there should exist some time t ≥ 0 such that ṗ(t ) ≤ 0. Now let us consider the following two individual channels: (1) The decay rates for both the individual channels in equation ( 28) are: for k = 1,2,3.Taking a convex combination of the individual channels in equation ( 28), we obtain the total channel: (T) where η ∈ [0, 1] and w(t) = ηq(t) + (1η)r(t).Now the decay rates for the total channel are: It is clear from equation ( 31) that it is possible for there to exist times t and t such that q(t ) < 0 and ṙ(t ) < 0 and w(t) ≥ 0 for all times t.This tells us that it is possible to use the depolarizing channel in equation ( 26) to show that we can convexly mix two non-Markovian channels to get a total Markovian channel.As a consequence of our calculations we can also see that it is impossible to have mixtures of two Markovian depolarising channels that yield a Non-Markovian depolarising channel.To see this consider the decay rates in equation ( 29), for the channels (1)  t and (2) t to be markovian we require that q(t) ≥ 0 and ṙ(t) ≥ 0 for all times t ≥ 0. However if this is true then the decay rates for total channel (T)   t in equation (31) will always be positive since, η q(t) + (1η)ṙ(t) ≥ 0 for all t ≥ 0. This tells us that mixing two Markovian depolarising channels always yields a Markovian depolarising channel.This result is useful as depolarising channels are the most commonly used noise models when studying the robustness of quantum algorithms [52,53] to noise.
The goal now would be to pick functions q(t) and r(t) such that (1)  t and (2)  t are non-Markovian and the total channel (T) t is Markovian.We observe that if we parameterize q(t) and r(t) as: Furthermore, setting η = 1 2 in equation (30), then using the bounds on q(t) and r(t), we get b(t) ≤ a(t) < 1b(t).Taking into consideration all the constraints outlined above, we choose the functions a(t) and b(t) as follows: Refer to Fig. 1 for the plots of the functions given in equation ( 33) above.From equation (29), we see that the decay rates for the individual channels (1)  t and (2)  t are non-Markovian.This is because for some time interval q < 0 and for some other interval of Figure 1 This is a plot of the functions q(t) and r(t) for the time interval t ∈ [0, 9].Both these functions satisfy the constraints that 0 ≤ q(t), r(t) < 1 and the functions are numerically close enough to zero that we can write: q(0) ≈ 0 and r(0) ≈ 0 Figure 2 This plot shows the derivatives of q(t) and r(t), we see that for some times both q(t) and ṙ(t) are negative which will lead to the decay rates in the generators of the individual channels to become negative.This implies that the individual channels are both non-Markovian Figure 3 Plot showing the function w(t) which parameterizes the total channel (T)  t , one can easily see that w(t) satisfies all the constraints above.This plot also shows ẇ(t), and it is clear that ẇ(t) ≥ 0 for all times t, which tells us that the decay rates of the total channel are all greater than or equal to zero which implies that the total channel (T)  t is Markovian time ṙ < 0 making the decay rates negative, leading to the channels being non-Markovian.Figure 2 shows the plots of q(t) and ṙ(t), which shows that both are negative for some times.The total channel (T)   t is parameterized by the function w(t) and it is clear from Fig. 3 that ẇ(t) ≥ 0 for all times t ≥ 0, so the total channel (T) t is Markovian.Refer to Appendix A for more intuition about how the functions q(t) and r(t) were chosen.Hence we have found an example of the convex sum of two non-Markovian channels (1)  t and (2) t yielding a Markovian total channel (T) t .

Simulation of quantum channels on a NISQ device
The simulation of the channels from the previous section on a NISQ device requires a twofold method.First, one needs to construct quantum circuits that implement the action of the quantum channel on a single qubit.Then we should perform a quantum pro-cess tomography (QPT) to reconstruct the quantum channel.The tomography will yield a non-physical quantum channel.Hence convex optimization must be used to construct the closest possible physical quantum channel to the channel obtained from the QPT.The following sub-sections shall detail the method above.

Quantum circuits for simulation of the channels on a NISQ device
To simulate the quantum channels in Sect.3, we must construct a quantum circuit that implements the channel.Usually, this is done by naively using the Stinespring representation of the channel [43], but this naive application of the Stinespring dilation may lead to a complicated quantum circuit unsuitable for the NISQ setting.To solve this problem we use the Stinespring representation in a way that takes into account the properties of the NISQ device, when we construct circuits for the channels.To design quantum circuits for our channels from Sect. 2 that are appropriate for the NISQ setting, we take advantage of the fact that in both cases (i.e.M + M = nM and nM + nM = M), the channels that were designed implement a Pauli operator with some probability.We can easily design circuits that can implement the channels by using ancilla qubits.This approach is inspired by [20] and allows us to construct NISQ-appropriate quantum circuits.

Circuits for Markovian channel addition
To implement the Markovian channel (1)  t with probability p(t) we use the circuit shown in Fig. 4 below.In the circuit the gate R y (θ ) is a rotation gate and its matrix representation is, where θ is the angle of rotation.The circuit is understood as follows it will apply the gate σ 1 = X to the input state ρ in with probability | sin(θ /2)| 2 and it will leave the state unchanged with probability | cos(θ /2)| 2 (i.e. when the ancilla is in the state |0 do not change the state ρ in and when the ancilla is in state |1 apply the X gate to the input state ρ in ).
From this, we can get the value of the angles in terms of the probability p(t) as: Since the angle θ is written in terms of p(t) the circuit in Fig. 4 will leave ρ in unchanged with probability p(t) and it will apply σ 1 to ρ in with probability 1p(t).Hence the circuit implements the channel (1)  t .Similarly, we can design a circuit that implements (2)   t by using the same ideas used to design the previous circuit.Figure 5 shows the quantum circuit that implements the channel (2)  t where the Y gate is applied to the input state with probability p(t).The angle θ is obtained using equation (35).Now that Figure 4 Quantum circuit implementing the Markovian channel (1)  t for probability p(t) Figure 5 Quantum circuit implementing the Markovian channel (2)  t for probability p(t) Figure 6 Quantum circuit implementing the total non-Markovian channel (T) t for probability p(t) we have designed circuits that implement the individual Markovian channels, we design a circuit that implements the total non-Markovian channel (T) t .In designing this circuit, we keep in mind that we want to leave the input state ρ in unchanged with probability p(t) and apply the noise operators σ 1 and σ 2 to the input state with probability 1-p(t) 2 .Refer to Fig. 6 for the quantum circuit that implements the total channel.This circuit implements the channel by first preparing the two ancilla qubits in the state: cos(θ 1 /2)|00 -sin(θ 1 /2) sin(θ 2 /2)|01 -sin(θ 1 /2) cos(θ 2 /2)|11 now using the ancilla qubits the circuit will apply X to the input state if the ancillae are in the state |01 , it will apply Y to the input state when the ancillae are in the state |11 (NB: the circuit applies XZ to the input state which is the same as applying Y up to a phase).The circuit leaves the input state unchanged when the ancillas are in the state |00 .From the state of the ancillae, we can determine the angles θ 1 and θ 2 : Now that we have the angles in terms of p(t), we have found the circuit that implements the channel (T) t .

Circuits for non-Markovian channel addition
To implement the channels that we have designed for the non-Markovian channel addition, we need a circuit that can implement the depolarizing channel that is parameterized by some function p(t) (Refer to equation ( 26)).We use the circuit simulating the depolarizing channel in [20] to do this.Refer to Fig. 7 for the circuit that simulates the depolarising channel.The three ancillea qubits are prepared in the state (cos(θ /2) |0 + sin(θ /2) |1 ) ⊗3 are used as control qubits for the controlled-X, controlled-Y and controlled-Z gates, so that the gates X, Y and Z are all applied to the input state ρ in with probability sin 2 (θ /2).To use this circuit to implement the depolarising channel in equation ( 26) the rotation angle θ needs to be in terms of the probability p(t) so that each gate is applied with probability p(t)/4.We observe that applying X and then Z to ρ in but not Y , is equivalent to just applying Y to ρ in and so on.The resulting equation that relates the probability p(t) to θ is, Solving this equation for θ in terms of p(t) we get, which will be the parameter in the R y gate in then circuit in Fig. 7.

Tomographic reconstruction of the channel
It is known that a quantum channel t has a Kraus representation [54]: where Kα are the Kraus operators that satisfy α K † α Kα = 1.In this work, we will consider the case of a single qubit channel then the Kraus operators are 2 × 2 matrices.If we choose a complete basis for the Kraus operators of a single qubit channel as {σ 0 = 1, σ 1 , σ 2 , σ 3 }, where σ i are the usual Pauli matrices.Then we can expand the Kraus operators in terms of this basis to get the process matrix representation of the quantum channel for a single qubit: Here χ mn is a positive and Hermitian 4 × 4 matrix called the process matrix and shall be determined using a quantum process tomography [55,56].Now if we know the process matrix, then we have a complete description of the channel t .
To determine the elements of the process matrix, we need to choose a complete set of input states.We choose the states The states from the set D form a complete set as the projectors constructed from each ket vector in this set can be used to construct the density operator of any physical single qubit state.The input states are sent through the channel, t .We can prepare the initial state of the qubit as each input state.Using quantum state tomography, we reconstruct the state after each input state is passed through the channel [57].Then the formulas from [55,56] are used to construct the χ matrix which allows us to reconstruct the channel t .
To perform a QPT on a quantum computer, we use the circuits from Sects.3.1.1and 3.1.2,and then the system qubit is prepared in one of the input states above and passed through the channel.We then perform the state tomography on the system qubit for each input state and get the corresponding counts.Using the counts obtained from the tomographic circuits, we can construct the process matrix [55].
We use the counts obtained from the tomographic circuits to construct an initial process matrix denoted by, χ in .The matrix χ in will not be positive and Hermitian.This is because we can only make a finite number of measurements on the system qubit.To correct this, we shall solve a convex optimization problem to find the closest possible process matrix, denoted χ c , to the ideal process matrix, which we denote as χ id [44].We shall use the matrix χ in as an initial condition when solving the optimization problem.To construct the optimization problem, we first need to parameterize χ c in terms of parameters that we can optimize.We also define χ c to be Hermitian and positive semi-definite for all parameter values.We now parameterise χ c as, where T = T(x 1 , x 2 , . . ., x 16 ) is a 4× 4 triangular matrix that is a function of 16 real variables x 1 , x 2 , . . ., x 16 and is shown below in matrix form: It is evident from this parameterization that χ c is positive semi-definite.Consider some arbitrary four-dimensional vector |ψ then, It is also evident that χ c is Hermitian.To find χ c with convex optimization, we need to define an objective function that will be minimized with respect to constraints.We define the objective function by the squared difference between the theoretical and experimental probability distributions for each of the counts obtained from the process tomography.
The following projective measurement operators are defined from the set D above, Next, we consider the input state ρ i = |ψ i ψ i | where |ψ i ∈ D. Then the theoretical probability of being in the state |ψ j after the application of the channel t to the initial state |ψ i , denoted p the ij , is, The experimentally obtained probability p exp ij of being the initial state ρ i and measuring in the state |ψ j is, where n ij is the counts obtained from the circuit with input state |ψ i and output state |ψ j and N being the total number of counts.Now we can define the objective function as, We should also consider that the channel t should be trace-preserving.This will be one of the constraints we define for our optimization problem.From equation (40) we see that for t to be trace-preserving we must have that, 4 m,n=1 We can weaken this constraint by just requiring that t be trace non-increasing, this leads to, 4 m,n=1 This constraint can be written as a positive semi-definite constraint i.e.
We can now state the optimization problem that we want to solve that will yield a positive semi-definite and Hermitian matrix χ c that is the closest to χ id , where we obtain an initial set of parameters (x (in) 1 , . . ., x (in) 16 ) from χ in .The solution to this problem will yield the optimal values for x 1 , . . ., x 16 so that χ c is as close to χ id as possible.
The following section will use the matrix χ c for the characterization.

Characterization of the channel as Markovian or non-Markovian
From the previous section, we know that the state of the system at some time t ≥ 0 is given in terms of the quantum channel t as, ρ(t) = t ρ(0), where ρ(0) is the initial state of the system at t = 0. We shall use the CP divisibility criteria to characterize the channel t as Markovian or non-Markovian.We start by writing the dynamical map t in the following way: where V t,s is called the intermediate map (IM) from time s to t.The maps V t,s form a family of two-parameter propagators.We say that t is CP divisible if V t,s is completely positive (CP) for all t ≥ s ≥ 0. The goal is to check that the map V t,s is CP.We will use the techniques in [38].From Sect. 3, we know that we can express a quantum channel in the following representation: where χ mn is a positive and Hermitian 4 × 4 matrix called the chi matrix, σ 0 = 1 and σ i for i = 1, 2, 3 are the Pauli matrices.Knowing the chi matrices for two different time durations, s and t, we can check whether the map V t,s that evolves the system from time s to time t is completely positive.To check this, we make use of the transfer matrices F(s) and F(t) for the maps s and t , respectively.The transfer matrix F(t) of a map t is a concrete matrix representation of the map in a given orthonormal basis [58].The transfer-matrix approach is useful as it allows us to represent the density matrix ρ as a stacked vector |ρ , now the evolution of the vector |ρ(0) can be written as: |ρ(t) = F(t)|ρ(0) which is nothing more than the action of the transfer matrix on the stacked vector [38].The elements of a transfer matrix F(t) are given explicitly as, where {G α } are a set of orthonormal operators with respect to the Hilbert-Schmidt inner product [58].We choose the set {G α } to be the standard matrix basis of M 2 (C), i.e.
, where {|0 , |1 } are the standard computational basis vectors of the single qubit.If we know the χ matrix for some time t, we can use this to calculate t G β and hence calculate F(t).Now using equation (58) and writing it in terms of transfer matrices, we arrive at F(t) = F(t, s)F(s).This tells us that if we have the transfer matrix for two times s and t, we can get the transfer matrix of the intermediate map, i.e.
F(t, s) = F(t)F -1 (s).For a given transfer matrix F(t) we can obtain the Choi matrix W (t).For a single qubit, this can be written as: This is derived by applying t to a single qubit of the maximally entangled state By the Choi-Jamiolkowski isomorphism, a dynamical map t is completely positive if and only if the corresponding Choi matrix of the map is positive [59,60].This tells us that the map t is CP if all the eigenvalues of the Choi matrix W (t) are non-negative.So for two times s and t such that t ≥ s ≥ 0, if the eigenvalues of the Choi matrix W (t, s), of the intermediate map V t,s , are non-negative then the intermediate map is completely positive and by our definition of CP divisibility, this tells us that the dynamical map t = V t,s s is Markovian.Any deviation from this leads to non-Markovian dynamics.

Results and discussion
This research aimed to demonstrate that more complex open quantum systems can be efficiently simulated on NISQ devices.In particular, we aimed to show that we can simulate convex mixtures of Markovian and Non-Markovian quantum channels on NISQ devices.In this section, we present quantum circuits based on the circuits in Sect.4.2 that account for the topology of the NISQ devices they will be run on.We also present the results of the simulations and evaluate the quality of the results using three metrics: fidelity of the process matrix (discussed in Sect.6.2), trace norm of the Choi matrix (discussed in Sect.6.3) and the minimum eigenvalue of W (t, s) in Sect.6.4.

Implementation of quantum circuits
For our experiments, the circuits are executed using the NISQ device provided by IBM Quantum Experience (IBMQE).Figure 8 gives an overview of the steps used to construct the quantum circuits, perform the tomography, and reconstruct the channel.The sevenqubit device called ibm_perth was used.Figure 9 gives a diagrammatic representation of the qubits in the quantum computer as well as the connections between the qubits.The The topology of the NISQ computer used for the experiments.We made use of the seven qubit device ibm perth for both cases M + M = nM and nM + nM = M.The circles represent the qubits, and the lines represent the connectivity between the qubits.For example, qubits three and five are connected, meaning one can perform a CNOT between them directly.The qubits that are coloured in red are used for the experiments, while the blue qubits are not.For the case M + M = nM, we use qubits one to three; for the case nM + nM = M, we use qubits zero to three Figure 10 Quantum circuit implementing the total non-Markovian channel (T)  t for probability p(t) with the addition of a swap gate to mitigate controlled operations between unconnected qubits connectivity of the qubits plays an essential role in constructing the circuits for the quantum channels.This is because, on all NISQ devices, the quality of the results depends on the number of controlled not gates (CNOTs) used.One should note that this is not the only factor but is most important to us in our simulations.The NISQ devices can only implement a small number of CNOTs before the qubit interactions, noise, and decoherence affect the results.Moreover, the implementation of CNOTs between qubits that are not directly connected requires a SWAP gate to swap the state of the qubits so that the CNOTs can be implemented between directly connected qubits, and then the states are swapped back, which is done automatically by the quantum computer.This leads to problems as SWAP gates are equivalent to three CNOTs, significantly increasing the number of CNOTs.
We have mitigated this issue by looking at the circuits designed in Sect. 3 and making some additions whenever we implement controlled operations between qubits that are not directly connected.Rather than letting the quantum computer perform SWAP operations, we manually add the SWAP gate to minimise the number of SWAPs needed, minimising the number of CNOTs used.
For the M + M = nM case, we use qubits one to three in Fig. 9. Qubit two is the system qubit, and qubits one and three correspond to the environment.For the Markovian channels (1)  t and (2)  t there are no controlled operations between unconnected qubits.However, for the non-Markovian total channel (T)  t , we see that we require a controlled operation between qubits three and two.
Figure 10 shows the new quantum circuit for the non-Markovian channel (T) t with the additional SWAP gate.
For the case nM + nM = M, we use qubits one as the system and qubits zero, two and three as the environment.Since there will be no controlled operations between unconnected qubits for the simulation of the channels in the non-Markovian channel addition case, we do not need to modify the circuits from Sect. 3 with additional SWAP gates.
Once the circuits are executed, a process tomography is performed on the channels for each time step t = 0.1s in the interval t ∈ [0, 3.9] for the Markovian channel addition (M + M = nM), and the interval t ∈ [0, 8.9] for the non-Markovian channel addition (nM + nM = M).The counts obtained from the results of these circuits are used to generate the process matrices χ in (t) for each time step.By adding Gaussian fluctuations, 100 process matrices for each time step were generated.Now, using convex optimisation, as outlined in Sect.4, the closest possible process matrices χ c (t) to the ideal process matrices χ id (t) were obtained.
We then pick an intermediate time s such that t ≥ s ≥ 0 and for each pair (t, s), 10,000 process matrices are obtained.These are used in the analysis to compute the three metrics.After optimisation, the first metric, the fidelity for the process matrices, measures how close these process matrices are to the ideal.The second metric, the trace norm of the Choi matrix of the channels at each time step, allows us to check if the simulated channels are trace-preserving.The last metric we considered is the minimum eigenvalue of the Choi matrix of the intermediate map W (t, s), which allows us to characterise the channel as either Markovian or Non-Markovian.These three metrics allow us to verify that the channels simulated are physical and satisfy the case of convexly mixing two Markovian channels to yield a non-Markovian channel and vice versa.

Fidelity of the process matrices
In both cases, we compute the process fidelity of the process matrix for each time t after optimisation, using the following formula [38,61]: We note that F p ∈ [0, 1] when F p = 1 this tells us that the process matrix is the same as the ideal, i.e. χ = χ id and when F p = 0 the process matrix is far from the ideal process matrix χ id .We shall compute the process fidelities for the channels (1)  t , (2)  t , (T) t for both the cases of Markovian channel addition and non-Markovian channel and plot them for each time step.
In Fig. 11, we plot the process fidelities for the Markovian channel addition (M + M = nM).We see that in Fig. 11 (a) and (b) the fidelities for the channel (1)  t and (2)  t are very close to 1.This tells us that our χ matrices are very close to the ideal case after MLE.In Fig. 11 (c), we see that the fidelities for the total channel (T) t while not as high as the other two channels are still relatively good enough for our experiment.It should be noted that the fidelity is lower in this case due to decoherence and dissipation in the quantum computer.
We plot the process fidelities for the channels in Fig. 12.We see that for Fig. 12 (a) and Fig. 12 (b) the fidelities for the individual channels are high and have a value of one for large parts of the time interval, indicating that the quality of the χ matrices is good.In Fig. 12 (c), the fidelity of the total channel is very good, although at time t = 1.8 s, the fidelity is low.This tells us that the χ matrices for the total channel will be close to ideal.For most of Figure 11 The process fidelities for the χ matrices of the implemented channels after MLE for times from 0 to 3.9 seconds with time step 0.1 seconds.(a) and (b) show the fidelities for the Markovian channels (1)  t and (2) t respectively.(c) shows the fidelity for the total non-Markovian channel (T) t Figure 12 The process fidelities for the χ matrices of the implemented channels after MLE for times from 0 to 8.9 seconds with time step 0.1 seconds.(a) and (b) show the fidelities for the non-Markovian channels (1)   t and (2)  t respectively, the fidelities for these channels are good as for (1)  t the fidelity is 1 for a large part of the time interval, for (2)  t the fidelities do fluctuate but are also very good.(c) shows the fidelity for the total Markovian channel (T) t the interval but at t = 1.8 s, the χ matrix will be far from the ideal matrix.These fidelities are high enough such that the accuracy for χ matrices used in our analysis is high.

Trace norm of Choi matrix
The trace norm for the Choi matrix at each time step t can be used to check if the simulated channel was trace-preserving.We calculate the Choi matrix W (t) for each time step by using equations ( 60)-(61) from Sect. 4. Once we have found the Choi matrices W (t) we find the trace norm, We know that for any quantum channel t to be trace preserving it must satisfy, The trace norm for the channels in both cases, M + M = nM and nM + nM = M, are calculated and plotted in Fig. 13 and Fig. 14, respectively.We see that in Fig. 13 that the trace norm of W (t) for the Markovian channel addition (M + M = nM) is a good enough approximation of 1 for the simulated channel to be considered trace-preserving.It should be noted that device noise and finite sampling are to blame for the deviations from 1.In Fig. 14, we see that for the case of the non-Markovian channel addition (nM + nM = M), the trace norm of W (t) is a better approximation of 1 as they cover the theoretical curve.  (1 t , (2)  t and (T) t respectively, for the Markovian channel addition (M + M = nM).For each case, the trace norm is a good enough approximation of 1 to state that the simulated channels are trace-preserving.The deviation and slightly long error bars result from device noise and decoherence Figure 14 (a)-(c) We plot the trace norm of the Choi matrices of the channels (1)  t , (2)  t and (T) t respectively, for the non-Markovian channel addition (nM + nM = M).The trace norm for these simulations covers the theoretical value and is a good enough approximation to say that the simulated channels were trace-preserving The same reasons for the deviations in the Markovian channel addition also cause the deviations seen here.In both cases, the trace norm is a good enough approximation to consider the simulated channels trace-preserving.

Minimum eigenvalue of W(t, s)
Now, we use the characterisation method outlined in Sect. 4 to classify the channels.By plotting the minimum eigenvalues of each Choi Matrix W (t, s) of the IM, we can easily classify the corresponding channel as Markovian if all the minimum eigenvalues are nonnegative or Non-Markovian if any minimum eigenvalue is negative.This is according to the CP-divisibility criteria.First, we look at the minimum eigenvalues for the Markovian channel (1)  t (s = 0.5) in Fig. 15 (a).We observe that the minimum eigenvalues are all negative.This differs from the theoretical minimum eigenvalues, which are all zero.This discrepancy is expected since it is difficult to measure zero in any experiment due to standard deviation.This discrepancy can also be seen in Fig. 15 (b) between the theoretical and experimental minimum eigenvalues for the Markovian channel (2)  t (s = 0.5).In both these cases, the minimum eigenvalues are very small negative values and can be accepted as approximations of zero.The minimum eigenvalues for the Non-Markovian channel (T) t (s = 0.5), as shown in Fig. 15 (c), show definitively that this channel is Non-Markovian as the minimum eigenvalues, within their standard deviation, are all negative.While the minimum eigenvalues do fluctuate due to noise, the eigenvalues follow the behaviour of the theoretical curve.15 The plot (a) shows the minimum eigenvalues of the IM for the Markovian channel (1)  t , it is clear that the experimental results are a good enough approximation of the zero eigenvalue.(b) Shows the minimum eigenvalues of the IM for the channel (2)  t .The experimental results in this case are better than for the first channel as the standard deviation is small for all the points.(c) shows the minimum eigenvalues of the IM for the total non-Markovian channel (T)  t , the points cover the theoretical curve, hence this channel is definitively non-Markovian by the CP divisibility criteria [29] Figure 16 The plot (a) shows the minimum eigenvalues of the IM for the non-Markovian channel (1)  t , it is clear that the experimental results are in agreement with the theoretical minimum eigenvalues of the IM.(b) Shows the minimum eigenvalues of the IM for the channel (2)  t , the experimental results in this case deviate from the theoretical curve but this channel is still non-Markovian by CP divisibility.(c) shows the minimum eigenvalues of the IM for the total Markovian channel (T)  t , the points cover the theoretical curve and the initial points, although negative, are still a good enough approximation of zero, this is a problem face when the minimum eigenvalue is zero.Hence, this channel is Markovian by the CP divisibility criteria [29] Next, we look at the minimum eigenvalues for the Non-Markovian channels (1)  t and (2) t (s = 3), shown in Fig. 16 (a) and (b).We see that the minimum eigenvalue for (1)   t is negative in the time interval t = 51 s to t = 5.8 s and for (2)  t the minimum eigenvalue is negative in the interval t = 3 s to t = 5.4 s.This shows that these channels are Non-Markovian.We can also see that the minimum eigenvalues follow the theoretical curves.While there are slight deviations from the curve, these deviations can be attributed to noise from the experiment.Lastly, in Fig. 16 (c) we look at the minimum eigenvalues for the Markovian channel (T)  t (s = 3).We see that the minimum eigenvalues between t = 0 and t = 6 are negative when the theoretical eigenvalues are zero.Once again, this discrepancy can be attributed to the difficulty of measuring zero in any experiment due to standard deviation.The minimum eigenvalues in this time interval are very small negative values and can be accepted as approximations of zero.Therefore, the channel can be classified as Markovian.
In each case, we see that the minimum eigenvalues follow the theoretical curves and that the channels can be successfully classified as either Markovian or Non-Markovian using the CP-divisibility criteria.

Conclusion
We have successfully performed simulations of mixtures of two Markovian single qubit Pauli channels that give rise to a non-Markovian single qubit Pauli channel and two non-Markovian single qubit depolarising channels that give rise to a Markovian single qubit depolarising channel.The success of these simulations has been verified in three ways.Firstly, we have shown by the CP-divisibility criteria that the channel reconstructed from the simulations of the M + M = nM mixtures and the nM + nM = M mixtures are non-Markovian and Markovian, respectively.We have also shown that the reconstructed channels have high fidelity to the theoretical channel.Lastly, we have shown that the Choi matrices have a required trace norm of one, indicating that the simulated channels are physical.This demonstrates the effectiveness of designing circuits with NISQ device topology in mind and provides more accurate results, as the quantities we use to benchmark our simulation are all excellent enough approximations of their theoretical values.This demonstrates that using the least squared objective function [44] and convex optimization is a better choice for reconstructing our channel than MLE.The success of these experiments shows that a NISQ computer can be used to simulate more complex open quantum systems.However, future work is needed to verify that the strategies used here generalise to more general quantum channels.In this work, we consider only single qubit Pauli channels.The simulation of many qubit and non-Pauli channels is not considered here but will be examined in future work.We would like to see if our effective strategies in simulating single qubit channels can aid in simulating single qubit non-Pauli channels and many qubit Pauli and non-Pauli channels.Future work could look at applying the developed experimental pipeline to different channels since the pipeline is not specific to the channels used in this work.Other future work includes testing and comparing different objective functions for optimization and changing the constraints.One could also look into using a semi-definite program for the optimization part of the reconstruction.

Figure 7
Figure 7 circuit implementing the depolarizing channel for a single system qubit, for the probability p(t).The angle θ is determined by the formula θ(t) = 1 2 arccos(1 -2p(t))

Figure 8 AFigure 9
Figure 8 A flow chart summarising the experimental procedure used in the simulation and reconstruction of the channels for both cases M + M = nM and nM + nM = M

Figure 13 (
Figure13(a)-(c) We plot the trace norm of the Choi matrices of the channels(1)  t ,(2)  t and(T)

Figure
Figure15 The plot (a) shows the minimum eigenvalues of the IM for the Markovian channel(1)  t , it is clear that the experimental results are a good enough approximation of the zero eigenvalue.(b) Shows the minimum eigenvalues of the IM for the channel(2)  t .The experimental results in this case are better than for the first channel as the standard deviation is small for all the points.(c) shows the minimum eigenvalues of the IM for the total non-Markovian channel(T)  t , the points cover the theoretical curve, hence this channel is definitively non-Markovian by the CP divisibility criteria[29]