Using a quantum work meter to test non-equilibrium fluctuation theorems

TL;DR

This paper introduces an experimental quantum work meter using cold atoms, capable of directly measuring and sampling work distributions to test non-equilibrium fluctuation theorems like the Jarzynski identity.

Contribution

The paper presents the first experimental implementation of a quantum work meter that directly measures work distributions in quantum systems.

Findings

Successful measurement of work distribution in a quantum system.

Verification of the Jarzynski fluctuation theorem in a quantum context.

Demonstration of a new tool for quantum thermodynamics experiments.

Abstract

Non-equilibrium quantum thermodynamics is essential to describe new devices that operate far from the regime where the usual thermodynamical laws are obeyed. When quantum fluctuations dominate, defining and measuring work and heat, two central concepts in classical thermodynamics, is non-trivial. For driven, but otherwise isolated, quantum systems, work is a random variable associated with the change in the internal energy, as the first law of thermodynamics indicates. In this paper we present the design and the experimental implementation of a quantum work meter (QWM) operating on an ensemble of cold atoms. Our QWM not only directly measures work but also directly samples its probability distribution. As the work probability distribution plays a central role in the fluctuation theorems of non-equilibrium quantum thermodynamics, the QWM is an ideal tool to test their validity. In particular, we use it to verify the Jarzynski identity.

Figures

Tables

Table 1. Table 1: Estimates of β ​ Δ ​ F 𝛽 Δ 𝐹 \beta\Delta F and Δ ​ F Δ 𝐹 \Delta F for three different temperatures. We show the estimation obtained using the Jarzynski identity (JI) and from a direct calculation of the partition function (PF).

$\beta E$	$\beta \Delta F$ (JI)	$\beta \Delta F$ (PF)	$\Delta F/E$ (JI)	$\Delta F/E$ (PF)
$0.58\pm 0.02$	$-0.36\pm 0.04$	$-0.35\pm 0.03$	$-0.62\pm 0.07$	$-0.60\pm 0.06$
$1.11\pm 0.02$	$-0.63\pm 0.05$	$-0.63\pm 0.04$	$-0.57\pm 0.05$	$-0.57\pm 0.04$
$1.75\pm 0.04$	$-0.92\pm 0.09$	$-0.93\pm 0.06$	$-0.53\pm 0.05$	$-0.53\pm 0.04$

Equations

$$⟨e^{-\beta w}⟩=e^{-\beta \Delta F},$$

$$P(w)=\sum _{n,m}{p}_n{p}_{m|n}\delta [w-({\tilde{E}}_m-E_n)].$$

$$P(w)=\sum _{n,m}{p}_n{p}_{m|n}\delta \left[w-({\tilde{E}}_m-E_n)\right],$$

$$p_n=\mathrm{Tr}(\rho {\Pi }_n),p_{m|n}=\frac{1}{p_n}\mathrm{Tr}({\tilde{\Pi }}_m{\mathcal{U}}_{\mathcal{S}}{\Pi }_n\rho {\Pi }_n{\mathcal{U}}_{\mathcal{S}}^{†}).$$

$$⟨e^{-\beta w}⟩$$

$$⟨e^{-\beta w}⟩=\frac{1}{Z}\sum _m\mathrm{Tr}({\tilde{\Pi }}_m)e^{-\beta {\tilde{E}}_m}=\frac{\tilde{Z}}{Z},$$

$$P(w_{nm})=p_{m|n}{p}_n.$$

$$P(w_{nm})=\mathrm{Tr}({\tilde{\Pi }}_m{\mathcal{U}}_{\mathcal{S}}{\Pi }_n\rho {\Pi }_n{\mathcal{U}}_{\mathcal{S}}^{†})=\mathrm{Tr}(\rho A_{nm}),$$

$$A_{nm}={\Pi }_n{\mathcal{U}}_{\mathcal{S}}^{†}{\tilde{\Pi }}_m{\mathcal{U}}_{\mathcal{S}}{\Pi }_n.$$

$$|\Phi (t_f)⟩=\tilde{U}(I_{\mathcal{A}}\otimes {\mathcal{U}}_{\mathcal{S}})U|\varphi ⟩\otimes |\xi ⟩$$

$$|\Phi (t_f)⟩=\sum _{n,m}{D}_{nm}|\varphi ⟩\otimes {\tilde{\Pi }}_m{\mathcal{U}}_{\mathcal{S}}{\Pi }_n|\xi ⟩$$

$${\rho }_{\mathcal{A}}(t_f)=\sum _{n,n^{\prime },m}\mathrm{Tr}\left({\tilde{\Pi }}_m{\mathcal{U}}_{\mathcal{S}}{\Pi }_n|\xi ⟩⟨\xi |{\Pi }_{n^{\prime }}{\mathcal{U}}_{\mathcal{S}}^{†}\right)D_{nm}|\varphi ⟩⟨\varphi |D_{n^{\prime }m}^{†}.$$

$${\rho }_{\mathcal{A}}(t_f)=\sum _{n,n^{\prime },m}\mathrm{Tr}\left({\tilde{\Pi }}_m{\mathcal{U}}_{\mathcal{S}}{\Pi }_n{\rho }_{\mathcal{S}}{\Pi }_{n^{\prime }}{\mathcal{U}}_{\mathcal{S}}^{†}\right)D_{nm}{\rho }_{\mathcal{A}}{D}_{n^{\prime }m}^{†}.$$

$$P_{\mathcal{A}}(p)=\sum _{n,n^{\prime },m}\mathrm{Tr}\left({\tilde{\Pi }}_m{\mathcal{U}}_{\mathcal{S}}{\Pi }_n{\rho }_{\mathcal{S}}{\Pi }_{n^{\prime }}{\mathcal{U}}_{\mathcal{S}}^{†}\right)⟨p-\lambda w_{nm}|{\rho }_{\mathcal{A}}|p-\lambda w_{n^{\prime }m}⟩.$$

