Paper deep dive
CN101 - A Digital Thermodynamic Computer for Generative AI
Lars Holdijk, Denis Melanson, Zier Mensch, Brandon Birchall, Vincent Cheung, Nicholas Lehrter, Maxwell Aifer, Samuel Duffield, Jan Ole Ernst, Rajath Salegame, Antonio J. Martinez, Gavin Crooks, Miranda Cheng, Zach Belateche, Marc Bright, Patrick J. Coles, Faris Sbahi
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 90%
Last extracted: 8/4/2026, 5:03:05 AM
Summary
The paper introduces CN101, a prototype digital thermodynamic computing chip that implements a substrate-independent formalization of equilibration-style computation. Unlike prior analogue implementations restricted to Langevin dynamics, CN101 uses discrete accumulator dynamics on standard CMOS hardware. The authors demonstrate that this approach supports 'anytime precision,' 'parallel sample aggregation,' and 'sequential parallelism,' validating the chip's effectiveness on generative AI workloads such as VAEs and flow matching for image generation and scientific problems.
Entities (9)
Relation Signals (7)
CN101 → developedby → Normal Computing Corporation
confidence 95% · 1 Normal Computing Corporation, New York, USA
CN101 → implements → Thermodynamic Computing
confidence 95% · CN101 implements the formulation through discrete accumulator dynamics on standard CMOS using stochastic computing principles.
CN101 → supports → Sequential Parallelism
confidence 90% · dependent stages of a computation operate concurrently rather than serially, a property we call sequential parallelism.
CN101 → uses → CMOS
confidence 90% · implements the formulation through discrete accumulator dynamics on standard CMOS
CN101 → appliedto → VAEs
confidence 88% · characterise CN101's success across conventional generative AI workloads in the form of VAEs
CN101 → appliedto → Flow Matching
confidence 88% · characterise CN101's success across conventional generative AI workloads in the form of VAEs and flow matching
Langevin dynamics → restrictedto → Analogue Substrates
confidence 85% · formulated exclusively through Langevin dynamics, restricting its implementations to analogue substrates
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Thermodynamic computing is an emerging hardware paradigm, in which stochastic physical dynamics serve as the direct computational primitive. The recent explosion of generative AI has only sharpened the search for alternative approaches to compute, and, as we show in this work, thermodynamic computing turns out to be well suited to this space. An important class of methods realises a function as the stationary expectation of an ergodic stochastic process: the answer is encoded in the time-averaged statistics of an equilibrating trajectory. To date, this equilibration-style class has been formulated exclusively through Langevin dynamics, restricting its implementations to analogue substrates and the engineering challenges those bring. In this work, we propose a substrate-independent formalisation of the equilibration-style formulation, in which the only object of design is the dynamical generator L* of an arbitrary ergodic process. The formalisation makes three hardware-level properties of the formulation explicit: the precision of a result is a knob set by how long the dynamics are run, sample averages decompose across independent trajectories, and dependent stages of a computation operate concurrently rather than serially, a property we call sequential parallelism. We instantiate the formalisation by fabricating a prototype digital thermodynamic computing chip, named CN101, that implements the formulation through discrete accumulator dynamics on standard CMOS using stochastic computing principles. We characterise CN101's success across conventional generative AI workloads in the form of VAEs and flow matching, applied to both image generation and scientific problems. Together, the formalisation and its digital instantiation show that the equilibration-style formulation is substrate-independent, and that its computational properties can be exploited on standard digital hardware.
Tags
Links
- Source: https://arxiv.org/abs/2608.00754v1
- Canonical: https://arxiv.org/abs/2608.00754v1
Trouble viewing inline? Open PDF directly →
Full Text
100,438 characters extracted from source content.
Expand or collapse full text
CN101 - A Digital Thermodynamic Computer for Generative AI Targeting the Sequential Bottleneck of Modern Generative AI Lars Holdijk1,2, Denis Melanson1, Zier Mensch3,4, Brandon Birchall1, Vincent Cheung1, Nicholas Lehrter1, Maxwell Aifer1, Samuel Duffield1, Jan Ole Ernst1, Rajath Salegame1, Antonio J. Martinez1, Gavin Crooks1, Miranda Cheng3,5, Zach Belateche1, Marc Bright1, Patrick J. Coles1, Faris Sbahi1 1 Normal Computing Corporation, New York, USA 2 University of Oxford, Oxford, United Kingdom 3 University of Amsterdam, Amsterdam, Netherlands 4 National Taiwan University, Taipei, Taiwan 5 Academia Sinica, Taipei, Taiwan Correspondence: research@normalcomputing.com, larsholdijk@gmail.com abstractbox Abstract Thermodynamic computing is an emerging hardware paradigm, in which stochastic physical dynamics serve as the direct computational primitive. The recent explosion of generative AI has only sharpened the search for alternative approaches to compute, and, as we show in this work, thermodynamic computing turns out to be well suited to this space. An important class of methods realises a function f(x)f(x) as the stationary expectation of an ergodic stochastic process: the answer is encoded in the time-averaged statistics of an equilibrating trajectory. To date, this equilibration-style class has been formulated exclusively through Langevin dynamics, restricting its implementations to analogue substrates and the engineering challenges those bring. In this work, we propose a substrate-independent formalisation of the equilibration-style formulation, in which the only object of design is the dynamical generator ℒL of an arbitrary ergodic process. The formalisation makes three hardware-level properties of the formulation explicit: the precision of a result is a knob set by how long the dynamics are run, sample averages decompose across independent trajectories, and dependent stages of a computation operate concurrently rather than serially, a property we call sequential parallelism. We instantiate the formalisation by fabricating a prototype digital thermodynamic computing chip, named CN101, that implements the formulation through discrete accumulator dynamics on standard CMOS using stochastic computing principles. We characterise CN101’s success across conventional generative AI workloads in the form of VAEs and flow matching, applied to both image generation and scientific problems. Together, the formalisation and its digital instantiation show that the equilibration-style formulation is substrate-independent, and that its computational properties can be exploited on standard digital hardware. 1. Introduction In recent years, there has been growing interest in thermodynamic computing as an alternative hardware paradigm[15, 14, 29, 2]. The explosion of generative AI has only sharpened this interest: as these workloads strain conventional accelerators, the search for alternative approaches to compute has intensified, and thermodynamic computing turns out to be well suited to them. The unifying premise across the field is that stochastic physical dynamics can serve as a direct computational primitive [15]. Rather than fighting noise to implement deterministic logic, thermodynamic computing uses stochasticity as part of the computational process. Within this broad concept, distinct computational styles have developed, each defining its own relationship between physical dynamics and computational output. Figure 1: The evaluation PCB carries CN101 together with an FPGA, power regulation, and a USB-C host interface. In this work, we focus on the class of methods that realises a function f(x)f(x) as the stationary expectation of an ergodic stochastic process[3, 46], an equilibration-style formulation in which the time-averaged statistics of a stochastic trajectory encode the answer. This class was established by Thermodynamic Linear Algebra (TLA)[3, 46] and validated in analogue hardware by the Stochastic Processing Unit (SPU)[46]; we describe the underlying dynamics in Section˜2.1. The same formulation has also been carried to further problem classes, such as Bayesian inference[4] and natural gradient descent[18]. Subsequent work has extended the formulation to nonlinear and non-equilibrium settings[64, 65], including generative modelling by reverse-time evolution of analogue Langevin systems[65]. Distinct from this line of work, a separate direction of thermodynamic computing has developed around discrete probabilistic bits sampling energy-based models[10, 33, 1, 14, 32]. Methods here are defined by a different computational formulation where the focus is on sampling from a programmable probability distribution rather than evaluating arbitrary functions using a thermodynamic process, as is done in the equilibration-style methods that motivate this work. We do not address the sampling focused direction here. Despite the rapid advances over the last years, a unifying definition for the equilibration-style formulation has not been written down independent of its analogue Langevin instances that have motivated the thermodynamic hardware development so far. Instead, each demonstration of the field carries its own implicit computational contract, with the consequence that each formulation has remained tied to a specific physical implementation. The most immediate practical consequence is that thermodynamic hardware has primarily been confined to analogue substrates, which bring well-known engineering challenges, such as device variability, limited dynamic range, and a manufacturing path that does not benefit from the established CMOS toolchain. To address this, we make three contributions. First, we give a substrate-independent formal definition of equilibration-based computation, in which the dynamical generator of an ergodic process is the only object of design. The analogue Langevin systems of prior work appear as one special case here (Section˜2), but the definition also allows for other formats such as discrete Markov chains. Second, we characterise three scaling properties of thermodynamic computation (anytime precision, parallel sample aggregation, and sequential parallelism) that are exploitable in hardware independently of the substrate. And third, to instantiate this formalisation in silicon, we present a digital thermodynamic computing chip, codenamed CN101, based on stochastic computing principles [22, 5]. CN101 implements the dynamical generator ℒL at the core of the presented thermodynamic computing formulation through accumulator-based finite-state machines updated by pseudo-random bit streams. As the chip is built on standard CMOS, its design directly benefits from established manufacturing and design tooling, in contrast to the analogue implementations of prior work. Notably, CN101 is a prototype built to prove out the substrate-independent formulation in silicon rather than as a finished or commercial design, and many of its blocks are deliberately first-generation, with substantial improvements left to subsequent chips in the Carnot Architecture. Fig.˜1 shows CN101 mounted on its evaluation PCB. To validate CN101’s operation we perform extensive testing across multiple workloads. First, we confirm that CN101 reproduces the formulation’s original problem class, solving linear systems, in fully digital form. We then turn to a conditional variational autoencoder generating MNIST digits, a nonlinear generative model run entirely on-chip and use it to characterise the discussed scaling properties of thermodynamic computing. Using multiple chips, we then shift our focus to modern generative AI methods and show CN101 operating on problems within image generation and the life sciences. Every model we deploy is a standard architecture of the kind trained and run on GPUs (a conditional variational autoencoder and flow-matching models) without specialised modifications. 2. Substrate Independent Thermodynamic Computing In the coming section we will develop the substrate independent formulation of thermodynamic computing, discuss the core challenges in implementing this framework, and outline three key characteristics and the benefits they bring when successful. First, however, we quickly discuss the fundamental direction in this branch of thermodynamic computing built around equilibration. 2.1. Background: Equilibrium Thermodynamic Computing Thermodynamic computing takes its name from a structural analogy with statistical mechanics. A physical system in contact with a heat bath relaxes to thermal equilibrium, where the Boltzmann distribution π(S)∝e−H(S)/kBTπ(S) e^-H(S)/k_BT assigns a probability to each configuration S and macroscopic observables are equilibrium expectations π[φ(S)]E_π[ (S)] of microscopic quantities φ . These averages are never evaluated by an explicit sum over configurations; they emerge, in the sense that once the system has equilibrated, a time average of φ along a single trajectory converges to the ensemble expectation. Thermodynamic computing inverts this relationship: rather than reading out the equilibrium properties of a given system, one engineers a system whose equilibrium expectation is a quantity one wishes to compute. A clear realisation of this idea is Thermodynamic Linear Algebra (TLA) [3, 46]. To solve a linear system As=bAs=b with symmetric positive-definite A, one runs the Ornstein–Uhlenbeck dynamics [62] dSt=−(ASt−b)dt+2dWt,dS_t=-(AS_t-b)\,dt+ 2\,dW_t, (1) whose stationary distribution is the Gaussian (A−1b,A−1)N(A^-1b,\,A^-1). Its two moments give two linear-algebra primitives: the mean A−1bA^-1b is the solution of the system, recovered as the time-averaged state 1T∫0TStdt 1T _0^TS_t\,dt, and the covariance A−1A^-1 is the matrix inverse, recovered from the sample covariance of the state. Neither requires a matrix factorisation. The Stochastic Processing Unit (SPU) [46] realised these dynamics in analogue hardware as a set of capacitively coupled RLC cells, with A encoded in the cells’ tunable capacitances and the cells driven by an injected current-noise source sampling the equilibrium voltages and forming their covariance; it inverted dense matrices in hardware. Two features of this construction are worth noting, as both carry over to the substrate-independent setting developed below. First, the computation is performed by relaxation rather than by an explicit algorithm: the SPU inverts a d×d× d matrix in a time that scales as (d2)O(d^2), against the (d3)O(d^3) of direct factorisation [3]. Second, the readout is available at any time. The running average is already a usable estimate at every T and improves as the process continues. We develop this property into a general scaling axis in Section˜2.3. What has so far been specific to TLA is the choice of Langevin dynamics on a continuous space, and it is precisely this choice that ties these realisations to analogue hardware. The formalism developed removes this restriction. 2.2. The Equilibration Formalism We now state this pattern as a definition that makes no commitment to any particular dynamics or physical substrate: A thermodynamic computation of f:ℝn→ℝmf:R^n ^m is specified by a tuple (S,ℒ,π,φ)(S,L,π, ), where S is a state space, ℒL is a dynamical generator parameterised by f and the input x, π is the unique stationary distribution satisfying ℒπ=0Lπ=0, and φ:S→ℝm :S ^m is a readout observable with π[‖φ(S)‖]<∞E_π[\| (S)\|]<∞. The computation is f(x)≜π[φ(S)],f(x)\; \;E_π[ (S)], (2) and the output is the time-averaged readout y^T=1T∑t=1Tφ(St), y_T= 1T _t=1^T (S_t), (3) which converges almost surely to f(x)f(x) as T→∞T→∞ whenever ℒL generates ergodic dynamics [48], regardless of initial state. To make the tuple concrete, take the linear system Ax=bAx=b of Section˜2.1. Here the state space S is that of the Ornstein–Uhlenbeck process, the generator ℒL is the operator of the dynamics dSt=−(ASt−b)dt+2dWtdS_t=-(AS_t-b)\,dt+ 2\,dW_t with A and b encoding the input, and the stationary distribution is the Gaussian π=(A−1b,A−1)π=N(A^-1b,\,A^-1). Taking the identity readout φ(S)=S (S)=S gives π[φ(S)]=A−1bE_π[ (S)]=A^-1b, so the time-average y^T y_T of the state is the running estimate of the solution, and reading the sample covariance instead recovers the inverse A−1A^-1. What the definition requires is only that the dynamics equilibrate to a unique stationary distribution. The vocabulary of energies and temperatures motivates the name and supplies the intuition, but it plays no role in the definition itself, and the generator need not correspond to any physical system. This is what makes the definition substrate-independent. The realisation of Section˜2.1 took ℒL to be a Langevin operator on a continuous state space, which requires continuous noise and so confines them to analogue hardware [46]. The definition given here, however, does not impose such a restriction: ℒL may be any ergodic generator of which the stationary distribution encodes the computation, including a discrete-state Markov chain realised entirely in digital logic (Fig.˜2). 02244π Langevin (continuous in time and value)02244π → Markov chain (discrete in time and value) Figure 2: Different forms of Equilibrium Thermodynamic Computing. Different stochastic processes can relax to the same π, and so compute the same f(x)=π[φ(S)]f(x)=E_π[ (S)] under the equilibrium definition of thermodynamic computing. Top: an analogue Langevin process wanders continuously, in both time and value. Bottom: a digital Markov chain takes discrete steps between discrete states. 2.2.1. Challenges Realising this computation for a given f in a hardware-native approach requires solving three coupled challenges related to the encoding, implementation, and efficiency. Encoding: The encoding problem is to find a generator ℒL and readout φ such that π[φ(S)]=f(x)E_π[ (S)]=f(x), demanding co-design of all three objects: generator, observable, and input encoding. For TLA this is known analytically: an Ornstein–Uhlenbeck process with A encoded in the drift has a Gaussian stationary distribution with mean A−1bA^-1b, and the identity readout suffices [3]. For general f, such as the input-output map of a neural network, no trivial recipe exists. Implementation: Beyond encoding, the generator must also be implementable as the native dynamics of a physical or digital substrate, not as a simulation running on top of a general-purpose processor. The SPU is the canonical example: its RLC circuit does not simulate the Ornstein–Uhlenbeck process but is that process, with the dynamics carried out by the device’s physics rather than by instructions on a host. This is what distinguishes thermodynamic computing from running Markov chain Monte Carlo on a GPU: on a GPU the substrate bears no structural relationship to the generator, whereas in thermodynamic computing the structure of the substrate is the structure of the generator. Efficiency: Even with the correct stationary distribution and a realisable generator, the dynamics must also reach and then explore that distribution fast enough for y^T y_T to converge in useful time. This is governed by two principles, burn-in and mixing time, as familiar from Markov chain Monte Carlo [8]. The first is burn-in: the chain must reach its stationary distribution at all, and the number of steps it takes to approach π is the relaxation time τrelax _relax. Below τrelax _relax the readout carries a systematic bias from the initial state that decays only as T grows past τrelax _relax. The second is mixing: once at equilibrium, successive samples along a trajectory must decorrelate for average to improve and reduce variance. The rate of this is set by the correlation time τcorr _corr. A hardware implementation that encodes the right distribution but either equilibrates or mixes slowly, with τrelax _relax or τcorr _corr large, is uncompetitive. 2.3. Scaling Axes of Thermodynamic Computing Because every output is a statistical average over a stochastic trajectory rather than the result of a fixed instruction sequence, an equilibration-based computation can be scaled along axes that have no clean analogue on conventional accelerators. We describe three. The first trades time for accuracy within a single run, the second trades hardware for accuracy across many runs, and the third removes the wall-clock cost of composing dependent computations. All three follow directly from the statistical nature of the readout, and each maps onto a concrete lever for scheduling work on real hardware. 2.3.1. Anytime precision. The readout y^T y_T is an average over T samples, and once past the burn-in for T beyond τrelax _relax, its error falls as the (T−1/2)O(T^-1/2) of the mixing regime [48]. Precision is therefore not a fixed property of the hardware but a knob set by how long the process is run past the initial burn-in: the computation returns a valid, if coarse, answer at every T and refines it monotonically in expectation as T grows [67]. A workload can be run only until it reaches the accuracy the downstream task requires and then terminated, rather than to some hardware-fixed bit width, so the energy that would have bought further precision is simply never spent. In a generative setting the same lever permits early exit: a draft can be inspected as it sharpens and abandoned before it is fully formed if it is not worth completing. Unlike a fixed-latency accelerator, which must finish all steps in computation before producing any output, an equilibration-based computation makes precision, time, and energy a tunable trade-off. 2.3.2. Parallel sample aggregation. A single run carries an intrinsic cost in the variance of the readout from one replica set by the variance of the observable under π and the correlation time τcorr _corr of the generator, Var(y^T)≈τcorrVarπ(φ)/TVar( y_T)≈ _corr\,Var_π( )\,/\,T [48]. How quickly this base variance falls is precisely the efficiency challenge of Section˜2.2.1. Fortunately, parallel sample aggregation provides a lever here, allowing one to trade off variance with extra hardware commitment. This lever comes from the observation that trajectories produced by independent runs of the same dynamics are statistically independent, so K replicas of the hardware running in parallel produce a combined estimator whose variance is K times smaller than a single replica of the same wall-clock duration, with no inter-replica coordination beyond a final averaging step [8]. The aggregate estimator y¯T(K)=1K∑k=1Ky^T(k) y_T^(K)= 1K _k=1^K y_T^(k) (4) satisfies Var(y¯T(K))=Var(y^T(1))/KVar( y_T^(K))=Var( y_T^(1))/K, so the rate of variance reduction scales with the aggregate throughput KTKT rather than either factor alone. This makes the precision of a computation a matter of resource allocation. In a datacentre, more chips or more power can be directed at a single latency-critical computation to reach a target accuracy sooner, and withdrawn again when that additional computation is no longer needed, with the variance of the result tracking the resources assigned to it. It needs to be noted here that there is a limit to what more hardware buys. Aggregation drives down the variance but leaves the burn-in bias untouched. At a fixed per-replica window T this can be averaged down only as far as the systematic floor that the burn-in leaves behind. Past that floor, precision is recovered not by adding replicas but by rolling each of them out beyond τrelax _relax. 2.3.3. Sequential parallelism. Of the three axes, this is perhaps the most consequential for generative AI, where the workloads of interest are long chains of dependent stages. Suppose the target function is realised not as a single thermodynamic computation but as a composition f=fL∘⋯∘f1f=f_L ·s f_1 of L stages, each fℓf_ itself a thermodynamic computation in the sense of Section˜2.2, with its own running estimate y^T(ℓ) y_T^( ). Each such estimate is available at every T and is a consistent estimator of its stage’s stationary expectation, as described by the anytime precision property. As such, a downstream stage fℓf_ can consume the running estimate y^T(ℓ−1) y_T^( -1) of the stage feeding it before that estimate has converged: early estimates are imprecise and improve over time, and as long as fℓf_ is continuous in its input, the downstream stage’s output tracks the converging upstream estimate and remains consistent with the composed function f. The continuous-mapping theorem [63] states that if Xn→X_n→ X almost surely and g is continuous at X, then g(Xn)→g(X)g(X_n)→ g(X) almost surely. Applied along an L-stage pipeline in which each stage stays consistent as the estimate feeding it converges, it gives y^T(L)→a.s.(fL∘⋯∘f1)(x)as T→∞, y_T^(L)\; a.s.\;(f_L ·s f_1)(x) T→∞, (5) provided each fℓf_ is continuous, regardless of whether the stages are run concurrently or sequentially. Viewed across all L stages at once, the pipeline is a single Markov chain over the joint state =(S(1),…,S(L))S=(S^(1),…,S^(L)) formed by concatenating the states of the L stages. Each stage advances its own component every cycle on the still-converging estimates of the stages feeding it, in the manner of a Jacobi-style fixed-point iteration over the whole unrolled chain. This is similar to the parallel-in-time relaxation of the Parareal family of evolution-equation solvers [41], which iterate the entire time domain to a fixed point rather than marching step by step, and the same Picard-style parallelisation has recently been shown to draw samples from diffusion models in far fewer sequential rounds than their nominal step count [55]. What relaxes is this joint state, and the end-to-end wall-clock is set by the time the joint chain takes to settle rather than by the sum of the per-stage latencies. The contrast with conventional execution is the essential point. In a standard pipeline each stage must produce a finished, exact output before the next can consume it, so an L-stage computation costs the sum of the per-stage latencies. Here the stages settle together rather than in turn, so that at the accuracy the output requires the joint relaxation time can fall well below L times a single stage. This is not to say that depth comes for free. Information must still propagate through the stages, so the advantage narrows as the target accuracy tightens. 3. CN101: A Digital Thermodynamic Computer Having established a formalism for equilibrium thermodynamic computing that is substrate-independent, we now get to the second core contribution of this work: realising it in digital silicon through the test chip CN101. The functions the chip is built to evaluate are deep nested compositions of affine maps and elementwise nonlinearities, x(L)=σ(W(L)σ(⋯σ(W(1)x+b(1))⋯)+b(L)), x^(L)=σ\! (W^(L)\,σ (·s\,σ(W^(1)x+b^(1))\,·s )+b^(L) ), (6) a class we return to in Section˜5. This nesting and how the output of each layer is the input of the next, such that execution runs as a deep sequential path, is the core property we consider for the chip’s functional design. 3.1. Architecture CN101 is the first silicon instantiation of the Carnot Architecture, a general design for realising the equilibration formalism of thermodynamic computing in digital logic through stochastic computing. In the following section we will discuss the basics of stochastic computing, the tile architecture representing individual stages of computation, how tiles are connected using a Stochastic Streaming Network-On-Chip (SSNoC) to compose stages, tighten the connection of the architecture and the equilibration formalism, and finally discuss how multiple chips can be combined to handle more compute intensive workloads. 3.1.1. Stochastic computing Rather than the binary word of a fixed-point representation, in which each bit carries a fixed place value, stochastic computing encodes a value in a random bitstream. A value x is carried by a sequence of independent Bernoulli bits b1,b2,…b_1,b_2,… with [bt]=xE[b_t]=x, and is recovered as the time-average x^T=1T∑t=1Tbt,x^T→a.s.x,Var(x^T)=x(1−x)T. x_T= 1T _t=1^Tb_t, x_T a.s.x, ( x_T)= x(1-x)T. (7) The encoded value is therefore an expectation and the readout an unbiased estimate of it. Reading a value as the time-average of an accumulator driven by a stochastic bitstream is the same stochastic-integration mechanism prior stochastic-computing hardware uses to solve differential equations [44, 43]. Precision, similar to the thermodynamic computing framework outlined before, is set by the length of the averaging window rather than by a fixed word width (Fig.˜3). Every bit in the stream carries equal weight, unlike the graded place values of a fixed-point or floating-point word, so that reaching a precision ϵε costs of order ϵ−2ε^-2 bits. As such, the representation is cheap when a coarse estimate suffices and expensive when an exact one is demanded. Similar to the anytime precision property discussed earlier. Specifically, CN101 supports unipolar as well as split-unipolar encoding, the latter of which is used for the experiments reported in the remainder of this paper. In split-unipolar encoding a single unipolar stream represents a value in [0,1][0,1] and a signed value is carried by a pair of such streams, so that signed arithmetic reduces to unipolar arithmetic on two channels [52]. (a)Fixed-point: each bit a fixed place value101112 1214 1418 18116 116= 0.6875=\;0.6875Stochastic: value is the fraction of onesx^T=612=0.5 x_T= 612=0.5(b)andppqqmean pqpq multiply, 0.5∧0.4→0.200.5 0.4→ 0.20orp1p_1p2p_2mean 1−∏i(1−pi)1- _i(1-p_i) accumulate Figure 3: Fixed-point versus stochastic representation. (a) A fixed-point word assigns each bit a fixed place value, whereas a stochastic bitstream encodes a value as its fraction of ones, recovered as the time-average x^T x_T. (b) In this representation arithmetic reduces to single logic gates: an and gate multiplies two independent streams to mean pqpq (worked example 0.5∧0.4→0.200.5 0.4→ 0.20), and an or gate accumulates a column of products. Note that this saturates to mean 1−∏i(1−pi)1- _i(1-p_i), which has to be corrected at readout. For two independent unipolar streams with means p and q, the bitwise and has mean pqpq, so that multiplication is a single and gate rather than a full multiplier. Similarly, summation is implemented using an or gate [22, 5]. As such, a matrix–vector product can be implemented efficiently using a grid of and gates feeding column-wise or accumulation. This forms the basis of the stochastic matrix multiplier (SMM) (Fig.˜4d) at the core of each tile. This and-multiply, or-accumulate datapath is often used in stochastic-computing neural-network accelerators [51, 34, 56]. Specifically, for the SMM on CN101, the weight operand W is held as a deterministic 8-bit value rather than as a second stream, so that only the input x adds variance. Each input bit is gated against the eight weight bits separately, so that every bit plane of W carries its own or accumulation column, and the eight columns are recombined by place value at readout. The and/or datapath is in this way kept intact while the weight retains its precision, which is why the array delivers a multi-bit result. (a) (b) (c) (d) Figure 4: CN101 stochastic digital thermodynamic computing. (a) Top-level architecture. Four tiles (each comprising an SMM, RNB, and OA) are connected by the SSNoC. The IO block provides host configuration access. (b) Die photograph of CN101. (c) Polysynchronous clocking. Each RNB is driven by either a global clock tree or a local ring oscillator, selected by a clock mux; two flip-flops handle metastability at the clock domain boundary. (d) OR-based accumulation in the stochastic matrix multiplier with AND gates computing bitstream products per element. 3.1.2. The tile One SMM together with a reconfigurable neuron bank (RNB) form the core computation unit on CN101, tiles. Each tile evaluates a single layer σ(Wx+b)σ(Wx+b) in full. The term tile is used in this specific sense throughout: a tile is a complete affine-plus-nonlinear step, not a sub-block of one matrix multiplication partitioned for data locality, as tile denotes on a GPU. The SMM stores the weight matrix W and forms WxWx through the and/or arithmetic above. The RNB adds the bias b and applies the nonlinearity σ. Each CN101 contains 4 of these tiles (Fig.˜4a). The nonlinearity is implemented along one of two paths. Where possible, we rely primarily on ReLU activations, which are trivial in the split-unipolar encoding since they act separately on the two sign channels. However, for generic nonlinearities the RNB also has access to a finite-state machine, the standard construction for nonlinear functions in stochastic computing [9, 40]. Within a tile the SMM passes its multi-bit result to the RNB over a direct, deterministic path, so that applying a layer adds no sampling noise of its own. The inputs x are themselves stochastic estimates, so that the tile’s output is still a stochastic quantity, with variance inherited from its inputs. What the deterministic intra-tile path avoids is the additional variance that re-encoding the value into a fresh bitstream would introduce. Beyond the multiplier and the neuron bank, each tile carries an output accumulator (OA), a readout channel whose role in composition we take up in Section˜3.3. 3.1.3. Composing tiles With a single tile evaluating a single layer, a deeper model can be assembled by composing multiple tiles. This is achieved by passing the output of one tile to the input of the next. CN101 carries this composition over its stochastic streaming network-on-chip (SSNoC), a reconfigurable Network-On-Chip that routes bitstreams between tiles. The bitstream representation required is generated within the RNB itself. Here each neuron carries a 65-bit linear-feedback shift register (LFSR) with a nonlinear output stage, of period exceeding 101910^19, so that this per-neuron pseudo-random number generator (PRNG) supplies the Bernoulli bits for that neuron’s streams [49]. As such, there exist two communication paths with complementary roles in the chip. Inside a tile, the SMM-to-RNB path is deterministic and multi-bit, so that a single layer is applied without added sampling noise. Conversely, between tiles, the SSNoC path carries bitstreams, is reconfigurable, and crosses between independent clock domains. The first preserves the accuracy of a layer and the second is what makes the composition of the chip into distinct equilibration units. One specific property that stochastic communication between tiles allows is a clean separation between clock domains that requires no synchronisation, resulting in a globally-asynchronous locally-synchronous design [11] at the level of a tile. A bitstream carries an expected value, and every bit of it carries the same expected value, so that the quantity a downstream tile estimates does not depend on which bits arrive or on exactly when they arrive. A delayed or dropped bit changes the number of samples averaged, not the value being averaged. As such, the SSNoC can move streams between tiles running on independent clocks without a global clock, each tile driven either by a global clock tree or by a local ring oscillator and left to settle at its own rate. Table˜1 lists the physical parameters of CN101 (Fig.˜4b). Table 1: Key physical parameters of CN101. Parameter Value Tile count 4 SMM array (per tile) 64×6464× 64 Neurons per RNB 64 Total MAC cells 16,384 Total neuron accumulators 256 Accumulator width 32 bits PRNG 65-bit LFSR + nonlinear output stage Operating frequency ≈ 500 MHz 3.2. CN101 and the equilibration formalism Having described how CN101 evaluates a function, we now make its connection to the equilibration formalism precise. Crucially, the chip does not compute its outputs in closed form: each value it carries is the expectation of a bitstream, obtained as a time-average over a window of T cycles, and the bitstreams passing through the active tiles together form a single ergodic process. The quantity the chip reports is the long-run behaviour of that process rather than the state it occupies on any one cycle, so that the chip realises a function in the precise sense of Section˜2, as the stationary expectation of its own dynamics. This joint relaxation is distinct from ancestral sampling. The tiles do not draw a finished sample x(ℓ)x^( ) from a conditional p(x(ℓ)∣x(ℓ−1))p(x^( ) x^( -1)) and hand it downstream. Instead, each updates every cycle on the current, still-converging estimate of the tile feeding it, realising the joint fixed-point relaxation of Section˜2.3 [41, 55]. What the chip reads out is the mean of that process, the layer output σ(Wx+b)σ(Wx+b) recovered as a time-average, and never a single draw from it. The four ingredients of the formalism each have a counterpart on the chip. The state StS_t is the joint configuration of the accumulators across the active tiles, an integer vector St∈ℤdS_t ^d whose dimension d is the number of active neurons in the RNBs. The generator ℒL advances this state by one cycle under the pseudo-random bits the neuron banks supply, and because the state space is discrete it is a Markov chain rather than the continuous Langevin operator that Section˜2 admits as its analogue counterpart. The weights, biases, and activations loaded onto the chip fix the stationary distribution π satisfying ℒπ=0Lπ=0, while the readout φ is taken from the output stream, whose time-average y^T=1T∑t=1Tφ(St) y_T= 1T _t=1^T (S_t) (8) converges to π[φ(S)]E_π[ (S)] as the window grows. As such, the function the chip computes is f(x)=π[φ(S)]f(x)=E_π[ (S)], the definition given in Section˜2. The same correspondence shows how CN101 meets the three challenges of Section˜2.2.1. Implementation is the most direct: the accumulator dynamics are the generator, carried out by the chip’s own logic rather than simulated on a host. Encoding is solved by the construction of tiles: instead of searching for a single generator whose stationary expectation is an arbitrary f, CN101 builds f layer by layer, each tile realising one affine-plus-nonlinear step as the stationary expectation of its own stream, with the weights supplied by ordinary training. What remains is efficiency, the requirement that the dynamics mix fast enough for y^T y_T to converge in useful time; this is a quantitative property of the chip rather than a question of principle, and we characterise it directly in Section˜4. 3.3. Composition across tiles and chips To realise more complex functions, we now turn to how the design of CN101 supports composition over multiple chips. Both directions of scaling, depth and width, exploit the same feature of how a tile computes: a tile does not deliver a value at a fixed moment, but produces an estimate that sharpens the longer it runs. A composition of tiles then depends only on these estimates, and not on the tiles sharing a synchronised state, with two consequences. First, tiles may run on independent clocks, since combining them requires no shared notion of when a value is ready. Second, and more important for scaling, a downstream tile may consume an upstream estimate that has not yet converged, so that all stages advance together rather than in turn. Taken together, these properties make the chip boundary a matter of implementation rather than of principle. Within a chip the tiles couple directly over the SSNoC, exchanging streams as they relax on independent clocks. Across a chip boundary the same coupling is realised by exchanging their running estimates, so that many chips act as a single system of tiles converging in unison to the collective result (Section˜B.1). We examine this scaling along the two axes in turn, each applying across chips as readily as within one. Depth: composing chips in sequence. The first axis adds depth by spreading a deep model across several chips. Writing the model as the composition f=fL∘⋯∘f1f=f_L ·s f_1 of its layers, we partition the layers into contiguous blocks, one block to a chip, so that the output of one chip is the input of the next. Each chip holds a running estimate of its block’s output, which the chip downstream consumes. That the chain converges to the intended result is the sequential-parallelism property of Section˜2.3: provided each layer map fℓf_ is continuous in its input, y^T(L)→T→∞(fL∘⋯∘f1)(x). y_T^(L)\; [T→∞]\; (f_L ·s f_1 )(x). (9) Width: composing tiles into a larger matrix. The second axis adds width by combining tiles, on one chip or across several, whose weight matrices together represent a matrix larger than a single 64×6464× 64 array can hold. This is similar to the tiling used for dense linear algebra on a GPU. Partitioning the weight matrix into blocks WijW_ij and the input correspondingly into blocks xjx_j, each tile evaluates one block product WijxjW_ij\,x_j, and the i-th block of the layer output is the sum (Wx)i=∑jWijxj (Wx)_i= _jW_ij\,x_j (10) of the partial products of the tiles assigned to it. Crucially, the tiles equilibrate independently, so that at any instant their partial products have converged to different degrees. This, however, does not bias the result, since each is an estimate of its block product and, expectation being linear, their sum is an estimate of (Wx)i(Wx)_i whatever the individual states. The unevenness in convergence across tiles enters only the variance of the combined estimate and not its expectation. 3.3.1. Output Accumulator Both axes rest on the same concept: a converging estimate is preserved under the summation that builds width and under the continuous composition that builds depth. What this earns in practice is freedom in how, and how often, the estimates are read from the tiles and passed on. To support this communication across chips each tile carries an Output Accumulator (OA). Crucially, the OA performs no part of the tile’s computation; it is a dedicated, non-disruptive channel that taps the tile’s output as the stream passes onto the SSNoC and integrates it into a running sum. The OA allows for the current time-averaged estimate to be read at any moment, without halting the neuron bank or perturbing the dynamics that generate the stream. A downstream stage, located on a different chip, can then draw its input from this accumulated value rather than from the live neuron bank. The reads from these accumulators can be staggered across tiles and taken at different rates. Because a downstream tile tolerates a stale input, a bounded delay between when an estimate is produced and when it is consumed changes only the transient accuracy of the downstream stage, and not the value it converges to. In place of a global barrier, each downstream input is refreshed over a window of cycles from the current reading of its source, so that no two inputs are updated at the same instant and no stage waits on another. Within a chip this exchange is continuous, the streams flowing between tiles as they run. Across chips it proceeds window by window, the window length setting how closely the assembly approaches the fully concurrent limit. 4. Single-chip validation: Digital TLA and Variational Autoencoders Having described how CN101 evaluates and composes functions, we now validate that it does so correctly using a single chip. We use two different workloads for this. The first is a linear system, the problem class for which the equilibration formalism was originally written down and on which the chip’s output can be checked against an exact solution. The second is a small generative model, a conditional variational autoencoder, whose output is checked against a floating-point reference of the same trained network while it executes the full multi-layer cascade. Figure 5: Thermodynamic linear algebra on CN101. (a) Parity plot of CN101 y^T y_T vs analytic A−1bA^-1b across dimensions d∈8,16,32,64d∈\8,16,32,64\ on a single tile and d=128d=128 using all four tiles in multi-tile composition (shown as ×). All points lie on the diagonal. (b) Variance (top) and relative ℓ2 _2 error (bottom) as functions of T, with T−1T^-1 and T−1/2T^-1/2 reference lines. Each curve averages 16 trajectories (4 system matrices × 4 PRNG seeds per matrix) per dimension. (c) Ground-truth (left) and CN101-reconstructed (right) A−1A^-1 at d=16d=16 for κ∈2,50,100κ∈\2,50,100\. 4.1. Digital thermodynamic linear algebra We begin with linear systems Ax=bAx=b with symmetric positive-definite A, realised on the chip as thermodynamic linear algebra [3, 46], for which the encoding is analytically controlled. It is a deliberately narrow validation rather than a representative run: it checks correctness against a known answer on a different datapath from the stochastic-computing equilibration the generative models use. The correct stationary distribution is known exactly, so any deviation between the chip’s output and A−1bA^-1b would directly be attributable to implementation error. For this workload every tile’s activation is set to linear integration, which implements the lattice random-walk discretisation of the Ornstein–Uhlenbeck dynamics [19, 47]: integer-valued accumulator updates whose stationary expectation is A−1bA^-1b to O(Δt)O( t). The weight matrices encode A and b rather than a trained network. Everything else (PRNG seeds, routing tables, and the window length T) is configured as for any other workload. 4.1.1. Results We characterise three properties in turn: correctness of the solution, T−1T^-1 variance scaling, and the operating envelope set by problem conditioning. Correctness. Fig.˜5a shows a parity plot of every component of the chip’s time-averaged readout y^T y_T at maximum T against the corresponding component of the analytic solution A−1bA^-1b, across dimensions d∈8,16,32,64d∈\8,16,32,64\ on a single tile and d=128d=128 using all four tiles composed via the SSNoC. All points lie on the diagonal across the full dynamic range. The multi-tile result is indistinguishable from the single-tile results, confirming that the SSNoC composition preserves the formulation without introducing systematic error. Convergence rate. Fig.˜5b shows variance and relative ℓ2 _2 error as functions of T over six orders of magnitude. The top panel plots Vars(y^T)Var_s( y_T) across 16 trajectories per dimension (4 system matrices × 4 independent PRNG seeds per matrix); the bottom panel plots the relative error for a representative seed. Both panels show a consistent three-regime structure across all dimensions. In the early regime (T≲105T 10^5), variance decays at slope ≈−1≈-1, matching the T−1T^-1 prediction, while the error decreases slowly as the chain is still shedding memory of its initial state; for T≲τrelaxT _relax this is the burn-in of Section˜2.2.1, in which the initial-state bias rather than the variance dominates the error. In the mid regime (105≲T≲10710^5 T 10^7) the variance decay visibly slows as the dynamics enter a transitional mixing phase; correspondingly, the error curve steepens towards its T−1/2T^-1/2 descent as variance becomes the dominant term in the MSE. In the late regime (T≳107T 10^7), variance resumes its T−1T^-1 slope before all curves flatten onto the bias floor at relative error ≈10−2≈ 10^-2, consistent with finite weight-register precision in the SMM. The convergence curves across dimensions are nearly coincident, indicating that the per-tile correlation time τcorr _corr, which sets the rate of the variance-limited mixing regime, is not strongly sensitive to problem dimension in this regime. Conditioning. Fig.˜5c shows the reconstructed inverse A−1A^-1 at d=16d=16 for three condition numbers, assembled column-by-column. At κ=2κ=2 the reconstruction is visually indistinguishable from the ground truth (relative error 0.009). At κ=50κ=50 the structure is accurately recovered (relative error 0.073). At κ=100κ=100 the main structure remains visible with some residual error (relative error 0.180), establishing the operating envelope of the current generation at this problem size. Figure 6: Conditional MNIST generation on a single chip. (a) Reconstructions of all ten MNIST classes at the concurrent operating point (T=500kT=500k cycles). Top row, CPU floating-point reference, bottom row, CN101. (b) Progression of a single digit’s reconstruction as the cycle budget T increases from 1010 to 500k500k; values below each panel are the scale-invariant RMSE against the CPU reference. (c) RMSE vs. T. Solid line: single-chip result at the per-T optimal operating point. Shaded lines: mean of K independent readouts (distinct RNG seeds), K∈2,4,8,16K∈\2,4,8,16\. Dashed line: bias floor ϵbiasVAE _bias^VAE. (d) Cycles required to reach each target RMSE τ for the sequential (∑iTi _iT_i) and concurrent (T) configurations. 4.2. Conditional MNIST generation with a single-chip VAE The second validation workload is a conditional variational autoencoder (VAE) [35, 59]. A conditional VAE on 8×88×8 MNIST fits on a single CN101, its four decoder layers mapping one-to-one onto the four tiles, and as such provides the within-chip demonstration of parallel sample aggregation and sequential parallelism. Unlike the linear system, the intended output here is not known in closed form; we validate instead against a floating-point reference of the same trained model. 4.2.1. Model architecture and training The decoder uses four fully-connected layers of width 64. The first layer concatenates a 16-dimensional latent code with a 10-way one-hot class label and projects to a 64-dimensional hidden representation. Two further 64×6464×64 hidden layers with ReLU activations follow, and the output layer projects to a 64-dimensional image vector reshaped to 8×88×8. The network is trained on 8×88×8 MNIST-like handwritten digits captured natively at 8×88×8 [6] rather than downsampled from the 28×2828×28 MNIST [38] they resemble, using the reparametrisation trick and a class-conditional evidence lower bound (ELBO) objective; training runs in floating point on a CPU and the weights are then loaded onto CN101. The four layers map one-to-one onto the four tiles: one weight matrix per tile, with the per-tile activation register set to ReLU for the three hidden layers and to LINEAR for the output layer. Two execution configurations are compared: concurrent and sequential. In the concurrent configuration all four tiles execute simultaneously and bitstreams flow between them through the network-on-chip. In the sequential configuration each layer runs in isolation to convergence and the integer output of each stage is re-encoded and fed into the next. The final tile’s output is the stream consumed by the output accumulator. 4.2.2. Results Fig.˜6a shows CN101 reconstructions of all ten digit classes alongside the floating-point reference at the concurrent operating point (T=500,000T=500,000 cycles per tile; hyperparameters in Appendix˜C). The reconstructed digits are recognisable across all ten classes with the scale-invariant RMSE averaged across 10 digits × 10 latent samples being 0.083±0.0280.083± 0.028. Fig.˜6b shows the evolution of a single digit’s reconstruction as T increases: at T=10T=10 the output is dominated by quantisation noise; the digit is clearly recognisable by T=10,000T=10,000 and converges onto the T=500,000T=500,000 limit by T≈105T≈ 10^5. Anytime-precision. Fig.˜6c plots the RMSE as a function of T, with per-T operating points selected by a search over the chip’s scale and precision knobs (details in Appendix˜C). The curve drops from RMSE≈0.79RMSE≈ 0.79 at T=100T=100 to a bias floor ϵbiasVAE≈0.09 _bias^VAE≈ 0.09 at T=106T=10^6, set by weight quantisation in the multiplier register banks and independent of T. At this operating point the reconstructed digits are perceptually indistinguishable from the floating-point reference at the output resolution of 8×88× 8 pixels. Reducing the floor further would require higher-precision weight registers. This sweep over T is an example of the anytime-precision axis discussed in Section˜2.3 on a generative workload. Parallel sample aggregation. Fig.˜6c overlays four curves obtained by averaging across K∈2,4,8,16K∈\2,4,8,16\ independent RNG seeds at each T. The K curves separate, approaching the 1/K1/ K spacing of the variance-limited regime (at T=1000T=1000 the 1→161→16 seed ratio is 2.3×2.3×, short of the ideal 4×4× because the bias floor still contributes). At small T the burn-in transient dominates the error and, being shared by replicas started alike, is not reduced by averaging (at T=100T=100 the ratio is 1.2×1.2×). At large T the bias floor limits further convergence. Sequential parallelism. Sequential parallelism Eq.˜5 predicts that the concurrent configuration delivers the same output as the sequential configuration, but with a joint relaxation time rather than the sum ∑iTi _iT_i of the per-stage times. Fig.˜6d quantifies this at four RMSE targets τ∈0.25,0.20,0.15,0.10τ∈\0.25,0.20,0.15,0.10\. For each target the smallest cycle budget is found for both execution forms: concurrent and sequential (search details in Appendix˜C). The concurrent budget grows from 7,3007,300 cycles at RMSE=0.25RMSE=0.25 to 175,300175,300 cycles at RMSE=0.10RMSE=0.10. The sequential total ∑iTi _iT_i tracks the same slope but is larger by a factor of 55–11×11×; at RMSE=0.10RMSE=0.10 the concurrent configuration reaches the target in 1.75×1051.75× 10^5 cycles against 8.8×1058.8× 10^5 cycles for sequential. 5. Thermodynamic computing for generative AI Having validated the single-chip operation of CN101, we now focus on its core target workload: generative AI. Modern Generative AI approaches are inherently sequential in nature and therefore benefit most from the scaling axis of thermodynamic computing implemented on CN101. Progress in deep learning has repeatedly come from depth rather than width: the breakthrough on ImageNet came not from wider networks but from very deep residual networks that made optimisation at depth stable [24]. The generative models that followed inherited this character by construction. Diffusion and flow-matching models are defined as continuous-time processes and generated by integrating a long chain of dependent steps, and large language models generate one token at a time, each conditioned on every token before it. Sequential depth is therefore not incidental to these models; it is where their capability comes from. The standard hardware these models run on, however, rewards the opposite. The parallelism of the GPU favours width over depth, larger batches and more parallel lanes over long dependent chains, and this “hardware lottery” [27] has shaped a decade of model design towards what existing accelerators execute cheaply, at the expense of the sequential depth on which the most capable generative models depend. Models have been made wider and shallower, and batched ever more aggressively, to suit the hardware rather than the structure of the problem. The scaling axes of Section˜2.3, and sequential parallelism in particular, invert this trade-off. Because a downstream stage can consume an upstream estimate while it is still converging, dependent stages relax concurrently rather than in sequence. Where the intermediate stages tolerate imprecision a deep inherently sequential computation can complete in close to the time of a single stage rather than in proportion to its depth, as we will experimentally validate in the coming sections. As such, thermodynamic computing matches the sequential structure of modern generative models rather than working against it. To be able to demonstrate this, each workload in this section runs across six CN101 chips that realise the multi-chip composition of Section˜3.3, each holding a block of the model and advancing on the still-converging estimates of the blocks upstream of it. The present chip is a prototype, and composing six of them is in part how we assemble the weight capacity these models need; but it is equally a benchmark of the composition itself, which is the mechanism by which larger models are assembled on future hardware. 5.1. Generative models as deep nested maps Before continuing, we now first make explicit that modern generative methods, such as diffusion [25, 58] and flow matching, can be described as a deep nesting of affine maps and nonlinearities of the form σ(Wx+b)σ(Wx+b) and therefore map well to the core operating model of thermodynamic computing. A diffusion or flow-matching model generates by integrating an ordinary differential equation x˙=vθ(x,t) x=v_θ(x,t) (or equivalent stochastic differential equation [20]). Its Euler discretisation xk+1=xk+Δtvθ(xk,tk) x_k+1=x_k+ t\,v_θ(x_k,t_k) (11) is a residual connection, so that the unrolled integration is a deep residual network in which the learned field vθv_θ is the per-step layer [24, 13, 21]. The number of these steps is the Number of Function Evaluations (NFE), which we denote N. A network of L layers unrolled over N steps is a nested computation L×NL× N layers deep, and it is this product, not L alone, that sets the sequential length of generation. The stochastic case extends this without leaving the nested form. A discretised stochastic differential equation, whether the noise enters additively, multiplicatively, or through the argument of the drift, takes the Euler–Maruyama form [36] xk+1=xk+f(xk,ΔWk),ΔWk∼(0,Δt). x_k+1=x_k+f(x_k,\, W_k), W_k (0, t). (12) Once the noise increments ΔWk W_k and the step schedule are fixed, the unrolled integration is again a deterministic deep nested map, differing from the ODE case only in the injected randomness [60]. 5.2. CIFAR-10 generation with a convolutional flow Flow matching [42] is one of the core methods of modern generative AI. It trains a velocity field so that integrating it transports a simple prior into the data distribution. Flow models built this way underlie much of image, video and molecular generation. It is therefore our first benchmark of what thermodynamic computing offers generative AI. We train a convolutional flow-matching model to generate CIFAR-10 images [37] and run it across the CN101 chips, asking two things: whether the substrate reproduces the same results in expectation as when run on a traditional accelerator, and how the thermodynamic computing scaling axes act on the sequential structure of the generation. 5.2.1. Setup The prototype chip’s only operation is a 64×6464× 64 matrix multiply with no direct native support for convolutions. As such, we express each convolution as a host-side shift [66] followed by an on-chip 1×11× 1 convolution. The shift re-indexes each channel’s spatial map so that a pixel’s neighbours are gathered into its own channel vector, a free memory operation, after which the 1×11× 1 convolution, a per-pixel Linear(C,C)Linear(C,C), is exactly the chip’s matmul, streamed over every pixel of the image (weight-stationary). The group normalisations and nonlinearities run on the host. The resulting model, ChipUNetTiny, is a ten-block U-Net [53] (a stem, two downsample–encoder stages, a bottleneck, two upsample–decoder stages with skip connections, and a head) on 32×32×332× 32× 3 images. Its ten 1×11× 1 blocks are pinned as ten weight-stationary tiles across the six-chip system (Section˜3.3), with two further tiles reallocated to share the busiest layers’ pixel streams, so every weight is resident throughout generation. The learned velocity field is integrated over N steps, each a full pass through the ten-layer U-Net, so the unrolled integration is 10×N10× N layers deep. The chip does not give each of those layers its own tile. The ten weight matrices are pinned one to a tile across the six CN101 chips, and because they recur at every integration step those ten tiles are reused, one evaluating each layer that shares its weights. This reuse bounds throughput, but it does not serialise the relaxation: sequential parallelism rests on the consistent-estimate property, so layers that share a tile still relax concurrently from one another’s running estimates. We train the flow-matching model on a set of ten CIFAR-10 images and apply reflow (rectified flow) [45] so that high-quality samples are produced in few integration steps. We use N=10N=10 throughout, so the unrolled integration is 100 layers deep. The model is deliberately narrow in scope: it is overfit to these ten images and does not generalise, the prototype’s on-chip weight memory being too small to hold a general CIFAR model. Every per-step matmul is read out as a time-average over a schedule rising to 5×1055× 10^5 cycles per tile. The first windows, over which the unrolled relaxation fills, are discarded from this average (a burn-in). We compare the chip-generated images against the floating-point reference of the same trained model. Figure 7: CIFAR-10 generation on chip. (a) Four generated samples: the floating-point reference (CPU), the chip output (CN101), and their per-pixel absolute difference (|Δ|| |). (b) Pixel-RMSE of a single generated image against the reference as the cumulative cycle budget grows to 10M cycles per pixel, with the chip image at cumulative budgets of 0.1, 1, 2, 5 and 10M cycles below. (c) Cumulative cycles per pixel to reach the same image quality (RMSE≈0.06RMSE≈ 0.06) for the 100-layer integration. 5.2.2. Results The model generates recognisable CIFAR images, and the chip reproduces the floating-point reference of the same model, sample for sample, at a pixel-RMSE of ≈0.06≈ 0.06, with the residual concentrated in high-frequency detail rather than the composition (Fig.˜7a). A sample begins as quantisation noise and sharpens into a recognisable car as cycles accumulate, the pixel-RMSE falling over the schedule until it settles onto the floor set by the precision of the weight registers (Fig.˜7b). This is an example of the anytime-precision axis discussed in Section˜2.3, realised across multiple chips. Comparing sequential and concurrent execution: run sequentially, each layer must converge before the next can use it, so the cost grows in proportion to the depth. Relaxed concurrently, the layers converge jointly, reaching the same image (pixel-RMSE ≈0.06≈ 0.06) in ≈10≈10M cycles per pixel against the sequential run’s ≈620≈620M, some 62×62× fewer at this target (Fig.˜7c). The sequential cost is a conservative bound: 60 of its 100 layers never meet the per-layer convergence target and exhaust the cycle cap. This points to a property of the sequential parallelism and depth-collapse worth stating in its own right, and to why the sequential baseline is so costly. The concurrent run does not spend cycles making the intermediate layers accurate. At the moment the image first becomes good (pixel-RMSE below 0.1, after only ≈3≈3M of the eventual cycles), the layers in the middle of the network still carry 20–70% per-tile error against the floating-point reference, and they remain inaccurate even at full convergence (Fig.˜8). The output tolerates this because the integration is residual: the image is x0+∑kdtvkx_0+ _kdt\,v_k, so per-step, per-layer errors average out and the residual dominates. Sequential parallelism thus allocates precision automatically to where it reaches the output and spends none where it does not, which is precisely the accuracy the layer-by-layer baseline wastefully demands of every intermediate tile. Figure 8: The chip spends no cycles on precision the output does not need. Per-tile error of the concurrent depth-collapse against the floating-point reference across the unrolled (step×layer) graph, at the window where the generated image first reaches pixel-RMSE <0.1<0.1. 5.3. Free-energy estimation on alanine dipeptide Generative modelling has increasingly turned to the molecular sciences, where diffusion and flow-matching models generate small molecules, protein structures and materials [50]. One problem in this space which has recently gained more interest is that of estimating free-energy difference [26], a thermodynamic quantity. We take this as our final demonstration, on alanine dipeptide, the standard small-molecule benchmark, and find that a thermodynamic substrate is especially well matched to this thermodynamic problem. The slow degrees of freedom of alanine dipeptide are the two backbone dihedrals (ϕ,ψ)(φ,ψ), whose Boltzmann distribution π(x)∝exp(−U(x)/kBT)π(x) (-U(x)/k_BT) concentrates in a small number of metastable basins (Fig.˜9). The quantity of interest is the free-energy difference between two basins i and j, ΔFij=Fj−Fi=−kBTlogZjZi, F_ij=F_j-F_i=-k_BT Z_jZ_i, (13) which fixes their relative equilibrium population through pj/pi=exp(−ΔFij/kBT)p_j/p_i= (- F_ij/k_BT). Free-energy differences of this kind set binding affinities, solubilities and conformational preferences, and their estimation is among the most demanding tasks in computational chemistry, because the basins are separated by barriers that equilibrium molecular dynamics crosses only rarely and so converges slowly. We stress that the task here is free-energy estimation rather than conformational sampling. Figure 9: Chip-generated conformations and the metastable basins. Left: chip-generated alanine-dipeptide conformations, one per metastable basin. Right: density of the reference molecular-dynamics ensemble in the backbone-dihedral plane (ϕ,ψ)(φ,ψ), with the six flat-bottom state-definition boxes that delimit the basins overlaid. Estimating free-energy differences. We estimate ΔF F by targeted free-energy perturbation [31, 30], using a trained Riemannian flow-matching model [42, 12] as the targeting map T that transports configurations from one basin to another. The work of transporting a configuration x0x_0 to x1=(x0)x_1=T(x_0) between basins with potentials UiU_i and UjU_j is W(x0)=Uj((x0))−Ui(x0)−kBTlog|det∇(x0)|,W(x_0)=U_j\! (T(x_0) )-U_i(x_0)-k_BT | (x_0) |, (14) where the log-determinant of the map is the time integral of the divergence of the velocity field vθv_θ along the trajectory [13], log|det∇(x0)|=∫01(∇⋅vθ)(xt,t)dt. | (x_0) |= _0^1 (∇\!· v_θ )(x_t,t)\,dt. (15) We train one such map between αR _R and each of the other five basins, and estimate the free-energy difference of each basin to αR _R with the Bennett acceptance ratio [7], applied to the forward and reverse works of that single pair (Appendix˜E). The remaining pairwise differences follow by subtraction, ΔFij=F^j−F^i F_ij= F_j- F_i. Figure 10: Free-energy estimation on alanine dipeptide. (a) On-chip free energies for each basin relative to αR _R against the umbrella-sampling reference (kJ/mol), with the CPU estimate of the same flows overlaid. (b) The free-energy estimate reconstructed at each cumulative cycle budget per tile, one curve per basin. (c) The running estimate over the chip works (solid) tracks the running estimate over the CPU works (dashed) as samples accumulate, for a tight basin (α′α ), a hard one (C5) and an intermediate one (αL _L). Why sequential parallelism helps. The cost of this estimator is set by the log-determinant Eq.˜15 rather than by the transport [23]. A configuration is carried faithfully in relatively few integration steps, but the divergence integral converges only as the number of steps grows, so it is the density estimation, not the sampling, that sets the depth of the computation. On CN101, using the sequential parallelism property, the divergence trace part of the computation can equilibrate at the same time as the transport rather than in a separate pass. 5.3.1. Setup We estimate ΔF F between all six metastable states of Alanine Dipeptide (αR _R, αL _L, αD _D, β, C5 and α′α ) through a star of flows centred on αR _R [26], so that one reference state anchors every estimate and the remaining pairwise differences follow by subtraction. For each basin we draw 6464 equilibrium configurations from a reference molecular-dynamics trajectory, transport them to αR _R and back in both directions using the learned flows, and combine the resulting forward and reverse works pairwise with the Bennett acceptance ratio. Each transport integrates the flow with a single Hutchinson probe and runs across six CN101 chips, its per-step matmul read out as a time-average over a schedule rising to 10710^7 cycles per tile. We compare against two references. The first is the umbrella-sampling estimate of the same free energies [61, 26], an independent physical ground truth. The second is the CPU estimate of the same flows, which isolates the contribution of the substrate from that of the model. The flows, the reference data and the training procedure are described in Appendix˜E. 5.3.2. Results Across the basins the chip reproduces the CPU estimate of the same flows to within 0.80.8 kJ/mol (Table˜2, Fig.˜10a), and, with C5 a noticeable exception, the free-energy differences track the umbrella-sampling reference. Where the flow is accurate the chip is accurate; the on-chip and CPU estimates agree closely. The transported configurations are themselves valid molecular structures, as illustrated in Fig.˜9. Table 2: On-chip free energies ΔF F (kJ/mol, relative to αR _R) against the CPU estimate of the same flows and the umbrella-sampling reference. Uncertainties on CN101 and CPU are bootstrap standard deviations over the 6464 configurations drawn from each basin. Basin CN101 CPU US ref. αL _L 7.4±0.87.4± 0.8 6.9±0.56.9± 0.5 7.47.4 αD _D 11.8±1.211.8± 1.2 12.5±0.912.5± 0.9 12.112.1 β −1.2±0.8-1.2± 0.8 −1.6±0.5-1.6± 0.5 −1.1-1.1 C5 −0.7±1.2-0.7± 1.2 −0.1±0.9-0.1± 0.9 1.41.4 α′α 6.4±0.56.4± 0.5 6.7±0.26.7± 0.2 6.66.6 Fig.˜10b shows how CN101 carries each value as a time-average, such that the free energy converges with the cycle budget rather than resolving at once. As the schedule rises to 10710^7 cycles per tile the substrate noise averages out and the estimate settles onto its converged value. However, the free energy is an exponential average of the per-sample works [68], so it does not simply average the substrate noise away as an ordinary mean would. Instead, residual noise can bias the estimate. Here that bias is however small and shrinks as the cycle budget grows, and combining forward and reverse transports suppresses it further, so the chip and CPU running estimates track each other as samples accumulate (Fig.˜10c). Discussion In this work we have given a substrate-independent formalisation of equilibration-style thermodynamic computing, and a digital instantiation of it in silicon as the CN101 prototype chip. Across four workloads, spanning linear systems, conditional MNIST generation, a convolutional flow on CIFAR-10, and molecular free-energy estimation, we have shown how the time-averaged readout y^T y_T converges to the intended function as the window grows, which is the anytime-precision axis realised in hardware. We characterise the other two axes directly: parallel sample aggregation on the VAE decoder, and sequential parallelism on both the within-chip VAE decoder and the cross-chip flow-matching pipeline. A recurring feature of the experiments is a small bias floor in the relative error of y^T y_T at large T, where the variance has been averaged away and a residual offset remains. On the linear-algebra workload this floor sits near one percent across the dimensions tested, and on the generative models it inherits a further contribution from the per-layer encoding error that accumulates through the network. An analogous floor was reported on the analogue Stochastic Processing Unit and attributed there to circuit non-idealities [46]. On CN101 the dominant contribution is instead the finite precision of the weight registers. Distinct from this equilibrium floor is a transient, or burn-in, bias. Unlike the variance, this transient is shared across independent replicas started from the same state, so parallel sample aggregation does not remove it. Both the equilibrium floor and this transient are properties of this instantiation rather than of the formulation. The principal open question these results raise is how to design generators that both equilibrate and mix quickly, with short τrelax _relax and τcorr _corr [39], on a digital substrate. One route is suggested by the field whose name the paradigm borrows. Rather than wait for a generator to equilibrate, the tools of stochastic and non-equilibrium thermodynamics drive it through a finite-time protocol and recover the equilibrium answer by reweighting the trajectories with the fluctuation theorems of Jarzynski and Crooks [31, 16, 54]. Finite-time thermodynamics bounds the dissipation such a protocol incurs through the geometry of optimal driving [57, 17], and recent work shows that the protocols can be learned [26]. A learned non-equilibrium protocol is, in our opinion, the most promising route past the mixing-time barrier. CN101 is a prototype, and the first chip in a broader architectural programme aimed at production-scale generative AI. The present paper characterises the behaviour of the formulation on this first generation, while the energy and latency targets that motivate the programme are deferred to later chips, on which the blocks left deliberately first-generation here are to be improved. Two choices in CN101 already anticipate that path. Its modular pseudo-random-number interface accepts a physical noise source in a future chip without a change to the datapath, and its polysynchronous clocking removes the global clock-distribution constraint that would otherwise bound how far a single chip, or an assembly of them, can scale. That a standard generative model, trained and run on GPUs and mapped onto the chip without custom layers, runs correctly on dynamics this far removed from a Langevin system is the result we find most telling: it indicates that thermodynamic computing in its equilibration-style formulation, and the hardware advantages that come with it, are within reach of the models the field already builds. Acknowledgements Normal Computing thanks the Advanced Research and Invention Agency’s (ARIA) Scaling Compute programme for funding this work. References [1] N. A. Aadit, A. Grimaldi, M. Carpentieri, L. Theogarajan, J. M. Martinis, G. Finocchio, and K. Y. Camsari (2022) Massively parallel probabilistic computing with sparse Ising machines. Nature Electronics 5 (7), p. 460–468. External Links: Document Cited by: §1. [2] M. Aifer, Z. Belateche, S. Bramhavar, K. Y. Camsari, P. J. Coles, G. Crooks, D. J. Durian, A. J. Liu, A. Marchenkova, A. J. Martinez, P. L. McMahon, F. Sbahi, B. Weiner, and L. G. Wright (2025) Solving the compute crisis with physics-based ASICs. External Links: 2507.10463, Link Cited by: §1. [3] M. Aifer, K. Donatella, M. H. Gordon, S. Duffield, T. Ahle, D. Simpson, G. Crooks, and P. J. Coles (2024) Thermodynamic linear algebra. npj Unconventional Computing 1 (1), p. 13. External Links: Document Cited by: §1, §2.1, §2.1, §2.2.1, §4.1. [4] M. Aifer, S. Duffield, K. Donatella, D. Melanson, P. Klett, Z. Belateche, G. Crooks, A. J. Martinez, and P. J. Coles (2024) Thermodynamic Bayesian inference. In 2024 IEEE International Conference on Rebooting Computing (ICRC), p. 1–20. External Links: Document Cited by: §1. [5] A. Alaghi and J. P. Hayes (2013) Survey of stochastic computing. ACM Transactions on Embedded Computing Systems 12 (2s), p. 92:1–92:19. External Links: Document Cited by: §1, §3.1.1. [6] E. Alpaydin and C. Kaynak (1998) Optical recognition of handwritten digits. Note: UCI Machine Learning Repository External Links: Document Cited by: §4.2.1. [7] C. H. Bennett (1976) Efficient estimation of free energy differences from Monte Carlo data. Journal of Computational Physics 22 (2), p. 245–268. External Links: Document Cited by: Appendix E, §5.3. [8] S. Brooks, A. Gelman, G. L. Jones, and X. Meng (Eds.) (2011) Handbook of Markov chain Monte Carlo. Handbooks of Modern Statistical Methods, Chapman & Hall / CRC. External Links: ISBN 978-1-4200-7941-8 Cited by: §2.2.1, §2.3.2. [9] B. D. Brown and H. C. Card (2001) Stochastic neural computation I: computational elements. IEEE Transactions on Computers 50 (9), p. 891–905. External Links: Document Cited by: §3.1.2. [10] K. Y. Camsari, R. Faria, B. M. Sutton, and S. Datta (2017) Stochastic p-bits for invertible logic. Physical Review X 7 (3), p. 031014. External Links: Document Cited by: §1. [11] D. M. Chapiro (1984) Globally-asynchronous locally-synchronous systems. Ph.D. Thesis, Stanford University. Note: Technical Report STAN-CS-84-1026 Cited by: §3.1.3. [12] R. T. Q. Chen and Y. Lipman (2024) Flow matching on general geometries. In International Conference on Learning Representations (ICLR), Cited by: Appendix E, §5.3. [13] R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. Duvenaud (2018) Neural ordinary differential equations. In Advances in Neural Information Processing Systems (NeurIPS), External Links: 1806.07366 Cited by: §5.1, §5.3. [14] P. J. Coles, C. Szczepanski, D. Melanson, K. Donatella, A. J. Martinez, and F. Sbahi (2023) Thermodynamic AI and the fluctuation frontier. In 2023 IEEE International Conference on Rebooting Computing (ICRC), p. 1–10. External Links: Document Cited by: §1, §1. [15] T. Conte, E. DeBenedictis, N. Ganesh, T. Hylton, J. P. Strachan, R. S. Williams, A. Alemi, L. Altenberg, G. Crooks, J. Crutchfield, L. del Rio, J. Deutsch, M. DeWeese, K. Douglas, M. Esposito, M. Frank, R. Fry, P. Harsha, M. Hill, C. Kello, J. Krichmar, S. Kumar, S. Liu, S. Lloyd, M. Marsili, I. Nemenman, A. Nugent, N. Packard, D. Randall, P. Sadowski, N. Santhanam, R. Shaw, A. Stieg, E. Stopnitzky, C. Teuscher, C. Watkins, D. Wolpert, J. Yang, and Y. Yufik (2019) Thermodynamic computing. External Links: 1911.01968, Link Cited by: §1. [16] G. E. Crooks (1999) Entropy production fluctuation theorem and the nonequilibrium work relation for free energy differences. Physical Review E 60 (3), p. 2721–2726. External Links: Document Cited by: Discussion. [17] G. E. Crooks (2007) Measuring thermodynamic length. Physical Review Letters 99 (10), p. 100602. External Links: Document Cited by: Discussion. [18] K. Donatella, S. Duffield, M. Aifer, D. Melanson, G. Crooks, and P. J. Coles (2026) Thermodynamic natural gradient descent. npj Unconventional Computing 3 (1), p. 5. External Links: Document Cited by: §1. [19] S. Duffield, M. Aifer, D. Melanson, Z. Belateche, and P. J. Coles (2025) Lattice random walk discretisations of stochastic differential equations. External Links: 2508.20883 Cited by: §4.1. [20] S. Duffield (2026) A complete decomposition of stochastic differential equations. External Links: 2601.07834 Cited by: §5.1. [21] W. E (2017) A proposal on machine learning via dynamical systems. Communications in Mathematics and Statistics 5 (1), p. 1–11. External Links: Document Cited by: §5.1. [22] B. R. Gaines (1969) Stochastic computing systems. In Advances in Information Systems Science, Vol. 2, J. T. Tou (Ed.), p. 37–172. Cited by: §1, §3.1.1. [23] W. Grathwohl, R. T. Q. Chen, J. Bettencourt, I. Sutskever, and D. Duvenaud (2019) FFJORD: free-form continuous dynamics for scalable reversible generative models. In International Conference on Learning Representations (ICLR), Cited by: §5.3. [24] K. He, X. Zhang, S. Ren, and J. Sun (2016) Deep residual learning for image recognition. In IEEE Conference on Computer Vision and Pattern Recognition (CVPR), External Links: 1512.03385 Cited by: §5.1, §5. [25] J. Ho, A. Jain, and P. Abbeel (2020) Denoising diffusion probabilistic models. In Advances in Neural Information Processing Systems (NeurIPS), Cited by: §5.1. [26] L. Holdijk, N. M. Anand, M. M. Bronstein, and M. Welling (2026) Learning escorted protocols for multistate free-energy estimation. In The Fourteenth International Conference on Learning Representations, External Links: Link Cited by: Appendix E, Appendix E, §5.3.1, §5.3, Discussion. [27] S. Hooker (2021) The hardware lottery. Communications of the ACM 64 (12), p. 58–65. External Links: Document Cited by: §5. [28] M. F. Hutchinson (1989) A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines. Communications in Statistics - Simulation and Computation 18 (3), p. 1059–1076. External Links: Document Cited by: Appendix E. [29] T. Hylton (2020) Thermodynamic computing: an intellectual and technological frontier. Proceedings 47 (1), p. 23. Note: MORCOM 2019 External Links: Document Cited by: §1. [30] C. Jarzynski (2002) Targeted free energy perturbation. Physical Review E 65 (4), p. 046122. External Links: Document Cited by: Appendix E, §5.3. [31] C. Jarzynski (1997) Nonequilibrium equality for free energy differences. Physical Review Letters 78 (14), p. 2690–2693. External Links: Document Cited by: §5.3, Discussion. [32] A. Jelinčič, O. Lockwood, A. Garlapati, P. Schillinger, I. Chuang, G. Verdon, and T. McCourt (2025) An efficient probabilistic hardware architecture for diffusion-like models. External Links: 2510.23972, Link Cited by: §1. [33] J. Kaiser and S. Datta (2021) Probabilistic computing with p-bits. Applied Physics Letters 119 (15), p. 150503. External Links: Document Cited by: §1. [34] K. Kim, J. Kim, J. Yu, J. Seo, J. Lee, and K. Choi (2016) Dynamic energy-accuracy trade-off using stochastic computing in deep neural networks. In Design Automation Conference (DAC), External Links: Document Cited by: §3.1.1. [35] D. P. Kingma and M. Welling (2014) Auto-encoding variational bayes. In International Conference on Learning Representations (ICLR), External Links: 1312.6114 Cited by: §4.2. [36] P. E. Kloeden and E. Platen (1992) Numerical solution of stochastic differential equations. Applications of Mathematics, Vol. 23, Springer. External Links: ISBN 3-540-54062-8 Cited by: §5.1. [37] A. Krizhevsky (2009) Learning multiple layers of features from tiny images. Technical report University of Toronto. Cited by: §5.2. [38] Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner (1998) Gradient-based learning applied to document recognition. Proceedings of the IEEE 86 (11), p. 2278–2324. External Links: Document Cited by: §4.2.1. [39] D. A. Levin, Y. Peres, and E. L. Wilmer (2017) Markov chains and mixing times. 2nd edition, American Mathematical Society. External Links: ISBN 978-1-4704-2962-1 Cited by: Appendix A, Discussion. [40] J. Li, Z. Yuan, Z. Li, C. Ding, A. Ren, Q. Qiu, J. Draper, and Y. Wang (2017) Hardware-driven nonlinear activation for stochastic computing based deep convolutional neural networks. In International Joint Conference on Neural Networks (IJCNN), p. 1230–1236. Cited by: §3.1.2. [41] J. Lions, Y. Maday, and G. Turinici (2001) A “parareal” in time discretization of PDE’s. Comptes Rendus de l’Académie des Sciences - Series I - Mathematics 332 (7), p. 661–668. External Links: Document Cited by: §2.3.3, §3.2. [42] Y. Lipman, R. T. Q. Chen, H. Ben-Hamu, M. Nickel, and M. Le (2023) Flow matching for generative modeling. In International Conference on Learning Representations (ICLR), Cited by: Appendix D, Appendix E, §5.2, §5.3. [43] S. Liu, W. J. Gross, and J. Han (2020) Introduction to dynamic stochastic computing. IEEE Circuits and Systems Magazine 20 (3), p. 19–33. External Links: Document Cited by: §3.1.1. [44] S. Liu and J. Han (2017) Hardware ODE solvers using stochastic circuits. In Design Automation Conference (DAC), External Links: Document Cited by: §3.1.1. [45] X. Liu, C. Gong, and Q. Liu (2023) Flow straight and fast: learning to generate and transfer data with rectified flow. In International Conference on Learning Representations (ICLR), Cited by: Appendix D, §5.2.1. [46] D. Melanson, M. A. Khater, M. Aifer, K. Donatella, M. H. Gordon, T. Ahle, G. E. Crooks, A. J. Martinez, F. Sbahi, and P. J. Coles (2025) Thermodynamic computing system for AI applications. Nature Communications 16 (1), p. 3757. External Links: Document Cited by: §1, §2.1, §2.1, §2.2, §4.1, Discussion. [47] Z. Mensch, L. Holdijk, S. Duffield, M. Aifer, P. J. Coles, M. Welling, and M. C. N. Cheng (2026) Robust stochastic gradient posterior sampling with lattice based discretisation. External Links: 2602.15925 Cited by: §4.1. [48] S. P. Meyn and R. L. Tweedie (2009) Markov chains and stochastic stability. 2nd edition, Cambridge University Press. External Links: ISBN 978-0-521-73182-9 Cited by: Appendix A, §2.2, §2.3.1, §2.3.2. [49] F. Neugebauer, I. Polian, and J. P. Hayes (2018) S-box-based random number generation for stochastic computing. Microprocessors and Microsystems 61, p. 316–326. External Links: Document Cited by: §3.1.3. [50] F. Noé, S. Olsson, J. Köhler, and H. Wu (2019) Boltzmann generators: sampling equilibrium states of many-body systems with deep learning. Science 365 (6457), p. eaaw1147. External Links: Document Cited by: Appendix E, §5.3. [51] A. Ren, Z. Li, C. Ding, Q. Qiu, Y. Wang, J. Li, X. Qian, and B. Yuan (2017) SC-DCNN: highly-scalable deep convolutional neural network using stochastic computing. In International Conference on Architectural Support for Programming Languages and Operating Systems (ASPLOS), p. 405–418. External Links: Document Cited by: §3.1.1. [52] W. Romaszkan, T. Li, T. Melton, S. Pamarti, and P. Gupta (2020) ACOUSTIC: accelerating convolutional neural networks through or-unipolar skipped stochastic computing. In Design, Automation & Test in Europe Conference & Exhibition (DATE), p. 768–773. External Links: Document Cited by: §3.1.1. [53] O. Ronneberger, P. Fischer, and T. Brox (2015) U-Net: convolutional networks for biomedical image segmentation. In Medical Image Computing and Computer-Assisted Intervention (MICCAI), p. 234–241. External Links: Document Cited by: §5.2.1. [54] U. Seifert (2012) Stochastic thermodynamics, fluctuation theorems and molecular machines. Reports on Progress in Physics 75 (12), p. 126001. External Links: Document Cited by: Discussion. [55] A. Shih, S. Belkhale, S. Ermon, D. Sadigh, and N. Anari (2023) Parallel sampling of diffusion models. In Advances in Neural Information Processing Systems (NeurIPS), Cited by: §2.3.3, §3.2. [56] H. Sim, D. Nguyen, J. Lee, and K. Choi (2017) Scalable stochastic-computing accelerator for convolutional neural networks. In Asia and South Pacific Design Automation Conference (ASP-DAC), p. 696–701. External Links: Document Cited by: §3.1.1. [57] D. A. Sivak and G. E. Crooks (2012) Thermodynamic metrics and optimal paths. Physical Review Letters 108 (19), p. 190602. External Links: Document Cited by: Discussion. [58] J. Sohl-Dickstein, E. A. Weiss, N. Maheswaranathan, and S. Ganguli (2015) Deep unsupervised learning using nonequilibrium thermodynamics. In International Conference on Machine Learning (ICML), Cited by: §5.1. [59] K. Sohn, H. Lee, and X. Yan (2015) Learning structured output representation using deep conditional generative models. In Advances in Neural Information Processing Systems (NeurIPS), p. 3483–3491. Cited by: Appendix C, §4.2. [60] Y. Song, J. Sohl-Dickstein, D. P. Kingma, A. Kumar, S. Ermon, and B. Poole (2021) Score-based generative modeling through stochastic differential equations. In International Conference on Learning Representations (ICLR), External Links: 2011.13456 Cited by: §5.1. [61] G. M. Torrie and J. P. Valleau (1977) Nonphysical sampling distributions in Monte Carlo free-energy estimation: umbrella sampling. Journal of Computational Physics 23 (2), p. 187–199. External Links: Document Cited by: §5.3.1. [62] G. E. Uhlenbeck and L. S. Ornstein (1930) On the theory of the Brownian motion. Physical Review 36 (5), p. 823–841. External Links: Document Cited by: §2.1. [63] A. W. van der Vaart (1998) Asymptotic statistics. Cambridge Series in Statistical and Probabilistic Mathematics, Cambridge University Press. External Links: Document Cited by: Appendix A, §2.3.3. [64] S. Whitelam and C. Casert (2026) Nonlinear thermodynamic computing out of equilibrium. Nature Communications 17 (1), p. 1189. External Links: Document Cited by: §1. [65] S. Whitelam (2026) Generative thermodynamic computing. Physical Review Letters 136 (3), p. 037101. External Links: Document Cited by: §1. [66] B. Wu, A. Wan, X. Yue, P. Jin, S. Zhao, N. Golmant, A. Gholaminejad, J. Gonzalez, and K. Keutzer (2018) Shift: a zero FLOP, zero parameter alternative to spatial convolutions. In IEEE Conference on Computer Vision and Pattern Recognition (CVPR), Cited by: §5.2.1. [67] S. Zilberstein (1996) Using anytime algorithms in intelligent systems. AI Magazine 17 (3), p. 73–83. External Links: Document Cited by: §2.3.1. [68] R. W. Zwanzig (1954) High-temperature equation of state by a perturbation method. I. nonpolar gases. The Journal of Chemical Physics 22 (8), p. 1420–1426. External Links: Document Cited by: §5.3.2. Appendix Appendix A Formal statement of the equilibration formalism This appendix provides formal statements of the claims made in Section˜2. We restate the formalism under formal assumptions and show that consistency is preserved under both parallel aggregation and sequential composition. Formalism and ergodicity. Let S be a measurable state space, and let (St)t≥0(S_t)_t≥ 0 be a discrete-time Markov chain on S with transition kernel PxP_x that depends parametrically on the input x to the computation. The dynamical generator ℒL is the forward operator associated with PxP_x, acting on distributions by ℒμ=μPx−μLμ=μ P_x-μ. Suppose PxP_x admits a unique stationary distribution πx _x satisfying ℒπx=0L _x=0, equivalently πxPx=πx _xP_x= _x, and suppose the chain is ergodic with respect to πx _x. Let φ:S→ℝm :S ^m be a measurable observable with π[‖φ(S)‖]<∞E_π[\| (S)\|]<∞. The computation is f(x)≜π[φ(S)],f(x)\; \;E_π[ (S)], and the chip’s output for a run of length T is y^T=1T∑t=1Tφ(St). y_T\;=\; 1T _t=1^T (S_t). By the ergodic theorem [48], y^T→f(x) y_T→ f(x) almost surely as T→∞T→∞, regardless of the distribution of S0S_0. As is standard for Markov chain averages, the variance of y^T y_T falls as τcorr/T _corr/T with the integrated correlation time τcorr _corr, while before stationarity the estimator carries a burn-in bias that decays once T exceeds the relaxation time τrelax _relax [39]. Parallel aggregation. Consider K independent runs of the same chain, each yielding a per-run estimator y^T(k) y_T^(k) as above. The aggregate estimator is y¯T(K)=1K∑k=1Ky^T(k), y_T^(K)\;=\; 1K _k=1^K y_T^(k), and independence of the runs gives [y¯T(K)]=[y^T(1)]E[ y_T^(K)]=E[ y_T^(1)] and Var(y¯T(K))=Var(y^T(1))/KVar( y_T^(K))=Var( y_T^(1))/K, so the rate of variance reduction is set by the aggregate sample throughput KTKT and the runs need no coordination beyond the final average. Sequential parallelism. The framework permits the dynamics to be split across subsystems that share no clock and exchange only aggregated outputs. Consider an L-stage pipeline in which stage ℓ implements a thermodynamic computation fℓ(xℓ)=πℓ[φℓ(S(ℓ))∣xℓ]f_ (x_ )=E_ _ [ _ (S^( )) x_ ], with xℓ+1x_ +1 produced from the output of stage ℓ . Let y^T(ℓ) y_T^( ) denote stage ℓ ’s time-averaged readout, run on the upstream estimate y^T(ℓ−1) y_T^( -1), and assume each stage is consistent in the sense that y^T(ℓ)→fℓ(xℓ) y_T^( )→ f_ (x_ ) almost surely whenever its input converges to xℓx_ . If each fℓf_ is continuous in its argument, applying this inductively along the pipeline gives y^T(L)⟶(fL∘⋯∘f1)(x1)a.s. y_T^(L)\; \;(f_L ·s f_1)(x_1) .s. as the per-stage windows grow, by the extended continuous-mapping theorem [63]. Consistency is therefore propagated through downstream stages, and the argument requires only stage-level consistency and continuity of the per-stage maps, not a shared clock between stages. Appendix B Compiling a neural network on CN101 Every workload is run on CN101 by the same recipe. The network is trained off-chip in floating point, without custom layers, and its weights are quantised only at deployment, to 88-bit integers in 64×6464× 64 blocks loaded into the stochastic matrix multipliers. A layer that fits a single tile occupies one block; a wider or deeper layer is partitioned into 64×6464× 64 blocks and assigned across tiles and chips, the partition being the only model-dependent choice. When a model fits the four tiles of a single chip, as the VAE decoder does, each neuron bank applies its layer’s bias and ReLU on-chip and the cascade runs as one dispatch; when a model is spread across tiles and chips, the chip evaluates the 64×6464× 64 matmuls and the host applies the bias, normalisation, activation and any integration step between windows. B.1. Composition and multi-chip execution Every concurrent configuration in the paper exercises the sequential parallelism of Section˜2.3: each stage advances on the still-converging estimate of the stage feeding it rather than waiting for it to settle. Within a chip this exchange is continuous. The four tiles run on independent ring-oscillator clocks, with cross-tile signals synchronised at bank boundaries, and exchange bitstreams over the stochastic streaming network. A four-layer cascade of tiles therefore runs as a single-chip dispatch with one final readout and no host intervention between layers. Across chips the same composition is realised window by window. Over a window of Δ cycles each chip advances on its current input; the host then reads each chip’s running estimate from its output accumulator, applies any host-side bias, normalisation or integration step, and refreshes the downstream inputs before the next window. The sync interval Δ can be chosen between fully sequential execution and the fully concurrent regime Δ=1 =1. Integrated multi-chip hardware, on which the estimate exchange is routed on-board rather than passed through the host, is left to a subsequent chip in the Carnot Architecture. Appendix C Conditional VAE: training and deployment The decoder architecture, dataset and tile mapping are given in Section˜4.2; here we record the training and deployment detail. The encoder, used only during training, mirrors the decoder and outputs the mean and log-variance of the 1616-dimensional latent. The network is trained end to end in floating point with the reparametrisation trick and a class-conditional evidence lower bound [59], the reconstruction and Kullback–Leibler terms weighted equally (β=1β=1). Training uses the Adam optimiser at learning rate 10−310^-3 and batch size 128128 for 100100 epochs; no CN101-specific regularisation is applied, and the weights are quantised to 88-bit integers only at deployment. For inference the weights are held stationary for the duration of a run. The per-T operating points of Fig.˜6c,d are selected by a search over the chip’s two configurable knobs that do not change the trained model: the input scale applied to each tile and the weight-register precision. For each cycle budget T the configuration minimising the scale-invariant RMSE against the reference is chosen, so that the reported curve traces the best operating point available at each T rather than a single fixed setting. Appendix D Convolutional flow on CIFAR-10: training and schedule The ChipUNetTiny architecture, the shift-then-1×11× 1 convolution mapping and the tile pinning are given in Section˜5.2; here we record the training and readout schedule. The velocity field is trained by flow matching [42] on ten CIFAR-10 images with the Adam optimiser at learning rate 2×10−42× 10^-4 and full-batch gradients for 20,00020,000 steps, and is then reflowed once [45]: samples drawn from the trained model are paired with their originating noise and the field is retrained on the straightened couplings for a further 20,00020,000 steps, so that generation needs only N=10N=10 integration steps. The model is deliberately overfit to these ten images, the prototype’s on-chip weight memory being too small to hold a general CIFAR generator, so the experiment establishes that the substrate runs the deep convolutional flow correctly rather than that this network generalises. At inference every per-step matmul is read out as a time-average over a schedule rising to 5×1055× 10^5 cycles per tile, with the first windows, over which the unrolled relaxation fills, discarded as a burn-in before averaging. Appendix E Alanine dipeptide: model, data, and free-energy protocol Reference data. We model alanine dipeptide in implicit solvent and describe each configuration in internal bond–angle–torsion coordinates [50]: 2121 bonds and 2020 angles on ℝ41R^41, and 1919 backbone and side-chain dihedrals on the torus 19T^19. The six metastable basins (αR _R, αL _L, αD _D, β, C5 and α′α ) are defined by flat-bottom boxes in the backbone-dihedral plane (ϕ,ψ)(φ,ψ) (Fig.˜9), following the state definitions of Holdijk et al. [26]. For each basin we draw equilibrium configurations from a restrained molecular-dynamics trajectory (2020 ns per basin, a flat-bottom restraint of strength k=100k=100 kJ mol-1 rad-2 confining the dynamics to the basin). These restrained ensembles are the reference distributions the flow transports between, and the unrestrained potential U(x)U(x) supplies the energies entering the work. Flow-matching model. For each pair (αR,X)( _R,X) we train a Riemannian flow-matching model [42, 12] on the product geometry 19×ℝ41T^19×R^41, with the dihedral block carried through a cos/sin / embedding so that the velocity field respects the periodicity of the torus. The velocity field is a residual multilayer perceptron of two blocks of width 128128 with layer normalisation and ReLU activations, conditioned on the integration time. It is trained to transport the restrained αR _R ensemble onto the restrained X ensemble, with the reverse map the same flow integrated backwards, using the Adam optimiser at learning rate 10−310^-3 and batch size 256256 for 200200 epochs. Free-energy estimation. We estimate free-energy differences by targeted free-energy perturbation [30], with the trained flow as the targeting map T and the work of a transport given by Eq.˜14 at temperature T=300T=300 K, its log-determinant estimated with a single Hutchinson probe per integration step [28] (rank one). We train a star of five maps, one between the reference basin αR _R and each of the other basins, following the construction of Holdijk et al. [26]. For each pair (αR,X)( _R,X) we estimate the free-energy difference ΔF(αR,X) F( _R,X) with the Bennett acceptance ratio [7], which combines the forward works of the transport αR→X _R\!→\!X with the reverse works of X→αRX\!→\! _R. With αR _R as the common reference this fixes each free energy F^i F_i relative to it, and the remaining pairwise differences follow by subtraction, ΔFij=F^j−F^i F_ij= F_j- F_i. Each free-energy difference draws 6464 configurations from each basin, and uncertainties are bootstrap standard deviations over those samples. Chip deployment. The flow is deployed on the CN101 cascade with its weights quantised in 64×6464× 64 blocks. Each transport integrates 3232 Euler steps, with the Hutchinson tangent carried in the same forward sweep rather than a separate pass, and every per-step matmul is read out as a time-average over a twenty-window schedule whose per-window length rises geometrically to 10710^7 cycles per tile (Fig.˜10b). Appendix F Reference computation, metrics and reproducibility Each chip result is compared against a reference that runs the same 88-bit-quantised model in exact floating-point arithmetic on the host. Because reference and chip share the model and differ only in the substrate, the comparison isolates the contribution of the substrate from that of the model. Errors are reported per workload in the natural metric. For the VAE the output is scored by a scale-invariant RMSE against the reference, computed after a single global scale and offset are removed, so that the metric reflects structural rather than overall-gain error. For the convolutional flow the metric is the pixel-RMSE of the generated image against the reference. For alanine dipeptide it is the free-energy difference recovered through the Bennett acceptance ratio. The randomness each run requires is drawn from the on-chip generators with seeds verified to produce non-overlapping sequences across neurons and across runs, and the independent replicas combined in parallel sample aggregation use distinct such seeds. Error bars are stated with each result: the VAE figures average over ten digit classes and ten latent draws per class, and the free-energy estimates use bootstrap standard deviations over the drawn configurations.