Experimental test of fluctuation relations for driven open quantum systems with an NV center

Quantum fluctuation relations [1–4] are exact relations in quantum statistical mechanics that allow for a compact and effective description of a plethora of non-equilibrium phenomena such as the statistical nature of the second law of thermodynamics and the emergence of an arrow of time [5], Onsager relations, Green–Kubo relations and high order relations between responses and multi-point correlation functions [6, 7], the microscopic mechanisms at the basis of Clausius inequality [8], and the Carnot theorem imposing the celebrated Carnot efficiency as the universal upper bound to efficiency of heat engines [9, 10]. Recently, they have been shown to be at the basis of thermodynamic uncertainty relations expressing lower bounds on entropy production in terms of the precision (signal-to-noise ratio) of fluctuating non-equilibrium charges (i.e., integrated non-equilibrium currents) [11, 12], or the asymmetry of their statistics [13].

Generally speaking, quantum fluctuation relations refer to a scenario where a quantum system, prepared in thermal equilibrium at some inverse temperature β₁, is then placed in contact (sequentially or at the same time) with several thermal baths at various temperatures β₂, ..., β_N while being possibly acted upon by a work source resulting in a time-dependent forcing [8]. In integral form, the quantum fluctuation relation reads [4]:

$$\begin{equation}\langle {\text{e}}^{-{\sum }_{i}{\beta }_{i}{\Delta}{E}_{i}}\rangle ={\text{e}}^{-{\sum }_{i}{\beta }_{i}{\Delta}{F}_{i}}\text{,}\end{equation} \tag{ 1 }$$

where ΔF_i is the free energy change between the (hypothetical) equilibrium state of component i of the global system + bath compound at the protocol's end time, and its initial free energy ⁷ , while ΔE_i is the according stochastic energy change as obtained by taking the difference of projectively measured eigen-energies of component i at the beginning and end of the driving protocol (here the so called two-point measurement scheme is implied [4]). The symbol ⟨⋅⟩ denotes average over the statistics of the joint probability distribution p(ΔE₁, ..., ΔE_N ) for the energy change of each component.

The above relation encompasses a number of typical situations encountered in non-equilibrium thermodynamics. For example, a system initially in equilibrium and in contact with a thermal bath at some inverse temperature β is acted upon by a time dependent force: in this case β₁ = β₂ = β, and since the work done by the external source is given by W = ΔE₁ + ΔE₂, the Jarzynski relation follows ⟨e^−βW ⟩ = e^−βΔF , with ΔF the global free energy variation [14]. This relation has been verified experimentally with closed quantum systems (namely in absence of a thermal bath) [15, 16], and for a quantum system in contact with a dephasing environment that does not modify the system energy [17]. However, the verification of this relation with an open quantum system that exchanges energy with its environment is still elusive, due to the fact that the measurement of total energy injected in system and bath is typically very challenging. Possible solutions to this challenge are offered by the interferometric method for the measurement of work statistics [18, 19] as discussed in [20], or with reference to superconducting circuit setups, by calorimetric measurements of energy exchanges between a qubit and a resistor [21].

A complementary situation is when two quantum systems prepared initially at different temperatures, are placed in thermal contact for some time τ. Provided the interaction energy is negligible, no work is associated to this process, and the general relation takes the form of the so called exchange fluctuation relation ⟨e^ΔβΔE ⟩ = 1 [22], where Δβ = β₂ − β₁ and ΔE = ΔE₁ = −ΔE₂. This relation has been recently verified in the quantum regime with NMR setup [23] and with nitrogen vacancy (NV) centers in diamond [24].

A situation of central interest in quantum thermodynamics combines the previous two scenarios, namely when not only a system exchanges energy with its surroundings because of a thermal gradient, but also with an external work source. That is indeed what happens during the operation of a heat engine. The experimental verification of equation (1) when there is a simultaneous exchange of heat and work poses significant challenges, as typical quantum platforms allow to experimentally access the system energy change ΔE = W − Q, but not the individual terms W and Q.

In this work, we address this challenge and experimentally explore the validity the quantum fluctuation relation for a driven-dissipative dynamics (equation (1)) realized by the ground state of a single negatively-charged NV center in diamond. The NV center is well established as a biocompatible quantum sensor [25–27], single-photon source [28], as well as building block for quantum registers [29–31], thanks to the coherent control of its ground state spin and optical addressability. The capability of coherently manipulating a two or three level spin system, and to engineer dissipative channels to mimic the thermal contact with a quantum reservoir, make the NV center also an excellent platform to explore quantum thermodynamics processes [24, 32]. Here, we use an NV center spin qubit to prove the validity of the fluctuation relation in equation (1) in two special cases of driven-dissipative dynamics, where the measurement of the spin energy change ΔE suffices: (i) when the environment with which the system exchanges energy has infinite temperature; (ii) when the duration of the process is chosen in such a way that the total work done on the system amounts to zero, although generally the power delivered to it is non-null.

The results obtained in this work expand the experimental basis over which quantum fluctuation relations have been tested, and prepare the ground for a fully fledged testing of equation (1) in the general case when energy is exchanged both in the form of heat and work.

To engineer a driven-dissipative dynamics, we use the electronic spin ground state of a single NV center in diamond at room temperature. An external magnetic field, B = 38.9 mT, is aligned along the NV spin quantization axis to remove the degeneracy between the spin projections m_S = ±1. The two states $\left\vert {m}_{\text{S}}=+1\right\rangle \equiv \left\vert 1\right\rangle$ and $\left\vert {m}_{\text{S}}=0\right\rangle \equiv \left\vert 0\right\rangle$ form the basis of a two-level system. The two-level system is optically initialized into $\left\vert 0\right\rangle$ before each experiment, and its state is read out by measuring the NV spin-selective photoluminescence intensity [33, 34].

We apply a continuous resonant microwave (mw) driving field, described by a Hamiltonian that, in the mw rotating frame, reads

$$\begin{equation}\mathcal{H}\left(t\right)=\frac{\hslash \omega \left(t\right)}{2}\left({\sigma }_{x}\enspace \mathrm{cos}\enspace \theta t+{\sigma }_{y}\enspace \mathrm{sin}\enspace \theta t\right),\end{equation} \tag{ 2 }$$

