Spectral Analysis of Current Fluctuations in Periodically Driven Stochastic Systems

Bertrand Lacroix-A-Chez-Toine Department of Physics of Complex Systems, Weizmann Institute of Science, Rehovot 7610001, Israel Oren Raz Department of Physics of Complex Systems, Weizmann Institute of Science, Rehovot 7610001, Israel

Abstract

Current fluctuations play an important role in non-equilibrium statistical mechanics, and are a key object of interest in both theoretical studies and in practical applications. So far, most of the studies were devoted to the fluctuations in the time-averaged current – the zero frequency Fourier component of the time dependent current. However, in many practical applications the fluctuations at other frequencies are of equal importance. Here we study the full frequency dependence of current statistics in periodically driven stochastic systems. First, we show a general method to calculate the current statistics, valid even when the current’s frequency is incommensurate with the driving frequency, breaking the time periodicity of the system. Somewhat surprisingly, we find that the cumulant generating function (CGF), that encodes all the statistics of the current, is composed of a continuous background at any frequency accompanied by either positive or negative peaks at current’s frequencies commensurate with the driving frequency. We show that cumulants of increasing orders display peaks at an increasing number of locations but with decreasing amplitudes that depend on the commensurate ratio of frequencies. All these peaks are then transcribed in the behaviour of the CGF. As the measurement time increases, these peaks become sharper but keep the same amplitude and eventually lead to discontinuities of the CGF at all the frequencies that are commensurate with the driving frequency in the limit of infinitely long measurement. We demonstrate our formalism and its consequences on three types of models: an underdamped Brownian particle in a periodically driven harmonic potential; a periodically driven run-and-tumble particle; and a two-state system.

I Introduction

Equilibrium statistical mechanics is a mature framework that provides a calculation scheme for the ensemble averages and the fluctuations of many quantities of interest, with only mild assumptions on the system, e.g. thermodynamic equilibrium and short range interactions. In practice, however, equilibrium systems are rather the exception than the norm. Despite many efforts, no equivalent general framework is known for systems which are far from thermal equilibrium. Characterising the mean value and fluctuations for most quantities of interest remains a challenging task, even in the relatively simple cases of a system in a non-equilibrium steady state Zia and Schmittmann [2007]; Ruelle [2003]; Zia and Schmittmann [2006] or a periodically driven system Brandner et al. [2015]; Raz et al. [2016]; Busiello et al. [2018].

Important progress has been achieved in the 50’s for systems which are close to thermal equilibrium, in the form of the fluctuation-dissipation relations Kubo [1966]. These relate the equilibrium fluctuations of a system to the rate of dissipation when the system is weakly driven away from equilibrium. These fluctuation-dissipation relations where generalized in the 90’s to the Gallavotti-Cohen fluctuation relations, which constitute rare exact results holding arbitrarily far from equilibrium. They imply a fundamental symmetry on the distributions of the entropy production Gallavotti and Cohen [1995]; Kurchan [1998] and of the currents Andrieux and Gaspard [2007a] in the system and naturally extend the fluctuation dissipation relations. They hold for non-equilibrium steady state systems Gallavotti and Cohen [1995]; Lebowitz and Spohn [1999]; Andrieux and Gaspard [2007b] as well as periodically driven systems Shargel and Chou [2009], although they are not expected to hold for currents in the latter case Harris and Schütz [2007].

The thermodynamic uncertainty relations, bounding the current fluctuations in the system using the entropy production Barato and Seifert [2015]; Gingrich et al. [2016]; Horowitz and Gingrich [2019] constitute another important progress that was recently achieved and holds arbitrarily far from equilibrium. They imply that decreasing the current fluctuations in the steady state comes at a cost in terms of the entropy production. These relations were extended recently to periodically driven systems in Proesmans and Horowitz [2019]; Barato et al. [2018]; Koyuk and Seifert [2019, 2020]; Koyuk et al. [2018]. Both the Gallavotti-Cohen and the thermodynamic uncertainty principle hold for arbitrary systems, but in practice they are of most interest in the case of small systems where large fluctuations occur with a non-vanishing probability.

In this study we consider the fluctuations of the current’s Fourier modes in stochastic non-equilibrium systems that are subject to time-periodic driving. The fluctuations of the zero frequency current (commonly referred to as the “DC current”) are tightly related to both the Gallavotti-Cohen and the thermodynamic uncertainty relations described above, and thus were studied intensively for systems which are subject to a time-independent thermodynamic forcing that breaks detailed balance, driving them into a non-equilibrium steady state Nyawo and Touchette [2016]; Nemoto and Sasa [2011a, b]; Proesmans and Derrida [2019]; Derrida and Sadhu [2019]; Nyawo and Touchette [2017]; Fischer et al. [2018]. However, only few works have been dedicated to other Fourier modes of the currents. To the best of our knowledge, only the zero-frequency current fluctuations Barato and Chetrite [2018]; Bertini et al. [2018] and the fluctuations of the current’s Fourier mode with the same frequency as that of the drive Chabane et al. [2020] have been investigated.

The motivation to study all of the current’s Fourier modes is natural from a theoretical perspective, but also from a practical one, as often these fluctuations have important implications. An interesting example is the Paul trap, commonly used to trap ions Blatt and Wineland [2008]. Ions in this trap display heating processes commonly orders of magnitude higher than expected Brownnutt et al. [2015] which are yet poorly understood. In this system ions are trapped by imposing an alternating current at high frequency giving also rise to fluctuations at all other frequencies. Although the driving frequency is far from the frequency that can heat the ions, the generated current fluctuation at the ion’s frequency might be an important component in the heating process.

In this manuscript we study the frequency dependence of current statistics when a stochastic Markovian system is periodically driven. Our main object of interest is the Cumulant generating function (CGF) of the current Fourier modes, which encodes the full information on the current statistics. First, we use the Oseledets theorem to extend the framework of the current cumulant generating function from the case where the driving frequency is commensurate with the Fourier mode of the current, to account for non-commensurate frequency ratios as well. Surprisingly, we find that this function has a structure similar to the famous Thomae’s function Beanland et al. [2009], namely it is composed of a continuous “background” for frequencies which are incommensurate with respect to the periodic driving, accompanied by (positive or negative) discontinuous “spikes” at commensurate frequencies. The integer denominator of the commensurate ratio between the frequencies at which the spikes appear can be directly related to the order of the cumulant that has a discontinuous behaviour at this frequency. Moreover, we explain how our results are modified for a finite time experiment, where the discontinuous spikes are smoothed into continuous peaks that gets sharper as the duration of the measurement increases. This general behavior is demonstrated on a Brownian particle in a periodically driven harmonic potential and a periodically driven trapped run-and-tumble particle, where analytical results for the first two cumulants can be obtained, and in a two-state system where the cumulant generating function can be numerically evaluate.

The setup studied in this manuscript is somewhat similar to the stochastic resonance setup, where a bi-stable system coupled to a noise source is weakly perturbed by time-periodic forcing (See Gammaitoni et al. [1998] for a review on stochastic resonance). However, the main object of interest in stochastic resonance is the response of the system to the driving as a function of the noise, where the response of the system is usually the transitions between the two meta-stable states at the driving frequency. In other words, the interest in stochastic resonance is in the driving frequency Fourier mode of the current between the two meta-stable state. In contrast, we consider a quite generic periodic driving, which is not necessarily weak and does not necessarily have two meta-stable states. We also consider all the Fourier modes of the current statistics, and not only the average current as in stochastic resonance. Furthermore, we are interested in the frequency dependence of the fluctuations.

The manuscript is organised as follows: In section II we detail the general setting that we consider and give the main technical background needed to obtain our results. In section III we provide a short summary of our main results and their implications. In section IV we derive our main results by obtaining general expressions of the cumulants. In section V we show how to apply our results to simple exactly solvable systems. Finally in section VI we conclude and give some insights on new unexplored directions. Some technical details of the computations and complicated expressions are relegated to the Appendices.

II Setup – Periodically driven Markovian systems

II.1 Markov propagator

In this manuscript we consider the statistics of the current’s Fourier modes in a periodically driven stochastic system. To this end, we assume that the evolution of the probability distribution associated with the system is Markovian, as detailed in what follows. The micro-state of the system, for example the position of a particle in space, is denoted by the vector ${\bf x}$ . We assume that the space of all ${\bf x}$ is simply connected ¹¹1This assumption implies some consequences on our results. Although a generalization to non-simply connected domains is possible, we do not consider these cases here. The probability distribution of all the micro-states at time $t$ is uniquely characterised by the initial state of the system, which for simplicity we assume to be a specific micro-state ${\bf x_{0}}$ , as well as by the propagator, $G({\bf x},t|{\bf x}_{0},t_{0})$ , which is the transition probability density to be in a final state ${\bf x}$ at time $t$ , provided that the system was at ${\bf x}_{0}$ at time $t_{0}$ . It satisfies the evolution equation,

	$\displaystyle\partial_{t}G={\cal L}({\bf x},\nabla_{\bf x},t)G\,,$		(1)
	$\displaystyle{\cal L}({\bf x},\nabla_{\bf x},t)=-\sum_{\beta}\partial_{x_{\beta}}F_{\beta}({\bf x},t)+\sum_{\beta,\gamma}\partial_{x_{\beta},x_{\gamma}}D_{\beta\gamma}({\bf x},t)\;,$

where here and in the following we use Greek letters for the indices of the micro-state space and with the initial condition

\,\,G({\bf x},t_{0}|{\bf x}_{0},t_{0})=\delta^{d}({\bf x}-{\bf x}_{0})\,,

(2)

where $\delta^{d}({\bf x})=\prod_{\beta=1}^{d}\delta(x_{\beta})$ is the $d$ -dimensional Dirac delta function. By virtue of probability conservation, the propagator is normalised to unity,

\forall t\geq t_{0}\;,\;\;\int d{\bf x}\,G({\bf x},t|{\bf x}_{0},t_{0})=1\;.

(3)

Throughout the manuscript we assume that the force $F_{\beta}({\bf x},t)$ and diffusion $D_{\beta\gamma}({\bf x},t)$ are both explicitly time dependent, non-singular and time-periodic with cycle time $T_{d}$ . We denote the corresponding fundamental angular frequency by $\omega_{d}=(2\pi)/T_{d}$ , and we refer to it as the driving frequency, as this time dependence is what drives the system out of equilibrium. We are interested in the cases where the system has a unique time-periodic stationary state in the limit $t\to\infty$ , and that it does not depend on the initial position ${\bf x}_{0}$ .

Equation (1) can be formally solved by

G({\bf x},t|{\bf x}_{0},t_{0})={\cal T}e^{\int_{t_{0}}^{t}d\tau{\cal L}({\bf x},\nabla_{\bf x},\tau)}\delta^{d}({\bf x}-{\bf x}_{0})\;,

(4)

where ${\cal T}$ is the time-ordered product. Exploiting the discrete time translation symmetry of the system, ${\cal L}({\bf x},\nabla_{\bf x},t+T_{d})={\cal L}({\bf x},\nabla_{\bf x},t)$ for any time $t$ , we expand the propagator using Floquet theory as Caceres and Lobos [2006]; Kim et al. [2010]

G({\bf x},t|{\bf x}_{0},t_{0})=\sum_{k=0}^{\infty}e^{-\lambda_{k}(t-t_{0})}f_{k}({\bf x},t)g_{k}({\bf x}_{0},t_{0})\;,

(5)

where the functions $f_{k}$ and $g_{k}$ are the eigenvectors $U_{T_{d}}f_{k}=e^{-\lambda_{k}T_{d}}f_{k}$ and $U_{T_{d}}^{\dagger}g_{k}=e^{\lambda_{k}T_{d}}g_{k}$ of the operators

	$\displaystyle U_{T_{d}}(t)={\cal T}e^{\int_{t}^{t+T_{d}}d\tau{\cal L}({\bf x},\nabla_{\bf x},\tau)}\;,$		(6)
	$\displaystyle U_{T_{d}}^{\dagger}(t_{0})={\cal T}e^{-\int_{t_{0}}^{t_{0}+T_{d}}d\tau{\cal L}^{\dagger}({\bf x}_{0},\nabla_{{\bf x}_{0}},\tau)}\;,$		(7)

where

{\cal L}^{\dagger}({\bf x},\nabla_{\bf x},t_{0})=\sum_{\beta}F_{\beta}({\bf x},t_{0})\partial_{x_{\beta}}+\sum_{\beta,\gamma}D_{\beta\gamma}({\bf x},t_{0})\partial_{x_{\beta},x_{\gamma}}

(8)

is the adjoint of the evolution operator ${\cal L}({\bf x},\nabla_{\bf x},t)$ . The Floquet eigenvalues $\lambda_{k}$ ’s are ordered in the following way,

\real(\lambda_{0})=0<\real(\lambda_{1})\leq\real(\lambda_{2})\leq\cdots

(9)

and they all have positive real part. Using our assumptions on the periodic solution being unique, we obtain that the eigenvalue $\lambda_{0}=0$ is non-degenerate. On the other hand, as $(t-t_{0})\to\infty$ , we obtain that

G({\bf x},t|{\bf x}_{0},t_{0})\xrightarrow[t\rightarrow\infty]{}f_{0}({\bf x},t)g_{0}({\bf x}_{0},t_{0})\;.

(10)

As this stationary state is independent of ${\bf x}_{0}$ , we must have that $g_{0}({\bf x}_{0},t_{0})$ is independent of both ${\bf x}_{0}$ and $t_{0}$ . It is convenient to choose $g_{0}({\bf x}_{0},t_{0})=1$ , so that the associated eigenfunction $f_{0}({\bf x},t)$ is positive and normalised to unity. Note that we assume this steady state to be non-singular, i.e. $f_{0}({\bf x},t)>0$ at least in a finite domain. More generally, exploiting the identities

	$\displaystyle G({\bf x},t\|{\bf x}_{0},t)=\delta^{d}({\bf x}-{\bf x}_{0})$		(11)
	$\displaystyle\int d{\bf y}\,G({\bf x},t\|{\bf y},\tau)G({\bf y},\tau\|{\bf x}_{0},t_{0})=G({\bf x},t\|{\bf x}_{0},t_{0})$		(12)

one can show that the eigenfunctions $f_{k}({\bf x},t)$ and $g_{k}({\bf x}_{0},t)$ are bi-orthogonal and satisfy a completeness relation

	$\displaystyle\sum_{k}f_{k}({\bf x},t)g_{k}({\bf x}_{0},t)=\delta^{d}({\bf x}-{\bf x}_{0})\;,$		(13)
	$\displaystyle\int d{\bf x}\,f_{k}({\bf x},t)g_{l}({\bf x},t)=\delta_{k,l}\;,$		(14)

where we remind that $\delta^{d}({\bf x})=\prod_{\beta=1}^{d}\delta(x_{\beta})$ is the $d$ -dimensional Dirac delta function while $\delta_{k,l}=1$ if $k=l$ and zero otherwise is the Kronecker delta.

II.2 Alternating charge

So far we discussed the probability propagator in periodically driven systems. However, in this manuscript we are interested in the fluctuations of the oscillating current at an arbitrary frequency $\omega_{c}$ . In practice, one needs to consider the associated time-integrated quantity which, in analogy to the DC case, we denote by the fluctuating charge. For a given realisation, the empirical $\omega_{c}$ -alternating charge is defined as

	$\displaystyle Q_{\omega_{c};\alpha}(t)=\int_{0}^{t}d\tau\,\partial_{\tau}x_{\alpha}(\tau)\,\cos(\omega_{c}\tau)$		(15)
	$\displaystyle=\left[x_{\alpha}(\tau)\cos(\omega_{c}\tau)\right]_{\tau=0}^{\tau=t}+\omega_{c}\int_{0}^{t}d\tau\,x_{\alpha}(\tau)\,\sin(\omega_{c}\tau)\;,$

