Bounds on current fluctuations in periodically driven systems

Small nonequilibrium systems in contact with a heat bath can be analyzed with the framework of stochastic thermodynamics. In such systems, fluctuations, which are not negligible, follow universal relations such as the fluctuation theorem. More recently, it has been found that, for nonequilibrium stationary states, the full spectrum of fluctuations of any thermodynamic current is bounded by the average rate of entropy production and the average current. However, this bound does not apply to periodically driven systems, such as heat engines driven by periodic variation of the temperature and artificial molecular pumps driven by an external protocol. We obtain a universal bound on current fluctuations for periodically driven systems. This bound is a generalization of the known bound for stationary states. In general, the average rate that bounds fluctuations in periodically driven systems is different from the rate of entropy production. We also obtain a local bound on fluctuations that leads to a trade-off relation between speed and precision in periodically driven systems, which constitutes a generalization to periodically driven systems of the so called thermodynamic uncertainty relation. From a technical perspective, our results are obtained with the use of a recently developed theory for 2.5 large deviations for Markov jump processes with time-periodic transition rates.


Introduction
Thermodynamics [1] is a major branch of physics concerned with the limits of operation of machines that transform heat into other forms of energy. This theory is limited to macroscopic systems such as a steam engine. However, the way heat and temperature relate to other forms of energy is also important for small nonequilibrium systems, such as molecular motors and colloidal heat engines. For such systems, thermal fluctuations are relatively large and they cannot be ignored.
Stochastic thermodynamics [2] generalizes thermodynamics to small nonequilibrium systems. A major question that arises within this theoretical framework that takes fluctuations into account is: What are the universal relations that rule fluctuations in small nonequilibrium systems? The fluctuation theorem is one such relation [3][4][5][6][7][8], it is a constraint on the probability distribution of entropy that generalizes the second law of thermodynamics.
A more recent universal relation associated with such fluctuations is the thermodynamic uncertainty relation from [9]. This relation establishes that precision of a thermodynamic current, such as the number of consumed ATP or the displacement of a molecular motor, has a minimal universal energetic cost. Possible applications of the thermodynamic uncertainty relation include the inference of enzymatic schemes in single molecule experiments [10], a bound on the efficiency of molecular motors that depends only on fluctuations of the displacement of the motor [11], a universal relation between power and efficiency for heat engines in a stationary state [12], and design principles in nonequilibrium self-assembly [13].
The thermodynamic uncertainty relation is a consequence of a more general bound on the full spectrum of current fluctuations [14,15]. Using large deviation theory [16][17][18][19], this bound is expressed as a parabola that is above the so called rate function, which quantifies the rate of exponentially rare events. A key feature of this parabolic bound is that it depends solely on the average entropy production, i.e., knowledge of the average entropy production implies a bound on arbitrary fluctuations of any thermodynamic current. There has been much recent work related to this universal principle about current fluctuations [20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37].
The parabolic bound applies to stationary states of Markov processes with timeindependent transition rates. Physically, this situation corresponds to systems that are driven by fixed thermodynamic forces, e.g., molecular motors driven by the free energy of ATP hydrolisis. Another major class of thermodynamic systems away from equilibrium is that of periodically driven systems, which can be described as Markov processes with time-periodic transition rates. Two experimental realizations of periodically driven systems are Brownian heat engines [38] and artificial molecular pumps [39].
There is a fundamental difference with respect to fluctuations between systems driven by a fixed thermodynamic force and periodically driven systems. As shown in [40], for a periodically driven system, the energetic cost of precision of a thermodynamic current can be arbitrarily small, in stark contrast to systems driven by a fixed thermodynamic force, for which this precision has a minimal universal cost, as determined by the thermodynamic uncertainty relation. Hence, the parabolic bound from [14,15] that depends on the average rate of entropy production does not apply to periodically driven systems. For the particular case of a time-symmetric protocol, a derivation of a thermodynamic uncertainty relation has been proposed in [29]. The relation between these two classes of nonequilibrium systems is also relevant for the mapping of artificial molecular machines, which are often driven by an external periodic protocol (see [41] for a counter-example), onto biological molecular motors, which are autonomous machines driven by ATP, as discussed in [42,43].
In this paper, we obtain a universal bound on current fluctuations in periodically driven systems that is also parabolic. For the particular case of a current with increments that do not depend on time, such as internal net motion in a molecular pump, our bound depends on a single average rate. However, this average rate is different from the entropy production. For a constant protocol that leads to time-independent transition rates, our bound becomes an even more general bound than the known bound for stationary states from [14,15]. A relevant technical aspect of our proof is as follows. The parabolic bound for stationary states has been proved in [15]. This proof uses a remarkable result for large deviations in Markov processes, i.e., the exact form of the rate function for 2.5 large deviations for stationary states [44][45][46][47]. More recently, the rate function of 2.5 large deviations for time-periodic transition rates has been obtained in [48]. We use this result to prove our bounds.
Similar to the parabolic bound for stationary states that implies the thermodynamic uncertainty relation, our global bound on large deviations leads to a trade-off relation between speed and precision in periodically driven systems. We obtain a tighter local bound on the rate function that leads to an improved trade-off relation between speed and precision. For the case of stationary states, this bound is also tighter then the bound determined by the thermodynamic uncertainty relation.
We also prove our results for the case of a cyclic stochastic protocol [40,49,50]. Such protocols are convenient to perform illustrative calculations with specific models. Furthermore, the proofs for stochastic protocols are a generalization of our results for deterministic protocols, since current fluctuations for a stochastic protocol with an infinitely large number of jumps are equivalent to current fluctuations for a deterministic protocol [50].
The paper is organized in the following way. In Sec. 2 we define the basic mathematical quantities and physical models. In Sec. 3, we introduce and illustrate our main results for the case of currents with time-independent increments. The bounds are derived in Sec. 4. We conclude in Sec. 5. Appendix A contains the proofs for the case of a stochastic protocol.