where σ_i are Pauli matrices, ω(t) is the Rabi frequency that is set by the mw amplitude, and θ is an angular frequency of the mw phase variation. We denote the instantaneous eigenstates of the Hamiltonian in equation (2) as ${\rho }_{{\pm}}\left(t\right)\equiv \left\vert {E}_{{\pm}}\left(t\right)\right\rangle \left\langle {E}_{{\pm}}\left(t\right)\right\vert$, and their corresponding eigenvalues as E_±(t) = ±ℏω(t)/2.

Besides the unitary evolution, the system interacts with a train of N_L temporally-equidistant short laser pulses, as depicted in figure 1. In a previous work [24], we have demonstrated that a photon absorption from each laser pulse results in the application of a quantum projective measurement and the opening of a dissipative channel. The application of a train of short laser pulses brings the spin into an asymptotic Gibbs state with inverse (pseudo)temperature β_R [24] and, as such, is analogous to the dissipation channel caused by the interaction with a heat reservoir R. The value of β_R is determined by a combination of factors, which are the relative orientation of the driving Hamiltonian operator with respect to the NV quantization axis (σ_z ), the Rabi frequency of the driving field, the laser absorption probability p_a, and the time between pulses [24].

Zoom In Zoom Out Reset image size

To measure the statistics of the spin energy change, and to verify the fluctuation relation in equation (1), we implement an experimental protocol consisting of three steps [24]:

• (a)  
  Initialization. The system is initialized into one of the initial Hamiltonian eigenstates {ρ_±(0)}, say {ρ₊(0)}, by applying a long laser pulse that optically pumps the spin in $\left\vert 0\right\rangle$, followed by an opportune spin-rotating mw gate.
• (b)  
  Evolution. The system evolves during a time t_f under the action of the mw driving, and is intermittently perturbed by short laser pulses.
• (c)  
  Readout. At time t_f, the spin is mapped back to the $\left\{\left\vert 0\right\rangle ,\left\vert 1\right\rangle \right\}$ basis and the NV photoluminescence intensity under laser illumination is recorded to measure the spin state.

The protocol is repeated several times and the integrated photoluminescence signal gives the probability for the spin to end up in state ρ₊(t). As a result, we experimentally construct the conditional probability P_+|+(t_f) that the state ρ₊(t) is detected at time t = t_f given that the two-level-system was initialized into ρ₊(0) at time t = 0. The same procedure is applied with the system initialized into ρ₋(0) to give P_+|−(t_f).

Formally, the conditional probability is expressed as

$$\begin{equation}{P}_{+\vert {\pm}}\left({t}_{\text{f}}\right)\equiv \mathrm{Tr}\left[{\rho }_{+}\left({t}_{\text{f}}\right)\enspace \mathcal{M}\left({t}_{\text{f}}\right)\left[{\rho }_{{\pm}}\left(0\right)\right]\right],\end{equation} \tag{ 3 }$$

where $\mathcal{M}\left({t}_{\text{f}}\right)\left[{\rho }_{{\pm}}\left(0\right)\right]$ is the density operator at time t_f, expressed as a quantum map $\mathcal{M}\left({t}_{\text{f}}\right)$ applied to the initial state ρ_±(0). Notice that, by definition, P_−|± = 1 − P_+|±. Therefore, measuring P_+|+(t_f) and P_+|−(t_f) is enough to completely reconstruct the energy jump probabilities.

Together with the initial probability P_i (0), the conditional probability P_j|i (t_f), with j = ± and i = ±, provides all necessary information to construct the two-level-system energy change probability distribution function:

$$\begin{equation}p\left({\Delta}E\right)=\sum\limits _{i,j}\delta \left({\Delta}E-{E}_{j}\left({t}_{\text{f}}\right)+{E}_{i}\left(0\right)\right){P}_{j\vert i}\left({t}_{\text{f}}\right){P}_{i}\left(0\right).\end{equation} \tag{ 4 }$$

During the statistical analysis of the experimental data, the results for each of the initial states are weighted with the according Gibbs probability ${P}_{{\pm}}\left(0\right)={\text{e}}^{-\beta {E}_{{\pm}}\left(0\right)}/{Z}_{0}$ hence emulating a projective measurement being performed onto a thermal mixture at inverse temperature β. Given that we measure the conditional probabilities P_j|i (t_f) for a high number of experimental realizations (∼10⁶) and for each Hamiltonian eigenstate, this protocol is equivalent to the two-time measurement protocol applied to a thermal mixed state with inverse temperature β.

The NV center is an open quantum system: during its evolution it exchanges energy in the form of heat with the laser field (mimicking a thermal reservoir) and in the form of work with the external time dependent driving. Using the symbol Q to denote the energy gained by the reservoir, equation (1) becomes in the present case:

$$\begin{equation}\langle {\text{e}}^{-\beta {\Delta}E-{\beta }_{R}Q}\rangle ={\text{e}}^{-\beta {\Delta}F}\text{,}\end{equation} \tag{ 5 }$$

where β_R denotes the reservoir inverse (pseudo)temperature, and ΔF = −β⁻¹ ln Z(t_f)/Z(0) is the free energy difference of the NV center, between its reference equilibrium state at time t_f and its initial equilibrium state. Note that the external driving only applies to the NV center and does not affect the reservoir, therefore the reservoir free energy is constant (ΔF_R = 0) and hence it does not appear in the above equation. In equation (5) the symbol ⟨⋅⟩ denotes the average over the joint probability distribution function p(ΔE, Q), while we only have experimental access to its marginal p(ΔE). So, in principle, in order to check the validity of equation (5) one needs to perform measurements of both Q and ΔE. The marginal distribution p(ΔE) is however sufficient to check the validity of equation (1) in two special cases, namely when the reservoir temperature is infinite (β_R = 0), and when the total work W = ΔE + Q performed on the system is null, and accordingly ΔE = −Q.

In the first case equation (5) reduces to:

$$\begin{equation}\langle {\text{e}}^{-\beta {\Delta}E}\rangle ={\text{e}}^{-\beta {\Delta}F}.\end{equation} \tag{ 6 }$$

Since the heat Q does not explicitly appear in equation (6), the average reduces to the one over the marginal p(ΔE). Note that although this expression appears formally similar to the Jarzynski relation [35], in an open system as considered here the energy exchange takes into account not only work but also heat flux. In the second case, equation (5) reduces to:

$$\begin{equation}\langle {\text{e}}^{-\left(\beta -{\beta }_{R}\right){\Delta}E}\rangle ={\text{e}}^{-\beta {\Delta}F}\text{,}\end{equation} \tag{ 7 }$$

where, for the same reason as above, the average can be understood as an average over the marginal p(ΔE). Note that in our experiment the power delivered to the NV spin is non-null at all times, namely the NV spin continuously exchanges energy with the work source. However, we will focus on specific times t_f at which the total work amounts to zero, in which case equation (7) holds.

In this section we report on the experimental verification of equation (6) as a special case of the general fluctuation relation in equation (1). To that end we set an infinite pseudo-temperature of the reservoir that mimics the thermal environment (β_R = 0). That is achieved by applying the following Hamiltonian:

$$\begin{equation}\mathcal{H}\left(t\right)=\frac{\hslash \omega \left(t\right)}{2}{\sigma }_{x}.\end{equation} \tag{ 8 }$$

Applying a sequence of short laser pulses while the system evolves under this Hamiltonian results in our NV center qubit reaching the completely mixed state, corresponding to β_R = 0, as we demonstrate in the experiment and with numerical simulations, as discussed below. Work is provided via a periodic temporal variation of the Hamiltonian amplitude ω(t), which in takes the form:

$$\begin{equation}\omega \left(t\right)=\frac{{\omega }_{0}}{2}\left(1+{\mathrm{cos}}^{2}\enspace \frac{\pi t}{{\tau }_{\text{A}}}\right),\end{equation} \tag{ 9 }$$

with period set to τ_A = 616 ns, and ω₀ = π/τ_A. Under these conditions, the eigenstates of the Hamiltonian remain constant in time, while the eigenvalues are time-dependent.

We applied the protocol described in the previous section to measure the conditional probabilities P_j|i (t_f), considered as functions of t_f. Since the Hamiltonian commutes with itself at different time values, P_j|i (t_f) remain constant during the time in-between laser pulses, while they present discontinuities caused by the dissipation at the times when laser pulses are applied, as shown in figure 2(a). Each experimental point shown in this figure was obtained by averaging over an ensemble of ∼10⁶ repetitions of the same protocol.

Zoom In Zoom Out Reset image size

The asymptotic behavior is such that, in the limit of t_f → ∞, the conditional probabilities P_j|i (t_f → ∞) → 0.5 for any i and j, meaning that our dissipative channel is for practical purposes analogous to a thermal reservoir of inverse temperature β_R = 0, as anticipated above. In figures 2(b) and (c), we show the simulated mean trajectory obtained by propagating the density operators $\left\vert {E}_{+}\right\rangle \left\langle {E}_{+}\right\vert$, and $\left\vert {E}_{-}\right\rangle \left\langle {E}_{-}\right\vert$, with the dissipative map $\mathcal{M}$ that describes the dynamics. See appendix A for details on the simulation. The concatenation of coherent driving and dissipation results in a map that intermittently pushes the state toward the center of the Bloch sphere.

Having measured the conditional probabilities of energy variation, we can obtain the energy variation probability distribution function in equation (4). This is enough to reconstruct the mean value of a generic function of ΔE.

The mean value of the energy variation ⟨ΔE⟩ is shown in figure 3. The probability associated with the initial Gibbs state was set to P₊(0) = 1/(1 + e²), hence the initial inverse temperature was β = 2/ℏω₀. The experimental values retrieved for ⟨ΔE⟩ are in agreement with the theoretical values of the sum of the mean work applied on the system plus the mean heat flux, ⟨W⟩ + ⟨Q⟩, as expected from the first law of thermodynamics. These mean values are calculated as $\langle W\rangle ={\int }_{0}^{{t}_{\text{f}}}\mathrm{d}{t}^{\prime }\enspace \mathrm{Tr}\left[\rho \left({t}^{\prime }\right)\frac{\mathrm{d}\mathcal{H}\left({t}^{\prime }\right)}{\mathrm{d}{t}^{\prime }}\right]$, and $\langle Q\rangle ={\int }_{0}^{{t}_{\text{f}}}\mathrm{d}{t}^{\prime }\enspace \mathrm{Tr}\left[\frac{\mathrm{d}\rho \left({t}^{\prime }\right)}{\mathrm{d}{t}^{\prime }}\mathcal{H}\left({t}^{\prime }\right)\right]$ [36]. More details on this calculation can be found in appendix B.

Zoom In Zoom Out Reset image size