where $x_{\alpha}(\tau)$ is the $\alpha^{\rm th}$ component of ${\bf x}(\tau)$ ²²2Note that a similar quantity can be defined with a $\sin$ function instead of the $\cos$ .. The random variable $Q_{\omega_{c};\alpha}(t)$ is a functional of the whole path realisation ${\bf x}(\tau)$ for all $\tau\in[0,t]$ . However, $Q_{\omega_{c};\alpha}(t)$ can be calculated for any specific realisation of the system using the above integral.

For $\omega_{c}=0$ , the empirical alternating charge coincides with the direct empirical charge. As we consider a simply-connected space, this direct charge is path independent, given the values of the initial $x_{\alpha}(0)$ and final $x_{\alpha}(t)$ positions, and is equal to $Q_{\omega_{c}=0;\alpha}(t)=x_{\alpha}(t)-x_{\alpha}(0)$ . The probability distribution function (PDF) $P_{\omega_{c}=0,\alpha}(Q;t)$ of this direct charge $Q_{0,\alpha}(t)$ is given, in the large time limit, by

P_{\omega_{c}=0,\alpha}(Q;t)\to\int d{\bf x}\,f_{0}({\bf x},t)\delta(x_{\alpha}-x_{0,\alpha}-Q)\;.

(16)

For a generic value of the frequency $\omega_{c}$ , the empirical alternating charge is path dependent, and therefore computing the probabilty distribution function $P_{\omega_{c},\alpha}(Q;t)$ of the alternating charge is a non-trivial task. It turns out to be simpler in this case to consider instead its moment generating function (MGF)

\displaystyle\Big{\langle}e^{\mu Q_{\omega_{c};\alpha}(t)}\Big{\rangle}_{{\bf x}_{0}}

\displaystyle\equiv\int_{-\infty}^{\infty}dQ\,e^{\mu Q}P_{\omega_{c},\alpha}(Q;t),

(17)

where the average $\langle\cdots\rangle_{{\bf x}_{0}}$ is taken over all the stochastic trajectories starting from the same initial state ${\bf x}(0)={\bf x}_{0}$ . This generating function can also be written as

\displaystyle\Big{\langle}e^{\mu Q_{\omega_{c};\alpha}(t)}\Big{\rangle}_{{\bf x}_{0}}

\displaystyle=\int d{\bf x}\,G_{\mu}({\bf x},t|{\bf x}_{0})\,,

(18)

where the function $G_{\mu}({\bf x},t|{\bf x}_{0})$ evolves in time according to the so-called ”tilted” operator Chetrite and Touchette [2013]

	$\displaystyle\partial_{t}G_{\mu}={\cal L}_{\mu}({\bf x},\nabla_{\bf x},t)G_{\mu}\,,$		(19)
	$\displaystyle G_{\mu}({\bf x},0\|{\bf x}_{0})=\delta^{d}({\bf x}-{\bf x}_{0})\,.$		(20)

The tilted operator can be written as

\displaystyle{\cal L}_{\mu}({\bf x},\nabla_{\bf x},t)=\sum_{k=0}^{2}\left[\mu\cos(\omega_{c}t)\right]^{k}{\cal L}^{k}({\bf x},\nabla_{\bf x},t)\;,

(21)

where

	$\displaystyle{\cal L}^{0}={\cal L}({\bf x},\nabla_{\bf x},t)\;,$		(22)
	$\displaystyle{\cal L}^{1}=F_{\alpha}({\bf x},t)-\sum_{\beta}\partial_{x_{\beta}}\left(D_{\alpha\beta}({\bf x},t)+D_{\beta\alpha}({\bf x},t)\right)\;,$
	$\displaystyle{\cal L}^{2}=D_{\alpha\alpha}({\bf x},t)\;.$

Similarly to the case of the propagator, one can formally express the solution as

G_{\mu}({\bf x},t|{\bf x}_{0})={\cal T}e^{\int_{0}^{t}d\tau{\cal L}_{\mu}({\bf x},\nabla_{\bf x},t)}\delta^{d}({\bf x}-{\bf x}_{0})\;.

(23)

We note, however, that in contrast with ${\cal L}({\bf x},\nabla_{\bf x},t)$ which is periodic by construction, ${\cal L}_{\mu}({\bf x},\nabla_{\bf x},t)$ is not necessarily time periodic, as it has both periodic components varying with frequency $\omega_{c}$ and frequency $\omega_{d}$ . Only in the case where the ratio between these frequencies is commensurate, namely when it is given by a rational number, does ${\cal L}_{\mu}({\bf x},\nabla_{\bf x},t)$ become periodic. In what follows, we first address the commensurate case, and then the more challenging case of non-commensurate frequencies.

II.3 Commensurate frequencies

Let us first consider the case where the frequencies of the drive $\omega_{d}$ and of the charge $\omega_{c}$ are commensurate. In this case there exist two natural numbers, $n,m\in\mathbb{N}$ , such that

T=nT_{d}=mT_{c}<\infty\;.

(24)

In this case the tilted operator is periodic with a finite period $T$ , namely ${\cal L}_{\mu}({\bf x},\nabla_{\bf x},t+T)={\cal L}_{\mu}({\bf x},\nabla_{\bf x},t)$ . One may then use the Floquet theory by defining the operators

	$\displaystyle U_{T;\mu}(t)={\cal T}e^{\int_{t}^{t+T}d\tau{\cal L}_{\mu}({\bf x},\nabla_{\bf x},\tau)}\;,$		(25)
	$\displaystyle U_{T;\mu}^{\dagger}(t_{0})={\cal T}e^{-\int_{t_{0}}^{t_{0}+T}d\tau{\cal L}_{\mu}^{\dagger}({\bf x}_{0},\nabla_{{\bf x}_{0}},\tau)}\;.$		(26)

Introducing the common ordered set of Floquet eigenvalues of these operators $\{\lambda_{k}(\mu)\,,\,\,\real[\lambda_{0}(\mu)]<\real[\lambda_{1}(\mu)]\leq\cdots\}$ and their associated respective eigenfunctions $f_{k}^{\mu}({\bf x},t)$ and $g_{k}^{\mu}({\bf x}_{0},t)$ , we may express the MGF as

\Big{\langle}e^{\mu Q_{\omega_{c};\alpha}(t)}\Big{\rangle}=\sum_{k}e^{-\lambda_{k}(\mu)t}a_{k}(t;{\bf x}_{0})\;,

(27)

where the functions $a_{k}(t;{\bf x}_{0})=\int d{\bf x}f_{k}^{\mu}({\bf x},t)g_{k}^{\mu}({\bf x}_{0},0)$ are periodic in $t$ with period $T$ . In contrast to the propagator $G({\bf x},t|{\bf x}_{0})$ , the tilted propagator $G_{\mu}({\bf x},t|{\bf x}_{0})$ does not propagate a probability distribution, and its largest eigenvalue $\lambda_{0}(\mu)$ is non zero in general, apart for $\mu=0$ where ${\cal L}_{\mu=0}({\bf x},\nabla_{\bf x},t)={\cal L}({\bf x},\nabla_{\bf x},t)$ . Using the Floquet theory, the cumulant generating function (CGF) in the large time limit reads

	$\displaystyle\chi_{\mu}(\omega_{c})$	$\displaystyle\equiv\lim_{t\to\infty}\frac{1}{t}\ln\Big{\langle}e^{\mu Q_{\omega_{c};\alpha}(t)}\Big{\rangle}_{\bf x_{0}}=-\lambda_{0}(\mu)$		(28)
		$\displaystyle=\sum_{p=1}^{\infty}{\cal Q}_{p}^{\alpha}(\omega_{c})\,\mu^{p}\;,$

where ${\cal Q}_{p}^{\alpha}(\omega_{c})$ corresponds to the large time limit of the rescaled cumulant of the alternating charge

{\cal Q}_{p}^{\alpha}(\omega_{c})\equiv\lim_{t\to\infty}\frac{\Big{\langle}Q_{\omega_{c};\alpha}^{p}(t)\Big{\rangle}_{\rm co}}{t}=-\left.\partial_{\mu}^{p}\lambda_{0}(\mu)\right|_{\mu=0}\;.

(29)

Here, the subscript _co refers to a connected correlation function, and $p$ is the order of the cumulant. The limit is well-defined as the cumulants scale at most linearly with time. The specific form of the CGF obtained here is consistent with the probability distribution function (PDF) $P_{\omega_{c}}(Q;t)$ of the alternating charge $Q_{\omega_{c};\alpha}(t)$ taking a large deviation scaling form in the large time limit. By definition, the MGF

\Big{\langle}e^{\mu Q_{\omega_{c};\alpha}(t)}\Big{\rangle}=\hat{P}_{\omega_{c}}(-i\mu;t)\;,

(30)

where $\hat{P}_{\omega_{c}}(k;t)$ is the Fourier transform of the PDF $P_{\omega_{c}}(Q;t)$ at the specific value $k=-i\mu$ . Taking the inverse Fourier transform, we obtain that the atypical fluctuations of the alternating charge are described by the large deviation form,

P_{\omega_{c}}(Q;t)=\int_{-i\infty}^{i\infty}\frac{d\mu}{2i\pi}e^{-\mu Q}\Big{\langle}e^{\mu Q_{\omega_{c};\alpha}(t)}\Big{\rangle}\asymp e^{-t\,\Phi_{\omega_{c}}\left(\frac{Q}{t}\right)}\;,

(31)

where we have used the notation (common in large deviation theory)

A(x,t)\asymp e^{-tB\left(\frac{x}{t}\right)}\Leftrightarrow-\lim_{t\to\infty}\frac{1}{t}\ln A(at,t)=B\left(a\right)\;.

The large deviation function $\Phi_{\omega_{c}}(J)$ is given by the Legendre transform of the CGF

\Phi_{\omega_{c}}(J)=\max_{\mu}\left[\mu\,J-\chi_{\omega_{c}}(\mu)\right]\;,

(32)

which is simply obtained via a saddle point approximation for large $t$ with $J=Q/t<\infty$ as $t\to\infty$ . Furthermore, as all the non-zero cumulants have the same scaling $O(t)$ , we also expect the central limit theorem to hold for the long-time typical fluctuations of the alternating charge. The corresponding distribution is Gaussian with fluctuations of order $\sqrt{t}$ around the average value

P_{\omega_{c},\alpha}(Q;t)\approx\frac{\displaystyle\exp\left(-\frac{(Q-\Big{\langle}Q_{\omega_{c};\alpha}(t)\Big{\rangle})^{2}}{2{\rm Var}(Q_{\omega_{c};\alpha}(t))}\right)}{\sqrt{2\pi{\rm Var}(Q_{\omega_{c};\alpha}(t))}}\;.

(33)

II.4 Incommensurate frequencies

We now consider the situation where the frequency $\omega_{d}$ of the drive and $\omega_{c}$ of the current are incommensurate. As the tilted operator ${\cal L}({\bf x},\nabla_{\bf x},t)$ is not periodic, we cannot use a Floquet expansion to describe the tilted propagator $G_{\mu}({\bf x},t|{\bf x}_{0})$ and it is not obvious how to compute the limit in equation (28). However, one can use Oseledets’ multiplicative ergodic theorem Oseledets [1968] to show, nevertheless, that this limit is well defined for any frequency $\omega_{c}$ . This implies in particular that all the cumulants scale at most linearly in time in this more generic case.

A naive way to compute the CGF in this case is to consider a series of frequencies $\omega_{c}(n)$ , commensurate with $\omega_{d}$ for any finite $n$ but that converges to the desired incommensurate value of $\omega_{c}$ as $n\to\infty$ . The CGF can be calculated for any $n$ using the technique described above for commensurate ratios, and in the limit $n\to\infty$ one can hope to estimate $\chi_{\mu}(\omega_{c})$ . However, this naive approach only works if the CGF $\chi_{\mu}(\omega_{c})$ is a continuous function of $\omega_{c}$ (at a fixed value of $\mu$ ), which we prove below quite surprisingly not to be the case. Therefore, at incommensurate frequency ratios the CGF can (to the best our knowledge) only be evaluated by taking the $t\to\infty$ limit in equation (28).

III Main results

Let us summarise our main findings and their consequences. Our first main result is that for a periodically driven system with driving frequency $\omega_{d}$ , the cumulant generating function (CGF) of the time-integrated alternating current $Q_{\omega_{c};\alpha}(t)$ , namely $\chi_{\mu}(\omega_{c})$ defined in in Eq. (28), is not a continuous function of the current’s frequency $\omega_{c}$ in the large time limit. The CGF varies smoothly for any current frequency $\omega_{c}$ incommensurate with the frequency $\omega_{d}$ of the drive but displays additional discontinuous spikes at commensurate frequencies. This feature can be explained by considering the behaviour of the rescaled cumulants ${\cal Q}_{p}^{\alpha}(\omega_{c})$ of order $p\geq 1$ as a function of the current frequency $\omega_{c}$ . The cumulants of even order, namely with $p=2r$ and $r\in\mathbb{N}^{*}$ , have a generically non-zero continuous background for irrational frequency ratios $\omega_{c}/\omega_{d}\notin\mathbb{Q}$ and they display discontinuous spikes at frequencies

\omega_{c}=\frac{n}{2z}\omega_{d}\;,\;\;n\in\mathbb{N}^{*}\;,\;\;1\leq z\in\mathbb{N}\leq r\;.

(34)

The cumulants of odd order, with $p=2r+1$ and $r\in\mathbb{N}$ , have zero continuous background for irrational ratios $\omega_{c}/\omega_{d}\notin\mathbb{Q}$ , but they do display discontinuous spikes at frequencies

\omega_{c}=\frac{n}{2z+1}\omega_{d}\;,\;\;n\in\mathbb{N}^{*}\;,\;\;0\leq z\in\mathbb{N}\leq r\;.

(35)

Note that the number of spikes increases with the order $p$ of the cumulant. However, the height of these spikes decreases with the order $p$ , and goes to zero at the $p\to\infty$ limit. A consequence of this result is that for any fixed $\omega_{c}$ which is incommensurate with $\omega_{d}$ , the CGF is symmetric $\chi_{\omega_{c}}(\mu)=\chi_{\omega_{c}}(-\mu)$ , implying an analogous symmetry of the PDF $P_{\omega_{c}}(Q,t)=P_{\omega_{c}}(-Q,t)$ in the large time limit. -

The discontiuous behavior of the cumulant generating function is a result of the long time limit: one clearly does not expect any discontinuity of the cumulant generating function at finite time. We find that at large, but finite time $t$ , the behavior of the cumulant generating function is qualitatively the same, but as expected the spikes in the spectrum of the cumulants are smeared over a typical frequency scale of order $\sim t^{-1}$ .

These results are illustrated by considering in details two exactly solvable examples and one example that can be calculated numerically. First, we consider an underdamped Brownian particle trapped in a periodically varying harmonic potential. Then, we obtain explicit results for a harmonically confined ”active” particle, i.e. subject to telegraphic noise. Lastly, we consider a time-periodic two-state system, where the cumulant generation function can be numerically evaluated.

IV Derivation of the cumulants