Markov processes with time-periodic transition rates and fluctuating observables
We consider a Markov jump process with finite number of states Ω. The space of states is written as {1, 2, . . . , Ω}. The transition rate from state i to state j at time t is denoted by w i j (t).
Since we are interested in periodically driven systems, these transition rates have a period τ, i.e., w i j (t) = w i j (t + τ). Furthermore, we assume that if w i j (t) 0 then w ji (t) 0.
The master equation that governs the time-evolution of P i (t), the probability to be in state i at time t, reads d dt In the long time limit, P i (t) tends to an invariant time-periodic distribution π i (t) = π i (t + τ). An important quantity in this paper is the average elementary current Fluctuations can be analyzed if we consider stochastic variables that are defined as functionals of a stochastic trajectory (a t ) 0≤t≤mτ , where mτ is the final time and m is an integer number. This trajectory is a sequence of jumps and waiting times. If a jump takes place at time t, the state of the system before and after the jump is denoted by a − t and a + t , respectively. Two basic fluctuating quantities are and where dt is an infinitesimal time-interval and t ∈ [0, τ]. The empirical density ρ i j (t) counts the number of jumps per period from i to j at time t. Even though both quantities are functionals of the stochastic trajectory, to simplify notation, we do not keep the explicit dependence on (a t ) 0≤t≤mτ . The fluctuating empirical current from state i to state j is given by The average in Eq. (2) is where the brackets denote an average over stochastic trajectories. A generic current J (m) α is defined by its periodic increments α i j (t), which are antisymmetric, i.e., α i j (t) = −α ji (t), as where i< j represents a sum over all pairs of states (i, j) with i < j and with non-zero transition rates. The current in Eq. (7) can also be written in the form In stochastic thermodynamics, physical observables such as heat fluxes and particle fluxes are expressed as currents J Furthermore, the diffusion coefficient associated with J (m) α is defined as An important current in stochastic thermodynamics is the entropy increase of the environment [2], which corresponds to the increments α i j (t) = ln w i j (t) w ji (t) . The average rate of entropy production is then given by The second equality follows from π i (t) = π i (t + τ) and from Eq. (1), which leads to ∂ t π i + j i J i j = 0.

