Rectification of energy and motion in non-equilibrium parity violating metamaterials

Uncovering new mechanisms for rectification of stochastic fluctuations has been a longstanding problem in non-equilibrium statistical mechanics. Here, using a model parity violating metamaterial that is allowed to interact with a bath of active energy consuming particles, we uncover new mechanisms for rectification of energy and motion. Our model active metamaterial can generate energy flows through an object in the absence of any temperature gradient. The nonreciprocal microscopic fluctuations responsible for generating the energy flows can further be used to power locomotion in, or exert forces on, a viscous fluid. Taken together, our analytical and numerical results elucidate how the geometry and inter-particle interactions of the parity violating material can couple with the non-equilibrium fluctuations of an active bath and enable rectification of energy and motion.


I. INTRODUCTION
Identifying mechanisms that rectify stochastic fluctuations is a longstanding problem in non-equilibrium statistical mechanics [1,2]. The Feynman Ratchet and pawl model, and its associated generalizations, have elucidated how systems can rectify stochastic fluctuations and act as microscopic engines that perform work and exert forces [3]. Indeed such models have provided a framework to understand how biological molecular motors can convert the energy derived from the hydrolysis of energy rich molecules into mechanical work [4][5][6][7]. While these advances provide powerful insights into how single-body systems can rectify fluctuations, much less is understood about rectification via collective effects in many body systems, despite recent advances [1,8]. Understanding many body rectification can lead to methods to manipulate the flow of energy across materials without any imposed temperature biases [9][10][11][12][13] and potentially facilitate the development of design principles for constructing synthetic molecular motor analogues. In this paper, we show how such rectification can be achieved in parity violating many body interacting systems. In particular, we show how a parity violating metamaterial [14] can spontaneously rectify energy and motion in the absence of any imposed gradients, when it is allowed to interact with a bath of active particles [15][16][17].
Our choice of parity violating metamaterial is inspired by a recent work where a metamaterial composed of interacting gyroscopes was introduced [14]. These gyroscopic metamaterials were shown to support chiral topological edge modes whose origin lies in a violation of time reversal symmetry (TRS) in the microscopic equations of motion of the interacting gyroscopes [14,18]. Crucially, the TRS violation is controlled by a combination of the spin of the gyroscopes and the geometry of the lattice.
In this manuscript, we consider the effect of immersing such a TRS-violating metamaterial in a bath of active particles that violate the fluctuation dissipation relation [19]. Our central result is a demonstration that this combination naturally rectifies stochastic fluctua-tions present in the bath. Crucially, the rectification mechanism relies on both the TRS-violating properties of the metamaterial and the TRS violations implicit in the single particle fluctuations of the bath. Indeed, simply coupling the chiral topological metamaterial to a heat bath does not result in any fluxes or symmetry breaking on account of the Bohr-van Leeuwen theorem [20] that forbids any non-equilibrium currents in thermalized magnetic systems. Similarly, active baths coupled to passive metamaterials do not in general lead to rectification of fluctuations.
Our model system generically supports a directed flux of energy across the metamaterial network. Unlike conventional energy flows, our energy flow does not require a temperature gradient. This flux can be routed through an otherwise isolated elastic object with the active network acting as a current source. Finally, we show that the microscopic mechanisms responsible for this energy flow can also potentially allow the parity violating active network to swim in and exert forces on a viscous fluid.
We analytically and numerically demonstrate these results for a variety of network geometries. In particular, we construct an intuitive diagrammatic approach for computing the energy flux (and consequently the swim speed) that shows how our results can readily be applied to tailor flows in arbitrarily complex networks. Taken together, our results establish a new mechanism for rectification of energy, motion and forces in non-equilibrium parity violating metamaterials. Unlike existing prescriptions, our mechanism exploits inherent asymmetries in the geometry and interactions of the material to achieve rectification.
The remainder of this paper is organized as follows: In Sec. II, we introduce our model parity violating metamaterial and provide a microscopic definition for the energy flux. In Sec. III-V we analytically identify the ingredients for rectification of energy fluxes and construct a diagrammatic approach that reveals a relationship between energy flux and network geometries. Finally in Sec. VI-VII we show that when the particles transmitting the energy flux are allowed to interact with a viscous low Reynolds Averaged energy flux from numerical calculations for a few typical networks. The flux direction and pattern can be controlled by the network geometry. In these figures, gray lines and dots represent the mechanical equilibrium structure of the network, and the size of the arrows is proportional to the magnitude of the averaged energy flux. The arrows are colored blue if it is counter-clockwise (CCW), red if clockwise (CW), and gray for fluxes not on the boundary. The numerical calculations were performed with all parameters set to 1.
number (Re) fluid, the nonreciprocal motions responsible for the energy flux can be utilized to pump the viscous fluid.

