Employing Circuit QED to Measure Nonequilibrium Work Fluctuations

We study an interferometric method for the measurement of the statistics of work performed on a driven quantum system, which has been put forward recently [Dorner et al., Phys. Rev. Lett. 110 230601 (2013), Mazzola et al., Phys. Rev. Lett. 110 230602 (2013)]. The method allows replacing two projective measurements of the energy of the driven system with qubit tomography of an ancilla that is appropriately coupled to it. We highlight that this method could be employed to obtain the work statistics of closed as well as open driven system, even in the strongly dissipative regime. We then illustrate an implementation of the method in a circuit QED set-up, which allows one to experimentally obtain the work statistics of a parametrically driven harmonic oscillator. Our implementation is an extension of the original method, in which two ancilla-qubits are employed and the work statistics is retrieved through two-qubit state tomography. Our simulations demonstrate the experimental feasibility.


Introduction
In the last two decades, the field of nonequilibrium statistical mechanics and thermodynamics has received a great momentum in its development due to the discovery of exact results, known by now as fluctuation relations.They characterize nonequilibrium phenomena in small systems well beyond the regime of linear response (in fact, to any order in the perturbative expansion) and pose stringent conditions on the form that the statistics of non-equilibrium fluctuating quantities, such as work and heat, can assume [1,2,3,4,5,6].For example, the statistics p[w, λ] of work w, performed by varying an external parameter in a time span [0, τ ] according to some pre specified protocol λ, is related to the statistics p[w, λ] performed when applying the time-reversed protocol λ, by the formula (Tasaki-Crooks relation) Here the initial state of the forward (backward) process is a thermal state of temperature 1/k B β and parameter value λ 0 ( λ 0 = λ τ ) where k B is Boltzmann's constant, β is the thermal energy and ∆F is the difference of free energy between the initial state of the backward protocol and the initial state of the forward protocol.We follow the notation of [3], where λ, without subscript, denotes the function specifying the value λ t of the parameter at each time t ∈ [0, τ ], and square brackets refer to the functional dependence.Similar expressions called exchange fluctuation relations [3,7,8] hold in transport scenarios, where one looks at the statistics of energy and/or particles transferred between two reservoirs at different temperature and/or chemical potential.Classically, the Tasaki-Crooks relation (1) has been tested in single molecule stretching experiments, where they have been used to obtain the free energy landscape from nonequilibrium work measurements [9,10,11].In contrast, the experimental verification in the quantum regime is very challenging.The problem lies in the fact that work is not an ordinary quantum mechanical observable [4,12].It cannot be obtained by a single projective measurement but rather by two projective measurements of the initial Hamiltonian H(λ 0 ) at time t = 0, and of the final Hamiltonian H(λ τ ) at time t = τ .The work is then given by the difference of the measured eigenenergies w = E λτ m − E λ 0 n .‡ Huber et al. [14] have proposed an experiment with trapped ions based on this two-measurement scheme but it has not been realized so far.It is worth emphasizing that such experiments would be very important especially because they will provide technological solutions to experimentally access the work statistics p[w, λ], which is a basic building block for the study of thermodynamics in the quantum regime.It is central for the investigation of, e.g., the thermodynamic cost of quantum operations, such as quantum gates, which form the basis of quantum computation and quantum information processing [15].‡ Recently, new quantum fluctuation relations have been found that do not involve projective measurements but focus instead on the change of the quantum expectation of the Hamiltonian [13].This type of relation is not investigated here.
One possible strategy to overcome the difficulties that the two measurements scheme poses has been proposed in Refs.[16,17].There the authors noted that intermediate quantum measurements of arbitrary observables do not alter the validity of the fluctuation relations, thus one might be able to retrieve the wanted information from continuously monitoring some properly chosen quantum observable representing the flux of the wanted quantity.As shown in Ref. [16] this is actually what one does in experiments of bi-directional counting statistics [18], where indeed the exchange fluctuation relation for electron transport has been verified experimentally by looking at the number of electrons crossing an interface [19,20].
Recently, Dorner et al. [21] and Mazzola et al. [22] have put forward a promising method for the measurement of work statistics that avoids the projective energy measurements, and replaces them with state tomography of a qubit (the "ancilla") that is appropriately coupled to the driven system.This possibility was anticipated by Silva, who first pointed out the formal equivalence between the work characteristic function (the Fourier transform of the work statistics) and the Loschmidt echo [23].The proposed implementations use trapped ions [21], and micro or nano-beams coupled to a qubit [22] while an experiment has just been performed using a nuclear magnetic resonance system [24].
With this contribution, we (i) review the method of Dorner et al. [21] and Mazzola et al. [22] (Sec.2), (ii) discuss important extensions thereof (Sec.3), and (iii) illustrate an implementation using a circuit QED setup (Sec.4).Most notably, as we shall discuss in Sec. 3, this new method offers a very promising tool for accessing the work statistics of systems that are strongly coupled to their environment [25].This is very important because so far no experimentally feasible method was known for this case.
The expression "circuit QED" refers to solid state devices that realize on a solidstate micro-chip [26,27] the physics of an atom interacting with a light mode in a cavity, a classic problem of quantum optics [28].Here the role of the atom is played by a superconducting qubit, and the cavity is formed by a planar wave-guide.Such devices have undergone a tremendous and fast development in the last decade [29,30,31], allowing for the experimental study of light-matter physics in parameter regimes that standard quantum-optics experiments cannot reach [32,33] and with an unprecedented flexibility.For example, in a circuit QED device, one can easily manipulate the level spacing of the qubit, which can span a whole range of values, from being resonant with the oscillator, to being far detuned from it.This allows for manipulation of the oscillator state and its read out.For example, following a theoretical proposal [34], Refs.[35,36] report on the qubit-assisted creation and read-out of Fock states and superpositions thereof.Importantly enough for this work, experiments have demonstrated full two qubit tomography [37].Moreover, circuit QED may be used for studying thermodynamic effects on the quantum scale [38,39].
Given the high flexibility offered by the state of the art in circuit QED, we believe that it constitutes a very promising tool-box not only for the development of quantum manipulation and read-out, but also for the study of its thermodynamic cost.The latter is an aspect which has been seldom addressed so far, but is important in order to achieve quantum computers that are not only efficient with respect to the accuracy of the logical quantum gates, but also with respect to avoiding detrimental heating.Here we suggest a circuit QED implementation for the measurement of work statistics, which constitutes a first step towards the development of this tool-box.