Large deviations
The rate function from large deviation theory quantifies exponentially rare events in the long time limit [16][17][18][19]. It is defined through the relation where the symbol ∼ means asymptotic equality in the limit m → ∞ and J (m) α lies in an infinitesimal interval around x. Our main result is a parabola that bounds I α (x), which is a convex function, from above. This parabola depends on an average rate. For the known parabolic bound for stationary states from [14,15], this rate is the average rate of entropy production σ in Eq. (11). In our bound for periodically driven systems, this rate is, in general, different from σ.
Current fluctuations can also be characterized by the scaled cumulant generating function where z is a real number. The cumulants associated with J (m) α can be obtained as derivatives of λ α (z) at z = 0. The scaled current generating function λ α (z) is a Legendre-Fenchel transform of the rate function I α (x), i.e.
If a parabola bounds I α (x) from above then a corresponding parabola, which can be determined from Eq. (14), bounds λ α (z) from below. For illustrations of our results we perform calculations of λ α (z) using known methods [40,50].

Stochastic protocol
We also consider the case of an external protocol that is stochastic [40,49,50]. In order to mimick a deterministic periodic protocol, this stochastic protocol is cyclic and has N states. The transition rate from state i to state j with the external protocol in state n = 0, 1, . . ., N − 1 is denoted by w n i j . The transition rate for a change in the external protocol from state n to state n + 1 mod N is γ, whereas the transition rate for the reversed transition is 0. Consider a deterministic periodic protocol characterized by the rates w i j (t) and the period τ. If the rates of the model with a stochastic protocol are w n i j = w i j (t = nτ/N) and γ = N/τ, then in the limit of N → ∞, current fluctuations for the stochastic protocol become equal to current fluctuations for the deterministic protocol [50]. Hence, the deterministic protocol corresponds to an asymptotic limit of a stochastic protocol.
In Appendix A, we derive bounds on current fluctuations for the case of a stochastic protocol. These derivations are similar to the derivation in Sec. 4 for a deterministic periodic protocol. An advantage of models with a stochastic protocol is that they are Markov processes with time-independent transition rates, which can simplify the exact evaluation of the scaled cumulant generating function in Eq. (13), as explained in [40]. Whereas the expressions in the main text are for the case of a deterministic protocol, the expressions for a stochastic protocol can be obtained from these expressions for a deterministic protocol by making the substitution

Case studies
2.4.1. Colloidal particle driven by a time-periodic field The first model in Fig. 1(a) is a biased random walk on a ring with Ω states driven by a time-periodic force F(t) ≡ F 0 cos(2πt/τ). A physical realization of this model is a charged colloid on a ring subjected to a time-periodic electrical field. We set Boltzmann constant k B and the temperature T to k B T = 1 throughout. The transition rate for a jump in the clockwise direction is k + (t) ≡ ke F(t)/Ω and the reversed transition rate is k − (t) ≡ k. These transition rates satisfy the generalized detailed balance relation [2]. The current we consider is the net number of jumps in the clockwise direction per unit time. For this model, the scaled cumulant generating function in Eq. (13) can be calculated exactly [50].

Molecular pump
The other two models are driven by a stochastic protocol. The model illustrated in Fig. 1(b) is a molecular pump with Ω = 3. This model has been introduced in [40]. The external protocol changes energies and energy barriers between states, which can lead to net rotation in the ring with three states. The number of states of the external protocol is N = 3. The states of the external protocol are denoted by 0, 1, 2, which correspond respectively to the top left circle, the top right circle and the bottom circle in Fig. 1(b). In this model, the energies and energy barriers are rotated in the clockwise direction by one step if a jump (with rate γ ) that changes the state of the protocol takes place. The energies are denoted by E 1 , E 2 , and E 3 , whereas the energy barriers are denoted by B 1 , B 2 , and B 3 . The internal transition rates are given by for j = i + 1, and for j = i − 1, where we assume periodic boundary conditions. An important property of molecular pumps is that the thermodynamic force is zero for any state n of the external protocol. This physical condition is manifested in the following restriction on the transition rates The current we consider is the net number of jumps in the clockwise direction per unit time.
The scaled cumulant generating function in Eq. (13) associated with this current can be calculated from the eigenvalue of a modified generator, as shown in [40].

