Measuring work and heat in ultracold quantum gases

The study of out-of-equilibrium thermodynamics has received a significant thrust thanks to the experimental advances in the control and manipulation of microscopic systems. From a fundamental point of view, these endeavours aim at clarifying the foundations of modern thermodynamics and its connection to information theory. From a more applied perspective, these studies aim at understanding limitations of microscopic engines and building more efficient ones. Heat and work, two ubiquitous concepts in traditional thermodynamics, assume in this context the role of stochastic variables whose fluctuations can be ingeniously related to equilibrium properties, as is the case of the celebrated Jarzynski equality [1]. Many physical systems have been realized to investigate non-equilibrium thermodynamics, including for instance strands of RNA [2], single electron boxes [3], levitated or trapped nanoparticles [4, 5], and colloidal particles trapped in optical potentials [6].

In the last decade, general interest has been directed towards the quantum regime of out-of-equilibrium thermodynamics. In this regime, the dynamics of a small quantum system is dominated by quantum rather than thermal fluctuations. Although many open questions remain unanswered, some of the concepts of non-equilibrium classical thermodynamics have been translated into the quantum domain (see for example [7, 8]). A measure of work based on a two-measurement scheme is now commonly accepted [9] and can be shown, for isolated systems, to fulfill a quantum extension of the Jarzynski equality [10]. For open systems the Jarzynski equality still holds if one considers changes of energy in the system and the environment together [11]. However if we consider energy changes in the system only, the fluctuation relation for the system energy ceases to work and contains a correction that depends on the properties of the environment [12, 13]. In fact Jarzynski equality is still valid if the corresponding evolution superoperator is unital ⁴ , i.e., if the completely mixed state (corresponding to infinite temperatures) remains unaltered after the open system evolution.

Although implementing directly the two-measurement scheme has proven to be challenging, alternative routes to measure work in quantum systems have been proposed. One of these employs a Ramsey scheme [14, 15] and has been experimentally implemented in a nuclear magnetic resonance setup [16]. These proposals have also been extended to the open system scenario [17, 18]. Other proposals to measure work in the quantum domain relies on counting phonon excitations in trapped ions [19, 20] or counting electrons in single electron boxes [3].

Recently, two of us proposed a different method to measure work which is based on the fact that, for quantum systems, work can always be measured by performing a positive operator valued measure (POVM) at a single time [21]. This simple observation, that remained unnoticed until recently, implies that work can be measured with a single projective measurement on an extended system. Thus, it is always possible to devise a measurement apparatus that yields the work value W which is a random variable distributed with the work probability P(W). In this paper we will generalize this method and show how to use it to measure work and heat in gases of cold atoms.

There has been a lot of interest in applying ideas of non-equilibrium quantum thermodynamics in the case of isolated quantum many-body systems [22–30]. Despite the experimental advances in the field of ultracold atoms, an ideal platform for the quantum simulation of many-body systems [31], an experimentally feasible proposal for measuring heat and work in these systems is still missing. The Ramsey scheme mentioned earlier is based on the global coupling of an auxiliary two-level system with the system under consideration and might not be well suited for a cold atomic system.

The proposal we present to measure work and heat in quantum gases generalizes the method proposed in [21] and consists in coupling the atoms with a continuous degree of freedom which can be realized by the light quadratures. The interaction will be chosen in such a way to induce a phase-space translation of the continuous variable position that is conditional on the value of the energy of the system under consideration. In short, the method consists of three steps: first we let atoms interact with light in such a way that correlations between them are established. Second, while light is stored in a quantum memory, we drive the atoms with the thermodynamic process we are interested in. Third, we retrieve the light beam from the memory and redirect it into the atomic ensemble enforcing a second interaction between them. After these three steps, a standard homodyne detection of the output light is performed. The key of the method is that the statistical distribution of work and heat on the atoms if fully encoded in the statistical distribution of the light quadratures.

The paper is organized as follows: in section 2 we present the key ingredients of the method, which generalizes the one presented in [21]. Then, in the following sections we showcase two cold atoms settings where our proposal can be implemented using a quantum non-demolition measurement based on the Faraday effect [32]. The first one is designed for measuring work in cold atomic ensembles and is described in sections 3 and 4; the second example is for ultracold atoms in optical lattices, described in section 5 where we show how to measure heat and work for the atoms. In the latter case, the measurement scheme allows us to discern if the open system dynamics is unital or not, by checking whether the Jarzynski equality is fulfilled. Finally in section 6 we summarize.

Let us consider a process where a quantum system with an initial state ρ is driven from an initial Hamiltonian H to a final one $\tilde{H}$. The work value W in each realization is defined as the energy difference $W={{\tilde{E}}_{m}}-{{E}_{n}}$, where E_n are the eigenvalues of H (i.e. $H|{{\phi }_{n}}\rangle ={{E}_{n}}|{{\phi }_{n}}\rangle$) and those of $\tilde{H}$ are denoted with ${{\tilde{E}}_{m}}$ (i.e., $\tilde{H}|{{\tilde{\phi }}_{m}}\rangle ={{\tilde{E}}_{m}}|{{\tilde{\phi }}_{m}}\rangle$). Thus, W is a random variable distributed according to the following probability distribution:

$$\begin{eqnarray*}&&P(W)=\mathop{\sum }\limits_{m,n}{{p}_{n}}{{p}_{m|n}}\;\delta [W-({{\tilde{E}}_{m}}-{{E}_{n}})],\end{eqnarray*}$$

where ${{p}_{n}}=\langle {{\phi }_{n}}|\rho |{{\phi }_{n}}\rangle$ and ${{p}_{m|n}}=|\langle {{\tilde{\phi }}_{m}}|{{U}_{E}}|{{\phi }_{n}}\rangle {{|}^{2}}$ and U_E is the unitary operation that represents the driving. As we mention in the introduction, there are many protocols that were proposed to experimentally reconstruct this probability distribution.

