Simulation of open quantum systems via low-depth convex unitary evolutions

Simulating physical systems on quantum devices is one of the most promising applications of quantum technology. Current quantum approaches to simulating open quantum systems are still practically challenging on NISQ-era devices, because they typically require ancilla qubits and extensive controlled sequences. In this work, we propose a hybrid quantum-classical approach for simulating a class of open system dynamics called random-unitary channels. These channels naturally decompose into a series of convex unitary evolutions, which can then be efficiently sampled and run as independent circuits. The method does not require deep ancilla frameworks and thus can be implemented with lower noise costs. We implement simulations of open quantum systems up to dozens of qubits and with large channel ranks.


I. INTRODUCTION
System-environment interactions play an important role in the dynamics of quantum systems, giving rise to phenomena such as dissipation, decoherence, relaxation, and particle or energy transfer processes [1,2].Ideal simulations of open quantum systems on classical computers suffer from exponentially scaling memory and runtime requirements, leaving complex physical systems largely out of reach without major approximations [3].Using quantum computers, several proposed techniques for simulating open systems exist, with a common approach being the encoding of nonunitary dynamics in a larger, dilated Hilbert space [4].Numerous techniques make use of Stinespring dilation [5,6], Sz.-Nagy dilation [7,8], and linear combination of unitaries [9][10][11].In practice, this dilation leads to long controlled gate sequences and large ancilla frameworks, both of which are practically at odds with NISQ-oriented applications.
Among a wide variety of open system processes, random-unitary channels form an important and practical subset.Also known as mixed-unitary, these channels admit an operator-sum decomposition [12,13] of the form where p i are probabilities and U i are unitary operators [14].Random-unitary channels naturally represent many probabilistic processes, including the well-known Pauli channels, and have found widespread utility across quantum information science.In fact, general Markovian processes can be mapped onto effective random-unitary channels through Pauli twirling [15,16] and randomized compiling techniques [17].They are particularly useful for noise modeling of quantum devices [18,19], as well as within error correction [20] and mitigation [21] schemes.
They have also been implemented in quantum simulations of thermal relaxation [22,23].As quantum computers continue to grow in size, efficient simulation of random-unitary channels thus becomes increasingly imperative.
To address this need, we propose a hybrid quantumclassical method for simulating random-unitary channels.Our approach involves classically sampling a channel's probability distribution {p i }, and then running a series of corresponding circuits on a quantum computer.Compared to ancilla-based methods, this reduces both the number of qubits and the depth of circuits required, thereby widening the class of simulations feasible on nearterm quantum computers.By outlining a practical and efficient algorithm, we seek to enhance both the understanding and utility of random-unitary channels.

II. RANDOM-UNITARY ESTIMATION
We present a hybrid quantum-classical method of simulating random-unitary channels for open quantum systems.For a state passing through a random-unitary channel of the form in Eq. ( 1), ρ → ρ = E(ρ), the expectation value of an observable Ô is Typically, ancilla-based methods attempt to faithfully produce ρ, which can then be accessed directly via measurement.
In our approach, we obtain the above expectation via an estimator of the random-unitary channel, following the formalism of Arrasmith et al. [24].Given the probability distribution {p i } and N total shots, we generate the following multinomial distribution: (3) arXiv:2307.14325v3[quant-ph] 6 Sep 2024 } FIG. 1. Schematic for simulating random-unitary channels via low-depth convex unitary evolutions.Open system dynamics are described mathematically by their operator-sum representation F, where Ki are Kraus operators.Processes which preserve the maximally mixed state are called unital channels, of which random-unitary channels E are the subset described by Eq. ( 1).An observable Ô of states evolved via random-unitary channels can be simulated via the low-depth, parallel quantum circuits depicted.Here, each circuit corresponding to the unitary Ui is simulated with E[si] shots, generated from the multinomial distribution Mult(N, {pi}).Each measurement is performed with respect to the observable Ô, and the results are then aggregated to compute the unbiased estimator Ê.
Here, the ŝi outcomes correspond to quantum circuits measuring the operator U i ÔU † i .Then, the estimator of the expectation values has the form where E[ŝ i ] = p i N .This matches the form of Eq. (2b) and is an unbiased estimator.Figure 1 represents this scheme pictorially.Because the probability distribution {p i } naturally sums to one, the variance of the estimator is simply the variance of the observable Ô divided by the number of shots N .Notably, it follows that the complexity and specific parameters {p i } of the randomunitary channel have no effect on the precision of the estimator.
We note that a similar heuristic approach was tried for low-dimensional systems, but associated with exponential scaling problems [22] and was limited to singlequbit applications [23].Crucially, this work demonstrates the feasibility and scalability of this method even for exponentially complex probability distributions {p i }, as long as this distribution can be sampled efficiently.We demonstrate these ideas for simple open systems which, even in lieu of error mitigation, show constant precision over exponentially scaling problems.