$$\sum _{n,m}\mathrm{Tr}\left({\tilde{\Pi }}_m{\mathcal{U}}_{\mathcal{S}}{\Pi }_n{\rho }_{\mathcal{S}}{\Pi }_n{\mathcal{U}}_{\mathcal{S}}^{†}\right)⟨p-\lambda w_{nm}|{\rho }_{\mathcal{A}}|p-\lambda w_{nm}⟩=\int dwP(w)f(p-\lambda w),$$

$$f(p)=⟨p|{\rho }_{\mathcal{A}}|p⟩.$$

$$\left|\sum _{n\ne n^{\prime },m}\mathrm{Tr}\left({\tilde{\Pi }}_m{\mathcal{U}}_{\mathcal{S}}{\Pi }_n{\rho }_{\mathcal{S}}{\Pi }_{n^{\prime }}{\mathcal{U}}_{\mathcal{S}}^{†}\right)⟨p-\lambda w_{nm}|\varphi ⟩⟨\varphi |p-\lambda w_{n^{\prime }m}⟩\right|\le \sum _{n\ne n^{\prime },m}\left|\mathrm{Tr}\left({\tilde{\Pi }}_m{\mathcal{U}}_{\mathcal{S}}{\Pi }_n{\rho }_{\mathcal{S}}{\Pi }_{n^{\prime }}{\mathcal{U}}_{\mathcal{S}}^{†}\right)\right|\frac{\sigma }{\hslash \sqrt{\pi }}{e}^{-\frac{{\sigma }^2{\lambda }^2}{4{\hslash }^2}(E_n-E_{n^{\prime }}{)}^2}$$

$$P_{\mathcal{A}}(p)=\int dwP(w)f(p-\lambda w),$$

$$z_{nm}(t)=-w_{nm}\frac{\lambda t}{m_a}+\frac{g}{2}{t}^2\mathrm{and}\Delta z_t=\frac{\sigma }{\sqrt{2}}\sqrt{1+{\left(\frac{\hslash t}{m_a{\sigma }^2}\right)}^2}.$$

Full text

I Introduction

Classical fluctuation theorems establish surprising relations between non-equilibrium and equilibrium concepts. In particular, the work performed on a system during non-equilibrium processes is connected with key concepts of equilibrium thermodynamics, such as the free-energy jarzynski1997nonequilibrium ; crooks1999entropy . These relations have been verified in various experiments involving microscopic thermodynamic systems hummer2001free ; liphardt2002equilibrium ; collin2005verification . Recent advances in quantum technologies enable the control of small quantum systems that can be manipulated far from the regime where the usual thermodynamical laws are obeyed. This triggered the development of the rapidly growing field of non-equilibrium quantum thermodynamics esposito2009nonequilibrium ; campisi2011colloquium ; seifert2012stochastic ; goold2016role .

When quantum fluctuations dominate, defining and measuring work and heat, two central concepts in classical thermodynamics, is non-trivial. For driven, but otherwise isolated, quantum systems, work $w$ is a random variable associated with the change in the internal energy talkner2007fluctuation , as the first law of thermodynamics indicates. Thus, the commonly accepted definition of quantum work requires a two-time measurement strategy, which consists of performing two projective energy measurements, one at the beginning and the other at the end of the process. Then, work is associated with the measured energy difference. However, implementing the two-time measurement is experimentally difficult huber2008employing ; watanabe2014generalized due to the fact that the two projective measurements are unavoidably disruptive (see Ref. an2015experimental for an ion trap implementation). Alternative methods to evaluate the work probability distribution that rely on the direct estimation of its Fourier transform were also proposed in Refs. mazzola2013measuring ; dorner2013extracting and later implemented in NMR experiments batalhao2014experimental .

In this paper we present the design and the experimental implementation of a “quantum work meter” (QWM) operating on an ensemble of cold atoms, combining the idea presented in Ref. roncaglia2014work and the experimental setup used in Ref. machluf2013coherent . Our QWM is conceptually different from previous work-measurement devices. Its main advantage is that the QWM efficiently samples $P(w)$, which is a direct observable in the experiment. Namely, our QWM not only directly measures work but also directly samples its probability distribution $P(w)$ [i.e. the outcome $w$ is obtained with probability $P(w)$]. As the work probability distribution plays a central role in the fluctuation theorems of non-equilibrium quantum thermodynamics, the QWM is an ideal tool to test their validity. In particular, we use it to verify the Jarzynski identity jarzynski1997nonequilibrium ; tasaki2000jarzynski ; kurchan2000quantum ; mukamel2003quantum ; talkner2007fluctuation .

II Results

Work measurement and the QWM. A QWM is an apparatus that measures the work performed on a driven quantum system whose Hamiltonian varies from an initial $H$ to a final $\tilde{H}$ with eigenvalues $E_{n}$ and $\tilde{E}_{m}$, respectively. For an isolated system $\mathcal{S}$, with a $D$-dimensional space of states, the number of different values of work is bounded by $D^{2}$. Therefore, the pointer of the QWM has $D^{2}$ distinct positions (one for each value of $w=w_{nm}=\tilde{E}_{m}-E_{n}$). The QWM presented here enables us to choose $H$ and $\tilde{H}$ (fixing the possible values of $w$) and to vary the intermediate driving (inducing different evolution operators denoted as $\mathcal{U}_{\mathcal{S}}$). In this way, we vary the probability $P(w)$, which depends on the intermediate driving $\mathcal{U}_{\mathcal{S}}$.

