Beating the House: Fast Simulation of Dissipative Quantum Systems with Ensemble Rank Truncation

We introduce a new technique for the simulation of dissipative quantum systems. This method is composed of an approximate decomposition of the Lindblad equation into a Kraus map, from which one can define an ensemble of wavefunctions. Using principal component analysis, this ensemble can be truncated to a manageable size without sacrificing numerical accuracy. We term this method \emph{Ensemble Rank Truncation} (ERT), and find that in the regime of weak coupling, this method is able to outperform existing wavefunction Monte-Carlo methods by an order of magnitude in both accuracy and speed. We also explore the possibility of combining ERT with approximate techniques for simulating large systems (such as Matrix Product States (MPS)), and show that in many cases this approach will be more efficient than directly expressing the density matrix in its MPS form. We expect the ERT technique to be of practical interest when simulating dissipative systems for quantum information, metrology and thermodynamics.


I. INTRODUCTION
In the emerging field of quantum technologies, environments play a crucial role. Their most obvious contribution is to introduce noise which destroys the delicate coherences necessary for truly quantum effects, but system-environment interactions can also harbour some surprising benefits when carefully controlled [1,2]. For example, one may use dissipation to prepare entangled states [3,4], induce nonreciprocal photon transmission [5] and amplification [6], exploit cooperative effects to reduce entropy [7], as well as engineer dynamics [8]. Given the manifold challenges and opportunities presented by open quantum systems, it is of vital importance to have a fast and scalable method for simulating them.
The first problem one encounters in modelling open systems is how to account for the fact that a full system +environment amalgam exists in a Hilbert space far too large to be simulated exactly. One must instead focus on an effective description of the system of interest, eliminating the environmental degrees of freedom. A number of approaches have been developed to tackle this problem, with one of the earliest examples being the Redfield equation [9][10][11]. While this equation provides a good approximation to dynamics at weak coupling, it does not guarantee the positivity of the density matrix [12], and as such has remained controversial [13].
While influence functional techniques (being nonperturbative) are particularly useful at strong environmental coupling, in the weak coupling regime it is far more common to employ the Lindblad master equation [39][40][41]. This equation represents the most general description for a Markovian process that preserves the essential features of the density matrix (complete positivity and trace preservation) [42,43]. It is this generality that has led to the Lindblad equation acquiring its status as the 'default' approach to simulating open systems [44], and has found application in a very broad range of physical settings [45][46][47][48][49][50]. Indeed, it has recently been shown that the secular approximation [51] used to derive the Lindblad equation is unnecessary, giving the equation a much broader regime of applicability than was previously supposed [13,52]. Recent work has also shown that the Lindblad equation can be used for systems undergoing various kinds of broadband control [53].
Given all this, an efficient scheme for computing Lindblad dynamics is of great importance across the many disciplines that use it, enabling the simulation of larger systems which may in turn lead to new discoveries. Achieving this however is problematic, as simulating this master equation requires the density matrix, which for a Hilbert space of size N H has N 2 H elements. This means that simulating the Lindblad  Given the exponential increase in N H that occurs as one adds particles to a many-body system, the difference between O(N 3 H ) and O(N 2 H ) translates into a tremendous disparity in simulation runtime between the Lindblad and Schrödinger equations for even a moderately large system. While in the last decade enormous strides have been made in the approximate simulation of large arXiv:2010.05399v1 [quant-ph] 12 Oct 2020 systems (see for example recent work simulating strongly interacting quantum thermal machines [54]), the methods employed do not address the fundamental difference in scaling between density matrices and wavefunctions. Instead, the density matrix is flattened into a vector [55][56][57] and the power of tensor network decompositions [58] are leveraged to efficiently evolve this large system. In this sense, one has simply transposed the density matrix into a wavefunction living in a Hilbert space of size N 2 H rather than N H . For this reason, there is a powerful motivation to express the Lindblad equation using wavefunctions of size N H rather than the full density matrix.
The standard solution to this issue has long been to employ a Wavefunction Monte-Carlo (WMC) method (although similar stochastic methods using classical noise have been proposed [59]). This technique exploits the fact that the evolution of the Lindblad equation is equivalent to averaging over an ensemble of stochastic wavefunction evolutions [60][61][62][63][64][65]. If the number of realisations n traj needed to converge the averaging to an accurate answer is much less than N H , then WMC provides an efficient alternative to the exact evolution of the Lindblad equation [66,67]. Of course, like all stochastic methods the number of trajectories needed to converge a WMC result will depend on the specifics of the model system, as well as the length of time being simulated. Consequently, there will be circumstances in which the n traj required will match or even exceed N H , negating the advantages of employing WMC.
In this paper we introduce a new wavefunction based method for describing Lindblad dynamics, which we term Ensemble Rank Truncation (ERT). This method is based on deriving a general expression for the Kraus map of an infinitesimal evolution of the Lindblad equation (up to the second order error in the time step, O(dt 2 )). This mapping provides a deterministic procedure for constructing an ensemble of wavefunctions which together accurately describe system expectations. The size of the ensemble generated through this method increases exponentially with time, necessitating a second procedure to counteract this. Once the size of the ensemble has exceeded a prespecified rank R, principal component analysis is performed to generate a truncated set of wavefunctions in the basis which most efficiently represents the density matrix. This method differs sharply from WMC in two important respects-first, its dynamics are entirely deterministic and do not require any stochastic averaging. This leads to the second distinguishing property, namely that the size of the ensemble produced in this method tends to be much smaller than that required by WMC, with R n traj . Our results show that these features allow ERT to 'beat the house', providing a great improvement over WMC for systems with weak environmental noise (damping and/or dephasing).
A full derivation of our ERT method is presented in Sec. II, together with a formal demonstration that using ERT with Matrix Product State (MPS) simulations of pure states is more efficient than directly using an MPS representation for a mixed-state evolution. Section III demonstrates the efficacy of ERT in three separate systems -a Heisenberg spin chain, a set of atoms in a cavity, and a dissipative Fermi-Hubbard model. In particular, we show here that ERT can provide order of magnitude improvements in both speed and accuracy simultaneously as compared to the standard WMF approach. We conclude in Sec. IV with a discussion of the results presented here, and suggest potential improvements and extensions to the ERT method.