III. DEMONSTRATIONS OF SCALABLE RANDOM-UNITARY SIMULATIONS
To highlight our approach, we present two cases, both of which lead to exponentially scaling problems.We first look at examples of the depolarizing channel, a Pauli channel with exponentially many operators with respect to the number of qubits.We then look at a discretized time-evolution process, whose operator complexity grows exponentially in time.In the first instance, due to the simplicity of the method and the lack of ancilla costs required, we are able to easily demonstrate this approach on superconducting transmon qubit devices.The second example allows for a straightforward demonstration on quantum devices as well, although it can also be classically simulated via exact storage of the two-qubit density matrix.

Let P
(n) i denote an n-qubit Pauli string, i.e.P (n) i ∈ {I, X, Y, Z} ⊗n .Then the n-qubit depolarizing channel [25] can be written as For n ∈ [1,27] (with ibmq montreal having 27 qubits), we simulate the depolarizing channel on ibmq montreal with strength parameter p = 0.5 and an initial state of we also demonstrate an ancillabased approach, using a simple linear combination of unitaries.While many such methods exist, we compare simply against the most common example [9], as generally all of these methods require ancilla circuits.The resource requirements of the two methods are compared in Table I.We measure the set of single-Z observables ⟨Z (i) ⟩, e.g., Z ⊗ I for two qubits or I ⊗ Z ⊗ I for three, and compare to the following analytic result: We plot the mean squared error, MSE = 1 2 , against the number of qubits for both the ancilla and low-depth methods.Despite the depolarizing channel's exponentially large sample space of 4 n operators, the probability distribution {p i } can be sampled efficiently.Because these operations belong to the Clifford group [26], this channel can be classically simulated in polynomial time according to the Gottesman-Knill theorem [27].Thus, we also include an instance of a simple non-Clifford state, where we prepare the qubits pairwise in the entangled state √ 2 |11⟩ and simulate the depolarizing channel for even values of n.We measure the set of pairwise observables X(2i) X(2i+1) and compare to the following analytic result: We plot the mean squared error against the number of simulated qubits using our low-depth hybrid method.
All results are combined in Fig. 2. As shown, our hybrid sampling method clearly outperforms the traditional ancilla-based approach.It also demonstrates consistently reliable results up to the maximum number of qubits.
In contrast, circuits to run the multiplexing operations for the ancilla-based method become infeasible after 3 qubits.
Despite the high numbers of qubits and the prevalence of errors, the simulations here still manage to faithfully capture the simulated state dynamics.In Fig. 3 pected bit-flip errors across each qubit skew the results away from the all-zero state, yet not in a way that corrupts the overall distribution.We model the overall distribution as the depolarizing channel coupled with bit-flip errors of p X = 4.7% per qubit.