By sampling $P(w)$, we use the QWM to verify a fundamental result in non-equilibrium quantum thermodynamics: the Jarzynski identity. This identity states that for any initial state with populations identical to the ones associated to a thermal Gibbs state and for any distribution $P(w)$, the linear combination ${\langle e^{-\beta w}\rangle=\sum_{w}e^{-\beta w}P(w)}$, where $\beta=1/k_{B}T$ is the inverse temperature of the system, is an equilibrium property (rather than a non-equilibrium one). The Jarzynski identity (see the Supplementary Note 1) reads

[TABLE]

where $\Delta F$ is the free energy difference between the thermal states associated with the Hamiltonians $H$ and $\tilde{H}$. In the absence of degeneracies, this implies that the vector formed by the $D^{2}-1$ measured probabilities belongs to a $D^{2}-2$ dimensional hyperplane: the “Jarzynski manifold” [as shown in the Supplementary Note 1, further constraints restrict this dimensionality to $(D-1)^{2}$]. With the QWM we measure $P(w)$ for different driving fields showing that all probability vectors belong to the same manifold. By characterising this manifold, we not only verify the identity but also independently estimate the free energy difference $\Delta F$ jarzynski1997nonequilibrium ; hummer2001free ; liphardt2002equilibrium ; collin2005verification .

The work distribution sampled by the QWM tasaki2000jarzynski ; kurchan2000quantum ; mukamel2003quantum ; talkner2007fluctuation is:

[TABLE]

Thus, $P(w)$ is the probability density of finding the energy difference $w$ after a measurement of $H$ followed by an intermediate driving $\mathcal{U}_{\mathcal{S}}$ and a final measurement of $\tilde{H}$. This is indeed the case if $p_{n}$ is the probability of obtaining $E_{n}$ when measuring $H$ and $p_{m|n}$ is the probability of obtaining $\tilde{E}_{m}$ when measuring $\tilde{H}$ given that $E_{n}$ was detected at the beginning. Equation  (3) defines a probability density that is independent of the initial coherences in the energy basis. For the discrete $D^{2}$ values of $w$ we will use $P(w)$ to denote the probability (not the density) of each $w$. The concept on which our QWM is based was first discussed in Refs. roncaglia2014work ; de2015measuring , where it was noticed that the work done on $\mathcal{S}$, can be detected by performing a generalised quantum measurement, which enables the number of outcomes to be larger than $D$. This can be done by entangling $\mathcal{S}$ with an ancilla $\mathcal{A}$ that stores a coherent record of $w$. Then a standard measurement on $\mathcal{A}$ can reveal $w$. Similar strategies have been later studied and extended to other contexts in Refs. de2015measuring ; talkner2016aspects ; anders2017work .

Design and operation of the QWM. A pictorial representation of the protocol we follow to operate the QWM is shown in Fig. 1a. The QWM is designed to measure the work done on a system $\mathcal{S}$ whose Hamiltonian changes from $H$ to $\tilde{H}$ and which is subjected to a driving $\mathcal{U}_{\mathcal{S}}$ in between. We couple $\mathcal{S}$ to a continuous variable system $\mathcal{A}$ and use $\hat{z}_{\mathcal{A}}$ to denote its position (the generator of translations along the momentum $p$). A coherent record of $w$ is created by an “entangling interaction” between $\mathcal{A}$ and $\mathcal{S}$ that must take place before and after the driving $\mathcal{U}_{\mathcal{S}}$. The unitary operators representing these interactions are: $U=e^{-i\lambda\,\hat{z}_{\mathcal{A}}\otimes{H}/\hbar}$ and $\tilde{U}=e^{i\lambda\,\hat{z}_{\mathcal{A}}\otimes{\tilde{H}}/\hbar}$, where $\lambda$ is a coupling parameter. Thus, $U$ and $\tilde{U}$ respectively translate $\mathcal{A}$ along $p$ by a displacement proportional to $(-\lambda H)$ and $\lambda\tilde{H}$. Then, as shown in detail in the Supplementary Note 3, the final measurement of $p$ on $\mathcal{A}$ yields a random result whose distribution $P_{\mathcal{A}}(p)$ is a smeared version of the true work distribution $P(w)$ defined in Eq. (3). In fact, outcome $p$ is obtained with a probability density $P_{\mathcal{A}}(p)=\int dwP(w)f(p-\lambda w)$, where the window function $f(p)=|\langle p|\phi\rangle|^{2}$ is fixed by $|{\phi}\rangle$, the initial state of $\mathcal{A}$ [thus, by localising $|{\phi}\rangle$ we improve the accuracy in the estimation of $P(w)$].

A “universal” QWM is an apparatus which can measure $w$ and sample $P(w)$ for any possible choice of $H$ and $\tilde{H}$. To build it, we need enough control to enforce the entangling operators $U$ and $\tilde{U}$ for any choice of $H$ and $\tilde{H}$. Remarkably, this is achieved for a $2$-level system by the atom chip implementation we describe below.

Experimental implementation of the QWM. To describe our QWM we should identify the physical systems representing $\mathcal{S}$ and $\mathcal{A}$, the way in which $H$ and $\tilde{H}$ can be chosen, and how the associated $U$ and $\tilde{U}$ are implemented. In our experiment we represent $\mathcal{S}$ by the subspace associated with the Zeeman sublevels ${|{1}\rangle\equiv|{F=2,m_{F}=1}\rangle}$ and $|{2}\rangle\equiv|{F=2,m_{F}=2}\rangle$ of a 87Rb atom that, as in Ref. machluf2013coherent , behaves as a two-level system (see below). The motional degree of freedom of the atom plays the role of $\mathcal{A}$.

