Capturing Non-Markovian Dynamics on Near-Term Quantum Computers

With the rapid progress in quantum hardware, there has been an increased interest in new quantum algorithms to describe complex many-body systems searching for the still-elusive goal of 'useful quantum advantage'. Surprisingly, quantum algorithms for the treatment of open quantum systems (OQSs) have remained under-explored, in part due to the inherent challenges of mapping non-unitary evolution into the framework of unitary gates. Evolving an open system unitarily necessitates dilation into a new effective system to incorporate critical environmental degrees of freedom. In this context, we present and validate a new quantum algorithm to treat non-Markovian dynamics in OQSs built on the Ensemble of Lindblad's Trajectories approach, invoking the Sz.-Nagy dilation theorem. Here we demonstrate our algorithm on the Jaynes-Cummings model in the strong coupling and detuned regimes, relevant in quantum optics and driven quantum systems studies. This algorithm, a key step towards generalized modeling of non-Markovian dynamics on a noisy-quantum device, captures a broad class of dynamics and opens up a new direction in OQS problems.


I. INTRODUCTION
Open quantum systems (OQSs), quantum systems that are coupled to their environment, are ubiquitous in the physical sciences [1][2][3][4][5][6][7][8][9][10] and many classical techniques exist to describe the dynamics of an OQS beyond the Markov approximation [11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29]. While the Lindblad formalism gives an efficient and accurate depiction of the dynamics in the weak coupling regime [30,31], the approach does not extend to systems that are strongly coupled to their environments [11,21]. Strong coupling can lead to non-Markovian effects such as recurrences of quantum properties, which are both important for a fundamental understanding of system dynamics and show promise for aiding in system control [32][33][34][35][36][37][38][39][40]. With the advent of quantum devices and the corresponding search for useful quantum advantage there has been an increased interest in algorithm development for physics and chemistry problems. Surprisingly, algorithm development for the treatment of open quantum systems has been limited to a few theoretical and experimental studies [3,[41][42][43][44][45][46]. While there have been impressive recent strides in the field, [44,46] the lag in the development of this field is in part due to the challenge of the non-unitary evolution of OQSs being cast into the framework of unitary quantum gates. To evolve an OQS unitarily, dilation methods must be used to incorporate the important environmental degrees of freedom into a new effective system. Early work faced computational scaling challenges with this dilation [42]; however, recent advances [47][48][49] have allowed accurate simulation of Lindbladian dynamics on a noisy-quantum device [44].
In this Article, we present a novel quantum algorithm to treat non-Markovian dynamics in open quantum systems by extending the Ensemble of Lindblad's Trajectories method onto a quantum computer (ELT-QC) by invoking the Sz.- * kheadmarsden@seas.harvard.edu † prineha@seas.harvard.edu Nagy dilation theorem [44,[47][48][49]. While previous work has looked at example systems for Markovian and non-Markovian dynamics on quantum computers, here we aim to provide the foundation for a universal and generalizable theory that allows for the investigation into open quantum system properties. We start with introducing the theory behind the ELT-QC method, followed by benchmarking the algorithm, and demonstrating the impact of this approach on problems in quantum optics.

A. Ensemble of Lindbladian Trajectories
Density matrix methods are a natural choice for modeling open quantum systems and there have been a plethora of methods to capture Markovian and non-Markovian dynamics, from perturbative to numerical techniques [11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][50][51][52][53][54]. However, in the non-Markovian regime many of these methods struggle from the same challenges such as maintaining the positivity, and therefore physical nature, of the system density matrix. The Ensemble of Lindblad's Trajectories Method is a recently developed, formally exact method depicted in Figure 1, where true density matrices and auxiliary density matrices at time t are represented as purple and gray respectively, and the variable τ represents the time lag [23]. The ELT method extends the Lindbladian formalism beyond the Markovian regime through the use of an ensemble of trajectories originating from different points in the system's history [23]. Due to the relationship between Lindblad operators and Kraus maps [55], the density matrix remains positive semidefinite for all time, and due to the ensemble average non-Markovian behavior is captured. This method was also generalized to treat systems of multiple fermions with accurate statistics [24]. In the simplest discrete form this is mathematically equivalent to writing the density arXiv:2005.00029v2 [quant-ph] 20 Sep 2021 1. An ensemble of Lindbladian trajectories whose weighted ensemble produces the true density matrix at time t.
matrix as, where T is the maximum memory, ω(τ i ) are the statistical weights of the i th trajectories, and e L(τi) are the propagators. Each trajectory is a Kraus map which we can represent by the following Lindbladian trajectory, where s represents an effective time within the mapping, and the Lindbladian superoperator L can be written in terms of the system Hamiltonian H and Lindbladian matrices C i , which account for the interaction of the reduced density matrix with its environment through M different channels [30]. From the properties of Kraus maps the trajectories produce positive semidefinite density matrices whose ensemble is also positive semidefinite [55].
In the continuous limit, this formulation can be related to a simplified form of the generalized master equation [23], where the traditional memory kernel is represented by the ensemble weights, the Lindbladian superoperator, and a timederivative. In this representation, just as in traditional use of Lindblad's equation, there are multiple avenues of approach to parameter selection, including but not limited to numerical optimization, extraction from short-time propagations, and use of experimental data.

B. Quantum Algorithm
Recent work by Hu et al. has been done to map the Lindblad equation into a dilated unitary evolution, then perform this unitary evolution on a quantum computer [44]. In this work, the Lindblad equation is first written in operator sum form, where the M i are Kraus maps corresponding to the Lindbladian channels represented by each C i in the original equation. Since Kraus maps are contraction mappings of a Hilbert space, the Sz.-Nagy dilation theorem guarantees that there exists a related unitary operator dilation in a larger Hilbert space [44]. While different orders of d-dilation exist and correspond to contraction mappings being applied d times to a Hilbert space, the 1-dilation is sufficient to consider either populations or coherence of a two-level system. Here, we focus on the 1-dilation which produces the unitary operator, where I is the identity matrix and D Mi = I − M † i M i . To perform the unitary evolution, the density matrix is reformulated in vector form and dilated to account for the added environmental degrees of freedom. The density matrix can be represented as an ensemble of vectors {v j } where the vectors need not be orthogonal. The {v j } vectors are then dilated by padding with zeroes to match the dimension of the unitary operator. In this work, only the population elements of the density matrix are considered. Given an ensemble of vectors {v j } and using this dilation, the population elements of the system are calculated as, where ρ k is the occupation number of the k th state and the summation is over the dilations of all the Kraus operators M i acting on all the vectors v j . Each term can be simulated in parallel, with circuits of gate count on the order of n 2 , where n is the dimension of the system [44]. Since the ELT method is an ensemble average of Lindbladian trajectories, the recent algorithm used to calculate Lindbladian trajectories on quantum devices can be used with the following generalization. In previous work when practically invoking the ELT method [23,24], each trajectory is written as,D where the Lindbladian term L is a constant in terms of time, and D(t − τ ) is the density matrix initialized from τ time steps in the system's history. This equation can be directly cast into the unitary dilation framework; however, it requires storage of each D(t − τ ) along the way and the production of a new circuit describing each D(t − τ ), or the production of an arbitrary state at each time point. While this is possible, it would be computationally taxing and preclude the possibility for a quantum speed up. Alternatively, using a variable change, we can shift the time dependence from the density matrix to the Lindbladian, where the Lindbladian term L depends on the time lag while D is a constant density matrix. Instead of requiring an arbitrary state preparation at every step, this formulation requires only variation in the time input to produce the Lindblad matrices and therefore the unitary dilations. The ELT increases the computational cost of a single Lindbladian trajectory in Ref. 44 by a multiplicative prefactor of T where T is the number of trajectories in the ensemble average.

A. Single Lindbladian Trajectory
The unitary evolution matrices of a two-level system in a single amplitude damping Lindbladian channel are given by [44], and where γ is the rate of decay. An initial density matrix of the form, is decomposed into an ensemble of dilated vectors [44], and With this decomposition, the four necessary terms for each Lindbladian trajectory are calculated through the circuits shown in Figure 2. The circuits were constructed and verified in both Qiskit [56] and QuTip [57,58]. It should be noted that the two qubit states follow the conventional notation where the vector (1, 0, 0, 0) represents both qubits in their ground states, (0, 1, 0, 0) represents q 0 in its excited state and q 1 in its ground state, (0, 0, 1, 0) represents q 0 in its ground state and q 1 in its excited state, and (0, 0, 0, 1) represents both qubits in their excited states. Using these circuits, we consider a two-level system in an amplitude damping channel with a system decay rate of γ = 1.52 · 10 9 s −1 . Using the solution from the Lindblad equation on a classical device, IBM's Qiskit simulator [56],  [56] and the x's are generated using IBM's London device [59].

B. Ensemble of Lindbladian Trajectories
Having verified this algorithm for a single trajectory of the Lindblad equation, we extend these circuits to treat the damped Jaynes-Cummings model, which consists of a single excitation in a two-level system coupled to a reservoir of harmonic oscillators [11,50,60,61]. This system is both exactly solvable and known to demonstrate non-Markovian behavior in the strong coupling and detuned regimes, making it an excellent benchmarking system for open quantum system methods [11,23,60]. The Hamiltonian is given by, where ω 0 is the system's transition frequency, λ is inversely proportional to the reservoir correlation time,â † andâ are the creation and annihilation operators, andσ x,y,z are the Pauli spin operators withσ ± =σ x ±σ y . In lieu of the cavity mode treatment, the bath spectral density is used [11], where ω 0 is the system transition frequency, ω the bath frequency, ∆ the detuning, and λ is related to the system-bath coupling strength.
First we consider the strong coupling case where the bath relaxation parameter is given by λ = 0.2γ. In this regime, a single Lindbladian trajectory fails to capture accurate dynamics and a more involved method is required [11,50]. The populations of the Jaynes-Cummings model in the strong coupling regime are shown in Figure 4 using the ELT-QC method with weights numerical optimized as compared to the exact solution using Maple [62], as was done in Ref. 23. The solid lines represent the exact solution while the dots and x's are results from quantum simulation using the circuits shown in Fig. 2 using IBM's simulator and London device respectively [56,59].  [56], and the x's are generated using IBM's London device [59].
The ELT on a quantum simulator and device agrees well with the exact solution, demonstrating ability of the ELT-QC algorithm to accurately capture dynamics in the strong coupling regime.
Next we consider the strong coupling and detuned case, where ∆ = ω − ω 0 = 0. In this regime, the bath relaxation parameter is given by λ = 0.3γ and the detuning by ∆ = 2.4γ. The time-evolution of the populations are shown in Figure 5 where lines represent the exact result and the dots and x's represent results from quantum simulation using the circuits shown in Fig. 2 [56], and the x's are generated using IBM's London device [59].
with the exact solution, demonstrating that the ELT-QC algorithm captures accurate dynamics in the detuned regime for the Jaynes-Cummings model.

IV. CONCLUSIONS AND OUTLOOK
Few algorithms exist for the treatment of open quantum systems on quantum devices, and existing algorithms are either restricted to Markovian systems or inconsistent with complete positivity. Our work presents an algorithm inspired by the Ensemble of Lindblad's Trajectories method [23] using an efficient dilation theorem [44] which allows for the treatment of non-Markovian dynamics of open quantum systems on a quantum device. The ELT-QC method is benchmarked on the Jaynes-Cummings model in the strong coupling and detuned cases on both a quantum simulator and real 5-qubit quantum computer, showing excellent agreement with the exact dynamics. While the Jaynes-Cummings model is a benchmark system, these results show our algorithm's ability to capture a wide variety of open quantum system dynamics accurately with tractable scaling. These regimes include strictly Markovian, or dissipative, dynamics and non-Markovian dynamics induced by either strong or detuned system-environment coupling.
Moreover, this quantum algorithm retains the significant advantages of its classical counterpart including an exact treatment of non-Markovian dynamics and complete positivity of the density matrices, meaning that the density matrices remain positive semidefinite with non-negative probabilities for all time. The ELT-QC algorithm offers a path towards polynomial time-evolution of an electronic system in the presence of a complex environment. An accurate yet tractable description of open quantum systems on quantum devices has a myriad of significant applications from catalytic chemistry and correlated materials physics to descriptions of hybrid quantum systems and spin systems.