{'fraction_non_alphanumeric': 0.036561495501643115, 'fraction_numerical': 0.010828380052795086, 'mean_word_length': 4.653030149223429, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Gigahertz Sub-Landauer Momentum Computing

We introduce a fast and highly-efficient physically-realizable bit swap. Employing readily available and scalable Josephson junction microtechnology, the design implements the recently introduced paradigm of momentum computing. Its nanosecond speeds and sub-Landauer thermodynamic efficiency arise from dynamically storing memory in momentum degrees of freedom. As such, during the swap, the microstate distribution is never near equilibrium and the memory-state dynamics fall far outside of stochastic thermodynamics that assumes detailed-balanced Markovian dynamics. The device implements a bit-swap operation -- a fundamental operation necessary to build reversible universal computing. Extensive, physically-calibrated simulations demonstrate that device performance is robust and that momentum computing can support thermodynamically-efficient, high-speed, large-scale general-purpose computing that circumvents Landauer's bound.


I. INTRODUCTION
Ever since the first exorcism of Maxwell's demon [1], determining how much energetic input a particular computation requires has been a broadly-appreciated theoretical question. In the current century, however, the question has taken on a markedly practical bent; a familiar example is the evolution of Moore's Law from initially provocative speculations decades ago to now addressing material, thermodynamic, and fabrication restrictions [2][3][4][5][6]. Transistor-based microprocessing presents fundamental scaling challenges that strictly limit potential directions for future optimization, and these challenges are no longer speculative. Clock speed, to take one example, has been essentially capped for two decades due to energy dissipation at high rates [7,8]. By some measures, Moore's law is already dead-as integrated circuit manufacturers go vertical, rather than face the expense of creating smaller transistors for 2D circuits that yield only marginal gains [9][10][11]. Given predicted explosive growth in societal demands for information processing and that digital microelectronics is now approaching the physical limits of available architectures [12], exploring alternative computing paradigms is not only prudent but necessary. One alluring vision for the future involves hybrid devices, composed of a suite of computing modules-classical/quantum, digital/analog, deterministic/thermal-each with its own architecture and function that operate in concert. A hybrid architecture allows dynamically harnessing the processing node best suited for the task at hand. The underlying insight is that a computing device's physical substrate should match its desired processing function [13]. In keeping with this, momentum computing demonstrated that low dissipation operations do not require quasi-static operation [14]. That is, energy-efficient computation can be fast in a low-dissipation device. Reference [14] introduced a design framework and theory for an arbitrarily low-cost, high-speed bit swap, a logicallyreversible gate (the only known logical framework with no nontrivial lower bound on its dissipation [12,15,16].) It demonstrated that a universal reversible gate-a Fredkin gate [17,18]-can be built by coupling three such devices together. However, any particular physically-instantiated implementation will come with its own restrictions and considerations that are likely to disallow performing the swap exactly as theorized. And so, an implementation linked to a particular substrate must be built and analyzed in its own right. We present a physically-realizable device and control protocols that implement a bit swap gate that operates in the sub-k B T energy regime using superconducting Josephson junctions (JJs)-a well-known and scalable microtechnology. We recently used this device to measure the thermodynamic performance of bit erasure [19,20]. That extensive experimental effort demonstrated in practical terms that the device proposed here is realizable with today's microfabrication technologies and allows for detailed studies of thermodynamic costs. And so, the device's design and control protocol open up exploring the energy scales of highly energy-efficient, high-speed, general-purpose computing. The Landauer While there are many different quantities one might wish to optimize, the perspective here sets the goal as minimizing the net work invested W when performing logical operations. It is well known that the most pressing physical limits on modern computation are power constraints [21], thus the measure is well suited to diagnose the problems with current devices as well as potential strengths of new ones. For over half a century now Landauer's Principle has arXiv:2202.07122v2 [cond-mat.stat-mech] 18 Nov 2022 exerted a major impact on the contemporary approach to thermodynamic costs of information processing [22,23]. Its lower bound of k B T ln 2 energy dissipated per bit erased has served as standard candle for energy use in physical information processing. To aid comparing other computing paradigms and protocols, we refer to this temperature-dependent information-processing energy scale as a Landauer: approximately a few zeptojoules at room temperature, and a few hundredths of a zeptojoule at liquid He temperatures. See Appendix A for further comparisons. To appreciate the potential benefits of momentum computing operating at sub-Landauer energies we ask where contemporary computing is on the energy scale. Consider recent stochastic thermodynamic analyses of single-electron transistor logic gates [24,25]-analogs to conventional CMOS technology. The upshot is that these technologies currently operate between 10 3 and 10 4 Landauers. More to the point, devices using CMOS-based technology will only ever be able to operate accurately above ≈ 10 2 Landauers [12,15]. In short, momentum computing promises substantial improvements in efficiency with no compromise in speed. Outline Here we provide a brief overview of each section and appendix in the text. Section II explains the importance of bit-swap operations and summarizes the protocol presented in Ref. [14]. Section III introduces the physical substrate, highlights why it is a good candidate, and addresses design restrictions. Section IV reports quantitative results on device performance as measured through detailed simulations of the microscopic degrees of freedom. Section V compares them to related results, both contemporary and foundational. Section VI concludes, summarizing the results and briefly outlining future directions and challenges for scaling up to general-purpose computing. Appendices include details necessary to understand the process by which the parameter space of control protocols was restricted and local work minima were found in simulation. Additionally, they also provide expository information that the interested reader might find relevant. In particular, Appendix A discusses the temperature dependent energy scale, the "Landauer". Appendix B outlines key physical differences between continuous-time Markov chains and hidden Markov chains. Appendix C presents the equations of motion of the bit-swap Josephson junction circuit in their dimensional form and their transformation to simulation-appropriate dimensionless equations. Appendix D details the process of algebraically eliminating large swaths of protocol parameter space. And, finally, Appendix E discusses the algorithmic details of the simulations.

II. BIT SWAP
The Landauer cost stood as a reference for so long since bit erasure is the dominant source of unavoidable dissipation when implementing universal computing with transistor logic gates. It is the elementary binary computation that most changes the Shannon entropy of the distribution over memory states. In this way, one sees k B T ln 2 not just as the cost of erasure, but as the cost of the maximally dissipative elementary operation on which conventional computing relies. And so, the Landauer naturally sets the energy scale for conventional computing.
Taking inspiration from Landauer's pioneering work, we investigate the cost of the most expensive operation necessary to physically implement universal momentum computing: a bit swap. The ideal bit swap has no error, but in the thermodynamic setting one is also interested in an implementation's fidelity. And so, we write a swap with error rates 0 and 1 as a stochastic mapping between memory states m ∈ {0, 1} from time 0 to time τ : The bit swap's dominance in the cost of universal momentum computing can be appreciated by considering the input-output mapping of the Fredkin gate-a 3-bit universal gate with memory states m x m y m z , m i ∈ {0, 1}. All inputs are preserved except for the exchange 101 ↔ 110. We can decompose the informational state space into two regions. If m x = 0, the operation is simply an identity, which trivially is costless. If m x = 1 and m y = m z , we once again have an identity. Thus, it is only the subspace of m x = 1, where m y = m z that a swap must take place. Reference [14] provides explicit potentials that impose effectively 1D swap potentials on a full 3-bit state space in order to implement the Fredkin gate, demonstrating that only 1D swap operations need contribute to the operation's thermodynamic cost.

A. Momentum Computing Realization
Storing information in a one-dimensional state space, it is not clear how to operate a thermodynamically-efficient bit swap with high accuracy. (In this, we recall the conventional interpretation of efficient to mean quasistatic or constantly-thermalizing Markovian dynamics [26,27].) At time t in the operation, the distribution of initial conditions corresponding to m(t = 0) = 0 must overlap with that corresponding to m(t = 0) = 1. And, from that point forward it is impossible to selectively separate them based on their initial positions. Information, and so reversibility, is lost.
Consider, instead, a computation that happens faster than the equilibration timescale of the physical substrate and its thermal environment. In this regime, a particle's instantaneous momentum can be commandeered to carry useful information about its future behavior. Our protocol operates on this timescale, using the full phase space of the underlying system's degrees of freedom to transiently store information in their momenta. Due to this, the instantaneous microstate distribution is necessarily far from equilibrium during the computation. Moreover, the coarse-grained memory-state dynamics during the swap are not Markovian; despite both the net transformation over the memory states and the microscopic phase space dynamics being Markovian. Nonetheless, the system operates orders of magnitude more efficiently than current CMOS but, competing with CMOS, the dynamics evolve nonadiabatically in finite time-on nanosecond timescales for our physical implementation below.
In this way, momentum computing offers up device designs and protocols that accomplish information processing that is at once fast, efficient, and low error. There is a trade-off-a loss of Markovianity in the memory-state dynamics. That noted, the dynamics of the memory states are faithfully described by continuous-time hidden Markov chains (CTHMCs) [28][29][30], rather than the continuous-time Markov chains (CTMCs) that are common in stochastic thermodynamics [26,27]. See Appendix. B for a brief review.

B. Idealized Protocol
Reference [14] describes a perfectly-efficient protocol for implementing a swap in finite time. The operation is straightforward. We begin with an ensemble of particles subject to a storage potential. The potential energy landscape V store (x) must contain at least two potential minima-positioned, say, at x = ±x 0 -with an associated energy barrier equal to max{V store (x), x ∈ (−x 0 , x 0 )} − V store (x 0 ). During storage, a particle's environment is a thermal bath at temperature T . As the height of the potential energy barrier rises relative to the bath energy scale k B T , the probability that the particle transitions between left (x < 0) and right (x ≥ 0) decreases exponentially. In this way, if we assign the left half of the position space to memory state 0 and the right half to memory state 1, the energy landscape is capable of metastably storing a bit m ∈ {0, 1}.
At the protocol's beginning, we instantaneously apply a new potential energy landscape V comp ≡ kx 2 /2. The system is then temporarily isolated from its thermal environment, resulting in the particles undergoing a simple harmonic oscillation. Waiting a time τ until the oscillation is only half completed, the potential is returned to V store . The initial conditions-for which x 0 < 0 (x 0 > 0)-have then been mapped to x τ > 0 (x τ < 0), achieving the desired swap computation. If V store is an even function of x, the computation requires zero invested work as well. This follows since the harmonic motion created a mirror image to the original distribution and the energy imparted to the system at t = 0 is completely offset by the energy extracted from the system turning off V comp at t = τ .  [32]. This is the origin of the variable subscripts to distinguish ϕ from ϕ dc and ϕx from ϕ xdc .

III. PHYSICAL INSTANTIATION
Due to its conceptual simplicity the protocol does not require any particular physical substrate. That said, the practical feasibility of performing such a computation must be addressed. One obvious point of practical concern is assuming the system can be isolated from its thermal environment during the computation. However, total isolation is not necessary. If τ τ R -the relaxation timescale associated with the energy flux rate between the system and its thermal bath-then the device performs close to the ideal case of zero coupling. As proof of concept, Ref. [14]'s simulations showed that this class of protocol is robust: thermodynamic performance persists in the presence of imperfect isolation from the thermal environment, albeit at an energetic cost. Thus, a system that obeys significantly-underdamped Langevin dynamics is an ideal candidate as the physical substrate for bit swap. We analyze in detail one physical instantiation-a gradiometric flux logic cell (Fig. 1), a mature technology for information processing. With suitable scale definitions, the effective degrees of freedom-Josephson phase sum ϕ and difference ϕ dc -follow a dimensionless Langevin equation [19,20,[31][32][33][34]: where x ≡ (ϕ, ϕ dc ) and v ≡ (φ,φ dc ) are vector representations of the dynamical coordinates. Enacting a control protocol on this system involves changing the parameters of the potential over time: = (ϕ − ϕ x (t )) 2 /2 + γ(ϕ dc − ϕ xdc (t )) 2 /2 + β cos ϕ cos(ϕ dc /2) − δβ sin ϕ sin(ϕ dc /2) .
The relationships between the circuit parameters and the parameters in the effective potential U are as follows. ϕ = (ϕ 1 + ϕ 2 )/2 − π and ϕ dc = (ϕ 2 − ϕ 1 ), where ϕ 1 and ϕ 2 are the phases across the two Josephson elements; ϕ x = 2πφ x /Φ 0 − π and ϕ xdc = 2πφ xdc /Φ 0 , where Φ 0 is the magnetic flux quantum and (φ x , φ xdc ) are external magnetic fluxes applied to the circuit; U 0 = (Φ 0 /2π) 2 /L, γ = L/2 , β = I + 2πL/Φ 0 , and δβ = I − 2πL/Φ 0 , where L and 2 are geometric inductances; and I ± ≡ I c1 ± I c2 are the sum and difference of the critical currents of the two Josephson junctions. All parameters are real and it is assumed that γ > β > 1 δβ. Some particularly important parameters of U are ϕ x and ϕ xdc , which control the potential's shape by where the the dynamical variables ϕ and ϕ dc localize in equilibrium, and γ, which controls how quickly ϕ dc localizes to the bottom of the quadratic well centered near ϕ dc = ϕ xdc . At certain control parameters (ϕ x , ϕ xdc ), the effective potential contains only two minima: one located at ϕ < 0 and one at ϕ > 0. So, the device is capable of metastably storing a bit, as described above. In point of fact, the logic cell has been often used as a double well in ϕ with a controllable tilt and barrier height [19,32,34].
The Langevin equation's coupling constants, λ and η, determine the rate of energy flow between the system and its thermal environment and the. They depend on the parameters L, R, and C. In the regimes at which one typically finds L, C, and R and with temperatures around 1 K, the system is very underdamped; ring-down times are O(10 3 ) oscillations about the local minima. (Notably, the device thermalizes at a rate proportional to R −1 . A tunable R allows the device to transition from the underdamped to overdamped regime, allowing for rapid thermalization, if desired.) Finally, θ is a dimensionless factor that depends on the relative inertia of the two degrees of freedom, it depends on the circuit architecture. Appendix C gives the equations of motion and thorough definitions of all parameters and variables in terms of dimensional quantities.

A. Realistic Protocol
With the device's physical substrate set, we now show how to design energy-efficient bit-swap control protocols. There are four parameters that depend primarily on device fabrication: I c1 , I c2 , R, and C. Two that depend on the circuit design: L and . And, four that allow external control: ϕ x , ϕ xdc , T (the environmental temperature), and τ (the computation time). Without additional circuit complexities to allow tunable L, R, and C, we assume that once a device is made, any given protocol can only manipulate ϕ x , ϕ xdc , T , and τ . A central assumption is that computation happens on a timescale over which the thermal environment has minimal effect on the dynamics, so the primary controls are ϕ x , ϕ xdc , and τ . ϕ x is associated with asymmetry in the informational subspace, and will only take a nonzero value to help offset asymmetry from the δβ term in U . Thus, ϕ xdc primarily controls the difference between V comp and V store , while τ governs how long we subject the system to V comp .
V store must be chosen to operate the device in a parameter regime admitting two minima on either side of ϕ = 0 as in Fig. 2. They must also be sufficiently separated so that they are distinct memory states when immersed in an environment of temperature T . In the ideal case, V comp is a quadratic well with an oscillation period τ = π m/k. However, U will never give an exact quadratic well unless β = δβ = 0. So, a suitable replacement is necessary. The closest approximate is at the relatively obvious choice ϕ xdc = −2π. In this case, the minima of both the quadratic and the periodic part of the potential lie on top of each other and the potential is well approximated by a quadratic function over most of the relevant position-domain. However, due to restrictions on V store , transitioning between V store and V comp may induce unnecessarily large dissipation since the oscillations in the ϕ dc dimension have a large amplitude. (See Appendix D for details.) Instead, to dissipate the minimum energy, the control parameters must balance placing the system as close as possible to the pitchfork bifurcation where the two wells merge, while still maintaining dynamics that induce the ϕ < 0 and ϕ > 0 informational states to swap places due to an approximately harmonic oscillation. Near this parameter value, one typically finds a "banana-harmonic" potential energy landscape. (See Fig. 2 for a comparison of the distinct potential profiles for storage and computation.)

B. Computation Time
The final design task determines the computation timescale τ . Under a perfect harmonic potential, the FIG. 3. A dynamic computation: 1, 500 trajectories from V store 's equilibrium distribution in the ϕ (top) and ϕ dc (bottom) dimensions. V comp is applied at t ∈ (1, 1 + τ ), denoted by heavy black lines. ϕ dc oscillations are several times faster than the others, as expected when γ 1. The work done on the system by the control apparatus, W0 = V comp (t = 1) − V store (t = 1), by its intervention at t = 1 is largely offset by the work absorbed into the apparatus by its intervention when V comp re-engages. Visually, we can track this energy flux by the nonequilibrium oscillations induced at t = 1 and the return to a near-equilibrium distribution at t = 1 + τ . Time is measured in units of √ LC, which is ≈ 2ns for the JJ device.
most energetically efficient τ is simply π m/k. This ensures that x(t = 0) = −x(t = τ ). Since the design has an additional degree of freedom beyond that necessary-the ϕ dc dimension-however, we must not only ensure our information-bearing degree of freedom switches signs, but also ensure that ϕ dc (t = 0) ≈ ϕ dc (t = τ ). This means that during time τ , the ϕ variables must undergo n + 1/2 oscillations and the ϕ dc variables must undergo an integer number of complete oscillations. (See Fig. 3.) Hence, τ must satisfy matching conditions for the periods of the oscillations in both ϕ and ϕ dc during the computation: showcases this by displaying the behavior observed during simulations near the ideal timescale. The local work minima coincide with local minima in the average kinetic energy, but not every kinetic energy minimum coincides with a work minimum. While there are kinetic energy minima every half-integer oscillation in ϕ dc , only integer multiples of ϕ dc oscillations yield minimum work.
The equations of motion governing the system are stochastic, dissipative, and nonlinear, so the frequencies of the different oscillations ω, ω dc are nontrivial nonlinear stochastic mappings of device parameters, initial positions, and protocol parameters. They are not easily determined analytically. However, they change smoothly with small changes in the parameters they depends on. Thus, we were able to use an algorithmic approach to find the timescales that yield local minima and explore the regions surrounding them.

C. Physically-Calibrated Bit Swap
We are most interested in the effect of parameters that are least constrained by fabrication. And so, all simulations assume constant fabrication parameters with I + , R, and C set to 2.0 µA, 371 Ω, and 4.0 nF, respectively. To explore how the I − asymmetry affects work cost, we simulated protocols with both a nearly-symmetric device (I − = 7 nA) and a moderately-asymmetric device (I − = 35 nA). Given devices with the parameters above, what values of the other parameters yield protocols with minimum work cost? This involves a twofold procedure. First, create a circuit architecture by setting L and γ, thus fully specifying the device; details in Appendix E. Second, determine the ideal protocols for that combination of device parameters.

D. Computational Fidelity
To determine the best successful protocol, we must define what a successful bit swap is. First, we set a lower bound for the fidelity f : f ≥ 0.99. We define f over an ensemble of N independent trials as: f = 1 − N e /N , with N e counting the number of failed trials, trials for which sign[ϕ(t = 0)] = sign[ϕ(t = τ )]. Second, the distribution over both ϕ(t = τ ) and ϕ(t = 0) must be bimodal with clear and separate informational states. The criteria used for this second condition is: were σ s and s are standard deviations and means of ϕ conditioned on statement s being true.
The final choice concerns the initial distribution from which to sample trial runs. For this, we used the equilibrium distribution associated with V store with the environmental temperature set to satisfy k B T = 0.05U 0 . Here, we ensure fair comparisons between different parameter settings by fixing a relationship between the potential's energy scale and that of thermal fluctuations. This resulted in temperatures from 400 − 1400 mK, though it is possible to create superconducting circuits at much higher temperatures [35][36][37][38] using alternative materials.
Sampling initial conditions from a thermal state assumes no special intervention created the system's initial distribution. We only need wait a suitably long time to reach it. Moreover, this choice is no more than an algorithmic way to select a starting distribution. It is not a limitation or restriction of the protocol. Indeed, if some intervention allowed sampling initial conditions from a lower-variance distribution, it could be leveraged into even higher performance.

IV. PERFORMANCE
Appendix E lays out the computational strategy used to find minimal W implementations among the protocols that satisfy the conditions above. Since the potential is held constant between t = 0 and t = τ , work is only done when turning V comp on at t = 0 and turning it off at ) . Thus, the mean net work cost is the sum W = W 0 + W τ . As we detail shortly, this yielded large regions of parameter space that implement bit swaps at sub-Landauer work cost. This result and others demonstrate the notable and desirable aspects of momentum computing: accuracy, low thermodynamic cost, and high speed. Let's recount these one by one.

