Mimicking Nonequilibrium Steady States with Stochastic Pumps

We establish a correspondence between two very general paradigms for systems that persist away from thermal equilibrium. In the first paradigm, a nonequilibrium steady state (NESS) is maintained by applying fixed thermodynamic forces that break detailed balance. In the second paradigm, known as a stochastic pump (SP), a time-periodic state is maintained by the periodic variation of a system's external parameters. In both cases, currents are generated and entropy is produced. Restricting ourselves to discrete-state systems, we establish a mapping between these scenarios. Given a NESS characterized by a particular set of stationary probabilities, currents and entropy production rates, we show how to construct a SP with exactly the same (time-averaged) values. The mapping works in the opposite direction as well. These results establish an equivalence between the two paradigms, by showing that stochastic pumps are able to mimic the behavior of nonequilibrium steady states, and vice-versa.


Introduction and Motivation
While there is currently no single theory that unifies all nonequilibrium phenomena, a number of useful paradigms of nonequilibrium behavior have emerged. These include: small perturbations near equilibrium; systems driven away from an initial state of equilibrium; spontaneous relaxation towards equilibrium; non-equilibrium steady states generated by fixed thermodynamic forces; and stochastic pumps driven by the time-periodic variation of external parameters. Theoretical frameworks developed within each paradigm -for instance, linear response theory to describe near-equilibrium perturbations -have contributed to a broader understanding of nonequilibrium processes.
In this work we focus on two of these paradigms: non-equilibrium steady states and stochastic pumps. These share certain features, notably the persistence of non-vanishing currents and entropy production, which invite a comparison between the two. As elaborated below, we will devise a mapping from one paradigm to the other: given a system in a nonequilibrium steady-state characterized by certain occupation probabilities, currents and entropy production rates, we will show how to construct a stochastic pump that exhibits the same (time-averaged) properties. The inverse direction, namely the construction of a nonequilibrium steady state with the same properties as a given, time-averaged stochastic pump, will also be discussed.
In the nonequilibrium steady-state (NESS) paradigm, a system driven by fixed thermodynamic forces -such as temperature gradients or chemical potential differences -reaches a steady state in which its statistical properties are stationary with time. Unlike an equilibrium state, a nonequilibrium steady state exhibits non-vanishing currents, reflecting the violation of detailed balance. In order to maintain such a state, a thermodynamic cost must be paid. This cost is measured by the continual depletion of a thermodynamic resource, such as a chemical fuel, resulting in the production of entropy in the system's thermal surroundings.
Biomolecular motors illustrate the NESS paradigm [7,8]: a reaction such as AT P hydrolysis (AT P → ADP + P i ) produces entropy in the surrounding solution, and the chemical potential difference between reactants and products provides the thermodynamic force. The "current" in this situation corresponds to the mechanical motion produced by the motor, for instance the systematic displacement of kinesin motor toward the positive end of a microtubule filament. For a recent review of the stochastic theory of nonequilibrium steady states, as applied to biochemical processes, see Ref. [5].
In the stochastic pump paradigm, a system is driven by the time-periodic variation of external parameters in the presence of a thermal reservoir. Typically, it is assumed that the dynamics satisfy detailed balance at every instant in time -in other words, if the parameters were suddenly frozen at their instantaneous values, the system would relax to an equilibrium state. Under suitable conditions, a periodically driven system reaches a time-periodic state with nonvanishing time-averaged currents. These currents are effectively "pumped" by the periodic variation of the parameters, and the cost associated with pumping these currents is the work invested in driving the parameters. Ultimately, the energy provided by this work is dissipated into the thermal reservoir, resulting in the production of entropy.
The study of stochastic pumps has been stimulated by experiments on artificial molecular machines [3,6], which are manipulated by the variation of external parameters to achieve some desired behavior. For instance, in experiments on catenanes -mechanically interlocked ring-like molecules -the aim was to produce unidirectional rotation of one ring around the other [9]. Theoretical investigations of SP have focused on slowly driven [2,17] as well as weakly driven [18] pumps, "no-pumping" theorems [12,14,4,10,11], the role of interactions [1] and fluctuations [15], and the ability to extract work from stochastic pumps [13].
Underlying both the experimental work on artificial molecular machines and the theoretical work on stochastic pumps is the broad goal of understanding how to achieve controlled motion at the molecular level, where thermal fluctuations are large. The focus on time-dependent driving is motivated in part by the difficulty of synthesizing artificial molecular systems that, like biological molecular motors, takes advantage of chemical potential differences to drive steady currents. It is often simpler to manipulate the system by varying external parameters such as temperature and the surrounding chemical environment.
In both nonequilbrium steady states and stochastic pumps, the generated currents can be viewed as desired outcomes, and entropy production as the cost of achieving them. In this perspective, fixed thermodynamic forces (NESS) and time-periodic external driving (SP) represent the tools at our disposal. It is then interesting to compare these two sets of tools with respect to the degree of control that can be achieved. In particular, in this paper we investigate whether time-periodic driving can always achieve the same time-averaged outcome (i.e. identical currents) as fixed thermodynamic forces, and at the same cost (identical time-averaged entropy production). In other words, can a stochastic pump "mimic" an arbitrary nonequilibrium steady state?
More precisely, we begin by considering a generic Markov process of random transitions among n states of a system. The transition rates are fixed in time and do not satisfy detailed balance, hence the dynamics lead to a NESS with nonvanishing currents and entropy production. We then show how to prescribe a stochastic pump that has, in the limit of many cycles, the same time averaged probabilities, currents and entropy production rates as the NESS. This prescription is constructive though not unique. Surprisingly, the construction does not require the solution of any differential equations, only linear, algebraic equations. By contrast, a mapping in the opposite direction (from SP to NESS) requires that we first determine the periodic state of the system, which involves solving a set of coupled ordinary differential equations with time-periodic parameters. Typically this can only be done numerically.
The rest of the manuscript is organized as follows: In section 2 we formulate the problem and review some useful known results. We then consider in section 3 the simpler direction, namely mapping a stochastic pump into a NESS. In section 4 we present two types of transformations for detailed balance matrices which play a key role in our construction of a SP that mimics NESS, and the "no current-loops" property that sets a constraint on pumping with time-dependent detailed balance matrix. The construction of a pumping protocol is described in section 5. We finish with few concluding remarks in section 6.