We consider in this section the cumulants of the alternating charge $Q_{\omega_{c};\alpha}(t)$ in the large time limit for a gaped system, i.e. $\real(\lambda_{1})>\lambda_{0}=0$ where $\lambda_{0},\lambda_{1}$ are the two lowest Floquet eigenvalues associated to the ”bare” evolution operator ${\cal L}({\bf x},\nabla_{\bf x},t)$ , i.e. in the absence of tilting $\mu=0$ (in this case the operator is periodic by definition and therefore the Floquet thery can be applied). As seen in section II, we expect that the rescaled cumulants

{\cal Q}_{p}^{\alpha}(\omega_{c})=\lim_{t\to\infty}\frac{\Big{\langle}Q_{\omega_{c};\alpha}^{p}(t)\Big{\rangle}_{\rm co}}{t}\;,

(36)

are of order $O(1)$ at most (we recall that a function of time $A(t)$ is defined to be of order $O(t^{\alpha})$ if $0<\lim_{t\to\infty}t^{-\alpha}A(t)<\infty$ ). In the following, we derive systematically for any value of $p$ its expression as a function of the parameters $\omega_{c},\omega_{d}$ and the Floquet spectrum $\{\lambda_{k},k\in\mathbb{N}\}$ . Before deriving this expression, we first use Eq. (15) to rewrite

	$\displaystyle Q_{\omega_{c};\alpha}(t)$	$\displaystyle=\left[x_{\alpha}(\tau)\,\cos(\omega_{c}\tau)\right]_{0}^{t}+q_{\omega_{c};\alpha}(t)$
	$\displaystyle q_{\omega_{c};\alpha}(t)$	$\displaystyle=\omega_{c}\int_{0}^{t}d\tau\,x_{\alpha}(\tau)\,\sin(\omega_{c}\tau)\;.$		(37)

For the systems that we are considering, we expect the average value and fluctuations of $x_{\alpha}(t)$ to remain of $O(1)$ as $t\to\infty$ . As shown in App. A, this simplifies the computation of the cumulants for any value of $p$ as

{\cal Q}_{p}^{\alpha}(\omega_{c})=\lim_{t\to\infty}\frac{\Big{\langle}q_{\omega_{c};\alpha}^{p}(t)\Big{\rangle}_{\rm co}}{t}\;.

(38)

From this identification, we already obtain that the system does not have a direct charge ${\cal Q}_{p}^{\alpha}(\omega_{c}=0)=0$ for any $p\geq 0$ . This is indeed expected, as non-zero cumulants of the charge would imply some fluctuations of the steady state current, which are incompatible with our assumption of simply-connected micro-state space. We note, however, that in a compact domain and a non-trivial topology, e.g. a particle on a ring, the empirical alternating and direct charges are not defined as per the second line of equation (15). Indeed, in that case, the fluctuation of the direct charge might be non-zero for such systems Busiello et al. [2018].

The cumulant of order $p$ of $q_{\omega_{c};\alpha}(t)$ reads

	$\displaystyle\Big{\langle}q_{\omega_{c};\alpha}^{p}(t)\Big{\rangle}_{\rm co}=\Big{\langle}\prod_{j=1}^{p}\int_{0}^{t}dt_{j}\,x_{\alpha}(t_{j})\,\omega_{c}\,\sin(\omega_{c}t_{j})\Big{\rangle}_{\rm co}$		(39)
	$\displaystyle=p!\,\omega_{c}^{p}\int_{0}^{t}dt_{p}\cdots\int_{0}^{t_{2}}dt_{1}\Big{\langle}\prod_{j=1}^{p}x_{\alpha}(t_{j})\Big{\rangle}_{\rm co}\prod_{j=1}^{p}\sin(\omega_{c}t_{j})\,,$

where the times $t\geq t_{p}\geq t_{p-1}\geq\cdots\geq t_{2}\geq t_{1}\geq 0$ in the second line are ordered. In this expression the $p$ -times connected correlation functions $\Big{\langle}\prod_{j=1}^{p}x_{\alpha}(t_{j})\Big{\rangle}_{\rm co}$ are conveniently expressed in terms of the lower order $n\leq p$ -times (disconnected) correlation functions $\Big{\langle}\prod_{j=1}^{n}x_{\alpha}(t_{j})\Big{\rangle}$ as follows Speed [1983]

	$\displaystyle\Big{\langle}\prod_{k=1}^{p}x_{\alpha}(t_{k})\Big{\rangle}_{\rm co}=$
	$\displaystyle\sum_{k=1}^{p}(k-1)!(-1)^{k+1}\sum_{\pi\in P_{k}(p)}\prod_{B\in\pi}\prod_{i=1}^{k}\Big{\langle}\prod_{j\in B_{i}}x(t_{j})\Big{\rangle}\;,$		(40)

where $\pi\in P_{k}(p)$ is an element of the groups of partitions of $\{1,\cdots,p\}$ into $k$ blocks $B_{i}$ ’s with $i=1,\cdots,k$ . We denote $n_{i}$ the number of elements in block $i$ of the partition $\pi$ , with $\sum_{i=1}^{k}n_{i}=p$ , and each element within a given block is ordered, i.e $B_{i}(1)<\cdots<B_{i}(n_{i})$ .

The $n$ -times (disconnected) correlation functions $\Big{\langle}\prod_{j=1}^{n}x_{\alpha}(t_{j})\Big{\rangle}$ can be computed explicitly by introducing the propagator between each of the times $t_{j}$ ’s, which results in

	$\displaystyle\Big{\langle}\prod_{j=1}^{p}x_{\alpha}(t_{j})\Big{\rangle}$		(41)
	$\displaystyle=\int d{\bf x}_{1}\cdots\int d{\bf x}_{p}\prod_{j=1}^{p}\left[x_{j,\alpha}\,G({\bf x}_{j},t_{j}\|{\bf x}_{j-1},t_{j-1})\right]$
	$\displaystyle=\sum_{l_{1},\cdots,l_{p}=0}^{\infty}\left[\prod_{j=1}^{p}e^{-\lambda_{l_{j}}(t_{j}-t_{j-1})}\right]C_{l_{1},\cdots,l_{n}}(t_{1},\cdots,t_{p})\;,$

where $t_{0}=0$ . The third line is obtained by inserting the Floquet expansion (5) of the propagator where $l_{1},\cdots,l_{p}$ refer to indices in the Floquet spectrum and the function

	$\displaystyle C_{l_{1},\cdots,l_{p}}(t_{1},\cdots,t_{n})=g_{l_{1}}({\bf x}_{0},0)\int d{\bf x}_{p}\,x_{p,\alpha}\,f_{l_{p}}({\bf x}_{p},t_{p})$
	$\displaystyle\times\prod_{j=1}^{p-1}\int d{\bf x}_{j}\,x_{j,\alpha}\,f_{l_{j}}({\bf x}_{j},t_{j})g_{l_{j+1}}({\bf x}_{j},t_{j})\;,$		(42)

is periodic with fundamental frequency $\omega_{d}$ in all its variables. Notice that as $g_{0}({\bf x},t)=1$ for all ${\bf x}$ and $t$ , the term $\int d{\bf x}_{j}\,x_{j,\alpha}f_{l_{j}}({\bf x}_{j},t_{j})g_{l_{j+1}}({\bf x}_{j},t_{j})$ connecting $l_{j}$ with $l_{j+1}$ effectively becomes independent of $l_{j+1}$ for $l_{j+1}=0$ . Thus this term simplifies as a product of smaller order correlation function

	$\displaystyle C_{l_{1},\cdots,l_{p}}(t_{1},\cdots,t_{n})=$		(43)
	$\displaystyle C_{l_{1},\cdots,l_{j}}(t_{1},\cdots,t_{j})\times C_{0,\cdots,l_{p}}(t_{j+1},\cdots,t_{n})\;,\;\;l_{j+1}=0\;.$

Note also that as $\lambda_{0}=0$ in equation (45) for $l_{j+1}=0$ , the product

	$\displaystyle\prod_{m=1}^{p}e^{-\lambda_{l_{m}}(t_{m}-t_{m-1})}=$		(44)
	$\displaystyle\prod_{m=1}^{j}e^{-\lambda_{l_{m}}(t_{m}-t_{m-1})}\prod_{n=j+2}^{p}e^{-\lambda_{l_{n}}(t_{n}-t_{n-1})}\;,\;\;l_{j+1}=0\;,$

effectively disconnecting the correlations before $t_{j+1}$ from the correlations after this time.

Using the above and inserting (41), (42) into equation (45), taking the limit where $t_{1}\gg\lambda_{1}^{-1}$ and $\delta t_{i}=t_{i}-t_{i-1}=O(1)$ for $i=2,\cdots,p$ , the connected correlation function simplifies to

	$\displaystyle\Big{\langle}\prod_{j=1}^{p}x_{\alpha}(t_{j})\Big{\rangle}_{\rm co}\approx\sum_{k=1}^{p}(k-1)!(-1)^{k+1}\sum_{\pi\in P_{k}(p)}\prod_{B\in\pi}$		(45)
	$\displaystyle\sum_{l_{1},\cdots,l_{p}=0}^{\infty}(1-\delta_{l_{p},0})\delta_{l_{1},0}\prod_{j=2}^{p}e^{-\mu_{j}(l_{2},\cdots,l_{p};\pi)\delta t_{j}}$
	$\displaystyle\times\prod_{i=1}^{k}\left[\tilde{C}_{l_{B_{i}(1)},\cdots,l_{B_{i}(n_{i})}}(t_{B_{i}(1)},\cdots,t_{B_{i}(n_{i})})\right]\;,$

where, making explicit use of relation (43), we have introduced the functions $\tilde{C}_{l_{1}}(t_{1})=\delta_{l_{1},0}C_{0}(t_{1})$ and for $p>1$

	$\displaystyle\tilde{C}_{l_{1},\cdots,l_{p}}(t_{1},\cdots,t_{p})=$		(46)
	$\displaystyle\delta_{l_{1},0}\prod_{k=2}^{p}(1-\delta_{l_{k},0})C_{l_{1},\cdots,l_{p}}(t_{1},\cdots,t_{p})\;.$

The values of the strictly positive rates $\mu_{j}(l_{2},\cdots,l_{p};\pi)>0$ for $2\leq j\leq p$ depend on the specific partition $\pi$ and read, supposing that $j\in B_{i}$ ,

	$\displaystyle\mu_{j}=\lambda_{l_{j}}+\lambda_{l_{j+1}}(1-\delta_{l_{j},0})$		(47)
	$\displaystyle+\sum_{r=j+1}^{p}\lambda_{l_{r}}\prod_{m\neq i}\left[\delta_{r\in B_{m}}\Theta(j-\max\{s\in B_{m},s<r\})\right]\;,$

with $\delta_{r\in B_{m}}=\sum_{s\in B_{m}}\delta_{r,s}$ . The values of $\lambda_{r}$ ’s for $r>j$ are included in $\mu_{j}$ if the times $t_{r}$ and $t_{j}$ do not belong to the same block and $t_{j}$ is larger than all the times $t_{s}<t_{r}$ belonging to the same block as $r$ . In particular, Eq. (45) shows that the $p$ -time connected correlation decays exponentially with the associated rate $\mu_{j}$ as soon as one of the time-difference $\delta t_{j}=t_{j}-t_{j-1}$ becomes large.

IV.1 Proof of the general result

We are now able to prove our main claim by computing the general expression of ${\cal Q}_{p}^{\alpha}(\omega_{c})$ . Changing variables in the integrals from the times $t_{1},\cdots,t_{p}$ into $\delta t_{2}=t_{2}-t_{1},\delta t_{3},\cdots,\delta t_{p}=t_{p}-t_{p-1}$ and $t_{p}$ , we can rewrite equation (39) as

	$\displaystyle\Big{\langle}q_{\omega_{c};\alpha}^{p}(t)\Big{\rangle}_{\rm co}=$		(48)
	$\displaystyle p!\,\omega_{c}^{p}\int_{0}^{t}dt_{p}\int_{0}^{t_{p}}d\delta t_{p}\cdots\int_{0}^{t_{p}-\sum_{k=3}^{p}\delta t_{k}}d\delta t_{2}$
	$\displaystyle\Big{\langle}\prod_{j=0}^{p-1}x_{\alpha}\left(t_{p}-\sum_{k=0}^{j-1}\delta t_{p-k}\right)\Big{\rangle}_{\rm co}$
	$\displaystyle\prod_{j=0}^{p-1}\sin\left(\omega_{c}\left(t_{p}-\sum_{k=0}^{j-1}\delta t_{p-k}\right)\right)\;.$

It is clear that the variable $t_{p}$ is of order $O(t)$ in the large $t$ limit. The $p$ -times connected correlation term in the integrand $\Big{\langle}\prod_{j=1}^{p}x_{\alpha}\left(t_{j}\right)\Big{\rangle}_{\rm co}=\Big{\langle}\prod_{j=0}^{p-1}x_{\alpha}\left(t_{p}-\sum_{k=0}^{j-1}\delta t_{p-k}\right)\Big{\rangle}_{\rm co}$ is given in Eq. (45). This function does not decay with $t_{p}$ for fixed values of the $\delta t_{j}$ ’s but rather is a periodic function of this variable. On the other hand, we have seen that this function decays exponentially with all the variables $\delta t_{j}$ ’s for $2\leq j\leq p$ . We thus expect that the dominating contribution to the integral will come from $\delta t_{j}=O(1)$ while $t_{p}=O(t)$ in the large $t$ limit and we may safely replace the upper bounds $t_{p}-\sum_{k=j+1}^{p}\delta t_{k}$ in the integral over $\delta t_{j}$ by $+\infty$ in this limit.

Introducing the Fourier expansion for the coefficients

C_{l_{1},\cdots,l_{p}}(t_{1},\cdots,t_{p})=\sum_{k_{1},\cdots,k_{p}=-\infty}^{\infty}C_{l_{1},\cdots,l_{p}}^{k_{1},\cdots,k_{p}}e^{i\sum_{j=1}^{n}k_{j}\omega_{d}t}\;,

(49)

as well as the identity

\sin(\omega_{c}t)=\sum_{\sigma=\pm 1}\frac{(-\sigma)}{2i}e^{-i\sigma\omega_{c}t}\;,

(50)

one can express equation (48) in the long time limit as

	$\displaystyle\Big{\langle}q_{\omega_{c};\alpha}^{p}(t)\Big{\rangle}_{\rm co}\approx$	$\displaystyle p!\,\omega_{c}^{p}\sum_{k=1}^{p}(k-1)!(-1)^{k+1}\sum_{\pi\in P_{k}(p)}\prod_{B\in\pi}\sum_{l_{1},\cdots,l_{p}=0}^{\infty}(1-\delta_{l_{p},0})\delta_{l_{1},0}\sum_{k_{1},\cdots,k_{p}=-\infty}^{\infty}\sum_{\sigma_{1},\cdots,\sigma_{p}=\pm 1}\prod_{m=1}^{p}\left(\frac{-\sigma_{m}}{2i}\right)$		(51)
		$\displaystyle\times C_{l_{1},\cdots,l_{p}}^{k_{1},\cdots,k_{p}}\int_{0}^{t}dt_{p}\int_{0}^{\infty}d\delta t_{p}\cdots\int_{0}^{\infty}d\delta t_{2}\prod_{j=2}^{p}e^{-\mu_{j}(l_{2},\cdots,l_{p};\pi)\delta t_{j}}\prod_{m=1}^{p}e^{i(k_{m}\omega_{d}-\omega_{c}\sigma_{m})\left[t_{p}-\sum_{n=0}^{m-1}\delta t_{p-n}\right]}\;.$

