Current fluctuations in periodically driven systems
Andre C Barato Raphael Chetrite





^1 Max Planck Institute for the Physics of Complex Systems,

Nöthnitzer Str. 38, 01187 Dresden, Germany

^2Laboratoire J A Dieudonné,

UMR CNRS 7351, Université de Nice Sophia Antipolis, Nice 06108, France






#1⟨#1 ⟩

Small nonequelibrium systems driven by an external periodic protocol can be described
by Markov processes with time-periodic transition rates. In general, current fluctuations 
in such small systems are large and may play a crucial role.  We develop a 
theoretical formalism to evaluate the rate of such large deviations in periodically driven systems. We show
that the scaled cumulant generating function that  characterizes current fluctuations is given 
by a maximal Floquet exponent. Comparing deterministic protocols with stochastic protocols, 
we show that, with respect to large deviations, systems driven by a stochastic protocol with an infinitely large 
number of jumps are equivalent to systems driven by deterministic protocols. 
Our results are illustrated with three case studies: a two-state model for a heat engine,
a three-state model for a molecular pump, and a biased random walk with a time-periodic affinity.     







§ INTRODUCTION



Periodic external control is used to operate a wide variety of thermodynamic machines that includes 
traditional idealized engines. Modern experimental examples of such machines are molecular pumps <cit.> 
and micro-sized heat engines <cit.>. For these small systems, thermodynamic currents, such as 
the work exerted on a molecular pump or the heat flow in a heat engine, display thermal fluctuations 
that can be relatively large. For instance, the prominent fluctuation theorem is a symmetry related 
to these fluctuations. 

Stochastic thermodynamics <cit.> is an emerging field  that applies to small systems with 
large fluctuations. Within this theory, periodically driven systems are modeled as Markov processes with time-dependent
transition rates that are periodic. Such modeling  has been used in several works that include:
models for stochastic resonance <cit.>, linear response theory for periodically driven systems <cit.>, 
theoretical studies for small heat engines far from equilibrium <cit.>, necessary conditions for 
the generation of a current in a molecular pump <cit.>, a mapping relating periodically driven systems 
with systems driven by a fixed thermodynamic affinity <cit.>, and 2.5 large deviations for Markov process with time-periodic generators <cit.>.
However, a generic formalism to evaluate the rate of large deviations of single 
currents is still not available.

Periodically driven systems reach a limiting periodic state that can be contrasted with nonequilibrium 
stationary states, which are described by Markov  processes with constant transition rates. Physically, this second 
case corresponds to a system driven by a fixed thermodynamic affinity. For nonequilibrium 
stationary states, a generic formalism to quantify large deviations of currents is available: The so called scaled cumulant generating function 
(SCGF) is determined by the maximum eigenvalue of a tilted generator <cit.>. This formalism can be used to 
calculate the SCGF of systems subjected to a potential that is periodic in space <cit.>, which is not the case for a 
time-periodic potential. 

In this paper, we develop a formalism to determine large deviations in periodically driven systems. 
We show that a fundamental (or monodromy) matrix from Floquet theory <cit.>,  which is related to a time-dependent 
tilted generator, quantifies current fluctuations.  Specifically, we show that this fundamental matrix 
is a Perron-Frobenius matrix and that its maximal eigenvalue gives the SCGF.


Deterministic protocols are commonly used in the study of periodically driven systems. Nevertheless, such systems can also be driven by a cyclic stochastic 
protocol that mimics the periodicity of a deterministic protocol <cit.>. In this case, the system and protocol together form a bipartite Markov process with time-independent 
transition rates. For stochastic protocols, current fluctuations can then be analyzed within the stationary state of this 
bipartite Markov process. We prove the equivalence of current fluctuations between systems driven by a 
deterministic protocol and systems driven by a stochastic protocol with an infinitely large number of jumps. Therefore, we show that 
a periodic protocol can be seen as a particular limit of a stochastic protocol.  

Illustrations of our results are performed with three models. An exactly solvable model 
for a heat engine, a model for molecular pump that we use to compare theory with numerical simulations, and 
one model for a biased random walk with a time-periodic affinity that has a particularly simple 
expression for the SCGF.

The paper is organized in the following way. In Section <ref>, we define fluctuating currents in Markov processes with time-periodic transition rates. 
The formalism for the calculation of the SCGF is developed in Section <ref>. In Section <ref>, we analyze
three case studies. Section <ref> contains the results for a stochastic protocol. We conclude in Section <ref>. 
In <ref>, we discuss fluctuating currents such as work, which can be written in a form that is apparently different 
from the generic currents we consider in the main text. We extend our results to diffusion processes in <ref>. 
An exact calculation of the SCGF for systems with two states and a piecewise constant protocol is presented in <ref>.    




§ GENERAL SETUP AND MATHEMATICAL DEFINITIONS




We consider a Markov process with a finite number of states 
Ω. The time-dependent transition rate from state i
to state j at time t is denoted by w_ij(t). These transition rates are 
time-periodic with a period τ, i.e., 

    w_ij(t+τ)= w_ij(t).

In the theoretical framework of stochastic thermodynamics, which we use in our illustrative examples, the transition rates 
fulfill the restriction that if w_ij(t)≠ 0 then w_ji(t)≠ 0. However, our mathematical results (in Section <ref> and Section <ref>)
do not rely on this assumption.  A physical interpretation for these transition rates based on the generalized detailed balance relation from stochastic thermodynamics 
can be found in <cit.>.

