Harnessing Fluctuations in Thermodynamic Computing via Time-Reversal Symmetries

TL;DR

This paper experimentally reveals structured work distributions during bit erasure, linked to time-reversal symmetries, and introduces fluctuation theorems that improve understanding of entropy production and dissipation in nanoscale thermodynamic computing.

Contribution

It introduces a new class of fluctuation theorems based on time-reversal symmetries, providing a diagnostic tool for thermodynamic computing at the nanoscale.

Findings

Structured work distributions observed during bit erasure.

Verification of new fluctuation theorems in superconducting circuits.

Enhanced estimates of dissipated heat and entropy production.

Abstract

We experimentally demonstrate that highly structured distributions of work emerge during even the simple task of erasing a single bit. These are signatures of a refined suite of time-reversal symmetries in distinct functional classes of microscopic trajectories. As a consequence, we introduce a broad family of conditional fluctuation theorems that the component work distributions must satisfy. Since they identify entropy production, the component work distributions encode both the frequency of various mechanisms of success and failure during computing, as well giving improved estimates of the total irreversibly-dissipated heat. This new diagnostic tool provides strong evidence that thermodynamic computing at the nanoscale can be constructively harnessed. We experimentally verify this functional decomposition and the new class of fluctuation theorems by measuring transitions between flux states in a superconducting circuit.

Tables

Table 1. Table S1 : Erasure protocol.

Stage		Drop Barrier		Tilt		Raise Barrier		Untilt
$t_s$	$t_0$		$t_1$		$t_2$		$t_3$		$t_4$
$\beta \left(t\right)$	$|$	$\frac{t_1-t}{t_1-t_0}$	$|$	0	$|$	$\frac{t-t_2}{t_3-t_2}$	$|$	$1$	$|$
$\delta \left(t\right)$	$|$	0	$|$	$\frac{t-t_1}{t_2-t_1}$	$|$	$1$	$|$	$\frac{t_4-t}{t_4-t_3}$	$|$

Table 2. Table S2 : (Top Left) Comparison of simulated class-average works and bounds for different trajectory classes: All trajectories 𝒵 → → 𝒵 {\overrightarrow{\mathcal{Z}}} , Success trajectories, Fail trajectories, and Transitional trajectories. These are identified in Fig. 2 (Four front panels). From left to right, columns give the estimated class-average work ⟨ W ⟩ C subscript delimited-⟨⟩ 𝑊 𝐶 \langle W\rangle_{C} , TCFT lower bound ⟨ W ⟩ C min superscript subscript delimited-⟨⟩ 𝑊 𝐶 min \langle W\rangle_{C}^{\text{min}} , and their difference E C subscript 𝐸 𝐶 E_{C} . 3.5 × 10 6 3.5 superscript 10 6 3.5\times 10^{6} simulations were run for each of the forward and reverse processes, with 96.2 % percent 96.2 96.2\% trajectories successfully ending in the R 𝑅 R informational state under the forward process. (Top Right) Comparison of ensemble-average work and bounds due to different partitions: Trivial partition; Untilt-Centric I partition, composed of Success, Fail, and Transitional; and Untilt-Centric II partition, described in follow-on work. From left to right, columns give the estimated ensemble-average work, the partition bound ⟨ W ⟩ Q min superscript subscript delimited-⟨⟩ 𝑊 𝑄 min \langle W\rangle_{Q}^{\text{min}} , and their difference E Q subscript 𝐸 𝑄 E_{Q} . All values in units of k B ​ T subscript 𝑘 B 𝑇 k_{\text{B}}T . (Bottom) Parallel results from the flux qubit experiment. All trajectories 𝒵 → → 𝒵 {\overrightarrow{\mathcal{Z}}} , Success , Fail , and Transitional trajectories identified in Fig. 3 (D). Data from 195 , 050 195 050 195,050 trajectories from the forward protocol and 250 , 000 250 000 250,000 trajectories from the reverse protocol.

Class $C$	${⟨W⟩}_C$	${⟨W⟩}_C^{\text{min}}$	$E_C$
All $\overrightarrow{𝒵}$	0.634	0.0	0.634
Success	0.713	0.683	0.030
Fail	-3.885	-3.951	0.066
Transitional	-0.546	-1.650	1.170

Equations

$$⟨W{⟩}_C$$

$$=⟨W{⟩}_C^{\text{min}},$$

$$\Delta F=-k_{\text{B}}T\ln \left(\frac{\mathcal{P}(C)}{\mathcal{R}(C^{\text{R}})}⟨e^{-W/k_{\text{B}}T}{⟩}_C\right).$$

$$V(x,t)=a{x}^4-b_0\beta (t)x^2-c_0\delta (t)x,$$

$$dx$$

$$dv$$

$$⟨e^{-W/k_{\text{B}}T}{⟩}_C=\frac{\mathcal{R}(C^{\text{R}})}{\mathcal{P}(C)}{e}^{-\Delta F/k_{\text{B}}T},$$

$$⟨W{⟩}_C$$

$$\Delta F=-k_{\text{B}}T\ln \left(\frac{\mathcal{P}(C)}{\mathcal{R}(C^{\text{R}})}⟨e^{-W/k_{\text{B}}T}{⟩}_C\right).$$

$$\Delta F=-k_{\text{B}}T\ln ⟨e^{-W/k_{\text{B}}T}{⟩}_{\overrightarrow{\mathcal{Z}}},$$

$$E_C$$

$$=k_{B}T{D}_{\text{KL}}\left[\mathcal{P}(\overrightarrow{\mathcal{Z}}|C)||\mathcal{R}({\overrightarrow{\mathcal{Z}}}^{\text{R}}|C^{\text{R}})\right]$$

$$\ge 0.$$

$$E_Q$$

$$=⟨W{⟩}_{\overrightarrow{\mathcal{Z}}}-\sum _{C\in Q}\mathcal{P}(C)⟨W{⟩}_C^{\text{min}}$$

$$=⟨W{⟩}_{\overrightarrow{\mathcal{Z}}}-⟨W{⟩}_Q^{\text{min}},$$

$$⟨W{⟩}_Q^{\text{min}}$$

$$⟨W⟩$$

$$\displaystyle=\Delta F+k_{\text{B}}T\bigg{(}\sum_{C\in Q}\mathcal{P}(C)D_{\text{KL}}\left[\mathcal{P}({\overrightarrow{\mathcal{Z}}}|C)||\mathcal{R}({\overrightarrow{\mathcal{Z}}}^{\texttt{R}}|C^{\texttt{R}})\right]+D_{\text{KL}}\left[\mathcal{P}(Q)||\mathcal{R}(Q^{\texttt{R}})\right]\bigg{)}~{},$$

$$\Delta S_{\text{sys}}=\sum _z\left[-{\pi }_{\tau }(z)\ln {\pi }_{\tau }(z)+{\pi }_0(z)\ln {\pi }_0(z)\right].$$

$$\Delta s_{\text{sys}}(\overrightarrow{z})=\ln \frac{{\pi }_0(z_0)}{{\pi }_{\tau }(z_{\tau })},$$

$$\Delta S_{\text{res}}(\overrightarrow{z})=-\beta Q(\overrightarrow{z})=\ln \frac{p(\overrightarrow{z}|z_0)}{r({\overrightarrow{z}}^{\text{R}}|z_{\tau }^{\text{R}})},$$

$$\Delta S_{\text{tot}}(\overrightarrow{z})$$

$$=\beta \left({\mathcal{H}}_{\tau }(z_{\tau })-F_{\tau }\right)-\beta \left({\mathcal{H}}_0(z_0)-F_0\right)-\beta Q(\overrightarrow{z})$$

$$=\beta \left(\Delta \mathcal{H}(\overrightarrow{z})-Q(\overrightarrow{z})-\Delta F\right)$$

$$=\beta \left(W(\overrightarrow{z})-\Delta F\right).$$

$$\Delta S_{\text{tot}}(\overrightarrow{z})$$

$$=\ln \frac{{\pi }_0(z_0)}{{\pi }_{\tau }(z_{\tau })}\frac{p(\overrightarrow{z}|z_0)}{r({\overrightarrow{z}}^{\text{R}}|z_{\tau }^{\text{R}})}$$

$$=\ln \frac{\mathcal{P}(\overrightarrow{z})}{\mathcal{R}({\overrightarrow{z}}^{\text{R}})}$$

$$\mathcal{P}(\overrightarrow{z})$$

$$\mathcal{R}({\overrightarrow{z}}^{\text{R}})$$