In particular, for fixed values of $\{k_{1},\cdots,k_{p}\}$ and $\{\sigma_{1},\cdots,\sigma_{p}\}$ , one can compute explicitly the integral over the variable $t_{p}$ , which is the only term still depending on $t$ in the limit $t\gg\lambda_{1}^{-1}$ . It yields

	$\displaystyle\int_{0}^{t}dt_{p}\;e^{i\Omega_{p}t_{p}}$	$\displaystyle=e^{i\frac{\Omega_{p}t}{2}}\frac{2}{\Omega_{p}}\sin\left(\frac{\Omega_{p}t}{2}\right)\;,$		(52)
	$\displaystyle\Omega_{p}$	$\displaystyle=\omega_{d}\sum_{j=1}^{p}k_{j}-\omega_{c}\sum_{j=1}^{p}\sigma_{j}\;.$		(53)

In particular, it is only of order $t$ in the special case where $\Omega_{p}=0$ , giving a resonance condition for each of the cumulants. The only terms that contribute to $\mathcal{Q}_{p}(\omega_{c})$ , are by Eq.(38) these specific resonant contributions which scale as $O(t)$ , namely for which the right hand side of Eq.(53) vanishes. There are several important consequences for this observation. Let us therefore focus on the equation

\displaystyle\omega_{d}\sum_{j=1}^{p}k_{j}-\omega_{c}\sum_{j=1}^{p}\sigma_{j}=0\;.

(54)

This equality holds when both sums are zero independently

	$\displaystyle\sum_{j=1}^{p}k_{j}$	$\displaystyle=0$
	$\displaystyle\sum_{j=1}^{p}\sigma_{j}$	$\displaystyle=0\;$		(55)

or when $\omega_{c}$ and $\omega_{d}$ achieve some commensurate values

\displaystyle\frac{\omega_{c}}{\omega_{d}}=\frac{\sum_{j=1}^{p}k_{j}}{\sum_{j=1}^{p}\sigma_{j}}\;.

(56)

For incommensurate $\omega_{d}$ and $\omega_{c}$ , Eq.(56) cannot hold, since both the numerator and denominator in the right hand side of the equation are integers. Therefore, the only valid resonance condition corresponds to the case where both quantities in Eq.(55) are zero. However, the second line of (55) cannot hold for odd $p$ , since $\sigma_{j}=\pm 1$ and in particular the second line of Eq. (55) cannot vanish. Thus, for incommensurate frequency ratios all the odd cumulants are exactly zero. For even values of $p$ , it is possible that each of the sums above is zero separately, and therefore generically even cumulants do not vanish.

Next, let us consider commensurate frequencies: in this case, there are two options for condition Eq.(53) to hold: either Eq.(55) holds, or instead Eq.(56) can hold. For even values of $p$ , the latter case is possible whenever the relations between $\omega_{c}$ and $\omega_{d}$ is of the form

\omega_{c}=\frac{n}{2z}\omega_{d}\;,\;\;n\in\mathbb{N}^{*}\;,\;\;1\leq z\in\mathbb{N}\leq\frac{p}{2}\;,

(57)

since $\sum_{j}k_{j}$ can take any integer value, whereas $\sum_{j}\sigma_{j}$ is bounded between $-p$ and $p$ . Similarly, for odd values of $p$ , condition Eq.(53) can hold when

\omega_{c}=\frac{n}{2z+1}\omega_{d}\;,\;\;n\in\mathbb{N}^{*}\;,\;\;0\leq z\in\mathbb{N}\leq\frac{p-1}{2}\;.

(58)

In both the even and odd cases we expect contributions to the cumulants in addition to those that appear in the non-commensurate cases. Therefore, generically they have discontinuities at these specific commensurate frequencies.

At this point, one can rightfully suspect that the cumulant generating function $\chi_{\mu}(\omega_{c})$ , which holds the discontinuity of every cumulant, is for fixed $\mu$ an everywhere discontinuous function of the alternating charge’s frequency $\omega_{c}$ . However, we next show that it is a continuous function over the frequencies $\omega_{c}$ incommensurate with $\omega_{d}$ , similarly to the Thomae’s function Beanland et al. [2009]. To this end, we next examine the structure of ${\cal Q}_{p}^{\alpha}(\omega_{c})$ as a function of $p$ .

IV.2 General expression of ${\cal Q}_{p}^{\alpha}(\omega_{c})$

Using the tools presented above it is possible to obtain the general expression for the cumulant of arbitrary order $p$ . It reads

	$\displaystyle{\cal Q}_{p}^{\alpha}(\omega_{c})=$	$\displaystyle p!\sum_{k=1}^{p}(k-1)!(-1)^{k+1}\sum_{\pi\in P_{k}(p)}\prod_{B\in\pi}\sum_{k_{1},\cdots,k_{p}=-\infty}^{\infty}\sum_{\sigma_{1},\cdots,\sigma_{p}=\pm 1}\sum_{l_{1},\cdots,l_{p}=0}^{\infty}(1-\delta_{l_{p},0})\delta_{l_{1},0}$		(59)
		$\displaystyle\times\prod_{j=1}^{p}\left(\frac{-\sigma_{j}}{2i}\right)\frac{\displaystyle\prod_{i=1}^{k}\left[\delta_{l_{B_{i}(1)},0}\prod_{k=2}^{n_{i}}(1-\delta_{l_{B_{i}(k)},0})C_{l_{B_{i}(1)},\cdots,l_{B_{i}(n_{i})}}^{k_{B_{i}(1)},\cdots,k_{B_{i}(n_{i})}}\right]}{\displaystyle\prod_{j=2}^{p}\left[\mu_{j}(l_{2},\cdots,l_{p};\pi)+i\sum_{l=j}^{p}\left(\sigma_{l}\omega_{c}-k_{l}\omega_{d}\right)\right]}\omega_{c}^{p}\,\delta_{\sum_{j=1}^{p}\left(\omega_{d}k_{j}-\omega_{c}\sigma_{j}\right),0}\;.$

The Kronecker delta function in the second line of this expression ensures that the resonance condition is fulfilled. For incommensurate frequencies, this happens only when Eq. (55) holds but for commensurate ratio, this also happens when Eq. (56) holds.

We now explore, for a fixed $\mu$ , the continuity of the cumulant generating function

\chi_{\mu}(\omega_{c})=\sum_{p=1}^{\infty}\frac{\mu^{p}}{p!}{\cal Q}_{p}^{\alpha}(\omega_{c})\;.

(60)

As discussed above, each term of this series has specific discontinuities as detailed in Eqs. (57) and (58). For a sequence of rational ratios that convergence to an irrational ratio, namely a sequence of $(n_{i},m_{i})$ and a corresponding $\omega_{c}^{i}$

\omega_{c}^{i}=\omega_{d}\frac{n_{i}}{m_{i}}\;,

(61)

such that $a=\lim_{i\to\infty}n_{i}/m_{i}$ is an irrational number, the resonance condition in Eq. (56) imposes that the cumulants which display a discontinuity must be of increasing order $p_{i}$ . To this end, let us examine carefully the denominator of equation (59). For any a value of $p_{i}$ , the resonance condition imposes

\frac{\sum_{j=1}^{p_{i}}k_{j}}{\sum_{j=1}^{p_{i}}\sigma_{j}}=\frac{n_{i}}{m_{i}}\;.

(62)

For a random strings $\{\sigma_{j},j=1,\cdots,p\}$ and $\{k_{j},j=1,\cdots,p\}$ with the constraint of fulfilling the latter equation, the resonance condition imposes the following equality:

\left|\sum_{l=j}^{p_{i}}\left(\sigma_{l}\omega_{c}-k_{l}\omega_{d}\right)\right|=\left|\sum_{l=1}^{j-1}\left(\sigma_{l}\omega_{c}-k_{l}\omega_{d}\right)\right|\;.

(63)

The left hand side of this equation scales, in the large $j$ limit, as $O(\sqrt{j})$ . Thus we expect the product appearing in the denominator of Eq. (59)

\prod_{j=2}^{p_{i}}\left[\mu_{j}+i\sum_{l=j}^{p_{i}}\left(\sigma_{l}\omega_{c}-k_{l}\omega_{d}\right)\right]

(64)

to grow in modulus rather rapidly as a function of $p_{i}$ (as we expect $\mu_{j}=O(1)$ , this term will be of order $\sqrt{p_{i}!}$ ). Therefore we expect that ${\cal Q}_{p_{i}}^{\alpha}(\omega_{c})/(p_{i}!)$ goes to zero as $p_{i}\to\infty$ . The discontinuous part of the CGF at frequency $\omega_{c}^{i}$ are provided by the terms ${\cal Q}_{p}^{\alpha}(\omega_{c})/(p!)$ of order $p>p_{i}$ and these rescaled cumulants are getting smaller and smaller as $p_{i}$ is increased, we expect that in the limit where $i\to\infty$ , the discontinuities vanish completely from the CGF.

To get some insight on equation (59), we next use it to extract the explicit expressions for the first two cumulants, which are the main objects of interest in most cases.

IV.3 Average alternating current

Rather than simply plugging $p=1$ in Eq. (59), it is both easy and insightful to rapidly re-derive the result in this particular case at finite but large time $t$ . To this end, we first compute the finite time first cumulant $\Big{\langle}q_{\omega_{c};\alpha}(t)\Big{\rangle}$ . It reads

\Big{\langle}q_{\omega_{c};\alpha}(t)\Big{\rangle}=\omega_{c}\sum_{l_{1}=0}^{\infty}\int_{0}^{t}dt_{1}\sin(\omega_{c}t_{1})e^{-\lambda_{1}t_{1}}C_{l_{1}}(t_{1})\;,

(65)

where we have used the Floquet expansion of the propagator, to obtain

\Big{\langle}x_{\alpha}(t)\Big{\rangle}=\sum_{l_{1}=0}^{\infty}e^{-\lambda_{1}t}C_{l_{1}}(t)\;,

(66)

and we remind that the function $C_{l_{1}}(t)$ is a periodic function of $t$ with fundamental frequency $\omega_{d}$ . We introduce the Fourier expansion

\displaystyle C_{l_{1}}(t)

\displaystyle=\sum_{k=-\infty}^{\infty}C_{l_{1}}^{k}e^{ik\omega_{d}t}\;.

(67)

Note that $\Big{\langle}x_{\alpha}(t)\Big{\rangle}$ is a real function of $t$ such that if the Floquet eigenvalue $\lambda_{l_{1}}$ is real, the associated function $C_{l_{1}}(t)$ must also be real and $C_{l_{1}}^{-n}={C_{l_{1}}^{n}}^{*}$ (we denote by $z^{*}$ here and in the following the complex conjugate of $z$ ). If the Floquet eigenvalue $\lambda_{l_{1}}$ is complex, it will come in a complex conjugate pair, say $\lambda_{l_{1}}^{*}=\lambda_{l_{1}+1}$ and one must have instead that $C_{l_{1}}(t)=C_{l_{1}+1}(t)^{*}$ , which yields $C_{l_{1}}^{-n}={C_{l_{1}+1}^{n}}^{*}$ . Inserting this expansion into Eq. (65), we obtain in the regime where $t\gg\lambda_{1}^{-1}$

\displaystyle\Big{\langle}q_{\omega_{c};\alpha}(t)\Big{\rangle}\approx

\displaystyle\sum_{k=-\infty}^{\infty}\frac{\omega_{c}\,C_{0}^{k}}{\omega_{c}^{2}-k^{2}\omega_{d}^{2}}\left[\omega_{c}(1-e^{ik\omega_{d}t}\cos(\omega_{c}t))+ik\omega_{d}e^{ik\omega_{d}t}\sin(\omega_{c}t)\right]+\sum_{k=-\infty}^{\infty}\sum_{l_{1}=1}^{\infty}\frac{\omega_{c}^{2}\,C_{l_{1}}^{k}}{\omega_{c}^{2}+(\lambda_{l_{1}}-ik\omega_{d})^{2}}\;.

(68)

Using that ${C_{0}^{k}}^{*}=C_{0}^{-k}$ since $\lambda_{0}=0$ , in the double scaling limit $t\to\infty$ , $\omega_{c}-n\omega_{d}\to 0$ with $\phi=(\omega_{c}-n\omega_{d})t=O(1)$ fixed, the rescaled alternating charge takes the following scaling form

	$\displaystyle\frac{\Big{\langle}q_{\omega_{c};\alpha}(t)\Big{\rangle}}{t}$	$\displaystyle\approx-\omega_{c}\,{\cal J}_{n}((\omega_{c}-n\omega_{d})t)$		(69)
	$\displaystyle{\cal J}_{n}(\phi)$	$\displaystyle=-\real\left[C_{0}^{n}\right]\frac{2\sin^{2}\left(\phi/2\right)}{\phi}+\imaginary\left[C_{0}^{n}\right]\frac{\sin\left(\phi\right)}{\phi}\;.$

We are mainly interested in the limit ${\cal Q}_{1}^{\alpha}(\omega_{c})=\lim_{t\to\infty}\Big{\langle}q_{\omega_{c};\alpha}(t)\Big{\rangle}/t$ , which can simply be obtained from this result as

	$\displaystyle{\cal Q}_{1}^{\alpha}(\omega_{c})$	$\displaystyle=-\omega_{c}\sum_{n=0}^{\infty}{\cal J}_{n}(0)\delta_{n\omega_{d},\omega_{c}}$
		$\displaystyle=-\omega_{c}\sum_{n=0}^{\infty}\imaginary\left[C_{0}^{n}\right]\delta_{\omega_{c},n\omega_{d}}\;.$		(70)

This structure is clearly discontinuous in the $t\to\infty$ limit as any term for $\omega_{c}\neq n\omega_{d}$ decays as $1/t$ . The average alternating charge ${\cal Q}_{1}^{\alpha}(\omega_{c})$ is zero for any frequency $\omega_{c}\neq n\omega_{d}$ that is not a harmonic of the driving frequency $\omega_{d}$ as one could have naively guessed: there is no average probability current in periodically driven systems, except at frequencies which are an integer multiple of the driving frequency.

IV.4 Variance of the alternating charge

We next turn to the computation of the variance of the alternating charge. In order to compute this object, we first compute the connected two-time correlation function for $t_{2}>t_{1}\gg\lambda_{1}^{-1}$

\displaystyle\Big{\langle}x_{\alpha}(t_{1})x_{\alpha}(t_{2})\Big{\rangle}_{\rm co}\approx\sum_{l_{2}=1}^{\infty}e^{-\lambda_{l_{2}}(t_{2}-t_{1})}C_{0,l_{2}}(t_{1},t_{2})\;,

(71)

where we have used that $C_{l_{1},l_{2}=0}(t_{1},t_{2})=C_{l_{1}}(t_{1})C_{0}(t_{2})$ . The expression of the second cumulant is again particularly simple as only the identity partition of $\{1,2\}$ yields a non-zero contribution in the connected $2$ -times correlation function. It reads

	$\displaystyle{\cal Q}_{2}^{\alpha}(\omega_{c})=$		(72)
	$\displaystyle\sum_{k_{1},k_{2}=-\infty}^{\infty}\sum_{\sigma_{1},\sigma_{2}=\pm 1}\sum_{l_{2}=1}^{\infty}\frac{-\omega_{c}^{2}\,\sigma_{1}\sigma_{2}\,C_{0,l_{2}}^{k_{1},k_{2}}\,\delta_{(k_{1}+k_{2})\omega_{d},(\sigma_{1}+\sigma_{2})\omega_{c}}}{2[\lambda_{l_{2}}+i(\sigma_{2}\omega_{c}-k_{2}\omega_{d})]}$
	$\displaystyle={\cal Q}_{2}^{{\rm b},\alpha}(\omega_{c})-\sum_{n=1}^{\infty}\sum_{l_{2}=1}^{\infty}\sum_{k=-\infty}^{\infty}\sum_{\sigma=\pm 1}\frac{\omega_{c}^{2}\,C_{0,l_{2}}^{\sigma n-k,k}\,\delta_{2\omega_{c},n\omega_{d}}}{2[\lambda_{l_{2}}+i(\sigma\omega_{c}-k\omega_{d})]}\,,$