A key element of the QWM presented here is the atom chip keil2016fifteen , which efficiently entangles the internal and motional degrees of freedom of an atom just $\sim 100\,\mu$m away from the chip surface, through short and strong Stern-Gerlach type magnetic gradient pulses. These pulses are generated using a $3$-current-carrying-wire setup on the chip surface (described in Ref. margalit2015self and Methods). A gradient pulse along the $z$ direction with amplitude $B^{\prime}$ and duration $\tau$, induces a momentum kick $m_{F}\,\delta p$ on an atom in the $m_{F}$ state ($\delta p\sim\mu_{B}g_{F}B^{\prime}\tau$, where $\mu_{B}$ and $g_{F}$ are, respectively, the Bohr magneton and the Landé factor machluf2013coherent ). The evolution of the state of the atom induced by such a pulse is described by the unitary operator $U_{p}=e^{i\,\delta p\,\hat{z}_{\mathcal{A}}\otimes\hat{\sigma}/\hbar}$, where the operator ${\hat{\sigma}=|{1}\rangle\langle{1}|+2|{2}\rangle\langle{2}|}$. This physical operation translates $\mathcal{A}$ along the momentum $p$ by a displacement $\delta p\,\hat{\sigma}$ (notice that the operator $\hat{\sigma}$ defines the magnetic dipole moment of the atom since $\hat{\sigma}=\sum_{m=1,2}m|{m}\rangle\langle{m}|$). As described below, we apply two gradient pulses with different amplitudes ($B^{\prime}$ and $\tilde{B}^{\prime}$) and different durations ($\tau$ and $\tilde{\tau}$). Thus, defining $H=E\hat{\sigma}$ and ${\tilde{H}}={\tilde{E}}\hat{\sigma}$, $U_{p}$ and $\tilde{U}_{p}$ implement the required entangling operation $U$ and $\tilde{U}$, respectively. In this implementation $\lambda$ is consequently replaced by $-\delta p/E$ and $\delta\tilde{p}/\tilde{E}$, enforcing $\tilde{E}/E=-\delta\tilde{p}/\delta p$. The momentum kicks induced by both pulses are controlled in the experiment, and consequently, by fixing their ratio, we can simulate an arbitrary system with initial and final Hamiltonians $H$ and $\tilde{H}$ which are characterised by $\tilde{E}/E$ having the same ratio. Finally, let us note that the two pulses utilise $B^{\prime}$ and $\tilde{B}^{\prime}$ with opposite signs to ensure that the sequence creates a record of work corresponding to $\tilde{E}_{m}-E_{n}$.

To achieve universality we only need to be able to fix the energy splitting $E$ and $\tilde{E}$ of $H$ and $\tilde{H}$, as well as their eigenbasis. The traces of $H$ and $\tilde{H}$ (the sum of their eigenvalues) do not affect $P(w)$ but only add a constant to all values of $w$. As arbitrary $E$ and $\tilde{E}$ can be simulated and any change of basis can be absorbed into $\mathcal{U}_{\mathcal{S}}$, we conclude that our atom chip QWM can sample $P(w)$ for an arbitrary 2-level system and is thus universal.

The 87Rb atoms are magnetically trapped in state $|{2}\rangle$ and evaporatively cooled to a Bose-Einstein condensation (BEC). The BEC is released from the trap and a radio-frequency (RF) pulse is used to prepare a superposition of $|{1}\rangle$ and $|{2}\rangle$. A strong homogeneous magnetic field (created by external coils) suppresses the transitions taking $|{1}\rangle$ into the $|{2,0}\rangle$ state (due to the non-linear Zeeman effect machluf2013coherent ). The initial populations ($p_{1}$ and $p_{2}$) are chosen so that $\beta E=\ln{(p_{1}/p_{2})}$. The initial motional state is a wave-packet localised in position and momentum.

It should be noted that the initial internal state of the atom, while having the same populations as defined by the temperature of a thermal state, is still a pure state. However, the quantum coherences of this initial state do not affect the results of the QWM. As explained in the Supplementary Note 3, the contribution of the initial coherences to the final probability is multiplied by the overlap between the motional states of the atom associated with the different values of work. Thus, when the atomic clouds associated with the different work values are well separated, the effect of initial coherences is negligible. In this regime, our experiment gives the same result as the one we would obtain by preparing an initial thermal state (with no coherences). The study of the importance of the initial coherences in the definition of work is an interesting topic in itself, which is beyond the scope of our paper (see, for example, Ref. solinas2015full ; kammerlander2016coherence ; perarnau2017no ; anders2017work ; sampaio2017impo ).

The experimental sequence, presented in Fig. 1b, is: i) prepare the initial state and release the cloud (which then freely falls along $z$, the direction of gravity), ii) apply the magnetic gradient $U$ along $z$, iii) apply the driving $\mathcal{U}_{\mathcal{S}}$ by exposing the atoms to a RF field resonant with the Zeeman splitting induced by the homogeneous bias field, iv) apply the gradient $\tilde{U}$, v) obtain an image of the four clouds after a time-of-flight and count the number of atoms in each cloud. More details of the experiment can be found in Methods and the Supplementary Note 3. For the experimental demonstration presented here we set the ratio between the measured momentum kicks induced by the two pulses to $-\delta\tilde{p}/\delta p=0.56\pm 0.02$. Hence, our realisation of the QWM samples the work distribution of a simulated system in which the energy splitting is reduced to $56\%$ of its original value, from $E$ to $\tilde{E}$, while driven by $\mathcal{U}_{\mathcal{S}}$.