II. ENSEMBLE RANK TRUNCATION
In this section, we derive the approximation we term Ensemble Rank Truncation (ERT). This approximation consists of two parts -an approximate Kraus map which allows the dissipative evolution of the system to be represented with an ensemble of wavefunctions, and a truncation of that ensemble to its principal components, which is equivalent to excluding components of the density matrix above a certain rank.

A. Approximate Kraus operators
It is well known that the time evolution of any system (closed or open) must be a Completely Positive Trace Preserving (CPTP) map [68], and that such maps on the density matrix can always be decomposed into a Kraus form [69], given by: whereρ(t) is the density matrix and theM k are Kraus operators which must collectively satisfy K k=1M † kM k = I. ( Expressing a dissipative evolution in this way has obvious advantages compared to the usual master equation, for which the Kraus form may be formally regarded as a solution. Furthermore, if one begins with a pure statê ρ(0) = |ψ ψ|, it is possible to express the expectation of an operatorÔ purely in terms of wavefunctions, rather than the full density matrix. To do so, one defines a set of wavefunctions {|ψ k } = {M k |ψ } K k=1 , with which expectations can be expressed as: The advantage of this form is that rather than an N H × N H density matrix, one instead requires only K wavefunctions of dimension N H . In the typical case N H K, it is far more efficient to calculate expectations via Eq. (3) than with the full density matrix. Of course, this approach is only possible if one knows thê M k that describe the system of interest.
Unfortunately, while there are many existence proofs for Kraus operators describing open system evolutions [70][71][72] (including for controllability of those systems [73]), the number and explicit form of these operators are known only for a few special cases [74]. This is a problem, as one often wishes to solve for the dynamics of a system described by a Lindblad master equation [41], given by whereĤ is the Hamiltonian and theÂ k are the environmental dissipators. As discussed previously, for large Hilbert space dimension N H , the calculation of the density matrix as compared to a set of K wavefunctions is extremely expensive. We would therefore like to find a way to express this equation in a Kraus form.
To do so, we first discretise Eq. (4): It is then possible to represent this expression as a Kraus map: where the operatorsÛ k andV k are defined aŝ Note that in the case of Hermitian dissipatorsÂ † k =Â k , the Kraus operators are unitary, and the infinitesimal evolution has the form of a random unitary channel [75].
In order to prove Eq. (6), we use the Zassenhaus [76] approximation, and after substituting ∆ = √ dt, X =Ĵ k √ dt and Y = ∓i √ KÂ k , we obtain the following representations for the operators to O dt 2 : These approximations can be used to expand a single term in the sum on the RHS of Eq. (6) to O(dt 2 ) : (13) and after some tedious algebra, one findŝ We now perform a final expansion of e dtĴ k to O(dt 2 ) and substitute both this andΓ k (t) back into Eq. (13) to obtain: Finally, inserting this equality into Eq. (6) recovers Eq. (5) to O(dt 2 ), proving thatÛ k andV k are the soughtfor approximate Kraus operators for an infinitesimal evolution. Note also that these operators approximately satisfy Eq. (2) with:

B. Truncating the Ensemble
For a single time-step, the computational complexity of performing the matrix-matrix calculation for the exact Lindblad equation is O(N 3 H ), but employing the approximate Kraus form, obtaining the set of wavefunctions {|ψ k } = {Û k |ψ ,V k |ψ } K k=1 needed to calculate expectations has a computational complexity of only O(KN 2 H ). Furthermore, using the approximate method one need only store KN H elements, whereas the density matrix requires N 2 H elements. Clearly for a single step, if N H K, both the speed and storage requirements for the calculation will be improved by approximately a factor of N H . Of course, the major problem with this method is that after each step, the number of wavefunctions increases by a factor of 2K, so that after M steps one must store a set of (2K) M wavefunctions.
To overcome this bottleneck, we employ a further approximation, aimed at limiting the size of the wavefunction ensemble. First, let R denote the pre-specified maximal rank of the density matrixρ. The choice of R is dictated by physics, the desired accuracy, and the available memory. For simplicity we assume that the initial density matrix is represented bŷ although the generalisation to a mixed density matrix is trivial. The procedure is then as follows -we propagate according to Eq. (6), saving the set {|ψ k } at each step until there are L = (2K) M > R wavefunctions. Now the size of this set exceeds R, we orthogonalise them such thatρ with {|ψ k } being orthogonal but unnormalised. This orthogonalisation is achieved by where U is a unitary matrix obtained from diagonalising the overlap matrix S ij = ψ i |ψ j , such that Here we assume that the eigenvalues are arranged in descending order, w 1 ≥ · · · ≥ w L . Note that obtaining the overlap matrix involves L 2 wavefunction dotproducts, meaning its computational cost is O(L 2 N H ).
For N H L 2 , this has a negligible cost relative to other operations. Furthermore, since the dimensionality of matrix S is L×L, its diagonalization is also computationally "cheap" and does not depend on the dimensionality N H of the Hilbert space.
Having calculated U , we can truncate it to an R × L matrix U R , using the eigenvectors associated with the R largest eigenvalues. U R can then be used to generate an orthogonal ensemble of wavefunctions: The generation of this truncated ensemble is akin to principal component analysis, a well-known statistical technique where R is the number of principal axes [77]. Using this procedure we need only store {|ψ R k } R k=1 , and after appropriate normalisation of this set we obtain the sought for ERT approximation, Propagating to later steps then simply repeats the same process of orthogonalisation followed by truncation to R wavefunctions, as shown in Fig. 1. Finally, we emphasize that stepping forward in time by generating the set {Û m |ψ R k ,V m |ψ R k } K m=1 is a fully parellelisable calculation, so that if 2KR does not exceed the number of threads available, the time required for stepping forward will be approximately independent of the rank. Furthermore, in the case of large systems where the calculation Schematic showing the ERT process for a single dissipator and R = 2. The initial wavefunction is propagated by the application of the U and V unitaries until the size of the set exceeds R. Each time this occurs the orthogonalisation procedure is performed and the R largest components of the set are retained for propagation to later steps.
of matrix exponents is prohibitively costly, the effect of applying any given Kraus operator to a member of the ensemble can instead be characterised as an ODE (with an explicit time-step of √ dt). A repeated propagation via Eq. (6) can also be represented as a non-commutative Newton binomial [78,79]. This insight offers a potential route to speed up calculations.

C. Combining Ensemble Rank-Truncation with Matrix Product States
It is worth taking a moment to place the ERT method in its proper context. It is the product of two separate approximations -an approximation for the infinitesimal Kraus operators which allows one to characterise the dynamics with a set of wavefunctions, and a truncation of that set to its principal components at each timestep. The former approximation is to the best of our knowledge a novel representation of the density matrix evolution, while the reduction of the set of wavefunctions to their principal components is designed to constrain the size of the ensemble to a manageable size.
It is interesting to note that both ERT and the matrixproduct-state (MPS) representation [80][81][82] make use of a singular-value decomposition to find the most efficient representation for a state. In the case of MPS it is used to find minimal-rank Schmidt decompositions to minimize the matrices that form the MPS representation for a pure state. Here it is used to minimize the number of pure states that are required to represent a mixed state.
In fact, there is no reason why ERT cannot be performed on sets of wavefunctions evolved using time-dependent MPS (time evolving block decimation (TEBD) [81]). One might ask whether combining ERT with an MPS representation for pure states provides any advantage versus directly using MPS to represent the evolving vectorized density matrix. To answer this, recall that in the MPS representation [80,81], a quantum state |Ψ is written as where {|i 1 ⊗ · · · ⊗ |i n } is a computational basis for an n-body system, and χ is the Schmidt rank quantifying the degree of entanglement. In other words, a state |Ψ is represented by n tensors Γ [1] , . . . , Γ [n] (each with χd elements) and n − 1 length χ vectors λ [1] , . . . , λ [n−1] . We may apply the MPS represention directly to a density matrixρ that has been flattened into a column vector |ρ [56]. Using Eq. (23) the vectorised density matrix is expressed as: where is a set of the vectorized matrices {|k n l n |} d−1 k,l=0 . In this case the tensors Γ [k] will each be composed of χd 2 elements. Consequently the vectorised density matrix in Eq. (24) will possess nd times more parameters than the wavefunction in Eq. (23). It therefore follows that so long as 2KR < nd, it will be more efficient to capture the system behaviour using ERT to construct an ensemble of R MPS wavefunctions rather than evolving a single vectorised density matrix with MPS. The reason for this boost is that with the help of principle component analysis ERT finds the optimal basis to represent a given density matrix; whereas, MPS explores the sparsity of the vectorized density matrix.

III. SIMULATION RESULTS
We now demonstrate the utility of ERT over WMC by applying it to two typical multi-body systems of enduring interest (a Heisenberg spin-chain and a collection of twolevel systems coupled to a cavity mode). We choose these systems to be small enough that the Linblad equation can be simulated directly for comparison, but large enough that the ensemble methods are significantly faster.
We also provide an example of applying ERT to fermionic systems by calculating the power spectrum generated by a driven dissipative Fermi-Hubbard system. The convergence of ERT to the exact result is first checked at a smaller system size, before being applied to a system too large to practically calculate (on the hardware used) the exact dynamics.
In order to assess the performance of the we introduce an integrated error E for a set of observables {Ô j }: where O j (t) is the observable expectation from the exact simulation and O A j (t) is the expectation calculated using an approximate method. In the example simulations considered, when the integrated error is on the order of 10 −1 the accuracy is more than adequate for most purposes, and an error of 10 −2 represents an excellent reproduction of the exact dynamics.

A. Heisenberg Spin-Chain
We first consider a simple N site spin-chain system, specifically a Heisenberg XXX model described by the Hamiltonian: Historically Hamiltonians of this type have been used to model magnetic systems [83] in order to calculate their critical points and phase transitions [84]. More recently, this class of models has also been used to study exotic phenomena such as many-body localisation [85]. Given our principal interest in this model is to assess the performance of ERT, in all cases we set h = J = 1, and initialise the system with all spins in the +x direction (i.e. σ (j) x = 1). We choose N = 12 sites corresponding to a Hilbert space dimension of N H = 2 12 = 4096.
In addition to the dynamics of the chain itself, we shall include two types of dissipator -an onsite dephasinĝ and a set of spin injection/absorption operators at the chain terminals: For both types of dissipator, γ and Γ characterise the strength of the environmental coupling to the spin system. For the latter, µ biases the driving (analogously to a chemical potential). Lindblad master equations are only appropriate when the damping is sufficiently slow compared to the internal dynamics. We explore the performance of ERT for damping rates Γ ∈ [10 −3 h, 10 −1 h], which covers essentially the full range over which Linblad equations are of interest. We explore the behavior both without dephasing and with a weak dephasing rate of γ = 10 −3 h.
The numerical simulations of this spin-chain are performed using the Quantum Toolbox in Python (QuTiP) [86,87], which includes optimised implementations of both the exact Lindblad master equation, and a parallelized WMC approximation. All simulations are performed using 4 cores of an Intel Xeon 8280 CPU (with 8 hardware threads), and the ERT method is implemented without any explicit parallelisation. Despite this, Fig. 2 & Fig. 3, show that even an unoptimized Python implementation of ERT is able to perfectly reproduce the exact Lindbladian dynamics, improving upon the exact simulation runtime by almost two orders of magnitude. More significantly, the same ERT implementation also significantly outperforms the optimised C implementation of the QuTiP WMC code in both accuracy and runtime.
A more systematic comparison of results is shown in Fig. 3, demonstrating the relative performance of the lowrank and WMC methods both with and without dephasing, while varying Γ between Γ = 10 −3 h and Γ = 10 −1 h. The first conclusion that can be drawn is that ERT yields larger advantages at weaker couplings. This fits with the natural intuition that at weaker couplings the set of wavefunctions generated by the unitaries in Eq. (6) will be well characterised by a relatively small set of orthogonal components. Nevertheless, even at the strongest damping considered, ERT offers comparable accuracy to WMC. In the case of non-zero dephasing, the presence of an additional 12 dissipators reduces the relative accuracy of ERT at lower ranks (although it remains more accurate than the equivalent WMC simulation for the majority of points). Significantly, runtimes are increased by the presence of additional dissipators, but this is to be expected when the calculation of the new wavefunctions at each timestep is run serially rather than in parallel. In this instance, the runtime for a single timestep will be roughly proportional to 2KR, as one must apply each of the 2K unitaries to the set of R wavefunctions. Since each of these operations is independent however, significant efficiencies could be achieved by parallelising this process when more cores are available.

B. Ensemble of two-level systems in a cavity
Our second example is a hybrid system consisting of a single cavity mode coupled to a number of otherwise independent two-level systems (which represent, e.g., cold atoms [88], color centers [89], or rare-earth dopants [90]). We assume that all the emitters are coupled to the same zero-temperature bath, and the cavity is coupled to a separate input/output port (e.g. a lossy end-mirror). In a recent development, a master equation in the Lindblad form has been obtained that accurately describes open systems regardless of the proximity (or otherwise) of the various transitions (and that is valid for all temperatures) [13]. Employing this master equation, placing the two-level systems on resonance with the cavity mode and working in the interaction picture, the effective system Hamiltonian is given bŷ Here Λ ij = ∆ i δ ij + λ i λ j includes the detuning between emitter i and the cavity, ∆ i , as well as the respective bath-induced Lamb shifts, λ i . Hereσ

(j)
− is the lowering operator for the jth atom from its excited to its ground state andâ is the annihilation for the cavity mode. The term proportional to g is the coupling between the emitters and the cavity mode, and the last term describes driving of the cavity with a coherent state at the rate of |β| 2 photons per second [91]. The Lindblad dissipators for the system are given bŷ which describes output coupling of the cavity mode at rate κ and collective decay of the exited states of all the emitters at the rate γ.
Once again we focus on the regime of weak coupling, setting g = 1, Λ ij = 20jgδ ij and κ = β For all simulations we take 8 atoms and initialise each atom in an equal superposition of its ground and excited states, along with an empty cavity. Simulations are run with the same resources as in the previous section, and as Fig. 4 demonstrates, ERT is able to outperform WMC in capturing the decay of the atomic ensemble at a smaller computational cost. When varying γ and κ, Fig. 5 shows once again that for weak couplings ERT provides a clear advantage over WMC in both efficiency and accuracy. Even at the largest γ considered, we see that at higher ranks the performance of ERT becomes comparable to WMC. This trend is repeated both when including off-diagonal elements in Λ ij , and initialising the system from different states (not shown).

C. Dissipative Fermi-Hubbard model
Finally we apply the ERT method to a fermionic system, specifically the Fermi-Hubbard model [92]. This model is a useful benchmark, given its status as one of the paradigmatic models for strongly correlated electronic systems. Its Hamiltonian is given bŷ whereĉ jσ is the fermionic annhilation operator for the jth site and spin σ, t 0 is the hopping parameter and U is the onsite potential. The Fermi-Hubbard model is of particular interest as under coherent external driving it exhibits a high degree of complexity in its dynamics [93,94]. When dissipation is added this model contains a number of intriguing features, including symmetry breaking phase transitions [95] and dynamic reversals of magnetic correlations [96].
To incorporate dissipation, we apply fermion injection/absorption operators for each species to the terminals of an N site chain and use these dissipators to drive the system from its (isolated) half-filled ground-state.
In order to simulate this model, we use the Python package QuSpin [97,98], which is capable of efficiently simulating fermionic systems. Note that as QuSpin lacks an efficient WMC solver for dissipative systems, in this example we forgo comparison to approximate methods.
In the case of a coherently driven, non-dissipative system, the electronic current may exhibit an optical phenomenon known as High Harmonic Generation (HHG). HHG occurs when the dipole acceleration a(t) = dJ(t) dt spectrum has a a highly nonlinear response to driving, generating frequencies many multiples larger than the driving field [99][100][101]. It is therefore natural to ask if similar non-linearities are present when one uses incoherent dissipative driving. As an initial test of the applicability of ERT to this model, we calculate the dipole acceleration a(t) in a U = t 0 , N = 6 (N H = 4 6 = 4096) system and compare it to the results of an exact simulation. As Fig. 6 shows, once again ERT is able to obtain good accuracy while still significantly improving on the time taken to run the calculation via exact methods.
Using ERT we are able to investigate whether or not non-linear effects are present in this dissipatively driven model at larger system sizes than can be practically calculated using exact methods. As an example, we consider the U = 0, N = 14 Fermi-Hubbard model with dissipative driving. In this free case, Hamiltonian symmetries reduce the effective size of the Hilbert space to N H = 16384, but this is still well beyond the dimension at which an exact master equation calculation is practical without significant computational resources. One can test the accuracy of the approximation in this scenario in a number of ways -the most straightforward being to check that as the rank R is increased, simulation results converge. A second check in this instance is utilising the fact that the Hubbard model is known to exhibit diffusive transport when dissipatively driven [102]. From this, one would expect that the steady state currents J f should depend linearly on the coupling strength. Figure 7 shows that at sufficiently large ranks this is the case, providing a further check for the accuracy of ERT. Figure 8 demonstrates that incoherent dissipative driving generates non-linear effects in the transient dipole acceleration, analogously to the coherent external driving which produces HHG. In this instance, without an external driving frequency to compare to, we instead define a 'fundamental frequency' ω 0 as the most strongly represented frequency by taking S(ω) = |F F T {a(t)}| 2 and setting max ω S(ω) = S(ω 0 ).
Examining the dipole acceleration spectrum, we find that the transient current contains number of distinct frequencies. This is particularly interesting when contrasted to the case of an isolated U = 0 system driven by an external transform-limited pulse, where overtones only appear at odd integer multiples of the driving frequency [103]. In the present case, while we observe a prominent peak in the spectrum at the third harmonic, One finds that at higher ranks the expected linear relationship is recovered, while lower ranks do not capture the final state sufficiently well to maintain this linearity. Note however that at weaker couplings results agree reasonably well between both R = 16 and R = 128, and it is only at stronger couplings that one is required to use a larger rank to obtain accurate results. there are also a number of off-integer and sub-linear frequencies present. The position of these peaks (including ω 0 ) remain constant over a wide range of Γ and µ, indicating that while the transient system response to constant incoherent driving is non-linear, the dissipative parameters meaningfully affect only the size of that response rather than its spectral character.

IV. DISCUSSION
We have introduced an approximate method for the fast simulation of dissipative quantum systems. This technique is based upon first representing an infinitesimal Lindblad evolution in a Kraus form, from which one can construct an ensemble of wavefunctions which capture the full dynamics of the system. When the size of the ensemble exceeds a prespecified rank R, it is truncated in the space of its principal components. We have termed this method Ensemble Rank Truncation (ERT), and have found a number of significant results. First, this method may employed in combination with Matrix Product State (MPS) wavefunction decompositions, and under certain easily satisfiable conditions ERT + MPS will in principle be more efficient than either MPS alone or MPS + WMC in simulating large dissipative systems.
The performance of ERT was investigated in a number of physically distinct systems, allowing for direct comparison to WMC simulations. Here it was found that in the regime of weak dissipative coupling, ERT offered order of magnitude improvements in both accuracy and computing time compared to WMC. This result is at least partially due to the difference between ERT and WMC in their rate of convergence to the true solution. In the latter case, WMC asymptotically converges at O( √ n traj ), while the convergence of ERT with R will depend upon the specifics of the system dynamics. While one cannot determine a priori the appropriate rank to use with ERT, it should be possible to use system identification techniques [104][105][106] to ascertain the effective rank of the system before beginning an ERT simulation. As we have seen in all the cases considered, the rank R required to accurately model a system with ERT is much lower than either the Hilbert space dimension N H or the number of trajectories n traj required by WMC. An important question is how ERT will scale in parallel implementations using many processors. For WMC each trajectory is independent, and the speedup from multiprocessing will be approximately linear [86]. For ERT, each infinitesimal Kraus map may also be independently applied to each element of the ensemble. This part of the scheme is therefore trivially parallelisable and we would expect a similar linear scaling. Although the orthogonalisation step requires some cross-communication to calculate the overlap matrix S, we found in Sec. II that the cost of this step compared to evolving an individual element of the ensemble will be small when N H L 2 . We therefore expect that in this scenario, a properly parallelised ERT implementation will scale with number of processors approximately as well as WMC.
While the advantages of ERT have already been demonstrated, there remains a good deal of scope for improving the method. For instance, one might choose to adopt an adaptive rank R at each step. As a simple example, one could improve accuracy by choosing R at each step such that the sum of discarded eigenvalues in Eq. (20) is below some tolerance. Furthermore, it should be possible to improve on the O(dt 2 ) error found in Eq. (6). This could be achieved either by re-derivinĝ U k ,V k using higher order Suzuki-Trotter decompositions [107], or through the application of Richardson extrapolation to Eq. (6). A final avenue for extension would be to improve the domain of applicability of the ERT, deriving the infinitesimal Kraus operators for the non-Markovian generalisation of the Lindblad equation [108].
In summary, ERT is a novel technique for the simulation of dissipative quantum systems that is capable of outperforming existing techniques, and can in principle be alloyed to state-of-the-art methods for approximating large quantum systems.