The method
In this section we shall briefly review the method for extracting the work statistics put forward by Dorner et al. [21] and Mazzola et al. [22].We shall follow primarily the presentation given in Ref. [21].
Given a driven quantum system described by the Hamiltonian the aim of the method is to provide an experimentally feasible prescription of how one can obtain its work statistics in an experiment.Here λ t and Q denote an externally applied generalized force and its conjugate displacement, respectively.Q is a quantum mechanical observable, whereas λ t is a classical quantity, whose evolution in time from t = 0 to t = τ is pre-specified [40].Prototypical examples are a forced oscillator [41,42], and a parametrically driven oscillator, i.e., an oscillator with a time dependent frequency [43,44].
The traditional prescription requires that the system is prepared at time t = 0 in the thermal state: where β is the inverse thermal energy and Z S (λ 0 ) = Tr e −βH(λ 0 ) .Projective measurements of H(λ 0 ) and H(λ τ ) are then performed at times t = 0 and t = τ , providing us with one eigenvalue of the initial Hamiltonian and one of the final Hamiltonian, E λ 0 n and E λτ m , respectively.The work w = E λτ m − E λ 0 n is then recorded, so that repeated measurements allow one to sample the work probability distribution function where p m|n [λ] denotes the transition probability from state n to state m induced by the protocol λ.
The major obstacle for implementing this prescription comes from the experimental difficulty to perform projective measurements on the system of interest.The method of Dorner et al. [21] and Mazzola et al. [22] circumvents this difficulty by coupling the system to an "ancilla", namely a qubit, which is used to read out the Fourier transform Here U S [λ] is the time evolution operator generated by the driving protocol λ, and • S denotes average over ρ S .Note that since P [w, λ] is a real function, the relation Following Dorner et al. [21], the system (S) is coupled to the ancilla (A) according to the Hamiltonian where σ z = Π + −Π − is the z-Pauli matrix, Π ± = |± ±| is the projector onto the ancilla states |± , and χ ± t are two independent driving protocols of duration T = τ + u which will be specified later.
Because the system-ancilla coupling commutes with the free ancilla Hamiltonian Choosing the time evolution operator reads It contains the operators U S [λ]e −iuH(λ 0 )/ and e −iuH(λτ )/ U S [λ] which appear in the expression of the characteristic function, Eq. ( 5).
Inspired by Ramsey interferometry, the idea is to prepare the system in a superposition of up and down states so that the two time evolutions interfere and the wanted information will be encoded in the state of the ancilla at the final time T = τ +u.This is achieved by the following protocol [21,22]: 1. Prepare the compound system in the state ρ 4. Perform a Hadamard operation σ H on the ancilla at t = T = τ + u.
After this protocol, the ancilla is described by the reduced density matrix where ℜ and ℑ denote real and imaginary parts, Tr S is the trace over the system Hilbert space, and The purpose of the first Hadamard transformation is to create a superposition of the up and down states.The second Hadamard recombines the entries of the ancilla density matrix at time T = τ + u and, hence it is not strictly necessary.
It is worth emphasizing that the diagonal coupling in Eq. ( 6) can often be realized only approximately.The implementation proposed in Sec. 4 is one example of an approximate diagonal coupling obtained by far detuning the ancilla with respect to the system's transition energies.With the perfectly diagonal coupling, the wavefunction is prepared in a superposition in which one component follows the forward protocol χ + , while the other follows the backward protocol χ − .While both components evolve independent of each other with the respective U S [χ ± ] of the isolated system, the ancilla collects the resulting phase difference.