Fig. 1c shows a typical image obtained by the QWM. Four clouds are visible. From the positions of the center of each cloud, $\bar{z}$, we infer the total momentum shift, $\bar{p}$, induced by the pulses on that cloud (we take into account both the free fall and the kicks induced by the pulses, see the Supplementary Note 3). Then, we obtain the corresponding value of work as $w=E\,\bar{p}/\delta p$ ($w$ is proportional to $E$, whose value, together with the experimental results, determine the work $w$). Furthermore, the probability $P(w)$ for each $w$ is directly measured by the number of atoms in each cloud. Notably, this experiment determines the entire $P(w)$ distribution in a single shot.

Testing the Jarzynski identity. We repeat the experiment fixing the timing, duration and pulse strength. We consider three initial $\beta$’s and vary the intermediate driving $\mathcal{U}_{\mathcal{S}}$ by changing the duration of the RF field. In this way, we obtain many sets of probability distributions, each of which defines a 3D–vector (as there are three independent probabilities). When we represent all these vectors in the same 3D-plot, we see that they all belong to the same $\beta$–dependent manifold. Fig. 2a shows that this manifold is a $\beta$–dependent line (the dimensionality of this “Jarzynski manifold” is $(D-1)^{2}$, which in this case equals $1$).

Using the measured work probabilities we calculate the exponential average of the work $\langle e^{-\beta w}\rangle$ for each driving field. Fig. 2b displays the value of ${G=-\ln[\langle e^{-\beta w}\rangle]=-\ln[\sum_{w}e^{-\beta w}P(w)]}$ as a function of the duration of the intermediate RF field, that parametrises $\mathcal{U}_{\mathcal{S}}$. As established by the Jarzynski identity, $G$ is independent of the driving field and only depends on $\beta$. The horizontal lines in Fig. 2b are the theoretically predicted values of $\beta\,\Delta F$, obtained from a direct calculation (with its own theoretical uncertainty, due to the error in the estimation of $\beta\,E$). This calculation simply involves computing the initial and final partition functions, respectively denoted as $Z$ and $\tilde{Z}$, and using the identity $\beta\,\Delta F=\ln(Z/\tilde{Z})$. We find that, as the Jarzynski identity establishes, $G=\beta\Delta F$. From Fig. 2b one can notice that the largest errors in the estimation of $\beta\Delta F$ appear for $\beta\,E=1.75$. In this case, $P(w)\lesssim 0.1$ for two values of $w$ and, therefore, the relative error in the atom number estimation is large, inducing a larger error in the estimation of $\beta\Delta F$.

In Table 1 we compare measured and estimated values of $\beta\Delta F$. The uncertainty in the estimation of $\beta\Delta F$ and $\Delta F$ is close to 10%, which is enough to distinguish the three values of $\beta\Delta F$. On the other hand, in the case of $\Delta F$, there is a significant overlap in the measured values which does not allow to properly distinguish between the three different cases due to the error in the estimation of $\beta E$.

III Discussion

We presented and implemented a QWM, a new device directly sampling the work distribution on an ensemble of cold atoms. Our QWM can be used to simulate the behaviour of an arbitrary $2$-level system. We implemented it with an atom chip and verified the Jarzynski identity over a wide range of non-equilibrium processes. This is the first experiment, and so far the only one, directly sampling $P(w)$ offering advantages and different perspectives over previous work measurement schemes. Remarkably, in this cold atom experiment, the QWM extracts full statistical information about the work distribution in a single shot.

IV Methods

Initial state preparation.  After preparing the BEC, a homogeneous magnetic field of $36.7\,$G ($25\,h\,\mathrm{MHz}/\mu_{B}$, where $h$ is Planck’s constant) is used to push the transition to $|{2,0}\rangle$ out of resonance by $\sim 180\,$kHz due to the non-linear Zeeman effect, which is larger than the power broadened driving RF field of $\mathcal{U}_{\mathcal{S}}$. This ensures that the atoms behave as $2$-level systems. The BEC is released from the trap and a RF pulse is used to prepare a superposition of $|{1}\rangle$ and $|{2}\rangle$. By varying the relative populations we consider three different pseudo-thermal states. The initial motional state is a wave packet $|{\phi}\rangle$, well localised at $z_{0}=91\pm 1.2\,\mu$m from the chip and momentum $\sim 0$.

Entangling operations and measurement. An inhomogeneous magnetic field is used to couple spin and motional degrees of freedom. This is generated by a current $I=0.85\,$A in the 3-wire setup during a time $\tau$. The three parallel gold wires lie on the $x$ direction of the chip surface (Fig. 1.b). They are $10$ mm long, $40\,\mu$m wide and $2\,\mu$m thick. Their centers are at $y=-100,0,100\,\mu$m and the same current run through them in alternating directions ($-I,I,-I$, respectively), creating a $2$D quadrupole field at $z=100\,\mu$m below the chip. After a time of flight of $2.4$ ms the atoms are at $z\sim 119\,\mu\mathrm{m}$. At this point the first gradient pulse implements $U$: $\tau=40\,\mu$s with an amplitude of $B^{\prime}\sim 95\,$G/mm, such that the momentum kick is along $+z$. Then, after 3.1 ms the atoms are at $\tilde{z}\sim 0.3\,$mm and the second gradient pulse implements $\tilde{U}$: $\tilde{\tau}=300\,\mu$s, $\tilde{B}^{\prime}\sim-7.5\,$G/mm, such that the momentum kick is along $-z$. The relative strengths of the spin-dependent forces sets the energy splitting of the Hamiltonians which in this case is on average $\tilde{E}/E=-\tilde{\delta}p/\delta p=0.56\pm 0.02$ (this is the measured value, that takes into account fluctuations in the initial position of the cloud and in the gradient pulses). In between the entangling operation $\mathcal{U}_{\mathcal{S}}$ is applied with a RF pulse. Finally, an image of the atomic clouds is obtained after a time-of-flight of $18.2\,$ms after the $2$nd gradient (the clouds are centered around $z\sim 3\,$mm). The position and number of atoms of each cloud is determined. The momentum shifts of each cloud (that codifies the value of $w$) are obtained from the difference in positions between the clouds, that follow approximately classical trajectories (see Supplementary Note 3).