Enzymatic reaction with stochastic substrate concentrations
The model illustrated in Fig. 1(c) is a model with Ω = 4 and two independent thermodynamic forces F n 1 and F n 2 , which depend on the state of the external protocol n. This model can be interpreted as a enzyme that can consume two different substrates and produces one product [9]. The two enzymatic cycles are E + S 1 → E S 1 → E P → E + P and E + S 2 → E S 2 → E P → E + P, where E is the enzyme, P is the product, S 1 is one substrate, and S 2 is another substrate. State 1 corresponds to the free enzyme E, state 2 corresponds to E S 1 , state 3 corresponds to E S 2 , and state 4 corresponds to E P. The external control of the concentrations of the substrates S 1 and S 2 generate thermodynamic forces that depend on n. The number of states of the external protocol is N = 2. The generalized detailed balance relation for this model reads The thermodynamic forces change between two values of the same modulus and different sign stochastically, i.e., F n 1 is given by F 0 1 = F 1 and F 1 1 = −F 1 , whereas F n 2 is given by F 0 2 = F 2 and The transition rate for a change of the external protocol is γ. The transitions rates are set to w n 12 = ke F n 1 /2 , w n 13 = ke F n 2 /2 , w n 14 = k, w n 21 = k, w n 24 = ke F n 1 /2 , w n 31 = k, w n 34 = ke F n 2 /2 , w n 41 = k, w n 42 = k, and w n 43 = k. The current we consider is the elementary current from state 1 to state 2, which corresponds to the net number of S 1 molecules that have been consumed per unit time. As is the case of the previous model, the scaled cumulant generating function in Eq. (13) can be calculated with the method explained in [40]. The parameters for the model represented in Fig 1(a) are set to F 0 /Ω = 2 and k = τ = 1. The parameters for the model represented in Fig 1(b) are set to E 1 = E 3 = B 1 = B 2 = 0, E 2 = 2, B 3 = 5, and γ = 1/10. The parameters for the model represented in Fig. 1(c) are set to F 1 = 2, F 2 = 1/2, k = 1, and γ = 1/10. (b) The ratio R α ≡ J 2 α /(2D α σ * ) ≤ 1/2 as a function of the rate γ for jump of the protocol. We have analyzed the model illustrated in Fig. 1(b) with parameters E 1 = 1, B 1 = 5, k = 1, and E 2 = E 3 = B 2 = B 3 = 0, and the model illustrated in Fig. 1 Fig.  1(c) with parameters k = γ = 1 and two values of F 2 .

Main results
In this Section we discuss our main results for currents with time-independent increments α i j (t) = α i j , which include the case of currents generated in a molecular pump. For time-independent increments, the results acquire a simpler form with a more direct physical interpretation. In Sec. 4, we present proofs of more general results, which, inter alia, also hold for currents with time-dependent increments. Physical examples of currents with timedependent increments include the heat and work currents in heat engines. The general features of our main results presented in this Section are the same irrespective of whether the protocol is deterministic or stochastic, which is discussed in Appendix A.

