Hidden Markov models for stochastic thermodynamics

The formalism of state estimation and hidden Markov models (HMMs) can simplify and clarify the discussion of stochastic thermodynamics in the presence of feedback and measurement errors. After reviewing the basic formalism, we use it to shed light on a recent discussion of phase transitions in the optimized response of an information engine, for which measurement noise serves as a control parameter. The HMM formalism also shows that the value of additional information shows a maximum at intermediate signal-to-noise ratios. Finally, we discuss how systems open to information flow can apparently violate causality; the HMM formalism can quantify the performance gains due to such violations.


Introduction
In 1867, at the dawn of statistical physics, Maxwell imagined a thought experiment that has both troubled and inspired physicists ever since [1]. In modern language, the issue is that traditional thermodynamics posits a strict separation between observable macroscopic motion (dynamical systems) and unobservable degrees of freedom (heat). But imagine-as can now be done experimentally on small systems where fluctuations are important-that it is possible to observe some of these hidden degrees of freedom. (Maxwell's thought experiment used a "demon" to accomplish the same task.) In any case, the entropy of the system is reduced, and one can use the lower entropy to extract work from the surrounding heat bath, in seeming violation of the Second Law of thermodynamics.
This blurring of macroscopic and microscopic degrees of freedom has led to a new field, stochastic thermodynamics, which clarifies how thermodynamics should be applied to small systems where fluctuations are observable and important [2]. As we will see below, the nature of information acquired about the fluctuations-especially the precision with which they are measured and the time they become available-is of great importance. Indeed, information is itself a thermodynamic resource, and stochastic thermodynamics can be extended to accommodate the acquisition, dissipation, flow, and feedback of information [3,4,5,6,7,8,9,10,11,12,13,14,15]. For a recent review, see [16].
The goal of the present contribution is to combine ideas from control theory (state estimation) [17,18] with ideas from computer science about hidden Markov models [19,20,21,22,23] in order to explain some recent surprising observations from stochastic thermodynamics about how Maxwell's demon operates in the presence of measurement errors [24]. As a bonus, the formalism we discuss suggests a number of interesting areas where the stochastic thermodynamics of information may be extended.

Coarse graining and discrete state spaces
In the simplest non-trivial example of a discrete state space, a state x can, at each discrete time point, take on one of two values, for example −1 and +1. While systems such as spin- 1 2 particles are inherently discrete, a broad range of physical systemseven classical, continuous state spaces-can often be well approximated by discrete systems after coarse graining. Figure 1(a) sketches such a system, a protein in solution that alternates between a loose unfolded (−1) and a compact folded (+1) state. Other biological examples of two-state systems include ion channels that can be open or closed, gene-transcription repressor sites that can be occupied or empty, and sensory receptors that can be active or silent (chapter 7 in [25]).  Figure 1 illustrates schematically how to coarse grain from a physical situation, such as a protein in water, to a discrete-time Markov model. In (a), we depict two states of the protein, labeled "unfolded" and "folded" or, equivalently, −1 and +1. The word "state" is here a shorthand for "macrostate" and is associated with many microstates, each of which corresponds to a slightly different protein conformation that preserves the general property in question. In (b), we project the full dynamics onto a one-dimensional subspace modeled by a double-well potential. States with x < 0 are classified as −1, and states with x > 0 are classified as +1. The symmetry of the potential implies that the protein spends equal time in the two states, which is a special situation. In (c), we show a graphical depiction of the discrete, two-state Markov chain dynamics, where in a time τ , states remain the same with probability 1 − a and hop to the other with probability a. In order for a two-state description to reasonably approximate the dynamics, the dwell time spent in each well must be much longer than the time scale for fast motion within a well. This holds when a single energy barrier E b separates two states and whose height is much larger than kT .
Why might we want to approximate physical systems by discrete state spaces?
• Clarity: We can isolate just the important degrees of freedom, letting the others be uncontrolled and even unobserved.
• Simplicity: The mathematical description is more straightforward.
• Generality: Any dynamics that can be modeled on a computer is necessarily discretized in both time and state.

Markov chains
Let us briefly recall the basics of discrete-state-space systems in discrete time. Consider a system described at time k by a state x k that can be in one of n possible states, indexed by the values 1 to n. The index is distinguished from its value, which, for a two-state system, might be {±1}, {0, 1}, or even {left, right}. Let P (x k = i) be the probability that, at time k, the system is in the state indexed by i. The distribution is normalized by enforcing n i=1 P (x k = i) = 1 or, more succinctly, x k P (x k ) = 1. For dynamics, we consider Markov chains, which are systems with discrete time and discrete states. The Markov property implies that the next state depends only on the current state, as illustrated graphically in figure 2, which may be compared with figure 1(c). states: x k-1 x k x k+1 For Markov chains, the dynamics are specified in terms of an n×n transition matrix A whose elements A ij ≡ P (x k+1 = i | x k = j) satisfy 0 ≤ A ij ≤ 1. That is, A ij gives the rate of j → i transitions. For example, a general two-state system has Notice that the columns of A sum to 1, as required by the normalization of probability distributions. In words, if you start in state j then you must end up in one of the n possible states, indexed by i. Figure 1(c) depicts (1) graphically, with a 0 = a 1 = a. A matrix with elements 0 ≤ A ij ≤ 1 and i A ij = 1 is a (left) stochastic matrix. Define the n-dimensional stochastic vector p k , whose elements p give the probability to be in state j at time k. Then 0 ≤ p (j) k ≤ 1 and j p (j) k = 1 and (2) More compactly, p k+1 = A p k , a linear difference equation with solution p k = A k p 0 known as the discrete-time master equation. Often, we seek the steady-state distribution, defined by p = Ap. One way to find p is to repeatedly iterate (2); another is to note that the steady-state distribution of probabilities corresponds to the eigenvector associated with an eigenvalue equal to 1. A stochastic matrix must have such an eigenvalue, since A − I is a matrix whose columns all sum to zero. They are then linearly dependent, with zero determinant.
For example, the two-state Markov model with transition matrix A given by (1) has eigenvalues λ = 1 and 1 − (a 0 + a 1 ). The normalized eigenvector corresponding to λ = 1 is For the symmetric case, a 0 = a 1 ≡ a and p * = 0.5 0.5 , independent of a. By symmetry, both states are a priori equally probable.

Hidden Markov models
Often, the states of a Markov chain are not directly observable; however, there may be measurements (or emitted symbols) that correlate with the underlying states. The combination is known as a hidden Markov model (HMM). The hidden states are also sometimes known as latent variables [22]. The observations are assumed to have no memory: what is measured depends only on the current state, and nothing else. The graphical structure of an HMM is illustrated in figure 4.
In the example of proteins that alternate between unfolded and folded states, the molecule itself is not directly observable. One way to observe the configuration is to attach a particle to one end of the protein and anchor the other end to a surface [26], as illustrated in figure 4(a). As the protein folds and unfolds, the particle moves up and down from the surface. We can illuminate the region near the surface using an evanescent wave via the technique known as total internal reflection microscopy. The intensity I(z) of light scattered by the bead at height z from the surface will decrease exponentially as I(z) ∝ e −z/z 0 , with z 0 ≈ 100 nm. The two states will then correspond to two different hidden: observed: x k-1 x k x k+1 y k-1 y k+1 y k Figure 3. HMM graphical structure. The states x k form a Markov process that is not directly observable. The observations y k depend only on x k . scattering intensities. The observation y k is the number of recorded photons, integrated over a time that is shorter than the dwell time in each local potential well.
Hidden Markov As with states, we can further simplify by discretizing the intensities, classifying as "dim" intensities below a given threshold and "bright" intensities above that threshold. "Dim" and "bright" then become two observation symbols. Because light scattering is itself a stochastic process, the protein can be in one state but emit the "wrong" symbol, as illustrated in figure 4(b). We can describe such a situation by defining the observations y k = ±1 and noting that they are related to the states probabilistically via an observation matrix B having components B ij ≡ P (y k = i|x k = j): where we suppose, for simplicity, that errors are symmetric. Because observations have no memory, the probability to observe y k depends only on the current state x k .
In words, the matrix B states that an observation is correct with probability 1 − b and wrong with probability b. Like the transition matrix A, the matrix B is stochastic, with columns that sum to 1. Its rows also sum to 1, but only because of the symmetry between states. Note that the number of observation symbols, m, need not equal the number of internal states, n. The m × n matrix B can have m bigger or smaller than n. The case of continuous observations (m → ∞) is also straightforward. Larger values of m increase knowledge of the underlying state somewhat.
One interesting feature of HMMs is that states x k follow a Markov process and so does the combined process for x k and y k , but not necessarily the observations y k . The analysis of HMMs is thus more difficult than for ordinary Markov processes.
The literature on HMMs is both vast and dispersed. For treatments of increasing complexity, see section 16.3 of Numerical Recipes [19], the bioinformatics book by Durbin et al. [20], a classic tutorial from the speech-recognition literature [21], the control-influenced book by Särkkä [27], and the mathematical treatment of Cappé et al. [23]. The tutorial by Rabiner [21] has been particularly influential; however, its notation and ways of deriving results are more complicated than need be, and some of its methods have been replaced by better algorithms. The discussion here is based largely on the cleaner derivations in [27].

State estimation
Hidden Markov models are specified by a transition matrix A and observation matrix B. Let us pose the following problem: Given the output of a hidden Markov model (HMM), what can be inferred about the states? The answer depends both on the information available and the exact quantity desired. Here, we focus on two cases: (i) Filtering, or P (x k |y k ). We estimate the probabilities for each state based on observations y k ≡ {y 1 , y 2 , . . . , y k } up to and including the present time k. Filtering is appropriate for real-time applications such as control. ‡ Another quantity of interest is the most likely path, defined as arg max x N P (x N |y N ), which may be found by an algorithm due to Viterbi [19]. For example, McKinney et al. study transitions between different configurations of a DNA Holliday junction, using fluorescence resonance energy transfer (FRET) to read out the states, and infer the most likely state sequence [28]. Since path estimates are less useful for feedback control, we will consider them only in passing, in section 8. We will also see that smoothing estimates provide a useful contrast with filter estimates.

Filtering
The filtering problem is to find the probability distribution of the state x k based on the past and current observations y k from time 1 to time k. We assume that the dynamics have been coarse grained to be Markov, so that the state x k+1 depends only on the state x k . Then P (x k+1 |x k , ✓ ✓ y k ) = P (x k+1 |x k ), where the "cancel" slash indicates conditional independence: conditioning on x k "blocks" the influence of all other variables. The x k−1 are blocked, too: the state at time k + 1 depends only on the state at time k.
From marginalization and the definition of conditional probability, we have Equation (5) predicts the state x k+1 on the basis of y k , assuming that the previous filter estimate, P (x k |y k ) is already known. Once the new observation y k+1 is available, we can use Bayes' Theorem and the memoryless property of observations, P (y k |x k , ✟ ✟ ✟ y k−1 ) = P (y k |x k ), to update the prediction (5) to incorporate the new observation. Then, where Z k+1 normalizes the distribution. Equations (5)-(6) constitute the Bayesian filtering equations [27,29]. Because of their importance, we collect them here: The normalization (partition function) Z k+1 is given by Note that the HMM literature, e.g., [19] and [20] , expresses (7) differently, using joint probabilities such as P (x k , y k ) rather than conditional probabilities such as P (x k |y k ). Using joint probabilities leads to the forward algorithm. Our notation emphasizes the similarities between HMM and state-space models of dynamics; the formulas of one apply mostly to the other, with x k ↔ dx k . For continuous state spaces with linear dynamics and Gaussian noise, (7) is equivalent to the Kalman filter [27]. Below, we will see that using conditional probabilities also has numerical advantages. Figure 5 shows filtering in action for a symmetric, two-state, two-symbol hidden Markov model. The time series of observations y k (markers) disagrees with the true state 30% of the time. The black line shows P (x k = 1|y k ). When that probability is below the dashed line at 0.5, the most likely state is 0. For the value of a used in the dynamic matrix (a = 0.2), the filter estimate x with the observation only 7% of the time, a noticeable improvement over the naive 30%. Notice that whenever the state changes, the filter probability responds, with a time constant set by both observational noise (b) and dynamics (a). A long string of identical observations causes filter confidence to saturate at p * (dashed line).
There is an advantage to recording the probability estimates (black line) rather than simply the MAP (maximum a posteriori) estimate, which here is just the more likely of the two possibilities. When the filter is wrong, the two probabilities are often not that different. An example is indicated by the arrow in figure 5. Thus, marginalizing (averaging) any prediction over all possibilities rather than just the most likely will improve estimates. Of course, a string of wrong symbols can fool the filter. See, in figure 5, the three wrong symbols just to the left of the arrow.
Below, we will see that the filtered estimate becomes significantly more reliable as a → 0. Intuitively, small a means that states have a long dwell time, so that averaging observations over times of the order of the dwell time can reduce the effect of the observational noise, which is quantified by the parameter b.

Smoothing
If we estimate the state x k after gathering N observations (N > k), we can use the "future" information to improve upon the filter estimate. In the control-theory literature, such estimates are called "smoother" estimates, as they further reduce the consequences of observation noise.
The smoother estimate has two stages. First, we use the filter algorithm (7) to calculate P (x k |y k ) and P (x k+1 |y k ) for each k ∈ [1, N]. Then we calculate P (x k |y N ) via a backward recursion relation from the final time N to the initial time 1.
The backwards recursion relation is initialized by P (x N , y N ), the last step of the forward filter recursion. To derive (9), we introduce the state x k+1 , which we will remove later by marginalization [27]. Thus, But using conditional probability and the Markov property. Substituting into (10), Summing both sides over x k+1 gives (9). The algorithm defined by (7) and (9) is equivalent to the Rauch-Tung-Striebel smoother from control theory when applied to continuous state spaces, linear dynamics, and white-noise inputs [27]. In the HMM literature, a close variant is the forwardbackward algorithm [20]. We can apply the smoother algorithm to the example of section 5.1 and obtain similar results. In figure 6, we plot the smoother estimate, with the filter estimate added as a light gray trace. Despite their similarity, the differences are instructive: The filter always lag (reacts) to observations, whereas the smoother curve is more symmetric in time. Flipping the direction of time alters the overall form of the filter plot but not the smoother. The smoother estimates are more confident than the filter estimates, as they use more information. Look at the time step indicated by the arrow. The filter estimate is just barely mistaken, but the smoother estimate makes the correct call, aided by the three correct observations that come before and the three after.
The phase lag apparent in the filter estimate is consistent with causality. Indeed, for continuous state spaces, the well-known Bode gain-phase relations-the "magnitudephase" equivalent of the Kramers-Kronig relations [30]-give the minimum phase lag for the output of a dynamical system that is consistent with causality. The smoother estimate in figure 6 has zero phase lag, as expected since it uses past and future information equally. Sudden jumps are anticipated by the smoother before they happen.
Intuitively, an estimator that uses more information should perform better. We can formalize this intuition via the notion of conditional Shannon entropy [31]. With where using a base-2 logarithm gives units of bits. For large-enough k, the average of H(x k |y k ) over y k becomes independent of k. Averaging over a single long time series of observations then leads to H(x k |y k ) = H(x| ← − y ), where ← − y denotes past and present observations. A similar definition holds for the smoother entropy, H(x k |y N ) and leads to a steady-state smoother entropy H(x| ← → y ), where ← → y includes both past and future observations. To characterize the performance of filtering and smoothing, we recall that for a two-state probability distribution, the entropy ranges from 1 bit (equal probabilities for each possibility) to 0 bits (certainty about each possibility).

Learning hidden Markov models
The state-estimation procedures described above assume that the transition matrix A, the emission matrix B, and initial probability P (x 1 ) are known. If not, they can be estimated from the observations y N . In the context of HMMs, the task is called, variously, parameter inference, learning, and training [19]. In the control-theory literature on continuous state spaces, it is known as system identification [32]. The general approach is to maximize the likelihood of the unknown quantities, grouped here into a single parameter vector θ. That is, we seek where it is better to compute L(θ) ≡ − ln P (y N |θ) because P (y N |θ) decreases exponentially with N, leading to numerical underflow. The negative sign is a convention from least-squares curve fitting, where χ 2 (θ) is also proportional to the negative log likelihood of the data [19]. We can find the total likelihood P (y N |θ) from the normalization condition in (7): where Z 1 ≡ P (y 1 ). Then where all right-hand-side terms depend also on θ. Since L(θ) is just a function of θ, we can use standard optimization routines to find the θ * that minimizes L.
In the HMM literature, an alternate approach to finding θ * is based on the Expectation Maximization (EM), or Baum-Welch algorithm [22,20]. In a two-step iteration, one finds θ by maximum likelihood assuming that the hidden states x N are known and then infers states x N from the smoother algorithm assuming θ is known. The algorithm converges locally but can be very slow. Indeed, the EM algorithm can seldom compete against the more sophisticated direct-optimization algorithms readily available in standard scientific programming languages. EM algorithms can, however, be the starting point for recursive variants that allow for adaptation [33]. A third approach to finding HMM parameters, based on finding the most likely (Viterbi) path, can also converge faster than EM and be more robust [34].

Control of discrete-state-space systems
We can now discuss the control of Markov models and HMMs. In the context of discrete state spaces, the control u k influences the transition probability, which becomes P (x k+1 |x k , u k ) and is described by a time-dependent transition matrix A k and a graphical structure illustrated in figure 8. Note that our previous discussion of state estimation (filtering) never assumed that the transition matrix is time independent. hidden: observed: x k-1 x k x k+1 y k-1 y k+1 y k The control of Markov chains is formally known as a Markov Decision Process (MDP), while that of HMMs is known as a Partially Observable Markov Decision Process (POMDP). Optimal-control protocols that minimize some cost function can be found using Bellman's dynamic programming, which is a general algorithm for problems involving sequential decisions [19,35]. In this setting, control is viewed as a blend of state estimation and decision theory [35,36]. The goal is to choose actions based on available information in order to minimize a cost function.
Here, we will present such ideas more informally, using a well-studied example: optimal work extraction from a two-state system with noisy observations and feedback. This problem is closely related to a famous thought experiment (recently realized experimentally [37]), Maxwell's demon.

Maxwell's demon
As discussed in the introduction, a Maxwell demon is a device where information about the state of a system is used to extract energy from a heat bath, in violation of the traditional form of the Second Law of thermodynamics. How is this possible? The catch is that we have assumed that information carries no cost. A first attempt at resolving the paradox hypothesized that energy is dissipated in acquiring information [1]. However, that turns out not to be true in general: one can sometimes acquire information without doing work. In its Kelvin-Planck formulation, the Second Law requires that no cyclic protocol of parameter variation can extract work from the heat bath of an equilibrium system held at constant temperature. Specifying a cyclic protocol can be subtle. Naively, a cyclic protocol requires that any potentials that are changed must be returned to their initial state; any mechanical part (pistons, etc.) that are moved must be moved back; and so on. But it also applies to information. In particular, any information acquired must be erased. In 1961, Landauer proposed that the erasure step necessarily required energy dissipation of at least kT ln 2 per bit, an amount that equals or exceeds the amount of work that can be extracted, thus saving (or extending) the Second Law [38]. Landauer's prediction has recently been confirmed experimentally [39,40], as has its converse, the Szilárd engine, which uses acquired information to extract work from a heat bath [37,41,42].  Figure 9. Converting information to work in a two-state system that hops back and forth between "left-well" and "right-well" states separated by a high energy barrier E b .
If the system is observed to be in its right-well state, then we can raise the left well without doing work. After a time τ , the well is lowered. If the left state is occupied, we extract an energy E that can be used to perform work.

A simple model, with fully observed states
We consider a particle in a fluid, subject to a double-well potential that may be manipulated by the experimenter ( figure 9). It is a useful setting for thinking about the issues raised by a Maxwell demon and is a situation that can now be realized experimentally [39,40]. We assume that the energy barrier is large (E b ≫ kT ), so that we can coarse grain to two-state Markov dynamics, as discussed in section 2. Henceforth, we set kT = 1. At intervals τ , we observe the state of the system and record which well the particle is in. For now, we assume this measurement is never wrong.
To extract work from a heat bath, we implement the following protocol: At t = 0, the potential is symmetric, with no energy-level difference between left and right wells. We then observe the particle. If we determine it to be in the right well, then, with no significant time delay, we quickly raise the left well to an energy E (and vice versa if in the left well). Raising the left well costs no work if we change the potential only where the particle is not present. From Sekimoto's formulation of stochastic energetics, the work done by an instantaneous change of potential is just ∆U, the change of potential evaluated at the position of the particle [43].
We then wait a time τ , keeping fixed the energy E of the left well. At some time, the particle may spontaneously hop to the left well, because of thermal fluctuations. At time τ , the left well is quickly lowered back to E = 0. If the particle happens to be in the left well, we extract an energy E from the heat bath. If not, no energy is extracted. Summarizing, the protocol is to measure the state; then raise the appropriate well by E and wait τ ; then lower the well back to 0.
Over many trials, the average extracted work W is given by E p τ , where p τ is the probability for the particle to be in the left well at time τ . But p τ also depends on E. To evaluate the relation, we consider the continuous time dynamics of the state of the system, allowing hops between states at arbitrary times t but still considering the hops themselves to be instantaneous. The discrete-time master equation p k+1 = A p k then becomesṗ = A p, where the matrix A has columns that sum to zero, to keep p normalized at all times. Normalization implies that a two-state system has but one independent evolution equation, p(t), which obeyṡ where ω − is the transition rate out from the left well and ω + is the transition rate into the left well. In equilibrium, detailed balance requires that ω + /ω − = e −E . Scaling time so that ω − = 1 then givesṗ = −p + e −E (1 − p) .
Settingṗ = 0 gives the steady-state solution p ∞ = 1/(e E + 1). Notice that E = 0 implies p ∞ = 1 2 , as expected for a symmetric double-well potential, and that E → ∞ implies that the particle is always in the right well (p ∞ → 0). For finite times, we solve (18) Note that we choose signs so that W > 0 corresponds to work extraction. Intuitively, for a given cycle time τ , an optimal energy E * maximizes the average work: if E is too small, you will extract work in many cycles, but the amount each time will be small. If E is too large, you will extract more work, but only very rarely, since the relative probability of being on the left side is e −E . For the quasistatic limit τ ≫ 1, W ≈ E/(e E + 1), whose maximum W * ≈ 0.28 for E * ≈ 1.28.
The second law of thermodynamics implies that W ≤ ∆F , where the free energy difference ∆F is just the difference in entropy ∆S, since the internal energy difference is zero for a cyclic process where the energies of both states are identical at beginning and end. The maximum entropy difference is ln 2 ≈ 0.69, which is considerably larger than the ≈ 0.28 found in the quasistatic limit of our protocol.
To achieve the ln 2 upper bound for extracted work per cycle, we need to allow E(t) to vary continuously in the interval 0 < t < τ (and to have jump discontinuities at the beginning and end of the interval). Such continuous-time protocols have been considered previously and lead to protocols that extract ln 2 of work in the quasistatic limit [44,45,24]. Nonetheless, we prefer our constant-E protocol: • The mathematics is simpler. The continuous version uses calculus of variations.
The discrete one requires only ordinary calculus.
• If implemented experimentally, the protocols would almost certainly be carried out digitally, with an output that is fixed between updates.
• When the goal is to optimize power extraction from the heat bath (rather than work per cycle), the constant-E and continuous protocols give identical results.
To explore this last point, we rewrite (19) for average power, P ≡ W /τ . Assuming, as a more careful analysis confirms, that maximum average power extraction occurs when τ ≪ 1, we have which has a maximum P * = 1/e ≈ 0.37 for E * = 1. The same result is found for the continuous protocol [24]. Since maximum energy extraction requires quasistatic, infinitely slow manipulations, the power at maximum energy tends to zero. Maximizing power extraction is arguably more interesting experimentally.

Hidden states
So far, we have assumed noise-free observations. If the observations are noisy, we have to infer the probability p(0) ≡ p 0 that the particle is in the left well. Assuming that the particle is likely in the right well (0 < p 0 < 1 2 ), then we should raise the left well. After a time τ has elapsed, (18) implies that with ω = e −E and ε ≡ e −(1+ω)τ . This expression is linear in p 0 , as the master equation (18) is linear. The discrete-time master equation for time step τ then is Matching terms with (21) gives β = ω(1−ε) 1+ω and α = ω+ε 1+ω . The complements are 1 − β = 1+ωε 1+ω and 1 − α = 1−ε 1+ω . Thus, when the left well is raised, the transition matrix A L (E, τ ) is Notice that the columns of A L sum to one, as they must and that the Markov transition matrix is no longer symmetric, as expected since we raise one of the wells. The novel aspect for us is that the transition matrix A L now depends on the energy level E, which can be set at each time step. When the right well is raised, matrix elements are switched, with left ↔ right. This amounts to swapping "across the diagonal" of the matrix. Thus, The previously analyzed case (19) for p 0 = 0 then represents the best-case scenario: the particle is definitely on the right, and there is never a penalty for raising the left well. For 0 < p 0 < 1 2 , we will occasionally do work in raising the well when the particle is present. Using (21) and maximizing over E, we can quickly calculate the maximum work extraction as a function of p 0 . Figure 10(a) shows that the maximum average extracted work decreases as the initial state becomes more uncertain. When p 0 = 1 2 , we have no information about the state of the system and cannot extract work from the heat bath, in accordance with the usual version of the second law. For p 0 > 1 2 , we would raise the right well, else we would be erasing information and heating the bath, rather than extracting energy from it. Figure 10(b) shows that the work extracted is nearly a linear function of the change in Shannon entropy between initial and final states. As in Szilárd's analysis, information was used to extract work from the heat bath. Here, the average slope (converted to nats) gives an efficiency of roughly 41%. Less than half the information gained is extracted as work by this particular protocol.

Two protocols
We have not yet specified how to estimate p 0 at the beginning of each time interval. We do so via the observations y k that are made at the beginning of each control period τ , before the choice of E. The observations have two symbols and are characterized by an observation matrix of the form of (4), with b the symbol error rate. We thus return to the formalism discussed in section 5, where p 0 → P (x k ), the state of the system at time k. Similarly, p τ → P (x k+1 ). The only difference is that we modify A by choosing E and which well to raise at each time step. Call the choice A k .
We can incorporate observations in two ways. One is to use only the observation y k to estimate P (x k ). Then Bayes' Theorem implies that P (x k |y k ) ∝ P (y k |x k ), where the prior P (x k ) = 1 2 , since left → right and right → left state transitions are equally likely. Although P (x k+1 |x k , u k ) does not satisfy this condition, the time-averaged sequence of transition matrices does: since left and right levels are raised at equal frequencies, the overall statistics are symmetric in the absence of other information. Here, u k is the control variable, a function of E.
The second way is to use the filtering formalism developed in section 5.1 to recursively compute P (x k |y k ). (Without information about the future, we cannot use smoothing.) We can say that the second strategy, which depends on past observations, uses memory whereas the first uses no memory. The procedure is then to • Measure y k .
• Update P (x k |y k ), based on {A k , B}, with the time-dependent transition matrix A k given by P (x k+1 |x k , u k ). The control u k is a function of E k .
• Determine E k+1 by minimizing W (E), the average work extracted in a cycle.
Iterated, the above algorithm leads to plots of the average extracted work as a function of the measurement-error probability b ( figure 11 ). In (a), the curve labeled memory, uses the Bayesian filter to estimate the state of the system. By "memory," we mean that the inference about which energy level to alter is based on all the observations y k up to time k. By contrast, in (b), the "no memory" curve uses only the current observation, y k . As before, the extra information from past states is most useful at intermediate values of error rate b. The difference curve, plotted at left below, resembles figure 7, which compared estimator entropies of the smoother and filter state estimates. The conclusion, again, is that extra information is most useful at intermediate signalto-noise ratios. Here, retaining a memory of past observations via the filter allows the Maxwell demon to extract more power from the heat bath.

Phase transition in a Maxwell demon
The continuous-protocol version of the Maxwell demon shows phase transitions in the behavior of the Maxwell demon as the symbol error rate b is varied [24]. To see that similar phenomena arise in the constant-E protocol discussed in this paper, compare the outcomes of the strategy that uses memory (y k ) with one using no memory (y k ). More precisely, we define a "discord" order parameter D, where y = ±1 represents the time series of observations andx = ±1 represents the state estimate, based in this case on the optimal filter. § If y andx always agree, D = 0. If y andx are uncorrelated, D = 1. Partial positive correlations imply 0 < D < 1. Put differently, D > 0 implies that there is value in having a memory, as the filter estimatex can differ from the observation. When D = 0, the filter always agrees with the observation, implying that there is no value in calculating the filter.
In figure 12(a), we plot the discord order parameter D against the symbol error rate b for three different cycle times, τ = 0.1, 1, and 10. There are many interesting features. For long cycle times, represented by τ = 10 and hollow markers, observations match the inferred state-defined here to be the more likely state, as determined by the probabilities from the filter algorithm. represented by τ = 0.1 and black markers, we observe two transitions that, upon closer inspection, are both discontinuous, corresponding to first-order phase transitions and marked by down-pointing black arrows. Finally, at b = 0.5, the order parameter D = 1, since there is no correlation between observation and the internal state (or its estimate). Interestingly, there is always a jump discontinuity in D at b = 0.5.

Phase transitions in state estimation
The phase transition observed in the Maxwell-demon model given in the previous section can also be seen in hidden Markov models that have nothing to do with thermodynamics. Figure 12(b) shows the discord order parameter D for a two-state, two-symbol HMM with x, y ∈ {−1, +1}, for three values of a. As in figure 12(a), there are first-order transitions for small values of a, continuous transitions for intermediate values, and no transitions for larger values. Intuitively, we need long dwell times in states (low values of a) so that we have time to average over (filter) the observation noise. If so, we may be confident in concluding the true state is different from the observed state. If the dwell time is short (high value of a), the best strategy is to trust the observations. Note that the values of a correspond roughly to the same regimes as implied by the values of τ ; however, we cannot make an exact mapping, since the Markov transition rate in the Maxwell-demon depends on the control u k , which depends on observation errors b.
As with the Maxwell-demon example, for given a there is a critical value of b, denoted b c . To calculate b c , we note that there is an upper limit to the confidence one can have in a given state estimate. As we can see in figure 5, this limit is achieved after a long string of identical observations, say y k = 1, that is {y 1 = 1, y 2 = 1, . . . , y k = 1}. See the string of eight +1 states in figure 5 as an example. More formally, we consider P (x k = 1|y k = 1). For k ≫ 1, the maximum value of the state probability approaches a fixed point p * at long times. The intuition is that even with a long string of +1 observations, you cannot be sure that there has not just been a transition and an accompanying observation error. We derive p * (a, b) in Appendix A and plot the results in figure 13(a). Let us denotex , the filter estimate of x k . To find conditions wherex (f) disagrees with y, we construct the extreme situation where a long string y k = 1 gives the greatest possible confidence that x k = 1. Then let y k+1 = −1. The discordant observation must lower the confidence in x k+1 to below 1 2 in order for the filter estimate and observation to disagree. Thus, the condition defining b c is Writing this condition out explicitly gives, after a calculation detailed in Appendix A, A similar calculation for the smoother, again detailed in Appendix A, leads to b smoother Figure 13(b) shows that the thresholds of simulated data agree with (27) and (28). Both filter and smoother estimates imply that there is a maximum value of a, call it a c , above which D = 0 for all b. For the filter a c = 1 4 , while for the smoother, a c = 1 3 . The higher value of a c reflects the greater value of smoother vs. filter inferences.

Mapping to Ising models
Although we have explained some features of figure 12, there is clearly more to understand. For example, there are both continuous and discontinuous transitions, as well as evidence for multiple transitions at fixed a. To begin to understand the reason for multiple phase transitions, we note that the two-state, two-symbol HMM can be mapped onto an Ising model [47,48]. Let us change variables: We use these definitions to formulate a "Hamiltonian" H = − ln P (x N , y N ) via where we have dropped constant terms that are independent of x k and y k . For a < 1 2 , the interaction term J > 0 is ferromagnetic: neighboring "spins" tend to align. The term h corresponds to an external field coupling constant. The field hy k is of constant strength and, for b < 1 2 , has a sign is equal to the observation y k . The picture is that a local, quenched field of strength hy k tries to align its local spin along the direction defined by y k . Notice that h = 0 for b = 1 2 : spins are independent of y k : observations and states decouple. A further change of variables (gauge transformation), z k = y k x k and τ k = y k y k+1 , gives which is a random-bond Ising model in a uniform external field h [49]. Starting in the late 1970s, both random-bond and random-field one-dimensional Ising chains were extensively studied as models of frustration in disordered systems such as spin glasses. In particular, Derrida et al. showed that the ground state at zero temperature has a countable infinity of transitions at h = 2J/m for m = 1, 2, . . . , ∞ [50]. Their transfer-matrix formalism is equivalent to the factorization of the partition function Z = k Z k given in (15).
The lowest-order transition, h = 2J, corresponds to a case where the external field at a site forces the local spin to align, because we are at zero temperature. In terms of the original HMM problem, the ground state corresponds to the most likely (Viterbi) path discussed briefly in section 5 [48]. While the Viterbi path differs from the filter estimate considered here, there may be a similar explanation for the multiple transitions apparent in figure 12.

Discussion
The formalism of hidden Markov models, or HMMs, can both simplify and clarify the discussion of stochastic thermodynamics of feedback using noisy measurements. Expressed in terms of the control-theory notation developed here, state estimation based on HMM formalism is an effective way to incorporate the effects of noisy measurements. As an application, we simplified a previous analysis of a Maxwell demon that uses observations to rectify thermal fluctuations. We saw that a surprising phase transition in the "discord" between observation and inferred state is also present in simple HMM models. At least in this case, the primary source of complexity seems to lie in the process of state estimation, rather than some feature of the thermodynamics.
Our study of phase transitions in the discord parameter follows the methods of Bauer et al. [24]; however, the mathematics is considerably more complicated in that case. We note that while Bauer et al. do observe a series of transitions in their numerics, they have not seen evidence for jump discontinuities (private communication). Perhaps the differences are also associated with the continuous protocol for varying E. More investigation is warranted.
Beyond simplifying specific calculations, the use of HMMs leads to other insights. For example, in figure 11, we saw that using a memory improves the performance of a Maxwell demon that extracts power from a heat bath. The greatest improvement was for intermediate values of the noise parameter b. Sivak and Thomson, studying a simple model of biological sensing, reached a similar conclusion [51].
The results presented here suggest a somewhat broader view. Figure 7 shows a similar result, where the smoother estimate outperforms the filter estimate. Here, performance is measured by the Shannon entropy of the estimated probability distribution. Again, we see that the best performance, relative to without memory, is at intermediate noise levels. Indeed, a variety of similar results can be obtained from many analogous quantities. For example, filter estimates based on continuous measurements with Gaussian noise also exceed those based on discrete observation measurements, with, again, a maximum at intermediate values of observation noise.
The common feature in all these different examples is that we compute some measure of performance-work extraction, Shannon entropy, etc.-as a function of added information. This added information can be previous observations ("memory"), offline observations, extra measurement precision, multiple measurements, and so on. In all cases, the greatest improvement is always at intermediate noise levels or, more precisely, at intermediate levels of signal-to-noise ratio. Intuitively, the observation makes sense: if information is perfect (zero noise), then more is superfluous. If information is worthless (zero signal), then more is again not better. But in intermediate cases, extra information adds value. Thus, Extra information is most useful at moderate signal-to-noise ratios.
It would be interesting to try to formalize these ideas further by defining a kind of "information susceptibility" in terms of a derivative of power extraction, etc. with respect to added information. In this context, it is worth noting the study by Rivoire and Leibler, who show that the value of information can be quantified by different information theoretic quantities, such as directed and mutual information, when the analysis is causal or acausal [52].
Finally, we note that while we have been careful to discuss the smoother as an offline analysis tool whereby data is analyzed after the fact, there are more interesting possibilities. As stochastic thermodynamics is generalized to accommodate information flows, we should also consider the equivalent to open systems. For quantities such as energy, we are used to the idea that a subsystem need not conserve energy and that we must account for both energy dissipation and energy pumping. Analogously, for information, we should consider both dissipation and the consequences of added information. Because such information comes from "outside" the system under direct study, causality need not be respected. For example, consider the problem of controlling the temperature of a house. A causal control system will simply respond to temperature perturbations after they occur. If it gets cold, the heater turns on. On the other hand, we know in advance that at night it gets cold, and we know, with effectively absolute certainty, the time the sun will set. Thus, we can anticipate the arrival of a cold perturbation and start to compensate for its effects before they occur. The resulting performance gain will be precisely analogous to the results shown in figure 6, where we compare filter and smoother estimates. (The quality of state estimates limits the quality of control.) The analysis of noisy discrete dynamics of HMMs is perhaps the simplest non-trivial setting where these ideas may be explored. More generally, outside influences will appear as additional inputs to a state node in a graphical representation. In this context, the Bayesian treatment of causality due to Pearl shows how to generalize inferences such as filtering and smoothing to Bayesian networks, which have a richer graphical structure than the chain-like Markov and HMMs sketched in figures 2, 4, and 8 [53,36]. Such techniques have been used in stochastic thermodynamics to study information thermodynamics on networks [54] and would seem to be the right approach to studying systems that are "causally open." In conclusion, we have introduced some of the properties of hidden Markov models that make them useful for simplifying the analysis of stochastic thermodynamics in the presence of feedback and noisy measurements, and we have seen how they suggest interesting areas for future research.
In the a-b parameter plane, the critical line b c (a) defines the border between the D = 0 and D > 0 phases. Informally, the line separates a region where there is no benefit to using the filter estimate from one where there is. We can use both filter and smoother state estimates to calculate D, giving two different critical lines.