Uncertainties. The main source of position error is the initial distance of the cloud from the atom chip, whose uncertainty is $\sim 1\%$. This error is later translated to momentum uncertainty, since the field gradients are position dependent. The field gradients have a fractional uncertainty of $10^{-3}$ due to current fluctuations machluf2013coherent . The central position of each cloud is estimated by fitting a Gaussian profile. Each work probability is estimated as a normalised sum of the measured optical density in a relevant region around the cloud, introducing probability uncertainty (due to atom numbers uncertainty). Our $\sim 5\,\mu\mathrm{m}$ optical resolution also induces an error in the determination of the position for each cloud. We perform three different runs for each combination of initial state population ratios and intermediate driving and use the average values of position and probability. This gives us a position uncertainty of $\sim 0.015\,\mathrm{mm}$ and a probability uncertainty of $\sim 0.015$ (standard error).

Acknowledgments. FC, AJR and JPP acknowledge financial support from ANPCyT (PICT 2013-0621 and PICT 2014-3711), CONICET and UBACyT. YM and RF gratefully acknowledge funding by the Israel Science Foundation, the EC Matter-Wave consortium [FP7-ICT-601180], and the German DFG through the DIP program [FO703/2-1]. SM acknowledges financial support by the Foundation for Fundamental Research on Matter (FOM), which is part of the Netherlands Organisation for Scientific Research (NWO). We also thank the BGU nano-fabrication team for making available the high quality chip, and Zina Binstock for helping with the electronics.

Appendix A Supplementary Material.

Supplementary Note 1 - Jarzynski identity and quantum work measurement

Work is the energy variation induced on a system $\mathcal{S}$ by a certain driving (the system is otherwise isolated). Suppose that the Hamiltonian of $\mathcal{S}$ changes from an initially $H$ to a final $\tilde{H}$ and that the system is driven by a certain unitary operator $\mathcal{U}_{\mathcal{S}}$ in between. We denote the eigenvalues and eigenvectors of $H$ as $E_{n}$ and $|{\varphi_{n,\alpha}}\rangle$ (where $\alpha$ labels different eigenstates with the same energy). In turn the eigenvalues and eigenvectors of $\tilde{H}$ are denoted as $\tilde{E}_{m}$ and $|{\tilde{\varphi}_{m,\gamma}}\rangle$ (again, $\gamma$ labels states with the same energy $\tilde{E}_{m}$). Suppose that we measure the energy at the initial and final times. The work probability distribution $P(w)$ is nothing but the probability density to obtain a value $w$ as the difference between the results of the two energy measurements. This can obviously be computed as

[TABLE]

where $p_{n}$ is the probability of initially measuring energy $E_{n}$ and $p_{m|n}$ is the probability of measuring energy $\tilde{E}_{m}$ at the end of the driving given that $E_{n}$ was detected at the beginning. If the initial state of $\mathcal{S}$ is $\rho$, the above probabilities can be simply written in terms of the projectors $\Pi_{n}=\sum_{\alpha}|{\varphi_{n,\alpha}}\rangle\langle{\varphi_{n,\alpha}}|$ and $\tilde{\Pi}_{n}=\sum_{\gamma}|{\tilde{\varphi}_{n,\gamma}}\rangle\langle{\tilde{\varphi}_{n,\gamma}}|$. Thus,

[TABLE]

Jarzynski identity follows immediately from the above definition of the work probability distribution. In fact, if we compute the exponential average of the work, we get

[TABLE]

If the initial state is thermal, then $\rho=\sum_{n}\Pi_{n}e^{-\beta E_{n}}/Z$ (where $Z$ is the partition function $Z=\sum_{n}g_{n}e^{-\beta E_{n}}$ with $g_{n}=\mathop{\mathrm{Tr}}\nolimits(\Pi_{n})$ being the degeneracy of level $E_{n}$). Then, we can replace $p_{n}=e^{-\beta E_{n}}/Z$ and perform the summation over the label $n$ by noticing that $\sum_{n}p_{m|n}=\mathop{\mathrm{Tr}}\nolimits(\tilde{\Pi}_{m})=\tilde{g}_{m}$. In this way we obtain

[TABLE]

where the partition function of the final Hamiltonian is $\tilde{Z}=\sum_{m}\tilde{g}_{m}e^{-\beta\tilde{E}_{m}}$, with $\tilde{g}_{m}=\mathop{\mathrm{Tr}}\nolimits(\tilde{\Pi}_{m})$ the degeneracy of level $\tilde{E}_{m}$. The final step that leads to the Jarzynski identity is simply to notice that for a Gibbs state we have $e^{-\beta\Delta F}=\tilde{Z}/Z$. Thus, the average exponential work becomes independent of the intermediate driving process $\mathcal{U}_{\mathcal{S}}$ and turns out to be determined by equilibrium properties. Thus, $\langle e^{-\beta w}\rangle=e^{-\beta\Delta F}$.