$$\int d\overrightarrow{z}[\overrightarrow{z}\in C]\mathcal{R}({\overrightarrow{z}}^{\text{R}})$$

$$=\int d{\overrightarrow{z}}^{\text{R}}[{\overrightarrow{z}}^{\text{R}}\in C^{\text{R}}]\mathcal{R}({\overrightarrow{z}}^{\text{R}})$$

$$=\int d\overrightarrow{z}[\overrightarrow{z}\in C^{\text{R}}]\mathcal{R}(\overrightarrow{z})$$

$$=\mathcal{R}(C^{\text{R}}).$$

$$\int d\overrightarrow{z}[\overrightarrow{z}\in C]\mathcal{P}(\overrightarrow{z})e^{-\beta \left(W(\overrightarrow{z})-\Delta F\right)}$$

$$=e^{\beta \Delta F}\int d\overrightarrow{z}\mathcal{P}(\overrightarrow{z},C)e^{-\beta W(\overrightarrow{z})}$$

Full text

**Jarzynski integral fluctuation theorem, Crooks detailed fluctuation theorem, Landauer’s Principle, thermodynamics of computation, information thermodynamics **

Physics dictates that all computing is subject to spontaneous error. These days, this truism repeatedly reveals itself: despite the once-predictable miniaturization of nanoscale electronics, computing performance increases have dramatically slowed in the last decade or so. In large measure, this is due to the concomitant rapid decrease in the number of information-bearing physical degrees of freedom, rendering information storage and processing increasingly susceptible to corruption by thermal fluctuations. All computing is thermodynamic. Controlling the production of fluctuations and removing heat pose key technological challenges to further progress. One comparable setting that gives some optimism, though, is the overtly functional behavior exhibited by biological cells—presumably functional information processing by small numbers of molecules subject to substantial thermal fluctuations. Computing technologies are very far away from this kind of robust information processing.

Only recently have tools appeared that precisely describe what trade-offs exist between thermodynamic resources and useful information processing—these are highly reminiscent of the centuries-old puzzle of how Maxwell’s “very observant and neat-fingered” demon uses its “intelligence” to convert disorganized heat energy to useful work [1]. In our modern era, his demon has led to the realization that information itself is physical [2, 3, 4]—or, most constructively, that information is a thermodynamic resource [5]. This has opened up the new paradigm of thermodynamic computing in which fluctuations play a positive role in efficient information processing on the nanoscale. We now conceptualize this via information engines: physical systems that are driven by, manipulate, store, and dissipate energy, but simultaneously generate, store, lose, communicate, and transform information. In short, information engines combine traditional engines comprised of heat, work, and other familiar reservoirs with, what we now call, information reservoirs [6, 7].

Reliable thermodynamic computing requires detecting and controlling fluctuations in informational and energetic resources and in engine functioning. For this, one appeals to fluctuation theorems that capture exact time-reversal symmetries and predict entropy production leading to irreversible dissipation [8, 9, 10, 11, 12, 13, 14]. We are now on the door-step of the very far-from-equilibrium thermodynamics needed to understand the physics of computing. And, in turn, this has started to reveal the physical principles of how nature processes information in the service of biological functioning and survival.

Proof-of-concept experimental tests have been carried out in several substrates: probing biomolecule free energies [15, 16, 17], work expended during elementary computing (bit erasure) [18, 19, 20, 21, 22, 23], and Maxwellian demons [24]. That said, the suite of contemporary principles (Supplementary Materials (SM) I) far outstrips experimental validation to date.

To close the gap, we show how to diagnose thermodynamic computing on the nanoscale by explaining the signature structures in work distributions generated during information processing. These structures track the mesoscale evolution of a system’s informational states and reveal classes of functional and nonfunctional microscopic trajectories. We show that the informational-state evolutions are identified by appropriate conditioning and that they obey a suite of trajectory-class fluctuation theorems, which give accurate bounds on work, entropy production, and dissipation. The result is a new tool that employs mesoscopic measurements to diagnose nanoscale thermodynamic computing. For simplicity and to make direct contact with previous efforts, we demonstrate the tools on Landauer erasure of a bit of information in a superconducting flux qubit.

As a reference, we first explore the thermodynamics of bit erasure in a simple model: a particle with position and momentum in a double-well potential $V(x,t)$ and in contact with a heat reservoir at temperature $T$. (Refer to Fig. 1.) An external controller adds or removes energy from a work reservoir to change the form of the potential $V(\cdot,t)$ via a predetermined erasure protocol $\{\left(\beta(t),\delta(t)\right):0\leq t\leq\tau\}$. $\beta(t)$ and $\delta(t)$ change one at a time piecewise-linearly through four protocol substages: (1) drop barrier, (2) tilt, (3) raise barrier, and (4) untilt. (See SM VI.) The system starts at time $t=0$ in the equilibrium distribution for a double-well $V(x,0)$ at temperature $T$. Being equiprobable, the informational states associated with each of the two wells thus contain $1$ bit of information [25].

The default potential, $V(\cdot,0)=V(\cdot,\tau)$, has two symmetric wells separated by a barrier. Following common practice we call the two wells, from negative to positive position, the Left ($L$) and Right ($R$) informational states, respectively.

The erasure protocol is designed so that the particle ends in the $R$ state with high probability, regardless of its initial state. Conducting our simulation $3.5\times 10^{6}$ times, $96.2\%$ of the particles were successfully erased into the $R$ state. Thus, as measured by the Shannon entropy, the initial $1$ bit of information was reduced to $0.231$ bits. Note that we chose the protocol to give partially inaccurate erasure in order to illustrate our main results on diagnosing success and failure.

At all other times $t$, $V(\cdot,t)$ has either one or two local minima, naturally defining metastable regions for a particle to be constrained and gradually evolve towards local equilibrium. We therefore define the informational states at time $0\leq t\leq\tau$ to be the metastable regions, labeling them $R$ and, if two exist, $L$ — from most positive to negative in position.

Since the protocol is composed of four simple substages, we coarse-grain the system’s response by its activity during each substage at the level of its informational state. Specifically, for each substage, we assign one of three substage trajectory classes: the system (i) was always in the $R$ state, (ii) was always in the $L$ state, or (iii) spent time in each. Sometimes there is only one informational state and so the latter two classes are not achievable for all substages.

We then focus on a single mesoscopic observable—the thermodynamic work expended during erasure. An individual realization generates a trajectory of system microstates, with $W(t,t^{\prime})$ being the work done on the system between times $0\leq t<t^{\prime}\leq\tau$; see SM VI. Let $W_{s}=W(t_{s-1},t_{s})$ denote the work generated during substage $s$ and $C_{s}$ the substage trajectory class. Figure 1 (Outer plot sequence) shows the corresponding substage work distributions $\Pr(W_{s},C_{s})$ obtained from our simulations. (See SM VII.)

The drop-barrier and tilt substage work distributions are rather simple, being narrow and unimodal. The raise-barrier distributions have some asymmetry, but are also similarly simple. However, the untilt work distributions (farthest right in Fig. 1) exhibit unusual features that are significant for understanding the intricacies of erasure. Trajectories that spend all of the untilt substage in either the $R$ state or $L$ state form peaks at the most positive (red) and negative (orange) work values, respectively. This is because the $R$-state well is always increasing in potential energy while the $L$-state well is always decreasing during untilt. In contrast, the other trajectories contribute a log-linear ramp of work values (blue) dependent on the time spent in each. The ramp’s positive slope signifies that more time is typically spent in the $R$ state.

Looking at the total work $W_{\text{total}}=W(0,\tau)$ generated for each trajectory over the course of the entire erasure protocol, we observe the strikingly complex and structured distribution $\Pr(W_{\text{total}})$ shown in Fig. 2(Rear). There are two clear peaks at the most positive and negative work values separated by a ramp. This highly structured work distribution, generated by bit erasure, contrasts sharply with the unimodal work distributions common in previous studies; see, for example, Fig. 2(inset) for the work distribution generated by a thermodynamically-driven simple harmonic oscillator translated in space or Fig. 2 in Ref. [12].