The associated master equation reads 

    d/dt P(i,t)= ∑_j[P(j,t)w_ji(t)-P(i,t)w_ij(t)],

where P(i,t) is the probability to be in state i at time t. 
This equation can be written in the vectorial form 

    d/dt|P_t⟩= ℒ_t|P_t⟩,

where |P_t⟩ is a vector with components P(i,t)  and ℒ_t is the stochastic matrix defined as 

    [ℒ_t]_ji≡ (1-δ_ij)w_ij(t) -δ_ij∑_kw_ik(t).


A stochastic trajectory from time 0 to time T=nτ is a sequence of jumps and waiting times, which is denoted by A_0^T. If a jump takes place at time t, the state before the jump 
is denoted a_t^- and the state after the jump is denoted a_t^+. A fluctuating current is a functional of the stochastic trajectory defined as 

    X[A_0^T]≡∑_0≤ t≤ Tθ_a_t^-,a_t^+(t).

For a current, the increments θ_i,j(t) are anti-symmetric, i.e., 

    θ_i,j(t)=-θ_j,i(t).

Fluctuations of this current in the long time limit are characterized by  the SCGF

    λ(z)≡lim_T→∞1/Tln⟨e^zX⟩=1/τlim_n→∞1/nln⟨e^zX⟩,

where the brackets mean an average over stochastic trajectories. 
The average current J and diffusion coefficient D are given by 

    J≡lim_T→∞1/T⟨ X⟩= λ'(0)

and

    D ≡lim_T→∞1/T(⟨ X^2⟩-⟨ X⟩^2)= λ”(0),

respectively. Similarly, higher order moments associated with X can be obtained by taking higher order derivatives
of λ(z) at z=0.  

In the long time limit, the system reaches an invariant limiting periodic distribution P^inv_i(t)=P^inv_i(t+τ).
The average current can be written in terms of this distribution as

    J=  1/τ∫_0^τ∑_i<jθ_ij(t)[P^inv_i(t)w_ij(t)-P^inv_j(t)w_ji(t)]dt,

where ∑_i<j means a sum over all links with non-zero transition rates in the network of states. 

An important observable in stochastic thermodynamics is the work done on the system due to the periodic variation of 
energy of the system. This fluctuating current is typically written as an integral of a function over the time interval T. However, 
as we show in appendix <ref>, observables such as work can also be written in the form given by Eq. (<ref>).


 




§ FLOQUET THEORY FOR THE SCGF






 §.§ General theory


The joint probability that the current is X and the system is in state i at time t is written as P(i,X,t), whereas the vector 
|P_t(X)⟩ has components P(i,X,t). The Laplace transform of |P_t(X)⟩ is given by 

    |P̃_t(z)⟩≡∑_X e^Xz|P_t(X)⟩.

The average ⟨e^zX⟩ in Eq. (<ref>) is related to this Laplace transform in the following way,  

    ⟨e^zX⟩= ∑_i=1^ΩP̃(i,z,T).

From the master equation (<ref>) and the tilted generator 

    [ℒ_t(z)]_ji≡ (1-δ_ij)w_ij(t)e^zθ_ijt -δ_ij∑_kw_ik(t),

we obtain

    d/dt|P̃_t(z)⟩=ℒ_t(z)|P̃_t(z)⟩.

Using the periodicity of ℒ_t(z), the formal solution of this equation at time T=nτ is

    |P̃_T(z)⟩=exp(∫_0^Tℒ_t(z)dt) |P̃_0(z)⟩= [exp(∫_0^τℒ_t(z)dt)]^n|P̃_0(z)⟩≡ℳ(z)^n|P̃_0(z)⟩,

where |P̃_0(z)⟩ is the initial condition and exp represents a time-reversed ordered exponential. 
This ordered exponential can be defined as the solution of the differential equation  

    d/dtexp(∫_0^tℒ_t'(z)dt')= ℒ_t(z)exp(∫_0^tℒ_t'(z)dt'),

where the initial condition is the identity matrix. 

The matrix 

    ℳ(z)≡exp(∫_0^τℒ_t'(z)dt')

is a central object that is known as fundamental matrix in Floquet theory <cit.>. 
The eigenvalues of this matrix are denoted by ρ_k(z), the right  eigenvectors by |r_k(z)⟩, and 
the left eigenvectors by ⟨l_k(z)|. The fundamental matrix can then be written as 

    ℳ(z)=   ∑_k=1^Ωρ_k(z)|r_k(z)⟩⟨l_k(z)|.

From Eq. (<ref>), by Setting T=τ, imposing the initial condition a_0=i, and 
restricting to trajectories that finish at state a_T=j, we obtain  

    ⟨e^zXδ_a_T,j|a_0=i⟩= [ℳ(z)]_ji.

This equation shows that all elements of the fundamental matrix ℳ(z) are positive. Hence, from the Perron-Frobenius theorem, the matrix ℳ(z) has 
a maximal real eigenvalue defined as  ρ_1(z). This fact together with the definition of the SCGF
in Eq. (<ref>), Eq. (<ref>), Eq.  (<ref>),  and Eq.  (<ref>) 
lead to the main result 

    λ(z)≃1/τlnρ_1(z),