Recently an alternative method that allows to sample the work probability distribution has been put forward in [21]. The method is based on the idea that work measurement is actually a POVM. As it is well known, any such generalized measurement can be implemented as a standard projective measurement on an enlarged system. A simple example of such strategy to implement the work POVM is depicted in figure 1. We assume that a system $\mathcal{S}$ is coupled to an auxiliary system $\mathcal{A}$ in such a way that $\mathcal{A}$ gets entangled with $\mathcal{S}$ keeping a coherent record of the energy at two times. In the simplest case (which will be generalized below), the interaction between $\mathcal{S}$ and $\mathcal{A}$ is such that it can be described by the unitary evolution operators ${{U}_{I}}={{{\rm e}}^{-{\rm i}\kappa PH}}$and ${{\tilde{U}}_{I}}={{{\rm e}}^{-{\rm i}\kappa P\tilde{H}}}$ .

Zoom In Zoom Out Reset image size

The auxiliary system $\mathcal{A}$ is a continuous degree of freedom and P is the generator of translations in the position quadrature. In between the two entangling operations the system is driven with the operator U_E. At the end, the ancillary system is measured in the X basis and the moments of the X variable can be estimated. The key of the method, as shown below, is that the distribution of results P(X) is a coarse-grained version of the full probability distribution of work.

To see how this method works we consider an initial thermal state for $\mathcal{S},$ i.e. ${{\rho }_{\beta }}={{{\rm e}}^{-\beta H}}/{{Z}_{\beta }}$ with ${{Z}_{\beta }}={\rm tr}\;{{{\rm e}}^{-\beta H}}$ the partition function. In turn, for $\mathcal{A}$ we consider a general state ${{\sigma }_{0}}$ (this is a generalization of the treatment presented in [21], where $\mathcal{A}$ was assumed to be a position eigenstate). The total state of the combined $\mathcal{S}-\mathcal{A}$ Universe can be obtained after the sequence of evolutions ${{U}_{T}}={{\tilde{U}}_{I}}{{U}_{E}}U_{I}^{\dagger }$. Now let us see the state after each step of the algorithm. Initially the state is ${{\sigma }_{0}}\;\otimes \;{{\rho }_{\beta }}$ and an entangling operation is applied, after this step the state can be written as

$$\begin{eqnarray*}&&U_{I}^{\dagger }({{\sigma }_{0}}\;\otimes \;{{\rho }_{\beta }}){{U}_{I}}=\mathop{\sum }\limits_{i}\frac{{{{\rm e}}^{-\beta {{E}_{i}}}}}{{{Z}_{\beta }}}{{{\rm e}}^{{\rm i}\kappa P{{E}_{i}}}}{{\sigma }_{0}}{{{\rm e}}^{-{\rm i}\kappa P{{E}_{i}}}}\;\otimes \;|{{\phi }_{i}}\rangle \langle {{\phi }_{i}}|.\end{eqnarray*}$$

Then, U_E is applied to system $\mathcal{S}$ and, after the last entangling operation, the final state ${{\rho }_{I}}={{U}_{T}}({{\sigma }_{0}}\;\otimes \;{{\rho }_{\beta }})U_{T}^{\dagger }$ can be written as

$$\begin{eqnarray*}&&{{\rho }_{I}}=\mathop{\sum }\limits_{i,m,n}\frac{{{{\rm e}}^{-\beta {{E}_{i}}}}}{{{Z}_{\beta }}}{{\chi }_{m,i}}\;\chi _{n,i}^{*}{{T}_{m,i}}{{\sigma }_{0}}T_{n,i}^{\dagger }\;\otimes \;|{{\tilde{\phi }}_{m}}\rangle \langle {{\tilde{\phi }}_{n}}|,\end{eqnarray*}$$

where the transition elements are ${{\chi }_{m,i}}=\langle {{\tilde{\phi }}_{m}}|{{U}_{E}}|{{\phi }_{i}}\rangle$. Moreover, the operators ${{T}_{m,i}}$ translate the state of $\mathcal{A}$ by an amount that depends on the energy difference $({{\tilde{E}}_{m}}-{{E}_{i}})$. Thus, they are defined as

$$\begin{eqnarray*}&&{{T}_{m,i}}={\rm exp} \left[ -i\kappa P({{\tilde{E}}_{m}}\ -{{E}_{i}}) \right].\end{eqnarray*}$$

From the total state ${{\rho }_{I}}$ we can compute the reduced density matrix of the auxiliary system ${{\rho }_{\mathcal{A}}}\equiv {\rm t}{{{\rm r}}_{\mathcal{S}}}{{\rho }_{I}}$. Thus,

$$\begin{eqnarray}&&{{\rho }_{\mathcal{A}}}=\mathop{\sum }\limits_{i,m}\frac{{{{\rm e}}^{-\beta {{E}_{i}}}}}{{{Z}_{\beta }}}|{{\chi }_{m,i}}{{|}^{2}}\ \ {{T}_{m,i}}\;{{\sigma }_{0}}\;T_{m,i}^{\dagger }.\end{eqnarray} \tag{ 1 }$$

From this expression we can compute the moments of the position variable of $\mathcal{A},$ defined as ${{\left\langle {{X}^{n}} \right\rangle }_{0}}={\rm tr}[{{X}^{n}}{{\rho }_{I}}]$. They turn out to be

$$\begin{eqnarray}\begin{array}{rcl} \left\langle {{X}^{n}} \right\rangle & = & \mathop{\sum }\limits_{i,m}\frac{{{{\rm e}}^{-\beta {{E}_{i}}}}}{{{Z}_{\beta }}}|{{\chi }_{m,i}}{{|}^{2}}{\rm tr}\left\{ {{\sigma }_{0}}{{\left[ X+\kappa ({{\tilde{E}}_{m}}-{{E}_{i}}) \right]}^{n}} \right\} \\ {} & = & \mathop{\sum }\limits_{k=0}^{n}\left( \begin{array}{ccccccccccccccc} n \\ k \\ \end{array} \right){{\kappa }^{n-k}}\;{\rm tr}({{\sigma }_{0}}{{X}^{k}})\left\langle {{W}^{n-k}} \right\rangle . \\ \end{array}\end{eqnarray} \tag{ 2 }$$

This equation establishes a simple relation between the moments of the position variable of $\mathcal{A}$ and those of the work distribution, which are defined as $\left\langle {{W}^{n}} \right\rangle ={{\sum }_{i,l}}\frac{{{{\rm e}}^{-\beta {{E}_{i}}}}}{{{Z}_{\beta }}}|{{\chi }_{l,i}}{{|}^{2}}{{({{\tilde{E}}_{l}}-{{E}_{i}})}^{n}}$. In particular, for the first two moments the equations are particularly simple. They read