where the continuous background ${\cal Q}_{2}^{{\rm b},\alpha}(\omega_{c})$ reads

{\cal Q}_{2}^{{\rm b},\alpha}(\omega_{c})=\sum_{l_{2}=1}^{\infty}\sum_{k=-\infty}^{\infty}\frac{\omega_{c}^{2}(\lambda_{l_{2}}-ik\omega_{d})\,C_{0,l_{2}}^{-k,k}}{(\lambda_{l_{2}}-ik\omega_{d})^{2}+\omega_{c}^{2}}\,.

(73)

In contrast to ${\cal Q}_{1}^{\alpha}(\omega_{c})$ , which is zero except for frequencies that satisfy $\omega_{c}=n\omega_{d}$ , the rescaled variance ${\cal Q}_{2}^{\alpha}(\omega_{c})$ is generically non-zero for any frequency $\omega_{c}$ and equal to ${\cal Q}_{2}^{{\rm b},\alpha}(\omega_{c})$ away from the integer and half-integer multiple of $\omega_{d}$ . Therefore, there are current fluctuations even at frequency with no average current, as expected.

It is instructive to compare these results with the power spectral density calculated for the stochastic resonance setup Shneidman et al. [1994]; Nikitin et al. [2007, 2005]. There, the peaks appear only at integer multiples of $\omega_{d}$ , whereas in our setup they appear at both integer and half integer multiples of $\omega_{d}$ .

The expressions for higher order cumulants can be obtained more explicitly but are quite cumbersome. In Appendices B and C, we respectively give explicit expressions for the third and fourth cumulant of the alternating charge.

IV.5 Large $\omega_{c}$ limit of the CGF

In the case $\omega_{c}\gg\omega_{d}$ , a simplified expression for the tilted operator ${\cal L}_{\mu}({\bf x},\nabla_{\bf x},t)$ can be obtained. The bare dynamics that evolves over the period $T_{d}$ is not fast enough to change over a period $T_{c}\ll T_{d}$ , such that the terms depending on $\omega_{c}$ can be replaced in practice by their cycle average value over the period $T_{c}$ , i.e.

	$\displaystyle{\cal L}_{\mu}({\bf x},\nabla_{\bf x},t)G_{\mu}({\bf x},t\|{\bf x}_{0})\approx\overline{{\cal L}_{\mu}({\bf x},\nabla_{\bf x},t)}^{\omega_{c}}G_{\mu}({\bf x},t\|{\bf x}_{0})\;,$
	$\displaystyle\overline{{\cal L}_{\mu}({\bf x},\nabla_{\bf x},t)}^{\omega_{c}}={\cal L}({\bf x},\nabla_{\bf x},t)+\frac{\mu^{2}}{2}D_{\alpha\alpha}({\bf x},t)\;,$		(74)

where $\overline{A}^{\omega_{c}}=T_{c}^{-1}\int_{0}^{T_{c}}A(t)dt$ . In the specific case where $D_{\alpha\alpha}({\bf x},t)=D_{\alpha\alpha}(t)$ is independent of ${\bf x}$ , one can show that the tilted propagator can be simply expressed as a function of the bare propagator as

G_{\mu}({\bf x},t)\approx G({\bf x},t|{\bf x}_{0},0)e^{\frac{\mu^{2}}{2}\int_{0}^{t}D_{\alpha\alpha}(\tau)d\tau}\;.

(75)

Form this result, one obtains that as $\omega_{c}\to\infty$

	$\displaystyle\chi_{\mu}(\omega_{c})$	$\displaystyle=\lim_{t\to\infty}\frac{1}{t}\ln\int d{\bf x}\,G_{\mu}({\bf x},t)\approx\frac{\mu^{2}}{2}\overline{D_{\alpha\alpha}}^{\omega_{d}}\;,$		(76)
	$\displaystyle{\cal Q}_{p}(\omega_{c})$	$\displaystyle=\left.\partial_{\mu}^{p}\chi_{\mu}(\omega_{c})\right\|_{\mu=0}\approx\delta_{p,2}\overline{D_{\alpha\alpha}}^{\omega_{d}}\;,$		(77)

where $\overline{D_{\alpha\alpha}}^{\omega_{d}}=T_{d}^{-1}\int_{0}^{T_{d}}D_{\alpha\alpha}(\tau)d\tau$ . In particular, it implies that both the typical and atypical fluctuations of the alternating charge are Gaussian in the $\omega_{c}\to\infty$ limit.

V Specific examples

We have derived general expressions for the first two cumulants of the alternating charge. Next, we demonstrate these results in specific examples, which are commonly used in statistical mechanics and where a full analytical computation is possible.

V.1 Underdamped Ornstein Uhlenbeck

As a first example, we consider the following dynamics

	$\displaystyle\dot{v}(t)$	$\displaystyle=\dot{\mu}(t)-\Gamma(t)\left[v(t)-\mu(t)+\Omega(t)x(t)-\sqrt{2D(t)}\eta(t)\right]\;,$
	$\displaystyle\dot{x}(t)$	$\displaystyle=v(t)\;,$		(78)

where $\eta(t)$ is a Gaussian white noise with zero mean $\Big{\langle}\eta(t)\Big{\rangle}=0$ and variance $\Big{\langle}\eta(t_{1})\eta(t_{2})\Big{\rangle}=\delta(t_{1}-t_{2})$ . The damping $\Gamma(t)$ , trapping frequency $\Omega(t)$ , diffusion coefficient $D(t)$ and drift $\mu(t)$ are all time-periodic functions with fundamental angular frequency $\omega_{d}$ . Keeping track of both $(x,v)$ at time $t$ , the system is Markovian and the general theory described in the sections above applies to this case.

By linearity of the above equation, both the position $x(t)$ and the speed $v(t)$ are expressed as (infinite) linear combinations of the Gaussian random variables $\eta(\tau)$ for $\tau\in[0,t]$ and thus also have Gaussian fluctuations. The same reasoning applies to the empirical alternating charge $Q_{\omega_{c};x}(t)$ , which is expressed as (infinite) sums of $x(\tau)$ for $\tau\in[0,t]$ . It yields that all the cumulants beyond the variance are zero, i.e.

	$\displaystyle{\cal Q}_{p}^{x}(\omega_{c})=0\;,\;\;p>2\;,$		(79)
	$\displaystyle\chi_{\mu}(\omega_{c})=\mu{\cal Q}_{1}^{x}(\omega_{c})+\frac{\mu^{2}}{2}{\cal Q}_{2}^{x}(\omega_{c})\;.$		(80)

A direct consequence is that the large deviation function takes the quadratic form

\Phi_{\omega_{c}}(J)=\frac{(J-{\cal Q}_{1}^{x}(\omega_{c}))^{2}}{2{\cal Q}_{2}^{x}(\omega_{c})}\;.

(81)

The joint bare propagator $G(v,x,t|v_{0},x_{0},t_{0})$ of the speed $v$ and position $x$ satisfies the equation

	$\displaystyle\partial_{t}G=$	$\displaystyle-(v+\mu(t))\partial_{x}G+\Gamma(t)\partial_{v}\left[(v+\Omega(t)x)G\right]$		(82)
		$\displaystyle+\Gamma(t)D(t)\partial_{v}^{2}G\;.$

We consider a simple example with constant damping $\Gamma(t)=\Gamma_{0}$ and trapping frequency $\Omega(t)=\Omega_{0}$ and introduce the Fourier series

\mu(t)=\sum_{k=-\infty}^{\infty}\mu_{k}\,e^{ik\omega_{d}t}\;,\;\;D(t)=\sum_{k=-\infty}^{\infty}D_{k}\,e^{ik\omega_{d}t}\;,

(83)

which correspond to periodic forcing of the drift and the diffusion coefficient. We suppose that the functions $\mu(t)$ and $D(t)$ are real, such that for all $k$ , $\mu_{k}=\mu_{R,|k|}+i\,\mbox{sgn}\left[k\right]\mu_{I,|k|}$ and similarly $D_{k}=D_{R,|k|}+i\,\mbox{sgn}\left[k\right]D_{I,|k|}$ . In this case, the position at time $t$ , starting with the initial condition $(x(0),v(0))=(0,0)$ is given by

	$\displaystyle x(t)=$	$\displaystyle\sqrt{\frac{\Gamma_{0}}{\Gamma_{0}-4\Omega_{0}}}\sum_{\alpha=\pm 1}\alpha$		(84)
		$\displaystyle\times\int_{0}^{t}d\tau\,e^{-\frac{(t-\tau)}{2}\left[\Gamma_{0}-\alpha\sqrt{\Gamma_{0}^{2}-4\Gamma_{0}\Omega_{0}}\right]}\,\left[\eta(\tau)+\mu(\tau)\right]\;.$

This expression allows to identify the Floquet spectrum $\lambda_{0}=0\leq\lambda_{1}\leq\lambda_{2}$ and $\lambda_{k}=+\infty$ for $k\geq 3$ with

	$\displaystyle\lambda_{1}$	$\displaystyle=\frac{\Gamma_{0}}{2}-\frac{\sqrt{\Gamma_{0}^{2}-4\Gamma_{0}\Omega_{0}}}{2}\;,$
	$\displaystyle\lambda_{2}$	$\displaystyle=\frac{\Gamma_{0}}{2}+\frac{\sqrt{\Gamma_{0}^{2}-4\Gamma_{0}\Omega_{0}}}{2}\;.$		(85)

Using the expression of the position, one can show that

	$\displaystyle\Big{\langle}x(t)\Big{\rangle}=$	$\displaystyle\sqrt{\frac{4\Gamma_{0}}{\Gamma_{0}-4\Omega_{0}}}$		(86)
		$\displaystyle\times\int_{0}^{t}d\tau\,e^{-\frac{\Gamma_{0}}{2}\tau}\sinh\left(\sqrt{\Gamma_{0}^{2}-4\Gamma_{0}\Omega_{0}}\frac{\tau}{2}\right)\,\mu(t-\tau)\;.$

Inserting the Fourier series of $\mu(t)$ and taking the long-time limit $t\gg\lambda_{1}^{-1}$ , one can identify this expression with

\displaystyle\Big{\langle}x(t)\Big{\rangle}\approx

\displaystyle\sum_{k=-\infty}^{\infty}C_{0}^{k}e^{ikt}\;,

(87)

which allows to obtain the expression of the coefficients

C_{0}^{k}=\frac{\Gamma_{0}\mu_{k}}{\Gamma_{0}(\Omega_{0}+ik\omega_{d})-k^{2}\omega_{d}^{2}}\;.

(88)

Using this expression, we obtain the expression of the first cumulant

	$\displaystyle{\cal Q}_{1}^{x}(\omega_{c})$	$\displaystyle=\sum_{n=1}^{\infty}Q_{1,n}\delta_{n\omega_{d},\omega_{c}}\;,$		(89)
	$\displaystyle Q_{1,n}=$	$\displaystyle\frac{\mu_{R,n}\Gamma_{0}^{2}n^{2}\omega_{d}^{2}+\mu_{I,n}\Gamma_{0}n\omega_{d}(n^{2}\omega_{d}^{2}-\Gamma_{0}\Omega_{0})}{\Gamma_{0}^{2}\Omega_{0}^{2}+\Gamma_{0}(\Gamma_{0}-2\Omega_{0})n^{2}\omega_{d}^{2}+n^{4}\omega_{d}^{4}}\;.$

Refer to caption — Figure 1: Comparison between the time-averaged alternating current $\Big{\langle}Q_{\omega_{c}}(t)\Big{\rangle}/t$ (in blue) computed numerically for a measurement time $t=10^{2}$ for a periodically driven underdamped Brownian particle with $\Gamma_{0}=2$ , $\Omega_{0}=1$ , $\omega_{d}=1$ and for $\mu(t)=3/2\cos(2\omega_{d}t)$ and the analytical expression at finite time (in orange) given by inserting the coefficient $C_{0}^{k}$ given in (88) into Eq. (69).

One can similarly compute the two-time correlation function

	$\displaystyle\Big{\langle}x(t_{1})x(t_{2})\Big{\rangle}_{\rm co}=\Big{\langle}x(t_{1})x(t_{2})\Big{\rangle}-\Big{\langle}x(t_{1})\Big{\rangle}\Big{\langle}x(t_{2})\Big{\rangle}$		(90)
	$\displaystyle=\sum_{\alpha_{1},\alpha_{2}=\pm 1}\frac{2\Gamma_{0}\alpha_{1}\alpha_{2}}{\Gamma_{0}-4\Omega_{0}}e^{-\frac{(t_{2}-t_{1})}{2}\left[\Gamma_{0}-\alpha_{2}\sqrt{\Gamma_{0}^{2}-4\Gamma_{0}\Omega_{0}}\right]}$
	$\displaystyle\times\int_{0}^{t_{1}}d\tau\,e^{-\frac{(t_{1}-\tau)}{2}\left[2\Gamma_{0}-(\alpha_{1}+\alpha_{2})\sqrt{\Gamma_{0}^{2}-4\Gamma_{0}\Omega_{0}}\right]}D(\tau)\;.$

From this expression, it is possible to identify the coefficients

C_{0,l_{2}}^{k_{1},k_{2}}=\frac{2\Gamma_{0}^{3/2}D_{k_{1}}}{\sqrt{\Gamma_{0}-4\Omega_{0}}}\frac{(\delta_{l_{2},1}-\delta_{l_{2},2})\delta_{k_{2},0}}{[\Gamma_{0}+ik_{1}\omega_{d}][2\lambda_{l_{2}}+ik_{1}\omega_{d}]}\;.

(91)

One can check that for $\Gamma_{0}<4\Omega_{0}$ , we have that $\imaginary[\lambda_{1}]=\imaginary[\lambda_{2}]=0$ and $\imaginary[C_{0,l_{2}}^{0,0}]=0$ while for $\Gamma_{0}>4\Omega_{0}$ we have $\lambda_{1}=\lambda_{2}^{*}$ such that $C_{0,1}^{0,0}={C_{0,2}^{0,0}}^{*}$ . After simplifications, one can obtain explicitly the value of the rescaled variance of the alternating charge by inserting the expression of the coefficients (91) and of the spectrum (85) into (72)

	$\displaystyle{\cal Q}_{2}^{x}(\omega_{c})=\frac{\Gamma_{0}^{2}\omega_{c}^{2}D_{0}}{\Gamma_{0}^{2}\Omega_{0}^{2}+\Gamma_{0}(\Gamma_{0}-2\Omega_{0})\omega_{c}^{2}+\omega_{c}^{4}}$		(92)
	$\displaystyle-\sum_{n=1}^{\infty}\frac{\Gamma_{0}^{2}\omega_{c}^{2}D_{R,n}(\Gamma_{0}^{2}\Omega_{0}^{2}-\Gamma_{0}(\Gamma_{0}+2\Omega_{0})\omega_{c}^{2}+\omega_{c}^{4})}{(\Gamma_{0}^{2}\Omega_{0}^{2}+\Gamma_{0}(\Gamma_{0}-2\Omega_{0})\omega_{c}^{2}+\omega_{c}^{4})^{2}}\delta_{n\omega_{d},2\omega_{c}}$
	$\displaystyle-\sum_{n=1}^{\infty}\frac{2\Gamma_{0}^{3}\omega_{c}^{3}D_{I,n}(\Gamma_{0}\Omega_{0}-\omega_{c}^{2})}{(\Gamma_{0}^{2}\Omega_{0}^{2}+\Gamma_{0}(\Gamma_{0}-2\Omega_{0})\omega_{c}^{2}+\omega_{c}^{4})^{2}}\delta_{n\omega_{d},2\omega_{c}}\;.$