where the symbol ≃ means asymptotic equality in the limit n→∞ and τ^-1lnρ_1(z) is the maximal Floquet exponent. 
The SCGF λ(z) can be evaluated by first 
solving Eq. (<ref>) and then calculating the maximal eigenvalue of ℳ(z). 
For z=0 the matrix elements in Eq. (<ref>) are transition probabilities, therefore, 
the maximum eigenvalue of the matrix is  ρ_1(0)=1, which implies λ(0)=0. It is worth noting that 
Eq. (<ref>) is not restricted to currents. This result is also valid 
for any functional of the stochastic trajectory with the form given in Eq. (<ref>) that has increments that do not
fulfill the anti-symmetry in Eq. (<ref>). We point out that 
the SCGF has been obtained as a maximum Floquet exponent for specific two-state models in <cit.>.

This relation between SCGF and maximal Floquet exponent is also valid for diffusion processes, as shown in <ref>. 
For piecewise protocols, transition rates w_ij(t) are piecewise. 
If the external protocol is piecewise constant, with the period 
divided into L pieces and τ=τ_0+τ_1+…+τ_L-1, the matrix  ℳ(z) 
defined in Eq. (<ref>) takes the form 

    ℳ(z)= exp(ℒ̃_L-1(z)τ_L-1)…exp(ℒ̃_1(z)τ_1)exp(ℒ̃_0(z)τ_0),

where ℒ̃_k(z) is the constant modified generator during the interval τ_k. In this last equation, it is assumed that 
the increments of the current are also piecewise constant. The SCGF
can then be obtained from the maximal eigenvalue of this matrix. For a piecewise protocol that is 
not constant, the expression of ℳ(z) becomes a product of ordered exponentials. 




 §.§ Expressions for average current and diffusion coefficient
  

The average current J, diffusion coefficient D, and higher order cumulants can be obtained without explicit 
evaluation of the maximum eigenvalue associated with ℳ(z) in the following way. A similar method for non-equilibrium stationary states 
has been introduced by Koza <cit.> (see also <cit.>). The characteristic polynomial associated with ℳ(z) is written as

    (x𝐈-ℳ(z))=∑_m=0^Ωc_m(z)x^m,

where 𝐈 is the identity matrix.
The maximum eigenvalue ρ_1(z) is a root of this polynomial, which leads to the equation

    ∑_m=0^Ωc_m(z)[ρ_1(z)]^m=0.

Taking a derivative with respect to z and setting z=0 we obtain 

    J= 1/τρ'_1(0)/ρ_1(0)=-1/τ∑_m=0^Ωc'_m(0)/∑_m=0^Ωmc_m(0),