Definitions
We limit ourselves to an ergodic, continuous-time Markovian system with n states. The evolution of the system consists of random, Poissonian transitions among these states, with transition rates that are governed either by a timeindependent rate matrix R (when analyzing nonequilibrium steady states) or a time-periodic rate matrix W(t) (for stochastic pumps), as described in more detail below. It is convenient to picture the system in terms of a "connectivity graph" with n nodes and a finite number of edges connecting given pairs of nodes. An edge between states i and j implies that the system can make transitions between these states (see Fig. 1). We will use p(t) to denote the vector of the probabilities whose i'th component p i (t) is the probability for the system to be in state i at time t.
In the NESS scenario, the evolution of the system obeys the master equation Figure 1: A four-state system described by a graph. Each node represents a state of the system. The edges represent non-vanishing transition rates between states. In this example, direct transitions between states 1 and 3 are not allowed.
For i = j, the matrix element R ij is the probability per unit time for a system in state j to make a transition to state i. The diagonal elements of R are negative, and are determined by conservation of probability: i R ij = 0. For the system to be ergodic, we demand that (i) if R ij = 0 then also R ji = 0, and (ii) the graph associated with R is connected. That is, for any pair of nodes (states) i and j, there exists a path from i to j, possibly through a sequence of intermediate nodes, along the edges of the graph. Under these conditions, Eq.(1) has a unique steady state solution, which we will denote by p ss , and any solution of Eq.(1) converges to this steady state in the long-time limit [16].
In addition to the steady state probabilities p ss , we will be interested in the steady state currents, defined by and the entropy production rates associated with these currents [16]: We will assume that some of these currents, and therefore the corresponding entropy production rates, are non-vanishing. In other words we assume the dynamics violate detailed balance, hence p ss is a genuinely nonequilibrium steady state.
In the stochastic pump scenario, the system obeys a master equation with a time-periodic rate matrix, where for some finite period T . If we momentarily treat t in Eq.(4) as a parameter of the rate matrix (rather than as the time variable), then for any fixed value of this parameter we will assume the rate matrix W(t) has a unique stationary solution π(t), and we further assume that In other words, the dynamics generated by W(t) (for fixed t) satisfy detailed balance. We will refer to π(t) as the equilibrium state of W(t). This is the state to which the system would relax if all the transition rates W ij were "frozen" in time.
Let us now return to thinking of t as time. For any solution of Eq.(4) the quantities (suppressing the argument t on the right side) represent instantaneous currents and entropy production rates. Under Eqs. (4,5) the statistical state of the system evolves asymptotically to a unique time periodic state, with currents and entropy production rates These are analogous to the quantities appearing in Eqs. (2,3), only J ps ij (t) and σ ps ij (t) are periodic with time, whereas J ss ij and σ ss ij do not vary with time. Throughout this paper we will be interested in the asymptotic properties of the system, therefore we will consider only the steady state p ss and the periodic state p ps (t), and not the process of relaxation to either of these states.
In order to compare the NESS and stochastic pump (SP) scenarios, let us define the time averaged quantities in the periodic state of the SP: The problem that we wish to study can now be formulated as follows.
Problem Formulation: Given a time-independent rate matrix R corresponding to steady state quantities p ss , J ss and σ ss , we want to construct a time periodic detailed balance rate matrix W(t) whose periodic state is described by the same quantities, after averaging over time: We denote the above problem as the "forward" problem. We will also be interested in the "inverse" problem: given a time-dependent detailed balance rate matrix W(t) corresponding to the time-averaged quantities p ps i , J ps ij and σ ps ij , we want to construct a time-independent rate matrix R such that Eq.(13) holds. As we discuss in more detail below, the solution of the inverse problem follows directly from known results, therefore we will focus mainly on the forward problem in this paper.
When R and W(t) give rise to dynamics that satisfy Eq.(13), we will say that the stochastic pump "mimics" the nonequilibrium steady state, and vice-versa.