Note that the continuous background,

{\cal Q}_{2}^{{\rm b},x}(\omega_{c})=\frac{\Gamma_{0}^{2}\omega_{c}^{2}D_{0}}{\Gamma_{0}^{2}\Omega_{0}^{2}+\Gamma_{0}(\Gamma_{0}-2\Omega_{0})\omega_{c}^{2}+\omega_{c}^{4}}\;,

(93)

is identical to the variance of the fluctuating charge for a system with constant diffusion coefficient $D_{0}$ , i.e. in the absence of any periodic drive. The only effect of the periodic drive is thus the emergence of the discontinuities at frequencies $\omega_{c}=n\omega_{d}/2$ for integer $n$ . The continuous part of the spectrum of ${\cal Q}_{2}^{{\rm b},x}(\omega_{c})$ presents a local maximum ${\cal Q}_{2}^{{\rm b},x}(\omega_{c}=\sqrt{\Gamma_{0}\Omega_{0}})=D_{0}$ for $\sqrt{\Gamma_{0}\Omega_{0}}$ away from integer and half-integer multiples of $\omega_{d}$ (unless $\omega_{d}=2\sqrt{\Gamma_{0}\Omega_{0}}/m$ for some $m\in\mathbb{N}^{*}$ ). In Fig. 2, we show a comparison between our analytical computation for the rescaled variance ${\cal Q}_{2}^{x}(\omega_{c})$ as a function of $\omega_{c}$ . The agreement is excellent for the background ${\cal Q}_{2}^{{\rm b},x}(\omega_{c})$ and relatively good for the spikes.

While the average alternating charge ${\cal Q}_{1}^{x}(\omega_{c})$ is independent of the diffusion coefficient $D(t)$ , the variance ${\cal Q}_{2}^{x}(\omega_{c})$ is independent of the drift $\mu(t)$ . For a finite value of $\Gamma_{0}$ , one has that ${\cal Q}_{2}^{x}(\omega_{c})\to 0$ as $\omega_{c}\to\infty$ which is consistent with the absence of a term of the form $\partial_{x}^{2}G$ in equation Eq. (82). Finally, the overdamped limit can easily be obtained in this example by taking the limit $\Gamma_{0}\to\infty$ . In this limit, the evolution equation for the effective propagator $G(x,t|x_{0},t_{0})$ reads

\partial_{t}G=\Omega(t)\partial_{x}\left((x-\mu(t))G\right)+D(t)\partial_{x}^{2}G\;.

(94)

The background variance reads ${\cal Q}_{2}^{\rm b}(\omega_{c})=D_{0}\omega_{c}^{2}/(\Omega_{0}^{2}+\omega_{c}^{2})$ . In the limit $\omega_{c}\to\infty$ , one has that ${\cal Q}_{2}(\omega_{c})\to D_{0}$ , which is consistant with the prefactor of $\partial_{x}^{2}G$ being $D(t)$ and $\overline{D(t)}^{\omega_{d}}=D_{0}$ .

V.2 Periodically driven run-and-tumble particle

We next consider a simple extension of the model discussed in the previous example to characterise a particle subject to telegraphic noise confined in a time varying potential, i.e.

\dot{x}(t)=-\Omega(t)x(t)+v(t)\sigma(t)\;,\;\;x(0)=0\;,

(95)

where $\sigma(t)=\pm 1$ is a telegraphic noise with $\Big{\langle}\sigma(t)\Big{\rangle}=0$ and $\Big{\langle}\sigma(t_{1})\sigma(t_{2})\Big{\rangle}=e^{-2\gamma|t_{2}-t_{1}|}$ . The noise $\sigma(t)$ is not Gaussian in this framework, which in turn renders the alternating charge $Q_{\omega_{c};x}(t)$ not Gaussian either. Hence

{\cal Q}_{p}^{x}(\omega_{c})\neq 0\;,\;\;p>2\;.

(96)

In Fig. 3, we have plotted the fourth order rescaled cumulant ${\cal Q}_{4}^{{\rm b},x}(\omega_{c})$ for this specific model, showing indeed that it converges to a non-zero value.

At any time $t$ , the process $(x,\sigma)$ is Markovian and this problem is therefore covered by the general theory exposed in the previous sections. The propagator $G_{\sigma,\sigma_{0}}(x,t|x_{0},t_{0})$ satisfies the equation

	$\displaystyle\partial_{t}G_{\sigma,\sigma_{0}}=$	$\displaystyle-\sigma v(t)\partial_{x}G_{\sigma,\sigma_{0}}+\Omega(t)\partial_{x}(xG_{\sigma,\sigma_{0}})$		(97)
		$\displaystyle+\gamma\left(G_{-\sigma,\sigma_{0}}-G_{\sigma,\sigma_{0}}\right)\;.$

One can again obtain the value of the position at time $t$ as

x(t)=\int_{0}^{t}d\tau\,e^{-\int_{\tau}^{t}\Omega(\tau^{\prime})d\tau^{\prime}}v(\tau)\sigma(\tau)\;.

(98)

Let us consider the special case where $\Omega(t)=\Omega_{0}$ . We introduce the Fourier series

v(t)=\sum_{k=-\infty}^{\infty}v_{k}\,e^{ik\omega_{d}t}\;.

(99)

As $\Big{\langle}x(t)\Big{\rangle}=0$ , it yields immediately ${\cal Q}_{1}(\omega_{c})=0$ . The two-time connected correlation function reads instead

	$\displaystyle\Big{\langle}x(t_{2})x(t_{1})\Big{\rangle}_{\rm co}=e^{-\Omega_{0}(t_{1}+t_{2})}$		(100)
	$\displaystyle\times\int_{0}^{t_{1}}d\tau_{1}\int_{0}^{t_{2}}d\tau_{2}\,e^{\Omega_{0}(\tau_{1}+\tau_{2})-2\gamma\|\tau_{1}-\tau_{2}\|}v(\tau_{1})v(\tau_{2})\;.$

Inserting the Fourier expansion of $v(t)$ result and computing explicitly the integrals, we identify this equation for large but finite time with

	$\displaystyle\Big{\langle}x(t_{2})x(t_{1})\Big{\rangle}_{\rm co}\approx$		(101)
	$\displaystyle\sum_{l_{2}=1}^{\infty}e^{-\lambda_{l_{2}}(t_{2}-t_{1})}\sum_{k_{1},k_{2}=-\infty}^{\infty}C_{0,l_{2}}^{k_{1},k_{2}}e^{i\omega_{d}(k_{1}t_{1}+k_{2}t_{2})}\;.$

We thus obtain the Floquet spectrum for this model with $\lambda_{0}=0$ and

\lambda_{1}=\min(\Omega_{0},2\gamma)\leq\lambda_{2}=\max(\Omega_{0},2\gamma)\;,

(102)

as well as the coefficients

$\displaystyle C_{0,2\gamma}^{k_{1},k_{2}}=$	$\displaystyle\frac{v_{k_{1}}v_{k_{2}}}{[\Omega_{0}+2\gamma+ik_{1}\omega_{d}][\Omega_{0}-2\gamma+ik_{2}\omega_{d}]}\;,$	(103)
$\displaystyle C_{0,\Omega_{0}}^{k_{1},k_{2}}=$	$\displaystyle\sum_{n=-\infty}^{\infty}\frac{v_{n}v_{k_{1}-n}\delta_{k_{2},0}}{\Omega_{0}+2\gamma+i(k_{1}-n)\omega_{d}}$	(104)
	$\displaystyle\times\left(\frac{2}{2\Omega_{0}+ik_{1}\omega_{d}}-\frac{1}{\Omega_{0}-2\gamma+in\omega_{d}}\right)\;.$

Note that for this model, we expect a dynamical phase transition for $\Omega_{0}=2\gamma$ where the two Floquet eigenvalues cross. In this example, the background rescaled variance of the charge ${\cal Q}_{2}^{\rm b}(\omega_{c})$ for $\omega_{c}$ different from integer and half-integer multiples of $\omega_{d}$ depends on the full time oscillations of $v(t)$ and reads

	$\displaystyle{\cal Q}_{2}^{{\rm b},x}(\omega_{c})=\frac{2\gamma\|v_{0}\|^{2}\omega_{c}^{2}}{(\Omega_{0}^{2}+\omega_{c}^{2})(4\gamma^{2}+\omega_{c}^{2})}$		(105)
	$\displaystyle+\frac{\omega_{c}^{2}}{\Omega_{0}^{2}+\omega_{c}^{2}}\sum_{n=1}^{\infty}\frac{4\gamma\|v_{n}\|^{2}(4\gamma^{2}+n^{2}\omega_{d}^{2}+\omega_{c}^{2})}{(4\gamma^{2}+n^{2}\omega_{d}^{2})^{2}+2(4\gamma^{2}-n^{2}\omega_{d}^{2})\omega_{c}^{2}+\omega_{c}^{4}}\;.$

Taking the limit $\omega_{c}\to\infty$ for finite $\gamma$ , one obtains that ${\cal Q}_{2}^{{\rm b},x}(\omega_{c})\to 0$ , which is consistent with the absence of term of the form $\partial_{x}^{2}G$ in Eq. (97). On the other hand, taking the limit $\gamma,v_{n}\to\infty$ with $|v_{n}|^{2}/\gamma=O(1)$ , one obtains the result for the overdamped Brownian motion ${\cal Q}_{2}^{\rm b}(\omega_{c})=D_{0}\omega_{c}^{2}/(\Omega_{0}^{2}+\omega_{c}^{2})$ where $D_{0}=(|v_{0}|^{2}+2\sum_{n>1}|v_{n}|^{2})/(2\gamma)$ .

V.3 Two-level system

In the previous two examples, we were able to compute the first few cumulants of the alternating charge. Achieving to obtain numerically the fourth order cumulant with enough precision was a hard task that required $\sim 10^{9}$ realisations. However, computing the cumulant generating function itself beyond Gaussian models, even numerically, seems a formidable task. We therefore consider next a different example, where the cumulant generating function can actually be evaluated numerically. Specifically, we consider a two level system, with a probability vector $|P(t)\rangle=(p_{1}(t),p_{2}(t))^{\rm T}$ , where $p_{i}(t)$ is the probability for the system to be in state $i$ at time $t$ . $|P(t)\rangle$ evolves in time according to a master equation,

|\dot{P}(t)\rangle=R(t)|P(t)\rangle\;,\;\;R(t)=\begin{pmatrix}-k_{1}(t)&k_{2}(t)\\ k_{1}(t)&-k_{2}(t)\end{pmatrix}\;,

(106)

where $R(t)$ is periodic with angular frequency $\omega_{d}$ . A given realisation of the process is a sequence of successive states of the system, $i_{1}\to i_{2}\to...\to i_{n}$ and the corresponding sequence of times $t_{k}$ at which the system jumped from state $i_{k}$ to $i_{k+1}$ . For any realisation, we define its empirical alternating charge at frequency $\omega_{c}$ as

Q_{\omega_{c}}(t)=\sum_{k}\Big{(}i_{k+1}-i_{k}\Big{)}\cos(\omega_{c}t_{k})\;,

(107)

where the sum is over the times $t_{k}<t$ where the system jumps from one state to the other of the specific realization. As in the continuous state systems, we are interested in the fluctuations of this quantity as $t\to\infty$ . As the system is only composed of two states, there is no cycle and therefore the direct current (DC) corresponding to $\omega_{c}=0$ is identically zero. Recent progress has been made on this particular system Takahashi et al. [2020] and in particular an exact representation of the CGF has been obtained for the DC currents. Following the method presented in Takahashi et al. [2020] but for alternating charge (see details in Appendix D), we obtain

\chi_{\mu}(\omega_{c})=-\lim_{t\to\infty}\int_{0}^{t}\frac{d\tau}{t}\left[k_{1}(\tau)+k_{2}(\tau)e^{-\mu\cos(\omega_{c}\tau)}y_{\mu}(\tau)\right]\;.

(108)

In that expression, the function $y_{\mu}(t)$ satisfies a first order non-linear differential equation for fixed $\mu$ , which reads Takahashi et al. [2020]

\dot{y}_{\mu}=(y_{\mu}-e^{\mu\cos(\omega_{c}t)})(k_{1}+k_{2}e^{-\mu\cos(\omega_{c}t)}y_{\mu})\;.

(109)

Using this particular representation, we show in Fig. (5) the CGF $\chi_{\mu}(\omega_{c})$ as a function of $\omega_{c}$ , calculated for $\mu=3$ and the following rates:

	$\displaystyle k_{1}(t)$	$\displaystyle=\exp\Big(\cos(\omega_{d}t)-1\Big{missing}),$
	$\displaystyle k_{2}(t)$	$\displaystyle=\exp\Big(-4\sin(5\omega_{d}t)+3\sin(7\omega_{d}t)-2\Big{missing})\;.$		(110)

The blue line was calculated for the series of frequencies $\omega_{c}=\frac{n}{10\pi}$ for integer values of $n$ , incommensurate with $\omega_{d}=1$ . The red dots were calculated for the series of frequencies $\omega_{c}=\frac{n}{m}$ for $m\in\mathbb{N}^{*}\leq 50$ and all $n\in\mathbb{N}\leq 6m$ , commensurate with $\omega_{d}$ . Both the commensurate and incommensurate cases were estimated by first solving numerically the differential equation (133) and then evaluating the integral in equation (108) up to $t=2\times 10^{4}$ .

VI Conclusion

In this article we have computed explicitly the large time fluctuations of the alternating charge $Q_{\omega_{c};\alpha}(t)$ (defined as the time integrated instantaneous alternating current in direction $\alpha$ ) for a periodically driven stochastic system. We have shown that the cumulant generating function of this quantity (and as a consequence its large deviation function) is not a continuous function of the frequency $\omega_{c}$ of the charge. In particular, we have shown that there exists both a continuous background for any value of the charge and driving frequency, respectively $\omega_{c}$ and $\omega_{d}$ , and additional peaks for commensurate pair of frequencies. This general result has been confirmed by considering some exactly solvable models.

It would be interesting to test these fairly general results experimentally, especially on systems where current fluctuations play an important role, e.g. trapped ions setups. For this, it is important to consider the finite time version of the results.

We have focused our interest in this article on systems with a discrete gaped Floquet spectrum on a simply-connected domain, where the position of the particle are confined to a finite portion of space. It would be interesting to consider similar results for systems on a finite loop, where there might be long-lasting average currents. Some simple exactly solvable cases can be considered, hinting that the cumulants of the alternating charge behave quite differently in this case, but obtaining a general result remains a challenge at the moment.

Acknowledgements : O.R. is the incumbent of the Shlomo and Michla Tomarin career development chair, and is supported by the Abramson Family Center for Young Scientists, the Minerva and by the Israel Science Foundation, Grant No. 950/19.

Appendix A Simplification of the cumulants

Let us first show that the CGF $\chi_{\mu}(\omega_{c})$ as defined in Eq. (28) is identical to the CGF

\tilde{\chi}_{\mu}(\omega_{c})=\lim_{t\to\infty}\frac{1}{t}\ln\Big{\langle}e^{\mu q_{\omega_{c};\alpha}(t)}\Big{\rangle}\;,

(111)

where we remind that