Global bound
The parabolic bound on the rate function is where The inequality σ * ≥ 0 comes from the fact that for fixed t every term in the sum i< j in Eq. (20) is not negative. In general, the average rate σ * is different from the thermodynamic rate of entropy production σ in Eq. (11). Furthermore, there is no simple inequality relating both quantities.
For the case of time-independent transition rates w i j (t) = w i j , σ * = σ and the bound (19) becomes This bound is the known parabolic bound for time-independent transition rates proved in [15]. Hence, Eq. (19) constitutes a generalization of this parabolic bound to periodically driven systems.
In terms of the scaled cumulant generating function, the bound in Eq. (19) is written as where we used Eq. (14). The universality of our result is illustrated in Fig. 2(a), where we compare the functioñ for the models in Fig. 1, with the lower boundz(1 +z). This bound, or the bound in Eq. (19), is a particular case of two bounds, one derived in Sec. 4.1 and the other derived in Sec. 4.5.

Trade-off between speed and precision
Taking the second derivative of I α (x) at x = J α , we obtain the diffusion coefficient D α defined in Eq. (10) as The inequality in Eq. (19) and the fact that this inequality is saturated at x = J α , leads to the following bound on D α , In Sec. 4.3, we derive a local quadratic bound on I α (x), which is valid for x close to the average J α . This local bound together with Eq. (25), gives a tighter bound on D α that reads whereσ The second inequality in Eq. (27) is a consequence of σ * ≥σ, which follows from the inequality where a and b are positive.
Rearranging the terms in Eq. (27), we write the following universal trade-off relation between speed and precision for periodically driven systems, where F α ≡ 2D α /J α is the Fano factor. The Fano factor characterizes the precision associated with J (m) α , whereas J α quantifies the speed. In periodically driven systems, a current with small fluctuations, as characterized by a small Fano factor F α , can only be as fast asσF α /2.
This trade-off relation is a generalization of the thermodynamic uncertainty relation to periodically driven systems. In particular, for the case of time-independent transition rates w i j (t) = w i j , inequality (30) implies the thermodynamic uncertainty relation F −1 α J α ≤ σ/2, since σ * = σ for this case. Furthermore, the inequality F −1 α J α ≤σ/2, for time-independent transition rates, provides an even tighter bound than the thermodynamic uncertainty relation.
This result is relevant for the mapping between an artificial molecular pump and a system driven by a fixed thermodynamic force such as a biological molecular motor, which is modeled with time-independent transition rates that lead to a nonequilibrium stationary state, proposed in [42]. With this mapping, one can construct a molecular pump that mimicks a stationary state and vice-versa, in the sense that both the average rate of entropy production and the average elementary currents between a pair of states are conserved. However, a mapping of a molecular pump onto a stationary state that also preserves fluctuations is not always possible, since a molecular pump may not fulfill the relation J 2 α /(2D α σ) ≤ 1/2, as shown in [40], whereas a system that reaches a nonequilibrium stationary state must fulfill this relation.

Discussion of the bounds
In Fig. 2(b), we show plots of R α ≡ J 2 α /(2D α σ * ) ≤ 1/2 as a function of the rate γ, which quantifies the speed of the protocol, for the models illustrated in Fig. 1(b) and in Fig. 1(c). For the first model, which is a molecular pump, we find that this bound is saturated if the transitions of the protocol are much slower than the internal transition rates associated with changes of the state of the system. For this model, in this limit the bound is saturated independent of the values of the energies and energy barriers. However, for the second model the bound is not saturated in this limit.
In Fig. 2(c), we show plots of R α ≡ J 2 α /(2D α σ * ) ≤ 1/2 for the model illustrated in Fig.  1(c). The quantity in the horizontal axis is the thermodynamic force F 1 . For this model, the bound is saturated for F 1 small and the other thermodynamic force F 2 = 0. This saturation of the bound is similar to the saturation of the bound for steady states known as thermodynamic uncertainty relation, which happens in the linear response regime [9].
Finally, let us comment on the rate σ * that we have introduced here. Its physical interpretation is that σ * , and not the rate of entropy production σ, provides a bound on the whole spectrum of fluctuations for any current (with time-independent increments) in a generic periodically driven system arbitrarily far from equilibrium. In terms of the trade-off relation from Eq. (30), σ * (and alsoσ) provides a limit on how precise and fast a thermodynamic current can be. The rate of entropy production σ quantifies the energetic cost of sustaining the operation of the nonequilibrium system. Interestingly, for time-independent transition rates corresponding to a system driven by a fixed thermodynamic force, σ * = σ is a rate that has both physical properties, i.e., it bounds current fluctuations and quantifies energetic cost.