Work statistics of arbitrary open quantum systems
As mentioned above the primary advantage of the interferometric scheme of Dorner et al. [21] and Mazzola et al. [22] is that it avoids projective measurements on the system of interest H S by replacing them with state tomography on the ancilla.This has a very important consequence in regard to the possibility of experimentally testing fluctuation theorems in open quantum systems with arbitrarily strong coupling to a thermal environment [25] Here H B is the thermal bath Hamiltonian and H SB is an arbitrarily strong coupling.Nevertheless, the fluctuation theorem continues to hold unaltered in this case [25], because when driving the system S, part of the injected energy may flow to the bath B and in the SB interaction.Thus the work spent to drive the system is given by the change in the S + B energy: w = E λτ S+B,m − E λ 0 S+B,n .But since one can see the S + B system as a closed system staying initially in a thermal state with (common) inverse temperature β, the ordinary fluctuation relation applies P [w, λ]/P [−w, λ] = e β(w−∆F S+B ) independent of the coupling strength.Using the expression of the free energy of an arbitrary open quantum system [45] Thus the fluctuation theorem remains unaltered in the case of arbitrary open quantum systems [25]: This result is the quantum version of a result obtained by Jarzynski for classical systems [46].The main difference between the classical and the quantum case is that while in the classical case one may obtain the work w performed on the S + B system by looking at the trajectory of S alone [46], in the quantum case this is impossible [3].In the quantum case, in principle, one should perform two projective measurements of the total Hamiltonian H S+B .Making a projective measurement on S alone is already a challenging task in many experimental setups; making a projective measurement of S + B seems much more difficult, if not impossible.The interferometric scheme may be effective in overcoming this problem.If now the open system is coupled to the ancilla which, in turn, has no direct contact to the environment, the S + B + A Hamiltonian reads Implementing the same interferometric scheme as in Sec. 2 with the initial state Π − ρ S+B results in the characteristic function of the open system H S+B (λ t ).Thus the interferometric approach provides, if the ancilla is well isolated, access to the work distribution of arbitrary open as well as closed nonequilibrium quantum systems.Most remarkably, our present discussion highlights that in the interferometric scheme of Dorner et al. [21] and Mazzola et al. [22], deviations from the fluctuation theorem are expected only as a consequence of thermal noise on the ancilla A. Thermal noise on the system S may affect the statistics of work itself, but not the validity of the fluctuation relations.
We emphasize that the fluctuation theorem for open quantum systems described in this section is fully general and exact.In particular it does not require the interaction H SB to be weak nor the initial S + B state to be uncorrelated.Quite on the contrary, in case of strong coupling the initial state ρ S+B contains correlations, and the subsequent evolution of the reduced system density matrix needs not be described by completely positive maps [47], nor has to be Markovian.In this regard newly introduced definitions of characteristic functions for open quantum systems in terms of the reduced system dynamics [48,49,50,51] must be regarded as approximate expressions whose validity is not guaranteed, and whose main object is generally not the work (i.e., the change in energy of S + B) but some other quantity that pertains to the system S only.