In the case considered here, the Hamiltonian time dependency is reflected in the oscillations of ⟨W⟩. On the other hand, the dissipation induced by the interaction with the laser pulses appears as discrete jumps in ⟨Q⟩. As mentioned before, the dissipation tends to equilibrate the population of the Hamiltonian eigenstates, thus the energy variation induced by the Hamiltonian driving is gradually reduced, which explains the damping on ⟨W⟩. For a time between laser pulses τ = τ_A carefully chosen to be equal to the Hamiltonian period (figure 3(b)), the flow of heat is perfectly synchronized to the times when work vanishes. Instead, when the synchronization is not perfect τ ≠ τ_A (figure 3(a)), the mean work accumulates and tends to a finite value. In these experiments, the Hamiltonian is such that $Z\left(t\right)\equiv \mathrm{Tr}\left[{\text{e}}^{-\beta \mathcal{H}\left(t\right)}\right]=2\enspace \mathrm{cosh}\left[\beta \enspace \hslash \omega \left(t\right)/2\right]$, with ω(t) defined in equation (9). Therefore, the values of ΔF = −β⁻¹ ln(Z(t_f)/Z(0)) can be computed analytically (gray lines in figure 3). The periodic time-variation of the Hamiltonian induces an oscillation with constant amplitude on ΔF. Notice that, before the very first laser pulse, when the system is closed, the difference ⟨W⟩_irr ≡ ⟨W⟩−ΔF ⩾ 0 is caused by the inner friction during the unitary evolution [37]. Given that the Hamiltonian eigenstates remain constant in time, the initial state remains unaltered by the unitary evolution. However, since the Hamiltonian eigenvalues change for a time 0 < t_f < τ_A, then the system is, on average, in a non-equilibrium state ${\rho }_{{t}_{\text{f}}}={\text{e}}^{-\beta \mathcal{H}\left(0\right)}/Z\left(0\right)$ (only describing an equilibrium thermal state for t_f = 0, τ_A, 2τ_A, ...). The quantum relative entropy between ${\rho }_{{t}_{\text{f}}}$ and the actual thermal state ${\rho }_{{t}_{\text{f}}}^{\text{th}}={\text{e}}^{-\beta \mathcal{H}\left({t}_{\text{f}}\right)}/Z\left({t}_{\text{f}}\right)$ defines the value of ⟨W⟩_irr [37, 38]. After the first laser pulse the system is no longer closed, and ΔF can exceed ⟨W⟩. In contrast, the free energy variation is always smaller or equal than the mean energy variation, as expected from Jensen's inequality applied to equation (6).

In order to effectively verify the validity of the FR in equation (6), we used the experimental values of the energy variation probability distribution to obtain ⟨e^−βΔE ⟩, and we compared it with the values of e^−βΔF = Z(t_f)/Z(0) computed analytically, as shown in figure 4. Notice that the oscillating e^−βΔF depends only on the time-varying Hamiltonian, therefore, its values are the same in figures 4(a) and (b). In contrast, ⟨e^−βΔE ⟩ depends on the time-varying Hamiltonian, and on the energy variation probability distribution function, which in turns depends on the time between laser pulses, and the strength of the dissipation, i.e., the mean amplitude of the energy jumps induced by the interaction with a single laser pulse. The fact that work and heat flux are perfectly compensated to allow equation (6) to be valid, is a consequence of the infinite pseudo-temperature of the environment. The agreement between the predicted values of e^−βΔF and the experimental values of ⟨e^−βΔE ⟩ verifies the validity of the FR in equation (6).

Zoom In Zoom Out Reset image size

Now we consider a driving Hamiltonian with fixed amplitude ω(t) ≡ ω₀ and a time varying phase such that,

$$\begin{equation}\mathcal{H}\left(t\right)=\frac{\hslash {\omega }_{0}}{2}\left({\sigma }_{x}\enspace \mathrm{cos}\enspace \theta t+{\sigma }_{y}\enspace \mathrm{sin}\enspace \theta t\right).\end{equation} \tag{ 10 }$$

In contrast to the case treated in the previous section, here the instantaneous eigenstates $\left\vert {E}_{{\pm}}\left(t\right)\right\rangle$ of $\mathcal{H}\left(t\right)$ are changing in time while the eigenvalues remain constant in time. The Hamiltonian, and accordingly its eigenvectors, are periodic functions of time, with period τ_θ = 2π/θ. As long as only the discrete times t_n = nτ_θ are concerned, the discrete-time dynamics of the system, due to its periodicity, can be conveniently expressed in terms of an effective Floquet-type Hamiltonian, reading, in this case

$$\begin{equation}{\mathcal{H}}_{\theta }=\frac{\hslash }{2}\left({\omega }_{0}{\sigma }_{x}+\theta {\sigma }_{z}\right).\end{equation} \tag{ 11 }$$

The according unitary ${U}_{\theta }={\text{e}}^{-\text{i}{\mathcal{H}}_{\theta }{\tau }_{\theta }/\hslash }$ advances the system from time nτ_θ to (n + 1)τ_θ . As opposed to the original Hamiltonian $\mathcal{H}\left(t\right)$, the discrete-time Hamiltonian ${\mathcal{H}}_{\theta }$ is time-independent. Accordingly, in the discrete-time picture no work source appears to be at play.

In our experiment the two-level system is prepared in one of the two eigenstates of ${\mathcal{H}}_{\theta }$, reading

$$\begin{equation}\left\vert {E}_{\theta ,{\uparrow}}\right\rangle \equiv \mathrm{sin}\left(\alpha /2\right)\left\vert 1\right\rangle +\mathrm{cos}\left(\alpha /2\right)\left\vert 0\right\rangle \end{equation} \tag{ 12 }$$

$$\begin{equation}\left\vert {E}_{\theta ,{\downarrow}}\right\rangle \equiv \mathrm{cos}\left(\alpha /2\right)\left\vert 1\right\rangle -\mathrm{sin}\left(\alpha /2\right)\left\vert 0\right\rangle ,\end{equation} \tag{ 13 }$$

where $\alpha =-\mathrm{arctan}\left({\omega }_{0}/\theta \right)$. The probability P_↑|↑ to start in the $\left\vert {\uparrow}\right\rangle$ state and find the system again in the $\left\vert {\uparrow}\right\rangle$ state after a time t_f reads ${P}_{{\uparrow}\vert {\uparrow}}=1-\frac{{\omega }_{0}^{2}}{{\omega }_{0}^{2}+{\theta }^{2}}\enspace {\mathrm{sin}}^{2}\enspace \frac{t\theta }{2}$. Similarly, it is ${P}_{{\uparrow}\vert {\downarrow}}=\frac{{\omega }_{0}^{2}}{{\omega }_{0}^{2}+{\theta }^{2}}\enspace {\mathrm{sin}}^{2}\enspace \frac{t\theta }{2}$. Note that these quantities are, as expected, periodic with period τ_θ . In addition, P_↑|↑ = 1 and P_↓|↑ = 0 at stroboscopic times t_n = nτ_θ , in accordance with the fact that $\left\vert {E}_{\theta ,{\uparrow}}\right\rangle$ and $\left\vert {E}_{\theta ,{\downarrow}}\right\rangle$ are eigenvectors of the discrete-time Hamiltonian ${\mathcal{H}}_{\theta }$. Figure 5(a) shows the conditional probabilities P_↑|↑ and P_↑|↓ measured in absence of laser pulses. Note the agreement with the analytic expressions, confirming that $\mathcal{H}\left(t\right)$ well describes the dynamics in absence of laser pulses.