General bounds
In this Section we derive the bounds that imply the results discussed in Sec. 3. We obtain two global bounds that imply the global bound in Eq. (19), the first one is given in Eq. (47) and the second one is given in Eq. (67). We also derive a local bound that leads to the inequality in Eq. (56), which generalizes the trade-off relation in Eq. (30).

First global bound
In our proof we use the theory for 2.5 large deviations for periodically driven systems developed in [48]. At the level 2.5 the joint distribution of all empirical densities defined in Eq. (3) and all empirical currents defined in Eq. (5) is considered. In our notation ρ(t) represents a vector with the empirical densities that has dimension Ω and J(t) is a vector with the empirical currents that has dimension M, where M is the number of unordered pairs of states with non-zero transition rates. The advantage of considering this level of large deviations is that the rate function can be calculated exactly as where and Note that the quantities G and a depend on the empirical density ρ. The empirical density and current in Eq. (31) fulfill the constraint d dt for all states i. To simplify the notation we write I cur 2.5 [J(t), ρ(t)] instead of the l.h.s. of Eq. (31).
The name level 2.5 large deviations can also refer to the rate function associated with the joint probability of the empirical density and the empirical flow defined in Eq. (4). The rate function with the empirical current can be obtained from the rate function with the empirical flow [48].
An important technique in large deviation theory is the so called contraction [16][17][18][19], for which the rate function associated with a coarse-graining of the number of variables can be obtained from the original rate function. Hence, the rate function for an arbitrary current J α can be obtained from a contraction of I cur 2.5 [J(t), ρ(t)], which leads to the expression where J(t) and ρ(t) are such that they fulfill Eq. (35) and the relation 1 τ In particular, this relation leads to the inequality whereG andã are functions ofρ as in (32) and (33). This inequality is valid for any pair of vectors that fulfill the constraints d dtρ for all states i, and 1 τ The inequality [15] ψ J i j , G i j , a i j ≤ 1 4 together with Eq. (38), leads to We are now left with the problem of finding a judicious choice of J (t),ρ(t) that fulfills the constraints in Eq. (39) and in Eq. (40). One such choice is The time-independent parameters K i j are antisymmetric, i.e., K i j = −K ji , and satisfy j i for all states i. Using this choice in Eq. (42), we obtain where and The global bound in Eq. (47), together with Eq. (25), leads to  (46), K i j can be seen as the current of some auxiliary Markov process with time-independent transition rates in the stationary state. A natural choice of K i j is to consider the time-integrated probability current, as defined in Eq. (21), i.e., For this choice and σ * K = σ * , where σ * is defined in Eq. (20). For currents with time-independent increments α i j (t) = α i j , we obtain i< jJi jᾱi j = J α , where J α is given by Eq. (9), and the bound in Eq. (47) becomes the bound in Eq. (19). For currents with time-dependent increments, which include the rate of extracted work and the rate of heat flow in a heat engine driven by periodic temperature variation, the rate J K in Eq. (52) is, in general, different from the average current J α .

Other possible choices for K
The freedom of choice for the parameter K depends on the network of states of the Markov process, with Eq. (46) limiting the number of independent currents K i j [51]. For instance, for the unicyclic model in Fig. 1(a), there is just one independent current and K i j is the same for all pairs of states. In this case, the ratio σ * K /J 2 K becomes independent of K and, therefore, there is only one bound in Eq. (47) regardless of the value of K i j . We note that the same argument about the freedom of choice for the parameter K applies to stochastic protocols, as is the case of the model in Fig.1(b) If we consider a model with the network of states shown in Fig. 1(c), then there are two independent K i j and different choices for these parameters can lead to different bounds in Eq. (47). Two particularly appealing choices for the parameter K are the choices that conserve the rate of entropy production or the average current in Eq. (47). The first choice corresponds to a K that fulfills the relation σ * K = σ and the second choice corresponds to a K that fulfills the relation J K = J α . Whether it is possible to set K in such a way that one of these relations is fulfilled is a question that depends on the model (or class of models) at hand.