B. Noisy Hamiltonian evolution
To look at the performance under a different exponentially growing scheme, we simulate the two-qubit transverse-field Ising model (TFIM) under a noise channel which occurs at each discretized step.
The TFIM Hamiltonian is z − h σ (1)  x + σ (2)   x , where J is the exchange interaction parameter and h quantifies the strength of the transverse magnetic field.For small time steps ∆t, we evolve this system as the Hamiltonian evolution operator composed with the random-unitary channel S. The resulting composed operator is itself a random-unitary channel, denoted as S.
We generate each time step recursively, by sampling a distribution as in Eq. ( 4): We simulate these dynamics on the ibmq kolkata device, with a depolarizing strength of p = 0.05 per time step, TFIM parameters J = 1 and h = 1, and the starting state |00⟩.Figure 4 shows the populations over time, both in the standard computational basis and in the eigenbasis of the Hamiltonian.Absent any noise, the eigenstates would be stabilized at constant populations.Fitting exponential decay functions to the eigenstate population curves yields an estimated decoherence time of T 1 = 4.84 ± 0.46 (a.u.), encompassing the predicted time of T 1 = 5.25 (a.u.).In this example, the possible number of circuits to sample from grows exponentially with the number of time steps.For our simulations, after 25 time steps there are 16 25 circuit permutations to consider, a seemingly intractable sample space.However, due to the constant norm of the probability distribution vector, the measurement statistics are agnostic to the complexity of the sample space.This enables our results to accurately capture the system-environment dynamics while using only 10 3 shots at each point in time.

IV. DISCUSSION AND CONCLUSION
These results demonstrate that for a certain class of quantum channels, the proposed hybrid quantumclassical approach offers significant efficiency advantages.By stochastically generating low-depth quantum circuits and then aggregating the measurement results, we enable the simulation of exponentially complex processes.
A key result in our work is that the variance of our estimator is unchanged by the specifics of the probability distribution.Because random-unitary channels consist of a convex combination of unitary operations, the set of sampling probabilities {p i } always sums to unity.While we chose the depolarizing channel to demonstrate scalability, any random-unitary channel can be simulated via this method, assuming only that its probability distribution can be sampled efficiently.This includes non uniform noise channels, in which each unitary operator U i has a distinct probability p i .
While our method leads to efficiency gains important for noisy quantum device applications, one limitation is that it does not create a coherent quantum state as an output.This means that any observables obtained via further evolutions of the state must also be processed through an estimator.For instance, if one wishes to use a mixed state as input to some broader algorithm, then any desired observables at the end of this algorithm must themselves be sampled using an estimator of the form in Eq. ( 4).However, if the channel's operator representation has exponentially many terms, implementation via ancilla-based methods will be significantly challenging as well.This problem is demonstrated in Fig. 2, where simulating a 3-qubit depolarizing channel with a linear dilation technique naively requires nearly 10 5 CX gates.
Other difficulties exist in practically implementing these simulations on quantum computers.For instance, random distributions that cannot be sampled in constant or polynomial time could themselves result in memory issues.Another problem exists for U i that cannot be efficiently implemented as products of two-qubit unitaries, i.e. unitaries with exponential complexity in their generators, although this does not exclude their generation via more sophisticated quantum simulators.These problems would also plague ancilla-based approaches.
Our method is related to other stochastic techniques within quantum simulation and open quantum systems simulation, where random effects can help reduce circuit requirements or classical run times [17][18][19][20][21][22][23].Stochastic compilation methods have also notably been used to con-struct first-order dynamics of Hamiltonian evolution [28] and twirled channels [29], and as a tool in implementing certain shallow channels [22].
While random-unitary channels do not encompass all open system dynamics, they are an important subset of quantum channels, and we envision applications of this protocol beyond the examples shown in this work.For example, noise on quantum computers can be modeled as Pauli channels, with empirically determined strength parameters p i [18].Further, more general quantum channels can be simulated using approximate random unitary forms [30].
Modeling environmental effects is essential in understanding the dynamics of many physical systems.The current work proposes a straightforward way to simulate random-unitary noise channels via classical probabilistic sampling of low-depth unitaries.For quantum simulations of random-unitary processes, the proposed hybrid approach has clear advantages over ancilla-based techniques.While the scope is limited to random-unitary channels, these represent an important class of quantum channels and open the door for a variety of other dynamics to be simulated on near-term quantum devices.The authors acknowledge the use of IBM Quantum services for this work.The views expressed are those of the authors and do not reflect the official policy or position of IBM or the IBM Quantum team.