Zoom In Zoom Out Reset image size

In the further experiments, short laser pulses are applied at the discrete times t_n = nτ_θ , and final energy measurements are performed at a final time t_f = N_L τ_θ . Accordingly, the spin free energy difference ΔF is null, and equation (7) takes on the simpler form

$$\begin{equation}\langle {\text{e}}^{-\left(\beta -{\beta }_{R}\right){\Delta}E}\rangle =1.\end{equation} \tag{ 14 }$$

The experimental data are shown together with a Monte-Carlo simulation (see appendix A). Figures 5(b)–(d) shows the impact of the laser pulses on the conditional probabilities P_↑|↑ and P_↑|↓ as a function of N_L. Note how these probabilities converge to the same asymptotic value ${P}_{{\uparrow}}^{\infty }$, meaning that, regardless of the initial condition, the system ends up in the $\left\vert {\uparrow}\right\rangle$ state with probability ${P}_{{\uparrow}}^{\infty }$. In other words, the system reaches a steady state which conforms to a Gibbs state of inverse temperature ${\beta }_{R}=-\frac{1}{\hslash \omega } \mathrm{ln} \left(\frac{{P}_{{\uparrow}}^{\infty }}{1-{P}_{{\uparrow}}^{\infty }}\right)$. The value of ${P}_{{\uparrow}}^{\infty }$, and therefore the value of β_R , depends on the ratio ω₀/θ. In the adiabatic regime, θ ≪ ω₀, the Hamiltonian is approximately proportional to σ_x , and ${P}_{{\uparrow}}^{\infty }\simeq 1/2$. Hence, in the adiabatic regime the value of β_R tends to zero. As the value of θ increases with respect to ω₀, the dynamics no longer belongs to the adiabatic regime, and the value ${P}_{{\uparrow}}^{\infty }$ approaches to zero, hence increasing the value of β_R . Note that, although the spin is driven by a time dependent Hamiltonian that delivers power continuously, at discrete times in which the energy is measured the system evolves as under the influence of an effective thermal reservoir of inverse temperature β_R and as if it were not in contact with any work source. Given these conditions, the mean energy variation is equal to the mean heat flux ⟨Q⟩, as demonstrated in figures 6(a)–(c) by directly comparing ⟨ΔE⟩, extracted from the experimental data, with the mean heat flux calculated with the theoretical model. See appendix B for details on how to calculate ⟨Q⟩ in this case. Remarkably, even when ⟨ΔE⟩ is different from zero, we demonstrate that the energy variation statistics follows the quantum fluctuation relation $\langle {\text{e}}^{-\left(\beta -{\beta }_{R}\right){\Delta}E}\rangle =1$ (equation (14)), as the general theory predicts. Figures 6(d)–(f) shows the experimental values of the quantity $\langle {\text{e}}^{-\left(\beta -{\beta }_{R}\right){\Delta}E}\rangle$, as obtained under various experimental conditions (specifically with various values of τ_θ and α). The experimental data confirm the validity of the quantum fluctuation relation as expressed in equation (7).

Zoom In Zoom Out Reset image size

The experimental verification of quantum fluctuation relations for open systems that simultaneously exchange heat—due to a thermal gradient—and work—due to a time-dependent driving—is a formidable challenge, owing to the general difficulty of measuring heat and work separately, while having typically experimental access to their difference ΔE.

With this work we establish a first milestone along the way toward a fully fledged experimental demonstration of fluctuation relations for driven-dissipative systems. The verification becomes feasible (i.e., it requires the statistics of ΔE only) in two special cases: (i) the external environment is at infinite temperature, and (ii) the system is subject to a periodic driving, so that, at multiples of the driving period the total work amounts to zero (while the power is generally non-null). In those cases the fluctuation relations boil down to the simplified expressions in equations (6) and (7), respectively, which we have experimentally verified.

Our experimental platform of choice was a single NV center in diamond. Due to its long coherence time, and the possibility to finely control its state and its interaction with engineered environments, this system naturally lend itself to study non-equilibrium thermodynamics in the quantum regime [24].

Our results represent the first experimental study of quantum fluctuation relation in driven open quantum systems.

The authors would like to thank P Cappellaro, F Poggiali, and S Gherardini for useful discussions, and F S Cataliotti and M Inguscio for continuous support.

All data that support the findings of this study are included within the article (and any supplementary files).

In order to model the qubit dynamics, we separate the evolution into a unitary part, describing the dynamics between consecutive short laser pulses, and the irreversible dynamics that describe the interaction with laser pulses.

The unitary operators of the map correspond to

$$\begin{equation}U\left({t}_{{N}_{\text{L}}-1},{t}_{{N}_{\text{L}}}\right)=\mathrm{exp}\left(-i{\int }_{{t}_{{N}_{\text{L}}-1}}^{{t}_{{N}_{\text{L}}}}\mathrm{d}{t}^{\prime }\enspace \mathcal{H}\left({t}^{\prime }\right)\right),\end{equation} \tag{ A.1 }$$

with t_N ≡ Nτ. We recall that τ is the time between laser pulses. In the case of the Hamiltonian in equation (8), considered in section 3, each unitary operator describes a rotation around the x axis, such that $U\left({t}_{{N}_{\text{L}}-1},{t}_{{N}_{\text{L}}}\right)=\mathrm{exp}\left(-i{\sigma }_{x}{\int }_{{t}_{{N}_{\text{L}}-1}}^{{t}_{{N}_{\text{L}}}}\mathrm{d}{t}^{\prime }\enspace \omega \left({t}^{\prime }\right)\right)$, with ω(t) defined by equation (9). In contrast, in the case studied in section 4, the unitary evolution is described by the Floquet Hamiltonian ${\mathcal{H}}_{\theta }=\frac{\hslash }{2}\left({\omega }_{0}{\sigma }_{x}+\theta {\sigma }_{z}\right)$ (equation (11)) only for discrete times t_n = nτ_θ , such that, for interpulse times τ = τ_θ , the unitary operator corresponds to $U\left({t}_{{N}_{\text{L}}-1},{t}_{{N}_{\text{L}}}\right)={U}_{\theta }=\mathrm{exp}\left(-i{\mathcal{H}}_{\theta }{\tau }_{\theta }/\hslash \right)$. Both of these unitary operators can be solved analytically and then included in a numerical model of the dynamics.