II. MODEL SYSTEMS AND ENERGY FLUX
Our model TRS-violating metamaterial is a springmass network in which each particle feels a Lorentz-like force in addition to spring-like interactions with neighbors and stochastic forcing from an active bath [19] (FIG. 1a). The equation of motion is: , where e ij is the unit vector from the equilibrium position of i to that of j, B = −Bẑ, and the lorentzforce's customary electric charge-like factor is absorbed in B. The construction of our model system is motivated by the recently constructed topological gyroscopic metamaterials [14] and in the linearized regime, the equations of our model system are equivalent to the equations of motion of the gyroscopic metamaterials [21]. The last two terms describe the active bath, which consists of friction −γv i and an Ornstein-Uhlenbeck (OU) colored noise η i [19]. The colored noise has finite correlation time τ with statistics where, for fixed τ , the parameter T a controls the variance of the colored noise, and I is the identity matrix with appropriate dimensions. The time evolution of the OU noise can be described according to the following equation [22], where ξ i is a delta function correlated white noise with unit variance. The friction −γv i and OU noise η i break fluctuation-dissipation relation, thus driving the system out of equilibrium [19]. The active bath reduces to a thermal equilibrium bath in the τ → 0 limit. Since the particles in our model are tethered to their lattice sites, rectification of fluctuations, if any, does not result in any particle flows. Rather, rectified fluctuations can affect the transport of energy or heat. To study such phenomena, the observable we mainly focus on is the time-averaged energy flux between particles at steady state. For a system with pairwise interactions and onsite potentials, the energy flux J ij from particle i to j reads To arrive at this formula, we first define the energy of a particle as the sum of its kinetic energy, on-site potential energy, and one half of the bond energies [23]. Then we write down the energy balance relations using ideas from stochastic energetics [24]. Finally we identify the energy exchanged due to particle-particle interactions as the energy flux J ij (derived in detail in the supplemental material [25]). We note that the energy flux can simply be interpreted as the rate at which work is done on particle j by particle i. Since this microscopic work is due to particles' stochastic motions, rather than due to an external control, the energy flux can also be interpreted as a heat flux [23,24]. The averaged energy fluxes Eq. (4) are identically equal to zero for a system at equilibrium with Boltzmann distribution. Starting from the linear equations Eq. (1), (3), the energy fluxes can be solved numerically using methods introduced in [ . 1b. We see nonzero energy rectification or energy fluxes can be generated in our chiral active system, and the flux direction or pattern changes with the network geometry. Using a linear response theory, we now develop analytical expressions for the energy flux that reveal how a combination of chirality, non-equilibrium activity, and network geometry is responsible for generating energy fluxes.

III. LINEAR RESPONSE THEORY FOR EN-ERGY FLUX
We begin by writing the equations of motion, Eq. (1), in frequency space,z Here, we have represented the displacement of all particles using a column vector z = i |i ⊗ z i , with |i denoting the 2D subspace of particle i.z(ω) andη(ω) denote the Fourier transform of z and the OU noise η, respectively. G + is the response matrix, in which matrix K encodes all on-site and spring forces F according to F = −Kz, and A is an antisymmetric matrix Eq. (5) describes how the displacement responds to the noise. Following the procedure in [28], the flux defined in Eq. (4) can be expressed using G + as a spectral integral where A as is an antisymmetric matrix A as = − |i j| ⊗ e ij e T ji + |j i| ⊗ e ji e T ij . The response function G + (ω) has no pole in the lower-half complex plane, but the colored noise introduces one pole at ω = −i/τ . Using the residue theorem we get a compact expression for the energy flux Eq. (7) and (10) will serve as our starting point to understand the energy flux. While they contain all the information required to compute energy fluxes, they have limited utility as design principles. Indeed, as written down, they require that the flux be recomputed de novo for each new network geometry and non-equilibrium bath activity. In the next sections, we show that it is possible to expand Eq. (7) and (10) in forms that reveals design principles for controlling energy fluxes.
As an aside, one important property that can be directly obtained from a similar linear response analysis is that the energy fluxes satisfy Kirchoff's law, i J ij = 0. The Kirchoff's law shows that on average there is no energy exchange between particles and the active bath. To derive the Kirchoff's law, we calculate the average heat exchange between particle i and the active bath v i · η i − γv i · v i , and following procedures in [28], this heat exchange can be shown to be zero (Supplemental

Material [25]
). The Kirchoff's law puts a strong constraint on possible energy flux patterns among particles, and some corollaries immediately follow, such as networks with no cycles cannot have nonzero flux, and fluxes of all bonds in a polygon network (as in FIG. 1b) are equal.

IV. INGREDIENTS FOR ENERGY RECTIFICA-TION AND THEIR ROLES
Compared with an ordinary thermal spring-mass network, that supports no energy fluxes in its equilibrium steady state, our parity violating metamaterial contains two extra components, the Lorentz force and the correlation in the noise. We first show that these two components provide two necessary ingredients required to ensure energy rectification in our model. Next, we clarify the role played by the geometry of the network in controlling the energy flux.
A. Lorentz force and non-equilibrium activity are necessary for the generation of an energy flux We begin with Eq. (7) that represents the averaged flux in terms of functions J F T (ω), h(ω). In particular, we note that the function J F T (ω) is proportional to the energy flux at Fourier frequency ω in an isolated damped variant of our network while the function h(ω) is proportional to the the noise spectrum, η(ω)η * (ω) = 2γT a h(ω)/t. To generate a nonzero flux, or equivalently make the integral nonzero in Eq. (7), we need two requirements (FIG. 2). Firstly, J F T (ω) should not be zero everywhere. In the absence of a magnetic field, B = 0, it can be easily shown from Eq. (7) that J F T (ω) = 0. Indeed, the response function G + is symmetric or reciprocal in this case, and since A as is antisymmetric, the trace tr G + (ω)A as = 0 at all values of ω. Nonzero B breaks the reciprocity of G + , and can thus generate a nonzero J F T (ω).
Nonzero J F T (ω) alone does not ensure a nonzero averaged energy flux, we further require that h(ω) not be constant. Indeed, noise characterized by a constant h(ω) function corresponds to fluctuation dissipation preserving white noise. In such cases, our model system would be in equilibrium even in the presence of a magnetic field according to the Bohr-van Leeuwen theorem [20]. A nonconstant h(ω) function corresponding to a colored noise source, such as the fluctuation dissipation violating OU noise, h(ω) = 1/(1 + ω 2 τ 2 ), can support a non-zero average energy flux.
In summary, we see that B-field and a colored noise are two necessary ingredients to generate flux in our model chiral systems. The role of the B-field is to break the reciprocity of response and generate Fourier modes such that J F T (ω) = 0. The role of the colored noise is to excite these modes in a weighted manner.
B. Energy flux can be tuned as a function of lattice geometry Apart from the these two ingredients, the geometry of the network also plays an important role. Indeed in the small γ regime, the existence of chiral modes with J F T (ω) = 0 can be heuristically explained by exploiting the connection between the slightly damped isolated variants of our system and the undamped isolated gyroscopic metamaterials [14,18]. Specifically, the slightly damped variant resonate near the eigen-frequencies of the undamped metamaterials, and hence exhibit Fourier modes that are close to the eigenmodes of the undamped system. Consequently, we infer that J F T (ω) = 0 as long as the corresponding eigenmodes in the undamped variant are chiral. As discussed in [14], the geometry of the network plays a crucial role in generating the chiral eigenmodes. At larger γ's, the Fourier modes of our damped isolated variant are no longer close to the chiral eigenmodes of gyroscopic metamaterials, but an emergent connection between eigenmodes and energy fluxes can be built (Supplemental Material [25]). In the next section, we further elaborate the role of geometry. Specifically, the central results of the next section provide compact expressions that elucidate the role played by geometry, Lorentz forces Each path is a pictorial representation of one term in the small-k expansion, and is depicted using green arrows. The magnitude of path of length n is on the order of k n . (b) Schematic of a polygon path and its outer angles θ1, θ2, . . . , θn. The flux of this path is simply Eq. (13). (c) For flux in complex networks, the leading order term is determined by the shortest cycles. Flux in the triangle part has order k 3 , and the prefactor c3 is the same as that in a standalone triangle network. Likewise for the pentagon part. As a result, the flux in a complex network can be viewed as a combination of fluxes in its constituent cycles. and non-equilibrium activity.

