Employing circuit QED to measure non-equilibrium work fluctuations

In the last two decades, the field of non-equilibrium 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 non-equilibrium 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–6]. For example, the statistics p[w,λ] of work w, performed by varying an external parameter in a time span [0,τ] according to some prespecified protocol λ, is related to the statistics $p[w,\widetilde \lambda ]$ performed when applying the time-reversed protocol $\widetilde \lambda$, by the formula (Tasaki–Crooks relation)

$$\begin{equation} \frac{p[w,\lambda]}{p[-w,\widetilde \lambda]}= \mathrm{e}^{\beta (w-\Delta F)}\, . \end{equation} \tag{ 1 }$$

Here the initial state of the forward (backward) process is a thermal state of temperature 1/k_Bβ and parameter value λ₀( $\widetilde \lambda _0 = \lambda _\tau$), 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 non-equilibrium work measurements [9–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λ0n.6 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 are 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].

One possible strategy to overcome the difficulties that the two-measurements scheme poses has been proposed in [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 [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 the 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] (section 2), (ii) discuss important extensions thereof (section 3) and (iii) illustrate an implementation using a circuit QED setup (section 4). Most notably, as we shall discuss in section 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 solid-state 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–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], Hofheinz et al [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.

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 [21].

Given a driven quantum system described by the Hamiltonian

$$\begin{equation} H_{\mathrm{S}}(\lambda_t) = H_0 - \lambda_t Q , \end{equation} \tag{ 2 }$$

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 τ 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:

$$\begin{equation} \rho_{\mathrm{S}}= \mathrm{e}^{- \beta H(\lambda_0)}/Z_{\mathrm{S}}(\lambda_0) , \end{equation} \tag{ 3 }$$

where β is the inverse thermal energy and ZS(λ0) = Tr e−βH(λ0). Projective measurements of H(λ0) and H(λτ) are then performed at times t = 0 and τ, providing us with one eigenvalue of the initial Hamiltonian and one of the final Hamiltonian, Eλ0n and Eλτm respectively. The work w = Eλτm − Eλ0n is then recorded, so that repeated measurements allow one to sample the work probability distribution function

$$\begin{equation} p[w;\lambda]= \sum_{m,n} \delta(w-(E_m^{\lambda_\tau}-E_n^{\lambda_0}))p_{m|n}[\lambda]\mathrm{e}^{-\beta E_n^{\lambda_0}}/Z_0 , \end{equation} \tag{ 4 }$$

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 in performing 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 G[u,λ] of the work statistics p[w,λ] [12]:

$$\begin{eqnarray} &&\fl G[u,\lambda] = \int \mathrm{d}w \mathrm{e}^{\mathrm{i}uw}p[w;\lambda] = \langle U_{\mathrm{S}}^\dagger[\lambda] \mathrm{e}^{\mathrm{i}u H(\lambda_\tau)/\hbar} U_{\mathrm{S}}[\lambda] \mathrm{e}^{-\mathrm{i}u H(\lambda_0)/\hbar} \rangle_{\mathrm{S}} \nonumber \\ &&= \langle (\mathrm{e}^{-\mathrm{i}u H(\lambda_\tau)/\hbar} U_{\mathrm{S}}[\lambda])^\dagger U_{\mathrm{S}}[\lambda] \mathrm{e}^{-\mathrm{i}u H(\lambda_0)/\hbar} \rangle_{\mathrm{S}} . \end{eqnarray} \tag{ 5 }$$

Here U_S[λ] is the time evolution operator generated by the driving protocol λ, and 〈·〉_S denotes the average over ρ_S. Note that since P[w,λ] is a real function, the relation G[ − u,λ] = G*[u,λ] holds.

Following Dorner et al [21], the system (S) is coupled to the ancilla (A) according to the Hamiltonian

$$\begin{equation} H_{\mathrm{S+A}} = \frac{\hbar\varepsilon}{2}\sigma^z + H_0 - (\chi^+_t \Pi_+ + \chi^-_t \Pi_-) Q , \end{equation} \tag{ 6 }$$

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 H_A = εσ^z/2, the evolution U_S+A[χ⁺,χ⁻] of S + A generated by the drivings χ⁺,χ⁻ is block diagonal in the basis {| + 〉,| − 〉}:

$$\begin{equation} U_{\mathrm{S+A}}[\chi^+,\chi^-] = \left(\begin{array}{@{}cc@{}}\mathrm{e}^{-\mathrm{i}\varepsilon T/2\hbar }U_{\mathrm{S}}[\chi^+] & 0 \\ 0 & \mathrm{e}^{\mathrm{i}\varepsilon T/2\hbar }U_{\mathrm{S}}[\chi^-]\end{array}\right) . \end{equation} \tag{ 7 }$$

Choosing

$$\begin{equation} \chi^+_t=\left\{ \begin{array}{@{}lll@{}} \lambda_t \quad {\rm for} \ t \in [0,\tau], \\ \lambda_\tau \quad {\rm for} \ t \in [\tau,\tau+u], \end{array}\right. \quad \chi^-_t=\left\{ \begin{array}{@{}lll@{}} \lambda_0 \quad {\rm for} \ t \in [0,u], \\ \lambda_{t-u} \quad {\rm for} \ t \in [u,\tau+u], \end{array}\right. \end{equation} \tag{ 8 }$$

the time evolution operator reads

$$\begin{equation} U_{\mathrm{S+A}}[\chi^+,\chi^-] = \left(\begin{array}{@{}cc@{}}\mathrm{e}^{-\mathrm{i}\varepsilon (\tau+u)/2\hbar } \mathrm{e}^{- \mathrm{i}u H(\lambda_\tau)/\hbar} U_{\mathrm{S}}[\lambda] & 0 \\ 0 & \mathrm{e}^{\mathrm{i}\varepsilon (\tau+u)/2\hbar }U_{\mathrm{S}}[\lambda] \mathrm{e}^{-\mathrm{i}u H(\lambda_0)/\hbar}\end{array}\right) . \end{equation} \tag{ 9 }$$

It contains the operators U_S[λ]e^{−iuH(λ₀)/ℏ} and e^{−iuH(λ_τ)/ℏ}U_S[λ], which appear in the expression of the characteristic function, equation (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 ρ_S+A = | − 〉〈 − |ρ_S at t < 0.
• 2.  
  Perform a Hadamard operation $\sigma ^{\mathrm {H}}=(\sigma ^x+\sigma ^z)/\sqrt {2}$ on the ancilla at t = 0.
• 3.  
  Let the S + A system evolve for a time τ + u according to U_S+A[χ⁺,χ⁻].
• 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

$$\begin{eqnarray} &&\rho_{\mathrm{A}}(u)= {\mathrm{Tr}}_{\mathrm{S}} \sigma^{\mathrm{H}} U_{\mathrm{S+A}} [\chi^+,\chi^-]\sigma^{\mathrm{H}} \rho_{\mathrm{S+A}} \sigma^{\mathrm{H}} U^\dagger_{\mathrm{S+A}}[\chi^+,\chi^-] \sigma^{\mathrm{H}} \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} &&= \left(1- [\Re L(u)] \sigma^z + [\Im L(u)]\sigma^y \right)/2 , \end{eqnarray} \tag{ 11 }$$

where ℜ and ℑ denote real and imaginary parts, Tr_S is the trace over the system Hilbert space, and

$$\begin{eqnarray} &&L(u)= \mathrm{e}^{-\mathrm{i}\varepsilon (\tau + u)/\hbar } G[u,\lambda] . \end{eqnarray} \tag{ 12 }$$

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.

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 equation (6) can often be realized only approximately. The implementation proposed in section 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.

3.1. 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]

$$\begin{equation} H_{\mathrm{S+B}}(\lambda_t) = H_{\mathrm{S}}(\lambda_t) + H_{\mathrm{B}} + H_{\mathrm{SB}} \, . \end{equation} \tag{ 13 }$$

Here HB is the thermal bath Hamiltonian and HSB 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 β,

$$\begin{equation} \rho_{\mathrm{S+B}} = \frac{\mathrm{e}^{-\beta H_{\mathrm{S+B}}(\lambda_0)}}{Z_{\mathrm{S+B}}(\lambda_0)}\, , \end{equation} \tag{ 14 }$$

the ordinary fluctuation relation $P[w,\lambda ]/P[-w,\widetilde {\lambda }]= \mathrm {e}^{\beta (w-\Delta F_{\mathrm {S+B}})}$ applies, independent of the coupling strength. Using the expression of the free energy of an arbitrary open quantum system [45] F_S(λ_t) = F_S+B(λ_t) − F⁰_B (where F⁰_B = −β⁻¹ ln Tr_B e^−βH_B, and F_S+B(λ_t) = −β⁻¹ ln Tr_S+B e^{−βH_S+B(λ_t)}), one immediately sees that ΔF_S+B = F_S+B(λ_t) − F_S+B(λ₀) = F_S(λ_t) − F_S(λ₀) = ΔF_S. Thus the fluctuation theorem remains unaltered in the case of arbitrary open quantum systems [25]:

$$\begin{equation} \frac{P[w,\lambda]}{P[-w,\widetilde{\lambda}]}= \mathrm{e}^{\beta(w-\Delta F_{\mathrm{S}})}\, . \end{equation} \tag{ 15 }$$

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

$$\begin{equation} H_{\mathrm{S+B+A}} = \frac{\hbar\varepsilon}{2}\sigma^z + H_0 + H_{\mathrm{B}} + H_{\mathrm{SB}} - (\chi^+_t \Pi_+ + \chi^-_t \Pi_-) Q . \end{equation} \tag{ 16 }$$

Implementing the same interferometric scheme as in section 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 non-equilibrium 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–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.

3.2. Exclusive versus inclusive work statistics

Fluctuation relations appear in the literature in two complementary ways, referred to as exclusive and inclusive viewpoints [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₀ − λ_tQ, i.e. including the driving term λ_tQ. One may want to look also at the statistics p₀[w₀,λ] of the change w₀ in the energy of the system H₀ excluding the driving energy λ_tQ.⁷ This is given by w₀ = e_m − e_n, where e_m and e_n are eigenvalues of H₀ obtained by projective measurements of H₀ at t = 0 and τ, respectively [53]. The exclusive fluctuation relation reads [1, 53]

$$\begin{equation} \frac{p_0[w_0,\lambda]}{p_0[-w_0,\widetilde{\lambda}]}= \mathrm{e}^{\beta w_0}, \end{equation} \tag{ 17 }$$

while the characteristic function of exclusive work is given by [53]

$$\begin{eqnarray} &&\fl G_0[u,\lambda] = \int \mathrm{d}w \mathrm{e}^{\mathrm{i}uw}p_0[w;\lambda] = \langle U_{\mathrm{S}}^\dagger[\lambda] \mathrm{e}^{\mathrm{i}u H_0/\hbar} U_{\mathrm{S}}[\lambda] \mathrm{e}^{-\mathrm{i}u H_0/\hbar} \rangle_{\mathrm{S}} \nonumber \\ &&= \langle (\mathrm{e}^{-\mathrm{i}u H_0/\hbar} U_{\mathrm{S}}[\lambda])^\dagger U_{\mathrm{S}}[\lambda] \mathrm{e}^{-\mathrm{i}u H_0/\hbar} \rangle_{\mathrm{S}} . \end{eqnarray} \tag{ 18 }$$

It differs from the exclusive work characteristic function in equation (5) in that e^{−iuH₀/ℏ} appears instead of e^{−iuH(λ_τ)/ℏ}. Accordingly, G₀[u,λ] can be accessed by means of the interferometric scheme described in section 2 by replacing in equation (8) the drivings χ^± by

$$\begin{equation} \chi^+_t=\left\{ \begin{array}{@{}lll@{}} \lambda_t \quad {\rm for} \ t \in [0,\tau], \\ 0 \quad {\rm for} \ t \in [\tau,\tau+u], \end{array}\right. \quad \chi^-_t=\left\{ \begin{array}{@{}lll@{}} 0 \quad {\rm for} \ t \in [0,u], \\ \lambda_{t-u} \quad {\rm for} \ t \in [u,\tau+u], \end{array}\right. \end{equation} \tag{ 19 }$$

so that the evolution takes on the form

$$\begin{equation} U_{\mathrm{S+A}}[\chi^+,\chi^-] = \left(\begin{array}{@{}cc@{}}\mathrm{e}^{-\mathrm{i}\varepsilon (\tau+u)/2\hbar } \mathrm{e}^{- \mathrm{i}u H_0/\hbar} U_{\mathrm{S}}[\lambda] & 0 \\ 0 & \mathrm{e}^{\mathrm{i}\varepsilon (\tau+u)/2\hbar }U_{\mathrm{S}}[\lambda] \mathrm{e}^{-\mathrm{i}u H_0/\hbar}\end{array} \right) , \end{equation} \tag{ 20 }$$

where we recognize the operators U_S[λ]e^{−iuH₀/ℏ} and e^{−iuH₀/ℏ}U_S[λ], appearing in the expression of the exclusive characteristic function, equation (18). Accordingly, the state of the ancilla at time τ + u encodes the information on the exclusive characteristic function G₀[u,λ].

We want to experimentally access the work statistics of a parametrically driven quantum oscillator, whose frequency changes in time according to ω²(t) = ω² − 4ωλ_t. Its Hamiltonian reads

$$\begin{eqnarray} H_{\mathrm{S}}(\lambda_t) &=& p^2/2m + m(\omega^2- 4 \omega \lambda_t) x^2/2 \nonumber \\ &=& \hbar\omega (a^\dagger a +1/2)- \hbar \lambda_t (a^\dagger +a)^2 , \end{eqnarray} \tag{ 21 }$$

where $a= x \sqrt {m \omega /2\hbar } + \mathrm {i} p \sqrt {1 /2m \omega \hbar }$ and $a^\dagger = x \sqrt {m \omega /2\hbar } - \mathrm {i} p \sqrt {1 /2m \omega \hbar }$ 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

$$\begin{eqnarray} H_{\mathrm{S+A}} &= \hbar \varepsilon \sigma^z/2 + \hbar \omega(a^\dagger a+1/2)+ \hbar g(a^\dagger +a)\sigma^x , \end{eqnarray} \tag{ 22 }$$

where σ^x, σ^z are Pauli matrices, ε is the qubit energy splitting and g is the qubit–oscillator interaction strength. Note that this Hamiltonian is generally not of the type of that in equation (6). Firstly, the qubit–system interaction does not commute with the free qubit Hamiltonian ℏεσ^z/2. Secondly, the interaction term is linear in a^† + a, whereas we aim at implementing an interaction quadratic in a^† + a. Thirdly, typical circuit QED setups do not provide the possibility of controlling 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 also contain 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:

$$\begin{equation} g \simeq \omega \ll \varepsilon_t \,. \end{equation} \tag{ 23 }$$

By applying the time-dependent unitary transformation

$$\begin{equation} \Omega_t= \mathrm{e}^{\mathrm{i} g(a^\dagger +a)\sigma_y/\varepsilon_t} \end{equation} \tag{ 24 }$$

and neglecting terms beyond second order in the small parameter g/ε_t, we obtain, up to a global energy shift, the Hamiltonian

$$\begin{eqnarray} &&\fl H_{\mathrm{S+A}}'(\varepsilon_t) = \Omega_t H_{\mathrm{S+A}} \Omega^ \dagger_t + \mathrm{i}\dot \Omega_t \Omega_t^\dagger \nonumber\\ &&= \frac{\hbar\varepsilon_t}{2} \sigma^z + \hbar \omega \!\left(a^ \dagger a+\frac{1}{2}\right)\! + \frac{\hbar g^2}{\varepsilon_t} (a^\dagger +a)^2 \sigma^z + \mathrm{i} \frac{\hbar \omega g}{\varepsilon_t}(a^\dagger-a)\sigma_y - \frac{ \hbar g \dot \varepsilon_t}{\varepsilon^2_t}(a^\dagger+a)\sigma_y . \nonumber\\ \end{eqnarray} \tag{ 25 }$$

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

$$\begin{equation} \dot \varepsilon_t / \varepsilon_t \simeq \omega \ll \varepsilon_t \,. \end{equation} \tag{ 26 }$$

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 \dot \varepsilon _t/\varepsilon ^2_t$ are comparable to the factor g²/ε_t appearing in the third term of equation (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 $\pm (\bar \varepsilon _t \pm \omega )\simeq \pm \bar \varepsilon _t=\pm t^{-1} \int _0^t \varepsilon _s\,\mathrm {d}s$ and the second term contains the much lower frequencies 0,±2ω. We thus conclude that

$$\begin{eqnarray} &&\fl H'_{\mathrm{S+A}} (\varepsilon_t) = \frac{\hbar \varepsilon_t}{2} \sigma^z + \hbar \omega \left(a^\dagger a +\frac{1}{2}\right)+ \frac{\hbar g^2}{\varepsilon_t} (a^\dagger +a)^2 \sigma^z \nonumber \\ &&= \frac{\hbar \varepsilon_t}{2} \sigma^z +\hbar \omega \left(a^\dagger a + \frac{1}{2}\right) + \left(\frac{\hbar g^2}{\varepsilon_t}\Pi_+ - \frac{\hbar g^2}{\varepsilon_t} \Pi_- \right) (a^\dagger +a)^2 \end{eqnarray} \tag{ 27 }$$

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 equation (6). The main difference is that in equation (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 equation (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.

Zoom In Zoom Out Reset image size

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, equation (27), the soft mode Hamiltonian. Like the dispersive Hamiltonian, the soft mode Hamiltonian is diagonal in the natural qubit–oscillator 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^†).

4.1. 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

$$\begin{eqnarray} &&\hspace{-1cm}H_{\mathrm{S+2A}} = \hbar \varepsilon_1 \sigma^z_1/2 + \hbar \varepsilon_2 \sigma^z_2/2 + \hbar \omega(a^\dagger a+1/2)+ \hbar(a^\dagger +a)(g_1\sigma^x_1+g_2\sigma^x_2) . \end{eqnarray} \tag{ 28 }$$

Following the derivation illustrated above, we shall work in the regime

$$\begin{equation} \omega, g_i \ll \varepsilon_{i,t}\, , \quad \dot \varepsilon_{i,t} / \varepsilon_{i,t} \simeq \omega \ll \varepsilon_{i,t}\, , \quad i=1,2 . \end{equation} \tag{ 29 }$$

By applying the transformation

$$\begin{equation} \Omega_t = \mathrm{e}^{\mathrm{i}g_1(a^\dagger +a)\sigma^y_1/\varepsilon_{1,t}}\mathrm{e}^{\mathrm{i}g_2(a^\dagger +a)\sigma^y_2/\varepsilon_{2,t}} \end{equation} \tag{ 30 }$$

and neglecting cubic or higher terms in g₁/ε₁,g₂/ε₂, as well as fast oscillating contributions, we arrive at

$$\begin{eqnarray} &&\fl H'_{\mathrm{S+2A}}(\varepsilon_{1,t},\varepsilon_{2,t}) = \frac{\hbar} {2}\varepsilon_{1,t} \sigma^z_1 + \frac{\hbar}{2} \varepsilon_{2,t} \sigma^z_2 + \hbar \omega \left(a^\dagger a + \frac{1}{2}\right)+ \hbar (a^\dagger +a)^2\left( \frac{g_1^2}{\varepsilon_{1,t} } \sigma^z_1+ \frac{g_2^2}{\varepsilon_{2,t}} \sigma^z_2\right). \end{eqnarray} \tag{ 31 }$$

Figure 2 illustrates how the original Rabi Hamiltonian (28) is well approximated by the diagonal Hamiltonian in equation (31). It is worthwhile rewriting H'_S+2A in terms of projectors Π_±,± onto the four states | ± ,±〉:

$$\begin{eqnarray} &&\fl H'_{\mathrm{S+2A}}(\varepsilon_{1,t},\varepsilon_{2,t}) =\hbar [\varepsilon^+_t + \omega (a^\dagger a+1/2) + \chi^+_t (a^\dagger +a)^2]\Pi_{++} \nonumber \\ &&+\, \hbar [\varepsilon^-_t + \omega (a^\dagger a+1/2) + \chi^-_t (a^\dagger +a)^2]\Pi_{+-} \nonumber \\ &&+\,\hbar [-\varepsilon^-_t + \omega (a^\dagger a+1/2) - \chi^-_t (a^\dagger +a)^2]\Pi_{-+} \nonumber \\ &&+\,\hbar [-\varepsilon^+_t + \omega (a^\dagger a+1/2) - \chi^+_t (a^\dagger +a)^2]\Pi_{--} , \end{eqnarray} \tag{ 32 }$$

where

$$\begin{eqnarray} &&\varepsilon^{\pm}_t = \frac{\varepsilon_{1,t} \pm \varepsilon_{2,t}}{2} , \end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray} &&\chi^{\pm}_t = \frac{g_1^2}{\varepsilon_{1,t}} \pm \frac{g_2^2}{\varepsilon_{2,t}}. \end{eqnarray} \tag{ 34 }$$

By focusing 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.

Zoom In Zoom Out Reset image size

4.2. 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 equation (8). This is achieved by solving equation (34) for ε_1,t,ε_2,t to obtain

$$\begin{equation} \varepsilon_{1,t} = \frac{2g_1^2}{ \chi^{+}_t+ \chi^{-}_t} \, , \quad \varepsilon_{2,t} = \frac{2g_2^2}{ \chi^{+}_t- \chi^{-}_t} . \end{equation} \tag{ 35 }$$

With this choice, the protocol goes as follows (see figure 3 top panel):

• 1.  
  Prepare the system at t < 0 in the state:

  $$\begin{equation} \rho_{\mathrm{S+2A}} = \frac{\mathrm{e}^{-\beta (\omega a^\dagger a - \lambda_0(a^\dagger+a)^2)}} {Z_{\mathrm{S+2A}}(\lambda_0)} \Pi_{--} . \end{equation} \tag{ 36 }$$
• 2.  
  Perform a Hadamard operation $\sigma ^{\mathrm {H}}_2= (\sigma ^x_2+\sigma ^z_2)/\sqrt {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₂ at time t = T = τ + u.

Zoom In Zoom Out Reset image size

This results in the two-qubit density matrix

$$\begin{eqnarray} \rho_{2\mathrm{A}}(u) &=& {\mathrm{Tr}}_{\mathrm{S}} \sigma_2^{\mathrm{H}} U_{\mathrm{S+A}} [\chi^+,\chi^-]\sigma_2^{\mathrm{H}} \rho_{\mathrm{S +2A}} \sigma^{\mathrm{H}}_2 U^\dagger_{\mathrm{S+A}}[\chi^+,\chi^-] \sigma^{\mathrm{H}}_2 \nonumber \nonumber \\ &=& \left(1- \Re L_2(u)\, \Sigma^z_2 - \Im L_2(u)\, \Sigma^y_2 \right)/2, \end{eqnarray} \tag{ 37 }$$

where

$$\begin{eqnarray} L_2(u) &=& \mathrm{e}^{(\mathrm{i}/\hbar) \int_0^{\tau+u}\varepsilon_{2,t} \mathrm{d}t} G[u,\lambda], \end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray}\Sigma^z_2 &=& \Pi_{-+}-\Pi_{--}, \end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray} \Sigma^y_2 &=& |{-}{+} \rangle \langle {-}{-}|\, +\, |{-}{-} \rangle \langle {-}{+}| \, . \end{eqnarray} \tag{ 40 }$$

Thus performing two-qubit state tomography gives the characteristic function G[u,λ] of the process in equation (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,0,ℏε_2,0.

Two qubit-state tomography [57] can be realized in this setup by means of a quantum non-demolition joint dispersive read-out [37]. This is possible due to the fact that the system and oscillator are far detuned. Noticing that only terms involving Σ^z₂ and Σ^y₂ appear in equation (37), the wanted information can be retrieved in the following way: (i) follow the protocol described above; (ii) at the end of the protocol, perform a measurement of the two-qubit observables Σ^z₂ and Σ^y₂. Repeat (i) and (ii) many times to obtain the expectation values 〈Σ^z₂〉 and 〈Σ^y₂〉. Then ℜL₂(u) = −〈Σ^z₂〉 ℑL₂(u) = −〈Σ^y₂〉.

4.3. Numerical solution

We numerically studied the case of a linear ramp in the protocol λ_t = λ₀ + vt using 1/ω as unit of time and the parameters g₁ = 2.5ω, g₂ = 0.5ω, λ₀ = 0.0625ω, v = 1.5λ₀ω/(2π), τ = 2π/ω. For ω = 100 MHz, this amounts to couplings g₁ = 250 and g₂ = 50 MHz, an initial qubit splitting ε_1,0 = 10 GHz and a velocity of v ≈ 150 (MHz)². 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₁ = 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 figure 3 (bottom left). The graph in the bottom right panel of figure 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 figure 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 equation (37), where the time evolution was obtained by numerical integration of the Liouville–von Neumann equation. The thermal energy β⁻¹ was chosen equal to ℏω. We then extracted the real and imaginary parts of the characteristic function G[u,λ] using equation (38). Figure 4 shows the work probability distribution obtained after the 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 $\widetilde \lambda _t = \lambda _\tau - v t$. The inset of figure 4 shows a good agreement between the logarithm of the ratio $p[w;\lambda ]/p[-w,\widetilde \lambda ]$ from our numerics, and the linear behavior expected from equation (1). The agreement is however not as good as one would expect from figure 2, showing very good agreement between the dynamics of the model Hamiltonian and actual Hamiltonian. The source of this error comes from the fast oscillating phase $\mathrm {e}^{(\mathrm {i}/\hbar )\int \varepsilon _{2,t} \mathrm {d}t}$ in equation (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.

Zoom In Zoom Out Reset image size

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 in 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)² 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 [54, 55].

This work was supported by the German Excellence Initiative 'Nanosystems Initiative Munich (NIM)' (MC, RB, PH), the Collaborative Research Center SFB 631 (RB, PH), the Spanish Government grants MAT2011-24331 (SK) and FIS2011-25167 (DZ), the DGA project FENOL (DZ), the EU project PROMISCE (DZ) and the COST Action MP1209 (MC, PH).

The following derivation of the soft mode Hamiltonian is along the lines of the derivation of the dispersive Hamiltonian presented by Schleich [28]. Our starting point is the Rabi Hamiltonian, equation (22). For simplicity we do not include the time dependence of ε. In the interaction picture, the qubit–oscillator coupling reads

$$\begin{equation} H^{\mathrm{I}}_{\mathrm{SA}}(t) = \hbar g(\sigma_+ a \mathrm{e}^{\mathrm{i}\Delta t}+ \sigma_- a^\dagger \mathrm{e}^{-\mathrm{i}\Delta t}+ \sigma_- a \mathrm{e}^{-\mathrm{i}\Gamma t}+ \sigma_+ a\mathrm{e}^{\mathrm{i}\Gamma t}) , \end{equation} \tag{ A.1 }$$

where Δ = ε − ω, Γ = ε + ω and σ_± are the qubit rising and lowering operators. In the interaction picture, H^I_SA(t) is the generator of the dynamics:

$$\begin{eqnarray} U_{t,0}&=& T\exp \left( -\frac{\mathrm{i}}{\hbar} \int_0^t \mathrm{d}t' H^{\mathrm{I}}_{\mathrm{SA}}(t')\right)\nonumber \\ &\simeq& 1 -\frac{\mathrm{i}}{\hbar} \int_0^t \mathrm{d}t' H^{\mathrm{I}}_{\mathrm{SA}}(t') - \frac{1}{\hbar^2} \int_0^t \mathrm{d}t' H^{\mathrm{I}}_{\mathrm{SA}}(t') \int_0^{t'} \mathrm{d}t'' H^{\mathrm{I}}_{\mathrm{SA}}(t'') . \end{eqnarray} \tag{ A.2 }$$

Plugging equation (A.1) into (A.2), the first order term reads

$$\begin{eqnarray} &&\int_0^t \mathrm{d}t' H^{\mathrm{I}}_{\mathrm{SA}}(t') = \hbar g\left( \sigma_+ a \frac{\mathrm{e}^{\mathrm{i}\Delta t} -1 }{\mathrm{i} \Delta} + \sigma_- a^\dagger \frac{\mathrm{e}^{-\mathrm{i}\Delta t} -1 }{-\mathrm{i} \Delta} + \sigma_- a \frac{\mathrm{e}^{-\mathrm{i}\Gamma t} -1 }{-\mathrm{i} \Gamma} + \sigma_+ a^\dagger \frac{\mathrm{e}^{\mathrm{i}\Gamma t} -1 }{\mathrm{i} \Gamma} \right) ,\nonumber \\ \end{eqnarray} \tag{ A.3 }$$

which can be used the to calculate the second order term

$$\begin{eqnarray} &&\fl \int_0^t \mathrm{d}t' H^{\mathrm{I}}_{\mathrm{SA}}(t') \int_0^{t'} \mathrm{d}t'' H^{\mathrm{I}}_{\mathrm{SA}}(t'') = \hbar ^2g^2\int_0^t \mathrm{d} t' \nonumber \\ &&\times \left( \sigma_+ \sigma_- a a^\dagger \frac{1-\mathrm{e}^{\mathrm{i}\Delta t'}}{-\mathrm{i} \Delta}+ \sigma_+ \sigma_- a^2 \frac{\mathrm{e}^{\mathrm{i}(\Delta-\Gamma) t'} -\mathrm{e}^{-\mathrm{i}\Delta t'}}{-\mathrm{i} \Gamma}\right. + \sigma_- \sigma_+ a^\dagger a \frac{1-\mathrm{e}^{-\mathrm{i}\Delta t'} }{\mathrm{i} \Delta}\nonumber \\ &&+\, \sigma_- \sigma_+ a^{\dagger 2} \frac{\mathrm{e}^{-\mathrm{i}(\Delta-\Gamma) t'}-\mathrm{e}^{-\mathrm{i}\Delta t'} }{\mathrm{i} \Gamma} + \sigma_-\sigma_+ a^2 \frac{\mathrm{e}^{\mathrm{i}(\Delta-\Gamma) t'} -\mathrm{e}^{-\mathrm{i}\Gamma t'}}{\mathrm{i} \Delta} + \sigma_- \sigma_+ a a^\dagger \frac{1-\mathrm{e}^{-\mathrm{i}\Gamma t'} }{\mathrm{i} \Gamma} \nonumber \\ &&\left .+\, \sigma_+ \sigma_- a^{\dagger 2} \frac{\mathrm{e}^{-\mathrm{i}(\Delta-\Gamma) t'} -\mathrm{e}^{\mathrm{i} \Gamma t'}}{-\mathrm{i} \Delta} + \sigma_+ \sigma_- a^\dagger a\frac{1-\mathrm{e}^{-\mathrm{i}\Gamma t'} }{-\mathrm{i} \Gamma} \right) . \end{eqnarray} \tag{ A.4 }$$

Recalling that the oscillator is slow compared to the qubit, for times which are short compared to the oscillator's period we can employ the approximation e^{±i(Δ−Γ)t'} = e^−±2iωt' ≃ 1. In performing the integration we neglect the fast oscillating terms of frequency Γ,Δ ≃ ε to obtain

$$\begin{eqnarray} &&\frac{1}{\hbar^2}\int_0^t \mathrm{d}t' H^{\mathrm{I}}_{\mathrm{SA}}(t') \int_0^{t'} \mathrm{d}t'' H^{\mathrm{I}}_{\mathrm{SA}}(t'') \simeq \frac{\mathrm{i} g^2}{\varepsilon} (a^\dagger + a)^2 \sigma^z t , \end{eqnarray} \tag{ A.5 }$$

where we used σ₊σ₋ − σ₋σ₊ = σ^z. Note that the first order term, equation (A.3), contains either fast oscillating contributions or non-relevant constant terms. Therefore it can be neglected at once, so that the propagator becomes

$$\begin{eqnarray} U_{t,0} \simeq& 1 -\frac{\mathrm{i} g^2}{\varepsilon} (a^\dagger + a)^2 \sigma^z t \simeq \exp\left(-\mathrm{i} \frac{ g^2}{\varepsilon} (a^\dagger + a)^2 \sigma^z t \right), \end{eqnarray} \tag{ A.6 }$$

which corresponds to the Hamiltonian $H^{\mathrm {I}}_{\mathrm {S+A}}(t) \simeq \frac {\hbar g^2}{\varepsilon } (a^\dagger + a)^2 \sigma ^z$. Going back to the Schrödinger picture, we finally arrive at H_S+A ≃ ℏεσ^z/2 + ℏ(ωa^†a + 1/2) + ℏg²(a^† + a)²σ^z/ε.