The dissipative dynamics associated with each laser pulse consists instead in an energy jump that is modeled as a POVM, with (1 − p_a) probability to leave the state unaffected, and a probability p_a to absorb photons. Upon absorption, the spin qubit is subject to a projective measurement of σ_z followed by an optical pumping of the spin modeled via Lindbladian master equation, with a Lindbladian jump operator $\sqrt{{\Gamma}}\left\vert 0\right\rangle \left\langle 1\right\vert$ that transfers spin population to the state $\left\vert 0\right\rangle$. The numerical simulation is then performed by using a Monte-Carlo simulation to determine the result of the POVM, hence determining if the laser pulse is absorbed or not. If the pulse is absorbed, Lindbladian operator is applied to the state of the system. In this way, we simulate single trajectories followed by the system, and we finally average over several realizations to obtain the mean trajectories shown in figures 2 and 5.

As described in the main text, the mean values of work and heat flux can be calculated as $\langle W\rangle ={\int }_{0}^{{t}_{\text{f}}}\mathrm{d}{t}^{\prime }\enspace \mathrm{Tr}\enspace \left[\rho \left({t}^{\prime }\right)\frac{\mathrm{d}\mathcal{H}\left({t}^{\prime }\right)}{\mathrm{d}{t}^{\prime }}\right]$, and $\langle Q\rangle ={\int }_{0}^{{t}_{\text{f}}}\mathrm{d}{t}^{\prime }\enspace \mathrm{Tr}\enspace \left[\frac{\mathrm{d}\rho \left({t}^{\prime }\right)}{\mathrm{d}{t}^{\prime }}\mathcal{H}\left({t}^{\prime }\right)\right]$ [36]. In this appendix we describe how to compute this quantities for the dynamics described in sections 3 and 4.

B.1. Mean work and mean heat flux for section 3

The spin qubit is driven by the Hamiltonian (equation (8)) $\mathcal{H}\left(t\right)=\frac{\hslash \omega \left(t\right)}{2}{\sigma }_{x}$, where ω(t) is a periodic function such that ω(t) = ω(t + τ_A). In addition the system interacts with short laser pulses applied at times t_f = N_L τ, where τ is the time between consecutive pulses and N_L is the number of pulses. Assuming that the duration of every short laser pulse interaction is given by δ, then the mean work and mean heat flux after the first laser pulse is given by

$$\begin{align}{\langle W\rangle }_{1}& ={\int }_{0}^{\tau }\mathrm{d}{t}^{\prime }\enspace \mathrm{Tr}\left[\rho \left({t}^{\prime }\right)\frac{\mathrm{d}\mathcal{H}\left({t}^{\prime }\right)}{\mathrm{d}{t}^{\prime }}\right]+{\int }_{\tau }^{\tau +\delta }\mathrm{d}{t}^{\prime }\enspace \mathrm{Tr}\left[\rho \left({t}^{\prime }\right)\frac{\mathrm{d}\mathcal{H}\left(\tau \right)}{\mathrm{d}{t}^{\prime }}\right]\\ & ={\int }_{0}^{\tau }\mathrm{d}{t}^{\prime }\enspace \mathrm{Tr}\left[\rho \left({t}^{\prime }\right)\frac{\mathrm{d}\mathcal{H}\left({t}^{\prime }\right)}{\mathrm{d}{t}^{\prime }}\right]\end{align} \tag{ B.1 }$$

$$\begin{align}{\langle Q\rangle }_{1}& ={\int }_{0}^{\tau }\mathrm{d}{t}^{\prime }\enspace \mathrm{Tr}\left[\frac{\mathrm{d}\rho \left({t}^{\prime }\right)}{\mathrm{d}{t}^{\prime }}\mathcal{H}\left({t}^{\prime }\right)\right]+{\int }_{\tau }^{\tau +\delta }\mathrm{d}{t}^{\prime }\enspace \mathrm{Tr}\left[\frac{\mathrm{d}\rho \left({t}^{\prime }\right)}{\mathrm{d}{t}^{\prime }}\mathcal{H}\left(\tau \right)\right]\\ & ={\int }_{\tau }^{\tau +\delta }\mathrm{d}{t}^{\prime }\enspace \mathrm{Tr}\left[\frac{\mathrm{d}\rho \left({t}^{\prime }\right)}{\mathrm{d}{t}^{\prime }}\mathcal{H}\left(\tau \right)\right].\end{align} \tag{ B.2 }$$

The second term in equation (B.1) vanishes because we assume the interaction with the laser pulses to be instantaneous, δ ≪ τ and δ ≪ τ_A, hence the Hamiltonian is kept constant during the interaction with a single laser pulse. Instead, the first term in equation (B.2) is equal to zero because the system evolves under unitary dynamics during the integration interval, meaning that $\frac{\mathrm{d}\rho }{\mathrm{d}t}=\frac{1}{i\hslash }\left[\mathcal{H},\rho \right]$. Then, the cyclic property of the trace implies that $\mathrm{Tr}\left[\frac{\mathrm{d}\rho }{\mathrm{d}t}\mathcal{H}\right]=\frac{1}{i\hslash }\left(\mathrm{Tr}\left[\mathcal{H}\rho \mathcal{H}\right]-\mathrm{Tr}\left[\rho {\mathcal{H}}^{2}\right]\right)=0$. Therefore, for a given final time t_f such that N_L τ < t_f < (N_L + 1)τ we can separate the mean work and heat flux as

