Robust measurements of n-point correlation functions of driven-dissipative quantum systems on a digital quantum computer

,

Introduction.-Openquantum systems, and in particular driven-dissipative systems, are among the most difficult problems to study in many-body physics, but also among the richest.The problem parameter space is vast; the bath as well as the system have their own inherent dynamics, and their interaction can be complex.Yet, in some sense there is a unification and emergent simplicity as the details often do not play a role when it comes to describing non-equilibrium steady (or periodic) states.These can be captured with a few parameters, have lost all knowledge of their history, and are stable to perturbations away from their fixed point.In other words, they are remarkably robust.
Their robustness has also recently been exploited in simulations on quantum computers, either relying the hardware intrinsic decoherence [1][2][3], by implementing Kraus maps and Lindblad operators [4][5][6][7][8][9][10][11][12][13][14], or by implementing non-Hermitian dynamics [15,16].In some cases, the existence of a fixed point has enabled quantum computers to perform the simulations far beyond the short coherence time of the qubits when the fidelity of one Trotter step is sufficiently high.It thus appears that driven-dissipative systems are promising problems for applications of near-term quantum computers.
Given this situation, it is critical to develop the tools and methodology to be able to interrogate these long-lived non-equilibrium states.For single-time operators this is not a problem, one simply measures the desired operator at some point during the evolution.However, there is a wide class of observables -correlation functions -that are of equal or greater importance as they describe the excitations of the system, and make a direct connection to experimental observables.The problem is that the typical protocols for the measurement of correlation functions [17][18][19][20][21][22] are based on the Hadamard test [23] [see Fig.
FIG. 1. Overview of the proposed new method and its region of applicability.Ancilla decoherence (blue) limits the region where the Hadamard test method (a,c) yields a result, whereas the reset-based method (b,d) has no such limitation.For closed system evolution (a,b) the system decoherence naturally limits the maximal time of interest, but for driven-dissipative systems (c,d) the reset method is necessary to go beyond the ancilla coherence time.

1(a)],
where the correlation function is measured with an ancilla qubit [19,[24][25][26][27][28][29].This approach does not robustly generalize to measuring correlation functions in open quantum systems, and we illustrate the issues with it in Fig. 1 for a 2-point correlation function.In short, the ancilla cannot capture the potentially long-time dynamics of the driven-dissipative open quantum systems because the ancilla has a short coherence time.For simulating closed quantum systems [Fig.1(a)] this is not a problem because the system has an equally short coherence time, and the region of inaccessibility by the Hadamard test approach has no information.On the other hand, for simulating open quantum systems [Fig.1(c)], this is a problem because the now-stable dynamics of the driven-dissipative system are inaccessible due to the ancilla decoherence.
In this Letter, we propose a strategy that (i) is a full framework for computing n-point correlation functions in open quantum systems, and (ii) is suitable for the near term where we cannot rely on long ancilla coherence lifetimes.Its crux is a measurement of the ancilla right after it is entangled with the system, and using the result of the measurement in post-processing to construct the desired n-point correlation function.We illustrate the idea in Fig. 1 (b) and (d).For simulating closed quantum systems this approach yields the same information as the standard Hadamard test; however, for open quantum systems the region of inaccessibility by ancilla decoherence is now fully accessible.
Our method is a simple strategy capable of measuring arbitrary unequal-time correlation functions between multi-qubit Pauli operators, and which works for both dissipative and unitary time evolution.As such, it subsumes and unifies the approaches of Ref. [30] (unequal-time commutators) and Ref. [31] (unequal-time anticommutators).It is hierarchical in the sense that extracting the information of an n-th order correlation function requires previous knowledge of lower-order correlation functions; but, it restores the robust nature of driven-dissipative systems, because it does not require system-ancilla entanglement to be maintained during the time evolution of the system.We verify the validity of our method by performing measurements of the single-particle Green's function of a driven-dissipative fermionic model using a Quantinuum quantum computer.Our results show excellent quantitative agreement between data and the theoretical predictions.
Target quantities.-Ourgoal is the calculation of correlation functions of a generic system (S) that can also dissipate energy through an interaction with a bath (E); so we employ the density matrix formalism, which is required to study open quantum systems [32].
The correlation functions are constructed as follows.Let {O i } be a set of operators in the Schrödinger representation acting on the system Hilbert space with i = 1, 2, • • • n, and let {t i } being a set of ordered time values such that t 0 < t 1 < t 2 < . . .t n−1 < t n , where t 0 is the initial time, then we define the n-th rank correlation function via Here, O i (t i ) is the operator O i in the Heisenberg representation, ρ(t 0 ) is the system density matrix evaluated at the initial time, V ti+1,ti is the time evolution superoperator that evolves the system from time t i to t i+1 (i.e.ρ(t i+1 ) = V ti+1,t i ρ(t i ) acting from left to right), and Tr S indicates a trace over the system subspace (mean-ing that the degrees of freedom of the bath have already been integrated out).For simplicity, and without loss of generality, we assume that the operators O i are unitary and Hermitian operators; addressing this case is sufficient to demonstrate the validity of our method, because a non-unitary operator can always be expanded as a linear combination of unitaries, chosen to also be Hermitian (e.g.Pauli strings).As will be shown in the next section, the correlation function in Eq. ( 1) can be extracted from the Hadamard test.
The alternative strategy that we propose will naturally yield correlation functions of nested commutators and anti-commutators of the form where [. , .] ± can be either commutators (-) or anticommutators (+), all chosen independently.The correlation function in Eq. ( 1) can be obtained from the one in Eq. ( 2) and vice versa by performing multiple measurements and then combining the different outcomes together.[33].We note that in the case of two-point functions, Eq. ( 1) corresponds to lesser or greater Green's functions while Eq. ( 2) to advanced, retarded, and Keldysh Green's functions [34], so both methods produce all the physical Green's functions needed to describe a time evolving quantum system.However, there are some limitations -for example one cannot directly calculate out-of-time-ordered correlation functions with the circuit in Fig. (3) and we leave possible generalizations of this method to future work.
n h 6 z J B L t j D + y R P X n 3 3 r P 3 4 r 1 + j Q 5 5 g 5 1 F 9 g P e 2 y f b y q I P < / l a t e x i t > U tn,tn 1 K < l a t e x i t s h a 1 _ b a s e 6 4 = " Z A a p n 9 r E L P n x o w 9 N 5 3 0 g I Z m 3 p u I / 3 m t i P p n n V j 6 Y U T o i 8 k h k g q n h 4 z Q M u 0 E e U 9 q J I J J c u T S 5 w I 0 E K G W H I R I x S g t K Z f 2 4 c x / v 0 j q 5 Z J z X C p f n R Q r 5 7 N m s m y f H b A j 5 r B T V m G X r M p q T L B 7 9 s S e 2 Y v 1 a L 1 a b 9 b 7 z 2 j G m u 3 s s T + w P r 4 B 8 O e Z A g = = < / l a t e x i t > FIG.
2. The standard interferometric scheme for measuring the n-time correlation function [23], as given in Eq. ( 1), for a dissipative circuit.Accurate results require that the ancilla register A1 maintain coherence over the entire duration of the circuit.
Hadamard test for driven-dissipative systems.-InFig. 2, we show how the interferometry scheme proposed in Ref. 23 generalizes to compute the n-time correlator defined in Eq. ( 1) for an open quantum system.In order to simulate dissipative dynamics, we need a generic k-qubit ancilla register (called A 2 ) that we take to be initialized into the state |0⟩ = |0⟩ ⊗k .A suitable unitary operation U t,t ′ K that entangles A 2 with the system register S followed by tracing out (ignoring) the state of the ancilla register can encode the non-unitary time evolution map V t,t ′ , which can be rewritten using the Kraus sum n h 6 z J B L t j D + y R P X n 3 3 r P 3 4 r 1 + j Q 5 5 g 5 1 F 9 g P e 2 y f b y q I P < / l a t e x i t > U tn,tn 1 K < l a t e x i t s h a 1 _ b a s e 6 4 = " Z A a p n 9 r E L P n / 3. Circuit to measure a generic n-time correlation function of the kind defined in Eq.( 2) using the robust strategy. representation: where K i are the so called Kraus operators satisfying the sum rule i K † i (t, t ′ )K i (t, t ′ ) = I.They are related to the unitary evolution of the system and ancilla bank by , with {|i⟩} being a complete basis for A 2 .In the interferometry scheme, we need an extra single-qubit ancilla register A 1 in which all the information about the correlation function (which is a complex number) will be stored.For example, in the case of n = 2, the final quantum state of the A 1 qubit reads: where C = ⟨O 2 (t 2 )O 1 (t 1 )⟩, and the binary variable α = {0, 1} indicates whether the S gate was applied or not.
Measuring the ancilla in the Z and Y bases determines the real and imaginary parts of the correlation function.This method is convenient because the complex information encoded in the correlation functions of a many-body system are found from single qubit measurements.However, this scheme requires maintaining the coherence of the A 1 ancilla (and thereby its entanglement with the system) for the full duration t n − t 1 .In the next section, we introduce an alternative robust scheme that does not require maintaining coherence of the A 1 ancilla, but at the cost of requiring a more complex measurement scheme.
Robust strategy.-InFig. 3, we show the alternative circuit to measure the correlation function defined in Eq. ( 1).This circuit is schematic, because it encodes all possible circuits that are employed to measure the set of correlation functions in Eq. ( 2).Here, each realization has chosen unitary operations acting on A 1 (selected from (S) αi H, where S and H are the phase gate and the Hadamard gate, respectively) for each time t i measured in the correlation function.The circuit shown in Fig. 3 naturally measures the set of correlation functions defined in Eq. ( 2) with the commutator or anti-commutator chosen from the n − 1 dimensional binary vector α = {α 1 , α 2 , ..., α n−1 }.It is important to note that after the S αi H operation is performed, the ancilla qubit A 1 is measured immediately afterwards and the measurement outcome (m i ) is stored; such a measurement destroys the entanglement between A 1 and the state encoded in the system and the A 2 ancilla bank.The state is then evolved to the next t i using the Kraus map decomposition defined in Eq. ( 3).The A 1 ancilla is then reset to its |+⟩ state and the process is repeated for each operator in the correlation function.In the last step, after the final time evolution from t n−1 to t n , the S register qubits will be in a final state ρ n and the operator O n is measured directly on the S register qubits, yielding results that depend on α.The correlation function is determined by classical post-processing of the accumulated results and the choice of α.
In general, the state of the system qubits at time t j+1 is obtained from the state at t j through the following map [35]: where the proportionality constant is given by tracing the RHS of the equation.Here, m j = {0, 1} is the result of the A 1 qubit measurement, and ρ j=1 = ρ(t 1 ) is given by the initial state of the system at time t 1 (see Fig. 3).
In order to show how this method works in practice, we discuss the two simplest cases: i.e. the two-point and the three-point correlation functions.For n = 2, the result of measuring O 2 directly on the system register will yield where ).Hence, when α 1 = 0, 1 the term in square brackets in Eq. ( 6) is proportional to ⟨[O 1 (t 1 ), O 2 (t 2 )] ∓ ⟩.This is precisely Eq. ( 2) when n = 2.In order to isolate our target quantity, i.e. the two-time correlation function ⟨[O 1 (t 1 ), O 2 (t 2 )] ∓ ⟩, we need to subtract the first two terms in Eq. ( 6) multiplied by N which are determined by equal-time averages that can be obtained by performing a set of simple extra measurements [35].It is worthwhile to note, that the measurement outcome of the ancilla m 1 = 0, 1 must be stored for extracting our target quantity.
For n = 3, measuring O 3 results in the following quantity: where C α t1,t2,t3 is a three-time correlation function that depends on the values of α.There are four possible values In addition, there are contributions denoted by R α , which is a remainder function.It is determined by performing additional measurements comparable to what is needed for lower-rank correlation functions (see supplemental information for details: [35] ).
We note that in the case of single-qubit [30,31,36] and two-qubit [37] correlators, there are alternative ways of measuring correlation functions that do not require the extra ancilla register A 1 , at the cost of performing more measurements on the system [35].
Circuit implementing the Trotterised time evolution UK (t, t + ∆t) of the model defined in Eq. ( 8).
Hardware implementation-In order to verify the validity of the protocol, we applied it to measure the Green's function of spinless free fermions in a lattice driven by a constant electric field that also dissipate energy through a coupling with a thermal bath.The Hamiltonian of this chosen system plus bath can be brought into a blockdiagonal form after performing a Fourier transform to momentum space as described in Ref. [7].Hence, the system's reduced density matrix factorizes as a tensor product in momentum space, i.e., k ρ k , and we can define a (diagonal in k) master equation for each 2×2 k-dependent density matrix ρ k , where the Lindblad operators are , with d k being the destruction operator of a lattice fermion with quasi-momentum k, ϵ k (t) = −2J cos (k + Ω t), with J being the hopping amplitude, Ω the amplitude of the applied DC field, and k the crystalline momentum.Γ sets the strength of the systemenvironment coupling and n F (x) = [1 + exp(β x)] −1 is the Fermi-Dirac distribution with β being the inverse of the bath temperature.In Fig. 4, we show the circuit implementing U K for the Kraus map related to Eq. ( 8).
The Lindblad operator L 1(2) encodes the physical process of a Bloch electron (hole) with momentum k to hop from the lattice to the bath with a probability given by Γn F [−ϵ k (t)] (Γn F [ϵ k (t)]).Such a decay process introduces a time dependence of the momentum distribution function of fermions and a damping of Bloch oscillations that eventually leads to a non-zero average of the DCcurrent [7,10,35,38].In Fig. 5, we show the retarded fermion Green's function + measured on Quantinuum's model H1 quantum computer.The retarded Green's function of the model can be computed exactly and its derivation and analytical form are given in the supplemental materials [35].
There is excellent quantitative agreement between the data produced by the quantum computer and the expected curves in presence of noise.It is worthwhile to note that in the presence of a driving field, the Green's function does not oscillate as a simple sinusoidal function and it presents extra features, such as the additional maxima and minima occurring between times 10 and time 19 [see Fig . 5], that are faithfully reproduced by the quantum computer data.