Exclusive vs. inclusive work statistics
Fluctuation relations appear in the literature in two complementary ways, referred to as exclusive and inclusive viewpoint [52,53].In our discussion so far we have been adopting the inclusive viewpoint, which addresses the probability of the change w in the system energy H 0 − λ t Q, i.e., including the driving term λ t Q.One may want to look also at the statistics p 0 [w 0 , λ] of the change w 0 in the energy of the system H 0 excluding the driving energy λ t Q §.This is given by w 0 = e m − e n , where e m and e n are eigenvalues of H 0 obtained by projective measurements of H 0 at t = 0 and at t = τ , respectively [53].The exclusive fluctuation relation reads [1,53] while the characteristic function of exclusive work is given by [53] G It differs from the exclusive work characteristic function in Eq. ( 5) in that e −iuH 0 / appears instead of e −iuH(λτ )/ .Accordingly, G 0 [u, λ] can be accessed by means of the interferometric scheme described in Sec. 2 by replacing in Eq. ( 8) the drivings χ ± by so that the evolution takes on the form where we recognize the operators U S [λ]e −iuH 0 / and e −iuH 0 / U S [λ], appearing in the expression of the exclusive characteristic function, Eq. ( 18).Accordingly the state of the ancilla at time τ + u encodes the information on the exclusive characteristic function G 0 [u, λ].