VI. APPENDIX A. Shot allocation
While the hybrid simulation method laid out in Sec.I is robust to different sampling schemes, we optimize the use of valuable computing resources with shot-frugal allocation methods.Practically, some of the operators within a random-unitary channel demand higher priority due to their different weightings p i in Eq. (2b).Throughout these results, we use weighted random sampling of the probability distributions {p i } to allocate shots, thus optimizing the use of computing time [24].More precisely, we sample the operator space {U i } with corresponding probabilities {p i } to determine the circuit run for each shot.This means that for a simulation using a total of N shots, on average p i N shots will be allocated to each operator U i within the sampled space.FIG. 5. Circuit schematic for ancilla-based simulation method [9].The state of the combined ancilla-target system is displayed above each red, dotted line.
To compute the variance of our estimator in Eq. ( 4), we follow the weighted random sampling calculations of Arrasmith et al. [24].This gives the following variance: where ρi := U i ρU † i .In comparison, when using Stinespring dilation, the variance of the estimator ÊS is simply 1 N times the variance of the expectation value ⟨O⟩ ρ: Analytically, these expressions are equivalent: Applying Eq. (2b) to O 2 yields O 2 ρ = i p i O 2 ρi .However, our hybrid sampling can only provide access to the individual expectation values O 2 ρi , making Eq.( 10) the practical means of calculating the variance.

B. Analytic expectation values
The n-qubit depolarizing channel in Eq. ( 5) can equivalently be written as follows [25]: where λ = ( 4 n 4 n −1 )p.This form greatly simplifies analytic calculations, and, in fact, the expectation value of a general observable Ô evolved under E can be calculated as follows: From this form, we computed the analytic expectation values ⟨ Ẑ(i) ⟩ and X(2i) X(2i+1) in Eqs. ( 6) and (7), respectively.

C. Ancilla-based algorithm
For the depolarizing channel simulations, we compare with results from an ancilla-based simulation algorithm [9], described as follows.Consider the randomunitary channel E(ρ) = m i p i U i ρU † i .First, encode the probabilities p i as amplitudes of ⌈log(m)⌉-many ancilla qubits.Suppose a gate P (θ) exists such that: Then, apply each unitary evolution U i to the target state |ψ⟩ conditioned on the ancilla state |i⟩.This accomplishes the desired effect of applying each unitary U i with probability p i .The overall algorithm is displayed in circuit form in Fig. 5.

FIG. 2 . 1 √ 2 |00⟩ + e iπ/ 4 √ 2
FIG.2.Simulation of the n-qubit depolarizing channel on the ibmq montreal device.The black and blue curves use the starting state |0 n ⟩ and compare the ancilla-based approach to the proposed hybrid sampling method.As suggested by the red × marks, the ancilla data is shown purely to contrast with the hybrid method.Based on the average errors, the number of CNOTs, even for n = 1, would result in a highly decohered state, so the specific MSE values are not particularly meaningful.The green curve performs the same hybrid simulation but with each pair of qubits prepared in the non-Clifford (NC) starting state 1 √ 2 |00⟩ + e iπ/4 √ 2 |11⟩, and with measurements in the X basis.Each point corresponds to 10 3 shots.

2 ρ
FIG. 3. Hamming weight distribution for the 27-qubit depolarizing channel, with starting state |0 n ⟩.From the 1000 total measurements, we plot the frequency of each Hamming weight.The blue bars are data from the ibmq montreal device, the red dots are analytically calculated ideal results, and the black diamonds are analytically calculated results that account for bit-flip errors.The bit-flip effects are most visible at the start of the distribution, as many bit-flip errors symmetrically cancel around the center.

FIG. 4 .
FIG. 4. Simulation of the two-qubit TFIM subject to depolarizing noise.The points represent data from ibmq kolkata, and the smooth curves represent the analytic results.The top figure displays populations over time in the standard basis, while the bottom figure displays populations in the eigenbasis of the TFIM Hamiltonian.Note that analytically p01(t) = p10(t), so we omit p01(t) for visual clarity.
This work is supported by the NSF RAISE-QAC-QSA, Grant No. DMR-2037783; the Department of Energy, Office of Basic Energy Sciences, Grant No. DE-SC0019215; the NSF DGE NRT-QISE, Grant No. 2125924; the NSF ECCS CAREER, Grant No. 1944085; and the NSF CNS, Grant No. 2247007.