Conclusion &
Outlook.-We have put forward a robust technique for the measurement of multi-point correlation functions of driven-dissipative quantum systems that can be applied in the realm of quantum simulations of complex models such as the Hubbard model.Unlike the Hadamard test, which requires us to keep the ancilla and system qubits coherently entangled, our new approach does not.This is advantageous for driven-dissipative systems, where the system is not coherent (see Fig. 1), although it comes at the cost of performing extra measurements, as well as requiring additional circuits of lower depth than the one needed to extract the target quantity.Our method naturally computes correlators of the form given in Eq. ( 2), which represent a myriad of response functions and experimental measurements.We applied our method to measuring the Green's function of free fermions driven out of equilibrium and interacting with a bath.The data obtained from the quantum computer are in an excellent agreement with the curves predicted by the theory.While this data constitutes an important proof of principle enabling the measurement of correlation functions on near-term quantum computers, further work needs to be done to use this approach to solve new problems in science.
Interestingly, given its generality, the Hadamard test has applications other than the measurement of correlation functions; for example, it has been proposed for determining important overlaps in the realm of variational quantum dynamics simulations [39,40] and also for the simulation of open quantum systems using quantum imaginary-time evolution [13].We therefore expect our robust alternative strategy to the Hadamard test to be suitable for these other applications as well.