Circuit QED implementation
We want to experimentally access the work statistics of a parametrically driven quantum oscillator, whose frequency changes in time according to ω 2 (t) = ω 2 − 4ωλ t .Its Hamiltonian reads where a = x mω/2 + ip 1/2mω and a † = x mω/2 − ip 1/2mω are the usual bosonic shift operators.To implement this Hamiltonian we consider a circuit QED setup, where a qubit is coupled to a single mode ω of a line resonator [27].The qubit-oscillator system is described by the Rabi Hamiltonian § For the sake of clarity we stress that the dynamics U [λ] is one and the same in both viewpoints.What changes is only what one looks at, namely w in the inclusive case and w 0 in the exclusive case.
where σ x , σ z are Pauli matrices, ε is the qubit energy splitting, and g is the qubitoscillator interaction strength.Note that this Hamiltonian is generally not of the type of that in Eq. ( 6).First, the qubit-system interaction does not commute with the free qubit Hamiltonian εσ z /2.Second, the interaction term is linear in a † + a, whereas we aim at implementing an interaction quadratic in a † + a. Third, typical circuit QED setups do not provide the possibility to control the interaction g in time, because g is fixed by the geometry of the device.With current technology [27] one can relatively easily control the qubit splitting ε, while setups allowing for the control of g so far have been studied only theoretically [54,55].A full description of the transmission line would contain also higher harmonics at multiples of the fundamental frequency ω.The protocol considered below, however, acts directly only upon the fundamental mode, while its harmonics are affected only indirectly via the qubit.We therefore assume that they remain close to their ground state and, thus, do not take them into account.
These issues can be partially solved by considering a time dependent qubit splitting ε t and working in a regime where the coupling g and the oscillator frequency ω are small: By applying the time-dependent unitary transformation and neglecting terms beyond second order in the small parameter g/ε t , we obtain, up to a global energy shift the Hamiltonian The last term comes from the explicit time dependence of the transformation Ω t .We shall consider a qubit driving that is slow compared to the qubit's own time scale while being comparable to the oscillator's time scale In this way the oscillator can be driven out of equilibrium while the qubit undergoes an adiabatic evolution.Note that the factors gω/ε t and g εt /ε 2 t are comparable to the factor g 2 /ε t appearing in the third term of Eq. ( 25).However the last two terms are oscillating much faster and can therefore be neglected within a rotating-wave approximation.This can be seen by going to the interaction picture with respect to ε t σ z /2 + (ωa † a + 1/2), where the last two terms contain the frequencies ±(ε t ± ω) ≃ ±ε t = ±t −1 t 0 ε s ds and the second term contains the much lower frequencies 0, ±2ω.We thus conclude that  λ0) .We used the driving ε t = g 2 /2λ t , where λ t = λ 0 + vt, and the following parameters: is a good approximation of H S+A in the chosen parameter regime.The form ( 27) is already rather close to the desired Hamiltonian in Eq. ( 6).The main difference is that in Eq. ( 6), one drives the two subspaces spanned by Π ± with two independent drivings χ ± t , whereas here we have only one driving parameter ε t that drives both subspaces at the same time.The other difference is that now the free qubit Hamiltonian is time dependent.This affects only an overall phase, which therefore is not our major concern here.Figure 1 illustrates how the original Rabi Hamiltonian ( 22) is well approximated by the diagonal Hamiltonian in Eq. ( 27).In the appendix we provide an alternative derivation of H ′ S+A based on the explicit calculation of the time-evolution generated by H S+A .
Note that the transformation ( 24) is similar but not quite the same as the transformation commonly employed in the dispersive regime [56].We might call the regime investigated here, where the oscillator is very slow, the soft mode regime, and the resulting effective Hamiltonian H ′ S+A , Eq. ( 27), the soft mode Hamiltonian.Like the dispersive Hamiltonian, the soft mode Hamiltonian is diagonal in the natural qubitoscillator basis.But while the dispersive Hamiltonian represents a qubit-oscillator coupling linear in (a + a † ), the soft mode Hamiltonian describes a coupling quadratic in (a + a † ).

Introducing a second qubit
In order to allow for the independent driving of two subspaces, we modify the method described above by introducing a second qubit.The work characteristic function measurement is thus assisted by two ancillae.Two-qubit state tomography has been reported recently in [37].Our starting Hamiltonian is: 2 ).( 28) Following the derivation illustrated above, we shall work in the regime By applying the transformation and neglecting cubic or higher terms in g 1 /ε 1 , g 2 /ε 2 , as well as fast oscillating contributions we arrive at Figure 2 illustrates how the original Rabi Hamiltonian ( 28) is well approximated by the diagonal Hamiltonian in Eq. (31).It is worthwhile rewriting H ′ S+2A in terms of projectors Π ±,± onto the four states |±, ± : where By focussing onto the subspace spanned by Π −+ and Π −− we see that by manipulating the two splittings ε 1,t and ε 2,t , one can realize two independent drivings χ + t and χ − t acting in the respective sub-subspace.This realizes all the ingredients that we need for implementing the characteristic function measurements protocol employing a circuit QED setup.