where we used Eq. (<ref>) and ρ_1(0)=1. Taking a 
second derivative with respect to z of Eq. (<ref>) and setting z=0 we obtain  

    D     = 1/τ{ρ”_1(0)-[ρ'_1(0)]^2}
        =-∑_m=0^Ωc”_m(0)+2ρ'_1(0)∑_m=0^Ωmc'_m(0)+[ρ'_1(0)]^2∑_m=0^Ωm^2c_m(0)/τ∑_m=0^Ωmc_m(0),

where we used Eq. (<ref>) and ρ_1(0)=1. Using these expressions,  J and D can be evaluated directly from the coefficients of the characteristic polynomial associated with 
ℳ(z). Taking higher order derivatives lead to similar expressions for higher order cumulants.   

























§ CASE STUDIES





























































 §.§ Heat engine



We introduce an exactly solvable two-state model for a heat engine with a piecewise protocol, which is illustrated in Fig. <ref>. This model is similar to a model for a heat 
engine with a stochastic protocol analyzed in <cit.>. One of the states has energy 0 and the other state 
has a higher energy that depends on time. The protocol is divided in four steps. First, the energy changes from E to E+Δ E at a cold temperature β_c^-1.
Second, the temperature changes from β_c^-1 to the hot temperature β_h^-1. Third, the energy changes back from E+Δ E to E. Fourth, the temperature 
changes back from β_h^-1 to β_c^-1. The inverse temperature takes the form 

    β^k= β_c[1-F_qh^k],

where F_q≡(β_c-β_h)/β_c≥ 0, h^0=h^1=0 and h^2=h^3=1. Physically, F_q is the thermodynamic affinity associated with the heat current <cit.>
and the temperature depends on time through the parameter h^k. The energy of the state with higher energy is given by   

    E^k= E+Δ E f^k,

where f^0=f^3=0, f^1=f^2=1. The parameter f^k gives the time-dependence of the energy and Δ E is the amplitude of the time-dependent part of the energy.

The time-intervals of the period are set to τ_0=τ_2=τ/2 and τ_1=τ_3=τ'→ 0, i.e.,
the energy changes happens after a time interval τ/2 and the temperature changes are instantaneous.
Hence, the number of pieces of the protocol is reduced from four to two. The transition rates for this model, which fulfill 
the generalized detailed balance relation <cit.>, are set to 

    w_+^0=we^-β_c E/2    and      w_-^0=we^β_c E/2,

for the first part of the period; 

    w_+^1=we^-β_h(E+Δ E)/2    and      w_-^1=we^β_h (E+Δ E)/2,

for the second part of the period.
The subscript + indicates a transition rate from the state with energy zero to the state with energy E^k, whereas the subscript - indicates the reversed  transition rate.
The superscript 0 indicates the first half of the the period and the superscript 1 indicates the second half of the period. 
  



The basic physics of the model is that part of the heat taken from the hot reservoir is transformed into extracted work, as explained in <cit.>.
Two currents of interest are the heat current X_q and the work current X_e. 
The piecewise version of the modified generator for a generic current reads  

    ℒ̃_k(z)=([       -w_+^k w_-^ke^-zθ^k;  w_+^ke^zθ^k       -w_-^k ]),

where k=0 for the first half of the period and k=1 for the second half of the period.
The increments for the heat current X_q are defined as θ_q^0≡ 0 and θ_q^1≡β_c(E+Δ E). The increments for the work current X_e 
are defined as θ_e^0≡ 0 and θ_e^1≡-1. 
Hence, we obtain the relation X_q=-β_c(E+Δ E)X_e, which with Eq. (<ref>) leads to 
λ_q(z)=λ_e(-β_c(E+Δ E)z) for the SCGF. In other words, for this simple model there is tight coupling
between the work and heat currents, as shown in Fig. <ref>. 

The exact calculation of the SCGF for a generic piecewise two-state model, is presented in  <ref>. 
For the present model, the SCGF associated with X_q is

    λ_q(z)=1/τln[f(z)+√([f(z)]^2-4e^-(w^0+w^1)/2)/2]

where

    f(z)≡    4e^-τ(w^0+w^1)/4sinh(w^0τ/4)sinh(w^1τ/4)[q_+^1q_-^0e^zβ_c(E+Δ E)+q_-^1q_+^0e^-zβ_c(E+Δ E)]
       +(e^-τ w^0/2q_+^0+q_-^0)(e^-τ w^1/2q_+^1+q_-^1)+(q_+^0+e^-τ w^0/2q_-^0)(q_+^1+q_-^1e^-τ w^1/2),

w^0≡ 2w cosh(β_c E/2), w^1≡ 2w cosh(β_h(E+Δ E)/2), q_±^0≡ w_±^0/w^0, and q_±^1≡ w_±^1/w^1.
The SCGF plotted  in Fig. <ref>  is a concave function of z, which is a generic property of a SCGF, and it becomes linear 
in z for large z, which is a peculiarity of two-state models. 



 §.§ Molecular pump



We now analyze a three-state model for a molecular pump with states i=1,2,3, which is similar to a model analyzed in <cit.>.
The key phenomena that happens in such pumps, is that even though there are no fixed thermodynamic affinities, a suitable time-periodic 
variation of energies E_i(t) and energy barriers B_i(t) can lead to net rotation in the three state system. Periodically driven 
molecular pumps can be realized experimentally with interlocked molecular rings <cit.>. The transition rates are set to   

    w_ii+1=we^E_i(t)-B_i+1(t)    and     w_i+1i=we^E_i+1(t)-B_i+1(t),

where for i=3 we have i+1=1 and the inverse temperature is β=1. The energies are given by 

    E_i(t)= -1+cos[2π (t+(i-1)/3)/τ],

and the energy barriers are given by 

    B_i(t)= cos[2π (t+(i-1)/3)/τ].





The increment of the fluctuating current of interest is such that it increases by one if there is a transition in the clockwise direction
(1→ 2, 2→ 3, and 3→ 1) and it decreases by one if there is a transition in the counter clockwise direction.
The time-dependent modified generator for this current reads  

    ℒ_t(z)=([ -w_12(t) - w_13(t)        w_21(t)e^-z         w_31(t)e^z;         w_12(t)e^z -w_21(t) - w_23(t)        w_32(t)e^-z;        w_13(t)e^-z         w_23(t)e^z -w_31(t) - w_32(t) ]).

We calculated the average current J and diffusion coefficient D using Eq. (<ref>) and Eq. (<ref>), respectively.  
The fundamental matrix ℳ(z) was evaluated with numerical solution of Eq. (<ref>). In Fig. <ref> we show that both quantities  show perfect 
agreement with results obtained from Monte Carlo simulations. Note that the sign of the current depends on the rate constant w that determines 
the speed of the transitions in relation to the period τ: for large enough w the current becomes negative. We have used a 
discrete-time algorithm for our simulations with a sufficiently small time-step,
a continuous-time algorithm for time-dependent transition rates can be found in <cit.>.







 §.§ Biased random walk with time-periodic affinity




In this Section, we consider a biased random walk on a ring with Ω states  driven by a time-periodic affinity 

    F(t)≡ F cos(2π t/τ).

The transition rate for a jump in the clockwise direction is given by 

    w_+(t)= wexp[ψ F(t)/Ω],

whereas the transition rate for 
a jump in the anti-clockwise direction is 

    w_-(t)= wexp[(ψ -1)F(t)/Ω],

where 0≤ψ≤ 1. The parameter ψ determines how F(t)  
influences forward and backward rates. We consider the current in the ring, which increases by one if a jump in the clockwise 
direction takes place and decreases by one if a jump in the anti-clockwise direction takes place. For example, for Ω=3 the modified 
generator   ℒ_t(z) takes the form in Eq. (<ref>) with w_12(t)=w_23(t)=w_31(t)= w_+(t) and 
w_21(t)=w_32(t)=w_13(t)= w_-(t).

This model has the peculiar property that the uniform vector ⟨1| is a left eigenvector of ℒ_t(z), for all t and z, with eigenvalue 

    λ_t(z)= -w_-(t)-w_+(t)+w_-(t)e^-z+w_+(t)e^z.

Hence, from the Dyson series of the ordered exponential and  Eq. (<ref>) for ℳ(z), we obtain

    ⟨1|ℳ(z)= ⟨1|exp(∫_0^τλ_t(z)dt).

Since ⟨1| is a positive vector, from the Perron-Forbenius  theorem, it must be the eigenvector associated with the maximum eigenvalue of ℳ(z), 
which leads to ρ_1(z)=exp(∫_0^τλ_t(z)dt). Using Eq. (<ref>), we then obtain  

    λ(z)= τ^-1∫_0^τλ_t(z)dt.

Explicit evaluation of the above integral leads to 

    λ(z)= w(1-e^-z) [e^zI_0(ψ F_0/Ω)-I_0( F_0(-1+ψ)/Ω)],

where I_0(x) is a modified Bessel function of the first kind. Interestingly, the average current in this model is given by 

    J=λ'(0)= w[I_0(ψ F/Ω)-I_0( F(-1+ψ)/Ω)].

Even though the thermodynamic affinity integrated over a period is zero, the average current can be non-zero: J is 
positive for ψ>1/2 and negative for ψ<1/2.     






The expression in Eq. (<ref>) for λ(z) in terms of an integral of the maximum eigenvalue associated 
with ℒ_t(z) is valid not only for the above model but for any current related to a modified generator that 
fulfils the property  

    ⟨1|ℒ_t(z)= ⟨1|λ_t(z),

for all z and t. Furthermore, from Eq. (<ref>), if the property (<ref>) is satisfied, the SCGF can be written as 

    λ(z)= τ^-1∑_ij∫_0^τ P^inv_i(t)[ℒ_t(z)]_jidt,

where this equation is valid not only for P^inv_i(t) but for an arbitrary probability distribution.
Expression (<ref>) for the evaluation of the SCGF has been proposed in <cit.> as a general expression for the SCGF. 
While it is correct for this peculiar case, in general, Eq. (<ref>) does not provide the correct SCGF. In Fig. <ref>, we show that 
the right hand side of Eq. (<ref>) is different from the SCGF for the model for a molecular pump from 
Section <ref>.   






§ STOCHASTIC PROTOCOL




Hitherto we have restricted to the case of a deterministic protocol. In this section, we consider a system driven by a stochastic
protocol, which is cyclic and, therefore, mimics periodicity. We show that a stochastic protocol with a infinitely large number of jumps is 
equivalent to a deterministic protocol with respect to the large deviations of fluctuating currents. 


For a stochastic protocol, the mathematical model is a bipartite Markov process with time-independent transition rates <cit.>.
The bipartite Markov process has Ω× N states, where N is the number of jumps of the external protocol. A state of the bipartite 
Markov process (i,n) is determined by the variable i that identifies the state of the system and the variable 
n=0,1,…,N-1 that identifies the state of the external protocol. The transition rates 
for a change in the state of the system are defined as 

    w_ij^n≡ w_ij(t=nτ/N),

where w_ij(t) is the transition rate of a corresponding Markov process with time-dependent transition rates that describes a 
deterministic protocol.  

Since the external protocol is stochastic, there is a transition rate associated with changes of the state of the protocol from 
n to n+1. The reversed transition rate is zero. 
From the periodicity of the protocol, for n=N-1 the protocol transitions back to n=0. The transition rate for a change in the protocol 
is set to N/τ, hence, the average time for the protocol to complete a cycle is τ. 

The stationary distribution of state (i,n) is denoted P_i^n, where the dependence on the total number of
jumps N is not shown for a compact notation. The conditional probability of state i given that the protocol 
is in state n is P(i|n)=P_i^n/P^n, where P^n≡∑_iP_i^n=1/N is the stationary probability that the protocol is 
at state n. As shown in <cit.>, in the limit of N→∞, the stationary distribution of the bipartite Markov process is equivalent 
to the invariant periodic distribution of the corresponding Markov process with time-periodic transition rates, i.e., 
P(i|n)→ P^inv_i(t), where τ n/N→ t. 

A stochastic trajectory  of the bipartite process, from time 0 to time T', is denoted by (A,Ξ)_0^T', where A represents the state of the system  
and Ξ represents the state of the protocol. A generic current, analogous to the current in Eq. (<ref>) for a deterministic 
protocol, is defined as

    X_N[(A,Ξ)_0^T']≡∑_0≤ t'≤ T'θ_a_t'^-,a_t'^+^ξ_t',

where θ_ij^n≡θ_ij(t=nτ/N). It can be shown that the SCGF associated with this current 
can be obtained from the tilted generator  ℒ(z) <cit.>, which is a matrix with dimension N×Ω given by
 
    ℒ(z)=([        ℒ_0(z)-𝐈N/τ                  0                  0                  …               𝐈N/τ;               𝐈N/τ      ℒ_τ/N(z)-𝐈N/τ                  …                  0                  0;                  0               𝐈N/τ     ℒ_2τ/N(z)-𝐈N/τ                  …                  0;                  ⋮                  ⋮                  ⋮                  ⋮                  ⋮;                  0                  0                  0                  … ℒ_(N-1)τ/N(z)-𝐈N/τ ]),

where 𝐈 is the identity matrix with dimension Ω and ℒ_nτ/N is equivalent to
the tilted generator from Eq. (<ref>) with θ^n_ij instead of θ_ij(t). 

The maximal eigenvalue associated with (<ref>) is written as  Λ_N(z). 
We now show that lim_N→∞Λ_N(z)= λ(z). The right eigenvector associated with Λ_N(z) is 
written as v⃗_N(z). From the equation ℒ(z)v⃗_N(z)=Λ_N(z)v⃗_N(z) and Eq. (<ref>) we obtain

    N/τ[v_i^n(z)-v_i^n-1(z)]= ∑_j[ℒ_nτ/N(z)]_ijv_j^n(z)-Λ_N(z)v_i^n(z),

where v_i^n(z) are the components of v⃗_N(z). We omit the dependence of these components on N for a compact notation. 
This equation can be written in the form

    N/τ[|v(z)⟩_n-|v(z)⟩_n-1]= ℒ_nτ/N(z)|v(z)⟩_n-Λ_N(z)|v(z)⟩_n,

where |v(z)⟩_n is vector with dimension Ω and components v_i^n(z). In the limit N→∞, we set nτ/N→ t' and 
|v(z)⟩_n→|v(z)⟩_t'. Since for n=N-1 the stochastic protocol jumps back to 
n=0, by construction |v(z)⟩_t'=|v(z)⟩_t'+τ. For N→∞, the vectorial form of Eq. (<ref>) 
then becomes 

    d/dt'|v(z)⟩_t'= [ℒ_t'(z)-lim_N→∞Λ_N(z)]|v(z)⟩_t'.

The formal solution of this equation reads 

    |v(z)⟩_t'= exp(-t'lim_N→∞Λ_N(z))exp(∫_0^t'ℒ_u(z)du) |v(z)⟩_0

Using the periodicity of the eigenvector, i.e., |v(z)⟩_τ=|v(z)⟩_0 and setting t'=τ, we obtain

    ℳ(z)|v(z)⟩_τ= exp(τlim_N→∞Λ_N(z))|v(z)⟩_τ,

where ℳ(z) is the fundamental matrix defined in Eq. (<ref>). Since  by construction |v(z)⟩_τ
is positive, from the Perron-Forbenius theorem 
exp(τlim_N→∞Λ_N(z)) is the maximal eigenvalue associated with ℳ(z), i.e., 

    λ(z)=lim_N→∞Λ_N(z).




In Fig. <ref> we provide two numerical illustrations, one with the model from Section <ref> and the other with the model 
Section <ref>, of the convergence of 
the SCGF for a stochastic protocol with increasing N towards the SCGF for a 
deterministic protocol. There is no generic inequality between Λ_N(z) and Λ(z).   








































§ DISCUSSION




For systems driven by an external periodic protocol, which are well described by Markov processes with time-periodic transition rates, 
large fluctuations of a thermodynamic current can now be determined with the formalism developed in Section <ref>. In particular, 
the SCGF can be evaluated by calculating the maximum Floquet exponent associated with the fundamental matrix ℳ(z).
Cumulants such as the average current J and the diffusion coefficient D can also be directly evaluated from this matrix, which provides 
a numerical method that can be more efficient than Monte Carlo simulations for small systems. Mathematically, beyond fluctuating currents, 
our formalism also applies to generic observables that count number of jumps with increments that do not have to be anti-symmetric. 
    

We have calculated analytically the SCGF in two models: A heat engine with piecewise constant protocol and 
a biased random walk with a time-periodic affinity. Furthermore, we have verified our theoretical results  
by showing agreement between results obtained from numerical evaluation  
of the maximal Floquet exponent and results obtained from Monte Carlo simulations for a model for a molecular pump.

The SCGF associated  with a stochastic protocol with an infinitely large number of jumps is equivalent to the SCGF associated with a deterministic 
protocol, as proved in Section <ref>. This proof provides a rigorous basis to the idea that periodically driven systems can be analyzed with 
the use of stochastic protocols, i.e., if a result about current fluctuations is valid for stochastic protocols than it should be valid for deterministic protocols, 
which can be obtained as a particular limit of a stochastic protocol. The advantage of working with stochastic protocols is that the system and protocol together 
form a bipartite Markov process with constant transition rates that reach a nonequlibrium stationary state, which are quite well known. 








§ WORK AND RELATED CURRENTS



A fluctuating current of interest in stochastic thermodynamics is the work done on the system 
due to the time variation of the energy levels. For jump processes, such current can be written in the form 

    X'[A_0^T]≡∫_0^Tf_a_t(t)dt,

where  f_i(t)= ∂_tg_i(t) and g_i(t) is periodic with period τ. For the case of 
work g_i(t) is the free energy of state i. 

The empirical density ρ_i(t) and the empirical flow C_ij(t) are the number of periods for 
which the trajectory is in state i at time t∈ [0,τ] and the number of transitions from 
i to j at time t, respectively. They are functionals of the stochastic trajectory 
from time 0 to time T=nτ that are defined as 

    ρ_i(t)≡1/n∑_k=0^nδ_a_kτ+t,i,

and

    C_ij(t)≡1/n∑_k=0^nδ_a_(kτ+t)^-,iδ_a_(kτ+t)^+,j.

As shown in <cit.>, in the large deviation regime, they fulfill the constraint

    d/dtρ(i,t)= ∑_j[C_ji(t)-C_ij(t)],

for every single stochastic trajectory. Moreover, both C_ji(t) and ρ_i(t) are periodic with 
period τ.  
  
Using the empirical density in Eq. (<ref>) we can rewrite Eq. (<ref>) as 

    X'[A_0^T]=n∑_i∫_0^τ∂_tg_i(t)ρ_i(t)dt,

where we used f_i(t)= ∂_tg_i(t). From the constraint in Eq. (<ref>) 
and the definition of the empirical flow in Eq. (<ref>), we obtain  

    X'[A_0^T]=∑_0≤ t≤ nτ(f_t(a_t^+)-f_t(a_t^-)),

which is a current of the form given in Eq. (<ref>). Hence, the current in Eq. (<ref>) can be 
written as in Eq. (<ref>) if f_i(t)= ∂_tg_i(t). This relation also holds for a stochastic 
protocol <cit.>. In general, if f_i(t) is not a derivative of a periodic function, this relation 
may not hold. Whereas extending the results from Sec. <ref> to such case should be straightforward, as we have done for 
diffusion processes in <ref>, proving the equivalence with stochastic protocols for this more general current remains 
an open challenge. 







§ DIFFUSION PROCESSES



We consider a d-dimensional diffusion process Y_t that follows the Fokker-Planck equation 

    ∂_tμ_t(y)=ℒ_t[μ_t](y)

where μ_t(y) is the probability density at time t and 
the adjoint of the generator is given by 

    ℒ^†_t=F_t∇+∇D_t/2∇.

where F_t is a time-periodic the drift vector, D_t is a time-periodic diffusion matrix, τ is the period, and ∇ is the nabla operator. 
This generator is then time-periodic. Hence, the invariant probability density related to the long time limit 
has the property μ^inv_t(y)=μ^inv_t+τ(y).

A stochastic current is a functional of the stochastic trajectory Y_0^nτ, from time 0 to time nτ,
defined as 

    X[Y_0^nτ]≡∫_0^nτ(f_t(Y_t)dt+g_t(Y_t)∘ dY_t),

where we use ∘ for the Stratonovich convention, f_t=f_t+τ is a scalar function, and g_t=g_t+τ is a vector field.
The Fokker-Planck current associated with ρ_t^inv is a vector  

    J_ρ_t^inv,t(y)≡F_t(y)ρ_t^inv(y)-1/2D_t(y)(∇ρ_t^inv)(y).

The typical behavior of X[Y_0^nτ] is given by 

    lim_n→∞X[Y_0^nτ]/nτ=1/τ∫_0^τdt(∫ρ_t^inv(y)f_t(y)dy +∫ J_ρ_t^inv,t(y)g_t(y)dy).





The SCGF is defined as 

    λ(z)≡lim_n→∞1/nτln(⟨exp[ zX_nτ]⟩).

Using the Feyman-Kac formula and the Grinasov lemma it can be shown that <cit.> 

    ⟨e^zXδ(Y_T-y')|Y_0=y⟩= [exp(∫_0^nτdtℒ_t(z))](y',y).

where the tilted generator ℒ_t(z) is the second order differential operator of the form

    ℒ^†_t(z)=zf_t+F_t(∇+zg_t)+(∇+zg_t)D_t/2(∇+zg_t).

Eq. (<ref>) is analogous to Eq. (<ref>) for jump processes. The fundamental operator is defined as 

    ℳ(z)≡exp(∫_0^τdtℒ_t(z)).

From Eq. (<ref>), this operator is also a Perron-Forbenius operator. Similar to the case for jump processes, ℳ(z) can 
be expanded in the form given by Eq. (<ref>). The real maximum eigenvalue of ℳ(z) is denoted ρ_1(z). Following the same 
procedure for a jump process one can define a Laplace transform of the probability of current, analogous to Eq. (<ref>) (with an integral instead 
of a sum). The time evolution of this Laplace transform follows Eq. (<ref>). Hence, the SCGF λ(z) can be expressed in terms of 
the maximum eigenvalue ρ_1(z) trough the formula (<ref>). 










§ EXACT EXPRESSIONS FOR TWO-STATE MODEL WITH PIECEWISE CONSTANT PROTOCOL



We consider a generic two-state model with a piecewise constant protocol (see <cit.> for similar 
calculations for a protocol with two pieces). The transition 
rates during the interval τ_k are denoted by w_+^k and w_-^k, where w_+^k is 
the transition rate from state 1 to state 2. The increment of 
the current from state 1 to state 2  during the interval τ_k is θ^k and the increment for 
the reversed transition is -θ^k. The piecewise version of the modified generator reads  

    ℒ̃_k(z)=([       -w_+^k w_-^ke^-zθ^k;  w_+^ke^zθ^k       -w_-^k ]).

From this matrix, we obtain 

    exp(-τ_kℒ̃_n(z))=1/w^k([      e^-τ_kw^kw_+^k+w_-^k e^-zθ^k(1-e^-τ_kw^k)w_-^k;  e^zθ^k(1-e^-τ_kw^k)w_+^k      w_+^k+e^-τ_kw^kw_-^k ])

where w^k≡ w_-^k+w_+^k.
The fundamental matrix for a piecewise protocol in Eq. (<ref>) is a product of the matrices in Eq. (<ref>). 
Expressing the maximum eigenvalue of the two by two fundamental matrix in terms of its trace and determinant, and 
using the  Abel-Jacobi-Liouville identity for the determinant

    (ℳ(z))= exp(∫_0^τ Tr(ℒ_t(z)) dt)=∏_n=0^N-1e^-τ_k w^k,

we obtain the following expression for the scaled cumulant generating function 

    λ(z)=1/τln[Tr(ℳ(z))+√([Tr(ℳ(z))]^2-4∏_n=0^N-1e^-τ_k w^k)/2].

Hence, we can calculate λ(z) from the trace of the fundamental matrix in Eq. (<ref>) with Eq. (<ref>).
In particular, for N=2  the trace of the fundamental matrix is given by  

    Tr(ℳ(z))  =     4e^-(τ_0 w^0+τ_1 w^1)/2sinh(w^0τ_0/2)sinh(w^1τ_1/2)[q_+^1q_-^0e^z(θ^1-θ^0)+q_-^1q_+^0e^z(θ^0-θ^1)]
       +(e^-τ_0 w^0q_+^0+q_-^0)(e^-τ_1 w^1q_+^1+q_-^1)+(q_+^0+e^-τ_0 w^0q_-^0)(q_+^1+q_-^1e^-τ_1 w^1),

where q_±^0≡ w_±^0/w^0 and q_±^1≡ w_±^1/w^1.

For the model for a heat engine in Section <ref>  the protocol has N=4 pieces. Since the changes in temperature are instantaneous, this number 
is reduced to N=2, with τ_0=τ_1=τ/2. The transition rates for this model w_+^0=we^-β_c E/2,  w_-^0=we^β_c E/2,
w_+^1=we^-β_h(E+Δ E)/2, and w_-^1=we^β_h (E+Δ E)/2. For the current X_q the increments,
are   θ_q^0=0 and θ_q^1=β_c(E+Δ E). For the current X_e the increments are θ_e^0=0 and θ_e^1=-1. 







































 




 















































§ REFERENCES


10
url<#>1#1urlprefixURL 





erba15
Erbas-Cakmak S, Leigh D A, McTernan C T and Nussbaumer A L 2015 Chem.
  Rev. 115 10081 
  

mart17
Martinez I A, Roldan E, Dinis L and Rica R A 2017 Soft Matter 13 22 

seif12
Seifert U 2012 Rep. Prog. Phys. 75 126001

gamm98
Gammaitoni L, Hänggi P, Jung P and Marchesoni F 1998 Rev. Mod. Phys.
  70 223

izum09
Izumida Y and Okuda K 2009 Phys. Rev. E 80 021121

izum10
Izumida Y and Okuda K 2010 Eur. Phys. J. B 77 499

izum15
Izumida Y and Okuda K 2015 New J. Phys. 17 085011

bran15
Brandner K, Saito K and Seifert U 2015 Phys. Rev. X 5 031019
  
proe15
Proesmans K and Van den Broeck C 2015 Phys. Rev. Lett. 115
  090601
  

proe16
Proesmans K, Cleuren B and Van den Broeck C 2016 J. Stat. Mech.  023202

ray17
Ray S and Barato A C 2017 Phys. Rev. E 96 052120
  

schm08
Schmiedl T and Seifert U 2008 EPL 81 20003

espo10
Esposito M, Kawai R, Lindenberg K and van den Broeck C 2010 Phys. Rev.
  E 81 041106

tu14
Tu Z C 2014 Phys. Rev. E 89 052148
  

raz16a
Raz O, Subaşı Y and
  Pugatch R 2016 Phys. Rev. Lett. 116 160601
  

raha08
Rahav S, Horowitz J and Jarzynski C 2008 Phys. Rev. Lett. 101
  140602

cher08
Chernyak V Y and Sinitsyn N A 2008 Phys. Rev. Lett. 101 160601

maes10
Maes C, Netocný K and Thomas S R 2010 J Chem. Phys. 132 234116
  
raha11
Rahav S 2011 J. Stat. Mech.: Theor. Exp.  P09020

mand14
Mandal D 2014 EPL 108 50001

raz16
Raz O, Subaşı Y and
  Jarzynski C 2016 Phys. Rev. X 6 021022


rots17
Rotskoff G M 2017 Phys. Rev. E 95 030101

sing10
Singh N and Wynants B 2010 J. Stat. Mech.  P03007

bert18
Bertini L, Chetrite R, Faggionato A and Gabrielli D 2018 arXiv:1710.08001

lebo99
Lebowitz J L and Spohn H 1999 J. Stat. Phys. 95 333

koza99
Koza Z 1999 J. Phys. A: Math. Gen. 32 7637

spec12
Speck T, Engel A, and Seifert U 2012 J. Stat. Mech.: Theor. Exp.  P12001


draz92
Drazin P G 1992 Nonlinear systems (Cambridge University Press, Cambridge)

grim90
Grimshaw R 1990 Nonlinear ordinary differential equations (Blackwell
  Scientific Publications, Oxford)


grif98
Grifoni M and Hänggi P 1998 Phys. Rep. 304 
  229

klau08
Klausmeier C A 2008 Theor. Ecol. 1 153
 

verl14
Verley G, den Broeck C V and Esposito M 2014 New J. Phys.
  16 095001
  

bara16
Barato A C and Seifert U 2016 Phys. Rev. X 6 041053
  
bara17b
Barato A C and Seifert U 2017 New J. Phys. 19 073021

verl13
Verley G, Van den Broeck C and Esposito M 2013 Phys. Rev. E 88 032137 

bara15
Barato A C and Seifert U 2015 J. Phys. Chem. B 119 6555

wach15 
Wachtel A, Vollmer J and Altaner B 2015 Phys. Rev. E 92 042132
  

prad97
Prados A, Brey J J and Sanchez-Rey B 1997 J. Stat. Phys. 89 709

chet15
Chetrite R and Touchette H 2015 Ann. Henri Poincaré 16
  2005