Two Useful Decompositions
The two well known decompositions described below, the first -an algebraic decomposition of rate matrices, and the second -a topological decomposition of the connectivity graph, will be extensively used in what follows.

Rate matrix decomposition
The following (unique) decomposition of any rate matrix R, obtained by Zia and Schmittmann [19], will prove to be useful: Here the multiplication is ordinary matrix multiplication, and S is a symmetric matrix whose elements in each column add up to zero, with negative entries only on the diagonal. J ss is the anti-symmetric current matrix defined in Eq.(2) and P = diag( p ss ) is the diagonal matrix with elements P ii = p ss i . An immediate corollary of Eq. (14) is the following statement: if a rate matrix R is the product of a symmetric rate matrix S and a diagonal matrix P (with positive diagonal entries summing to unity), then R satisfies detailed balance, i.e. there are no currents in the stationary state.

Cycle Decomposition
The currents that characterize a NESS are in general not independent of one another, as they must satisfy constraints arising from the conservation of probability. These constraints embody Kirchoff's law of currents. The cycle decomposition method provides a convenient tool to account for these constraints [16]. Briefly, in a connected network with N nodes and E edges, the conservation of probability imposes N − 1 constraints among the E currents (one current per edge). It is then convenient to identify C = E − N + 1 fundamental currents, using the following procedure. First, we build a connectivity graph for the system, as described above and illustrated in Fig.(2) for the case of four nodes and six edges (hence C = 3). We then construct a maximal spanning tree, by removing C edges without breaking the connectivity of the graph; this tree, illustrated by the red dash-dotted lines in Fig.(2), has no cycles. In the context of the original graph, the C edges that are removed to form the spanning tree are called fundamental edges. The currents along these edges are the fundamental currents (the black arrows in the figure), and the currents along the edges of the tree are the spanning tree currents (the red arrows in the figure).
For NESS, the steady state currents along the C fundamental edges can take on any values, independently of one another. However, once these fundamental currents are set, the spanning tree currents are then uniquely determined by conservation of probability: the total sum of incoming and outgoing currents at each state must vanishes in the steady state. Therefore, the number of degrees of freedom in the matrix J ss is not (n 2 −n)/2 as for an arbitrary anti-symmetric matrix, but is determined by the graph topology.
Unlike NESS, for stochastic pumps the fundamental currents at each moment do not fix the currents on the spanning tree edges since probabilities can temporarily accumulate on the vertices (Kirchoff's current law does not apply at any instant of time). However, J ps on the fundamental edges do dictate J ps on the spanning tree, since the average probabilities are conserved.

Mapping SP to NESS
In this section we consider the inverse problem defined at the end of Sec. 2.1, which conceptually is the simpler direction: given a time-dependent periodic rate matrix W(t), how do we construct a time-independent rate matrix R whose steady state properties satisfy Eq. (13)?
To construct R, we first solve for the periodic averages p ps i , J ps ij and σ ps ij associated with W(t). These can be calculated by obtaining the periodic solution of the master equation, p ps i (t), and then plugging this solution into Eq. (12). In most cases, however, finding the periodic solution must be done numerically.
Once p ps i , J ps ij and σ ps ij are known, we next have to build a rate matrix R whose steady state properties, p ss , J ss and σ ss , satisfy Eq. (13). A nice consequence of Eq. (14) is that if all the currents along the graph edges are nonzero, then the quantities p ss , J ss and σ ss uniquely determine R. We express this relationship by the shorthand notation To see this, we use Eq. (14) to write the entropy production rates as By this equation, the elements of J ss and σ ss uniquely determine S, if the currents are non-zero: If some of the currents are zero, then the corresponding elements of S are not uniquely determined from the currents and entropy production alone, and they can be arbitrarily chosen. Once S has been obtained consistently with J ss and σ ss , it can be combined with P and J ss , via Eq. (14), to give R.
In the remainder of this paper, we will address the forward problem, namely how to construct, for a given NESS, a mimicking SP protocol. While this problem is conceptually more complicated than the inverse problem discussed above, it turns out that it is computationally much simpler and does not require any solution of differential equations.

Key Ideas
Here we establish three technical results that will play a crucial role in the construction of W(t).

No Current-Loops in Detail Balance Systems
A system satisfying detailed balance has non-vanishing currents when the instantaneously probability distribution differs from the equilibrium state of the instantaneous rate matrix. These currents, however, cannot form a current loop. That is, no loop i, j, k, · · · m, i on the graph associated with the system can have all the currents oriented in the same direction around the loop. As an example of a current loop, consider the system described in Fig. (2). No detailed balance system can have instantaneous currents equal to the currents in the loop 1 → 3 → 4 → 1 (or in the loop 2 → 3 → 4 → 2), since they all have the same orientation.
To see why a current loop is inconsistent with detailed balance, let us consider the currents along the edges of a loop i, j, k, ...m, i. Using Eq. (14) to decompose the detailed balance rate matrix W = SΠ −1 , the currents generated by a distribution q satisfy . . .
Summing these equations we get that (J /S) = 0 around the loop. This means that not all the currents can have the same sign, since each S ij > 0. Therefore there are no current loops. (A similar argument was used in [12]).
By contrast, the steady state currents of a NESS must form at least one current-loop. To see this, just choose a site (denoted by i) through which some of the steady state currents flow. Conservation of probability implies that at least one of these currents is going out of the site i, so we choose such a current say from i into site j. Now from the site j again, there is at least one current going out, say to to site k. Following the same argument, from any site we can "flow" with a current into a new site, but since the number of sites is finite, after no more than n such steps we must come back to a site we already visited. Therefore, there must exist at least one current-loop.

Transformation for the diagonal part of the decomposition of W
Suppose we have a detailed balance rate matrixŴ and two instantaneous probability distributions q and p, neither of which necessarily correspond to a stationary distribution. We would like to transformŴ into a different detailed balance rate matrix, W, such thatŴ q = W p. This is achieved by the following transformation where P and Q are the diagonal matrices corresponding to p and q, respectively. Using indices, this reads The transformation given by Eq.(21a) affects only the diagonal part of the decomposition and has the following properties: 1. IfŴ satisfies the detailed balance condition, then so does W. This can be seen by using the decomposition ofŴ as in Eq. (14) in the above transformation and noting that W is a symmetric rate matrix times a diagonal matrix, and therefore it has no currents in its steady state.
2. The instantaneous currents of a system described byŴ with probabilities q are the same as those of a system described by W with the probabilities p. This follows by substituting Eq.(21b) into Eq. (7).
3. From Eqs.(21b, 8) it follows that the instantaneous entropy production rates along each edge (σ ij ) for a system described byŴ with probabilities q and for a system described by W with probabilities p are the same.
The significance of this transformation can be stated as follows. If we have a rate matrixŴ and probabilities q, which produce instantaneous currents J and entropy production rates σ, then for any other probability distribution p we can construct the rate matrix W that generates the same J and σ.

Transformation for the symmetric part of the decomposition of W
Currents arise in a system described by a detailed balance W when the instantaneous probability distribution p differs from the equilibrium distribution, π. We next show how to vary the magnitudes of these currents (but not their directions) while keeping p and π fixed. As the directions of the currents do not vary under this transformation, no loops can be formed in accordance with the "No loop condition" in Sec.(4.1). Let us use the decomposition W = SΠ −1 where S is symmetric and Π = diag( π). The currents are then given by We see that by varying S ij we vary the magnitude of the current J ij , but not its sign, since S ij ≥ 0. This transformation is complementary to the one given by Eq. (21): it enables us to change the currents (and entropy production rates) while keeping the probabilities p and π fixed, by tuning the symmetric part of W in the decomposition given by Eq. (14). By contrast, with the previous transformation we can vary the probabilities p and π at fixed currents and entropy production rates, by tuning the diagonal part of W.

Construction of the pumping protocol
Given a rate matrix R one can calculate its steady state p ss (the null eigenvector of R) and thus J ss and σ ss using Eqs. (2,3). Our goal is to construct a periodic pumping protocol -a time-dependent detailed balance rate matrix W(t) with some period T -such that Eqs. (13) hold, or in other words a SP that mimics the NESS. For simplicity, we assume that there are no edges along which the steady state currents are zero. The case involving zero currents along some edges is analyzed in the appendix.
We note that in general there are many SP's that mimic any specific NESS -the mapping is not one-to-one. Out of the many SP that mimic the NESS, we wish to choose one using a relatively simple construction, yet generic enough to mimic any NESS. Naively we would like to have a SP whose periodic state gives rise to time-independent quantities, p ps i (t) = p ss i , J ps ij (t) = J ss ij and σ ps ij (t) = σ ss ij . Such a construction is, unfortunately, impossible. This can be seen from the result of Sec. (4.1), which states that for all t, J ps ij (t) cannot have any current loops, whereas J ss ij must have at least one current loop. Thus, we can only hope to achieve a mapping between the time-averaged quantities associated with the SP and those of the NESS, as in Eq. (13). This also implies that at least some of the currents of the SP must be time-dependent.
The construction described below, though the simplest we could find, is nevertheless somewhat convoluted. We therefore first give an overview before proceeding to the detailed description. The main reason for the complication is the fact that the periodic solution, p ps i (t), is a highly non-trivial function of the pumping protocol W(t). To avoid this complication, we construct simultaneously both the pumping protocol W(t) and its periodic solution p ps i (t). This is achieved in 5 steps.
In the first step, described in Sec.(5.1), we divide the pumping protocol temporal interval of duration T into two equal half-intervals, designated as a and b. We then assign an arbitrary, fixed detailed balanced matrixW a and an arbitrary, fixed probability distribution q a , to be associated with the first half-period of driving. We will refer toW a and q a together as the seed for that half-cycle, and these quantities will be used to set the current directions during that time interval. Next, the seed and current directions for the first half-cycle are used to assign a seed (W b , q b ) and current directions for the second halfcycle. The current directions during the first half-cycle are opposite of those of the second half-cycle, and both sets are, by construction, consistent with the "no current loop" condition in Sec. (4.1).
In the next two steps, we use the seeds to construct fixed sets of currents (J a ij and J b ij ) and entropy production rates (σ a ij and σ b ij ) for the first and second halves of the cycle, whose averages over the two halves are equal to the steadystate values that we wish to mimic: for all i = j. This is done first for the fundamental currents in Sec. (5.2) and then for the spanning tree currents in Sec.(5.3). Up to this point, the currents and entropy production rates for the two half-cycles have been constructed from the initial, arbitrary seeds, but the corresponding time-periodic rate matrices W a (t) and W b (t) that actually generate these currents and entropies are not yet known. In the fourth step, described in Sec. (5.4), we use the transformation of Sec. (4.3) to adjust the symmetric part of the seed matrices,W a,b , arriving at new rate matricesŴ a,b that produce the desired currents J a,b , for the seed probabilities q a,b . These currents, together with the desired averaged probabilities p ps i , fix p ps i (t). In the last step (Sec. 5.5) we use the transformation of Sec. (4.2), together with the symmetric parts ofŴ a,b , to construct rate matrices W a (t) and W b (t) for which p ps i (t) is the periodic solution of the master equation.
In the specific protocol described below, the entries of the matrix W(t) are not continuous functions of time, as they have discontinuities between the two T /2 intervals. These discontinuities are not essential, and can be removed at the expense of making the construction less transparent.
To improve the clarity of presentation, some of the formal definitions of the construction are followed by a concrete application to the example of a 4-state system described in Fig.(2). In this example the NESS system has 4 states with p ss = (0.1, 0.2, 0.3, 0.4). The fundamental currents were chosen such that they form a loop, and their values are J ss 31 = 3, J ss 43 = 2 and J ss 14 = 1. The currents for the spanning tree edges are then dictated by Kirchoff's law -the sum of currents in each vertex must be zero. The corresponding current-matrix is Finally, we choose the entropy production rate to be 1 along all the edges. Using Eq.(16) the matrix S can be calculated, and is given in Eq.(42) in the Appendix. The matrix R giving rise to this particular NESS can be constructed using Eq.(14).

Step 1-choosing the seed
In what follows, superscripts a and b stand for quantities associated with the first and second halves of the period, respectively. Our first step is to choose an arbitrary detailed balance matrix on the graph. This is done by choosing an equilibrium state for the first half period π a and a symmetric rate matrixS from which we composeW a =S(Π a ) −1 where Π a = diag( π a ). Next, we choose a fixed probability distribution q a = π a that satisfies log for any i = j. It is always possible to satisfy this condition, by choosing q a close enough to π a . As we will see in the next section, this condition is necessary for the consistency of our construction. For the second half of the period we replace π a and q a by vectors with components π b i = 1/π a i -from which we constructW b =S(Π b ) −1 -and q b i = 1/q a i . The two probability vectors π b and q b are not normalized, but as will become clear, this normalization does not play any role, and it will prove to be simpler to work with these unnormalized vectors. We note that the currents generated byW a and q a ,J a ij =S ij π a j q a j − π a i q a i , have, on each edge, opposite signs to those generated byW b and q b , given bỹ The values ofJ a,b ij will not explicitly be used in what follows -only their directions, which by construction are consistent with detail balance. We additionally note that To illustrate this part of the construction with our four-state example, we chooseS BothS and π a are, in fact, arbitrary. We next note that min ij [log R ij p ss j − log R ji p ss i ] = 1/3. If we therefore choose q a such that |log π a i − log q a i | < 1/6,then Eq.(25) is satisfied. Any vector close enough to π a will do. As an example we use q a = (0.23, 0.24, 0.26, 0.27). In the second half of the period these correspond to π b = (4, 4, 4, 4) and q b = (4.3478, 4.1667, 3.8462, 3.7037).

Step 2 -fundamental edge currents
In this step we set the currents along the fundamental edges to be constant during each of the two half cycles, such that their time averages (the average between the first and second halves) is equal to J ss , and the time average of the entropy production rates is equal to σ ss . Moreover, we will choose these fundamental currents to have the same directions as the fundamental currents iñ J a during the first half-cycle, and to be reversed in direction during the second half-cycle.
We first make sure that for each fundamental edge, the average of J a ij and J b ij , and therefore the time averaged current, is exactly J ss ij . To that effect we introduce the following rule: • If the direction of the NESS current, J ss ij , is the same as the direction ofJ a ij defined above, then we set in the first half of the period J a ij = (2+α ij )J ss ij and in the second half of the period J b ij = −α ij J ss ij , where α ij are positive and will be determined below.
• If the direction of the NESS current, J ss ij , is not the same as the direction ofJ a ij , then we set in the first half of the period J a ij = −α ij J ss ij and in the second half J b ij = (2 + α ij )J ss ij .
According to this rule, the directions of the currents during the first half of the period are the same as that ofJ a ij and are opposite to those in the second half. Moreover, by the above construction the time averaged currents on the fundamental edges have the required values: Next we determine the values of the α ij 's so as to satisfy the requirement on the entropy production rates. Assuming for the moment that the probability distributions in the first and second halves of the period are given by q a and q b and that the equilibrium distributions of the detailed balance matrices during the first and second halves of the period are given by π a and π b respectively, then the entropy production rates with the currents J a ij and J b ij during the two halves of the cycle are given by Substituting in these equations J a ij and J b ij in terms of α ij and demanding that 1 2 (σ a ij + σ a ij ) = σ ss ij , gives an equation for α ij whose solution is Eq.(25) ensures that indeed α ij > 0.
To illustrate this step, we examine the signs of the currents generated by the matrixW a =S(Π a ) −1 and the probability q a , denoted byJ a , on the fundamental edges: sign J a 13 = +1, sign J a 14 = +1, sign J a 34 = +1.
While the direction of the current along the 1-4 edge is the same as the current orientation of J ss , for the other two fundamental edges the directions of J a and J ss are not the same. We next solve Eq.(32) for α ij . The explicit expression, as well as the currents, are given in the appendix. Fig(3) shows the currents on the first half period (left) and second half period (right). Note that (i) the direction of the currents along each edge are opposite in the two halves of the period (ii) as discussed above, in the first half period the direction of the current along the 1-4 edge is the same as that of Fig(2), but the direction of currents along the other fundamental edges are different from that of Fig(2), and (iii) there are no current loops in Fig(3).

Step 3 -Spanning tree currents
So far we have shown how to construct the fundamental edge currents. Next we discuss the currents along the edges of the spanning tree. For both the first and the second half-periods we impose the following two constraints: 1. The sum of currents feeding into any site i during the first half-period must be equal to minus the same quantity during the next half-period: These constraints ensure that we indeed have a periodic time evolution: Note that these are only n − 1 independent equations, since conservation of probability adds the constraint ij J a,b ij = 0 . 2. For each spanning tree edge we demand that The number of edges in the spanning tree is n − 1, and therefore these are n − 1 additional equations. They ensure that the time-averaged entropy production rate along the i, j edge are the same as σ ss ij . All together Eqs. (34, 35) are 2(n − 1) linear equations for 2(n − 1) unknowns: J a,b on the spanning tree, which has n − 1 edges. Moreover, the directions J a,b ij that solve these equations are the same as that ofJ a,b ij respectively. To see this, let us use the definition of q b and π b in the second condition above, together with the definition of σ ss ij (Eq.3): Taking the absolute value of both sides in the above equation and using Eq.(25) implies that However, by the construction of J a ij and J b ij their average is equal to the steady state current, J ss ij , therefore Eqs.(38, 37) are consistent with each other only if the sign of J a ij is opposite to that of J b ij . Moreover, as σ ss ij > 0 (this follows from Eq.(16)), the signs of J a ij and J b ij must be the same as that of log π a i q a j /q a i π a j and log π b i q b j /q b i π b j , respectively, otherwise the left hand side of Eq.(35) will be negative. Therefore it is also the same as the sign ofJ a ij .

Step 4 -Transforming the symmetric part for the currents
At this stage we have constructed the currents for both half cycles with the same directions as the currents ofW a,b with q a,b . We can now use the transformation of the symmetric part of the decomposition of the rate matrix (defined in section 4.3) to change the symmetric part ofW such that the currents generated by the transformed matrix and the seed probabilities q a,b have the constructed values: Importantly, all these off-diagonal elements are positive. This follows from the fact that the denominators are justJ a,b ij /S ij , but by our construction the signs of J a,b ij is the same as that ofJ a,b ij . The diagonal elements of S a,b are now determined by the requirement that the sum of each of the columns is zero.

Step 5 -Forming a solution to the master equation
We have now arrived at the matricesŴ a,b = S a,b (Π a,b ) −1 . These are the transformed seed matrices, which have been constructed (by tuning the elements of S a,b ) so as to produce the desired currents and entropy production rates when the probability vectors are q a,b . However, the time average of q a,b is not p ss , and in fact q a,b are not solutions of the master equation: ∂ t q a,b =Ŵ a,b q a,b .
To remedy this situation we use the transformation of the diagonal part of the decomposition ofŴ, defined in section 4.2.
First, we want the time-averaged probabilities to be equal to p ss . Second, we already know what ∂ t p i should be, in terms of the desired currents J a,b . Namely, ∂ t p i = j J ij . Therefore we do the following: 1. Calculate m a,b = ∂ t p =Ŵ a,b q a,b . These are the temporal slopes of the probabilities that solve the master equation, during the first and second halves of the period.
2. Choose T such that 0 < p ss i ±(T /4)m a,b i < 1 for any i. This choice ensures that the probabilities stay bounded between 0 and 1, and is always possible by taking T to be small enough.
The matrix W(t) has all the periodic state averages we demand, and its periodic state solutions in the two halves of the period are, by construction, p a,b (t).
For our example, we calculate the slopes ∂ t p on the two half cycles by proper summation of the currents: ∂ t p a = (37.89, −2.51, 5.85, −41.22) and   Fig.(4). Plugging these into Eq.(40) gives W(t). To verify that the solution of the master equation with the constructed W(t) has the required properties we solve this system numerically. The numerical results for p(t), shown as red circles in Fig(4), agree (up to numerical error) with the analytical solution.

Conclusions
We conclude with a few comments on our construction. First, the protocol presented above is clearly not unique. For example, different seeds or choices of a spanning tree result in different protocols. The non-uniqueness might be used, in principle, to match additional quantities, e.g. fluctuations around the average or the rate of decay towards the steady and periodic states. In addition, as mentioned previously, in our construction W(t) has discontinuities between the two halves of the period. This results from discontinuities in S, Π and Q. The discontinuities in S can be avoided if the currents are not taken to be fixed during the two halves of the period, but changing with time, such that the currents at t = T /2 vanish. The discontinuities in Π and Q can be avoided if we change q as a function of time, crossing π at t = T /2. This, however, makes the construction more cumbersome. Next, we note that in the above construction both the symmetric and the diagonal parts of the decomposition of W change with time. This is known to be an essential feature of all pumping protocols [14].
Lastly, we note that for NESS there is no minimal entropy production rate associated with a given set of currents, as is evident from Eq.(16): for any J the entropy production rates of the NESS can be made arbitrarily small by taking S to be large enough [19]. The mapping presented here shows this is also the case for stochastic pumps: finite currents can be pumped with arbitrarily small values of dissipation, which is somewhat surprising. Stated more generally, for any connected graph, both a NESS and a SP can be constructed to have any desired set of (time-averaged) probabilities, non-zero currents and positive entropy production rates, provided the currents obey Kirchhoff's law.
Let us now consider the construction of a stochastic pump with small entropy production rates but large currents. For the entropy production rates of the NESS to be very small with non-vanishing currents, log R ij p ss j /R ji p ss i must be very small, say of order ε. Eq.(25) then implies that log π a i q a j /π a j q a i ∼ O(ε) as well. Exponentiating this gives π a i q a j /π a j q a i = 1 + O(ε) thus (π a j ) −1 q a j − (π a i ) −1 q a i ∼ O(ε). By Eq.(39) this implies S ∼ O(ε −1 ). Therefore, in stochastic pumps -as in NESS -small entropy production rates with non-vanishing currents come at the cost of large values in the symmetric part of the rate matrix. But when the elements of S are large, the corresponding transitions occur very rapidly. Thus, we obtain finite currents at arbitrarily low dissipation when there is a large separation of timescales, with transitions among the n states of the system -and therefore relaxation to equilibrium -occurring much more rapidly than the external driving of parameters. The greater the separation of timescales, the more the stochastic pump approaches the adiabatic (quasistatic) limit [2], in which the system remains in thermal equilibrium at all times and there is no entropy production.
spanning tree edges are dictated by Kirchoff's law -the sum of currents in each vertex must be zero. The corresponding currents matrix is To set S, we choose the entropy production rate to be 1 along all the edges. Using Eq. (16) We next note that min ij [log R ij p ss j − log R ji p ss i ] = 1/3. If we therefore choose q a such that |log π a i − log q a i | < 1/6, condition Eq.(25) is automatically satisfied. Any vector close enough to π a will do. As an example we use q a = (0.23, 0.24, 0.26, 0.27). In the second half of the period these correspond to π b = (4, 4, 4, 4) and q b = (4.3478, 4.1667, 3.8462, 3.7037).
Note that the direction of the current in the 1−4 edge is the same as the desired current orientation, but for the other two edges the directions are different from the desired currents.
Together with the directions the fundamental currents in the two half cycles are given by Note that the averages are the same as the desired values.
In the third step we solve for the spanning tree currents, by solving Eqs.(34,35) for J 1,2 12 , J 1,2 23 and J 1,2 24 . They are given by In the fourth step we calculate S a and S b using Eq.(39). We present here only the upper part of them as they are symmetric: In the last step, we first calculate ∂ t p on the two half cycles by proper summation of the currents: ∂ t p a = (37.89, −2.51, 5.85, −41.22) and ∂ t p b = (−37.89, 2.51, −5.85, 41.22), which as expected cancel each other. These are the slopes of the p(t)'s in the first and second halves of the period. To keep 0 < p(t) < 1, we need to choose T small enough, say for simplicity T = 0.01. Using T we can calculate p a,b (t) which are the actual probability distributions in the two halves. These linear functions are plotted in Fig.(4). Plugging these into Eq.(40) gives W(t).

B Zero currents in the NESS
If J ss ij = 0 for some edges, then clearly σ ss ij = 0 for the same edges as well. However, since σ ps ij (t) ≥ 0 identically (this follows from Eqs. (10) and (11)), we must have σ ps ij (t) = 0 for all t in order for the time-averaged entropy production rate of the periodic state to be zero. This can only happen if the currents along these edges are zero at all time, J ps ij (t) = 0. A simple prescription to set these currents to zero is to use the following modified rate matrix: where P is the matrix with the steady state of R on its diagonal, J ss is steady state current matrix of R (but with different factor in front of it), and S ij = 0 J ij = 0 and i = j with the diagonal elements ofS changed to make the sum of columns is zero. This modified rate matrix, in which the edges with zero steady state currents have been "removed" (R ij = 0 on these edges), has the same steady state as R, but its steady state currents are 1.5 larger then those of R. From Eq.(16) it also follows thatσ ss ij = 1.5σ ss ij .
Note that forcing some of theR ij to be zero might make the time-dependent system non-ergodic, since it might disconnect some of the states from the others at all times. This, however, can be overcome by dividing the time interval into three equal intervals rather then two. In the first two parts we repeat the construction as before, but usingR instead of R. In the last interval, we choose W(t) = S c (Π c ) −1 , with S c ij = 1 for any i = j and Π c the diagonal matrix with p ps (t = 0) on its diagonal. With this construction, there are no currents in the periodic solution during the last time interval, the average currents and entropy production rates are the required ones, and the last interval ensures that the system is ergodic.