Local bound
We now derive a local quadratic bound on I α (x) that leads to the first inequality in Eq. (27). For a and G fixed, a Taylor expansion of the function ψ(J, G, a) for J around the value G, leads to Applying this Taylor expansion to Eq. (38) withρ andJ given by (43) and (44), respectively, we obtain the local bound where J K is defined in Eq. (48) and .
The local bound in Eq. (54) together with Eq. (25) leads to A generic model-independent choice for K is the one given in Eq. (51), i.e., K i j =J i j . If, in addition, the increments are time-independent, the bound in Eq. (56) becomes the trade-off relation between speed and precision in Eq. (30). We recall that from Eq. (29), σ * K ≥σ K , thus, the bound in Eq. (56) is stronger than the bound in Eq. (50).

Bounds for time-independent transition rates
Here, we stress that the bounds for time-periodic transition rates derived above imply new bounds for the case of time-independent transition rates that lead to a non-equilibrium stationary state. For time-independent transition rates, and for currents with time-independent increments, the terms in Eq. (47) become and Hence, from Eq. (47) we have the bound For K i j = J i j , Eq. (59) becomes the known parabolic bound for stationary states from [14,15]. Furthermore, for time-independent transition rates Eq. (56) becomes This bound is tighter then the bound on the diffusion coefficient that follows from Eq. (59). For the case K i j = J i j , Eq. (60) becomes an even stronger bound than the thermodynamic uncertainty relation, as discussed in Sec. 3.

Second global bound
We can obtain a bound different from the global bound in Eq. (47) by considering a choice for J i j (t) that is different from the one in Eq. (44). We write the stationary distribution of a master equation with frozen transition rates w i j (t) as µ i (t). This quantity is known as accompanying density [52]. Due to the periodicity of w i j (t) we have µ i (t) = µ i (t + τ). We consider the bound in Eq. (42) withρ i (t) = π i (t) and where c 1 (t) and c 2 (t) are time-periodic functions, K i, j is antisymmetric and fulfill the relation in Eq. (46), and Since j i M i j (t) = 0, which comes from the definition of the accompanying density µ i (t), this choice fulfills the constraint in Eq. (39). Setting K i j =J i j , c 1 (t) = c 1 , and c 2 (t) = c 2 , the constraint in Eq. (40) applied to the choice in Eq. (62), leads to where q is an arbitrary real number and The bound in Eq. (42) then becomes Minimization over the single parameter q gives the tightest bound on the large deviation function. For q = 0 we obtain the bound in Eq. (47) with K i j =J i j . However, for q = 1 we obtain a bound that cannot be obtained from Eq. (47), which reads This bound with σ * µ that depends on the accompanying density µ i (t) seems particular appealing for the case of time-scale separation, with the transition rates much faster than the inverse of the period τ −1 .

Conclusion
The thermodynamic uncertainty relation and the parabolic bound on current fluctuations that generalizes it, constitute major recent developments in stochastic thermodynamics that are valid for Markov processes with time-independent transition rates that reach a stationary state, which describes a system driven by fixed thermodynamic forces. We have generalized these bounds to periodically driven systems. Similar to the bound for stationary states, we obtained a bound that depends on the single average rate σ * . However, for periodically driven systems this average rate is, in general, different from the thermodynamic entropy production σ. These rates have two essential physical properties: while σ quantifies the energetic cost of maintaining the system out of equilibrium, σ * provides a generic limit to current fluctuations.
The quite high degree of universality of our results are encouraging with respect to possible applications. For instance, we have found a trade-off relation between speed and precision in periodically driven systems for currents that have time-independent increments. Physically, such relation tells us that if one wants to generate net motion in a artificial molecular pump driven by an external periodic protocol, there is a universal limit on how fast and precise this net motion can be.
For the case of the thermodynamic uncertainty relation for stationary states, several applications have been proposed [10][11][12][13]. Figuring out how to extend these applications to periodically driven systems is an interesting direction for future work. One particular instance would be to extend the universal relation between power, efficiency and fluctuations from [12] to periodically driven heat engines. The more general bounds derived in Sec. 4 that apply to time-dependent increments, might be important for these applications. Finally, our results could be verified in experimental systems such as periodically driven colloidal particles and artificial molecular pumps.