The protocol
First, the two drivings ε 1,t and ε 2,t are chosen in such a way as to realize the protocols χ + t , χ − t in Eq. ( 8).This is achieved by solving Eq. ( 34) for ε 1,t , ε 2,t to obtain: With this choice, the protocol goes as follows (see Fig. 3 top panel): 1. Prepare the system at t < 0 in the state:  31), red line.The plot shows the evolution of the population of the first three eigenstates of the oscillator.The inset shows the corresponding evolution of the first qubit population.The initial state was ρ S+2A (see Eq. 36), and ε 1,t , ε 2,t were chosen as in Fig. 3, bottom right panel, as to realize the drivings χ ± t shown in Fig. 3, bottom left panel, corresponding to a linear ramp λ t = λ 0 + vt.We used the following parameters: g 1 = 2.5 ω, g 2 = 0.5 ω, λ 0 = 0.0625 ω, v = 1.5λ 0 ω/2π, 1/β = ω.

Perform a Hadamard operation σ H
2 = (σ x 2 + σ z 2 )/ √ 2 on the second qubit at time t = 0. 3. Let the S + 2A system evolve for a time τ + u according to H S+2A (ε 1,t , ε 2,t ). 4. Perform a Hadamard operation σ H 2 at time t = T = τ + u.This results in the two-qubit density matrix where Σ Σ Thus performing two-qubit state tomography gives the characteristic function G[u, λ] of the process in Eq. ( 21) at the point u apart from a known phase factor.The state ρ S+2A can be prepared by thermalizing the S + 2A system at a temperature such that β −1 ≃ ω ≪ ε 1,0 , ε 2,0 .Two qubit-state tomography [57] can be realized in this setup by means of quantum non-demolition joint dispersive read-out [37].This is possible due to the fact that system and oscillator are far detuned.Noticing that only terms involving Σ z 2 and Σ y 2 appear in Eq. ( 37), the wanted information can be retrieved in the following way.(i) Follow the protocol describe above.(ii) At the end of the protocol, perform a measurement

Numerical solution
We numerically studied the case of a linear ramp in the protocol λ t = λ 0 + vt using 1/ω as unit of time and the parameters g 1 = 2.5ω, g 2 = 0.5ω, λ 0 = 0.0625ω, v = 1.5λ 0 ω/2π, τ = 2π/ω.For ω = 100 MHz this amounts to couplings g 1 = 250 MHz and g 2 = 50 MHz, an initial qubit splitting ε 1,0 = 10 GHz and a velocity of v ≈ 150 (MHz) 2 .The level splitting of the second qubit goes to infinity at the beginning and at the end of the protocol, corresponding to a complete decoupling.The cutoff we introduced to handle this divergence is equivalent to ε 2,0 = 10 GHz, which is within current experimental reach.Fabrication of a slow oscillator with a slow frequency of 100 MHz does not seem to pose any special technological challenge.A slow mode oscillator can be made by increasing the resonator's length [58].Perhaps more challenging is reaching the strong coupling g 1 = 2.5ω.There is currently a strong interest in this regime of ultra-strong coupling, and we are optimistic that it will be soon reached [54,59,60].The time development of the two drivings χ ± t , is illustrated in Fig. 3 bottom left.The graph in the bottom right panel of Fig. 3, shows the corresponding time evolution of the two qubit energy splittings ε i,t , i = 1, 2. Note that ε 2,t diverges for t → 0 and for t → τ + u.In our simulation, ε 2,t was cut at the value of 100ω.This results in a deviation of the actual drivings χ ± t from those reported in Fig. 3 bottom left panel, for those values of t where the two χ's approach.For small u (as compared to τ ), this deviation becomes more relevant.With the so chosen parameters, the condition ( 29) was obeyed at all times t ∈ [0, τ + u].
We computed ρ 2A (u) according to Eq. ( 37) where the time evolution was obtained by numerical integration of the Liouville-von Neumann equation.The thermal energy β −1 was chosen equal to ω.We then extracted the real and imaginary parts of the characteristic function G[u, λ] using Eq.(38). Figure 4 shows the work probability distribution obtained after inverse Fourier transform of the so-obtained G[u, λ].The blue dots show the values of the work probability density function as obtained by integrating the model Hamiltonian (21) directly.The approximations introduced by our implementation result in a spread of the peaks, as compared to the expected ones, and to the emergence of further peaks in the work probability at high w (not shown).Because of normalization these effects lower the height of the relevant peaks.We repeated the same procedure for the time reversed protocol λ t = λ τ − vt.The inset of Fig. 4 shows a good agreement between the logarithm of the ratio p[w; λ]/p[−w, λ] as from our numerics, and the linear behavior expected from Eq. ( 1).The agreement is however not as good as one would expect from Fig. 2 showing very good agreement between the dynamics of the model Hamiltonian and actual Hamiltonian.The source of error is coming from the fast oscillating phase e (i/ ) ε 2,t dt in Eq. (38), which has to be taken away before the inverse Fourier transformation is applied.This may pose an issue at the experimental level as well.