$$\begin{eqnarray}&&\left\langle X \right\rangle ={{\left\langle X \right\rangle }_{0}}+\kappa \left\langle W \right\rangle ,\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\left\langle {{X}^{2}} \right\rangle ={{\left\langle {{X}^{2}} \right\rangle }_{0}}+2\kappa {{\left\langle X \right\rangle }_{0}}\left\langle W \right\rangle +{{\kappa }^{2}}\left\langle {{W}^{2}} \right\rangle ,\end{eqnarray} \tag{ 4 }$$

where ${{\langle \;\rangle }_{0}}$ denotes average on the initial state ${{\sigma }_{0}}$. These equations can be used to obtain simple relations between the dispersion (defined as $\Delta {{X}^{2}}={{\left\langle {{(X-{{\left\langle X \right\rangle }_{{}}})}^{2}} \right\rangle }_{{}}}$) and the skewness (defined as $\Delta {{X}^{3}}={{\left\langle {{(X-{{\left\langle X \right\rangle }_{{}}})}^{3}} \right\rangle }_{{}}}$) of the X coordinate, and those of the work distribution. Thus,

$$\begin{eqnarray}\begin{array}{rcl} \Delta {{X}^{2}} & = & \Delta X_{0}^{2}+{{\kappa }^{2}}\Delta {{W}^{2}}, \\ \Delta {{X}^{3}} & = & \Delta X_{0}^{3}+{{\kappa }^{3}}\Delta {{W}^{3}}. \\ \end{array}\end{eqnarray} \tag{ 5 }$$

This also shows that the scheme can also be used to test linear response results which relate the dissipated energy to the variance of the work distribution [1]. The above equations are worth analyzing: it is clear that the choice of the initial state ${{\sigma }_{0}}$ imposes strong constraints on the accuracy of the estimation of the properties of the work distribution. In fact, it is clear that in order to estimate $\Delta {{W}^{n}}$ by measuring $\Delta {{X}^{n}}$, it is better to choose initial states with small dispersions. The only states for which such dispersions vanish are the position eigenstates, which were considered in [21]. However, for a continuous variable system such as the one we are considering here, these states are unphysical. Instead, in this paper we will consider realistic scenarios for which the initial state is, typically, a coherent state (or a squeezed one). If instead of pure states we use mixed ones, it is obvious that we lose accuracy. In fact, if the initial state is thermal (for a harmonic oscillator with frequency ω) we have $\langle {{X}^{2}}\rangle \propto {\rm coth} \left[ \omega \beta /2 \right]$. Therefore, the precision of the estimate of work dispersion decreases with the temperature (or, equivalently, to achieve the same precision in the estimate of the work dispersion, we would need to measure the dispersion in X with much higher accuracy).

There is another generalization of the method presented in [21] that turns out to be useful for our purpose here. In fact, we will consider a more general interaction Hamiltonians between $\mathcal{S}$ and $\mathcal{A}$. As we will show, if the Hamiltonian is nonlinear in the momentum of $\mathcal{A}$ then the estimate of the moments of the work distribution may be simpler, and even more precise. To see this we consider an interaction Hamiltonian which induces an evolution operators given as ${{U}_{I}}={{{\rm e}}^{{\rm i}\kappa {{P}^{\alpha }}H}}$, for integer values of $\alpha$. In this case, it is simple to extend the previous results and to obtain an analytic expression for the moments of the work distribution. In fact, we find that

$$\begin{eqnarray}&&\Delta {{X}^{n}}=\Delta X_{0}^{n}+{{\kappa }^{n}}\Delta {{(WP_{0}^{\alpha -1})}^{n}},\end{eqnarray} \tag{ 6 }$$

a formula which is valid for $n=1,2,3$. A particularly simple case is attained for $\alpha =2$. Then, the second moment satisfy

$$\begin{eqnarray*}&&\left\langle {{X}^{2}} \right\rangle ={{\left\langle {{X}^{2}} \right\rangle }_{0}}+{{\kappa }^{2}}\left\langle {{W}^{2}} \right\rangle {{\left\langle {{P}^{2}} \right\rangle }_{0}}.\end{eqnarray*}$$

where we assumed ${{\langle X\rangle }_{0}}=0$ as is the case of a thermal symmetric state. This has an obvious interpretation: by considering an initial state which is squeezed in position we reduce $\Delta {{X}_{0}}$. Then, the estimate of $\Delta X$ (for fixed accuracy in the measurement of $\Delta X$) is higher than in the linear case. Again, all these results are independent of the initial state of the apparatus and will be useful in what follows.

In this section we start by explaining a scheme to reconstruct the probability distribution of the work done on or extracted from a cold atomic ensemble. The state of the ensemble, composed by N2-level atoms, can be described in terms of the collective angular momentum ${\boldsymbol{J}}$ that is the sum of the atomic spins. The components of the angular momentum operator fulfil the usual commutation relations (assuming throughout the paper that $\hbar =1$): $[{{J}_{x}},{{J}_{y}}]=i{{J}_{z}}$ and all the cyclic permutations. The ensemble is subject, as in previous experiments, to a magnetic field ${\boldsymbol{B}} (t)$ that can be continuously changed in time along any direction. The Hamiltonian governing the dynamics of the ensemble is therefore:

$$\begin{eqnarray}&&H(t)=-\gamma {\boldsymbol{B}} (t)\cdot {\boldsymbol{J}} ,\end{eqnarray} \tag{ 7 }$$

where γ is the gyromagnetic ratio and ${\boldsymbol{B}} (t)\equiv |{\boldsymbol{B}} |({{n}_{x}}(t),{{n}_{y}}(t),{{n}_{z}}(t))$ thus we are assuming that only the direction ${\boldsymbol{n}} (t)$ of the magnetic field and not its magnitude changes in time. The instantaneous eigenstates of H(t) coincide with those of the projection of ${\boldsymbol{J}}$ along the magnetic field direction ${\boldsymbol{n}} (t)$ and we label them as $|{{m}_{{\boldsymbol{n}} (t)}}\rangle$ with eigenvalue ${{E}_{m}}(t)=-\gamma |{\boldsymbol{B}} |{{m}_{{\boldsymbol{n}} (t)}}$.