When $\mathcal{S}$ has a finite dimensional Hilbert space, $w$ can only take discrete values. In this case, instead of using the probability density we use directly the probability for each value of $w$. Thus, we can arrange $P(w)$ as a vector with $D^{2}$ components. Taking into account the normalisation condition, the vector of independent probabilities is $D^{2}-1$ dimensional. In turn, being the Jarzynski identity a linear equation in terms of $P(w)$ it constraints the probability vector to belong to a $D^{2}-2$ dimensional hyperplane. But, of course, there are further constrains reducing the number of independent probabilities. The simplest way to obtain this number is to go back to the definition of $P(w)$ in Supplementary Equation (3) and count the number of free parameters we have. In this equation, the probabilities $p_{n}$ are fixed by the initial state. Then, the number of independent probabilities is determined by the number of independent parameters in $p_{m|n}$. For the non-degenerate case we consider here (where all the values of work are different), the calculation is simple. In fact, the coefficients $p_{m|n}$ form a doubly stochastic matrix (since they are all positive numbers such that $\sum_{n}p_{m|n}=\sum_{m}p_{m|n}=1$). For such square matrix of dimension $D$, there is always $(D-1)^{2}$ free parameters. Indeed, this is the dimensionality of the manifold of where the probability vector lies. Jarzynski identity establishes that this manifold is a $\beta$-dependent hyperplane (a line in our case, where $D=2$).

Supplementary Note 2 - Work measurement as a POVM

Let us consider a system $\mathcal{S}$ with a $D$-dimensional space of states and show, in a simple way, that the work measurement can be viewed as a generalised quantum measurement. For this, we start by writing the probability for a given value of work $w_{nm}=\tilde{E}_{m}-E_{n}$ as

[TABLE]

Using the formula for the transition probability, we find that

[TABLE]

where

[TABLE]

It is simple to verify that the operators $A_{nm}$ expand the identity as $\mathbb{I}=\sum_{n,m}A_{nm}$ and that they are positive semi-definite (i.e., that for any state $|{\chi}\rangle$ we have $\langle{\chi}|A_{nm}|{\chi}\rangle\geq 0$). Therefore, the operators $A_{nm}$ define a positive operator valued measure (POVM), which is the most general type of measurement one can perform in quantum mechanics. Therefore, work can be measured in the same way as any POVM can: A powerful result (Neumark’s theorem) establishes that any POVM can be realised by coupling the system $\mathcal{S}$ with an ancillary system $\mathcal{A}$ and then performing a standard projective measurement on $\mathcal{A}$. This measurement can be performed at a single time. Thus, surprisingly, the two-time work measurement strategy can be replaced by a single-time strategy (which is the basic idea exploited by our QWM). In the following section we show how one can construct an approximation for that ideal apparatus, that we call Quantum Work Meter (QWM).

Supplementary Note 3 - Probability distribution for the outcome of a QWM

Here we compute the probability distribution for the result of the measurement of the auxiliary register of a general Quantum Work Meter. The protocol defining the apparatus is shown in Figure 1 of the main text. A system $\mathcal{S}$ is coupled to an ancillary one $\mathcal{A}$. This ancilla is a continuous variable system (of course this can be relaxed). The system $\mathcal{S}$ and the ancilla $\mathcal{A}$ are subject to the following evolution: $i)$ an entangling interaction $U$ is applied (which correlates $\mathcal{S}$ and $\mathcal{A}$), $ii)$ the evolution $\mathcal{U}_{\mathcal{S}}$ is applied on the system $\mathcal{S}$, $iii)$ a second entangling interaction $\tilde{U}$ is applied. Finally, after this sequence $\mathcal{A}$ is measured. The initial state of the system formed by $\mathcal{S}$ will be assumed to be a product state. For simplicity, we first assume that the states are pure and denote them as $|{\xi}\rangle$ (the state of $\mathcal{S}$) and $|{\phi}\rangle$ (the state of $\mathcal{A}$). We will later generalise the result for an initial state which is a tensor product of arbitrarily mixed states. After the sequence of operations we described above, the total final state is:

[TABLE]

The nature of the two entangling operations, that was discussed in the main text, is such that they both induce translations of $\mathcal{A}$ which depend on the state of $\mathcal{S}$. More specifically, $U=e^{-\frac{i}{\hbar}\lambda\hat{z}_{\mathcal{A}}\otimes H}$ and $\tilde{U}=e^{\frac{i}{\hbar}\lambda\hat{z}_{\mathcal{A}}\otimes\tilde{H}}$, where we use $\hat{z}_{\mathcal{A}}$ to denote the generator of translations of $\mathcal{A}$ along a certain variable $p$. Using this, it is simple to rewrite the final state as

[TABLE]

In this equation the displacement operators $D_{nm}$ act on the states of $\mathcal{A}$ and are defined as $D_{nm}=e^{\frac{i}{\hbar}\lambda w_{nm}\hat{z}_{\mathcal{A}}}$.

The interpretation of the above equation is simple: After the sequence of operations, $\mathcal{S}$ and $\mathcal{A}$ become entangled in such a way that a record of $w_{nm}$ is stored in $\mathcal{A}$. The states $|{F_{nm}}\rangle\equiv D_{nm}|{\phi}\rangle$ are “flag states” associated with the different values of work. When these states are orthogonal, they can be unambiguously distinguished and the value of work can be retrieved. Below, we will consider a more realistic scenario where the initial state of $\mathcal{A}$ is a localised coherent state. In that case, the flag states are displaced coherent states (which are simply translated along the variable $p$ direction by an amount that is proportional to $w_{nm}$). These states are not strictly orthogonal, but have a finite overlap. This induces an error in the work estimation protocol. However, the error can be exponentially reduced by simply increasing the interaction strength $\lambda$ (as the overlap exponentially decreases with $\lambda$).

From the above expression it is simple to obtain the quantum state of $\mathcal{A}$ by computing its reduced density matrix (which is obtained from the total state by tracing out the system $\mathcal{S}$). Thus,

[TABLE]

This expression can be generalised to the case where the initial states of $\mathcal{S}$ and $\mathcal{A}$ are initially mixed. In fact, if $\rho_{\mathcal{S}}$ and $\rho_{\mathcal{A}}$ respectively denote the initial density matrices of $\mathcal{S}$ and $\mathcal{A}$, the final state of $\mathcal{A}$ is