q_{\omega_{c};\alpha}(t)=\omega_{c}\int_{0}^{t}d\tau\,x_{\alpha}(\tau)\sin(\omega_{c}\tau)\;.

(112)

To show this identity, we use that

Q_{\omega_{c};\alpha}(t)=q_{\omega_{c};\alpha}(t)+x_{\alpha}(t)\cos(\omega_{c}t)-x_{\alpha}(0),

(113)

where $x_{\alpha}(0)=x_{0}$ is fixed. We introduce the joint PDF $P_{t}(q,{\bf x})$ for the respective random variables $q_{\omega_{c}}(t)$ and the final position ${\bf x}(t)$ . We use Bayes’ theorem as

P_{t}(q,{\bf x})=P_{t}(q|{\bf x})P_{t}({\bf x})\;,

(114)

where $P_{t}({\bf x})=G({\bf x},t|{\bf x}_{0},0)$ is the PDF of the final position ${\bf x}(t)$ and $P_{t}(q|{\bf x})$ is the PDF of $q_{\omega_{c};\alpha}(t)$ for a fixed final position ${\bf x}(t)={\bf x}$ . In the large time limit, the final position ${\bf x}(t)$ is of order $O(1)$ for the trapped systems that we consider while one naturally expects that the random variable $q_{\omega_{c};\alpha}(t)=O(t)$ . We thus expect the large deviation form

P_{t}(q|{\bf x})\asymp e^{-t\varphi\left(\frac{q}{t},{\bf x}\right)}\;,

(115)

while there is no large deviation form for the final position

-\lim_{t\to\infty}\frac{1}{t}\ln G({\bf x},t|{\bf x}_{0},0)=\lambda_{0}=0\;.

(116)

We may then show that on the one hand the MGF for $Q_{\omega_{c};\alpha}(t)$ reads in the large time limit

	$\displaystyle\Big{\langle}e^{\mu Q_{\omega_{c};\alpha}(t)}\Big{\rangle}$		(117)
	$\displaystyle=e^{-\mu x_{0}}\int d{\bf x}\,G({\bf x},t\|{\bf x}_{0},0)e^{\mu x_{\alpha}\cos(\omega_{c}t)}\int dq\,e^{\mu q}P_{t}(q\|{\bf x})$
	$\displaystyle\asymp e^{-\mu x_{0}}\int d{\bf x}\,f_{0}({\bf x},t)e^{\mu x_{\alpha}\cos(\omega_{c}t)}\int\frac{dQ}{t}\,e^{t\left[\mu Q-\varphi(Q,{\bf x})\right]}\;,$

while on the other hand the MGF for $q_{\omega_{c};\alpha}(t)$ reads

	$\displaystyle\Big{\langle}e^{\mu\,q_{\omega_{c};\alpha}(t)}\Big{\rangle}$		(118)
	$\displaystyle=e^{-\mu x_{0}}\int d{\bf x}\,G({\bf x},t\|{\bf x}_{0},0)\int dq\,e^{\mu q}P_{t}(q\|{\bf x})$
	$\displaystyle\asymp e^{-\mu x_{0}}\int d{\bf x}\,f_{0}({\bf x},t)\int\frac{dQ}{t}\,e^{t\left[\mu Q-\varphi(Q,{\bf x})\right]}\;.$

Using a saddle point-approximation then naturally yields the result

\chi_{\mu}(\omega_{c})=\tilde{\chi}_{\mu}(\omega_{c})=\max_{Q,{\bf x}}\left[\mu Q-\varphi(Q,{\bf x})\right]\;.

(119)

As the CGF for the two random variables $Q_{\omega_{c};\alpha}(t)$ and $q_{\omega_{c};\alpha}(t)$ are identical, it yields immediately the identity for the rescaled cumulants

{\cal Q}_{p}^{\alpha}(\omega_{c})=\lim_{t\to\infty}\frac{\Big{\langle}Q_{\omega_{c};\alpha}^{p}(t)\Big{\rangle}_{\rm co}}{t}=\lim_{t\to\infty}\frac{\Big{\langle}q_{\omega_{c};\alpha}^{p}(t)\Big{\rangle}_{\rm co}}{t}

(120)

of arbitrary order $p$ .

Appendix B Expression of the third order cumulant ${\cal Q}_{3}^{\alpha}(\omega_{c})$

Let us now consider the cumulant of order $p=3$ . To obtain its value, we first consider the three times connected correlation function for $t_{3}>t_{2}>t_{1}\gg 1$

\displaystyle\Big{\langle}x_{\alpha}(t_{1})x_{\alpha}(t_{2})x_{\alpha}(t_{3})\Big{\rangle}_{\rm co}=\sum_{l_{2},l_{3}=1}^{\infty}\prod_{j=2}^{3}e^{-\lambda_{l_{j}}(t_{j}-t_{j-1})}C_{0,l_{2},l_{3}}(t_{1},t_{2},t_{3})-\sum_{l_{3}=1}^{\infty}\prod_{j=2}^{3}e^{-\lambda_{l_{3}}(t_{j}-t_{j-1})}C_{0,l_{3}}(t_{1},t_{3})C_{0}(t_{2})\;.

(121)

Note that $p=3$ is the first term for which there are more than one partition of $\{1,\cdots,p\}$ that give a non-zero contribution to this connected correlation, namely $\{1,2,3\}$ and $\{1,3\}\{2\}$ .

	$\displaystyle{\cal Q}_{3}^{\alpha}(\omega_{c})=$	$\displaystyle\frac{3\omega_{c}^{3}}{4i}\sum_{k_{1},k_{2},k_{3}=-\infty}^{\infty}\sum_{\sigma_{1},\sigma_{2},\sigma_{3}=\pm 1}\sum_{l_{3}=1}^{\infty}\frac{\sigma_{1}\sigma_{2}\sigma_{3}}{\lambda_{l_{3}}+i(\sigma_{3}\omega_{c}-k_{3}\omega_{d})}\delta_{\sum_{j=1}^{3}k_{j}\omega_{d},\sum_{j=1}^{3}\sigma_{j}\omega_{c}}$		(122)
		$\displaystyle\times\left[\sum_{l_{2}=1}^{\infty}\frac{C_{0,l_{2},l_{3}}^{k_{1},k_{2},k_{3}}}{\lambda_{l_{2}}+i\left[(\sigma_{2}+\sigma_{3})\omega_{c}-(k_{2}+k_{3})\omega_{d}\right]}-\frac{C_{0,l_{3}}^{k_{1},k_{3}}C_{0}^{k_{2}}}{\lambda_{l_{3}}+i\left[(\sigma_{2}+\sigma_{3})\omega_{c}-(k_{2}+k_{3})\omega_{d}\right]}\right]\;.$

This third order cumulant is always zero if the frequencies $\omega_{c}$ and $\omega_{d}$ are incommensurate. It is only non-zero for frequencies of the form $\omega_{c}=n\omega_{d}$ and $\omega_{c}=n\omega_{d}/3$ for $n\in\mathbb{N}^{*}$ .

Appendix C Expression of the fourth order cumulant ${\cal Q}_{4}^{\alpha}(\omega_{c})$

One can use the expression of the connected four point functions for $t_{4}>t_{3}>t_{2}>t_{1}\gg 1$ with $\delta t_{j}=t_{j}-t_{j-1}$ for $j=2,3,4$ , which reads

	$\displaystyle\Big{\langle}x_{\alpha}(t_{1})x_{\alpha}(t_{2})x_{\alpha}(t_{3})x_{\alpha}(t_{4})\Big{\rangle}_{\rm co}=$		(123)
	$\displaystyle\sum_{l_{2},l_{3},l_{4}=1}^{\infty}\prod_{j=2}^{4}e^{-\lambda_{l_{j}}\delta t_{j}}C_{0,l_{2},l_{3},l_{4}}(t_{1},t_{2},t_{3},t_{4})-\sum_{l_{2},l_{4}=1}^{\infty}e^{-\lambda_{l_{4}}(\delta t_{4}+\delta t_{3})-\lambda_{l_{2}}\delta t_{2}}C_{0,l_{2},l_{4}}(t_{1},t_{2},t_{4})C_{0}(t_{3})$
	$\displaystyle-\sum_{l_{3},l_{4}=1}^{\infty}e^{-\lambda_{l_{4}}\delta t_{4}-\lambda_{l_{3}}\delta t_{3}-\lambda_{l_{3}}\delta t_{2}}C_{0,l_{3},l_{4}}(t_{1},t_{3},t_{4})C_{0}(t_{2})$
	$\displaystyle-\sum_{l_{3},l_{4}=1}^{\infty}e^{-\lambda_{l_{4}}\delta t_{4}-(\lambda_{l_{4}}+\lambda_{l_{3}})\delta t_{3}-\lambda_{l_{3}}\delta t_{2}}\left[C_{0,l_{4}}(t_{1},t_{4})C_{0,l_{3}}(t_{2},t_{3})+C_{0,l_{4}}(t_{2},t_{4})C_{0,l_{3}}(t_{1},t_{3})\right]$
	$\displaystyle+\sum_{l_{4}=1}^{\infty}\prod_{j=2}^{4}e^{-\lambda_{l_{4}}\delta t_{j}}C_{0,l_{4}}(t_{1},t_{4})C_{0}(t_{2})C_{0}(t_{3})\;.$		(124)

The partition of $\{1,2,3,4\}$ that contribute besides the identity are for two blocks $\{1,2,4\}\{3\}$ , $\{1,3,4\}\{2\}$ , $\{1,3\}\{2,4\}$ , $\{2,4\}\{1,3\}$ and for three blocks $\{1,4\}\{2\}\{3\}$ . Using the explicit expressions for the rates $\mu_{j}$ ’s corresponding to each partition, we obtain

	$\displaystyle{\cal Q}_{4}^{\alpha}(\omega_{c})=\frac{3\omega_{c}^{4}}{2}\sum_{k_{1},\cdots,k_{4}=-\infty}^{\infty}\sum_{\sigma_{1},\cdots,\sigma_{4}=\pm 1}\sum_{l_{4}=1}^{\infty}\frac{\prod_{j=1}^{4}\sigma_{j}}{\lambda_{l_{4}}+i(\sigma_{4}\omega_{c}-k_{4}\omega_{d})}\delta_{\sum_{j=1}^{4}k_{j}\omega_{d},\sum_{j=1}^{4}\sigma_{j}\omega_{c}}$		(125)
	$\displaystyle\times\left(\sum_{l_{2},l_{3}=1}^{\infty}\frac{C_{0,l_{2},l_{3},l_{4}}^{k_{1},k_{2},k_{3},k_{4}}}{\prod_{n=2}^{3}\left[\lambda_{l_{n}}+i\sum_{j=n}^{4}(\sigma_{j}\omega_{c}-k_{j}\omega_{d})\right]}+\frac{C_{0,l_{4}}^{k_{1},k_{4}}C_{0}^{k_{2}}C_{0}^{k_{3}}}{\prod_{n=2}^{3}\left[\lambda_{l_{4}}+i\sum_{j=n}^{4}(\sigma_{j}\omega_{c}-k_{j}\omega_{d})\right]}\right.$
	$\displaystyle-\sum_{l_{3}=1}^{\infty}\frac{C_{0,l_{3},l_{4}}^{k_{1},k_{3},k_{4}}C_{0}^{k_{2}}}{\prod_{n=2}^{3}\left[\lambda_{l_{3}}+i\sum_{j=n}^{4}(\sigma_{j}\omega_{c}-k_{j}\omega_{d})\right]}-\sum_{l_{2}=1}^{\infty}\frac{C_{0,l_{2},l_{4}}^{k_{1},k_{2},k_{4}}C_{0}^{k_{3}}}{\left[\lambda_{l_{4}}+i\sum_{j=3}^{4}(\sigma_{j}\omega_{c}-k_{j}\omega_{d})\right]\left[\lambda_{l_{2}}+i\sum_{j=2}^{4}(\sigma_{j}\omega_{c}-k_{j}\omega_{d})\right]}$
	$\displaystyle\left.-\sum_{l_{3}=1}^{\infty}\frac{1}{\sum_{j=3}^{4}\left[\lambda_{l_{j}}+i(\sigma_{j}\omega_{c}-k_{j}\omega_{d})\right]}\left[\frac{C_{0,l_{4}}^{k_{1},k_{4}}C_{0,l_{3}}^{k_{2},k_{3}}}{\lambda_{l_{4}}+i\sum_{j=2}^{4}(\sigma_{j}\omega_{c}-k_{j}\omega_{d})}+\frac{C_{0,l_{4}}^{k_{2},k_{4}}C_{0,l_{3}}^{k_{1},k_{3}}}{\lambda_{l_{3}}+i\sum_{j=2}^{4}(\sigma_{j}\omega_{c}-k_{j}\omega_{d})}\right]\right)\;.$

Appendix D Cumulant generating function for the alternating charge of the 2 level system

A convenient way to compute the CGF $\chi_{\mu}(\omega_{c})$ is to introduce a counting vector $|P_{s}(t)\rangle$ (see Touchette [2009] for details about this method), such that the $i^{\rm th}$ component of $|P_{s}(t)\rangle$ is the ensemble average of $e^{\mu Q_{\omega_{c}}(t)}$ on trajectories that are in the state $i$ at time $t$ . The initial condition is taken such that $|P_{\mu}(0)\rangle=|P(0)\rangle$ at $t=0$ , and it evolves with time according to

	$\displaystyle\|\dot{P}_{\mu}(t)\rangle$	$\displaystyle=R_{\mu}(t)\|P_{\mu}(t)\rangle\;,$		(126)
	$\displaystyle R_{\mu}(t)$	$\displaystyle=\begin{pmatrix}-k_{1}(t)&k_{2}(t)e^{-\mu\cos(\omega_{c}t)}\\ k_{1}(t)e^{\mu\cos(\omega_{c}t)}&-k_{2}(t)\end{pmatrix}\;,$		(127)

and $R_{\mu}(t)$ is called the “tilted rate matrix”. The vector $|P_{\mu}(t)\rangle$ can be (formally) expressed at any time $t$ as

|P_{\mu}(t)\rangle=\mathcal{T}e^{\int_{0}^{t}R_{\mu}(\tau)d\tau}|P_{\mu}(0)\rangle\;,

(128)

where $\mathcal{T}$ is the time ordered product. Using this counting vector method, the cumulant generating function is given by Lebowitz and Spohn [1999]

	$\displaystyle\chi_{\mu}(\omega_{c})$	$\displaystyle=$	$\displaystyle\ln\langle e^{\mu Q_{\omega_{c}}(t)}\rangle=\lim_{t\to\infty}\frac{1}{t}\ln\langle 1\|P_{\mu}(t)\rangle$		(129)
		$\displaystyle=$	$\displaystyle\lim_{t\to\infty}\frac{1}{t}\ln\langle 1\|\mathcal{T}e^{\int_{0}^{t}R_{\mu}(\tau)d\tau}\|P(0)\rangle\;,$		(129)

where $\langle 1|=(1,1)$ . While this analytical formula is exact, it is not very convenient as the time ordered exponential in Eq. (128) is a complicated object and hard to compute in practice. To simplify the computation of the counting vector, we first introduce a dynamical invariant of the process. We then obtain an explicit expression for the counting vector $|P_{\mu}(t)\rangle$ in terms of this dynamical invariant and finally compute the CGF.

The derivation below follows the derivation provided in Takahashi et al. [2020], but generalizes it for all the current’s Fourier components.

D.1 Dynamical invariant

A dynamical invariant corresponds to a matrix $F_{\mu}(t)$ which is diagonalisable, has time-independent eigenvalues that we set here to $\pm 1$ for convenience and satisfies for any time the differential equation