Supplemental material for: Robust measurements of n-point correlation functions of driven-dissipative quantum systems on a digital quantum computer
Lorenzo Del Re,

CORRELATION FUNCTIONS ON THE KELDYSH CONTOUR
We discuss more in detail what kind of correlation functions can be evaluated by the circuit in Fig. (3) in the main text and what are those for which a generalization of our scheme is needed.In general, since the operators O i (t i ) do not commute with each other, one can construct n! correlation function out of all possible permutations of the operators in Eq. ( 1) in the main text.However calculating all possible nested commutators/anti-commutators of Eq.( 2) gives the possibility of isolating only 2 n−1 of these permutations.
Here we argue that the permutations we have access to correspond to components of time-ordered correlation functions on a simple two-branch Keldysh contour [see Instead of giving a rigorous proof of our statement, we will consider the specific example for n = 3.In this case, the permutations we have access to through Eq. ( 2) are: where C ± are the forward and backward branches.In the same figure we show that the particular permutation ⟨O 2 O 3 O 1 ⟩ can be seen as a time-ordered function on the contour C where t 1 = t 1+ and t 3 = t 3+ lie on the forward branch while t 2 = t 2− lie on the backward branch.
Conversely, the permutation ⟨O 3 O 1 O 2 ⟩ cannot be interpreted as a time-ordered correlation function on C. Instead, it could be seen as a time-ordered correlation function on the more complex contour C 2 shown in Fig. (S1b) that contains two additional branches.The same contour has to be considered for the calculation of out-of-time correlators [2], and therefore a generalization of our scheme is needed in order to have access to those quantities.
Detail of the circuit presented in Figure 3 in the main text that is responsible for a single time step evolution from time tj to time tj+1.

TIME STEP UPDATE
In this section, we derive Eq. ( 5) of the main text, which describes how the density matrix is updated from the j-th step to (j + 1)-th one.In Figure S2, we show the section of the circuit presented in Fig. 3 of the main text which accounts for the single step time evolution.We note that we have two ancilla registers A 1 and A 2 .A 2 is a single ancilla register given by k-qubits needed to implement the dissipative time evolution (i.e. the Kraus map) of the system density matrix.This register is initialized in the computational basis as |0⟩ ⊗k .A 1 is a single-qubit ancilla register needed for the correlation function measurement protocol.We marked three intermediate points in the circuits (i), (ii) and (iii) which correspond respectively to the density matrices of system plus ancilla register A 1 ρ (i) , ρ (ii) and ρ (iii) .
After the first control operation on the system we have: After the action of the phase-gate S αj and Hadamardgate H on the ancilla A 1 , the density matrix is expressed as follows: After the ancilla A 1 is measured (with result m j = {0, 1}), the density matrix then becomes where If we then perform the dissipative time evolution using register A 2 from time t j to time t j+1 and trace out the bath degrees of freedom, we obtain the following density matrix for the system register: which coincides with Eq. ( 5) in the main text.

EXTRA MEASUREMENTS FOR THE TWO-TIME CORRELATION FUNCTION
In Figure S3, we show the circuit needed to obtain the information about the quantity ⟨O 2 (t 2 )⟩ O1 entering Eq. ( 6) in the main text.We note that the extra-ancilla register A 1 is not needed in this case, because this quantity represents a one-time average.This quantity plus the average ⟨O 2 (t 2 )⟩ multiplied times the normalisation factor N = {2 + (−1) m1 i α1 ⟨O 1 (t 1 )⟩ + h.c.} −1 must be subtracted from the right hand side of Eq. ( 6) in order to isolate the two-time correlation function Hence, in order to compute the two-time correlator, four measurements must be performed, namely the one described in Eq. ( 6), and those yielding A2:

EXPLICIT CALCULATIONS FOR THE THREE-POINT CORRELATION FUNCTION
Here we derive Eq. ( 7) in the main text and give some insights about the remainder function and the additional measurements that are needed for its computation.In the derivation, we use the following abstract notation for the time evolution: V tj+1,tj [ρ(t j )] = ρ(t j+1 ), instead of its sum representation.
The final state before the measurement in the system register in the case of n = 3 reads: where: Because the time evolution preserves the trace, the normalization factor will be given by: There are a total of 16 terms on the RHS of Eq. (S5): only four of them, those that are contained in ρ (III) 3 and its conjugate transpose, contain the information about the three-point correlation function.The other 12 terms give rise to the remainder function R α when O 3 is finally measured.As we shall see, these 12 spurious terms contain information about two-point and one-point (averages) correlation functions.Therefore, in order to extract the three-particle correlation function we will need to perform preliminary measurements of a few averages and two-point correlation functions that can be obtained using similar circuits as that one shown in the main text for n = 2.
If we substitute ρ 2 from Eq. ( 7) into Eq.(S6) we obtain: Analogously we obtain the explicit expression for Eq.(S7), that reads: It is worthwhile to note that all the terms in Eqs.(S10) and (S11) contribute solely to the remainder function.
For example, let us consider the following average that emerges from the last two terms in the RHS of Eq. (S11), i.e.
+h.c.: Such a quantity can be measured using the circuit shown in Fig. S4, where the unitary operation O 2 acts solely on the system qubits after the time evolution from t 1 to t 2 .
|0i A1: A2: n h 6 z J B L t j D + y R P X n 3 3 r P 3 4 r 1 + j Q 5 5 g 5 1 F 9 g P e 2 y f b y q I P < / l a t e x i t > x E / w K r 3 r y J l 7 9 B g / + i 7 s x B z X W q a j q p q s r S J Q 0 Z N s f V m V i c m p 6 p j p b m 5 t f W F y q L 6 + c m T j A c E 4 u c k h q 5 J n X S I J w 8 k m f y Q l 6 N J + P N e D c + Z q M 5 I 9 v Z J 3 9 g f P 4 A 2 + 6 f e g = = < / l a t e x i t > H |0i < l a t e x i t s h a 1 _ b a s e 6 4 = " < l a t e x i t s h a 1 _ b a s e 6 4 = " g 8 2 A D 1 i s 6 / e 8 O k P c r n D t V B i h r 1 a a l 6 k T W T Y w f s k B 0 z h 5 2 x K r t i N V Z n g o 3 Z E 3 t m L 9 a j 9 W q 9 W e 8 / o 0 t W t l N k f 2 B 9 f A M q E J m m < / l a t e x i t > FIG.S4.One of the preliminary measurements needed in order to compute the remainder function Ra in Eq. ( 7).
The explicit expressions of ρ (III) 3 reads: We notice that the first two terms in the RHS of the last equation contribute to the remainder function.For example, the second term in the RHS of Eq. (S12) would give rise to the following average . that could be measure using the circuit displayed in Fig. S5.
n h 6 z J B L t j D + y R P X n 3 3 r P 3 4 r 1 + j Q 5 5 g 5 1 F 9 g P e 2 y f b y q I P < / l a t e x i t > ⇢(t 1 ) < l a t e x i t s h a 1 _ b a s e 6 4 = " x E / w K r 3 r y J l 7 9 B g / + i 7 s x B z X W q a j q p q s r S J Q 0 A c E 4 u c k h q 5 J n X S I J w 8 k m f y Q l 6 N J + P N e D c + Z q M 5 I 9 v Z J 3 9 g f P 4 A 2 + 6 f e g = = < / l a t e x i t > |+i < l a t e x i t s h a 1 _ b a s e 6 4 = " F f f that is the one appearing in Eq. (7) in the main text.
The remainder function R α can be obtained in a similar way but adding up all the remaining 12 terms appearing in Eq. (S5) multiplying them times O 3 and computing the trace.In order to extract the information about C α t1,t2,t3 , we will need to perform additional measurements similar to those shown in Figs.S4 and S5.

SPECIFICS OF THE MODEL SYSTEM IMPLEMENTED ON THE QUANTUM COMPUTER
We model spinless electrons hopping on a onedimensional lattice with nearest-neighbor hopping.The electrons are placed in an electric field and are connected to a fermionic bath that provides dissipation ton the system.See Refs. 3 and 4 for more details.

Time evolution
The time evolution of such a model from time t to t + ∆t is governed by the infinitesimal Kraus map that can be derived from the master equation and it is given Analogously we obtain the following expression for the greater Green's function: ))e −Γ(t−t ′ ) e −if k (t,t ′ ) .(S20) Let us notice that only two Kraus operators, namely K 1 and K 2 were involved in the time evolution of the operator ρ d k (d † k ρ) in the expression of the lesser (greater) Green's function.From Eqs. (S19) and (S20), we obtain the retarded Green's function as following:

HARDWARE IMPLEMENTATION
The data was taken on Quantinuum's model H1-2 quantum computer, a QCCD trapped ion quantum computer that (at the time of taking this data) operated with 12 qubits.Each qubit is comprised of the |F = 0, m F = 0⟩ and |F = 1, m F = 1⟩ hyperfine states of a Yb 171 ion.Memory errors are only routinely measured and reported at the depth-1 circuit time of 20ms (the depth-1 circuit time is the average time required to execute a dense depth-1 circuit with random gate connectivity), at which time they are in the low 10 4 's per qubit.Given that two qubit gate errors were approximately 2 − 3 × 10 3 , we are confident that memory errors made a marginal contribution to the overall circuit error.

RESOURCE ESTIMATES ON THE NUMBER OF MEASUREMENTS
The number of measurements needed to attain a given precision ϵ is proportional to ϵ 2 or equivalently, the precision of our measurements scales as 1/ √ n, which is a consequence of the Hoeffding's inequality.The total computational cost will then be given by n multiplied times the number of gates needed for state preparation (which can be roughly approximated by those contained in the time evolution).In other methods, where more than one observable must be measured [6] this must also be taken into account (as a multiplicative factor).Let us note that the cost of adding the A 1 ancilla is negligible in comparison to the cost of performing time evolution for long times.Moreover, assuming that O i are not entangling gates acting on ℓ qubits, which is a sensible choice given that many correlation functions can be expressed as a linear combination of correlation functions of strings of Pauli operators, entangling the A 1 ancilla used in our scheme introduces a cost of roughly 2ℓ extra CNOTS.However, if we assume that time-evolution is the bottleneck of the quantum circuit and contains far more CNOTS than ℓ, then the extra complexity introduced by entangling the ancilla with the system is negligible.
t e x i t s h a 1 _ b a s e 6 4 = " J Z I S x a F D D i G F 1 u m k m w P L y 4 k A r w F r 9 f J u 1 6 z T u r 1 W / P K 4 2 r v J k i O S L H p E o 8 c k E a 5 I Y 0 S Y t w o s g T e S Y v z t x 5 d d 6 c 9 5 / R g p P v H J I / c D 6 + A a X 0 l F c = < / l a t e x i t > t e x i t s h a 1 _ b a s e 6 4 = " S A Y v 1 Z n A P e b a P I l A 7 b k

FIG. 5 .
FIG. 5.Imaginary (upper panel) and real (lower panel) parts of the retarded fermion Green's function as a function of time (parameters are wavevector k = −0.5/a0where a0 is the lattice constant, dimensionless electric field is Ω = 1, the dissipation rate to the fermionic bath is Γ = 1/16 and the bath temperature is 0.01 in units of the hopping).Circles represent data from a Quantinuum model-H1-2 quantum computer, with error bars representing 2σ confidence intervals.The primary source of error in the implemented circuits is due to noise on the two-qubit gates [(2 − 3) × 10 −3 average infidelity].However, the resulting noise model leads to results that are barely distinguishable by eye from exact circuit simulations (black lines) and are omitted from the figure.
Fig.(S1a) ], while the remaining ones are components of time-ordered correlation functions on a more complicated contour with multiple branches [see Fig.(S1b)].
t e x i t s h a 1 _ b a s e 6 4 = " g u 8 9 O 9 u 3 2 Y u q 0 FIG. S1.(a) Simple Keldysh contour C made up of the two forward (C + ) and backward (C − ) branches.Black dots specify on which branch the three different operators appearing in the permutation ⟨O2O3O1⟩ are evaluated.(b) Multiple-branch Keldysh contour C2 needed for the evaluation of the permutation ⟨O3O1O2⟩.In this case, one cannot evaluate O3 on one of the first two branches and keep the product O3O1O2 ordered on the contour at the same time.Time-ordering along the Keldysh contour, regardless of the number of branches (or whether time increases or decreases), simply proceeds by ordering along the arc length of the contour from start to finish.
r w F r 9 f J u 1 6 z T u r 1 W / P K 4 2 r v J k i O S L H p E o 8 c k E a 5 I Y 0 S Y t w o s g T e S Y v z t x 5 d d 6 c 9 5 / R g p P v H J I / c D 6 + A a X 0 l F c = < / l a t e x i t > U t3,t2 K < l a t e x i t s h a 1 _ b a s e 6 4 = " U i R P U M P H E K B Q x H W e 9 l 0 x r H p Y w u w = " > A A A C D X i c b V D L S g N B E J y N r x h f U U / i Z T A I H i T s J o I e R S + C l w i u C s m 6 9 I 5 t H D L 7 Y K Z X k G X 6 W n a l A M q k d m n P V 9 i l e z + I V M c O w 6 y k 7 s g 3 s t c n 3 4 D u K W y X y m 7 F H Y P P E m 9 C y m y C q 3 b p y + 9 E I g a Z 7 b J 9 d s g 8 d s L O 2 C W 7 Y j U m 2 A N 7 Y s / s x X l 0 X p 0 3 5 / 2 n t e B M Z n b Y H z g f 3 + N a o A s = < / l a t e x i t > |+i < l a t e x i t s h a 1 _ b a s e 6 4 = " t e x i t s h a 1 _ b a s e 6 4 = " J Z I S x a F D D i G F 1 u m k m w P L y 4 k A r w F r 9 f J u 1 6 z T u r 1 W / P K 4 2 r v J k i O S L H p E o 8 c k E a 5 I Y 0 S Y t w o s g T e S Y v z t x 5 d d 6 c 9 5 / R g p P v H J I / c D 6 + A a X 0 l F c = < / l a t e x i t > U t3,t2 K < l a t e x i t s h a 1 _ b a s e 6 4 = " U i R P U M P H E K B Q x H W e 9 l 0 x r H p Y w u w = " > A A A C D X i c b V D L S g N B E J y N r x h f U U / i Z T A I H i T s J o I e R S + C l w i u C s m 6 9 I 5 t H D L 7 Y K Z X k G X

|0i< l a t e x i t s h a 1 _H S ↵ 2 <
FIG. S5.Quantum circuit needed in order to compute a contribution to the remainder function Rα arising from ρ (III) 3 and its transpose conjugate [see Eq. (S12) and the text beneath].