V. RELATIONSHIP BETWEEN FLUX AND NETWORK GEOMETRY
One might expect, as is the case in ordinary Newtonian mechanical metamaterials, that the energy flux in our active system be simply given by considering the eigenmodes of the damped isolated system and adding the flux carried by each of these eigenmodes, weighted by their degree of excitation by the colored noise. This is however not the case. In our system, flux carried by a pair of excited eigenmodes is not equal to the sum of the fluxes carried by each eigenmode individually because the friction and the random forcing can induce cross-couplings between the different eigenvectors. While the flux can be described in terms of the eigenmodes of a modified isolated system, as described in the previous section, such a connection still does not provide an efficient framework for controlling and manipulating the energy fluxes using the network geometry and activity.
In this section we develop a diagrammatic technique, which provides a simple intuitive method to compute energy fluxes and elucidates the relationship between flux and network geometry. The diagrammatic technique is constructed by expanding the expressions for the energy flux Eq. (10) in small-k regime and shows how the energy flux across a bond can be expressed as a sum over paths traversed along the network (Eq. (11)). Our perturbation theory assigns a geometry dependent pre-factor for each path, thus elucidating the role played by network geometry in ensuring rectified energy fluxes (Eq. (13)). Together, the central results of this section, summarized in Eq. (11), (13), provide compact expressions that elucidate how geometry, B-field, and correlation time τ of the colored noise can combine to generate energy flows in networks with arbitrarily complex geometry and topologies. Crucially through this diagrammatic approach we demonstrate that the energy flux is controlled by the local geometry of the lattice. This result is surprising given the non-local dependence suggested by Eq. (10). The diagrammatic techniques developed here enable the efficient control energy flux patterns in arbitrarily complex networks.