We now compute the work done on the atomic ensemble, initially in the state $\rho (0)$, due to the variation of the magnetic field from ${\boldsymbol{B}} (0)$ to ${\boldsymbol{B}} (\tau )$ in a time τ. The ensemble state at any time can be calculated as $\rho (t)=U(t)\rho (0){{U}^{\dagger }}(t)$ where we have defined the unitary evolution operator which fulfills Schrödinger equation:

$$\begin{eqnarray}&&i\frac{\partial }{\partial t}U(t)=H(t)U(t)\end{eqnarray} \tag{ 8 }$$

with the initial condition $U(0)=\mathbb{I}$. We would like to stress here that we are not making any assumption on the time variation, slow or fast, of the direction of the magnetic field.

Taking as a definition the two time protocol, work W is a classical stochastic variable with probability distribution:

$$\begin{eqnarray}&&P(W)=\mathop{\sum }\limits_{m,m^{\prime} }{{p}_{m}}{{p}_{m^{\prime} |m}}\delta \left[ W-{{E}_{{{m}^{\prime }}}}(\tau )+{{E}_{m}}(0) \right],\end{eqnarray} \tag{ 9 }$$

where ${{p}_{m}}=\langle {{m}_{{\boldsymbol{n}} (0)}}|\rho (0)|{{m}_{{\boldsymbol{n}} (0)}}\rangle$ is the probability to find the initial state in the initial Hamiltonian eigenstate $|{{m}_{{\boldsymbol{n}} (0)}}\rangle$ and

$$\begin{eqnarray*}&&{{p}_{m^{\prime} |m}}={{\left| \ \left\langle m_{{\boldsymbol{n}} (\tau )}^{\prime }\left| U(\tau ) \right|{{m}_{{\boldsymbol{n}} (0)}} \right\rangle \ \right|}^{2}}\end{eqnarray*}$$

is the conditional probability that evolving with the evolution operator $U(\tau )$ the initial Hamiltonian eigenstate $|{{m}_{{\boldsymbol{n}} (0)}}\rangle$ the state of the system is found, at time τ, in the final Hamiltonian eigenstate $|m_{{\boldsymbol{n}} (\tau )}^{\prime }\rangle$.

We start with a simple case where we assume that the initial magnetic field is pointing along the z direction and, at t = 0, is instantaneously rotated to the y axis, thus $U(\tau )=\mathbb{I}$. We assume that the ensemble is initially in thermal equilibrium with inverse temperature β (assuming the Boltzmann constant k_B = 1) so that its state is:

$$\begin{eqnarray}&&\rho (0)=\frac{1}{Z(0)}{{{\rm e}}^{-\beta H(0)}},\end{eqnarray} \tag{ 10 }$$

where $Z(0)={\rm tr}\;{{{\rm e}}^{-\beta H(0)}}$ is the initial partition function ensuring the normalization of the state density matrix. In this case the work probability distribution (WPD) depends on the overlaps $|\langle m_{y}^{\prime }|{{m}_{z}}\rangle {{|}^{2}}$ between the angular momentum eigenstates along the z and y directions. These can be calculated in terms of the Wigner D-matrix but the result is cumbersome and will not be reported here. The results for the WPD can be found in figure 2. It can be observed that for very low temperatures the probability distribution resembles a Gaussian function. This can be explained as follows. The initial state is polarized along the z direction, so each spin is in a superposition of the up and down states along the y axis. As the total state is the tensor product of each spin wavefunction, the resulting distribution is binomial, thus approaching a Gaussian shape for large number of atoms. More precisely, for a large number of particles, and using Holstein–Primakoff approximation, the atomic state can be regarded as a coherent state. For the instantaneous quench we are considering, the WPD depends only on the transition probabilities ${{p}_{m^{\prime} |m}}$ which, in the Holstein–Primakoff picture, represents the wave function squared of such coherent state, therefore a Gaussian function, its position-like operator being proportional to the angular momentum J_y along the final magnetic field.

Zoom In Zoom Out Reset image size

For large temperatures this is not true anymore, and other transitions from initial excited states acquire a higher weight. These give rise to many more peaks distorting the WPD to a skewed function. We have studied the normalized skewness ${\rm Skew}[{\rm W}]$ of the WPD as a function of temperature. The normalized skewness is defined as:

$$\begin{eqnarray}&&{\rm Skew}[{\rm W}]=\left\langle {{\left( \frac{{\rm W}-\langle {\rm W}\rangle }{{{\sigma }_{W}}} \right)}^{3}} \right\rangle ,\end{eqnarray} \tag{ 11 }$$

where ${{\sigma }_{W}}$ is the work standard deviation.

The results, reported in figure 3, show that the skewness is always negative meaning that, although most of the probability is located to the right of the maximum of the distribution, there is a long tail of small probabilities to the left of the maximum. This is not uncommon for the WPD [19, 28] and sometimes it gives rise to non-zero probability for negative work values. Figure 3 shows also an interesting result: the skewness approaches zero for very small, as we said earlier, or very large temperatures. In the large temperature limit the skewness also approaches zero because the initial state is proportional to the identity meaning that all energy eigenstates are equally probable. This leads to a symmetric distribution as the transitions probabilities are symmetric: $|\langle {{m}_{z}}|{{m}_{y}}\rangle {{|}^{2}}=|\langle -{{m}_{z}}|{{m}_{y}}\rangle {{|}^{2}}=|\langle {{m}_{z}}|-{{m}_{y}}\rangle {{|}^{2}}=|\langle -{{m}_{z}}|-{{m}_{y}}\rangle {{|}^{2}}$. Thus P(W) is symmetric around zero and the skewness reduces to zero. For intermediate temperatures ${{(\beta \gamma |B|)}^{-1}}\sim 5$ there is a maximum of the absolute value of the skewness.

Zoom In Zoom Out Reset image size

We now consider a slow quench of the magnetic field and calculate the work done on the ensemble for different speeds $\omega =1/\tau$. We therefore assume that the magnetic field rotates at constant angular speed ω as:

$$\begin{eqnarray}&&{\boldsymbol{B}} (t)=|{\boldsymbol{B}} |({\rm cos} \omega t{\boldsymbol{\hat{k}}} +{\rm sin} \omega t{\boldsymbol{\hat{j}}} ).\end{eqnarray} \tag{ 12 }$$