$$\begin{equation}\langle W\rangle =\sum\limits _{n=1}^{{N}_{\text{L}}}{\langle W\rangle }_{n}+{\int }_{\tau {N}_{\text{L}}}^{{t}_{\text{f}}}\mathrm{d}{t}^{\prime }\enspace \mathrm{Tr}\left[\rho \left({t}^{\prime }\right)\frac{\mathrm{d}\mathcal{H}\left({t}^{\prime }\right)}{\mathrm{d}{t}^{\prime }}\right]\end{equation} \tag{ B.3 }$$

$$\begin{equation}\langle Q\rangle =\sum\limits _{n=1}^{{N}_{\text{L}}}{\langle Q\rangle }_{n},\end{equation} \tag{ B.4 }$$

where

$$\begin{align}{\langle W\rangle }_{n}& \equiv {\int }_{\left(n-1\right)\tau }^{n\tau }\mathrm{d}{t}^{\prime }\enspace \mathrm{Tr}\left[\rho \left({t}^{\prime }\right)\frac{\mathrm{d}\mathcal{H}\left({t}^{\prime }\right)}{\mathrm{d}{t}^{\prime }}\right]\\ & =\frac{\hslash }{2}{\int }_{\left(n-1\right)\tau }^{n\tau }\mathrm{d}{t}^{\prime }\frac{\mathrm{d}\omega \left({t}^{\prime }\right)}{\mathrm{d}{t}^{\prime }}\enspace \mathrm{Tr}\left[\rho \left({t}^{\prime }\right){\sigma }_{x}\right]\end{align} \tag{ B.5 }$$

$$\begin{align}{\langle Q\rangle }_{n}& \equiv {\int }_{n\tau }^{n\tau +\delta }\mathrm{d}{t}^{\prime }\enspace \mathrm{Tr}\left[\frac{\mathrm{d}\rho \left({t}^{\prime }\right)}{\mathrm{d}{t}^{\prime }}\mathcal{H}\left(\tau \right)\right]\\ & =\mathrm{Tr}\left[\left(\rho \left(n\tau +\delta \right)-\rho \left(n\tau \right)\right)\mathcal{H}\left(\tau \right)\right].\end{align} \tag{ B.6 }$$

Note that the density operator at any given time t can be written as

$$\begin{equation}\rho \left(t\right)={P}_{+}\left(t\right)\left\vert {E}_{+}\right\rangle \left\langle {E}_{+}\right\vert +\left(1-{P}_{+}\left(t\right)\right)\left\vert {E}_{-}\right\rangle \left\langle {E}_{-}\right\vert +\left(q\left(t\right)\left\vert {E}_{-}\right\rangle \left\langle {E}_{+}\right\vert +\mathrm{c}.\mathrm{c}.\right),\end{equation} \tag{ B.7 }$$

where 0 ⩽ P₊(t) ⩽ 1 is the population on the eigenstate $\left\vert {E}_{+}\right\rangle$, and q(t) is a complex number describing the coherence of the state in the Hamiltonian basis. We recall that $\left\vert {E}_{{\pm}}\right\rangle$ are eigenstates of σ_x with eigenvalues ±1, which implies that Tr[ρ(t)σ_x ] = 2P₊(t) − 1. Therefore, equation (B.6) can be rewritten as

$$\begin{equation}{\langle Q\rangle }_{n}=\hslash \omega \left(n\tau \right)\left({P}_{+}\left(n\tau \right)-{P}_{+}\left(\left(n-1\right)\tau \right)\right).\end{equation} \tag{ B.8 }$$

In addition, since the Hamiltonian is always parallel to σ_x , then the population P₊(t) during the unitary evolution is kept constant. Therefore, equation (B.5) can be recast as

$$\begin{equation}{\langle W\rangle }_{n}=\frac{\hslash }{2}\left(2{P}_{+}\left(\left(n-1\right)\tau \right)-1\right)\left(\omega \left(n\tau \right)-\omega \left(\left(n-1\right)\tau \right)\right).\end{equation} \tag{ B.9 }$$

This means that knowing the population at times t_f = nτ is enough to calculate ⟨W⟩_n and ⟨Q⟩_n . The population, as we just mentioned, are unaltered by the unitary evolution, so we need to understand how they are affected by the short laser pulses. As described in appendix A, the interaction with a single laser pulse has a probability 1 − p_a of leaving the system unaltered, and a probability p_a to project the state of the system into one of the eigenstates of σ_z , and then applying a Lindbladian evolution that pumps population from $\left\vert 1\right\rangle \left\langle 1\right\vert$ to $\left\vert 0\right\rangle \left\langle 0\right\vert$. Hence, the state of the system after absorbing a laser pulse is diagonal in the σ_z basis, meaning that its populations in the σ_x basis must be equal to 1/2. Therefore, we can write a recursive equation for the population in the σ_x basis

$$\begin{equation}{P}_{+}\left(n\tau \right)=\left(1-{p}_{\text{a}}\right){P}_{+}\left(\left(n-1\right)\tau \right)+{p}_{\text{a}}/2.\end{equation} \tag{ B.10 }$$

Notice that P₊(0) = 1/(1 + e), as mentioned in the main text. Therefore, equation (B.10) can be used recursively to obtain

$$\begin{equation}{P}_{+}\left(n\tau \right)=\frac{1}{2}\left(1-{\left(1-{p}_{\text{a}}\right)}^{n}\left(1-2{P}_{+}\left(0\right)\right)\right).\end{equation} \tag{ B.11 }$$

Finally, we can recast equations (B.9) and (B.8) as

$$\begin{equation}{\langle W\rangle }_{n}=\frac{\hslash }{2}\left(\omega \left(\left(n-1\right)\tau \right)-\omega \left(n\tau \right)\right){\left(1-{p}_{\text{a}}\right)}^{n-1}\left(1-2{P}_{+}\left(0\right)\right)\end{equation} \tag{ B.12 }$$

$$\begin{equation}{\langle Q\rangle }_{n}=\frac{\hslash }{2}\omega \left(n\tau \right){p}_{\text{a}}{\left(1-{p}_{\text{a}}\right)}^{n-1}\left(1-2{P}_{+}\left(0\right)\right).\end{equation} \tag{ B.13 }$$

