Quantum many-body Jarzynski equality and dissipative noise on a digital quantum computer

The quantum Jarzynski equality and the Crooks relation are fundamental laws connecting equilibrium processes with nonequilibrium fluctuations. They are promising tools to benchmark quantum devices and measure free energy differences. While they are well established theoretically and also experimental realizations for few-body systems already exist, their experimental validity in the quantum many-body regime has not been observed so far. Here, we present results for nonequilibrium protocols in systems with up to sixteen interacting degrees of freedom obtained on trapped ion and superconducting qubit quantum computers, which test the quantum Jarzynski equality and the Crooks relation in the many-body regime. To achieve this, we overcome present-day limitations in the preparation of thermal ensembles and in the measurement of work distributions on noisy intermediate-scale quantum devices. We discuss the accuracy to which the Jarzynski equality holds on different quantum computing platforms subject to platform-specific errors. The analysis reveals the validity of Jarzynski's equality in a regime with energy dissipation, compensated for by a fast unitary drive. This provides new insights for analyzing errors in many-body quantum simulators.


I. INTRODUCTION
More than a century after the foundations of thermodynamics and equilibrium statistical mechanics have been laid, the field of nonequilibrium physics remains an area of active research.The study of how macroscopic properties arise from microscopic quantum dynamics, initiated the development of quantum thermodynamics [1][2][3].Significant theorerical breakthroughs include the widely-applicable eigenstate thermalization hypothesis [4][5][6][7][8], and a handful of interesting cases of its violation [9][10][11]; at the same time, ongoing experimental progress made it possible to cool down highly controlled systems to temperatures dominated by quantum rather than thermal fluctuations [12][13][14][15].Despite this progress, we still lack a comprehensive theoretical framework or a complete set of principles to describe macroscopic phenomena in nonequilibrium physics.
A remarkable achievement in the field is the Jarzynski equality [16,17].It establishes a mathematical relation between the work  applied in a time-dependent process, and the free energy difference Δ between the initial thermal ensemble and the thermal ensemble associated with the final Hamiltonian: = 0, (1) where  =  −1 is the inverse temperature (we work in units  B = 1, ℏ = 1), and  diss =  −Δ is the dissipated work done during the process; the average ⟨•⟩  ( ) is performed over the work distribution () obtained from repeated applications of the protocol.Remarkably, Eq. (1) holds for any nonequilibrium protocol without restrictions; for close-to-equilibrium FIG. 1. Schematics of our protocol to test the quantum Jarzynski equality, cf.Eq. (1), in the quantum many-body regime.We simulate the dynamics of a system of interacting qubits initiated in a thermal state of the transverse-field Ising chain  i .The qubits then evolve under a non-equilibrium process on a quantum computer affected by energy dissipation.Finally, we extract the work distribution for the quantum circuit w.r.t. the final Hamiltonian  f .At the same time, we independently compute the exact theory prediction for the free energy difference.We compare both results against each other to test the validity of Jarzynski's equality [red box].
processes it reduces to the well-known fluctuation-dissipation theorem [16].
In the following, we focus on a formulation of the quantum Jarzynski equality valid for unital quantum channels, where work is defined by means of a two-point measurement scheme [28][29][30]: Let us denote the initial energy levels by  i  , and final energy levels -by  f  .In the following, we refer to a projective measurement in an eigenstate of the initial or final Hamiltonian, as an energy measurement; the measurement output is the corresponding eigenenergy  i  or  f  , respectively.The work is then defined as the difference between the measured final and initial energies, While this definition does not generalize to the most general open systems, together with recent progress in quantum simulation, it allowed for the experimental test of Eq. ( 1) using trapped ions [12,38], cold atoms [39], nuclear magnetic resonance (NMR) experiments [40], nitrogen-vacancy (NV) centers [41], and superconducting qubits [42].
Besides its fundamental importance in quantum statistical physics, testing the quantum Jarzynski equality is also of practical interest, as it allows us to measure free energy differences, which can be used, e.g., to characterize the onset of chemical reactions.In a recent study, Jarzynski's equality was used to extract approximate free energy differences in two-or three-qubit systems using minimally entangled typical thermal states [43].
In this work, we propose the idea that the quantum Jarzynski equality using a two-point measurement scheme provides a valuable benchmark for the performance of quantum devices.The only requirement for the equality to hold is the doubly-stochastic property of the underlying dynamics (see Appendix A and [29,38] for more details).This property is fulfilled for all unital channels, including pure unitary dynamics and dephasing noise.However, it does not hold for processes with energy dissipation.By experimentally testing Eq. ( 1), it is in principle possible to isolate contributions of processes violating double-stochasticity present during the dynamics and the measurement, which is of fundamental interest for improving current quantum technologies.
In the many-body regime (i.e., for many interacting degrees of freedom), the work distribution broadens with the square root of the system size [6], which requires an exponentially large number of projective measurements in order to estimate the left-hand-side of Eq. (1).Moreover, for quantum systems, further challenges arise: First, work fluctuations are not a measurable observable in quantum mechanics [25].Hence, testing the quantum Jarzynski equality presumes the ability to measure in the energy eigenbasis.This is notoriously difficult in practice since many-body energy eigenstates are typically volume-law entangled in real space.Second, preparing a quantum many-body system in a close-to-perfect thermal state can be demanding, and often comes with a substantial overhead of resource costs [46,47].For these and related reasons, the quantum many-body regime provides formidable challenges.
In the present study, we use classical presampling or midcircuit measurements to prepare thermal ensembles.In particular, mid-circuit measurements allow us to prepare a thermal ensemble of the transverse field Ising model [48] for up to  = 16 qubits.In contrast to many common approaches [46,47,[49][50][51], this enables us to prepare the thermal initial state of the transverse field Ising model without the overhead of using ancilla qubits.This innovation allows us to reach system sizes one order of magnitude larger compared to previous studies, which opens the possibility of using the Jarzynski equality as a benchmark for many-body quantum devices.
In a simulation on digital quantum computers, we test Jarzynski's equality for up to  = 16 qubits in a system of strongly interacting spins subject to a nonequilibrium protocol.We also analyze the quantum Crooks relation [52] [cf.Eq. ( 18)] for  = 8 qubits -an infinitesimal version of Jarzynski's equality -which had hitherto only been tested for a single two-level system [40,53].Unlike earlier experiments [12,38,40,43], we do not use a parametric quench, but a protocol of random entangling gates.
Moreover, we compare the accuracy of our results with the relaxation times on the quantum processor, and identify a novel feature in this regime of nonequilibrium qubit dynamics: we show that Jarzynski's equality holds approximately even in the presence of accumulating dissipation effects, so long as the execution time of gates is short compared to the thermal qubit relaxation time [so-called  1 ].In addition, we show that our nonequilibrium protocol improves the results in comparison to a pure energy dissipation process.
This paper is organized as follows: After introducing the challenges and their resolutions to test the quantum Jarzynski equality on a digital quantum computer in Sec.II, we provide results for the Jarzynski equality and the Crooks relation in Sec.III.In Sec.IV, we develop a theory to understand the results in the presence of energy dissipation.Finally, we compare our results with previous experiments in Sec.V and discuss the potentials of our work for benchmarking quantum devices and measuring free energy differences.

II. SIMULATION ON QUANTUM COMPUTERS
As already pointed out in the introduction, a quantum simulation of Jarzynski's equality in the many-body regime faces some restrictive challenges.
First and foremost, the work distribution requires the measurement of initial and final eigenenergies [25], cf.Eq. ( 2).At the moment, this is only feasible for a suitable choice of FIG. 2. The protocol of our quantum experiments (A) The unitary rotation gate defined in Eq. ( 13) and subsequent measurements prepare a thermal ensemble of the transverse field Ising model.The measurement has two effects: it prepares a thermal ensemble and is used as the initial energy measurement to determine the work distribution.(B) Transformation from the initial energy to the computational basis using a Bogoliubov (red) and a Fourier transformation (light blue).The gates are defined in Appendix C, the fermionic SWAP operators (dark green) keep track of the correct sign structure under permutation.(C) Now we apply the actual nonequilibrium protocol.In our case we execute a "random" circuit with  = 3 blocks.Each block consists of a layer of single qubit Haar random unitary gates (different colors) and a sequential CNOT layer.(D) Final measurement of the qubits in the computational basis corresponding to the final energy measurement of of Eq. (7).More details about the thermal state preparation and its accuracy in the quantum simulation can be found in App. C. the initial and final Hamiltonians, and requires the ability to measure in their respective eigenbases.
We emphasize that the ability to apply general unitary transformations is a distinctive feature of digital quantum computers.In contrast, this is currently not possible, in general, for analog quantum simulators, such as cold atom systems [54]; at present, this renders obtaining quantum work distributions in the many-body regime elusive on such platforms.
Second, the verification of Eq. ( 1) requires the preparation of an initial thermal ensemble.Most approaches for Gibbs state preparation require an overhead of ancilla qubits [46,47,[49][50][51].This uses up valuable qubits and makes the study of thermal states in the many-body regime difficult for the current generation of noisy intermediate-scale quantum (NISQ) devices.
Finally the number of required measurements increases with system size: According to [55,56], the number of shots  scales approximately at least as To get an estimate for ⟨ diss ⟩  ( ) , we can use the fluctuationdissipation relation [7] to obtain the realation Recall that the variance of generic non-adiabatic work distribution scales linearly with the system size  [6].As a consequence, the number of required shots scales at least exponentially in the system size: We stress that this scaling can even become worse in the case of generic nonequilibrium protocols which can create long tails in the work distribution [55,56].
In the following, we describe in detail how we tackle each of these challenges.In Sec.II A, we motivate our choice for the initial and final Hamiltonians.Section II B elucidates our choice of nonequilibrium protocol, which is summarized diagrammatically in Fig. 2.

A. Choice of initial and final Hamiltonian
Jarzynski's equality depends neither on the specific form of the initial or final Hamiltonians, nor on the protocol we evolve the state with.However, a suitable choice of these ingredients can help to address the above-mentioned challenges, and enables its verification through a simulation on a quantum device operated in the many-body regime.

Initial and final Hamiltonian
Notice that determining the exponentiated work distribution requires measuring the initial and final energy, i.e., two-point measurements [25].In general, this can be achieved by applying a unitary transform to switchfrom the computational (i.e., the Pauli -basis) into the energy eigenbasis, which is equivalent to diagonalizing the initial and final Hamiltonian, respectively.For generic systems, this requires a circuit of at least polynomial depth in the number of qubits , using a quantum phase estimation algorithm [57].However, for systems equivalent to free fermionic models, circuits increasing logarithmically with system size suffice [48].
Apart from this practical restriction, the validity of Eq. ( 1) imposes no further restrictions on the choice of initial and final Hamiltonians.Thus, we choose the initial Hamiltonian to be the transverse field Ising model with periodic boundary conditions [58], This system is integrable and can be mapped to free fermions through a Jordan-Wigner transformation.It can be diagonalized using a shallow circuit as described in Ref. 48.
Although  i is equivalent to a free-fermion model, its manyparticle eigenstates are entangled and, therefore, feature genuine quantum correlations; thus, the circuit dynamics in the many-body regime goes beyond existing work on single-and few-particle models previously analyzed in the context of verifying Jarzynski's equality [12,38,40,42,43,59].
As a final Hamiltonian, we chose the simple Hamiltonian which is already diagonal in the computational basis, such that measurements in the eigenbasis of  f are straightforward.
Overall, the choice of the initial and final Hamiltonians, which are either already diagonal in the computational basis, or where a circuit is known to diagonalize them in practice, enables us to perform a two-point measurement protocol to determine the work; this is essential to measure the left-hand side of the Jarzynski equality, cf.Eq. ( 1).The partition functions for the initial and final ensemble (and thus the free energy difference) can be computed exactly for our choice of Hamiltonians.The transverse field Ising model is equivalent to a system of non-interacting fermions with the energies   defined in App.C; the energy offset is given by  c = 1 − √ 2 (As mentioned above, we we work in units  B = 1, ℏ = 1).The free energy therefore reads as Since the energy offset   enters as a constant additive term both in the definition of the work  [Eq.(2)] and in the free energy  i above, and because the latter appear in the exponents on both sides of Eq. (1), we can ignore  c in the following discussion.The free energy of the final Hamiltonian is then In order to validate Eq. (1), it is therefore sufficient to focus on measuring the average exponentiated work.

Gibbs ensemble preparation: midcircuit measurements vs. classical presampling
To test the Jarzynski equality, we need to prepare a thermal ensemble w.r.t. the initial Hamiltonian on the quantum simulator.This is a difficult task on quantum devices in general since energy eigenstates are not necessarily accessible; however, the mapping of the Ising model to a model of non-interacting fermions allows us to use the initial energy measurement to prepare the thermal ensemble using a shallow circuit, as we now explain [see also parts (A) and (B) of the circuit in Fig. 2].
The thermal density matrix for a system of non-interacting fermions can be written as where   denotes again the energy of the single-particle eigenmodes.In the case of free-fermions, each energy eigenstate |⟩ is uniquely defined by the collection of excited singleparticle eigenmodes; This allows us to introduce an equivalence between bitstrings and eigenstates: We can map an occupied single-particle eigenmode  to the -th qubit in the excited state, written as |1⟩  .In the initial state, |0⟩ ⊗ , all qubits are in the |0⟩-state; thus, an excitation of the -th mode is equivalent to applying an -gate to qubit .Midcircuit measurements: starting from the product state |0⟩ ⊗ , we can implement the Gibbs ensemble in two steps: (i) to prepare a thermal ensemble in the free fermion basis, we start from the state |0⟩ and first apply a rotation gate   to each qubit, of the form with   implicitly defined by (ii) we then perform a subsequent measurement of all qubits.
The angles   are chosen such that this projective measurement collapses the state with probability equal to that of the Gibbs ensemble; hence, the measurement statistics correspond to sampling from the Gibbs state in Eq. ( 12).These two steps constitute part (A) of our circuit protocol, see Fig. 2.
To complete the procedure, we have to apply a transformation from the energy to the computational eigenbasis.This transformation consists of a fermionic Bogoliubov transform and a fermionic Fourier transform [cf.App.C]; these constitute part (B) of our circuit, see Fig. 2. In the case of  = 2  spins ( ∈ N), the Fourier transform can be decomposed into fermionic SWAP gates and two-body Fourier transform gates, which can be efficiently implemented on a quantum computer [48].
We emphasize that, besides its use for Gibbs state preparation, the mid-circuit measurement simultaneously serves as a measurement of  i  , which is necessary to determine the work distribution.
Classical presampling: unfortunately, mid-circuit measurements are currently not feasible on all present-day quantum computing platforms.Wherever they are not available, we use classical presampling [41,48,60].To this end, we prepare a randomly chosen many-body eigenstate |⟩ with probability equal to its Boltzmann weight, following Eq.(11).This is equivalent to sampling a bitstring directly from the Boltzmann distribution provided bitstring preparation can be performed with unit fidelity.
Since, at the end of the day, we want to perform a quantum simulation end-to-end, we avoid classical presampling whenever possible, and stick to midcircuit measurements for the largest system sizes  = 16.Furthermore, classical presampling requires preparing different circuits corresponding to each different initial bitstring, while the use of mid-circuit measurements requires only one circuit for all possible initial bitstrings.As a consequence, this results in a speed-up for the execution of the simulation.
After having discussed the difficulties related to measuring work and thermal state preparation, we now move on to address the remaining challenges concerning the exponential scaling of measurement shots and free energy differences.

B. Nonequilibrium protocol
In contrast to previous experiments, the local control over qubit interactions offers significant freedom in the choice of nonequilibrium protocol, cf.part (C) of our circuit in Fig. 2. Our goal is to devise a protocol that makes use of the intrinsic features of quantum computers in order to address the exponential scaling of the number of required measurements with the system size.
Because of the local nature of gates on digital quantum computers, we choose a circuit whose dynamics does not describe a parametric deformation between the initial and final Hamiltonians.This allows us to explore nonequilibrium protocols distinct from previous experiments [12,40].
Instead, in our digital quantum simulations, we apply a protocol made out of  sequential blocks, as shown in Fig. 2(C).Each block consists of a layer of single-qubit Haar random unitary gates, followed by a sequential layer of CNOT-gates.In comparison to pure Haar random circuits, our circuit allows for a more native and shallow implementation on NISQ devices; moreover,  =  − 1 blocks are sufficient to build up bipartite von Neumann entanglement in the system close to the Page value [61], as we demonstrate in a classical emulation in App.D.
Our choice of a non-equilibrium protocol results in a work distribution (), that approaches a Gaussian for system sizes  ≳ 8 [cf.App.D].Thus, the tails of the work distribution exhibit a Gaussian decay and do not dominate the average of the exponentiated work [55,56].
However, the number of measurements to determine the distribution of the exponentiated work still scales exponentially with the system size  and thus eventually gives rise to a bottleneck in testing Eq. ( 1).In our case, we find that for the case of  = 16, 2 16 measurements are sufficient.As we show below, this is feasible on current devices.
In order to determine the final energy (and with it, the work), we perform up to 2 24 measurements on superconducting quantum computing architectures, and between 2 12 (Quantinuum) and 2 16 (IonQ) measurements on trapped ion platforms.The reason for the smaller number of measurements on trapped ion devices is the higher simulation cost of the quantum simulation due to lower protocol repetition rates, which restricts us to use less data points.For further details, see App.B.
Last, it is curious to note that the separation of parts (B) and (C) in our protocol [Fig.2] is somewhat arbitrary.On the one hand, part (B) belongs naturally together with part (A) in a thermal state preparation subprotocol.On the other hand, we can interpret part (B) as part of the nonequilibrium protocol (C) applied to the system.Since Jarzynski's equality holds for arbitrary protocols, equilibrium or nonequilibrium, the accuracy of implementation of part (B) is not crucial for the accuracy to which we verify Jarzynski's equality, provided the errors are systematic and are, thus, identical across trials.Note that this does not compromise the many-body character of our circuit: indeed, the entanglement created by the circuit reaches the maximal Page curve, cf.Fig. 9 in App.D.

III. JARZYNSKI EQUALITY AND CROOKS RELATION ON A NOISY QUANTUM DEVICE
After having discussed the details of the implementation of a test for the many-body quantum Jarzynski equality on quantum devices, we now provide simulation results for the test of Jarzynski's equality and Crooks' relation on presentday quantum computing platforms.

A. Jarzynski relation in the few-and many-body regime
The results presented below are obtained in experiments on the ibmq_guadalupe superconducting quantum processor [62].A detailed comparison with other digital quantum architectures is included in App.G.For the technical details of the various devices including gate fidelities and relaxation times  1 and  2 , we refer the interested reader to App.B 3.
In Fig. 3 we display the deviations from Jarzynski's equality as a function of the inverse temperature.The three panels show the same data, but viewed through the lens of different quantifiers of the deviation.
Let us open up the discussion in Fig. 3 (a) by introducing the plain deviation quantifier ln e − diss  ( ) , which vanishes whenever Jarzynski's equality holds.The curves show an approximately linear scaling with the inverse temperature  in the high-temperature limit.The dominant contributions to this deviation quantifier originate from violations of doublestochasticity with equal contributions from all final energy eigenstates.By contrast, for inverse temperatures  ≳ 1, the deviations are determined by processes involving the ground state of the final system.As a consequence, the deviations for large  ≳ 1 will converge to a constant value, which depends on the concrete effect of energy dissipation on the process.For a detailed discussion, see App.E 2. The comparison between the three system sizes in Fig. 3 (a) reveals increasing absolute deviations from Jarzynski's equality with increasing .We suspect that this is due to the increasing number of qubits to be read out: a primary contribution to the overall error arises from measurement errors, and hence the error in the work measurement scales linearly with system size.We provide a more detailed analysis of various errors in Sec.IV and App.E. For a meaningful and systematic comparison of the results for different system sizes we, therefore, consider different quantifiers for relative deviations in panels (b) and (c).
Figure 3 (b) shows the ratio Δ sim /Δ, introduced in previous work [12,42].Here, is the free energy difference obtained from the distribution average of the exponentiated work and Δ is the theoretical FIG. 3. Testing Jarzynski's equality on a digital quantum computer using midcircuit measurements as a function of inverse temperature , for three different systems sizes  = 4, 8, 16.(a) ⟨e − diss ⟩  ( ) as a function of inverse temperature .The unnormalized data depend on temperature and sytem size: deviations from Jarzynski's equality for a dissipation-free system are sensitive to the system size  and grow with increasing inverse temperature.(b) Δ sim /Δ as a function of inverse temperature .The normalization largely removes the dependence on temperature and sytem size.(c) The same data, but now normalized by the maximum possible deviation [see Eq. ( 17)].The dashed red line shows the prediction from our theory derived in Sec.IV.Jarzynski's equality holds better than one part in ten, irrespective of the system size; the deviation agrees well with our theoretical prediction.We took 2 16 measurements on ibmq_guadalupe; the number of blocks in the circuit [cf.Fig. 2] is seven for  = 4, 8, and three for  = 16.Further technical details of the device can be found in App.B 3.
prediction for the free energy difference.Any deviations from the ideal, dissipation-free case, are indicated by deviations of the ratio from unity.We observe Δ sim /Δ > 90% for all inverse temperatures and only a weak system-size dependence, since the latter is absorbed by the scaling of the denominator.
Since Δ is a protocol-dependent quantity and not immediately accessible for a given setting, we introduce another normalization, which facilitates a straightforward quantitative comparison with previous experiments and which is tailored to quantify the amount of energy dissipation in the system.For this purpose, we consider the process of dissipative decay to the ground state of the final system as the natural reference for deviations from the Jarzynski equality, because it constitutes its worst possible violation.The resulting work obtained by the two-point measurement scheme when starting from an initial energy  i is  =  f 0 −  i .Accordingly, we introduce where the right-hand side is the ratio of the exponentiated energy differences for the purely dissipative process, and Δ sim for the true process simulated on the quantum computer.The bar (•) denotes an average over the initial thermal ensemble.This allows us to define the relative deviation from the theoretical prediction, This quantity is bounded from below and above.Whenever Jarzynski's equality holds, we have E () = 0. On the other hand, for a purely dissipative process, E () = 1.In contrast to previous deviation quantifiers, the normalization in E () by the worst-case scenario allows us to directly compare the amount of dissipation in our simulations and previous Relative deviation from Jarzynski's equality as a function of the waiting time  idle (solid lines), for four different inverse temperatures , at  = 8.The -axis is normalized by the average waiting time  1 = 118 s [63].For comparison, the diamonds denote the results for our protocol shown in Fig. 3(b), placed at their respective execution time [dashed vertical line].The application of gates between the two measurements in the protocol [cf.Fig. 2] improves the accuracy of our results by up to more than one order of magnitude.experiments, cf.Table I in the Discussion.Moreover, small deviations from the equality can be resolved logarithmically, while the upper bound simultaneously justifies a quantitative interpretation.
The relative deviation from Jarzynski's equality as defined in Eq. ( 17) is plotted in Fig. 3 (c).Consistently with Fig. 3 (b) the results are largely independent of the inverse temperature  and the relative deviations do not exceed 10%.As for the previous normalization Δ sim /Δ, E () exhibits only a weak system-size dependence.Because increasing the system size requires a larger circuit depth via the sub-circuits (B) and (C) in Fig. 2, one could naïvely expect that scaling to large system sizes is quickly hampered by dissipation effects which accumulate with increasing circuit size.However, this appears not to be the case in the observed behavior in Fig. 3, where we show data up to  = 16 qubits.Therefore, we now briefly investigate our circuit's susceptibility to energy dissipation.
The purely dissipative process leading to the right-hand-side of Eq. ( 16) can be emulated on a NISQ device by applying a so-called "idle process" [64]: this is a similar process as the one described in Fig. 2, but with the circuit parts (B) and (C) replaced by free evolution for a variable duration  idle .In other words, the two measurements at the end of circuit part (A) and in part (D) are separated by the idle time  idle .Hence, to access the regime of validity of Eq. ( 16), we have to apply an idle process of time  idle ≫  1 , with  1 being the dissipation time.
Figure 4 shows the relative deviation from Jarzynski's equality for an idle process, as a function of the waiting time  idle for  = 8 qubits.As expected, for large waiting times  idle ≳  1 , E → 1 reaches the limit of a pure dissipative process.We now want to compare the time required to reach the purely dissipative regime with the execution time of the circuit from Fig. 2. Before we do this, notice first that the average dissipation time of ibmq_guadalupe is  1 ≈ 118 s, while the execution time for our circuit [Fig.2, (A) to (D)] is   = 57 s for  = 8 qubits; hence   / 1 ≈ 0.48 [cf.Diamonds on dashed vertical line in Fig. 4] [for comparison, for 16 qubits we have   / 1 ≈ 2.95].Furthermore, the normalized deviation of the quantum Jarzynski equality Eq. ( 17) does not depend on the number of entangling blocks, as is shown in Fig. 12 in App.D 4.
Comparing the deviation values with the idle process, the accuracy of the process from Fig. 3 appears quite striking, since the deviation there is up to more than one order of magnitude smaller than an idle protocol of the same waiting time,  idle =   .In fact, the deviation values are comparable to those of the idle process with the smallest idle time investigated, from which we can deduce that the main source of deviations from Eq. (1) originates in part (D) of the protocol.
To sum up, our nonequilibrium circuit has the property of preventing energy dissipation effects from accumulating, despite increasing circuit execution time; this is, in turn, reflected in Jarzynski's equality being obeyed to a high degree of accuracy even in the many-body regime of  = 16 qubits.We will discuss this observation in more detail in Sec.IV.

B. Crooks relation in the many-body regime
While Jarzynski's equality is a statement about workaveraged processes, we can also check to what extent its infinitesimal [or in-sample] version, the Crooks relation [29,52] ln holds [see App.A for a derivation].To the best of our knowledge, the latter was only tested in two-level quantum systems so far [40,53].For a process described by a sequence of unitary gates, and ignoring noise for the time being, the backward process can be implemented in a straightforward way by reversing the gate order and taking the inverse of each individual gate.
We tested the Crooks relation for  = 8 qubits, using 2 24 ≈ 1.7 × 10 7 measurements to reduce statistical fluctuations due to a limited number of measurements.The corresponding data are shown in Fig. 5.The validity of the Crooks relation requires that all data points lie on the black line, apart from statistical fluctuations arising from a finite number of measurements.Indeed, the data follow the trend predicted by the Crooks relation.However, the deviation between the theory prediction and the measurement outcomes cannot solely be explained by the statistical uncertainty due to a limited number of measurements: The ratio between occurrences of forward and backward processes in Eq. ( A3) is up to one order of magnitude larger than expected due to statistical errors.These deviations are a consequence of the presence of energy dissipation in the device.

IV. THEORETICAL ANALYSIS
In order to understand the small deviations in the validation measurement for the Jarzynski equality and the strong fluctuations in the Crooks relation, we analyze a simplified single-qubit toy model subject to a periodic application of a single-qubit rotation gate, instead of a random sequence of gates.As we discuss below, the main results of this analysis extrapolate also to the protocol from Fig. 2.
Consider the single-qubit density matrix : We analyze the repeated [periodic] application to the qubit of a unitary gate  of the form The density matrix vector evolves according to Let us assume that the excited state has a finite lifetime  1 ; thus, we can use a single amplitude damping channel [65] to describe this effect.The time evolution over one period is then given by Here F  is a completely positive trace-preserving map on the space of density matrices with a unique largest eigenvalue  = 1 for  > 0. The Kraus operators K  describing the damping process are given by The parameter  is the fraction of the excited state population which decays to the ground state after one application of F .It is related to the relaxation time  1 and the gate time Repeating the process  times, the density matrix evolves in the long time limit as lim A simple calculation gives, for  > 0, with and Jarzynski's equality holds for  ( , ) = 0, which is indeed the case at  = 0. Our quantum simulations testing Jarzynski's equality are in the regime   / 1 ≪ 1.In order to understand this regime, we perform an expansion of  ( , ) around  = 0; together with Eq. ( 24) this gives This regime is approached in the case of an infinitely fast drive, as can be seen from the relation in Eq. ( 29).Note that the actual time of the protocol does not matter in this case, since this analysis holds in the limit of infinitely many periods  → ∞.Indeed, there is no obstruction if the protocol time greatly exceeds the relaxation time scale  1 .

B. Extension to multiple qubits
Although we restricted the analysis to the case of a single qubit with a periodic protocol, we can apply it to random gates and multi-qubit systems.For simplicity, let us only take into account the independent decay of qubits in the excited state to their ground state, and ignore many-particle effects for a moment.Although this is a crude simplification, it will turn out that it gives quantitative good results, such that this approximation is well justified a posteriori.The relevant physical parameter is , as defined in Eq. (24).The application of single-qubit gates can be interpreted as a drive applied to the system, which repopulates excited states and thus compensates for the deviations from Jarzynski's equality caused by the energy dissipation process.
To make a quantitative estimate for our experiment, we determine the parameter  and the measurement error of ibmq_guadalupe.The execution time of single-qubit gates is negligible in comparison to two-qubit gates   , so  is given by: see App.B 3 for the characteristic time scales of the device.Thus, in our approximation, the effect of 2-qubit entangling gates only enters via the gate time   .Furthermore, averaging over the angle  for each single rotation gate,  ( , ) reduces to In addition, we have to factor in the dissipation during the measurement process, i.e., excited qubits which are misidentified to be in the ground state.In the following, we neglect the opposite case, i.e., a qubit in the ground state misidentified to be in the excited state.For the device ibm_guadalupe, we get an approximate measurement error of  0 ≈ 0.04.The steadystate probability for each single qubit to be in the excited state is then given by (32) Since  0 ( ) ≪  0 , the model predicts that the deviations from our experiment are almost exclusively due to measurement imperfections.The analytically estimated deviation from a dissipation-free evolution is indicated as a red dashed line in Fig. 3. Taking the simplifications of our analysis into account, the result agrees quantitatively well with our simulation results and suggests that measurement errors are indeed the dominant source for deviations from the quantum Jarzynski equality.

C. Relation to DiVincenzo's third criterion
The compensation for energy dissipation is reminiscent of DiVicenzo's criteria for quantum computing [66].To be more concrete, the third criterion states that scalable quantum computing requires long decoherence times in comparison to the time scale of operational gates.
In our case, we are not interested in scalable fault-tolerant quantum computing, but in a weaker question, namely under which conditions the quantum Jarzysnki equality Eq. ( 1) is fulfilled.This only requires long dissipation timescales.In contrast to fault-tolerant quantum error correction, it is not necessary to actively correct for errors; the dynamics itself already compensates for the energy dissipation.Put differently, the application of single-qubit gates can be interpreted as a re-population of excited states.The quantum Jarzynski equality, thus, holds in an effective way, although every single component of the dynamics is spoiled by energy dissipation.

V. DISCUSSION & OUTLOOK
We proposed a protocol to test the quantum Jarzynski equality and the Crooks relation in the many-body regime on nearterm quantum computing devices, in the presence of different errors.We pushed the state of the art for the quantum simulation of Jarzynski's equality up to 16 qubits, and for the Crooks relation -to 8 qubits, respectively.We identified the implementation of two-point work measurements, the preparation of a thermal ensemble on a digital quantum simulator, and the exponential growth of the required number of measurements as the major practical challenges to reach the many-body regime.
To address the first two challenges, we developed a protocol sequence that prepares a canonical ensemble by using mid-circuit measurements.While for small system sizes classical presampling is still feasible, mid-circuit measurements allow us to prepare thermal ensembles for up to  = 16 qubits without running multiple independent experiments and without any overhead incurred by using ancilla qubits.Taking the transverse field Ising model and an   -model as exemplary initial and final Hamiltonians, respectively, we performed energy measurements with circuits scaling at most logarithmically in the system size .This approach is exact for Hamiltonians equivalent to single-particle systems; it is currently an open question whether, to what accuracy, and under which conditions, one could prepare thermal ensembles for more general Hamiltonians using similar protocols.A possible route could be to use variation quantum eigensolvers for the preparation of thermal ensembles for non-integrable systems [67].
While the exponential (in the system size) number of projective measurements still appears as the dominant bottleneck in testing the quantum Jarzynski equality when increasing the number of degrees of freedom, we were able to collect enough measurement data to test the latter for up to  = 16 qubits.As a side product, our protocol in Fig. 2 reveals the ingredients of Eq. ( 1) in a simple way: the quantum Jarzynski equality is a relation between initial and final eigenenergies, and a doublestochastic transfer matrix connecting them.Any other terms, including the transformation to a computational basis, can be absorbed into the transition process itself.The quantum Jarzynski equality thus probes only the double-stochasticity of the process and not the accuracy of different parts of the protocol.
Let us also compare our results to previous experiments on the quantum Jarzynski equality, cf.Table I.To this end, we first extract the results for the left-hand side of Eq. ( 1), and the theory prediction of the free energy of the corresponding experiments; then we compute the maximum deviation of an idle process, cf.Eq. ( 16), which gives us the normalized deviation defined in Eq. ( 17) for each experiment [cf.Appendix G for details].The comparison shown in Table I clearly suggests that the accuracy of our results is comparable with most of the earlier experiments; however, we reach an order of magnitude larger system sizes, where quantum many-body effects become pronounced.
Moreover, in contrast to previous experiments, the protocol duration of our circuits is comparable to, or even exceeds, the average dissipation time  1 of the NISQ devices; therefore, energy dissipation is no longer negligible.We checked that our results do not depend on the specific choice of randomness in our protocol.Furthermore, we demonstrated that the relative accuracy of our results is almost independent of system size and circuit depth.We also developed a theoretical model which predicts the empirical observations in a quantitative manner.By employing a fast drive that compensates for energy dissipation, we thus found the Jarzynski equality to be effectively valid in this regime, even though the dynamics is not doubly stochastic.
While these observations are interesting on their own, our work demonstrates two promising practical applications: First, testing the equality can be used to investigate errors on NISQ devices in new ways: Since Jarzynski's equality is only sensitive to processes violating double-stochasticity, it can be used to quantify and single out this effect.Even though simpler protocols to determine the dissipation time  1 of devices exist [65], our approach can be used as the generalization of these approaches to the many-qubit regime: Jarzynski's equality can not only be used to detect the decay of excited states as shown in our analysis in Sec.IV, but it is, more generally, sensitive to any violation of double-stochasticity.It remains an open question for future studies to investigate to what extent this allows us to refine our understanding of correlated error processes on modern quantum devices.For instance, the commonly used protocol for measuring  1 can be interpreted as special cases of testing the Jarzynski equality for a single spin system, as we show explicitly in App.F. Second, our analysis emphasizes important limitations concerning the measurement of free energy differences -one of the promising applications of the Jarzynski equality: in generic non-integrable systems, the exact work distributions can only be extracted by diagonalizing the circuits, which requires at least polynomially deep circuits using quantum phase estimation [68] and is thus infeasible in the many-body regime.It remains an exciting question for future research to find approximations to Eq. ( 1), which allow for a scalable method to extract free energies using the quantum Jarzynski equality in the many-body regime [43].

Appendix A: Quantum Jarzynski Equality and Crooks Relation
Let us recall the derivation of the quantum Jarzynski equality [28,30].Consider a system described by a Hamiltonian  i with eigenstates  i | i ⟩= i  | i ⟩, coupled to a thermal reservoir of inverse temperature .The system is thus described by a thermal ensemble with the density matrix is the partition function.We now decouple the system from the reservoir, and let it evolve according to a dynamical process  (not necessarily unitary, see below).At the end of this process, the instantaneous final Hamiltonian  f of the system has eigenstates We denote by  → the transition probability from the initial eigenstate | i ⟩ to the final eigenstate | f ⟩.
While the dynamical process  need not be unitary, we require that the transition probabilities satisfy the following two sum rules: The left-hand equality reflects the conservation of probability.The right-hand sum rule is less obvious and implies the socalled double-stochasticity of the matrix   ; this condition is fulfilled for unitary dynamics  → = |⟨ f || i ⟩| 2 , but is also conserved throughout evolution in the presence of additional decoherence noise [38].By contrast, energy dissipation violates the right-hand sum rule condition [38,69].
Using these definitions, we can now prove the quantum Jarzynski equality for the process introduced above: where we used   = f  − i  according to Eq. ( 2) and the double-stochasticity of  → in the second line, and the definition of free energy  i = − ln  i in the last line.We note that Jarzynski's equality imposes no restrictions on the initial and final Hamiltonians; in particular, they need not be identical.

NISQ architecture characteristics
In order to test the effect of different noise types, we run our circuits on five different devices using two different architectures: superconducting qubits (ibm_perth, ibmq_guadalupe and Rigetti Aspen-11), and trapped ion platforms (Quantinuum H1 and an 11-qubit system of IonQ).We extract the exponential of the work, Eq. ( 2), to test the Jarzynski equality, cf.Eq. ( 1).As mentioned above the latter is valid also in the presence of noise which does not violate double-stochasticity, and is only sensitive to errors that violate the second sum rule in Eq. (A1).The size of the deviation from the theoretical prediction Eq. ( 1), valid for an ideal dissipationless device, gives us therefore information about the amount processes violating double-stochasticity during the simulation.
The quality of qubits is often measured by means of the average gate times   and the relaxation times  1 and  2 .The timescale of dephasing errors,  2 , is not relevant for our purposes, since depolarizing errors do not violate doublestochasticity [30,38].Only the thermal relaxation time  1 or, more concretely, the ratio  =   / 1 matters, as it sets the decay rate for excited states, cf.Sec.III.While the two-qubit fidelities for all architectures fall between 95% and 99.5%, the -factor depends strongly on the underlying architecture.
As discussed in Sec.IV, the accuracy to which the Jarzynski equality Eq. ( 1) holds in experimental setups, depends on (a) the ratio between gate time   and the relaxation time  1 , (b) measurement errors, and (c) statistical errors due to a finite number of measurements.The statistical error reduces with the square root of the number of measurements.Thus, by testing to what accuracy Eq. ( 1) holds, we gain information about processes violating double-stochasticity in the quantum device.FIG. 6. Validation of the Jarzynski equality Eq. ( 1) for  = 4 qubits as a function of inverse temperature  on ibm_perth (blue), Rigetti Aspen-11 (orange) and IonQ (green), 65, 536 shots with classical sampling.Red: comparison with a circuit executed on Quantinuum H1 with a slightly different sampling, see text, and 4, 000 shots.The data for superconducting platforms are denoted by circles, for trapped ions -by squares.The trapped ion platforms show a better performance than the superconducting qubit architectures, since ln⟨e − diss ⟩ is closest to zero -the theory prediction value.
The two-qubit gate time for IonQ-devices is   ∼ 200 s, and the relaxation time  1 ∼ 10 7 s, resulting in a factor  IonQ ∼ O (10 5 ).On the other hand, the timescale for IBMplatforms are   ∼ 400 ns,  1 ∼ 160 s, yielding  IBM ∼ 400.Based on these estimates, we anticipate obtaining significantly smaller deviations from Jarzynski's equality Eq. ( 1) on trapped ion compared to superconducting platforms.Furthermore, due to the high-quality factor  IonQ we can assign deviations from the Jarzynski equality on trapped ion platforms solely to measurement errors.

Comparison on different devices
Besides presenting results for the largest system sizes and the most resource-intensive simulations (see main text), it is worth comparing the performance of the different devices on smaller systems.This is necessary to obtain an overview of the scaling with system size  of the simulations, and their expected resource costs.Furthermore, it allows us to get valuable insights about the quality of the experimental platforms.
For a direct comparison of the performance on the five different devices, we estimate the validity of Jarzynski's equality using 2 16 measurement shots and  = 4, 8 qubits.In the case of Quantinuum, we are restricted by the smaller number of measurements to  = 4 qubits and to a different sampling scheme due to the higher financial resource costs for circuit evaluation on these devices: We sample each eigenstate 200 times, and weigh the results by the Boltzmann distribution during classical postprocessing.In this way, we are able to extract data for different inverse temperatures with a total of 4, 000 shots.II.Accuracy of the Jarzynski equality for  = 0.7, executed on different devices for an extension of the protocol from Fig. 2 to  = 8 qubits.The case for mid-circuit measurements is shaded in grey.The values for the runtimes on IonQ and Rigetti are extrapolated, using information about the compiled circuit and gate times.The accuracy of the results on ibmq_guadalupe is comparable to that on the trapped ion device IonQ.
Because of hardware limitations on Rigetti and IonQ, that prevent the use of mid-circuit measurements, in order to make a fair comparison between the various devices, in the following we use classical presampling for the preparation of the Gibbs ensemble.
We emphasize that the comparison in this section is not a statement about the quality of the different devices for an accurate simulation of a quantum circuit: as explained in Sec.A, the quantum Jarzynski equality is only sensitive to noise channels violating double-stochasticity, and thus a good accuracy for the validation of the circuit does not imply an accurate evaluation of the corresponding quantum circuit.
Accuracy of the results.-Theresults for  = 4 qubits are shown in Fig. 6 for five different values of the inverse temperature ; for  = 8 and a fixed  = 0.7, see Table II.The error bars indicate the statistical uncertainty due to a finite number of measurements.
We apply no error mitigation; instead, we consider the measurement process itself as part of the protocol with its own errors, which, depending on their origin, may or may not violate the double stochasticity condition.Moreover, in the case of IonQ, error mitigation was infeasible, since the data collection continued over a few days.However, we can still obtain rough estimates for the measurement errors and their impact on the validity of Jarzynski's equality.
The errors for the IonQ, Rigetti and IBM devices are within 2% to 5%, even for large values of  and  = 8 qubits.Note that the theory predicted value from Eq. ( 1), ln ⟨e − diss ⟩  ( ) = 0, does not fall within the error bars for the superconducting architectures-a direct manifestation of violations of doublestochasticity.
Even if deviations from Jarzynski's equality are not detectable with the current experiment, a closer look reveals that the measurement process itself is an energy-dissipative process in this case: our quantum simulations reveal that excited qubits are detected in the ground state with a probability of roughly 1%.This process violates the second sum rule in Eq. (A1), and thus leads to a weak violation of the quantum Jarzynski equality.As discussed in App.E, our simulations reveal a means to detect the readout error as a residual deviation in the limit of infinite measurements.
It is also insightful to compare the deviation from Eq. ( 1) with and without error-prone mid-circuit measurements on a IBM quantum device.In the case of  = 8 qubits, the error is increased by almost a factor of 2 due to the first measurement, as is shown in Table II.This already shows that most of the deviation from the theoretical prediction is due to the measurement process.
On Rigetti's Aspen-11, we found it was essential to disable "fencing" to obtain quality results on par with other devices.Fencing makes two-qubit gates executed sequentially, even when acting on different qubits within the same circuit layer.While this reduces crosstalks and increases the fidelity of the individual operations, it leads to circuits with a longer execution time, making relaxation effects through  1 more prominent.
Runtime and resource cost.-Regarding the execution time, note that different physical platforms have different run times.In the case of superconducting qubits, gates are implemented as a sequence of microwave pulses, where the average duration of each such gate is of the order of 10 ns to 100 ns [73]; thus, the overall circuit duration for  = 8 qubits is of the order of 50 s.In the case of trapped ions, gates are implemented via twophoton Raman processes, with gate times ranging from 10 s to 100 s [70].As a consequence, the execution time for a circuit is also three orders of magnitude longer, which presents a relevant bottleneck for us when increasing the system size or the number of measurement shots.
The use of different architectures also affects the simulation cost.The simulation cost on trapped ion systems is at least 30 times higher than on superconducting qubits, which is caused by the longer circuit evaluation times.This is the major bottleneck we encounter for simulations on trapped ion devices: it limits the number of circuit evaluations to 2 16 and the system size to  = 8 qubits on IonQ, and to 4, 000 evaluations on  = 4 qubit systems on Quantinuum H1, respectively.
Ability to perform midcircuit measurements-In the end, we want to perform a quantum simulation without any classical presampling.The thermal state preparation introduced in Sec.II A 2 requires the ability to perform midcircuit measurements.From the above platforms, only the IBM and Quantinuum devices are currently capable of performing this task.

Technical data of the various quantum devices
The technical details, including coherence times, thermalization times, and gate times, as well as the average two-qubit fidelities on the different devices we used, are shown in Table III.

Appendix C: Preparation of a thermal distribution for the transverse field Ising model
In the following section, it is explained how the eigenstates of the transverse field Ising model are mapped to computational basis states of the quantum simulator by a shallow circuit.This step is crucial to test the quantum Jarzynski equality on current devices.denotes the relaxation rate of off-diagonal elements [65].  denotes the average execution time of the native two-qubit gate on each of these devices.The average two-qubit fidelity is the average fidelity of the system-specific two-qubit gate. denotes the ratio between  1 and gate time   , as introduced in App.B. While the qubit fidelity is comparable for all devices,  1 varies over several orders of magnitude, depending on the underlying architecture.

Theory
As we discussed in the main text, the transverse field Ising model can be mapped to a noninteracting fermionic Hamiltonian.It is therefore possible to prepare it in a Gibbs ensemble, using the protocol described in Sec.II A 2, and with the additional help of a unitary transformation between the energy and spin bases of the system.The different transformation steps are explained here in detail, following Ref.48.
The transverse field Ising model is given by the Hamiltonian The transverse field here is chosen of the same strength as the Ising interaction, although this is not a necessary requirement for testing Jarzynski's equality using our protocol.In order to simplify the realization on a quantum computer, we impose periodic boundary conditions, and add an additional Pauli string to eliminate unwanted terms that appear in the Jordan-Wigner transformation: Note that the multi-body term becomes negligible in the thermodynamic limit.As a first step, we transform the Hamiltonian into fermionic modes using a Jordan-Wigner transform: This gives the following fermionic Hamiltonian: The wave function can be expressed as  13) and ( 14).The ensemble arising from repeated measurements is described by a Gibbs state.After the measurement, the circuit realizes a transformation from the energy eigenbasis of the transverse field Ising model to the computational basis.Here   are Fourier gates, and   -Boguliobov gates.Note that fermionic SWAP gates are required, that contain an extra sign compared to qubit swap gates.
Here |Ω  ⟩ is the vacuum state, i.e.,   |Ω  ⟩ = 0. Note that the coefficients of the wave function do not change; thus, the Jordan-Wigner transformation does not add any additional gates to the quantum circuit.However, we have to keep track of fermionic signs when swapping fermionic modes.The next step is to apply a Fourier transform: .

(C6)
The Fourier transform for  = 2  ( ∈ N) qubits can be implemented with a quantum circuit of depth log().To see how, we split the Fourier transform into even and odd sites: The two terms on the right-hand-side independently represent a Fourier transform for /2 fermions.The case  = 2  is particularly appealing, since we can keep iterating this step until we end with a Fourier transform of only two fermions, which can be easily implemented using two-qubit gates.
To do so, we need the fermionic swap gate (note the sign  structure, green gates in Fig. 7) and the Fourier gates (light blue gates in Fig. 7) In our case, we have to restrict to  ≤ 16, due to a limited number of available qubits.The above transformations lead to the Hamiltonian Finally, in order to diagonalize the Hamiltonian, we have to apply a Bogoliubov transformation: This transformation can be achieved with gates (red gate in Fig. 7) of the form where [48] The transformation steps for  ≥ 4 qubits are analogous to the case described above.This casts the Hamiltonian in diagonal form: with eigenenergies Here,   is a constant energy offset   = 1 − √ 2. Note that we ignore this term in the work distributions and free energy computations, since it only appears as a constant factor in the partition function  i , which is not relevant for testing the validity of Eq. (1).
In order to prepare the Gibbs state, we consider the diagonalized Hamiltonian.In this case a thermal ensemble can be prepared using projective measurements, as explained in Sec.II A 2. To apply a transformation back into the computational basis, we have to reverse the unitary Fourier and Bogoliubov transformations described above.The corresponding circuit is shown in Fig. 7.The data indicate that our chosen protocol is away from the adiabatic regime, for which ⟨ diss ⟩=0.

Accuracy of the thermal state preparation using midcircuit measurement
Let us now investigate the accuracy of the mid-circuit measurement state preparation.To do so, we detect the probability to prepare a state | i ⟩ after part (A) of the protocol displayed in Fig. 2, and compare it with the probability distribution of the canonical ensemble.The results are shown in Fig. 8, using the data for  = 8 qubits from Fig. 3.The quantum simulation prepares the correct thermal ensemble for different inverse temperatures  and the entire range of initial energies.Note that we observe a slight temperature dependence in the fluctuations: Low energy states are prepared more often than predicted by theory; higher energy states are slightly underrepresented.As we discussed in Sec.IV, the measurement itself is a dissipative process, causing decay from the excited to the ground state of the qubits.
In general, device imperfections do not allow us to perfectly prepare a given target state.In our case, we obtain instead of an eigenstate |  ⟩ of the transverse field Ising model, a density matrix   sim .To determine the quality of the state preparation process, we obtain the form of the density matrix for  = 4 qubits using quantum state tomography [68].This allows us to compute the average single particle fidelity As before we apply no error mitigation in this case.

Appendix D: Analysis of the nonequilibrium protocols
This section is dedicated to an analysis of the specific circuit protocols we selected to use in this study.We show here numerically that we are operating in a nontrivial quantum many-body regime, by computing the operator entanglement entropy of our circuits [74], and the work distribution they give rise to.Finally, we present results for different choices of the one-body random unitary gates and show that the specific choice of random gates has only a minor impact on the accuracy for validation of Eq. (1).

Operator entanglement entropy
We analyze our circuits from an entanglement perspective and show that our ideal circuits are sufficient to create entan-glement close to the Page value.This, in turn, demonstrates that we operate in the quantum many-body regime.
For a given system, we can choose complete basis sets of operators {   } and {  } which are orthonormal and have only support on subsystem  and its complement  =   , respectively.An In the following, we consider subsets of the form  = {0, . . .}, where  is the position of the last qubit included in the subset.The operator entanglement entropy for such a partition is defined as [74] We now compute the operator entanglement entropy of the unitary operator of our protocol for different partitions.The results are shown in Fig. 9.It is clear that our chosen operators with  − 1 blocks already exhibit an entanglement (operator) entropy close to the Page-value [61],  Page =  ln 2 − 1/2.

Work distribution
The work distributions   () for two different inverse temperatures  = 0.1 and  = 1.0 are shown in Fig. 10, for the protocol from Fig. 2. We obtain an approximately Gaussian distribution for both temperatures.
At this point it is interesting to compare the work distribution with the free energy difference.The free energy difference is marked in Fig. 10 by a dashed black line, the average work by a dashed magenta line.Since Δ < ⟨⟩   , the second law of thermodynamics holds, as expected.However, it is also visible that with a finite probability the extracted work ⟨⟩   for single realizations is smaller than the free energy difference; such datapoints are known as "microscopic violations" of the second law [75].
Furthermore, in order to demonstrate that our circuit operates away from the adiabatic regime, we compute ⟨ diss ⟩ and compare it to Eq. ( 1), see Table IV.Since  diss ≫ 0, it follows that the adiabatic approximation does not hold for our chosen protocols.

Different circuit realizations
Let us also check the effect of different single-qubit random unitaries on our results.To do so, we repeat our experiment for 4 qubits and 3 layers on ibm_perth for different circuit realizations, by choosing different random unitary gates.As we can see in Fig. 11, the deviation from Jarzynski's equality does not depend strongly on the particular choice of random gates in our circuits.2 with seven blocks and  = 8 qubits: theory prediction (red) vs. experimental simulations (blue).The free energy difference Δ is indicated by a dashed black line, the average work by a magenta line.Although the measured work can be smaller than the free energy for some shots, the average work is larger than the free energy difference and the second law of thermodynamics thus holds.The work distributions have a nearly Gaussian shape.The absence of strong tails in the work distribution reduces the number of required measurements.The accuracy of the result is in all cases comparable, i.e., it does not depend strongly on the choice of the random single-qubit gates in the protocol of Fig. 2.

Dependence on the number of circuit blocks
Finally, we investigate the deviations from Eq. (1) as a function of the non-equilibrium-protocol length.To do this, we consider the blocks introduced in part (C) of the circuit in Fig. 2; we can then stack multiple such blocks with different unitaries one after the other to increase the circuit depth.As is shown in Fig. 12, the deviations barely scale with the size of the non-equilibrium protocol for more than two-three blocks.In this case, most of the errors accumulated during the circuit are compensated for by applications of single-qubit random gates, as is explained in Sec.IV.  1) as a function of the circuit depth of the non-equilibrium protocol.The legend shows different inverse temperatures  for a system of  = 8 qubits.A single block here is defined in part (C) of the circuit introduced in Fig. 2. The deviations from Eq. ( 1) are almost independent of the number of blocks in the non-equilibrium protocol.
Here   denotes the probability to prepare a given initial eigenstate, and  → is the transition matrix of the process, defined in Sec. A. In the optimal case of infinitely many measurements and error-free preparation, we have   = e − i  / i .There are two generic ways for errors to occur.The first one is as a statistical error in   , The first term on the right-hand side in the above equation is the probability distribution of the canonical ensemble.The second term denotes the statistical error for the measurement probability of each state | i ⟩, which satisfies the sum rule    = 0.For a large enough number of states, we can assume the errors   to be independent of one another and neglect this constraint.The statistical error of measurement probability for each state | i ⟩ can be modeled by a binomial distribution: For each shot, we obtain this state with probability e − i  / i and measure another state with probability 1 − e − i  / i ; then the statistical uncertainty is given by the variance of the binomial distribution, and scales as where  denotes the number of shots.This expression is the variance of a binomial distribution for an event occurring with probability e − i  / i .The second way an error can occur is through a measurement of  → .We can write where K→ comprises any unitary or doubly-stochastic contributions that Jarzynski's equality is insensitive to.Using the (simplifying) assumption that the correct size of each matrix element is roughly 1/ with  = 2  (this is justified by using a scrambling circuit, and becomes more accurate with increasing system size  ≫ 1), the error of  → scales with the number of measurement shots as This is again the standard deviation for a binomial process of an event with probability 1/ and an effective number of repetitions e − i  / i .We used here (1 − 1/) ≈ 1. Taking these considerations into account, we can now divide the contributions of the fluctuations into four parts: Δ = 1 (E7) Equation (E7) is Jarzynski's equality.Assuming that all terms  → ∼ 1/ (as a consequence of using random circuits, which should not prefer any particular eigenstate transitions), and invoking the central limit theorem (the variance of a sum is the sum of the variances), we find that the second term, Eq. (E8), scales as

(E13)
Let us now test these expressions and the underlying assumptions we made above.Figure 13 shows the scaling of the errors for a noise-free quantum simulation (red) and a comparison with the exact error bars using the exact size of the matrix elements (blue).The data show that the simplification we used to compute error bars is well justified.

FIG. 5 .
FIG. 5. Test of the Crooks relation ln   ( )   (− ) as a function of the work .The black line denotes the theoretical prediction for a noiseless device, Eq. (18).Although the results from the quantum simulation follow the theory prediction, deviations indicate the violation of double-stochasticity.The data shown here are for  = 1.0 , 2 24 ∼ 1.7 × 10 7 shots and  = 8 qubits on ibmq_guadalupe.

FIG. 9 .
FIG.9.Classically computed operator entanglement entropy  1 (), Eq. (D3) of the circuits as a function of the cut-position x, for a different number of blocks of the non-equilibrium protocol (legend), for  = 8 qubits.We chose  − 1 blocks for our simulations, see Fig. 2. The dashed horizontal lines indicate the Page curve  Page = 2 ln 2 − 2 2−−1 , with  = min(,  − ). − 1 blocks are sufficient to generate an entanglement curve close to the Page curve.

FIG. 10 .
FIG.10.Histograms of the work distribution   () for the protocol in Fig.2with seven blocks and  = 8 qubits: theory prediction (red) vs. experimental simulations (blue).The free energy difference Δ is indicated by a dashed black line, the average work by a magenta line.Although the measured work can be smaller than the free energy for some shots, the average work is larger than the free energy difference and the second law of thermodynamics thus holds.The work distributions have a nearly Gaussian shape.The absence of strong tails in the work distribution reduces the number of required measurements.

FIG. 11 .
FIG. 11.Validation of the Jarzynski equality  = 4 qubits as a function of inverse temperature  for different choices of random unitaries in the nonequilibrium protocol, simulated on ibm_perth.The accuracy of the result is in all cases comparable, i.e., it does not depend strongly on the choice of the random single-qubit gates in the protocol of Fig.2.

7 FIG. 12 .
FIG.12.Test of the Jarzynski equality Eq. (1) as a function of the circuit depth of the non-equilibrium protocol.The legend shows different inverse temperatures  for a system of  = 8 qubits.A single block here is defined in part (C) of the circuit introduced in Fig.2.The deviations from Eq. (1) are almost independent of the number of blocks in the non-equilibrium protocol.

FIG. 13 .
FIG.13.Classically emulated validation of the Jarzynski equality Eq. (1) for  = 8 qubits as a function of the number of shots, noisefree results at inverse temperature  = 0.1.The green error bars indicate the error due to a finite number of measurements in the emulation.The orange error bars are computed using the simplification of roughly equal-size transition matrix elements [see text].Blue: Classical emulation, but now with thermal noise of similar size as on ibmq_guadalupe.The impact of shot noise becomes negligible with an increasing number of measurements (error bars decrease), such that the effects of violating double-stochasticity become visible (remaining finite plateau value).

TABLE I .
Comparison of different experiments for the validation of the quantum Jarzynski equality, using the normalized deviation defined in Eq. (17).The relative accuracy E of our simulations is comparable with previous experiments on few-qubit systems.More details for the extraction of the data and the determination of  are provided in the appendix G.
operator  can now be decomposed as