For this particular choice the eigenenergies E_m do not depend on time and there are no degeneracies. We thus expect that for sufficiently small angular speed ω the evolution to be (quantum) adiabatic: since there are no transitions induced by the time variation of the Hamiltonian, the state populations do not change in time and the state at all times remains in thermal equilibrium. In this regime we expect the average work $\langle W\rangle$ to approach the free energy difference $\Delta F$ which, for the process we consider, is null. For higher speed ω we expect the process to excite the system and bring it out of equilibrium. This in turn produces irreversible work defined as:

$$\begin{eqnarray}&&{{W}_{{\rm irr}}}=\langle W\rangle -\Delta F=\langle W\rangle ,\end{eqnarray} \tag{ 13 }$$

where the last equality follows from our assumptions that the modulus of the magnetic field does not change.

The results for $\langle W\rangle$ are shown in figure 4. As we expected, for very small ω the average work tends to zero while growing and approaching a limiting value for very fast quenches. This value coincides with the average work calculated assuming instantaneous quenches $U(\tau )=\mathbb{I}$. The figure also shows the dependence of the average work for different temperatures. For high temperatures the average work reduces as the system initially occupies many excited states. In the limit of infinite temperature, the initial state of the system is the unitary invariant completely mixed state proportional to the identity. In this limit, the average work is zero because any transformation leaves the state unaltered.

Zoom In Zoom Out Reset image size

4.1. The scheme

We propose a scheme, following [21], to experimentally reconstruct the probability distribution of the work done on an atomic ensemble when varying the applied magnetic field. To this end we use a light-matter interface based on the Faraday rotation [33]. If light polarized along the x axis propagates along the YZ plane and illuminates the atomic ensemble at an angle α with the z axis, the interaction Hamiltonian reads:

$$\begin{eqnarray}&&{{H}_{I}}(\alpha )=\frac{a}{T}{{S}_{z}}({\rm cos} \alpha {{J}_{z}}+{\rm sin} \alpha {{J}_{y}}),\end{eqnarray} \tag{ 14 }$$

where a is the coupling constant and T is the duration of the pulse. The Stokes operators are defined as:

$$\begin{eqnarray}&&{{S}_{x}}=\frac{1}{2}(a_{x}^{\dagger }{{a}_{x}}-a_{y}^{\dagger }{{a}_{y}}),\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{{S}_{y}}=\frac{1}{2}(a_{y}^{\dagger }{{a}_{x}}+a_{x}^{\dagger }{{a}_{y}}),\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{{S}_{z}}=\frac{1}{2i}(a_{y}^{\dagger }{{a}_{x}}-a_{x}^{\dagger }{{a}_{y}}),\end{eqnarray} \tag{ 17 }$$

where the operators a_x and a_y annihilates a photon with polarization along x and y, respectively. We assume that the light pulse is strongly polarized along the x axis: ${{S}_{x}}\approx \langle {{S}_{x}}\rangle ={{N}_{{\rm ph}}}/2$ where N_ph is the number of photons. Within this approximation, we can treat the Stokes operators in the two perpendicular directions as conjugated variables: ${{S}_{z}}=P\sqrt{{{N}_{{\rm ph}}}/2}$ and ${{S}_{y}}=X\sqrt{{{N}_{{\rm ph}}}/2}$, so that $[X,P]=i$.

Using these assumptions the evolution operator corresponding to a pulse with Hamiltonian (14) is:

$$\begin{eqnarray}&&{{U}_{I}}(\alpha )={\rm exp} \left[ -i\kappa PJ(\alpha ) \right],\end{eqnarray} \tag{ 18 }$$

where $\kappa =a\sqrt{{{N}_{{\rm ph}}}/2}$ and $J(\alpha )=({\rm cos} \alpha {{J}_{z}}+{\rm sin} \alpha {{J}_{y}})$. With atomic ensemble at room temperatures the coefficient κ could be very small for our purposes, as we would need a value $\kappa \approx 1$. For ultracold atoms the optical depth, and therefore κ, could be made larger although results in this direction have not yet been demonstrated. We could also write the transformation ${{U}_{I}}(\alpha )$ as:

$$\begin{eqnarray}&&{{U}_{I}}(\alpha )={\rm exp} \left[ i\tilde{\kappa }P{{H}_{A}}(\alpha ) \right],\end{eqnarray} \tag{ 19 }$$

where ${{H}_{A}}(\alpha )=\gamma |{\boldsymbol{B}} |J(\alpha )$, which is equivalent to H(t) in equation (7), and we set $\tilde{\kappa }=\kappa /(\gamma |{\boldsymbol{B}} |)$. Thus it is clear that transformation ${{U}_{I}}(\alpha )$ is a spatial translation of the continuous state of light conditional on the atomic ensemble energy. It is this conditional interaction that makes it possible to read the WPD from the state of the light.

We now follow the idea from [21]. Initially the polarization fluctuation state of the light is assumed to be characterized by a Gaussian wave function, $|X=0{{\rangle }_{L}}$ centred in zero with variance ${{\sigma }^{2}}$. Although the vacuum would correspond to ${{\sigma }^{2}}=1/2$ we carry on our analysis for generic σ thus encompassing also squeezed states. For the sake of simplicity we assume that the atoms are initially in a pure state $|\psi {{\rangle }_{A}}$, but the same exact scheme works also for mixed states.

As illustrated in figure 5, the protocol consists in shining the atoms with a laser beam, strongly polarized along x and propagating along a direction on the yz plane and forming an angle ${{\alpha }_{0}}+\pi$ with the z axis. During this first step, light and atoms interact with a Hamiltonian proportional to ${{H}_{A}}({{\alpha }_{0}})$. While the beam is stored in a quantum memory, the atoms undergo the process during which the magnetic field is rotated eventually pointing to the direction in the yz plane forming an angle ${{\alpha }_{1}}$ with the z axis. The atomic state is evolved with evolution operator U(t) fulfilling equation (8). Finally, the light beam is retrieved from the quantum memory and let pass through the atoms along a direction forming an angle ${{\alpha }_{1}}$ with the z axis. During this step light and atoms interact with a Hamiltonian proportional to ${{H}_{A}}({{\alpha }_{1}})$. Thus, at the end the state of the light encodes the difference between the final and initial energy for each posible quantum trajectory. It is at the very end, when the measurement is performed, that coherence is destroyed. In this way the method samples W with probability P(W).