Conclusions
We have extended the interferometric scheme of Dorner et al. [21] and Mazzola et al. [22] for the measurement of work distributions.The method lends itself straightforwardly to the application to open quantum systems, even in the regime of strong dissipation, which represents a crucial advantage beyond the works by Dorner et al. [21] and Mazzola et al. [22].We further showed how it can be modified to address the exclusive work fluctuation theorem of Bochkov and Kuzovlev [1].
Our central contribution is the illustration of a realistic implementation of the method with current circuit QED technology.A new feature of the proposed implementation is the introduction of a second ancilla qubit and the use of two-qubit state tomography.This technique may prove useful in all experimental scenarios where, as in the present case, two independent drivings might not be easily achieved with a single qubit.Our numerical calculations show the experimental feasibility.In the proposed implementation, the driving λ t (a † + a) 2 is achieved indirectly by driving the qubit splitting, and working in the soft mode regime (slow oscillator).As an alternative to the proposed implementation one could control λ t directly.This can be implemented by coupling a flux qubit to a SQUID as illustrated in Ref. [54,55].
Thus, by state tomography of the ancilla density matrix at time t = T = τ + u one can recover the value of the characteristic function G[u, λ] for a given u.By repeating the whole procedure for various values of u ∈ (0, ∞), one obtains G[u, λ] on the positive real axis.Using G[−u, λ] = G * [u, λ] one obtains G[u, λ] on the whole real axis, and then, by inverse Fourier transform the work statistics p[w, λ].In practice one can sample the characteristic function only on a finite domain.This, in turn, limits the accuracy with which the work probability distribution function can be resolved.

Figure 2 .
Figure 2. Comparison between the dynamics generated by the Tavis-Cummings Hamiltonian in Eq. (28), black line, and the dynamics generated by the diagonal Hamiltonian in Eq. (31), red line.The plot shows the evolution of the population of the first three eigenstates of the oscillator.The inset shows the corresponding evolution of the first qubit population.The initial state was ρ S+2A (see Eq. 36), and ε 1,t , ε 2,t were chosen as in Fig.3, bottom right panel, as to realize the drivings χ ± t shown in Fig.3, bottom left panel, corresponding to a linear ramp λ t = λ 0 + vt.We used the following parameters: g 1 = 2.5 ω, g 2 = 0.5 ω, λ 0 = 0.0625 ω, v = 1.5λ 0 ω/2π, 1/β = ω.
Figure 3. Top: Schematics of the two-qubit protocol.Bottom left: time evolution of the two driving parameters χ ± t , Eq. (8), for a linear ramp of λ t .Bottom right: the time evolution of the two qubit splittings ε i,t that realize the drivings χ ± t , see Eq. (35).