A. Path summation and its rules
Starting from the flux formula Eq. (10), we expand the flux to different orders in the spring constant k. Then for each order, we further collect contributions from various paths. In the final result, we write the total flux as a sum over the flux of paths (FIG. 3a), (Supplemental Material The path rules are as follows. For the flux from i to j, valid paths are l = i → j → l 3 → l 4 → · · · → l n → i, where l a and l b either has to be bonded or l a = l b . Paths that contain equal numbers of i → j and j → i do not contribute (e.g. path i → j → i), because either the path itself vanishes or it cancels with another path. As a result, paths appear as cycles. The term S l is defined as We note that all parameters (m, k g , k, B, γ, τ ) are condensed into this single angle α. (K s ) l b la ≡ l b | (K − k g I) |l a /k is a 2 × 2 submatrix used to calculate the non-dimensionalised spring force on l b due to the displacement of l a . −l means l in the reversed order. The interval of convergence depends on the geometry of the whole network as well as the parameter α. The typical value of the upper bound of k/k 0 ranges between 0.3 and 0.6. The paths can be represented using diagrams, from which the flux J path l can be calculated easily. For instance, the first diagram in FIG. 3a represents the path 1 → 2 → 3 → 1.
To calculate S l , one writes (−K s ) l b la for each arrow l a → l b , R α for each node l a , then multiply these matrices in the reversed order, and calculate the trace, e.g.
To get S −l , one takes the result of S l and replace α by −α. Finally, J path l can be calculated from the difference between S l and S −l .

B. Contributions to energy flux from leadingorder paths
Paths of length n show up as terms of order (k/k 0 ) n in the energy flux expression. In the small k regime, the main contribution to the flux comes from the lowestorder paths. The usual lowest-order paths are polygonal cycles with no loops (loops are self-connecting edges like l a → l a ). For these polygonal paths, the flux formula Eq. (11) reduces to a simple form (Supplemental Material where α is defined in Eq. (12), n is the number of nodes and θ i 's are outer angles (FIG. 3b). Eq. (13) illustrates how geometry of the network, as characterized by the angles θ i , together with the condensed parameter α that encodes the nonreciprocity due to the B field and the violation of fluctuation dissipation due to the colored noise, combine to generate energy fluxes.
For polygon networks, Eq. (13) gives a direct relationship between the lowest-order flux and the network geometry. As an example, flux in an arbitrary triangle is J ∝ k 3 sin θ 1 sin θ 2 sin θ 3 + O(k 4 ), whose k 3 term is always positive or CCW. For complex networks, Eq. (13) implys that its lowest-order flux can be viewed as a result of combining the flux of its constituent polygons, as illustrated in FIG. 3c. This is because the flux of a polygonal path Eq. (13) is not affected by any side chains on the nodes, and J path polygon for a polygon in a complex network is the same as J path polygon for the polygon when standalone. Starting from the diagrammatic expansion, the presence of localized energy fluxes in some networks (FIG. 1b) can be readily understood. Consider for instance the paths contributing to flux along an edge in a honeycomblike network. Away from the boundary, the lower order path contributions to the energy flux cancel each other, and higher order path contributions become dominant. The size of the leading order path that contributes to the energy flux along a bond increases as a function of the distance of the bond from boundary. Such a scaling results in an exponential localization of the energy flux at the boundary of the network (Supplemental Material [25]).
At the outset, given the long-ranged nature of elastic fluctuations, one might expect that the fluxes depend on topology in a highly non-local manner. Indeed, the expressions for flux in Eq. (10) suggest a complex non-local connection between the flux and the topology of the network. The results of this section show, however, that the flux can be in fact be controlled using effectively local rules. Hence, using the diagrammatic decomposition introduced here, it is easy to program energy flux patterns in networks of arbitrary complexity. Our diagrammatic expansion hence allows us to go beyond the need for de novo calculations suggested by our previous equations.

VI. ENERGY FLUX IN A PASSIVE SEGMENT COUPLED TO AN ACTIVE NETWORK
A canonical setup for the study of energy transport is a passive material bar placed between tho heat reservoirs held at constant temperature (FIG. 4a). The generic result is that energy flows from the 'hot' reservoir to the 'cold' reservoir, and in the absence of a temperature difference there can be no net energy flux through the bar.
Placing a passive material bar -in this case three masses not influenced by any magnetic fields and connected by springs in a linear geometry -across a gap in our activated metamaterial, as illustrated in FIG. 4a, reveals a very different behavior: despite the absence of temperature gradients a persistent energy flux is measured through the passive material bar. From numerical calculations, the magnitude of the energy flux decays exponentially with the bar length. In the small k limit, the flux in a bar with n sites can be evaluated using the diagrammatic techniques developed in the previous section and is equal to J n = J 1 (k/(k g + m/τ 2 )) (n−1) . Using numerical simulations, we plot the instantaneous flux transmitted across the bonds on the passive segment in  (FIG. 4b). This is reflected in the successive peaks in the instantaneous flux profile across the bonds of the passive segment. The spacing between the peaks matches the sound speed in the passive chain. This result, in combination with the results of the previous sections, shows how one design active parity violating metamaterials that can act as energy pumps and support energy transport in passive materials even in the absence of any temperature gradients.
Crucially, these results demonstrate how parity violating metamaterials can rectify non-equilibrium fluctuations. In the following section, we consider whether these rectified fluctuations when placed in contact with a vis- cous fluid, can act as low Re swimmers or fluid pumps [29][30][31].

VII. NON-RECIPROCAL MOTIONS RESPONSI-BLE FOR ENERGY FLUXES CAN BE USED TO GENERATE FORCES
The rectification of energy has been our main focus so far. In this section, we show that it is possible to exploit the energy flux to rectify motions when our model systems are allowed to interact with a viscous fluid (FIG. 5a). We begin by considering the motion of the three masses in the passive material bar discussed in the previous section and in particular, consider the effect that their motion would have on a viscous fluid. We first do so by taking their recorded trajectories and asking whether three particles following these trajectories would 'swim' in an external fluid.This calculation ignores any back action from the fluid on the dynamics of the segment. We then consider the effect of these forces in Sec. VII B and discuss regimes in which our energy conducting passive segment can generate forces when immersed in a viscous fluid. Together, these results demonstrate how a parity violating active metamaterial can be manipulated to exert forces and power motion in nanoscale materials.
FIG. 5. Utilization of non-equilibrium parity violating dynamics to power swimming and pumping in low Reynolds number (Re) media. (a) The energy flux in the active parity violating network is accompanied by nonreciprocal motions of the particles (i,ii). The nonreciprocal motion is a schematic for illustration purposes, and the real data is much more noisy. Using these nonreciprocal motions as input, a three bead linear object can be made to swim through a low Reynolds number (Re) medium (iii). Our analytical results in (b) below show that the swim speed Vs is in fact proportional to the energy flux J . (iv) Finally by immersing the passive segment into a low Reynolds fluid, the nonreciprocal motions can be used to pump the fluid. In this manner, the energy fluxes can be rectified for locomotion and force generation. In this section we consider a system of three spheres arranged in a linear configuration [31] (FIG. 5a(iii)), placed in a viscous fluid. This system is a minimal model for low Reynolds number swimming/pumping action. If the lengths of the two springs connecting the spheres, L 1 (t) = L + ∆L 1 (t), L 2 (t) = L + ∆L 2 (t) are varied according to some prescribed protocol, the time-averaged swim speed is (Eq. (12) in [31]) where a is the radius of the bead. Assumptions for this equation are a/L 1, ∆L i /L 1, and total external force on the swimmer is zero.
We now imagine recording the motions of the passive segment when it is connected to our parity violating metamaterial as in Sec. VI (FIG. 5a) and not coupled to a viscous fluid. This recorded motion can be used a protocol for modulating the configuration of an equivalent swimmer passive segment that is placed in a viscous fluid. We compute the swim speed V s of the swimmer using Eq. (14) and find that it is in fact proportional to the energy flux, J , conducted through the passive segment when it is coupled to the chiral active network (FIG. 5b). The nonreciprocal motions that is responsible for energy fluxes can also be used as a protocol to generate motion in a low Re fluid. The proportionality constant between the swim speed, V s and the energy flux, J , can be calculated using a modified diagrammatic technique (Supplemental Material [25]), where k 0 = k g +m/τ 2 (B, γ = 0 for the passive segment). This result Eq. (15) holds beyond small-k regime because all orders of paths are considered. FIG. 5b and Eq. (15) together establish that one can relate the swim speed to the flux of energy in the parity violating metamaterial. Similar proportionality between V s and J can be expected for other types of three-sphere swimmers, such as one where one sphere is much larger than the other two [32]. This is because the swim speed is generically proportional to the area enclosed in the ∆L i space [33]. This area is also proportional to the energy flux J (Supplemental Material [25]).
B. Force generation in a viscous medium using the rectified energy fluxes Finally, we now consider the scenario where the passive segment, which is connected to the active network, is immersed in a fluid (FIG. 5a(iv)). Since the passive segment is tethered due to k g and its connection to the active network, it cannot swim indefinitely. However, a stalled swimmer can potentially pump the fluid [34]. We now analyse the regime for parameters where this pumping is possible. In order to utilize the connection between the swim speed and energy flux, we require that the fluid minimally perturbs the dynamics of the passive segment. Such a regime requires the dissipation rate due to the viscous fluid to be much smaller than the energy flux through the segment. This condition can be expressed as η f av 2 J, where v is the characteristic velocity of a bead in the passive segment, η f is the dynamic viscosity of the fluid, and J is the energy flux in the absence of the fluid. In addition, the constraint of low Reynolds number requires that Re = ρ f av/η f 1, where ρ f is fluid's density. Writing these two conditions together, the requirement reads J η f av 2 ρ f a 2 v 3 . Here, we have ignored hydrodynamic interactions between the beads, because it is a higher-order perturbation with the order of a/L.
As a numerical example, this condition can be satisfied by setting k = 5 × 10 −5 kg/s 2 , k g = 1 × 10 −6 kg/s 2 for springs (values of optical traps [34]), a = 10 −6 m for all beads (size used in [34]), T a = 10 −18 J, τ = 1s for the active bath [35], ρ f = 10 3 kg/m 3 , η f = 10 −3 kg/(m · s) for liquid (water), and B = 10 −5 kg/s for the B-field. For these parameters, we obtain v = 3.8 × 10 −6 m/s, and the three relevant scales J = 5 × 10 −19 J/s, η f av 2 = 1.4 × 10 −22 J/s, a 2 v 3 ρ f = 5.4 × 10 −29 J/s, which indeed satisfy J η f av 2 ρ f a 2 v 3 . If the separation between beads is L = 10 −7 m, the swim speed calculated from Eq. (15) is 1 × 10 −8 m/s. The swim speed can provide an estimate of the potential speed at which the fluid can be pumped due to the energy fluxes. This speed can be scaled up by increasing T a . We note the extremely high value of the B field required to generate a sizable averaged flux J . To make practical use of this model, we anticipate that it will become necessary to instead consider models with interacting gyroscopes [14] or Coriolis forces [36]. Large B-field or Lorentz force analogues are easier to achieve in these cases.
In summary, the results of this section show that the energy fluxes generated by our active parity violating metamaterial can be used to rectify motion and generate forces on the nanoscale. Our calculations also show how the forces generated by our metamaterial are in fact proportional to the energy fluxes. Hence, together with the results from the previous sections that show how energy fluxes in arbitrarily complex active metamaterial networks can be controlled, our results provide a broad framework to generate and modulate forces using active parity violating metamaterials.

VIII. CONCLUSION
In conclusion, we have established a general set of design principles for rectifying energy and motion in nonequilibrium parity violating metamaterials. In particular, our central results show how a combination of time reversal symmetry violation due to the interactions and Lorentz forces in the metamaterial, and due to the nonequilibrium fluctuations of the active bath, can result in a general strategy for rectification of energy and motion. Extending these ideas to non-equilibrium parity violating metamaterials with non-linear interactions or materials composed of active chiral particles can potentially lead to new strategies for the construction of synthetic molecular motor analogues. These ideas will be considered in future work. In this section, we derive the formula for energy flux, Eq. (4) in the main text. The force F in this section is a general conservative force, which does not need to be linear. We use the following strategy to determine the energy flux. First we define the energy E i of particle i, and then write down an energy balance relation, that expresses the infinitesimal energy change dE i using stochastic calculus. Finally, we collect terms in dE i that couple neighboring particles together and identify this as the energy transfer among particles.
The energy of particle i is defined as where the first term is the kinetic energy, the second term denotes the on-site potential, and the last term is the shared spring energy between the particle and its neighbors. We use Ito's formula to calculate dE i . Ito calculus provides the advantage that the stochastic terms in dE i vanish under time-averaging. For a stochastic differential equation (SDE) of variable X(vector) with drift µ(vector) and diffusion σ(matrix) where dW is a vector consisting of standard Wiener processes, Ito's formula gives the SDE of function f (X) where ∇ X denotes the gradient with respect to X, the superscript T denotes the transpose, and tr denotes the trace. We begin by writing the equation of motion of our system Eq. (1) in the form of a stochastic differential equation Eq. (S2). We represent N particles' position by a column vector z = N i=1 |i ⊗ z i , where |i denotes the 2D subspace of particle i. Similar representations are also applied to v and η. Then we get where U is the total energy of the system, A is an antisymmetric matrix A = i |i i| ⊗ 0 1 −1 0 , and diag() means a block-diagonal matrix.

arXiv:1907.05938v1 [cond-mat.soft] 12 Jul 2019
Now we apply Ito's formula Eq. (S3) to our system by associating with the function f (X), the energy of particle i, E i (X). The nonzero terms in the gradient of E i are The term (∇ T X E i )µ reads where we used F ji = −F ij and v T i Av i = 0. The term 1 2 tr σσ T ∇ X ∇ T X f and ∇ T X f are zero. Finally, the energy change can be written as J ij is identified as the energy transferred per unit time from particle i to j, and h i is identified as the energy transferred from the bath to particle i. As for the steady-state average of J ij , we use dU ij /dt = 0 and the chain rule to simplify Eq. (S12) and arrive at the expresion (S15)

II. NUMERICAL METHOD FOR SOLVING THE ENERGY FLUX
One straightforward numerical method to compute the flux J is as follows. Our system is determined by the network geometry and parameters m, k g , k, B, γ, τ, T a . Given the equation of motion Eq. (S2),(S4)-(S6), one can numerically solve for the covariance C = XX T from the matrix equation µC + Cµ T = σσ T [S1, S2]. Finally, the flux Eq. (S12), which is bilinear in x and v, can be extracted from the covariance C. Numerical calculations of J are performed using Mathematica with custom code [S3] III. ENERGY FLUX FROM LINEAR RESPONSE THEORY Following [S4], we derive the expression of the energy flux (Eq. (7) and (10) in the main text) using a spectral linear response theory.

A. Fourier modes for energy flux
We define Fourier transform (FT) of a function f (t) as (S16) The FT of the equation of motion Eq. (S2),(S4)-(S6) reads v(ω) = iωz(ω), (S18) where G + is the response function The energy flux J ij from Eq. (S12) can be expressed as a bilinear function in z and v, by writing the linearized force F in terms of z, This bilinear form enables us to write the time integral of energy flux Q = t 0 dt J(t ) as a sum of Fourier modes q ω using Parseval's theorem, where the superscript * denotes the complex conjugate.
Pairingq ω and its conjugateq −ω gives a real function, which would be beneficial for subsequent derivations.
(S30) B. Integrating over the Fourier modes In long time limit, the sum can be approximated by an integral Eq. (S26) and (S30) can then be turned to an integral expression of the flux In the next steps, we simplify this integral with the help of the property [S4] Using this property, the trace of A q ω becomes Plugging this trace into Eq. (S32), we get the integral form for the flux Eq. (7) This integral can be calculated using the residue theorem. Since Im G + (−ω) = − Im G + (ω), Im tr G + A as 1+ω t τ 2 is an odd function of ω, and its line integral vanishes.
The integrand vanishes at ω → ∞, so the line integral can be converted to a contour integral along the counterclockwise semicircle R in the lower-half plane The noise correlation τ introduces a pole of the integrand at ω = −i/τ , thus the contour integral can be evaluated as and the response function at ω = −i/τ reads In theory, the equation Eq. (S38) provides the analytical solution of the flux, because the inverse matrix Eq. (S39) can be expressed analytically. In practice, analytical solutions can be easily calculated for small networks, but are hard for large networks. Nevertheless, some general properties of the flux can be obtained from Eq. (S38) after some algebra. For network with only horizontal and vertical bonds (FIG. 1b), all fluxes are zero. For two networks whose slanted bonds have opposite angles (FIG. 1b), their fluxes are opposite. Changing B to −B would change the flux J to −J.

IV. KIRCHOFF'S LAW
The derivation of the Kirchoff's law is similar to the derivation of the energy flux, except that we use the energy from the bath to the particle h i in Eq. (S13) instead of J ij in Eq. (S12).
Following the procedure in the last section from Eq. (S26) to (S32), we arrive at an integral expression for the h i with a different A q Using the property of G ± Eq. (S33), we get so the trace of A q ω vanishes From Eq. (S40), h i is also zero, so on average there is no energy exchange between the particle and the bath.
Because the average change of E i is zero, and Ė i = − j J ij + h i , we obtain the Kirchoff's law

CONNECTION TO ISOLATED GYROSCOPIC NETWORKS
Since our model is built upon the well-studied isolated system in Refs [S5-S8], we would like to build a connection between our energy flux in the active system and eigenmodes in those studies. In this section, we show that the flux formula Eq. (10) can be decomposed to a weighted sum over eigenmodes Eq. (S47). Then we apply this result to a honeycomb network as an example.
The Fourier analysis from Sec. IV in the main text is not suitable for this connection, because Fourier modes and eigenmodes are related only at small γ's (FIG. S1a and b), but they become dissimilar at larger γ's (FIG. S1a and  c). The underlying discrepancy between Fourier modes and eigenmodes comes from the fact that eigenmodes are a natural basis for the isolated network, whereas Fourier modes have an extra factor of friction or damping. In addition to this extra factor γ, the active system also has extra factors of m and τ . The factor m comes from the order of dynamics: the active system is second order in time, while the gyroscopic dynamics in [S5] is first order, which corresponds to the m → 0 limit.
Our starting point is Eq. (10). The key point is that, in the function G + (−i/τ ) from the equation, γ, m, τ are not independent factors. Rather, they act collectively through In effect, the extra factors m, γ, τ only add a modification to k g . Following these ideas, we imagine a reference isolated system with modified on-site spring constant k g,τ . Then after some algebra, the flux J in active system can be written as a weighted sum of the flux of each eigenmode J eig ωe in the reference system (see the next section for the derivation), Here ω e is the discrete eigen-frequency of the reference system, not to be confused with the continuous Fourier frequency ω. The amplitude of eigenmode is set such that its energy is T a , and J eig ωe is the time-averaged energy flux. A related result is a so called "sum rule", namely, the unweighted sum of all modes is zero, This "sum rule" can be derived from direct calculations (see the next section). From this eigenmode decomposition, the discussion of time-reversal symmetry in the isolated system [S5] immediately carries over to the active system. For network geometries that satisfy time-reversal symmetries, the energy flux of eigenmodes are zero. Thus through Eq. (S47), the flux in active system is also zero. This result can alternatively be obtained from Eq. (10) through some linear algebra.
As an application, we will analyze the flux in the honeycomb network using the eigenmode decomposition Eq. (S47) and the "sum rule" Eq. (S48). The flux pattern in the active honeycomb network displays CCW flux localized on the boundary (FIG. 1b). This localization is reminiscent of the edgemode in [S5] (FIG. S1b), however, their directions are opposite. From the decomposition Eq. (S47), the edgemodes should contribute a large CW flux in the active system, but somewhat surprisingly, the net flux is CCW. To better analyze the contribution from each eigenmode, we look at a simple honeycomb lattice with only one layer (FIG. S2a). This lattice has four bands (FIG. S2b), two bulk bands (blue, red) and two edge bands (green, yellow). The weighted flux of each band is plotted in FIG. S2c. We see that the CW edge band does contribute a large CW flux (green curve in FIG. S2c), however, due to the "sum rule", the unweighted sum of other bands has to be CCW. In the honeycomb lattice, it happens that many of this CCW fluxes are contained in the lower bulk band (blue curve in FIG. S2c and example mode in FIG. S2f). When the flux gets weighted, the CCW flux from lower bulk band outweighs CW flux from the edgemodes, the other two bands (yellow and red curve in FIG. S2c) also contribute to CCW flux, although relatively small. As a result, the net flux is CCW, which is opposite to the flux of the edgemode.

VI. DERIVATION OF EIGENMODE DECOMPOSITION OF ENERGY FLUX
To derive the eigenmode decomposition Eq. (S47), we first look at the reference isolated system, write down its eigenmodes Eq. (S52) and time-averaged energy flux Eq. (S56). Then we turn to the active system and decompose the flux Eq. (S38) using the eigenmodes to get Eq. (S64). Finally we show that the flux from these two sides are actually related in Eq. (S67). Lastly we also derive the "sum rule" Eq. (S69).
A. Reference isolated system The reference isolated system has first-order gyroscopic dynamics as in [S5]. In the setup with Lorentz force, the dynamical equation can be obtained by setting the mass to zero, and replacing the force matrix K by Following [S5], we convert to complex representation with z c ≡ x + iy x − iy where z 0 is the amplitude, and it will be specified shortly. The eigenmode needs a combination of ω e and −ω e to ensure that the motion of x and y is real-valued. Mathematically, this combination is possible because of a symmetry in this eigenvalue problem, when there is ω e , there is also solution −ω e with u −ωe = 0 I I 0 u * ωe .
A related property we need later is that, the left eigenvector v ωe can be expressed as v ωe = c ωe −I 0 0 I u ωe , where c ωe is a real prefactor to ensure normalization v T ωe u ωe = 1. If there are degenerate eigenvectors (like v 1 ωe , v 2 ωe , . . . ), we choose an orthonormal basis set, i.e. v i,T ωe u j ωe = 0 for i = j. With the introduction of c ωe , we now set the amplitude z 0 to z 2 0 = −2c ωe T a /ω e B, such that the energy of the eigenmode is T a . The instantaneous energy flux J ωe of mode z c ωe writes From the expression of mode Eq. (S52), When averaging over time, terms like e ±2iωet vanish, so we get Plugging in z 2 0 = −2c ωe T a /ω e B, the time-averaged flux of the eigenmode J eig ωe reads Now we turn to the active system, and the starting point is Eq. (S38). We need to decompose G τ ≡ G + (−i/τ ) into modes as below, The averaged flux J reads The second term can be shown to be zero, tr OAA as O −1 (u ωe v T ωe − u −ωe v T −ωe ) = 0. So the mode decomposition in its preliminary form reads C. The relationship between isolated system and active system Now we need to find relationship between these two fluxes J ωe Eq. (S64) and J eig ωe Eq. (S56). We will write J eig ωe in a form that looks similar to J ωe . Converting A J to A as using A as = −(A J − A T J ), and u ωe to v ωe using v ωe = c ωe −I 0 0 I u ωe , we get From direct calculation, AO −1,T = i 2 OA, and J eig ωe becomes the same as J ωe apart from a factor Comparing Eq. (S66) with (S64), we derived the relationship between flux from active system and isolated system Lastly, we show that the unweighted sum of J eig ωe is zero. This unweighted sum reads where U is the collection of all right eigenvectors U = u ωe,1 u ωe,2 · · · , and likewise for V . Since U V T = I from orthonormality, this sum vanishes

VII. PATH SUMMATION OF ENERGY FLUX
To derive the path summation formula Eq. (11) and the path rules, we start from Eq. (S38), expand around a noninteracting reference system with k = 0 to get Eq. (S74), discuss the convergence radius in Eq. (S76), then insert resolution of identity to make each term representable by a path as in Eq. (S79), and arrive at the path summation formula in (S80). We also provide a convenient way to calculate S −l in Eq. (S81), and a heuristic interpretation of S l in Eq. (S82).
We now turn on k. We denote the inter-particle part of the force matrix K as kK s , where the factor k is extracted so that the matrix K s is dimensionless. The blocks of K s read Then G τ reads In small k/k 0 regime, this matrix inversion can be expanded as To find the convergence radius, we can write the eigen-decomposition of the matrix (−K s )(k 0 G τ 0 ) as (−K s )(k 0 G τ 0 ) = W ΛW −1 , where Λ is the diagonal matrix that contains all eigenvalues λ i 's, then the flux becomes For terms in the series to be convergent, k k0 should satisfy Before inserting resolution of identity to make paths, we note that the matrix A as and K s have common blocks, A as = − |i j| ⊗ e ij e T ji + |j i| ⊗ e ji e T ij and i|K s |j = e ij e T ji , so that A as can merge with the series of G τ . (S77) Now we use the expansion Eq. (S74), and look at the contribution of its (n − 1)'th-order term to the first term of the flux, k tr i| 1 If n − 1 = 0, this term vanishes, so we only need to consider n − 1 ≥ 1 case. Insert n − 1 resolution of identity I = N la=1 |l a l a |, and plug in k 0 G τ 0 Eq. (S71) and K s Eq. (S72), we get where (−K s ) l b la ≡ l b |−K s |l a . We will denote path l = i → j → l 1 → l 2 → · · · → l n−2 → i, and its corresponding term in the above summation as S l The second term of the flux in (S77) can be treated similarly, and it results in S −l , where −l means path l in its reversed order. Combining Eq. (S79) and (S77), we get the path summation formula of the flux

B. Path rules and discussions
The path rules can be extracted from the expression of S l and J path . From the element (−K) l b la in S l , we see that either l a , l b are bonded, or l a = l b , otherwise (−K) l b la = 0. So the path has to be a closed walk along the edges of the network. From J path l for flux from i to j, we see that if the path contains equal numbers of i → j and j → i, the net contribution is zero. Because, either l = −l, so J path l ∝ S l − S −l = 0, or l ≡ −l is another path, and J path l + J path l = 0. To calculate S −l , there is a convenient way given that S l is known. Based on the transformation below, S −l can be obtained by taking the result of S l , then replacing α by −α.
To interpret S l in a more heuristic way, we insert I = e ij e T ij + e ij,⊥ e T ij,⊥ to the trace in Eq. (S79), where e ij,⊥ denotes the unit direction perpendicular to e ij . Because (−K s ) ji e ij,⊥ = 0, the trace reduces to a matrix product This expression means the following operations: starting from a unit displacement of i along e ij , j would be displaced according to the force (−K s ) ji e ij , after which j is rotated by angle α; then start from j and perform similar operations for (−K s ) l1j and R α ; finally, the transmission goes back to i; we project the displacement onto e ij , and this value is S l (apart from the prefactor ( k k0 ) n ).

C. Flux of polygon paths
Here we write down the flux formula for a polygon path without loops. It is easier to work in local coordinates, where each node has its own coordinate system. For node i in the path, let the outer angle from i to i − 1 be π, and the angle from i to i + 1 be θ i . Then the matrix (−K s ) i+1,i reads The trace in S l becomes S l /( k k 0 ) n = tr In some situations, the contribution of polygonal paths vanish, and higher order paths with loops become dominant. Unlike in polygonal paths, paths with loops are affected by side-chains.
One situation is when the polygon path itself vanishes. In FIG. S3a, the flux of lowest order path, square (FIG. S3b), is zero, so the main contribution comes from the path with length 5 (FIG. S3b). Through the loop in this path, the orientation of the side-chain controls the flux direction in the main square, without changing the geometry of the main cycle (as seen in FIG. 1b).
Another situation is that two polygon paths cancel each other, which happens in honeycomb-like networks away from the boundary (FIG. S4a). With careful calculations, the fluxes for n l ≥ 1 are not zero, they rather appear as an exponential decay (FIG. S4b). By changing the geometric angle θ, the decay length varies non-monotonically, and has a cusp at θ = α (FIG. S4c). This decay and its relationship with θ can be explained by considering the paths. While the hexagon path constitutes the lowest-order path at the boundary, it vanishes for n l ≥ 1 due to cancellations. Parameters used for numerical calculations are θ = π/6, k/k0 = 0.01, α = π/4. (c) Decay length d l changes with the network angle θ non-monotonically, and the curve has a cusp at θ = α = π/4. At small k/k0, perturbation theory results agree with numerical calculations. (d) The first non-vanishing path pair for n l = 1 has length 7. The two paths do not cancel, because the loop in the bulk and at the boundary have different values.
The first non-vanishing pair of paths for n l = 1 is shown in FIG. S4d, in which the loop exploits the asymmetry between the bulk side (with a vertical bond at the blue loop) and the boundary side (with no vertical bonds at the red loop). For every increment of one layer, the length of paths increases by 4. So the flux at layer n l is on the order of k 4n l +3 , which exhibits an exponential decay e −(n l −1)/d l . Through the calculation of these paths, we get the decay length d l = −1/ log[4(k/k 0 ) 4 (sin(θ + α) sin(θ − α)) 2 ]. From this result, we see that the cusp at θ = α in FIG. S4c is due to the term sin(θ − α). In fact, at the special point θ = α, paths like FIG. S4d vanish, and we need to consider even higher-order paths. First we rewrite V s as the following where we have defined J s ij ≡ (x i − x j )(v i + v j ) . J s ij is proportional to the energy flux via J ij = kij 2 J s ij , where k 12 = k 23 = k, and k 31 = 0 (because J 31 = 0, there is no energy flux from 3 to 1). We see that both J s 12 and J s 23 are proportional to the flux J apart from a factor k, so the remaining task is to find the relationship between J s 31 and J or J s 12 . We use a diagrammatic technique with the modification that the paths should contain only one 3 → 1 segment. This modification is a consequence of the fact that particle 3 and 1 are not bonded. We now illustrate the correspondence between the paths for J s 31 and for J s 12 . For each path l for J s 31 , we can construct n paths for J s 12 by reversing l then replacing 1 → 3 by 1 → 2(→ 2) n → 3, where n = 0, 1, . . . . An example construction of paths is shown in FIG. S5. For J s 12 , all its paths can be constructed in this way. As a result, there is a 1 to n correspondence between the paths for J s 31 and for J s 12 , which leads to the relationship where k 0 = k g + m/τ 2 (B, γ = 0 for the passive part), and the factor −2 k k0 comes from the loop 2 → 2. Plugging Eq. (S98) to the expression of V s Eq. (S97), we obtain the proportionality which is Eq. (15) in the main text. Since we have considered all the paths, this result can be analytically continued to arbitrarily large k.
From this diagrammatic technique we also see that, the proportionality constant is independent of the geometry of the active part of the network. This is because the paths through the active part for J s 31 and for J s 12 are identical.