We can understand the mechanisms behind this structure when decomposing Fig. 2 (Rear)’s total work distribution under the untilt substage trajectory classes $C_{4}$. We label trajectories that spend all of the untilting substage in the $R$ state Success since, via the previous substages, they reach the intended $R$ state by the untilting substage and remain there until the protocol’s end. Similarly, trajectories that spend all of the untilt substage in the $L$ state are labeled Fail. The remaining trajectories are labeled Transitional, since they transition between the two informational states during untilt, potentially succeeding or failing to end in the $R$ state. Figure 2 (Three front plots) shows the work distribution for each of these three trajectory classes. Together they recover the total work distribution over all trajectories shown in Fig. 2(Rear). Though, now the thermodynamic contributions to the total from the functionally-distinct component trajectories are made apparent.

Exploring the mesoscale dynamics of erasure revealed signatures of a “thermodynamics” for each trajectory that is closely associated with successful or failed information processing. We now introduce the underlying fluctuation theory from which the trajectory thermodynamics follow. Key to this is comparing system behaviors in both forward and reverse time [8, 9, 10, 11, 12, 13, 14]. (See SM III and IV.)

This suite of trajectory-class fluctuation theorems (TCFTs) applies to arbitrary classes of system microstate trajectories obtainable during a thermodynamic transformation. Importantly, they interpolate between Jarzynki’s equality [10] and Crooks’ detailed fluctuation theorem [12], as the trajectory class varies. This lower bounds the average work $\langle W\rangle_{C}$ over any measurable subset $C$ of the ensemble of system microstate trajectories ${\overrightarrow{\mathcal{Z}}}$, where $W$ is the total work for a trajectory:

[TABLE]

with $\Delta F$ the change in equilibrium free energy over the protocol, $\mathcal{P}(C)$ the probability of realizing the class $C$ during the protocol, and $\mathcal{R}(C^{\texttt{R}})$ the probability of obtaining the time reverse of class $C$ under the time-reverse protocol. $k_{\text{B}}$ is Boltzmann’s constant.

The TCFTs lead to several consequences. First, well-formulated trajectory classes allow accurate estimates of the works for their trajectories, even with limited knowledge of system response under the protocol and its time reverse. Second, they strictly and more strongly bound the average work over all trajectories compared to the equilibrium free energy change $\Delta F$. Third, they provide a new expression for obtaining equilibrium free-energy changes:

[TABLE]

Remarkably, this only requires statistics for a particular class $C$ and its reverse $C^{\texttt{R}}$ to produce the system’s free energy change. Since rare microstate trajectories may generate sufficiently negative works that dominate the average exponential work, this leads to a substantial statistical advantage over direct use of Jaryznski’s equality $\Delta F=-k_{\text{B}}T\ln\langle e^{-W/k_{\text{B}}T}\rangle_{\overrightarrow{\mathcal{Z}}}$ for estimating free energies [26].

To explore these predictions, we selected a superconducting flux qubit composed of paired Josephson junctions (Fig. 3(A)), resulting in a double-well nonlinear potential that supports information storage and processing (Fig. 3(B)). SM VIII.1 explains the physics underlying their nonlinear equations of motion, comparing the similarities and differences with our model’s idealized Langevin dynamics.

Despite control protocols for double-well potentials that perform accurate and efficient bit erasure [27], we run the flux qubit in a mode that yields imperfect erasure (Fig. 3(C)). As with the simulations, our intention is to illustrate how trajectory classes and the TCFTs can be used to diagnose and interpret success and failure in microscopic information processing using only mesoscopic measurements of work.