Zoom In Zoom Out Reset image size

Mathematically the state of atoms and light before light is measured is

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} |\Psi {{\rangle }_{AL}} & = & {{U}_{I}}({{\alpha }_{1}})U(t)U_{I}^{\dagger }({{\alpha }_{0}})|\psi {{\rangle }_{A}}|X=0{{\rangle }_{L}} \\ {} & = & \mathop{\sum }\limits_{m,m^{\prime} }{{c}_{m}}{{c}_{m^{\prime} |m}}\ \ |{{E}_{{{m}^{\prime }}}}{{\rangle }_{A}}|X=\kappa ({{E}_{{{m}^{\prime }}}}-{{E}_{m}}){{\rangle }_{L}} \\ \end{array},\end{eqnarray} \tag{ 20 }$$

where ${{c}_{m}}{{=}_{A}}{{\langle {{E}_{m}}|\psi \rangle }_{A}}$ and ${{c}_{m^{\prime} |m}}{{=}_{A}}{{\langle {{E}_{{{m}^{\prime }}}}|U(t)|{{E}_{m}}\rangle }_{A}}$ and where states like $|X=\kappa ({{E}_{{{m}^{\prime }}}}-{{E}_{m}}){{\rangle }_{L}}$ represent the initial state of light rigidly translated by the quantity $\kappa ({{E}_{{{m}^{\prime }}}}-{{E}_{m}})$.

The reconstructed work distribution can be found from the probability density distribution of the X quadrature of light (assuming no degeneracies):

$$\begin{eqnarray}&&{{P}_{L}}(X)=\mathop{\sum }\limits_{m,m^{\prime} }{{p}_{m}}{{p}_{m^{\prime} |m}}\frac{1}{\sqrt{2\pi {{\sigma }^{2}}}}{\rm exp} \left[ -\frac{{{(X-\kappa ({{E}_{{{m}^{\prime }}}}-{{E}_{m}}))}^{2}}}{2{{\sigma }^{2}}} \right],\end{eqnarray} \tag{ 21 }$$

where we have identified ${{p}_{m}}=|{{c}_{m}}{{|}^{2}}$ and ${{p}_{m^{\prime} |m}}=|{{c}_{m^{\prime} |m}}{{|}^{2}}$. Notice the difference with the work distribution in equation (9): apart from the conversion factor κ the light distribution corresponds to a coarse grained version of P(W) where Dirac delta functions have been replaced by Gaussians with width σ. We therefore expect a faithful reconstruction of the WPD when $\sigma /\kappa$ is sufficiently smaller than the energy change ${{E}_{{{m}^{\prime }}}}-{{E}_{m}}$.

Using equations (3) and (4) we can obtain the first two moments of the light distribution:

$$\begin{eqnarray}&&\langle X\rangle =\kappa \langle W\rangle .\end{eqnarray} \tag{ 22 }$$

And for the second moment:

$$\begin{eqnarray}&&\langle {{X}^{2}}\rangle ={{\sigma }^{2}}+{{\kappa }^{2}}\langle {{W}^{2}}\rangle \end{eqnarray} \tag{ 23 }$$

so that the variance of the light distribution is:

$$\begin{eqnarray}&&\Delta {{X}^{2}}={{\sigma }^{2}}+{{\kappa }^{2}}\Delta {{W}^{2}}\end{eqnarray} \tag{ 24 }$$

Therefore provided that κ is sufficiently strong we can estimate the first two moments of the work distribution by measuring the light fluctuations. A similar two- or multiple-passage protocol has been previously discussed in [34] for the implementation of a quantum memory.

4.2. An example

To showcase our proposal we consider the process described in section 3. The atoms are initially in thermal equilibrium and subject to a magnetic field along the z direction. The magnetic field is suddenly rotated to the y direction and we want to reconstruct the WPD of this process and compare it with the exact one calculated in section 3.

The distribution for the light quadrature X can be found by inserting in equation (21) the expressions for p_m and ${{p}_{m^{\prime} |m}}$ used in section 3. The light probability distribution is shown in figure 6 for zero temperature and for a temperature ${{\beta }^{-1}}=1$. The light distribution is the sum of narrow Gaussians at each of the red points of figure 2. Therefore it represents a coarse grained version of it. Nevertheless, all the important features such as the first few moments and the overall shape agree with the exact result. So even with modest resources like using coherent states ${{\sigma }^{2}}=1/2$ and a coupling $\kappa =2$ it is possible to reconstruct quite faithfully the work probability distribution.

Zoom In Zoom Out Reset image size

The ultimate test of our reconstructed WPD is Jarzynski equality [1]

$$\begin{eqnarray}&&\langle {{{\rm e}}^{-\beta W}}\rangle ={{{\rm e}}^{-\beta \Delta F}}=1,\end{eqnarray} \tag{ 25 }$$

where the last equality follows from the fact that for us the free energy difference $\Delta F$ is zero.

Since the work variable W corresponds to the renormalized quadrature $X/\kappa$ we compute:

$$\begin{eqnarray}\begin{array}{rcl} \langle {{{\rm e}}^{-\beta X/\kappa }}\rangle & = & \int _{-\infty }^{\infty }{{{\rm e}}^{-\beta X/\kappa }}{{P}_{L}}(X){\rm d}X \\ {} & = & \mathop{\sum }\limits_{m,m^{\prime} }{{p}_{m}}{{p}_{m^{\prime} |m}}{{{\rm e}}^{-\beta ({{E}_{{{m}^{\prime }}}}-{{E}_{m}})}}{\rm exp} \left[ \frac{{{\sigma }^{2}}{{\beta }^{2}}}{2{{\kappa }^{2}}} \right] \\ {} & = & \langle {{{\rm e}}^{-\beta W}}\rangle {\rm exp} \left[ \frac{{{\sigma }^{2}}{{\beta }^{2}}}{2{{\kappa }^{2}}} \right]={\rm exp} \left[ \frac{{{\sigma }^{2}}{{\beta }^{2}}}{2{{\kappa }^{2}}} \right], \\ \end{array}\end{eqnarray} \tag{ 26 }$$

where in the last equality we used Jarzynski relation (25). Thus, Jarskynski equality is estimated with a correction that decreases with the coupling κ and the temperature and decreases with the width σ of the initial light polarization state. A similar results was found for generalized energy measurements [35].