A. Accuracy
Tradeoffs between a computation's fidelity and its thermodynamic cost are now familiar-an increase in accuracy comes at the cost of increased W or computation time [39][40][41][42][43][44]. These analyses conclude that accuracy generally raises computation costs. Momentum computing does not work this way. In fact, it works in the opposite way. The low cost of a momentum computing protocol comes from controlling the distribution over the computing system's final state. Due to this, fidelity and low operation cost are not in opposition, but go hand in hand, as Figs. 4 and 5 demonstrate.

B. Low Thermodynamic Cost
Conventional computing, based on transistor-network steady-state currents, operates nowhere near the theoretical limit of efficiency for logical gates. Even gates in Application Specific Integrated Circuits (ASICs) designed for maximal efficiency operate on the scale of 10 4 − 10 6 Landauers [45,46]. The physically-calibrated simulations described above achieved average costs well below a Landauer for a wide range of parameter values with an absolute minimum of W min = 0.43 Landauers, as shown in Figure 6 (left). For the less-ideal asymmetric critical-current device (right panel), the cost increases to only W min = 0.60 Landauers. And, the bulk of the protocols we explored operated at < 10 Landauers. Altogether, the momentum computing devices operated many orders of magnitude lower than the status quo. Moreover, the wide basins reveal robustness in the device's performance: an important feature for practical optimization and implementation.