Equations (B.12) and (B.13) can be used to calculate the values of ⟨W⟩ and ⟨Q⟩ (equations (B.3) and (B.4)). The mean values of work and heat flux are therefore completely defined by the parameters N_L, P₊(0), t_f, τ/τ_A, and p_a, the latter being the only free parameter fitted to the experimental data in figure 2.

B.2. Mean heat flux for section 4

For the stroboscopic times described in section 4, the unitary evolution of the system is described by ${U}_{\theta }={\text{e}}^{-\text{i}{\mathcal{H}}_{\theta }{\tau }_{\theta }/\hslash }$, where ${\mathcal{H}}_{\theta }$ is the time independent Floquet Hamiltonian in equation (11), and τ_θ is the time between two consecutive laser pulses, and the time at which the energy of the system is measured. Given that ${\mathcal{H}}_{\theta }$ is time independent, then the mean work measured at stroboscopic times t_f = τN_L is equal to zero. On the other hand, in analogy to appendix B.1, the mean heat flux measured at those stroboscopic times is

$$\begin{align}\langle Q\rangle & =\sum\limits _{n=1}^{{N}_{\text{L}}}{\langle Q\rangle }_{n}=\sum\limits _{n=1}^{{N}_{\text{L}}}{\int }_{n\tau }^{n\tau +\delta }\mathrm{d}{t}^{\prime }\enspace \mathrm{Tr}\left[\frac{\mathrm{d}\rho \left({t}^{\prime }\right)}{\mathrm{d}{t}^{\prime }}{\mathcal{H}}_{\theta }\right]\\ & =\sum\limits _{n=1}^{{N}_{\text{L}}}2{E}_{\theta }\left({P}_{{\uparrow}}\left(n\tau \right)-{P}_{{\uparrow}}\left(\left(n-1\right)\tau \right)\right)\\ & =2{E}_{\theta }\left({P}_{{\uparrow}}\left({N}_{\text{L}}\tau \right)-{P}_{{\uparrow}}\left(0\right)\right),\end{align} \tag{ B.14 }$$

where ${E}_{\theta }=\hslash \sqrt{{\omega }_{0}^{2}+{\theta }^{2}}/2$.

Let us assume that the density operator ρ(τ) represents the state of the system after applying a single short laser pulse to the state ρ(0). As explained in appendix A, the qubit interaction with a short laser pulse is modeled as a POVM with probability (1 − p_a) to leave unaltered the state of the qubit and a probability p_a to apply a projective measurement of σ_z followed by Lindbladian dynamics under the jump operator $\sqrt{{\Gamma}}\left\vert 0\right\rangle \left\langle 1\right\vert$. Under this model, the probability to obtain E_θ,↑ = E_θ when measuring the energy of ρ(τ) is

$$\begin{align}{P}_{{\uparrow}}\left(\tau \right)& =\mathrm{Tr}\left[\left\vert {E}_{\theta ,{\uparrow}}\right\rangle \left\langle {E}_{\theta ,{\uparrow}}\right\vert \rho \left(\tau \right)\right]\\ & =\left(1-{p}_{\text{a}}k\right){P}_{{\uparrow}}\left(0\right)+\frac{1}{2}{p}_{\text{a}}\left(k-{p}_{\mathrm{d}}\enspace \mathrm{cos}\enspace \alpha \right),\end{align} \tag{ B.15 }$$

where p_d ≡ 1 − e^−Γδ is the optical pumping probability associated with the Lindbladian dynamics, $\alpha =-\mathrm{arctan}\left({\omega }_{0}/\theta \right)$ is the angle that defines the ratio between the terms in the Floquet Hamiltonian (equation (11)), ${P}_{{\uparrow}}\left(0\right)=\mathrm{Tr}\left[\left\vert {E}_{\theta ,{\uparrow}}\right\rangle \left\langle {E}_{\theta ,{\uparrow}}\right\vert \rho \left(0\right)\right]$ is the probability associated with the initial state, and where we have defined, to simplify notation, a factor

$$\begin{equation}k\equiv 1+\left(1-{p}_{\mathrm{d}}\right)\mathrm{cos}\enspace {\alpha }^{2}.\end{equation} \tag{ B.16 }$$

The relation in equation (B.15) can be recursively applied in order to obtain

$$\begin{equation}{P}_{{\uparrow}}\left({N}_{\text{L}}\tau \right)={\left(1-{p}_{\text{a}}k\right)}^{{N}_{\text{L}}}{P}_{{\uparrow}}\left(0\right)+\frac{1}{2}\left(1-{\left(1-{p}_{\text{a}}k\right)}^{{N}_{\text{L}}}\right)\left(1-\frac{{p}_{\mathrm{d}}}{k}\enspace \mathrm{cos}\enspace \alpha \right).\end{equation} \tag{ B.17 }$$

Therefore, the mean heat flux in equation (B.14) can be recast as

$$\begin{equation}\langle Q\rangle ={E}_{\theta }\left(1-{\left(1-{p}_{\text{a}}k\right)}^{{N}_{\text{L}}}\right)\left(1-\frac{{p}_{\mathrm{d}}}{k}\enspace \mathrm{cos}\enspace \alpha -2{P}_{{\uparrow}}\left(0\right)\right).\end{equation} \tag{ B.18 }$$

While p_a and p_d depend on the laser pulse duration and intensity, all the other parameters are precisely set in experiment, with negligible uncertainty. Note that p_d can be extracted from the asymptotic state probability ${P}_{{\uparrow}}^{\infty }={\mathrm{lim}}_{{N}_{\text{L}}\to \infty }\enspace {P}_{{\uparrow}}\left({N}_{\text{L}}\tau \right)=\frac{1}{2}\left(1-\frac{{p}_{\mathrm{d}}}{k}\enspace \mathrm{cos}\enspace \alpha \right)$. Thus, the only free parameter in equation (B.18) is the laser absorption probability p_a, which we fit to the experimental data (see figure 5).