Appendix A. Stochastic Protocol
Appendix A. 1

. Mathematical definitions
The master equation for the model with a stochastic protocol reads d dt P n i = j i P n j w n ji − P n i w n i j + γ(P n−1 i − P n i ), where n − 1 = N − 1 for n = 0 and P n i is the time-dependent distribution. The stationary distribution of state (i, n) is denoted by π n i . The stationary distribution of the state n of the protocol is given by π n ≡ i π n i = 1/N, which comes from the solution of the master equation (A.1) for the stationary distribution. The conditional probability for the system to be in state i given that the protocol is in state n is written as π(i|n) = π n i /π n = Nπ n i . Consider a timeperiodic Markov process with rates w i j (t) and period τ. If the transition rates fulfill the relation w n i j = w i j (t = nτ/N) and γ = N/τ, then, in the limit N → ∞, π(i|n) → π i (t) [40], where n = [tN/τ] and [·] denotes the integer part. Therefore, if we consider the average elementary current J n i j ≡ π n i w n i j − π n j w n ji in the limit of N → ∞, we obtain NJ n i j → J i j (t), where n = [tN/τ]. This relation is important for the connection between the cases of a deterministic and stochastic protocols. A stochastic trajectory is denoted by (b t ) 0≤t≤t f , where t f is the final time. Note that a state of the Markov process here is specified by the variable that determines the state of the system i and the variable that determines the state of the protocol n. The stochastic trajectory has a fluctuating number of jumps N f , the time interval between two jumps is denoted ∆t k , with k = 0, 1, . . ., N f , and the state of the Markov process during the time-interval ∆t k is denoted b k .
The empirical density of state (i, n), which is the fraction of time spent in this state, is defined as δ b k ,(i,n) is the Kronecker delta between the state of the trajectory b k and the state (i, n). The notation here in the appendix is different from the notation in the main text for the case of a deterministic protocol. If we compare Eq. (A.3) with Eq. (3), we see that here the upper index in ρ n i refers to the state of the stochastic protocol and is equivalent to t in ρ (m) i (t), for which the upper index m refers to the time interval of the stochastic trajectory. For a more compact notation we do not keep the dependence of the fluctuating quantities on the time interval t f .
The empirical current from state (i, n) to state ( j, n) reads For the case of a stochastic protocol, we also consider the empirical flow (or unidirectional current) from state (i, n) to state (i, n + 1), where n + 1 = 0 for n = N − 1, which is defined as The average of this empirical flow in the stationary state is C n i ≡ C n i = γπ n i . A generic fluctuating current is written as where α n i j = −α n ji are the increments. If we compare this expression with Eq. (7), which is the expression for a deterministic protocol, we see that an integral over a period divided by the period τ for a deterministic protocol becomes a sum over n divided by the total number of states of the protocol N for a stochastic protocol. Note that the factor 1/N does not appear in front of the sum in the r.h.s of Eq. (A.6) due to Eq. (A.2). The average current in the stationary state reads The rate function associated with J α is defined as where ∼ means asymptotic equality in the limit t f → ∞. The scaled cumulant generating function for a stochastic protocol is defined as These two quantities are related by a Legendre-Fenchel transform, as in Eq. (14). 2(K i j ) 2 π n i w n i j + π n j w n ji , (A.14) which is equivalent to Eq. (55). The parameter K i j in these equations is anti-symmetric, i.e., K i j = K ji , and thus fulfill j i K i j = 0 for all i.