Using the following parameters: $N=40;\kappa =2;{{\sigma }^{2}}=1/2;\beta =1$, we obtain

$$\begin{eqnarray}&&\langle {{{\rm e}}^{-\beta X/\kappa }}\rangle =1.06\end{eqnarray} \tag{ 27 }$$

which is only 6 $\%$ from the expected result. A plot of the correcting factor ${\rm exp} \left[ \frac{{{\sigma }^{2}}{{\beta }^{2}}}{2{{\kappa }^{2}}} \right]$ is shown in figure 7 where it is clear that values of κ above 5 gives a negligible correction to the Jarzynski equality.

Zoom In Zoom Out Reset image size

5.1. Generalities on fluctuation relations in open quantum system

So far we have discussed a method to reconstruct the probability distribution of work done or extracted from an isolated system. We now extend the method to a non-unitary evolution in which the system S is coupled to an environment E during the process. Relaxing the assumptions of unitary processes, we have to be careful when talking about work. The system in fact exchanges energy, which we may well call heat, with the environment. Thus it is more accurate to talk about the energy change $\mathcal{E}$ of the system and its fluctuation relations [36]. We are assuming as in the previous section that we perform a two-time energy measurement on the system, the only difference is that now the system evolution is not unitary.

In the open quantum system scenario, it is common to introduce a complete positive trace preserving (CPTP) map Φ acting on the initial density matrix ${{\rho }_{S}}(0)$ (assuming no initial system-environment correlations) to obtain the evolved density matrix at time t. If the evolution of the combined system plus environment is unitary and governed by the operator U_SE the map can be expressed as:

$$\begin{eqnarray}&&{{\rho }_{S}}(t)=\Phi [{{\rho }_{S}}(0)]={\rm T}{{{\rm r}}_{E}}\left[ {{U}_{SE}}{{\rho }_{S}}(0)\;\otimes \;{{\rho }_{E}}(0)U_{SE}^{\dagger } \right].\end{eqnarray} \tag{ 28 }$$

The map can be conveniently cast in terms of Kraus operators:

$$\begin{eqnarray}&&\Phi [{{\rho }_{S}}(0)]=\mathop{\sum }\limits_{k}{{A}_{k}}{{\rho }_{S}}(0)A_{k}^{\dagger },\end{eqnarray} \tag{ 29 }$$

where, due to the trace preserving nature of the map, the Kraus operators fulfil ${{\sum }_{k}}A_{k}^{\dagger }{{A}_{k}}=\mathbb{I}$. It is possible to define a dual map ${{\Phi }^{*}}$ as

$$\begin{eqnarray}&&{{\Phi }^{*}}[{{\rho }_{S}}(0)]=\mathop{\sum }\limits_{k}A_{k}^{\dagger }{{\rho }_{S}}(0){{A}_{k}}\end{eqnarray} \tag{ 30 }$$

which however is not in general trace preserving. A map Φ is called unital if the corresponding dual map ${{\Phi }^{*}}$ is trace preserving. This condition is equivalent to requiring that Φ maps the completely mixed state ${{\mathbb{I}}_{S}}$ into itself: $\Phi ({{\mathbb{I}}_{S}})={{\mathbb{I}}_{S}}$.

It has been shown before [12, 13] that when calculating an analogous relation to Jarzynki's one obtains a result that depends on the dual map:

$$\begin{eqnarray}&&\langle {{{\rm e}}^{-\beta \mathcal{E}}}\rangle ={\rm Tr}\left[ {{\rho }_{S}}(0)\Phi (\mathbb{I}) \right],\end{eqnarray} \tag{ 31 }$$

where ${{\rho }_{S}}(0)=\frac{{{{\rm e}}^{-\beta {{H}_{S}}}}}{{\rm tr}\;{{{\rm e}}^{-\beta {{H}_{S}}}}}$ and the quantity on the right-hand side has been called efficacy of the process. Thus, Jarzynski equality for the energy change is fulfilled, i.e. the right hand side is 1 as we are considering zero free energy change, if and only if the map Φ is unital.

5.2. An example with atomic spins in optical lattices

To test the ideas discussed in the previous paragraphs, we consider the setup sketched in figure 8. A superlattice potential of double wells is created with the aid of two standing waves with wave vectors having a ratio of 2. For large enough intensities, and assuming no vacancies, each well will contain exactly one atom, i.e. the system is in a Mott insulator with unit filling. Probing ultracold atoms in superlattice potentials has been proposed in [37]. We assume the atom sitting in the left well to be the system and the atom in the right well to be the environment. In this limit tunnelling is suppressed and a super-exchange interaction between the pseudo-spin internal levels can be induced by lowering the barrier between the two wells. The spins are initially in thermal equilibrium at the same temperature:

$$\begin{eqnarray}&&{{\rho }_{R}}(0)=\frac{{{{\rm e}}^{-\beta {{H}_{R}}}}}{{\rm Tr}{{{\rm e}}^{-\beta {{H}_{R}}}}},\quad R=S,E\end{eqnarray} \tag{ 32 }$$

with ${{H}_{R}}=B\sigma _{x}^{R}$ and $S,E$ indicates the system and environment spins, respectively.

Zoom In Zoom Out Reset image size

To measure energy change $\mathcal{E}$, we first measure the initial energy of the system by projecting the initial density matrix on the eigenstates $|\rangle +$ and $|\rangle -$ of H_S. Then the thermodynamic process consists in coupling system and environment with the XXZ interaction:

$$\begin{eqnarray}&&{{H}_{SE}}=\sigma _{x}^{S}\sigma _{x}^{E}+\sigma _{y}^{S}\sigma _{y}^{E}+\Delta \sigma _{z}^{S}\sigma _{z}^{E}\end{eqnarray} \tag{ 33 }$$

and evolving it in time for a time t with the evolution operator ${{U}_{SE}}={\rm exp} [-i{{H}_{SE}}t]$. In equation (33), Δ is the interaction anisotropy which can be tuned by accurately changing the atoms scattering length near a Feshbach resonance [31].

We then consider the reduced density matrix of the system only and measure again the energy H_S. The probability distribution of the energy change $\mathcal{E}$ is:

$$\begin{eqnarray}\begin{array}{rcl} P(\mathcal{E}) & = & \langle +|{{\rho }_{S}}(0)|+\rangle \ \langle +|{{\rho }_{+}}(t)|+\rangle \delta (\mathcal{E}) \\ {} & {} & +\langle +|{{\rho }_{S}}(0)|+\rangle \ \langle -|{{\rho }_{+}}(t)|-\rangle \delta (\mathcal{E}-2B) \\ {} & {} & +\langle -|{{\rho }_{S}}(0)|-\rangle \ \langle +|{{\rho }_{-}}(t)|+\rangle \delta (\mathcal{E}+2B) \\ {} & {} & +\langle -|{{\rho }_{S}}(0)|-\rangle \ \langle -|{{\rho }_{-}}(t)|-\rangle \delta (\mathcal{E}), \\ \end{array}\end{eqnarray} \tag{ 34 }$$

where

$$\begin{eqnarray}&&{{\rho }_{\pm }}(t)={\rm T}{{{\rm r}}_{E}}\left[ {{{\rm e}}^{-i{{H}_{SE}}t}}|\pm \rangle \langle \pm |\;\otimes \;{{\rho }_{E}}(0){{{\rm e}}^{i{{H}_{SE}}t}} \right].\end{eqnarray} \tag{ 35 }$$

It is easy to check that the probability distribution equation (34) is normalized. From equation (34) is it possible to write an analogous of Jarzinsky equality:

$$\begin{eqnarray}&&\langle {{{\rm e}}^{-\beta \mathcal{E}}}\rangle =1+{{{\rm tanh} }^{2}}(\beta B){\rm sin} (2t){\rm sin} (2\Delta t).\end{eqnarray} \tag{ 36 }$$

Thus, not surprisingly we get a time-dependent correction to Jarzinsky equality due to the openness of the system evolution.

However, notice that for $\Delta =0$, corresponding to the XXZ model, the Jarzynski equality is fulfilled at all times. This is because the CPTP map Φ that evolves the system becomes unital for $\Delta =0$. In fact, applying the map to the completely mixed state, we obtain:

$$\begin{eqnarray}\Phi \left[ \frac{{{\mathbb{I}}_{S}}}{2} \right]=\frac{1}{2}\left( \begin{array}{ccccccccccccccc} 1 & g(t) \\ g(t) & 1 \\ \end{array} \right),\end{eqnarray} \tag{ 37 }$$

where $g(t)={\rm tanh} (\beta B){\rm sin} (2t){\rm sin} (2\Delta t)$. Thus the 1-norm of the difference of $\Phi [{{\mathbb{I}}_{S}}/2]$ from the identity is related to the violation of the Jarzynski equality:

$$\begin{eqnarray}&&\left\|\Phi \left[ \frac{{{\mathbb{I}}_{S}}}{2} \right]-\frac{{{\mathbb{I}}_{S}}}{2}\right\|{_{1}}=\frac{1}{2}|g(t)|.\end{eqnarray} \tag{ 38 }$$

5.3. Scheme to reconstruct energy change in optical lattices

We finally consider the reconstruction of the dissipated energy distribution with the Faraday rotation scheme. As shown in figure 8 the scheme is very similar to the one with an atomic ensemble. This time, the pulse produces a standing wave with a double period with respect to the optical lattice. This means that only the left-most spin in each double well is strongly illuminated by the light probe. In this way we can measure the total energy of all the identical system spins in the lattice.

The pulse is first sent through the atomic array and then stored in a quantum memory, as before. In this first stage the light polarization fluctuation contains information about the initial energy of the system. Then, while the light is stored, the atoms interact according to the XXZ interaction described before. After this, the light pulse is retrieved from the memory, its polarization is rotated by 180 degrees by a half-wave plate, and let pass through the atoms again, thus reading the final energy of the system spins. The pulse is finally analysed with a homodyne detection.

As in section 4, the reconstructed distribution is obtained from equation (34) with the substitution:

$$\begin{eqnarray}&&\delta (\mathcal{E})\to \frac{1}{\sqrt{2\pi }\sigma }{\rm exp} \left[ -\frac{1}{2{{\sigma }^{2}}}{{(X-\kappa \mathcal{E})}^{2}} \right].\end{eqnarray} \tag{ 39 }$$

Thus the probability distribution of the light quadrature X is a coarse-grained version of the true energy change probability distribution. The average exponentiated energy for a single spin is corrected by a factor:

$$\begin{eqnarray}&&\langle {{{\rm e}}^{-\beta X/\kappa }}\rangle ={{{\rm e}}^{\frac{{{\sigma }^{2}}{{\beta }^{2}}}{2{{\kappa }^{2}}}}}\langle {{{\rm e}}^{-\beta \mathcal{E}}}\rangle \end{eqnarray} \tag{ 40 }$$

therefore if $\sigma \ll \kappa$ the reconstruction is possible. Notice that when the map is unital ($\Delta =0$) the reconstructed Jarzynski quantity is time independent. Therefore even if the reconstructed result differs from the correct one, from its time dependence, it unambiguously signals the unitality of the map.

So far, we have calculated the dissipated energy distribution for a pair of spins. As the light interacts with all the N pairs of atoms in the lattice, the total dissipated energy is the sum of all the energies of each system atom. As these behave independently the joint probability distribution is factorized, so that the expectation value of the exponential becomes the Nth power of the results in (36) and (40).

In summary, we have proposed an experimentally feasible method to reconstruct the full distribution of the energy change, specifically work and heat, of ultracold atomic gases. Although our proposal employs a light-matter interface based on the quantum Faraday rotation, we stress that it could be adapted to other similar setups, for example a Bose–Einstein condensate in a cavity. Finally, our proposal is able to reveal fundamental properties of non-unitary evolutions that can be exploited for quantum environment engineering.

We thank S Deffner, J Goold, J Sherson for fruitful discussions. G D C acknowledges hospitality at Universidad de Buenos Aires where part of this work was carried on. G D C acknowledges support by the UK EPSRC (EP/L005026/1), the John Templeton Foundation (grant ID 43467), the EU Collaborative Project TherMiQ (Grant Agreement 618074) and the COST Action MP1209. A J R and J P P acknowledge support from ANPCyT (PICT-2010–02483, PICT-2013–0621), CONICET and UBACyT.