C. High Speed
Paralleling accuracy, the now-conventional belief is that computational work generally scales inversely with the computation time: W ∼ 1/τ [39,[47][48][49]. Again, this is not the case for momentum computing, as Figs. 4 and 5 demonstrate. Instead, there are optimal times τ * that give local work minima and around which the work cost increases.
Optimal τ * s are upper bounded: the devices must operate faster than particular timescales-timescales determined by the substrate physics. The bit swap's low work cost requires operating on a timescale faster than the rates at which the system exchanges energy and information with the environment. Thus, momentum computing protocols have a speed floor rather than a speed limit. However, even assuming perfect thermal isolation there is a second bound on τ * . The computation must terminate before the initially localized ensemble-storing the memory-decoheres in position space due to dispersion. For our JJ device this is the more restrictive timescale. Due to local curvature differences in the potential, the initially compact state-space regions corresponding to peaks of the storage potential's equilibrium distribution begin to decohere after only one or two oscillations. Once they have spread to cover both memory states, the stored information is lost. This means it is most effective to limit the duration of the swap to just a half-oscillation of the ϕ coordinate. For our devices, this typically corresponds to operating on timescales < 15 ns.

V. RELATED WORK
Reversible computing implementations of various operations have been proposed many times over many decades. Perhaps the most famous is the Fredkin billiards implementation [18]. While ingenious, it suffers from inherent dynamical instability (deterministic chaos) and cannot abide any interactions with the environment. At the other end of the spectrum is a family of superconducting adiabatic implementations [50][51][52][53][54][55][56][57]. These are low cost Note that when the asymmetry is low, it can effectively be offset by the parameter ϕx, but for higher asymmetry, protocols that cost less than 1 Landauer are less common.
in terms of dissipation and are stable, but they suffer from fundamental speed limits due to the adiabaticity requirement: W ∝ 1/τ .
Other recent implementations [58][59][60] of reversible logic using JJs are more akin to the proposal at hand, in that they require nearly-ballistic dynamics and attempt to recapture the energy used in a swap at the final step. While these implementations are markedly different, their motivation follows similar principles. Particularly, the framework for asynchronous reversible computing proposed in [59,60] might serve as a testbed for momentum computing elements.
Another distinguishing feature of the present design is that the phenomenon supporting the computing is inherently linked to microscopic degrees of freedom evolving in the device's phase space. This moves one closer to the ultimate goal of using reversible nanoscale phenomena as the primitives for reversible computing-a goal whose importance and difficulty were recognized by Ref. [61].
Working directly with the underlying phase space also allows incorporating the thermal environment. And, this facilitates characterizing the effect of (inevitable) imperfect isolation from the environment.
It is worth noting the similarity between the optimal timescales τ * and the principal result in Ref. [62] in which a similar local minima emerges when comparing thermodynamic dissipation to computation time. These minima also come from certain matching conditions between the rate of thermalization and the system's response time to its control device. Another qualitatively similar result [63] found faster operation could lead to reduced errors in overdamped JJs under periodic driving. These similarities could point to a more general principle at play.

VI. CONCLUSION
Our detailed, thermodynamically-calibrated simulation of microscopic trajectories demonstrated that momentum computing can reliably (i) implement a bit swap at sub-Landauer work costs at (ii) nanosecond timescales in (iii) a well-characterized superconducting circuit. These simulations served two main purposes. The first highlights momentum computing's advantages. The proposed framework uses the continuum of momentum states to serve as the auxiliary system that allows a swap. In doing so, it eliminates the associated tradeoffs between energetic, temporal, or accuracy costs that are commonly emphasized in thermodynamic control analyses [39][40][41][42].
Momentum computing protocols are holistic in that low energy cost, high fidelity, and fast operation times all come from matching parallel constraints rather than competing ones.
The second purpose points out key aspects of the proposed JJ circuit's physics. The simulations reveal several guiding principles-those that contribute most to decreasing work costs for the proposed protocols. The system is so underdamped that thermal agitation is not the primary cause of inefficiency. The two main contributors are (i) the appearance of dispersive behavior in the dynamics of an initially-coherent region of state space and (ii) asymmetries inherent to the device that arise from differing critical currents in the component superconducting JJ elements. Notably, if the elements are very close to each other in I c , then symmetry can be effectively restored by setting the control parameter ϕ x to counteract the difference. However, the more asymmetry, the harder it is to find ultra low-cost protocols; cf. Fig. 6 left and right panels. Note, too, that initial-state dispersion can be ameliorated by using a V comp that is as harmonic (quadratic) as possible. However, this typically requires lower inductance L, possibly complicating circuit fabrication. Additionally, the potential-well separation parameter β's linear dependence on L hinders the system's ability to create two distinct states during information storage. Though these tradeoffs are complicated, our simulations suggest that dispersion can be controlled, yielding swap protocols with even lower work costs. Since the protocol search space is quite high-dimensional and contains many local-minima, we offer no proof that the protocols found give the global work minimum. Very likely, the thermodynamic costs and operation speed of our proposed JJ momentum computing device can be substantially improved using more sophisticated parameter optimization and alternative materials. Even with the work cost as it stands, though, sub-Landauer operation represents a radical change from transistor-based architectures. One calibration for this is given in the recent stochastic thermodynamic analysis of a NOT gate composed of single-electron-state transistors [24] that found work costs 10 4 times larger. Note, too, that running at low temperatures requires significant off-board cooling costs, as required in superconducting quantum computing. Our current flux qubit implementation requires operating at liquid He temperatures [19,20]. However, there are also JJs that operate at N 2 temperatures, promising system cooling costs that would be 2 to 3 orders of magnitude lower [35][36][37][38]. Additionally, the physics necessary to build a momentum computing swap-underdamped behavior and controllable multiwell dynamics-is far from unique to superconducting circuits. As an example, nanoelectromechanical systems (NEMS) are another well-known technology that is scalable with modern microfabrication techniques. NEMS provide the needed nonlinearity for multiple-well potentials, are extremely energy efficient, and have high Q factors even while operating at room temperature [64][65][66]. Momentum computing implemented with NEMS rather than superconductors completely obviates the cooling infrastructure and so may be better suited for large-scale implementations.
That said, the JJ implementation at low temperatures augmented with appropriate calorimetry will provide a key experimental platform for careful, controlled, and detailed study of the physical limits of the thermodynamic costs of information processing. Thus, these devices are necessary to fully understand the physics of thermodynamic efficiency. And so, beyond technology impacts, the proposed device and protocols provide a fascinating experimental opportunity to measure energy flows that fluctuate at GHz timescales and at energy scales below thermal fluctuations. Success in these will open the way to theoretical investigations of the fundamental physics of information storage and manipulation, time symmetries, and fluctuation theorems [43,67].

Appendix A: The Landauer: A Standard Candle for Thermodynamic Computation
A long and checkered history underlies the physics of information and energy, arguably originating in the paradox of Maxwell's Demon [68]. Most recently, though, the paradigm of thermodynamic computing emerged to frame probing their limits [69]. In this setting, Landauer's Principle says that k B T ln 2 energy units must be expended to erase a single bit of information. Beyond erasure, though, his Principle also stands as a challenge-Can conventional computing paradigms operate at sub-Landauer scales? It seems not. Landauer's theory and follow-on results [70][71][72] and recent experiments [42,73] verified the lower bound.
To apply more broadly, Landauer's Principle generalizes to W ≥ k B T ∆H, where ∆H is the change in Shannon entropy between a computational system's initial and final information-bearing states [71,74,75]. Despite the Principle's generalization beyond bit erasure, the Landauer scale remains a familiar reference point for the energy costs of binary operations; its familiar use coming at the expense of ignoring specifics of any given logical operation [44].
An efficient bit-swap operation, for example, has zero generalized Landauer cost, as it is logically reversible. However, since many thermodynamic computing architectures do not have access to dynamics that can accomplish reversible computing efficiently, the Landauer scale provides a common reference to compare gate performance across physical substrates and design paradigms. It also facilitates comparing across substrates that operate at different temperatures. Table I lists thermodynamic energies for a range of physical environments. Table II gives Landauer work energies for various information processing operations in environments and at temperatures where thermodynamic computers operate.

Appendix B: Limits of Stochastic Thermodynamics for Information Processing
Stochastic thermodynamics [26,27] f ( p, t).
The resulting dynamics are necessarily represented by a master equation over the memory-state distributioṅ p(t) = A(t) p(t) [76]. This framework is powerful, yielding great insight into physical processes when its assumptions are met. The framework, however, does not apply to momentum computing. To appreciate why, consider justifying Markovian dynamics over memory states. Assume a microscopic physical system S that serves as a computational substrate. While allowing the universe to be deterministic, S can exhibit stochastic dynamics since it represents only a portion of the partially-observed universe. The very typical assumption that S's local environment acts as a large weakly-coupled heat bath with quickly relaxing degrees of freedom yields dynamics on S that are also Markovian and, therefore, can be represented by CTMCs. However, computationally-useful memory states are not the CTMC-obeying microstates of S, but a set M of mesostates that represent coarse-graining over S. It is possible, depending on the variables or timescales of interest, that this coarse-graining ignores only rapidly-relaxing subsystems of S. Then M inherits the Markov property that governs the microstates [26]. This strategy-coarse graining over physical degrees of freedom irrelevant to the dynamics-is analogous to establishing S as a stochastic, Markovian subsystem of the universe. A straightforward example of this case is when M consists of positional degrees of freedom and S evolves by overdamped Langevin dynamics. When implementing momentum computing, however, the coarse-graining yielding M is applied over hidden microstates that contain dynamically relevant information not determined from M's instantaneous realizations. As a consequence, the dynamic over the coarse-grained states is not Markovian. CTMCs cannot be used. A straightforward example of this arises when M consists of positional degrees of freedom and S evolves by underdamped Langevin dynamics. On the downside, a general analytical treatment of such partially-observed systems (continuoustime hidden Markov chains) is highly nontrivial [27,[77][78][79]. On the upside, the possibility of hidden states allows for substantially more general forms of computation. As the results here showed, the benefits of this expanded space are quite substantial. Note that the bit swap computation is, in general, problematic to implement using CTMCs since input-output mappings whose determinants are negative are disallowed when memory-state dynamics are restricted to obey CTMCs. Formally, auxiliary systems can be added to the set of memory states. Done correctly this again permits using CTMCs in the augmented state space to accomplish the computation [76]. However, physically-embedded computations do not generally allow the required perfect control over the system Hamiltonian. Indeed, one need look no further than the present work to see how nontrivial it is to implement an operation as simple as a harmonic oscillation in a physically-realistic device. Moreover, adding auxiliary subsystems increases statespace dimension and complicates control apparatus and   control protocols. Due to the increased complication, in many settings, adding auxiliary dimensions is simply not physically possible. On top of this, the timescale of these augmented computations must be longer than the equilibration time of the auxiliary systems and thermal environment. In this way, adding auxiliary systems imposes additional speed limits to computations. In short, adding auxiliary subsystems addresses the shortfalls of CTMCs, but does not sidestep their fundamental limitations.
We illustrate this by considering an efficient bit swap implemented via a Markovian embedding. First, it augments the system with an unoccupied auxiliary state A to serve as a transient memory. It then quasistatically translates memory state 0 to A, while memory state 1 is translated to 0. Finally, it quasistatically translates A to 1.
Quasistatic processes cost arbitrarily little work, but they take arbitrarily-long times. To compute faster (τ → 0), the work cost will diverge as 1/τ [39,[47][48][49]. Increasing fidelity requires raising the scale of the barrier separating the states. Doing so, though, increases the energetic cost at a given computational speed; maintaining the same work cost, then, requires slowing the operation. In short, the trade-offs in Markovian embedding complicate design and, more to the point, reduce performance.

Appendix C: Flux Qubit Dimensionless Equations of Motion
In terms of the dimensional degrees of freedom, the flux qubit equations of motion are: where the dimensional ϕs are related to the main text's dimensionless fluxes and phases by the magnetic flux quantum 2π/Φ 0 . With the addition of thermal noise, the Langevin equation is: where κ ≡ k B T . Matching these variables to the equations of motion yields: where subscript i has been dropped in favor of a vector representation.
The task is to write each physical quantity z in terms of a dimensional constant and dimensionless variable by defining scaling factors according to the following prescription: z ≡ z z c , where z c is a dimensionful constant.
Setting m c = C and ν c = 1/R are obvious choices. Additionally, since the potential factors into U = U 0 ×U ( 2π Φ0 ·x), a good choice for positional scaling is It is advantageous to write nondimensional kinetic energies as 1 2 m v 2 without additional scaling factors. This means setting the energy scaling as: This does not uniquely determine the energetic scale, since t c is still free. The two obvious choices are to scale to the temperature-KE = 1 corresponds to k B T units of dimensional energy-or to the potential energy scale-KE = 1 corresponds to U 0 units of dimensional energy. Choosing the latter yields: Evidently, the timescale is t c = √ LC, which is a workable timescale for our purposes given that the dynamics of interest happen on the scale of τ ≈ ω LC . Setting the timescale to the potential energy rather than the thermal energy may well become common practice in simulating momentum computation, since protocols must be timed precisely with respect to the dynamics of the potential energy surface.
The Langevin equation, in terms of the nondimensional quantities defined above, becomes: Simplifying algebra then yields: 2dt .
Finally, we define λ, θ, and η as nondimensional parameters that serve as our dimensionless Langevin coefficients. This yields the Langevin equation for the simulations detailed in Appendix E: with: where:

Symmetric Approximation
We can obtain reasonable estimates for good ϕ xdc values by assuming a perfectly symmetric device δβ = 0. Furthermore, we also set ϕ x = 0 for all cases. This allows two symmetric wells on either side of ϕ = 0. In practice, since δβ = 0 in a real device, ϕ x would be calibrated to compensate for the asymmetry; see Sec. D 2.
In the symmetric case, the potential splits into two components-periodic and quadratic: The periodic term allows for multiple minima, while the quadratic terms force the dynamical variables to stay close to their respective parameters. This localization means we focus only on the the area near ϕ = ϕ x and ϕ dc = ϕ xdc .
To employ the potential most flexibly, we must characterize the relevant fixed points that occur in this region. Following Refs. [19,33], we choose to search in the domain −π < ϕ < π and −2π < ϕ dc < 0. Fixed points occur when all components of the gradient vanish: The first condition is met whenever ϕ = 0 and, also, when ϕ β sin ϕ = cos 1 2 ϕ dc . Consider the case where ϕ = 0-the "central" fixed point. To find the ϕ dc location of the fixed point ϕ 0 dc , we look to the gradient's second term. This yields the condition: The central fixed point occurs close to the parameter ϕ xdc , but is offset by a value ≤ β/2γ. The equation above can be solved numerically with ease to find the location of the central fixed point. To classify the fixed point, we look at the Hessian. While the general expression for the eigenvalues is rather verbose, the case where ϕ = 0 simplifies to: λ 2 > 0 as long as γ > β/4. And, since we assume γ > β, this condition is always met. Thus, this fixed point is either a saddle point or a minimum based on whether ϕ 0 dc is greater or less than ϕ c dc ≡ −2 cos −1 1 β , respectively. (We only use the negative branch of cos −1 due to the domain of ϕ dc .) See Fig. 7 for an example of the behavior of the central fixed point for typical parameters. We can also find an expression for ϕ c xdc (β, γ) ≡ F 0 (ϕ dc = ϕ c dc ), the critical value of the control parameter at which the central fixed point transitions between a saddle point and a minimum: Naively, the best strategy to form a low cost protocol is to take values of ϕ xdc just above and below ϕ c xdc . However, there are several factors that introduce complications. For one, the energy scale separating the two wells when ϕ xdc ≈ ϕ c xdc is very small and it will typically be overwhelmed by thermal energy at the temperatures of interest (400−1400 mK). A second is that the approximation of δβ = 0 actually has a most pernicious effect near to a subcritical pitchfork bifurcation-yielding a potential with a third (undesirable) minimum rather than a single one.
When ϕ = 0 we can rewrite Eqs. (D2) and (D3): The potential is symmetric, so these fixed points come in pairs ϕ ± . Substituting ϕ dc /2 = − cos −1 (ϕ ± /β sin ϕ ± ) into the second equation yields the following for ϕ ± : Note that the sign changes due to the domain restriction of ϕ dc . Figure 8 shows how these fixed points behave as β, γ, and ϕ xdc change. The value of ϕ xdc tangent to the curve when ϕ = 0 corresponds to the critical control parameter value ϕ c xdc , which can be seen by verifying lim ϕ→0 F ± (ϕ) = ϕ c xdc . As a last note, different values of β and γ have qualita- Once again, the full derivative is quite verbose. However, taking the limit ϕ → 0 gives: Interestingly, when β > β * , there is always a parameter space region with three distinct minima. This might be useful, in fact, for single-bit computations that require more states. For bit swap, though, the goal is for the system to jump between a V store with 2 minima and a V comp with a single minimum (see Figure 10). And so, care must be taken to avoid the three-minima regions when β > β * . It is likely that more evolved solution strategies will improve results.

δβ = 0
The device just considered is ideal. In reality δβ = 0, and exact analytic work is much less fruitful. Introducing the asymmetric terms augments the potential: In short, one must vary ϕ x to offset the effect of δβ, provided a symmetric potential is preferred. There are two obvious strategies to minimize the effects of asymmetry. Either a strategy that minimizes the effect of U asym at the central fixed point-the "min of mid" strategy-or at the fixed points at ϕ ± -the "min of max" strategy. It stands to reason that one uses the former to set ϕ x for V comp and the latter for V store . The "min of mid" strategy is easy to implement. Simply set the derivative of ∂ ϕ U asym | ϕ=0 = 0, with the intent of having the asymmetrical part of the potential be as flat as possible near ϕ = 0. Simple algebra yields: ϕ x = −δβ sin ϕ dc /2. The "min of max" strategy requires numerical solution. First, note that the maximum value of U asym occurs when ϕ = ϕ max = arccos( ϕx δβ sin .5ϕ dc ). Then, use a symbolic solver (e.g., SymPy's nsolve function) to find the value of ϕ x that minimizes U asym (ϕ max , ϕ x , δβ, ϕ dc ). Figure 9 shows that the effect of δβ = 0 is, unsurprisingly, the most noticeable near the bifurcation of the central fixed point. For the bit swap, as described in Sec. II B, we need only two different profiles for the potential: one in which we have two symmetric wells and one in which we have a single well placed midway between them. Thus, we must keep the ϕ xdc parameter sufficiently far away from ϕ c xdc . The strategy employed in the simulations described below always involves setting a minimum distance that ϕ xdc must be from ϕ c xdc , in order to avoid falling into the pitfalls described here.

Appendix E: Searching for Minimal-Work Bit Swaps
The following lays out the computational strategy to find low work-cost implementations.
We are most interested in the effect of parameters that are the most removed from fabrication, so all simulations assume JJ elements with I + , R, and C set to 2.0 µA, 371 Ω, and 4.0 nF, respectively. To explore how asymmetry affects work cost, we simulated protocols with a nearly-symmetric device with I − = 7 nA, a moderatelysymmetric device with I − = 35 nA, and an asymmetric device with I − = 60 nA. Additionally, k B T is always scaled to U 0 , so that κ ≡ k B T /U 0 = 0.05.
Given devices with the parameters above, what values of the remaining parameters yield protocols with minimum work cost? This involves a twofold procedure. First, create the circuit architecture by setting L and γ by hand; thus, fully specifying the device. Second, determine the ideal protocols for that combination of device parameters through simulation.
L's order of magnitude was chosen from previous results [19,20,[31][32][33][34] to be 10 −9 H. Noting that a lower L results in a more harmonic potential during computation, we set a minimum L to be 0.3nH. This is in order to stay within the parameter range for which β > 1 and we can still use the analytic expressions derived above. To assure γ > β, γ values were tested in the range [3.0, 20.0].
After choosing a pair of circuit parameters L and γ, we turn to simulation. First, V store must be chosen by setting ϕ store x and ϕ store xdc . This is done by calculating ϕ store xdc ≡ ϕ c xdc + ∆S, where ϕ c xdc (γ, β) is from Eq. (D7). The parameter ∆S is initialized manually to a value ∆S * when starting a new round of simulations. (∆S * = 0.16 was used in the heatmaps shown in Fig. 6.) Then, using the "min of max" method (Sec. D 2), we set ϕ store Finally, V store is tested by sampling 50,000 states from V store 's equilibrium distribution using a Monte Carlo algorithm. The resulting ensemble is verified by determining that it contains two well-separated informational states by asserting that: where s and σ s are means and standard deviations of ϕ conditioned on s being true. If the ensemble fails the test, ∆S is incremented and the process is repeated. If the ensemble succeeds, we have found a viable V store . Then, we move on to establish V comp by choosing ϕ comp x and ϕ comp xdc . Similar to ϕ store xdc , ϕ comp xdc ≡ ϕ c xdc −∆C with ∆C manually set. The value of ∆C does effect the eventual work cost, but the work costs vary smoothly, and a single value of ∆C tends to work well over a large parameter range. Manually setting a single value for ∆C, rather than allowing it to adjust itself to fall into a local minimum, substantially reduces simulation run time. However, we expect that given more compute resources a wider range of sub-Landauer protocols will be discovered. Figure 11 shows the effect of changing ∆C for three different devices. Once ∆C is chosen, we use the "min of mid" (Sec. D 2) method to set ϕ comp x and fully determine V comp . Next, a preliminary simulation is run to identify an approximate value of the computation time τ . To make the simulation run quickly, the ensemble above is coarsegrained into two partitions based on whether ϕ > 0 or ϕ < 0. Then, each partition is coarse-grained again into ≈ 250 representative points through histogramming. A Langevin simulation is run over the histogram data, exposing it to V comp for a time O(10) √ LC. This ensures capturing the time with the best bit swap. Next, weighting the simulation results by histogram counts within each partition, we obtain conditional averages for an approximation of the behavior over the entire ensemble. These averages are parsed for a set of times at which there are indications of a successful and low-cost bit swap: ϕ(t = 0) < 0 > 0, ϕ(t = 0) > 0 < 0, and values of φ and φ dc that are close to zero. See, for example, the blue highlighted portion on the top panel of Fig. 4. In this way, a range (τ min , τ max ) is determined for τ . Now, a larger simulation is completed to determine τ that give the lowest work value. Another 40,000 samples are generated from V store 's equilibrium distribution, and a Langevin simulation is run on the full ensemble by exposing it to V comp for τ max time units. Since the potential is held constant between t = 0 and t = τ , work is only done when turning V comp on at t = 0 and turning it off at t = τ . The average work done at t = 0 is W 0 ≡ V comp (ϕ(0), ϕ dc (0))−V store (ϕ(0), ϕ dc (0)) and returning to V comp at time t costs W t ≡ V store (ϕ(t), ϕ dc (t)) − V comp (ϕ(t), ϕ dc (t)) . Thus, the mean net work cost at time t is the sum W (t) = W 0 + W t . Additionally, for each t ∈ (τ min , τ max ) we calculate the fidelity f (t) and whether the final states are well-separated informational states, s(t): Finally, we choose the minimum work protocol via inf (W (t) : f (t) ≥ 0.99, s(t) = True). After this, we move on to the next pair of L and γ. Typically, these are chosen to be individually close to the last pair. And, and instead of re-initializing ∆S to its initial value by hand, we decrement ∆S from its current value by a small amount if ∆S > ∆S * , using this value as the starting point for the next L and γ pair. This allows the value of ∆S to drift from its starting point towards more favorable values as the parameters change, while still preferring to be close to the known well-behaved parameter value ∆S * . Setting a new initial value for ∆S goes full circle, to find the next minimum work protocol by repeating the procedure.
This procedure yielded rather large ranges of parameter space over which we found very low work-cost bit swap protocols. Here, we offer no proof that the protocols found achieve the global minimum work, since the protocol space is high dimensional and contains many local minima. That said, improved algorithms and a larger parameter-range search should result in even lower work costs.
Langevin simulations of the dimensionless equations of motion employed a fourth-order Runge-Kutta method for the deterministic portion and Euler's method for the stochastic portion of the integration with dt set to 0.005 √ LC. (Python NumPy's Gaussian number generator was used to generate the memoryless Gaussian variable r(t).)