Interplay between the geometric, linear magnetic inductances and the nonlinear Josephson inductances gives rise to a potential landscape that can be controlled with external bias fluxes. It is natural to call the $\phi_{\text{x}}$ and $\phi_{\text{xdc}}$ fluxes, threading the differential mode and the small SQUID loop, respectively, the tilt and barrier controls. (See (Fig. 3(A) caption.) SM VIII presents a derivation of the flux qubit potential and details its calibration. All experiments presented here were carried out at a temperature of $500$ mK.

To execute an erasure protocol, we first choose an information-storage state with a tall barrier and two equal-depth wells. The two-dimensional potential for this at the calibrated device parameters is depicted in Fig. 3(B). We implement the bit erasure protocol as a time-domain deformation imposed by the two control fluxes that starts and ends at the storage configuration. The amplitudes of the control waveforms in reduced units are small; see Fig. 3(C). Hence, the microscopic energetics change linearly as a function of the control fluxes.

We use a local dc-SQUID magnetometer to continuously monitor the trapped flux state in the device—Readout 1 in Fig. 3(A). The digitized signal has a rise time of $100~{}\mu$s, after which the two logical states are discriminated virtually without error. A typical magnetometer trace $V(t)$ acquired during the execution of the erasure protocol is shown in Fig. 3(C). We operate the magnetometer with a low-amplitude AC current bias at 10 MHz to avoid an increase in the effective temperature during continuous readout of the flux state due to wideband electromagnetic interference.

To collect work statistics, we repeat the erasure protocol $10^{5}$ times. We identify the logical-state transitions from the magnetometer traces as zero-crossings, recording the direction $\delta_{i}$—sign convention: +1 ($-$1) for a L-to-R (R-to-L) transition—and the time $t_{i}$ relative to the start of the protocol. We evaluate a single-shot work estimate $W=\sum_{i}\delta_{i}U_{LR}(t_{i})$, where $U_{LR}(t)=U_{R}(t)-U_{L}(t)$ is the biasing of the potential minima at time $t_{i}$. Making use of the linearity of the system energetics and the choice of offsets and compensation coefficients, we find $U_{LR}(t)=A\left(\phi_{x}(t)-\phi_{x}(0)\right)$, with the coefficient $A=210\text{K}\times k_{\text{B}}$ evaluated from the calibrated potential. The above work estimate based on the logical-state transitions is an accurate estimate of the true microscopic work assuming that the timescales for the state transitions and for changes in the control parameters are much slower than the intra-well equilibration. (See SM VI.)

The total work distribution estimated from the flux qubit experiments is shown as the rear-most distribution in Fig. 3(D). Using the previous microstate trajectory partitioning into the Success, Fail, and Transitional trajectory classes reveals a decomposition of the total work distribution given by Fig. 3(D)(Three front panels). The close similarity with our simulations (Fig. 2) is notable. Especially so, given the rather substantial differences between the simulated system (idealized double-well potential and thermal noise, exactly one-dimension system, …) and the experimental system (complex potential in two dimensions, nonideal fluctuations, …). A priori it is not clear that the theoretical predictions of the informational classes should apply so directly and immediately to the real-world qubit. In point of fact, these differences serve to emphasize the descriptive power of the mesoscopic work fluctuation theorems: despite substantial differences in system detail they successfully diagnose the information-processing classes of microscopic trajectories. This robustness will be especially helpful in monitoring thermodynamic computing in biological systems, where, in many cases, information-bearing degrees-of-freedom cannot be precisely modeled.

We experimentally demonstrated that work fluctuations generated by information engines are highly structured. Nonetheless, they strictly obeyed a suite of time-reversal symmetries—the trajectory-class fluctuation theorems introduced here. The latter are direct signatures of how a system’s informational states evolve and they identify functional and nonfunctional microscopic trajectory bundles. We showed that the trajectory-class fluctuation theorems naturally interpolate between Jarzynski’s integral and Crooks’ detailed fluctuation theorems, providing a unified diagnostic probe of nonequilibrium thermodynamic transformations that support information processing.

Using them we gave a detailed mechanistic analysis of the thermodynamics of the now-common example of erasing a bit of information as an external protocol manipulated a stochastic particle in a double-well potential (simulation) and the stochastic state of a flux qubit (experiment). To give insight into the new level of mechanistic analysis possible, we briefly discussed the untilt trajectory-class partitioning. Though ignoring other protocol stages, this was sufficient to capture the basic trajectory classes that generate the overall work distribution’s features. Partitioning on informational-state occupation times during barrier raising and untilting—an alternative used in follow-on studies—yields an even more incisive decomposition of the work distributions and diagnosis of informational functioning. The corresponding bounds on thermodynamic resources obtained via the TCFTs also improve on current estimation methods. The net result is that trajectory-class analysis can be readily applied to debug thermodynamic computing by engineered or biological systems.

Supplementary Materials

Materials and Methods: Derivation of trajectory-class fluctuation theorems, further discussion and interpretation, and experiment implementation, calibration, and work estimation methods.

**Supplementary Materials

*Harnessing Fluctuations in Thermodynamic Computing

via

Time-Reversal Symmetries

One Sentence Summary:

Distinct distributions of thermodynamic work identify signatures of successful and failed information processing in a microscale flux qubit

Gregory Wimsatt, Olli-Pentti Saira, Alexander B. Boyd, Matthew H. Matheny,

Siyuan Han, Michael L. Roukes, and James P. Crutchfield **

Materials and Methods

The following presents details, derivations, and explanations for the theoretical claims and experimental results. First, we give a synopsis of recently developed principles of thermodynamic computing. Then, we introduce the model microscopic stochastic thermodynamical system, including its equations of motion and its physical calibration. The following sections explain how to use and interpret the trajectory-class fluctuation theorems and provide their derivation. Their practical application to work estimation for trajectories and classes, in light of the alternative kinds of work (inclusive and exclusive) and methods for experimental estimation, is laid out. A brief commentary on the substage work distributions then follows. Finally, we turn to describe the flux qubit, its equations of motion, implementation, calibration, and measurement.

I Principles of Thermodynamic Computing: A recent synopsis

A number of closely-related thermodynamic costs of computing have been identified, above and beyond the house-keeping heat that maintains a system’s overall nonequilibrium dynamical state. First, there is the information-processing Second Law ****[28]**** that extends Landauer’s original bound on erasure ****[2]**** to dissipation in general computing and properly highlights the central role of information generation measured via the physical substrate’s dynamical Kolmogorov-Sinai entropy. It specifies the minimum amount of energy that must be supplied to drive a given amount of computation forward. Second, when coupling thermodynamic systems together, even a single system and a complex environment, there are transient costs as the system synchronizes to, predicts, and then adapts to errors in its environment ****[29, 30, 31]****. Third, the very modularity of a system’s organization imposes thermodynamic costs ****[32]****. Fourth, since computing is necessarily far out of equilibrium and nonsteady state, there are costs due to driving transitions between information-storage states ****[33]****. Fifth, there are costs to generating randomness ****[34]****, which is itself a widely useful resource. Finally, by way of harnessing these principles, new strategies for optimally controlling nonequilibrium transformations have been introduced ****[35, 36, 37, 27]****.

II Microscopic Stochastic Thermodynamical System

For concreteness, we concentrate on a one-dimensional system: a particle with position and momentum in an external potential $V(x,t)$ and in contact with a heat reservoir at temperature $T$. An external controller adds or removes energy from a work reservoir to change the form of the potential $V(\cdot,t)$ via a predetermined erasure protocol $\{\left(\beta(t),\delta(t)\right):0\leq t\leq\tau\}$. (See Supplementary Materials (SM) VI for details on the alternative definitions of work.) The potential takes the form:

[TABLE]

with constants $a,b_{0},c_{0}>0$. During the erasure protocol, $\beta(t)$ and $\delta(t)$ change one at a time piecewise-linearly through four protocol substages: (1) drop barrier, (2) tilt, (3) raise barrier, and (4) untilt, as shown in Table S1. The system starts at time $t=0$ in the equilibrium distribution for a double-well $V(x,0)$ at temperature $T$. Being equiprobable, the informational states associated with each of the two wells thus contain $1$ bit of information ****[25]****. The effect of the control protocol on the system potential and system response is graphically displayed in Fig. 1.

We model the erasure physical information processing with underdamped Langevin dynamics:

[TABLE]

where $k_{\text{B}}$ is Boltzmann’s constant, $\gamma$ is the coupling between the heat reservoir and system, $m$ is the particle’s mass, and $r(t)$ is a memoryless Gaussian random variable with $\langle r(t)\rangle=0$ and $\langle r(t)r(t^{\prime})\rangle=\delta(t-t^{\prime})$.

For comparison to experiment, we simulated erasure with the following parameters, sufficient to fully specify the dynamics: $\gamma\tau/m=500$, $2k_{\text{B}}T\tau^{2}a/(mb_{0})=2.5\times 10^{5}$, $b_{0}^{2}/(4ak_{\text{B}}T)=7$, and $\sqrt{8a/b_{0}^{3}}c_{0}=0.4$. The resulting potential, snapshotted at times during the erasure substages, is shown in Fig. 1(Inner plot sequence).

Reliable information processing dictates that we set time scales so that the system temporarily, but stably, stores information. To support metastable-quasistatic behavior at all times the relaxation rates of the informational states are much faster than the rate of change of the potential, keeping the system near metastable equilibrium throughout. The entropy production for such protocols tends to be minimized.

III Trajectory-Class Fluctuation Theorems: Use and Interpretation

Here, we describe the trajectory-class fluctuation theorems, explaining several of their possible implications and exploring their application to both the simulations and flux qubit experiment. Their derivations are given in the section following.

First, consider a forward process distribution $\mathcal{P}$, defined by the probabilities of the system microstate trajectories ${\overrightarrow{\mathcal{Z}}}$ due to an initial equilibrium microstate distribution evolving forward in time under a control protocol. Then, the reverse process distribution $\mathcal{R}$ is determined by preparing the system in equilibrium in the final protocol configuration and running the reverse protocol. The reverse protocol is the original protocol conducted in reverse order but also with objects that are odd under time reversal, like magnetic fields, negated. The time-reversal of a trajectory ${\overrightarrow{z}}=(z_{0},z_{1},\ldots z_{\ell})$ is ${\overrightarrow{z}}^{\texttt{R}}=(z_{\ell}^{\texttt{R}},\ldots,z_{1}^{\texttt{R}},z_{0}^{\texttt{R}})$, where $z_{i}^{\texttt{R}}=-z_{i}$ if $z_{i}$ is odd under time-reversal (e.g., momentum or spin), otherwise $z_{i}^{\texttt{R}}=z_{i}$. For a measurable subset of trajectories $C\subset{\overrightarrow{\mathcal{Z}}}$, let $\langle\cdot\rangle_{C}$ denote an average over the ensemble of forward process trajectories conditioned on the trajectory class $C$. Let $\mathcal{P}(C)$ and $\mathcal{R}(C^{\texttt{R}})$ denote the probabilities of observing the class $C$ in the forward process and the reverse class $C^{\texttt{R}}=\{{\overrightarrow{z}}^{\texttt{R}}|{\overrightarrow{z}}\in C\}$ in the reverse process, respectively.

We first introduce a trajectory-class fluctuation theorem (TCFT) for the class-averaged exponential work $\langle e^{-W/k_{\text{B}}T}\rangle_{C}$:

[TABLE]

with $\Delta F$ the system equilibrium free energy change. We also introduce a class-averaged work TCFT:

[TABLE]

This employs the Kullback-Liebler divergence $D_{\text{KL}}[~{}\cdot~{}]$ taken between forward and reverse process distributions over all class trajectories ${\overrightarrow{z}}\in C$, conditioned on the forward class $C$ and reverse class $C^{\texttt{R}}$, respectively. If we disregard this divergence, which is nonnegative and would generally be difficult to obtain experimentally, we then find the lower bound $\langle W\rangle_{C}^{\text{min}}$ on the class-averaged work of Eq. (1).

In the limit of class $C$ possessing only a single trajectory, we recover detailed fluctuation theorems as in Ref. ****[12]****. If, however, we take $C$ to be the entire set of trajectories ${\overrightarrow{\mathcal{Z}}}$, we recover integral fluctuation theorems as in Jarzynski’s equality ****[10]****. Thus, the TCFTs are a suite that spans the space of fluctuation theorems between the extreme of the detailed theorems, that require very precise information about an individual trajectory, and the integral theorems, that describe the system’s entire trajectory ensemble. SM IV below provides proofs for both TCFTs.

We can rearrange Eq. (S2) to obtain Eq. (2)—an expression for estimating equilibrium free energy changes:

[TABLE]

Thus, to estimate free energy one sees that statistics are needed for only one particular class and its reverse. Generally, this gives a substantial statistical advantage over direct use of Jaryznski’s equality:

[TABLE]

since rare microstate trajectories may generate negative work values that dominate the average exponential work ****[26]****. The problem is clear in the case of erasure. Recall from Fig. 2(Three front panels) that Fail trajectories generate the most-negative work values. In the limit of higher success-rate protocols that maintain low entropy production, failures generate more and more negative works, leading them to dominate when estimating average exponential works.

In contrast, to efficiently determine the change in equilibrium free energy from Eq. (2), its form indicates that one should choose a class that (i) is common in the forward process, (ii) has a reverse class that is common in the reverse process, and (iii) generates a narrow work distribution. This maximizes the accuracy of statistical estimates for the three factors on the RHS. For example, while the equilibrium free energy change in the case of our erasure protocol is theoretically simple (zero); the Success class fits the criteria.

We can then monitor the class-averaged work in excess of its bound:

[TABLE]

The inequality in Eq. (1) is a refinement of the equilibrium Second Law and therefore the bound $\langle W\rangle_{C}^{\text{min}}$ generally provides a more accurate estimate of the average work of trajectories in a class compared to the equilibrium free energy change $\Delta F$. More precisely, as we will see below, an average of the excess $E_{C}$ over all classes $C$ in a partition of trajectories must be smaller than the dissipated work $\langle W\rangle-\Delta F$. For trajectory classes with narrow work distributions, this can be a significant improvement. We can see this by Taylor expanding the LHS of Eq. (S2) about the mean dimensionless work $\langle W/k_{\text{B}}T\rangle_{C}$. This shows that Eq. (1) becomes an equality when the variance and higher moments vanish. SM V below delves more into moment approximations. In any case, trajectory classes with narrow work distributions have small excess works $E_{C}$.

To estimate $\mathcal{R}(C^{\texttt{R}})$, we ran $3.5\times 10^{6}$ simulations of the reverse process. Table S2 shows that the Success and Fail classes have small excesses and, as seen in Fig. 2(Three front panels), these classes indeed have narrow work distributions. Elsewhere we explore these and additional partition schemes, finding that the Transitional trajectories can be further partitioned to yield narrow work distributions so that all trajectory classes have small excesses $E_{C}$. In short, this demonstrates how well-formulated trajectory classes allow accurate estimates on the works for all trajectories.

To measure the efficacy of a given partition $Q$ of trajectories into classes, we ask what the ensemble-average of class-average excess works is:

[TABLE]

with $\langle W\rangle_{Q}^{\text{min}}=\sum_{C\in Q}\mathcal{P}(C)\langle W\rangle_{C}^{\text{min}}$.

From Eq. (1), we see that $\langle W\rangle_{Q}^{\text{min}}$ is the coarse-grained lower bound on ensemble-average dissipation from Ref. ****[38]****:

[TABLE]

where $D_{\text{KL}}\left[~{}\cdot~{}\right]$ is the Kullback-Liebler divergence between forward and reverse process distributions over the trajectory classes $C\in Q$. Since Kullback-Liebler divergences are nonnegative, such a bound always provides an improvement over the equilibrium Second Law. Table S2 shows both $\langle W\rangle_{Q}^{\text{min}}$ and $E_{Q}$ for the trivial partition $\{{\overrightarrow{\mathcal{Z}}}\}$, our three-class partition, labeled Untilt-Centric I, and the improved partition described in follow-on work, labeled Untilt-Centric II. In this case, the latter two also provide an improvement on the nonequilibrium Second Law which, assuming metastable starting and ending distributions, provides a lower bound on the average work equal to $0.533$, the change in nonequilibrium free energy.

We can appeal to Landauer’s erasure bound—$k_{B}T\ln 2\approx 0.693\,k_{B}T$—to calibrate the excesses $E_{C}$ and $E_{Q}$. We see for the simulation data that our three-class partition Untilt-Centric I provides class-average work bounds that, on average, are only about $11\%$ of $k_{B}T\ln 2$ from the actual class-average works. The more refined Untilt-Centric II partition reduces this excess to about $5\%$ while the trivial partition fails by about $91\%$ of $k_{B}T\ln 2$.

The experimental data matches these results, with the largest discrepancy occurring for the class-average excess for the Fail class. This is not wholly coincidence, since we determined the parameters of the experimental protocol by adjusting the parameters of a simple two-state Markov simulation to obtain a work distribution similar to that obtained by our Langevin simulations of the Duffing potential system described in the main text. However, it is interesting that this was sufficient to provide matches in both the clean decomposition of the total work distribution by trajectory class and the quantities of Table S2.

We also recover the equality of Ref. ****[38]**** for the ensemble-average work by averaging Eq. (S3) over each class:

[TABLE]

which of course is lower bounded by $\langle W\rangle_{Q}^{\text{min}}$.

These results suggest the criterion for optimal trajectory partitions: Select a partition sufficiently refined to yield tight bounds on class-average works, but no finer. Machine learning methods for model order selection will provide a basis for a natural classification scheme for trajectories that captures all relevant thermodynamics and information processing.

By changing our forward and reverse processes $\mathcal{P}$ and $\mathcal{R}$ to begin in system microstate distributions other than equilibrium, a yet-broader class of TCFTs emerge. We can then find analogous results for heats and comparisons with works and nonequilibrium free-energy changes. We explore these in depth elsewhere.

IV TCFT Derivations

We now present derivations for the two TCFTs introduced in Eqs. (S2) and (S3).

Assume that the system dynamics is described by a Hamiltonian specified in part by an external control protocol, as well as by a weak coupling to a thermal environment that induces steady relaxation to canonical equilibrium.

Start the system in equilibrium distribution $\pi_{0}$ for Hamiltonian $\mathcal{H}_{0}$ and run a protocol until time $\tau$, causing the system Hamiltonian to evolve to $\mathcal{H}_{\tau}$. If we then hold the Hamiltonian at $\mathcal{H}_{\tau}$ for a long time, the system relaxes into the equilibrium distribution $\pi_{\tau}$. The system’s ensemble entropy change from $t=0$ to $t=\infty$ is then:

[TABLE]

The trajectory-wise system entropy difference is defined to be:

[TABLE]

where ${z}_{0}$ and ${z}_{\tau}$ are the initial and final microstates of system microstate trajectory ${\overrightarrow{z}}$, respectively. Averaged over all trajectories ${\overrightarrow{z}}\in{\overrightarrow{\mathcal{Z}}}$, this then becomes the ensemble entropy change.

Let $p({\overrightarrow{z}}|{z}_{0})$ denote the probability of obtaining system microstate trajectory ${\overrightarrow{z}}$ via the protocol conditioned on starting the system in state ${z}_{0}={\overrightarrow{z}}(0)$.

Now, start the system Hamiltonian at $\mathcal{H}_{\tau}$ and run the reverse protocol, ending the Hamiltonian at $\mathcal{H}_{0}$. We then obtain the trajectory ${\overrightarrow{z}}$ with a different conditional probability: $r({\overrightarrow{z}}|{z}_{0})$.

Assuming microscopic reversibility and given a system trajectory ${\overrightarrow{z}}$, the change in the heat bath’s entropy is:

[TABLE]

where $\beta=1/k_{\text{B}}T$, $Q({\overrightarrow{z}})$ is the net energy that flows out of the heat bath into the system given the trajectory ${\overrightarrow{z}}$, and $(\cdot)^{\texttt{R}}$ denotes time-reversal. This holds for systems with strictly finite energies and Markov dynamics that induce the equilibrium distribution when control parameters are held fixed ****[39]****. Both our simulated Duffing potential system and flux qubit obey these requirements at sufficiently short time scales. Then we can express the total trajectory-wise change in entropy production due to a trajectory ${\overrightarrow{z}}$ as the sum of system and heat reservoir entropy changes:

[TABLE]

Since $\pi_{t}({z})=e^{-\beta(\mathcal{H}_{t}({z})-F_{t})}$, with $F_{t}$ the system’s equilibrium free energy at time $t$, we can write:

[TABLE]

Using Eq. (S5), we also have:

[TABLE]

Combining, we obtain a detailed fluctuation theorem:

[TABLE]

From here, we derive our first TCFT by integrating each side of Eq. (S6) over all trajectories ${\overrightarrow{z}}$ in a measurable set $C\subset{\overrightarrow{\mathcal{Z}}}$. Starting with the LHS and recalling the Iverson bracket $[\cdot]$, which is $1$ when the interior expression is true and [math] when false, we have:

[TABLE]

The first three steps used the unity of the Jacobian in reversing a microstate, the definition $C^{\texttt{R}}=\{{\overrightarrow{z}}^{\texttt{R}}|{\overrightarrow{z}}\in C\}$, and swapping all instances of ${\overrightarrow{z}}^{\texttt{R}}$ with ${\overrightarrow{z}}$. Integrating the RHS of Eq. (S6) then gives:

[TABLE]

Combining, we have our first TCFT, Eq. (S2).

To obtain the second TCFT, we first change the form of Eq. (S6):

[TABLE]

Then we calculate the class-average. The equilibrium free energy change is unaffected while the rightmost term becomes:

[TABLE]

which gives Eq. (S3)’s TCFT.

V Class-Averaged Work Approximation for Narrow Distributions

Here, we demonstrate that the class-averaged work $\langle W\rangle_{C}$ approaches its bound $\langle W\rangle_{C}^{\text{min}}$ when the variance and higher moments of the class’ distribution of works vanish. One concludes that $\langle W\rangle_{C}^{\text{min}}$ is a good approximation for $\langle W\rangle_{C}$ when the class’ work distribution is narrow.

We first express the LHS of Eq. (S2) in terms of the unitless distance of work from its class-average:

[TABLE]

with $x=\beta(W-\langle W\rangle_{C})$. Then, we Taylor expand the exponential inside the class-average:

[TABLE]

with $a=\sum_{n=2}^{\infty}\frac{(-1)^{n}}{n!}\langle x^{n}\rangle_{C}$. Equation (S2) then gives:

[TABLE]

Since $e^{-x}$ is convex,

[TABLE]

so $a\geq 0$. Then:

[TABLE]

The second line becomes an equality when $a$ goes to zero, which occurs as the variance and higher moments vanish.

VI Work Definitions and Experimental Estimation

Properly estimating the required works and devolved heats from experimental devices undergoing cyclic control protocols requires explicitly and consistently accounting for energy and information flows between the system, its environment, and the controlling laboratory apparatus. To this end, we construct a model Hamiltonian universe for common processes involving small systems interacting with laboratory apparatus and a thermal environment. After deriving key equalities for two definitions of work, the inclusive and exclusive works, we define a method of approximating them in appropriate cyclic protocols.

VI.1 The Model Universe and Hamiltonian

To study a small system that exchanges energy with its environment in the forms of heat and work, we introduce a model universe: a system of interest, a heat bath, and a lab (laboratory apparatus) that controls the system and derives any needed energy from a work reservoir. The system directly interacts with both the heat bath and the lab, but the heat bath and lab are not directly coupled.

We assume that a Hamiltonian $\mathcal{H}$ describes the universe’s evolution and that there is a set of generalized coordinates which can be sensibly partitioned into those for the system, heat bath, and lab. Then, we decompose the universe Hamiltonian into the following form:

[TABLE]

where $s$, $b$, and $l$ denote both the generalized coordinates and conjugate momenta for the system, bath, and lab, respectively. For any universe Hamiltonian $\mathcal{H}$, there can be many choices for this decomposition.

We also define the system Hamiltonian $\mathcal{H}^{\prime}$ as the three components that depend on the system coordinates:

[TABLE]

First, consider the subset of lab coordinates $l$ for which $h_{\text{S,L}}$ has nontrivial dependence. These so-called protocol parameters $\lambda$ are often simple and much fewer than the entire set of $l$. We often assume that we have total control of their evolution. More precisely, under an appropriate preparation for the lab at time $t=0$, a specific trajectory for the protocol parameters $\{\lambda(t)\}_{t}$ for $0\leq t\leq\tau$ is guaranteed for all preparations of the heat bath and system coordinates. We refer to the parameter trajectory as the protocol.

Suppose the heat-bath degrees of freedom that interact with the system change much faster than the system’s. We can assume that the system response to the bath resembles Brownian motion. On the time scale of changes in the system coordinates, then, we ignore the system-bath interaction term $h_{\text{S,B}}$ in writing the system Hamiltonian $\mathcal{H}^{\prime}$:

[TABLE]

The latter decomposition into kinetic energy $T$ and potential energy $V$ can be used to write Langevin equations of motion for the system. Furthermore, if the heat bath has a relaxation time sufficiently short that it is roughly in equilibrium at all times with fixed temperature, then its influence on the system will be memoryless.

VI.2 Inclusive and Exclusive Works and Heats

The basic scenario for executing a protocol is as follows. The universe coordinates begin according to a given initial distribution $\Pr(s)$ at time $t=0$ and they evolve in isolation until $t=\tau$. As above, we assume that a well-defined protocol $\{\lambda_{t}\}_{t}$ emerges due to our preparation of the lab coordinates.

We label all energy exchanged between the system and lab as work and all energy exchanged between the system and heat bath as heat. Since the lab is directly coupled only to the system, the work they exchange is given by the change in energy of the lab’s work reservoir. Similarly, since the heat bath is directly coupled only to the system, the heat exchanged is given by the change in the heat bath’s energy.

Note that this requires choices as to what constitute the energies of the three universe subsystems. While $H_{\text{B}}$, $H_{\text{S}}$, and $H_{\text{L}}$ define energies for the heat bath, system, and work reservoir, respectively, what of $h_{\text{S,L}}$ and $h_{\text{S,B}}$? If all subsystems were macroscopic, these interaction terms would be negligible. While it may be desirable to assume that the system is only weakly coupled to the heat bath—so that $h_{\text{S,B}}$ can be ignored—$h_{\text{S,L}}$** can be significant in many important small systems.**

And so, in general, we define the system energy to be $H_{\text{S}}$ plus any portions of $h_{\text{S,L}}$ and $h_{\text{S,B}}$. Then the work reservoir energy is $H_{\text{L}}$ plus the rest of $h_{\text{S,L}}$, while the heat bath energy is $H_{\text{B}}$ plus the rest of $h_{\text{S,B}}$. To make these distinctions clear we label two types of works, each corresponding to the two extremes for allocation of $h_{\text{S,L}}$ between the system and work reservoir: the inclusive work $W$ and the exclusive work $W_{0}$ ****[40]****. Specifically:

[TABLE]

We can similarly define the inclusive heat $Q$ and exclusive heat $Q_{0}$ depending on how we allocate $h_{\text{S,B}}$ between the system and heat bath:

[TABLE]

The inclusive work corresponds to fully including $h_{\text{S,L}}$ in the system energy, while the exclusive work corresponds to excluding it. Inclusive and exclusive heat correspond similarly with respect to $h_{\text{S,B}}$.

There is a key relation between the inclusive and exclusive works:

[TABLE]

That is, the inclusive work for an interval of time equals the sum of the exclusive work and the change in the system-lab interaction term $h_{\text{S,L}}$.

**In the above expressions, calculating the rate of change of a work or heat requires the time derivative of one or more of $H_{\text{L}}$ and $H_{\text{B}}$. This can be problematic. Fortunately, there are alternate forms that are amenable. One can show that the inclusive work rate is given by: **

[TABLE]

This is a more common definition for the work rate in small-system nonequilibrium thermodynamics. And, it allows the work to be calculated as:

[TABLE]

The exclusive work $W_{0}$ has a corresponding form:

[TABLE]

For the case where $h_{\text{S,L}}$ is a scalar potential for $s$, this is the product of the corresponding force with velocity. This makes the exclusive work equal to a familiar mechanics definition of work as the integral of the dot product of force and displacement:

[TABLE]

In this way, we write the inclusive and exclusive work rates in terms of the rates of change of the system and work-reservoir interaction term $h_{\text{S,L}}$ with respect to either the system or work reservoir coordinates.

VI.3 Approximating Inclusive Work Experimentally

For the flux qubit experimental system investigated here, we assume the following:

[TABLE]

That is, as far as the flux qubit and work reservoir are concerned, the only relevant energies at least partially ascribable to the flux qubit are its kinetic energy and the potential energy with the work reservoir. $h_{\text{S,L}}$ must then capture the change in the potential $V$ due to changes in the protocol parameters. We could simply define $h_{\text{S,L}}(s,\lambda)=V(s,\lambda)$ so that $H_{\text{S}}(s)=T(s)$. However, it is more useful to allocate the initial potential energy to $H_{\text{S}}$. That is:

[TABLE]

For cyclic protocols where $V(\cdot,\lambda_{0})=V(\cdot,\lambda_{\tau})$, such as in our erasure operation, $h_{\text{S,L}}(s(t),\lambda(t))$ vanishes for all trajectories at $t=0,\tau$. By Eq. (S7) we then have the useful equality $W=W_{0}$ between inclusive and exclusive works taken over the entire protocol.

Estimating $W$ for a system trajectory is then equivalent to estimating $W_{0}$ for the cyclic protocols we consider. In the flux qubit, the form of $h_{\text{S,L}}$ is known and the specific protocol $\{\lambda_{t}\}_{t\in[0,\tau]}$ is known. Unfortunately, we lack sufficient information about its instantaneous state $s$ at all times, since the device’s physics precludes precise measurements of system flux $\phi$—the relevant part of $s$ for determining the potential $h_{\text{S,L}}$. Instead, we do have reliable measurement of large and stable changes in the flux $\phi$. This specifically monitors when the system moves between wells in a double-well potential $V(\cdot,\lambda(t))$, if the rate of transition between wells is sufficiently slow.

And so, we can use information about the flux $\phi$ to approximate the exclusive work contribution at each moment in time. Then, adding up these contributions yields an approximation to the total exclusive work $W_{0}$ over the entire protocol and therefore of the inclusive work $W$ over the entire protocol. Note that the protocols used here maintain two wells at all times for the system flux $\phi$. We develop the approximation in two steps.

VI.3.1 First-Order Approximation

We first partition the potential in flux space into three segments. Two segments constitute the wells for the flux in which that state spends all its time except for very brief transitions between wells. Then, the third segment connects the two wells, capturing the dynamics arising from crossing the barrier that separates them.

We require that the partitioning allows the following two approximations. First, the particle spends negligible total duration in between the two wells. Second, the wells do not change shape over the protocol, but instead simply raise or lower in potential at different times, if they change at all. This means that the shape of the system-lab interaction term $h_{\text{S,L}}(\cdot,\lambda(t))$ at any time $t$ is very simple in the two wells—flat.

The result is that the exclusive work over any time duration is easily calculable from the experimental data. During times when the flux remains in a well, the exclusive work must be zero, since $h_{\text{S,L}}$ does not change with $s$. During a transition, the shape of $h_{\text{S,L}}$ does not change due to the first approximation. Then, the exclusive work is the difference in heights of the two wells as measured by $h_{\text{S,L}}$:

[TABLE]

where $\lambda(t)$ is the protocol parameter setting at any time during the transition and $w_{0}$ and $w_{1}$ are arbitrary flux values in the starting and ending wells, respectively.

Thus, the total inclusive work $W$ over the protocol for a trajectory is simply the sum of the jump contributions above for each transition.

VI.3.2 Second-Order Approximation

In point of fact, the potential wells do change shape. Fortunately, our method for calculating the inclusive work over the protocol remains valid under weaker constraints on the protocol.

We first require the protocol to maintain two metastable regions, the informational states, at all times; each possessing a unique local potential minimum continuously in time. We denote the flux value at the potential minima of informational state $i$ at time $t$ as $\phi_{i}^{t}$. The protocol must also evolve slowly enough so that the potential landscape changes slowly compared to the system’s relaxation rate in each metastable region. Both of these criteria are met by our erasure protocol.

Consider a short duration $\Delta t$ during which the potential $V(\cdot,t)$ changes little but long enough compared to the relaxation rates of the informational states. Consider two cases: either the system crosses the barrier between the two informational states during this time or it remains in one informational state.

First, suppose that the system transitions from one informational state $i$ to the other $j$. Denote the system flux at the beginning of the transition as $\phi^{t}$ and at the end as $\phi^{t+\Delta t}$. By Eq. (S7), the exclusive work contribution $\Delta W_{0}^{\text{trans}}$ is the difference of the inclusive work contribution and the change in system-lab interaction term $\Delta h_{\text{S,L}}$. The change $\Delta h_{\text{S,L}}$ can itself be broken down into two terms, one for the difference in $h_{\text{S,L}}$ between the informational-state minima and the other for the change in $h_{\text{S,L}}$ local to the respective minima. In other words:

[TABLE]

where $\Delta m^{t}$, Eq. (S12)’s first term, is the change in $h_{\text{S,L}}$ at the informational-state minima and $\Delta l^{t}$, Eq. (S12)’s second term is the change in $h_{\text{S,L}}$ of the system with respect to the informational-state minima. Our protocol ensures that the total number of transitions is so small and the time durations so narrow that we can ignore the total contributions of inclusive works $\Delta W^{\text{trans}}$ due to these transition durations. Then, we approximate the exclusive work contribution during a transition via:

[TABLE]

Suppose, now, that the system remains in one informational state $i$ during a time interval $\Delta t$. Since the relaxation rate is fast compared to the duration $\Delta t$, we assume that the system visits all microstates in the informational state roughly in proportion to the local equilibrium distribution. Then, the inclusive work contribution $\Delta W^{\text{stay}}$ is approximately independent of the specific system trajectory during this time and, instead, is determined by the time duration and the informational state $i$. If during this time we simultaneously shift the entire potential up by a given amount, we add an inclusive work contribution equal to the potential shift but the system trajectory is unchanged. Thus, the actual inclusive work contribution is equal to an amount due to the change in the system-lab interaction term at the informational-state minimum plus an amount due solely to the change in potential shape at the informational state with respect to its minimum. That is:

[TABLE]

where $\Delta W_{s}$ is the inclusive work contribution due to the change in potential shape at the informational state. Equation (S13) applies equally well here in describing the change in system-lab interaction term. Thus, the exclusive work contribution $\Delta W_{0}^{\text{stay}}$ for this time interval is:

[TABLE]

The result is that we have exclusive work contributions for both durations when the system transitions between informational states and when it remains in one.

To find the total exclusive work over the protocol for a given trajectory we add up the contributions. The sum of all local $h_{\text{S,L}}$ changes $\Delta l^{t}$ over all durations is the net local change in $h_{\text{S,L}}$. Recall, though, that the minima of the informational states begin and end at the same values. And so, the total local change in $h_{\text{S,L}}$ reduces to the absolute change in $h_{\text{S,L}}$. However, since we chose $h_{\text{S,L}}(\cdot,0)=h_{\text{S,L}}(\cdot,\lambda(t))=0$, this must vanish:

[TABLE]

We can now specify our final approximation: At any time $t$, the inclusive work contribution $\Delta W_{s}$ due to the change in potential shape is independent of the informational state. This is reasonable for our erasure protocol since the asymmetric contribution to the change in potential—the tilt—is slight. While it clearly breaks the symmetry of the double-well potential by changing the well heights, it has less effect on the well shapes and even less in making those shapes distinct.

Then, we can assume that the sum of $\Delta W_{s}$ for a trajectory is the same as that staying in one informational state the entire time. Since the protocol is very slow and cyclic, though, a particle that stays in one informational state the entire time must receive approximately zero inclusive work $W$. Given that the sum over all $\Delta W_{s}$ must be equal to $W$ for such a trajectory, it must also be negligible.

Altogether, the total exclusive work is approximately given by the sum over all transitions between informational states of the difference in potential at the informational-state minima:

[TABLE]

To reiterate, since $\Delta h_{\text{S,L}}=0$, this is also the total inclusive work $W(0,\tau)$ for a trajectory over the entire protocol.

VII Substage Work Distributions Commentary

Here, we briefly interpret several features of the substage work distributions observed in Fig. 1(Outer left plots).

The distributions for barrier dropping and tilting are narrow, symmetric peaks; see Fig. 1 (Outer left plots). Barrier raising also has a rather narrow peak, composed primarily of trajectories always in the $R$ state, but also exhibits a bulge toward positive work; see Fig. 1 (Top right). Note that the $L$ state is created mid-way through barrier raising, allowing for trajectories that spend some time in either informational state, but disallowing trajectories that spend all time in the $L$ state. The former induce the positive work bulge toward less negative works, which while notable will not be further explored here.

The substage work distributions for untilting presents the most striking picture; see Fig. 1 (Bottom right). Always-$R$** trajectories induce a large positive work peak (red), always-$L$ trajectories induce a large negative work peak (orange), and all other trajectories induce a ramp between them (blue).**

These features can be directly interpreted by following the locations of the potential minima over time and noting how the shifting potential adds or removes energy from a particle. During barrier dropping, to take one example, the protocol raises both minima by over $7~{}k_{\text{B}}T$, resulting in a narrow, peaked work distribution with a mean near $7~{}k_{\text{B}}T$.

Most interesting is the untilt substage. Since most particles start and then stay in the $R$ state for this substage, a large positive work is probable, due to the rising $R$-state well. However, it is also possible for the system to start in and then get stuck in the $L$-state well, resulting in a large negative work. The final possibility is transitioning between states during untilting, resulting in an intermediate range of less-likely work values. For trajectories that do transition between states during untilting, it is more likely to spend more time in the $R$ state, since it is energetically favored, resulting in the rising probability with increasing work in their work distribution—giving rise to the log-linear ramp in the work distribution.

Note that there are small peaks on each end of this third class’ distributions that require a more nuanced explanation. When a particle crosses a barrier—due to random thermal excitation—the surplus energy may quickly send the particle back to the previous well before it can be dissipated. Such particles then spend almost all of the substage in this first well, generating a work value accordingly. Statistics of the ramp proper are due to particles that have time to locally equilibrate before crossing any barriers.

Follow-on work develops the theory underlying this detailed mechanistic analysis and analyzes similar behavior in all metastable-quasistatic processes.

VIII Flux Qubit Device, Calibration, and Measurement

The benefits of the flux qubit device are several-fold. First, their physics provide a genuine two-degree of freedom dynamics, while other comparable experiments on Maxwellian demons and bit erasure are very high dimensional, only indirectly providing an effectively few-degree of freedom dynamics ****[41, 19, 21]****. Second, they operate at very high frequency and so one readily captures the substantial amounts of data required to accurately estimate rare-fluctuation statistics. Third, they leverage recent advances in manufacturing technology led by efforts in quantum computing. Fourth, being constructed via modern integrated circuit technology they form the basis of a technology that will scale to large, multicomponent circuit devices for more sophisticated thermodynamic computing. And, finally, in the near future flux qubits will facilitate experiments that probe the thermodynamics of the transition to quantum information processing.

At the microscopic level, a fraction of the electrons in a superconducting metal form bosonic Cooper pairs—a quantum-coherent condensate. For designing superconducting electronic circuits, though, one can forgo the microscopic description and work with higher-level phenomena, such as flux quantization and the Josephson relations for weak links. Importantly, the circuit-level degrees of freedom are not coarse-grained quantities, but display a full range of quantum behavior, including quantized excitations, coherent superpositions, and entangled states in such circuits. For our purposes here, however, we run the device so that it exhibits only classical stochastic dynamics, reserving quantum information thermodynamic explorations for the future.

This section lays out the basic physics of the flux qubit device and details of the experimental implementation.

VIII.1 Flux qubit physics

Our experimental information processor is a special type of superconducting quantum interference device (SQUID) with two degrees of freedom—a gradiometric flux qubit or the variable-$I_{c}$** rf SQUID introduced by Ref. [42]. Notably, the energies associated with the motion perpendicular to and along the escape direction differ substantially by about a factor of $12$. Practically, this asymmetry reduces the two-dimensional potential to one dimension. The net result is a device with an effective double-well potential with barriers as low as $\Delta U\sim k_{\text{B}}T$ that operates at frequencies in the GHz range. The potential shape is controlled by fluxes that are readily controlled by currents. SQUID device parameters, used to determine the potential shape and energy scales, were all independently determined.**

The variable-$I_{c}$** rf SQUID replaces the single Josephson junction in a standard rf SQUID with a symmetric dc SQUID with small inductance $\beta_{\text{dc}}=2\pi\ell I_{c0}/\Phi_{0}\ll 1$, where $2\ell$ is the loop inductance, $I_{c0}=i_{c1}+I_{c2}$ is the sum of critical currents of the two junctions, and $\Phi_{0}$ is the flux quantum $h/2e$. This architecture gives a device whose parameters can be accurately measured and that can be selected to exhibit a range of phenomena including thermal activation, macroscopic quantum tunneling, incoherent relaxation, photon-induced transitions, and macroscopic quantum coherence. It also allows us to perform, as we demonstrate, nanoscale thermodynamic computing.**

Its macroscopic dynamical variables are the magnetic flux $\Phi$ through the rf SQUID loop and $\Phi_{\text{dc}}$ through the dc SQUID loop. Based on the resistively-capacitively-shunted junction model of Josephson junctions, in the classical limit the variable-$I_{c}$** rf SQUID’s deterministic equations of motion are [42]:**

[TABLE]

In units of $\Phi_{0}/2\pi$, the 2D potential for the variable-$I_{c}$** rf SQUID is $U(\phi,\phi_{\text{dc}})=U_{0}f(\phi,\phi_{\text{dc}})$ with:**

[TABLE]

where $U_{0}=\Phi^{2}_{0}/(4\pi^{2}L)$. Here, $\gamma=L/(2\ell)$ is the ratio of rf and dc SQUID inductances; $\phi_{x}$ ($\phi_{\text{xdc}}$) is the external flux applied to the rf (dc) SQUID loop; $\phi$ ($\phi_{\text{dc}}$) is the flux enclosed in the rf (dc) SQUID loop; $\beta_{0}=2\pi LI_{c0}/\Phi_{0}$; and $\delta\beta=2\pi L(I_{c2}-I_{c1})/\Phi_{0}$.

For large-amplitude tuning of the external controls, the system response to $\phi_{x}$ ($\phi_{\text{xdc}}$) is $2\pi$ ($4\pi$) periodic. We make use of the global features to accurately determine the coefficients of the potential.

In the experiment, cross-coupling between the barrier and tilt controls was canceled by an affine transformation $(\phi_{\text{x}},\phi_{\text{xdc}})\rightarrow(\phi_{\text{x}}+\alpha\phi_{\text{xdc}},\phi_{\text{xdc}})$, with the coefficient $\alpha$ chosen such that the equilibrium population of the left and right wells was unaffected to first order by the barrier control $\phi_{\text{xdc}}$.

Operating the magnetometer generates wide-band local electromagnetic interference that can affect the dynamics of the flux qubit. A careful study of the back-action indicates that low-amplitude operation of the magnetometer can induce transitions in a manner that corresponds to a shift in the effective tilt and flux controls. Importantly, the effective temperature under magnetometer operations was not elevated from $500$ mK.

The dynamical variable $\phi$ describes the in-phase motion of the two junctions that results in a current circulating in the rf SQUID loop. The dynamical variable $\phi_{\text{dc}}$ describes the out-of-phase motion, resulting in a current circulating in the dc SQUID loop. The shape of the effective potential is completely determined by the dimensionless function $f(\phi,\phi_{\text{dc}})$ and the energy scale of the potential is determined by $U_{0}$. With suitable device parameters and applied fluxes ($\phi_{\text{x}}$** and $\phi_{\text{xdc}}$) one obtains a double-well potential. The barrier height $\Delta U$ separating the two wells is readily adjusted by varying $\phi_{\text{xdc}}$. The effective potential is plotted in Fig. 3(B) with parameters: $\beta_{0}=6.2$, $\gamma=12$ and $\delta_{b}=0.2$.**

VIII.2 Experimental implementation

The junctions were $1\times 1\mu\text{m}^{2}\text{Nb/Al}_{2}\text{O}_{3}/\text{Nb}$ tunnel junctions of very low subgap leakage, typically having a quality factor of $V_{m}\approx 70$ mV at $4.2$ K.

We followed a standard procedure (see, e.g., Ref. ****[42]****) for calibrating the flux qubit parameters. An outline of the steps is given below. A complete description of the measurements is presented elsewhere.

First, by executing wide-range sweeps of the coil currents $I_{\mathrm{tilt}}$ and $I_{\mathrm{barrier}}$, parameter values corresponding to single-valued and bistable potential landscapes are recorded. A linear transformation from $I_{\mathrm{tilt}}$ and $I_{\mathrm{barrier}}$ to $(\phi_{\text{x}},\phi_{\text{xdc}})$ is established by matching the experimental periodicity with the theoretical one $(2\pi,4\pi)$. Linear cross-talk from $I_{\mathrm{barrier}}$ to $I_{\mathrm{tilt}}$ is calibrated by orthogonalizing the global response. Cross-talk from $I_{\mathrm{tilt}}$ to $I_{\mathrm{barrier}}$ can be assumed to be small due to the symmetry of the on-chip flux lines and is taken to be zero.

The parameter values $\beta_{0}=6.2$, $\gamma=12$ and $\delta_{b}=0.2$ are determined by equating the observed extent of hysteresis at $\phi_{\text{xdc}}=0$ and the differential flux response $d\left<\phi\right>/d\phi_{\text{x}}$ at $\phi_{\text{xdc}}=2\pi$ to theoretical predictions. The prefactor $U_{0}=56.3$ K is determined by equating the observed escape energy for inter-well transitions at high temperatures with $k_{\text{B}}T$. The plasma frequency $\omega_{p}=1/\sqrt{LC}=2\pi\times 13.7$ GHz is determined from the observed low-temperature cross-over temperature $T_{c}r=103$ mK to macroscopic quantum tunneling (MQT) dominated dynamics. We obtain an upper bound $Q=\omega_{p}RC<130$ from the coupling to the passive shunt resistor of the magnetometer. Parameter calibration measurements are performed in such a way that the effect of magnetometer back-action is nulled through pulsing of the readout or otherwise minimized. The effective temperature under continuous magnetometer operation was determined by repeating the measurement for escape energy for interwell transitions and comparing the result to that obtained under pulsed magnetometer operation.