[TABLE]

From the above equation we obtain the probability density for detecting the value $p$ in a measurement of $\mathcal{A}$. Thus,

[TABLE]

The contribution of the diagonal ($n=n^{\prime}$) and off-diagonal ($n\neq n^{\prime}$) terms play a different role in the above expression. In fact, it is simple to show that the diagonal contribution is a smeared version of the true work distribution. Thus,

[TABLE]

where the window function is

[TABLE]

Therefore, the off-diagonal terms of $\rho_{\mathcal{S}}$ are responsible for the error in the work estimation and should be made small for it to be accurate. It is simple to show that if the momentum wave function of the initial state $|{\phi}\rangle$ is a Gaussian with a momentum dispersion $1/\sigma$ ($\sigma$ is the position dispersion), then the off-diagonal terms are bounded by:

[TABLE]

Thus, by increasing $\lambda$ (the interaction strength) or $\sigma$ (the position dispersion of the initial state) we exponentially reduce the error in the work estimation. It is worth noting that by increasing $\sigma$ we reduce the momentum uncertainty and localise the initial state in momentum. Naturally, the method becomes precise when the initial localisation in momentum is much smaller than the difference between the first momentum kicks (which is fixed by the product $\lambda(E_{n}-E_{n^{\prime}})$). When these conditions are satisfied, the off-diagonal terms can be neglected and the probability density to detect $p$ is

[TABLE]

where the window function $f(p)$ is defined as $f(p)=|\langle{p}|{\phi}\rangle|^{2}$ (which, for a coherent state is simply ${f(p)=e^{-\frac{\sigma^{2}}{\hbar^{2}}p^{2}}\sigma/\hbar\sqrt{\pi}}\;$).

QWM using an atom chip

Here we consider the implementation of the QWM using a cloud of atoms in a chip. As described in the main text, $\mathcal{S}$ is the pseudo spin 1/2 associated with the $F=2$, $m_{F}=1,2$ hyperfine states of a ${}^{87}\rm{Rb}$ atom (which, as discussed in the main text, behaves as a two level atom). $\mathcal{A}$ is encoded in the motional degrees of freedom of each atom. There are two subtle differences between the ideal protocol for a QWM described in the previous section and the implementation in an atom chip. The first different concerns the final measurement. Thus, in the previous section we computed the probability for a final momentum measurement but in the real experiment we are forced to measure the atomic position by taking an image of the atomic clouds. Therefore we will show below that the position measurement enables us to determine work and sample $P(w)$. The second difference is that in the real experiment we should take into account the fact that the ancilla $\mathcal{A}$ evolves during the whole process because the atoms actually move (they freely fall along the vertical direction). Thus, the real protocol describing the experiment is shown in Fig. 3.

Let us now analyse this process. We can first neglect the free fall taking place between the two entangling operations (we take this into account later) and assume the initial state of $\mathcal{A}$ is a coherent state localised around initial values of position and momentum which we arbitrarily take as $z=0$ and $p=0$. We denote this state as $|{\phi}\rangle=|{0,0}\rangle$. The calculation presented in the above section should be slightly modified. In this case the flag states are $|{F_{nm}}\rangle=u^{(1)}_{fall}D_{nm}|{0,0}\rangle$. Taking into account that $u_{fall}=e^{-\frac{i}{\hbar}t(\hat{p}_{\mathcal{A}}^{2}/2m_{a}-m_{a}g\hat{z}_{\mathcal{A}})}$ (where $t$ is the duration of the free-fall, $m_{a}$ is the mass of the atoms, and $g$ the gravity acceleration) we can easily compute the expectation value of the position for each flag state (as well as the corresponding position dispersion). In fact, we have

[TABLE]

Therefore, we can notice that by measuring the final position of the atoms we can infer the value of the momentum before the fall and thus acquire information about work $w$. In fact, the difference between the positions of the clouds is proportional to the difference in the values of work. The price we have to pay, is that the spread of the wave packets increases during the free-fall.

In a more realistic description of the experiment, we need to include also the free fall between the entangling gates ($u_{\rm fall}^{(1)}$ in the Supplementary Figure 3). It is easy to verify that since $u_{\rm fall}^{(1)}\,\hat{z}_{\mathcal{A}}\,u_{\rm fall}^{(1)\dagger}=\hat{z}_{\mathcal{A}}+\hat{p}_{\mathcal{A}}\,t/m_{a}+\mathbb{I}_{\mathcal{A}}\,g\,t^{2}/2$ (where $t$ is the duration of the free fall), then $u_{\rm fall}^{(1)}\,e^{\frac{i}{\hbar}\lambda\,\hat{z}_{\mathcal{A}}\otimes\tilde{H}}\,u_{\rm fall}^{(1)\dagger}=e^{\frac{i}{\hbar}\lambda\,\hat{z}_{\mathcal{A}}\otimes\tilde{H}}e^{\frac{i}{\hbar}\frac{t}{m}\lambda\,\hat{p}_{\mathcal{A}}\otimes\tilde{H}}e^{\frac{i}{\hbar}\theta\,\mathbb{I}_{\mathcal{A}}\otimes\tilde{H}}$. Thus, the first term is the usual entangling operation and the last term is just a phase depending on the value of the energy. In turn, the second term is an entangling operation where the atom is displaced along position (instead of momentum) depending on the state of $\mathcal{S}$. In summary, the free-fall in between the entangling gates simply induces an extra translation of the atoms by an amount that depends on the final value of the energy.

Finally, we show in the Supplementary Figure 4 the distribution probability that we obtain in each experiment. In the plots we show the projections of the Jarzynski manifold, that appears in Fig. 2 (a) of the main text, onto the different axes.