\dot{F}_{\mu}(t)=R_{\mu}(t)F_{\mu}(t)-F_{\mu}(t)R_{\mu}(t)\;.

(130)

A convenient way to parametrise a matrix satisfying the first two properties is Takahashi et al. [2020]

F_{\mu}(t)=\begin{pmatrix}z_{\mu}(t)&(1+z_{\mu}(t))y_{\mu}(t)^{-1}\\ (1-z_{\mu}(t))y_{\mu}(t)&-z_{\mu}(t)\end{pmatrix}\;,

(131)

with left and right time-dependent orthonormal eigenvectors $\langle r_{\sigma}(t)|l_{\sigma}^{\prime}(t)\rangle=\delta_{\sigma,\sigma^{\prime}}$ with $\sigma,\sigma^{\prime}=\pm$ that read Takahashi et al. [2020]

	$\displaystyle\|l_{+}(t)\rangle=\begin{pmatrix}\displaystyle 1\\ \displaystyle\frac{1-z_{\mu}(t)}{1+z_{\mu}(t)}y_{\mu}(t)\end{pmatrix}\;,$		(132)
	$\displaystyle\langle r_{+}(t)\|=\frac{1}{2}\begin{pmatrix}\displaystyle 1+z_{\mu}(t)&\displaystyle\frac{1+z_{\mu}(t)}{y_{\mu}(t)}\end{pmatrix}\;,$
	$\displaystyle\|l_{-}(t)\rangle=\frac{1}{2}\begin{pmatrix}\displaystyle 1-z_{\mu}(t)\\ \displaystyle-(1-z_{\mu}(t))y_{\mu}(t)\end{pmatrix}\;,$
	$\displaystyle\langle r_{-}(t)\|=\begin{pmatrix}\displaystyle 1&\displaystyle-\frac{1+z_{\mu}(t)}{1-z_{\mu}(t)}\frac{1}{y_{\mu}(t)}\end{pmatrix}\;.$

For the matrix $F_{\mu}(t)$ to satisfy Eq. (130), the functions $y_{\mu}(t)$ and $z_{\mu}(t)$ need to satisfy the differential equations

	$\displaystyle\dot{y}_{\mu}=$	$\displaystyle(y_{\mu}-e^{\mu\cos(\omega_{c}t)})(k_{1}+k_{2}e^{-\mu\cos(\omega_{c}t)}y_{\mu})\;,$		(133)
	$\displaystyle\dot{z}_{\mu}=$	$\displaystyle k_{2}e^{-\mu\cos(\omega_{c}t)}(1-z_{\mu})y_{\mu}-k_{1}e^{\mu\cos(\omega_{c}t)}\frac{(1+z_{\mu})}{y_{\mu}}\;.$

Any pair of initial values $(y_{\mu}(0),z_{\mu}(0))$ gives a valid dynamical invariant. Thus, choosing such a pair defines a valid dynamical invariant. We focus here on the long-time behaviour such that the initial condition is not relevant and choose for convenience $y_{\mu}(0)=0$ and $z_{\mu}(0)=-1$ . The equation for $y_{\mu}(t)$ has two fixed point: a stable fixed point $-k_{1}(t)/k_{2}(t)e^{\mu\cos(\omega_{c}t)}<0$ and an unstable fixed point $e^{\mu\cos(\omega_{c}t)}>0$ . Note that starting from $y_{\mu}(0)\leq 0$ , the solution $y_{\mu}(t)$ oscillates around the stable fixed point and remains negative at all time $t$ . The equation for $z_{\mu}(t)$ is linear and can be solved exactly. In particular, it is easy to realise that the function $z_{\mu}(t)$ grows exponentially with time $t$ as $y_{\mu}(\tau)<0$ for any $\tau\in[0,t]$

z_{\mu}(t)\asymp e^{-\int_{0}^{t}d\tau\left[y_{\mu}(\tau)k_{2}(\tau)e^{-\mu\cos(\omega_{c}\tau)}+\frac{k_{1}(\tau)}{y_{\mu}(\tau)}e^{\mu\cos(\omega_{c}\tau)}\right]}\;.

(134)

In the large time limit, the expressions of the eigenvectors in Eq. (135) simplify

	$\displaystyle\|l_{+}(t)\rangle\approx\begin{pmatrix}\displaystyle 1\\ \displaystyle y_{\mu}(t)\end{pmatrix}\;,\langle r_{+}(t)\|\approx\frac{z_{\mu}(t)}{2}\begin{pmatrix}\displaystyle 1&\displaystyle y_{\mu}^{-1}(t)\end{pmatrix}\;,$
	$\displaystyle\|l_{-}(t)\rangle\approx\frac{z_{\mu}(t)}{2}\begin{pmatrix}\displaystyle-1\\ \displaystyle y_{\mu}(t)\end{pmatrix}\;,\langle r_{-}(t)\|\approx\begin{pmatrix}\displaystyle 1&\displaystyle-y_{\mu}^{-1}(t)\end{pmatrix}\,.$		(135)

D.2 Expressing the counting vector and the CGF

Coming back to the problem of finding an expression for the counting vector $|P_{\mu}(t)\rangle$ , one can check by using Eq. (126) together with Eq. (130) that the following identity holds

\partial_{t}(F_{\mu}(t)|P_{\mu}(t)\rangle)=R_{\mu}(t)F_{\mu}(t)|P_{\mu}(t)\rangle\;.

(136)

It is then possible to express the counting vector $|P_{\mu}(t)\rangle$ in the basis of the eigenvectors $|l_{\pm}(t)\rangle$ of $F_{\mu}(t)$ . In this basis, one has that

|P_{\mu}(t)\rangle=\sum_{\sigma=\pm}a_{\sigma}(t)|l_{\sigma}(t)\rangle\,,\,\,a_{\sigma}(t)=\langle r_{\sigma}(t)|P_{\mu}(t)\rangle\,.

(137)

Inserting this equation into Eqs. (126) and (136), one obtains that for $\sigma,\sigma^{\prime}=\pm$

(\langle r_{\sigma^{\prime}}(t)|R_{\mu}(t)|l_{\sigma}(t)\rangle-\langle r_{\sigma^{\prime}}(t)|\dot{l}_{\sigma}(t)\rangle)a_{\sigma}(t)=\delta_{\sigma,\sigma^{\prime}}\dot{a}_{\sigma}(t)\,.

(138)

Inserting the expressions of the eigenvectors in Eq. (132), using the expression of the tilted rate matrix and the differential equations (133), one can then check that $\langle r_{\sigma^{\prime}}(t)|R_{\mu}(t)|l_{\sigma}(t)\rangle=\langle r_{\sigma^{\prime}}(t)|\dot{l}_{\sigma}(t)\rangle$ for $\sigma\neq\sigma^{\prime}$ . On the other hand, Eq. (138) yields

	$\displaystyle\|P_{\mu}(t)\rangle=$	$\displaystyle\sum_{\sigma=\pm}a_{\sigma}(0)e^{\int_{0}^{t}d\tau\phi_{\sigma}(\tau)}\|l_{\sigma}(t)\rangle\;,$		(139)
	$\displaystyle\phi_{\pm}(\tau)=$	$\displaystyle\langle r_{\pm}(\tau)\|R_{\mu}(\tau)\|l_{\pm}(\tau)\rangle-\langle r_{\pm}(\tau)\|\dot{l}_{\pm}(\tau)\rangle\;.$		(140)

We can now express $|P_{\mu}(t)\rangle$ in the large time limit as a function of the two functions $z_{\mu}(t)$ and $y_{\mu}(t)$ . Inserting the expressions of the eigenvectors in Eq. (135), using the expression of the tilted matrix $R_{\mu}(t)$ and the differential equations in Eq (133), one obtains at large time

|P_{\mu}(t)\rangle\asymp e^{-\int_{0}^{t}d\tau\left[k_{1}(\tau)+k_{2}(\tau)e^{-\mu\cos(\omega_{c}\tau)}y_{\mu}(\tau)\right]}\;.

(141)

Finally, using Eq. (129), we obtain the explicit expression of the CGF as given in Eq. (108) of the main text.

References

Zia and Schmittmann [2007] R. Zia and B. Schmittmann, Journal of Statistical Mechanics: Theory and Experiment 2007, P07012 (2007).
Ruelle [2003] D. P. Ruelle, Proceedings of the National Academy of Sciences 100, 3054 (2003).
Zia and Schmittmann [2006] R. K. Zia and B. Schmittmann, Journal of physics A: Mathematical and general 39, L407 (2006).
Brandner et al. [2015] K. Brandner, K. Saito, and U. Seifert, Physical review X 5, 031019 (2015).
Raz et al. [2016] O. Raz, Y. Subasi, and C. Jarzynski, Physical Review X 6, 021022 (2016).
Busiello et al. [2018] D. M. Busiello, C. Jarzynski, and O. Raz, New Journal of Physics 20, 093015 (2018).
Kubo [1966] R. Kubo, Reports on progress in physics 29, 255 (1966).
Gallavotti and Cohen [1995] G. Gallavotti and E. G. D. Cohen, Physical review letters 74, 2694 (1995).
Kurchan [1998] J. Kurchan, Journal of Physics A: Mathematical and General 31, 3719 (1998).
Andrieux and Gaspard [2007a] D. Andrieux and P. Gaspard, Journal of statistical physics 127, 107 (2007a).
Lebowitz and Spohn [1999] J. L. Lebowitz and H. Spohn, Journal of Statistical Physics 95, 333 (1999).
Andrieux and Gaspard [2007b] D. Andrieux and P. Gaspard, Journal of statistical physics 127, 107 (2007b).
Shargel and Chou [2009] B. H. Shargel and T. Chou, Journal of Statistical Physics 137, 165 (2009).
Harris and Schütz [2007] R. J. Harris and G. M. Schütz, Journal of Statistical Mechanics: Theory and Experiment 2007, P07020 (2007).
Barato and Seifert [2015] A. C. Barato and U. Seifert, Physical review letters 114, 158101 (2015).
Gingrich et al. [2016] T. R. Gingrich, J. M. Horowitz, N. Perunov, and J. L. England, Physical review letters 116, 120601 (2016).
Horowitz and Gingrich [2019] J. M. Horowitz and T. R. Gingrich, Nature Physics , 1 (2019).
Proesmans and Horowitz [2019] K. Proesmans and J. M. Horowitz, Journal of Statistical Mechanics: Theory and Experiment 2019, 054005 (2019).
Barato et al. [2018] A. C. Barato, R. Chetrite, A. Faggionato, and D. Gabrielli, New Journal of Physics 20, 103023 (2018).
Koyuk and Seifert [2019] T. Koyuk and U. Seifert, Physical review letters 122, 230601 (2019).
Koyuk and Seifert [2020] T. Koyuk and U. Seifert, Phys. Rev. Lett. 125, 260604 (2020).
Koyuk et al. [2018] T. Koyuk, U. Seifert, and P. Pietzonka, Journal of Physics A: Mathematical and Theoretical 52, 02LT02 (2018).
Nyawo and Touchette [2016] P. T. Nyawo and H. Touchette, Physical Review E 94, 032101 (2016).
Nemoto and Sasa [2011a] T. Nemoto and S.-i. Sasa, Physical Review E 83, 030105 (2011a).
Nemoto and Sasa [2011b] T. Nemoto and S.-i. Sasa, Physical Review E 84, 061113 (2011b).
Proesmans and Derrida [2019] K. Proesmans and B. Derrida, Journal of Statistical Mechanics: Theory and Experiment 2019, 023201 (2019).
Derrida and Sadhu [2019] B. Derrida and T. Sadhu, Journal of Statistical Physics 176, 773 (2019).
Nyawo and Touchette [2017] P. T. Nyawo and H. Touchette, EPL (Europhysics Letters) 116, 50009 (2017).
Fischer et al. [2018] L. P. Fischer, P. Pietzonka, and U. Seifert, Physical Review E 97, 022143 (2018).
Barato and Chetrite [2018] A. C. Barato and R. Chetrite, Journal of Statistical Mechanics: Theory and Experiment 2018, 053207 (2018).
Bertini et al. [2018] L. Bertini, R. Chetrite, A. Faggionato, and D. Gabrielli, in Annales Henri Poincaré, Vol. 19 (Springer, 2018) pp. 3197–3238.
Chabane et al. [2020] L. Chabane, R. Chetrite, and G. Verley, Journal of Statistical Mechanics: Theory and Experiment 2020, 033208 (2020).
Blatt and Wineland [2008] R. Blatt and D. Wineland, Nature 453, 1008 (2008).
Brownnutt et al. [2015] M. Brownnutt, M. Kumph, P. Rabl, and R. Blatt, Reviews of modern Physics 87, 1419 (2015).
Beanland et al. [2009] K. Beanland, J. W. Roberts, and C. Stevenson, The American Mathematical Monthly 116, 531 (2009).
Gammaitoni et al. [1998] L. Gammaitoni, P. Hänggi, P. Jung, and F. Marchesoni, Reviews of modern physics 70, 223 (1998).
Note [1] This assumption implies some consequences on our results. Although a generalization to non-simply connected domains is possible, we do not consider these cases here.
Caceres and Lobos [2006] M. O. Caceres and A. M. Lobos, Journal of Physics A: Mathematical and General 39, 1547 (2006).
Kim et al. [2010] C. Kim, P. Talkner, E. K. Lee, and P. Hänggi, Chemical Physics 370, 277 (2010).
Note [2] Note that a similar quantity can be defined with a $\sin$ function instead of the $\cos$ .
Chetrite and Touchette [2013] R. Chetrite and H. Touchette, Physical review letters 111, 120601 (2013).
Oseledets [1968] V. I. Oseledets, Trudy Moskovskogo Matematicheskogo Obshchestva 19, 179 (1968).
Speed [1983] T. Speed, Australian Journal of Statistics 25, 378 (1983).
Shneidman et al. [1994] V. A. Shneidman, P. Jung, and P. Hänggi, EPL (Europhysics Letters) 26, 571 (1994).
Nikitin et al. [2007] A. Nikitin, N. Stocks, and A. Bulsara, Physical Review E 76, 041138 (2007).
Nikitin et al. [2005] A. Nikitin, N. Stocks, and A. Bulsara, Physics Letters A 334, 12 (2005).
Takahashi et al. [2020] K. Takahashi, Y. Hino, K. Fujii, and H. Hayakawa, Journal of Statistical Physics 181, 2206 (2020).
Touchette [2009] H. Touchette, Physics Reports 478, 1 (2009).

	$\displaystyle\|l_{+}(t)\rangle=\begin{pmatrix}\displaystyle 1\\ \displaystyle\frac{1-z_{\mu}(t)}{1+z_{\mu}(t)}y_{\mu}(t)\end{pmatrix}\;,$		(132)
	$\displaystyle\langle r_{+}(t)\|=\frac{1}{2}\begin{pmatrix}\displaystyle 1+z_{\mu}(t)&\displaystyle\frac{1+z_{\mu}(t)}{y_{\mu}(t)}\end{pmatrix}\;,$
	$\displaystyle\|l_{-}(t)\rangle=\frac{1}{2}\begin{pmatrix}\displaystyle 1-z_{\mu}(t)\\ \displaystyle-(1-z_{\mu}(t))y_{\mu}(t)\end{pmatrix}\;,$
	$\displaystyle\langle r_{-}(t)\|=\begin{pmatrix}\displaystyle 1&\displaystyle-\frac{1+z_{\mu}(t)}{1-z_{\mu}(t)}\frac{1}{y_{\mu}(t)}\end{pmatrix}\;.$