Spectral theory of fluctuations in time-average statistical mechanics of reversible and driven systems

We present a spectral-theoretic approach to time-average statistical mechanics for general, non-equilibrium initial conditions. We consider the statistics of bounded, local additive functionals of reversible as well as irreversible ergodic stochastic dynamics with continuous or discrete state-space. We derive exact results for the mean, fluctuations and correlations of time average observables from the eigenspectrum of the underlying generator of Fokker-Planck or master equation dynamics, and discuss the results from a physical perspective. Feynman-Kac formulas are re-derived using It\^o calculus and combined with non-Hermitian perturbation theory. The emergence of the universal central limit law in a spectral representation is shown explicitly on large deviation time-scales. For reversible dynamics with equilibrated initial conditions we derive a general upper bound to fluctuations in terms of an integral of the return probability. Simple, exactly solvable examples are analyzed to demonstrate how to apply the theory. As a biophysical example we revisit the Berg-Purcell problem on the precision of concentration measurements by a single receptor. Our results are directly applicable to a diverse range of phenomena underpinned by time-average observables and additive functionals in physical, chemical, biological, and economical systems.


I. INTRODUCTION
Many experiments on soft and biological matter, such as single-particle tracking [1][2][3][4] and single-molecule spectroscopy [5][6][7][8][9][10][11][12][13], probe individual trajectories. It is typically not feasible to repeat these experiments sufficiently many times in order to allow for ensemble-averaging. It is, however, straightforward to analyze the data by means of time-averaging along individual realizations, which naturally leads to the study of time-averages, that is, functionals of stochastic processes.
From a theoretical point of view analytical results were obtained for the occupation time statistics for discretestate Markov switching [47][48][49][50], for the local time at zero and occupation time above zero of a Brownian particle diffusing in a simple one-dimensional potential [46,52,53], the occupation time inside a spherical domain of a Brownian particle moving in free space [51] and for a free, uniformly biased and harmonically bound particle undergoing subdiffusion [54,55]. Exact results were also obtained for occupation time statistics for a general class of Markov processes [56] and a discrete stationary non-Markovian sequence [57]. Large deviation functions * agodec@mpibpc.mpg.de for various non-linear functionals of a class of Gaussian stationary Markov processes were studied in [45]. Numerous important results on functionals have also been obtained in the context of persistence in spatially extended non-equilibrium systems [58]. Exact results were recently obtained on local times for projected observables in stochastic many-body systems [59,60], which provided insight into the emergence of memory on the level of individual non-Markovian trajectories. Notwithstanding, a general approach to fluctuations of time-average statistical mechanics for arbitrary initial conditions remains elusive.
Here, we present a spectral-theoretic approach to finite time-average statistical mechanics of ergodic systems. In mathematical terms we focus on the statistics of bounded, local additive, and potentially lowerdimensional, functionals of normal ergodic Markovian stochastic processes with continuous and discrete statespaces. The paper is organized as follows. We first provide in Sec. II a brief introduction into time-average statistical mechanics. In Sec. III we re-derive the wellknown Feynman-Kac formulas for Markovian diffusion with Itô calculus. In Sec. IV.A spectral theory combined with non-Hermitian perturbation theory is then applied to obtain our main result -exact expressions for the mean, variance and correlations of time-average observables for any non-stationary preparation of the system expressed explicitly in terms of the eigenspectrum of the underlying forward and backward generator of the dynamics, which may correspond to Fokker-Planck diffusion or Markovian jump dynamics governed by master equation. We demonstrate explicitly the emergence of a central limit law in a spectral representation on large deviation time-scales. We particularly derive in Sec. IV.C a general upper bound on fluctuations in terms of an integral of the, generally non-Markovian, return probability that is valid for generators of overdamped dynamics obeying detailed balance, which is the second main re- sult of this article. Finally, simple analytically solvable examples are provided in Sec. V to demonstrate how to apply the theory. We conclude in Sec. VI.

II. TIME-AVERAGE STATISTICAL MECHANICS
Traditional (classical) ensemble statistical mechanics describes physical observations as averages over individual realizations of the dynamics at single (or multiple) pre-determined times. For example, the ensemble average of an observable V (x t ) at a time t for an ergodic stochastic process x τ (0 ≤ τ ≤ t) starting from some non-stationary initial condition p 0 (x 0 ) is defined by where Ω is the state space of the process and P t (x|x 0 ) is the so-called propagator, i.e. P t (x|x 0 )dx (upper line) is the probability that the process is found in x ∈ Ω within the increment dx at time t given that it started at t = 0 at x 0 . Note that if x is continuous valued (x ∈ Ω ⊂ R d ; see Fig. 1 solid red line) then P t (x|x 0 ) is a probability density, whereas if x is from a discrete state space Ω ( Fig. 1 solid blue line) the integral in Eq. (1) (upper line) becomes a sum Ω dx → x∈Ω and P t (x|x 0 ) becomes a plain probability as shown in the lower line of Eq. (1). If x 0 is sampled and averaged over a stationary (invariant) measure p(x 0 ) = P inv (x 0 ) or t becomes sufficiently (i.e. ergodically) large P ∞ (x|x 0 ) ≡ P inv (x) the ensemble average becomes time-independent Conversely, in single-molecule dynamics, singleparticle tracking and other related experiments one probes individual realizations of x τ within the interval 0 ≤ τ ≤ t and instead analyzes the observation by taking a time-average. Such time-average observables are in general random, fluctuating quantities with non-trivial statistics. For example, for a physical observable V (x t ), which may correspond e.g. to the squared displacement [61,62] or local time [51,59,60,63] in single particle tracking or the FRET efficiency [5][6][7][8] or distance between two optical traps [9][10][11][12][13] in single-molecule fluorescence and force spectroscopy, respectively, the (local) time-average is defined as and depends on the entire history of x τ until time t (see also dotted lines in Fig. 1). The statistical evolution of V t is therefore a non-Markovian process with a probability density function that V t = ψ [46,52,53,59,60] where δ(z) is the Dirac delta function and · x0 denotes the average over all paths starting at x 0 , i.e. p 0 (x) = δ(x − x 0 ), and propagating until time t. The corresponding result for arbitrary initial conditions p 0 (x 0 ), i.e. P t (ψ|p 0 ), can be obtained from P t (ψ|x 0 ) as shown in Sec. IVb. The empirical "probability density function" θ x (t) that estimates P inv (x) along a single trajectory in time-average statistical mechanics is the so-called local time fraction defined as [59,60] which allows to rewrite the time average (3) in the form where δ(x − x t ) denotes the Dirac delta function if x ∈ Ω is continuous, whereas δ(x − x t ) denotes the Kronecker delta if x ∈ Ω is integer valued. Because the dynamics x t is ergodic we have lim t→∞ θ x (t) = P inv (x) [52,59,60] and Moreover, on the so-called large-deviation time-scale, i.e. on the time-scale that is finite but longer that the longest relaxation time of x t , we have convergence in the mean, V tLD = V inv , and Gaussian fluctuations around the mean value [45,59,60,64,65]. For finite, and in particular sub-ergodic (i.e. supra-large deviation) times the statistics of V t is, however, non-trivial. In order to exploit the role of the local time fraction θ x (t) as a "propagator" in time-average statistical mechanics (via Eq. (6)) we now relate the statistics of θ x (t) to the probability density of a general time-averaged physical observable V t defined in Eq. (3). The characteristic function of V t , i.e. if V t ≥ 0 is non-negative the Laplace transform is defined as where we inserted Eq. (6) in the last. Eq. (8) relates the statistics of the time average V t to the statistics of all local time fractions θ x (t) (x ∈ Ω). We note that Eq. (8) must be modified if V t can also be negative such that a Fourier transform is required instead, which amounts to replacing u → iω with ω ∈ R in Eq. (8). The probability density is obtained from Eq. (8) by Laplace inversion where c ∈ R lies to the right of all singularities ofP t (u|x 0 ) and we have assumed that P t (ψ|x 0 ) is of exponential order for sufficiently large ψ (in the following section the conditions on V (x) will be made more precise). In case of negative support Eq. (9) becomes the inverse Fourier transform. Here we are particularly interested in fluctuations of time averages V t and W t of different physical observables V (x τ ) and W (x τ ) which are quantified by [59,60] where σ 2 V (t) = C V V (t) denotes the variance of V t and C V W (t) the covariance between V t and W t . Using Eq. (8) we obtain more generally where x i , y j ∈ Ω, respectively. Thereby, the ensemble average corresponding to the last term in Eq. (11) is obtained by differentiating the characteristic function with respect to the Laplace (or Fourier) variable.
It therefore follows, that the fluctuations and (linear) correlations of general time-average observables are fully specified by multi-point correlations functions of the local time fraction. These are derived on the basis of the Feynman-Kac formalism, which is presented below for continuous diffusion processes. Extension to discrete state dynamics is discussed afterwards.

III. FLUCTUATIONS OF ADDITIVE FUNCTIONALS
A. Itô approach to Feynman-Kac theory of additive functionals It seems to be customary in the physics literature to start from a path integral approach to Feynman-Kac theory [46] and in the following to derive a backward Fokker-Planck equation for the characteristic function [52,53].
Here we provide a simple derivation of the "forward" Feynman-Kac theory based on Itô calculus.
We consider a d-dimensional Markovian diffusion x t ∈ Ω ⊂ R d process in the presence of a drift F(x) driven by a d-dimensional Gaussian white noise governed by the Itô equation where σ is a d×d noise matrix such that D = σσ T /2 becomes a symmetric positive (semi)definite diffusion matrix, dW t is an increment of a d-dimensional Wiener process, such that W t = 0 and dW t,i dW t ,j = δ(t − t )δ ij dt. We assume throughout that F(x) is sufficiently confining to assure that the process x t is ergodic with a steady state probability density P inv (x). Multiplying the time average in Eq. (3) by the trajectory length t we transform the time average to the additive functional We now consider the joint process of x t and ψ t . According to Itô's lemma any twice differentiable function f (x, ψ) with x t and ψ t defined by Eqs. (12) and (13), respectively, satisfies where we inserted the diffusion matrix D = σσ T /2 and for the last term we used dψ t = V (x t )dt that follows from Eq. (13). Using Itô's lemma (14) we derive in the following the time evolution of the joint probability density Q t (x, ψ|x 0 ) to find the system in state x and the functional value ψ at time t given that the system started from x 0 . For convenience we first focus on positive functionals (ψ ≥ 0). Using a test function that vanishes at the boundary f (x, 0) = 0 we obtain after some calculations [66] which is obtained as follows.
In the first line of Eq. (15) we differentiate both sides of the identity To obtain the second line we inserted Itô's Lemma (14) and finally performed an integration by parts. Since Eq. (15) holds for any function f that vanishes at the boundary ψ = 0 we obtain where we have defined the forward generatorL = ∇ x · D∇ x − ∇ x · F i (x) and further introduced a boundary term V (x)δ(ψ)Q t (x, 0|x 0 ) that vanishes for ψ > 0 and is derived in the following two steps. First, Eq. (15) holds for all functions f (x, ψ) with f (x, 0) = 0, which immediately gives Eq. (16) without the last term for ψ > 0 (see also [66]). Second, the last term in Eq. (16) is required to guarantee the conservation of probability ∞ 0 dψ Ω dx∂ t Q t (x, ψ|x 0 ) = 0, i.e., to correct for the fact that there is a non-zero probability that the functional has a vanishing value ψ = 0. Finally, performing a Laplace transform of Eq. (16)Q t (x, u|x 0 ) ≡ ∞ 0 e −uψ Q t (x, ψ|x 0 )dψ we obtain the forward Feynman-Kac partial differential equation for the characteristic function of the joint density of position and ψ where Q t (x, ψ|x 0 ) is the central object of the 'forward' Feynman-Kac approach [15,46].
We now relax the assumption by allowing for a negative support of ψ t , which we denote explicitly by ψ t → Ψ t . In this case we need not to make any additional assumptions on f (x t , Ψ t ) for Ψ t = 0 because naturally lim |Ψ|→∞ Q t (x, ψ|x 0 ) = 0. The lower boundary of integration over Ψ in Eq. (15) is extended to −∞ and boundary terms therefrom resulting from the partial integration vanish. Due to the boundary conditions on Q t (x, Ψ|x 0 ) the resulting Eq. (15) implies which upon taking a Fourier transformQ t (x, ω|x 0 ) ≡ ∞ −∞ e −iωΨ Q t (x, Ψ|x 0 )dΨ leads to the forward Feynman-Kac partial differential equation Note that the generalization to the case of a joint problem for multiple functionals ψ i t with i = 1, . . . , p is straightforward. Introducing the vectorial notation ψ t ≡ t 0 V(x τ )dτ with the corresponding Laplace images u the resulting equation reads which allows for the computation of higher-order correlation functions (11). What we actually seek for is the marginal probability density of ψ given x 0 defined by which is the statement of the Feynman-Kac theorem [15]. Note that the corresponding characteristic func-tionP t (u|x 0 ) is the solution of the 'backward' Feynman-Kac problem [52,53]. Moreover, the marginal probability density of x given x 0 being the the plain propagator To convert the probability density of ψ t (Ψ t , respectively) to the density of the time-average observable ψ t we perform a straightforward change of variables. Therefore, if ψ t is the local time tθ x (t) then the density of the local-time fraction is P t (θ|x 0 ) ≡ t Ω Q t (x, tψ|x 0 )dx and its characteristic function is given byP t (u|x 0 ) =

B. From the forward to the backward Feynman-Kac equation
For convenience we henceforth adopt the bra-ket notation, where the "ket" |h denotes a vector, the "bra" the integral operator g| ≡ Ω dxg † (x), and the scalar product is defined as g|h ≡ Ω dxg † (x)h(x). Introducing, moreover, the "flat" state |− ≡ Ω dx|x and −| ≡ Ω dx x|, Eqs. (19) and (20) for a general initial condition p 0 (x), i.e. |p 0 = Ω dx 0 p 0 (x 0 )|x 0 and To arrive at the second equality we have used Green's identity, introduced the adjoint (or backward) Fokker- and used that the Laplace transform of a real function f (t) transforms asf (s † ) =f † (s) under complex conjugation. In the following subsection we show that for Markov jump processes (22) the theory can be adopted one-by-one.

C. Markov-jump dynamics and additive functionals
Markov jump processes correspond to a discrete state space Ω in which the system jumps with a constant rate w xy from state x ∈ Ω to another state y ∈ Ω, such that the propagator P t (x|x 0 ) satisfies the master equation where x|L|y = w yx if x = y and x|L|x = − y =x y|L|x . According to the celebrated Gillespie algorithm [67] a single trajectory x τ (0 ≤ τ ≤ t) consist of a sequence exponentially distributed local waiting time in a given state x ∈ Ω before transiting instantaneously to another state. After the jump has occurred the process starts a new. Denoting the number of transitions from state x to state y by n xy (t) and identifying the sum of all local waiting times in x up to time t by tθ x (t) the path probability (or path weight) (0 ≤ τ ≤ t) generated by the Markov dynamics (23) can be written as [37] Furthermore, translating the integral into a sum in the time average from Eq.
where in the last term truncated the diagonal of the generator x|L(u)|x ≡ x|L|x − uV (x)/t, whereas kept the off diagonal elements unchanged, that is, x|L(u)|y ≡ x|L|y if x = y. SinceL and V (x) are real, Eq. (22) therefore holds also for Markov-jump processes. In the following we derive a spectral theory, which unifies diffusion processes and Markov jump processes.

D. Spectral theory of non-Hermitian generators
We henceforth employ a spectral-theoretic approach and are thus required to make some more specific assumptions about the underlying dynamics in order to assure that the generatorL is diagonalizable. An excellent account of the theory for Markov-jump dynamics gov-erned by a discrete-state master equation can be found in [68]. In the case of Fokker-Planck dynamics we consider that x t is an ergodic Markovian diffusion evolving according to Eq. (12) with the drift field F(x) not necessarily corresponding to a potential field (which thus include systems with a broken detailed balance) but at the same time requiring that it is sufficiently confining, that is, it grows sufficiently fast as |x| → ∞ to assure thatL has a pure point-spectrum. Moreover, we require thatL is diagonalizable [69]. A more detailed mathematical exposé of the requirements for, and properties of,L can be found in [60]. In all practical examples presented below we will in fact assume that the dynamics is overdamped and thatL obeys detailed balance [66,70] implying that it is orthogonally equivalent to a self-adjoint operator and hence automatically diagonalizable.
Let −λ k (Re(λ k ) ≥ 0), L k |, and |R k denote the eigenvalue and orthonormal left and right eigenstates ofL, and −λ † k , R k | and |L k the corresponding orthonormal eigen- Then written in the respective eigenbasesL andL † read with the ground-state eigenvalue λ 0 = 0 and corresponding null-space |R 0 ≡ |P inv and L 0 | ≡ −|. In the respective dual eigenbasis the propagator where

IV. CHARACTERISTIC FUNCTION NEAR ZERO VIA NON-HERMITIAN PERTURBATION THEORY
BecauseL is in general not self-adjoint we need to separately perturb left and right eigenstates. Based on Eqs. (8) and (11) we only require the moment-generating function (22) in the limit |u| → 0 to calculate moments of arbitrary order. To keep the treatment generally applicable we want to utilize the spectral expansion ofL (L † , respectively). We therefore aim to diagonalize the "tilted" propagator in Eq. (22) in the limit when u vanishes.
First we must confirm that the tilted propagator The singularities of Eq. (28) correspond to the perturbed eigenvalue spectrum {−λ k (u)} ofL(u)) and diagonalizability is broken whenever one or more singularities are not simple poles (see e.g. [72]). Eq. (22) shows that |u| = 0 is not an accumulation point. Moreover,L has a pure point spectrum therefore an arbitrarily small |u| cannot cause the emergence of poles of second order in Eq. (28) that would break diagonalizability, akin to the "avoided crossing theorem". Therefore, in the limit |u| → 0 the tilted generatorL(u) is diagonalizable, u can be taken as real [73], and the eigenspectrum ofL(u) corresponds to a regular perturbation of the original eigenvalue problemL|R k = −λ k |R k (L † |L k = −λ † k |L k ) and we seek for a perturbative expansion of the tilted eigenspectrum, e.g. ofL † (u): Note that while the spectra ofL andL † are complex conjugates (see Eq. (26)) the perturbation is in fact symmetric (see Eq. (28)). Therefore, the spectra ofL(u) and L † (u) are not complex conjugates except for the unperturbed part, which we denoted in Eq. (29) by parenthesis λ † . According to Eq. (29) multiple functionals yield additive perturbations. It thus suffices to carry out the calculations for u → u and write the corresponding general result by inspection. We are interested in up to second order moments (11) and therefore need to evaluate the perturbation up to second order in u: The 0-th order of the expansion gives the solution of the unperturbed system. For the higher orders we need to solve Eqs. (31) matching terms of equal order. Introducing the coupling elements V lk ≡ R l |V |L k we obtain (details of the calculation are shown in Appendix A) However, while they are orthogonal by construction, the resulting perturbed eigenstates are not normalized anymore, i.e. R k (u)|L k (u) = 1. Hence, we need to postnormalize them such that where from it follows that We now use the perturbative 2 nd -order eigenspectrum to diagonalize the moment-generating function in Eq. (22) where moreover (39) and p 0 |L k (u) = p 0 |L k +u p 0 |L 1 k +u 2 p 0 |L 2 k +O(u 3 ), where |L 1 k and |L 2 k are given by Eq. (34). In the special case when initial conditions x 0 are drawn from the steady state probability density, i.e. p 0 (x) = P inv (x), these equations simplify to Turning now to the case where V (x) extends negative values we must simply replace u by iω. In order to obtain moments up to second order (i.e. Eq. (10) we are now left with evaluating derivatives with respect to u.

A. Mean value, fluctuations and correlation functions for general initial conditions
We first focus on the case of a single time-average observable ψ t ≡ ψ t /t (see Eq. (13)). We will only present the result in terms of the spectrum of the backward gen-eratorL † since those corresponding toL follow trivially from Eq. (22). The mean value for an arbitrary initial condition p 0 (x) (note that p 0 |− = 1 for any normalized initial condition), ψ t p0 = −∂ uPt (u|p 0 )| u=0 is found to be given by with the anticipated ergodic limit V 00 = ψ t→∞ p0 = P inv |V |− ≡ Ω V (x)P inv (x)dx. The result (41) is equally valid in cases where V (x) can become negative (as long as its is bounded). Moreover, from Eq. (41) is easy to discern the large deviation asymptotic and in the special case of steady-state initial conditions p 0 (x 0 ) = P inv (x 0 ) Eq. (41) reduces to the timeindependent ergodic result ψ t inv = V 00 . The second moment follows from ψ 2 t p0 = ∂ 2 uPt (u|p 0 )| u=0 and reads yielding the variance and we further introduce the following notational convention for localized initial conditions p 0 ( . From Eqs. (41), (43), and (44) follows the anticipated ergodic result ψ 2 t→∞ p0 = V 2 00 , proving that in the ergodic limit ψ t becomes deterministic (i.e. the variance vanishes, σ 2 ψ,p0 (t → ∞) = ψ 2 t→∞ p0 − ψ t→∞ 2 p0 = 0). Conversely the largedeviation asymptotic reads, independent of p 0 (x) embodying the emergence of the central limit theorem.
One can further show that all higher cumulants decay to zero faster that 1/t, and since the only distribution with a finite number of non-zero cumulants is the Gaussian distribution [74], the large deviation mean value (42) and variance (46) specify the entire asymptotic probability density for tReλ † 1 1. In the special case of steadystate initial conditions, p 0 (x) = P inv (x), we instead find the variance (see also [59]) Note that for overdamped systems in detailed balance we have λ † k ∈ R and moreover V 0k = V k0 ∈ R, which implies (compare Eqs. (46) and (47)) that the fluctuations for stationary initial conditions are bounded from above by We now inspect the correlation between two functionals ψ 1,t and ψ 2,t (see Eq. (10)) defined as C ψ 1 ψ 2 (t) = ψ 1,t ψ 2,t p0 − ψ 1,t p0 ψ 2,t p0 . The mean values were derived in Eq. (41) so we only require the mixed second moment ψ 1,t ψ 2,t p0 , which is obtained from the joint moment-generating function (i.e. generalization of Eq. (37) to two variables) as ψ 1,t ψ 2,t p0 = ∂ 2 uvPt (u, v|p 0 )| u=v=0 . A lengthy calculation leads, upon introducing the coupling elements U kl ≡ R k |U (x)|L l and the shorthand notation W klmn = V kl U mn + U kl V mn , to the exact result ψ 1,t ψ 2,t p0 = W 0000 2 and we note that ψ 1,t inv ψ 2,t inv = W 0000 /2 implying that for an ergodic system any two functionals asymptotically decorrelate, lim t→∞ C ψ 1 ψ 2 (t) = 0. Equations In the case of stationary initial conditions p 0 (x) = P inv (x) we have that P inv |L k = δ k0 and as a result of Eq. (48) the covariance reduces to (note that W k00k = W 0kk0 and W 0000 /2 = ψ 1,t inv ψ 2,t inv ) Finally, in the large deviation regime we recover lim t 1/λ1 the 1/t scaling reflecting the emergence of the central limit theorem. Therefore, it follows that an arbitrary set of m functionals ψ t in the large deviation limit exhibits Gaussian statistics. If we denote the vector of mean values as ψ inv and introduce the symmetric covariance matrix C with diagonal elements C ii = t lim t 1/Reλ1 σ 2 ψ i (t) (see Es. (46)) and off-diagonal elements C ij = t lim t 1/Reλ1 C ψ i ψ j (t) (see Es. (50)) then ψ t obeys the asymptotic Gaussian limit We now introduce the rescaled variablesψ ≡ ψ t √ t and scaled mean µ ≡ ψ t / √ t inv , which upon renormalization lead to a time-independent density. Moreover, we define with the shifted mean and stretched variancẽ Then the limit law (51) implies that the univariate large deviations Q LD t (ψ) and conditional bivariate large deviations Q LD t (ψ 1 |ψ 2 ) ≡ Q LD t (ψ 1 ,ψ 2 )/Q LD t (ψ 2 ) collapse, upon rescaling, onto a universal Gaussian master curve where N x (0, 1) denotes the Gaussian probability density with zero mean and unit variance. The explicit rescaling (i.e. Eq. (54) combined with Eq. (46) and Eq. (50)) leading to the collapse onto a master unit normal density in the large deviation limit are the main practical consequence of our large deviation result.

B. Degenerate eigenspectra
Note that if the spectrum ofL has degenerate eigenstates (such as e.g. in single-file diffusion [59,60]) special care is required for initial conditions that do not correspond to the steady state, i.e. p 0 (x) = P inv (x), as a result of the singularities the degeneracy causes in Eqs. (43) and (48). As it is customary in regular perturbation theory (see e.g. [75]) one must first post-diagonalize all the respective degenerate sub-spaces prior to using Eqs. (43) and (48). Once this has been taken care of (in any of the many possible methods [76]) and the degenerate eigenstates are replaced by their appropriate linear combinations, Eqs. (43) and (48) can be used as they stand.

C. A general upper bound for overdamped reversible dynamics
WhenL corresponds to reversible overdamped dynamics (i.e. D −1 F(x) = −β∇ x U (x) is a gradient field), or to a reversible Markov jump process (i.e. the transition matrix elements in Eq. (23) satisfy the symmetry y|L|x / x|L|y = e βU (x)−βU (y) ) the large deviation asymptotic (46) provides an upper bound for fluctuations at any time. Let us define the projection operator Γ x (ψ; V ) ≡ Ω dxδ(y − x) [60] as well as the (generally non-Markovian) joint probability density of starting from ψ and returning to the initial value ψ in an ensemble of trajectories x t starting from the equilibrium probability density P inv (x 0 ) = P eq (x 0 ) as [60] G eq t (ψ, ψ) ≡Γ x (ψ; V )Γ x0 (ψ; V )P t (x|x 0 )P eq (x 0 ). (55) Then Eq. (46) and the spectral decomposition of P t (x|x 0 ) in Eq. (27) imply the general upper bound where the inequality saturates as t → ∞. Eq. (56) enables to obtain an upper bound on fluctuations of ψ t from the integral over the return probability (55) and thus bounds the fluctuations of time-averaged observables ψ t (with ψ t of type (13)) by means of the ensemble propagator. Eq. (56) is the main practical result of this work. Interestingly, we found a similar bound involving the integral of the return probability in the study of large deviation asymptotics of the first passage functional [77].

D. Physical interpretation of the results
We now provide some intuition about the developed theory. As time evolves the value of ψ t for an ergodic process eventually becomes only weakly correlated, and the statistics of ψ t passes first through the large deviation regime (51) where the central limit theorem kicks in with Gaussian statistics, and finally ends up in Khinchin's law of large numbers where it becomes deterministic and equal to ψ t inv [78]. For simplicity we start in the large deviation regime (46).
By using spectral theory we map fluctuations of ψ t onto the eigenmodes ofL (and/orL † respectively), with the 'similarity' to a given eigenmode reflected by the overlap elements V 0k , V k0 . Since on these time-scales all memory of the initial condition is lost, which is equivalent to imposing stationary initial conditions, only overlaps from and to the ground state are relevant. Moreover, due to the orthogonality of eigenmodes these projections are statistically independent. Each eigenmode has a finite lifetime or correlation time 1/λ k , therefore, in a time t λ −1 k any k-th projection acts as shot-noise and there will be tλ k independent realizations of such a projection reducing the (co)variance by a factor 1/tλ k (see Eqs. (46) and (50)). In the limit t → ∞ the Gaussian converges to a Dirac delta lim t→∞ Q LD t (ψ) = δ(ψ − µ). At shorter times non-trivial corrections arise to these large deviation results due to strong correlations between the values of ψ t at different times t. As a result of these correlations the 'completely decorrelated' large deviation results in Eqs. (46) and (50)) becomes reduced by a term that seems to reflect the 'effective probability of mode k to persist until t', t −1 t 0 e −λ k τ dτ = (1 − e −λ k t )/λ k t (see Eqs. (47) and (49) as well as Eqs. (B3)and (B4)). In the case of general initial conditions p 0 (x) further terms arise (see Eqs. (45) and (48)) that reflect the memory of the initial condition. These terms, however, are difficult to interpret beyond the point that they reflect projections coupling different excited eigenstates and thus describe fluctuation-modes that are more complicated than simple excursions starting and ending in the steady state.

V. APPLICATIONS OF THE THEORY
We now apply the theory to a collection of simple illustrative problems. Due to the fundamental role played by the local time fraction θ x (t) and because it determines the dynamics of other time-average observables (see Eq. (6)) we focus on θ x (t) alone. The coupling elements are therefore simply given by V lk = Ω dyR l (y)δ(x−y)L k (y) ≡ R l (x)L k (x). We first present explicit results for local times for continuous space-time Markovian diffusion processes and an irreversible (i.e. driven) three-state uni-cyclic network. Next, we apply the theory to a simple two-state Markov model of the Berg-Purcell problem [28][29][30]36], i.e. the physical limit to the precision of receptor-mediated measurement of the concentration of ligand molecules.
As minimal, exactly solvable models of continuousspace Markovian diffusion we consider a Wiener process confined to a unit interval with reflective boundaries and the Ornstein-Uhlenbeck process. To demonstrate the theory for Markov-jump dynamics we consider a random walk in a finite harmonic potential and a simple 3-state uni-cyclic network.

A. Local time fraction of the Wiener process in the unit interval
The propagator of the Wiener process confined to a unit interval (i.e. L = 1) is the solution of x in a unit interval are given by λ W k = k 2 π 2 (we express time is expressed in units of τ = L 2 /D) and the eigenvectors read [66,71], x is self-adjoint. The mean local time-fraction, θ x (t) , the variance σ 2 θ (t) and covariance C θ1θ2 (t) for the confined Wiener process are shown in Fig. 2. In the case of equilibrium initial conditions θ x (t) is constant and equal to P inv (x), and the fluctuations of θ x (t) are largest at the boundaries as a result of repeated collisions with the walls. Notably, starting from localized conditions θ x (t) a a function of x, in contrast to the ensemble propagator P W t (x|x 0 ), displays a pre-existent cusp located at the initial condition x 0 (see Fig. 2 c). The fluctuations of θ x (t) are larger near the initial condition and at the boundaries. Note that the fluctuations are always larger for equilibrium initial conditions (compare Fig. 2b and Fig. 2d).

B. Local time fraction of the Ornstein-Uhlenbeck process
Trajectories of the one-dimensional Ornstein-Uhlenbeck process are solutions of the Itô equation and on the level or probability density correspond to the Fokker-Planck equation (∂ t − D[∂ 2 x + γ∂ x x])P OU t (x|x 0 ) = 0 with initial condition P OU 0 (x|x 0 ) = δ(x − x 0 ) and natural boundary conditions lim |x|→∞ P OU t (x|x 0 ) = 0. To connect continuous processes to discrete ones we translate the Fokker-Planck equation of the Ornstein-Uhlenbeck process into a random walk on a lattice with spacing ∆x and the harmonic potential γx 2 entering transition rates according to in a confined domain Ω conf = {−l, −l + ∆x, . . . , l − ∆x, l} ⊂ Ω. The matrixL is tri-diagonal and satisfies y y|L|x = 0 for all x ∈ Ω and x ∈ Z∆x. We diagonal-izedL numerically using the library [79]. The mean, variance and correlation function for the continuous-space Ornstein-Uhlenbeck process (59) obtained from Brownian dynamics simulations are depicted in Fig. 3 (symbols) and are in excellent agreement with the spectraltheoretic results for the corresponding lattice random walk approximation (60) (lines). In Fig. 3f we also investigate the full probability density function of the fraction of occupation time in the interval x ∈ [0, 0.01], i.e. ψ t = t −1 t 0 1 [0,0.01] [x(τ )]dτ , on large deviation time scales and compare it the the theoretical Gaussian prediction Eq. (46). Note that while the eigenspectrum of the generator of continuous Ornstein-Uhlenbeck dynamics is unbounded implying that the spectral-theoretic result would require the summation of a large number of terms, the summation in the lattice approximation is limited by the number of lattice points. Therefore, except for very short times, where the lattice approximation naturally breaks down, this example demonstrates that our formalism applies equally well for Markov jump processes and diffusion dynamics.

C. Local time fraction in a driven uni-cyclic network
Let us in the following address a simple 3-state model with broken detailed balance to also address driven systems. The model corresponds to a simple cycle with states 1, 2 and 3, where all rates in a given direction are equal but each has the same forward/backward asymmetry. The model may represent, for example, a molecular motor such as the F1-ATPase driven by ATP hydrolysis [80]. The corresponding transition matrix of the model readsL and has eigenvalues λ 0 = 0, λ 1,2 = −9/2 ± i √ 3/2 and eigenvectors |R 0 = 3 −1 (1, 1, 1) T and As a result of broken detailed balance the eigenspectrum is complex. In Fig. 4 we analyze the mean (panel a), fluctuations (panel b) and correlation function (panel c) of the local time fraction θ x (t) in the various states for non-equilibrium steady-state initial conditions (light blue lines) and conditions initially localized in state x 0 = 2, i.e. |p 0 = (0, 1, 0) T . The theoretical result (lines) show an excellent agreement with simulations (symbols) carried out using the Gillespie algorithm [67]. We also confirm the Gaussian statistics of the local time fraction ψ t from Eq. (51) in Fig. 4d.

Generic behavior of local time fraction in ergodic systems
Note that an exhaustive study of the statistics of the local-time fraction is beyond the scope of this work. Nevertheless, we discuss some here general features of θ x (t). The manner in which θ x (t) , starting from some nonequilibrium initial conditions, approaches the ergodic invariant measure P inv (x) can be highly non-trivial and even non-monotonic (see e.g. Figs. 5a). Even when the ergodic limit is reached, where the variance ceases to depend on time, i.e. tσ 2 θx (t) = f (t), the fluctuations display a non-trivial behavior (see Fig. 5). For example, in the case of the Wiener process fluctuations are enhanced close to the boundaries, while for the Ornstein-Uhlenbeck process they become depressed near the minimum. Both results may be interpreted in terms of random 'oscillations' around a typical position and confined by a boundary amplifying fluctuations.
Moreover, the time-dependence of θ x (t) for nonstationary initial conditions is often non-monotonic or has a non-monotonic derivative (see Fig. 6a and c). A comparison between θ x (t) starting from stationary (dashed lines) and localized (full lines) initial conditions illustrates the two coexisting decorrelation mechanisms of θ x (t) at different times, one corresponding to selfaveraging and emergence of the central limit theorem (compare dashed and dotted lines), the other in addition to forgetting the memory of the initial condition (full lines). Stationary initial conditions always give rise to larger fluctuations than non-stationary initial conditions (compare dashed and full lines in Fig. 6b and d), and in the particular case of equilibrium initial conditions for systems obeying detailed balance σ 2 θ (t) is a monotonically decaying function of time t (see Eq. (47) and Fig. 6b and d) with an upper bound given by the large deviation asymptotic (see Eq. (56) with a ∝ 1/t scaling dictated by the central limit theorem (dotted lines in Fig. 6b and  d).
The covariance of the local time fraction between a pair of points x 1 and x 2 in the continuous setting (Figs. 2e and 3e) or x = 1 and x = 2 in the discrete setting (Fig. 4c), C θ1θ2 (t), displays a similarly non-trivial and non-monotonic dependence on time and initial conditions  p 0 (x 0 ) as shown in Figs. 2e, 3e, and 4c. The striking dependence on the tagged points reflects a directional persistence of individual trajectories in-between said points and can therefore be used as a robust indicator of directional persistence and thus 'temporally ordered exploration' on the level of a single trajectory.

Universal asymptotic Gaussian limit law for time-average physical observables
Finally, we comment on the universal asymptotic Gaussian limit law Eq. (51) for Markovian as well as non-Markovian time-average physical observables of type (3) of ergodic stochastic dynamics of the form given in Eqs. (12) and (23). Namely, using the asymptotic results (42), (46), and (50) in the large-deviation probability density function (51), and rescaling to the centered and time-independent variables Ξ and Ξ i|j defined in Eqs. (52) and (53), we can rescale the probability density of any time-average physical observable Q LD t (ψ i ), and the conditional probability density of a time-average physical observable given another time-average physical observable Q LD t (ψ i |ψ j ), to collapse at long times onto a unit normal probability density (54). For the three models studied here, Figs. 2f, 3f, and 4d, and additionally for the conditional probability density function of occupation time fraction in x ∈ [0.1, 0.4] given occupation time fraction in y ∈ [0.6, 0.9] for the Wiener process, we demonstrate this collapse explicitly in Fig. 7.

D. Precision limit of concentration measurement by a single receptor
Let us now investigate the physical limit to the precision of concentrations measurements by means of the simplest two-state Markov jump process with states Ω = {0, 1} [28]. The receptor can either be occupied by a ligand (x = 1) or be empty (x = 0). Let the background ligand concentration be c and assume that the ligand binds with a rate kc and unbinds with rate k, ignoring for simplicity any spatial variations of concentration. The generator and its eigenvectors are given bŷ being the only non-zero eigenvalue. The left eigenvectors corresponding to Eq. (63) are L 0 | = −| = (1, 1) and L 1 | = (c, −1). Moreover, since the entire state space has only two states we have θ 0 (t) = 1 − θ 1 (t), which implies that the mean values of the respective local time fractions are given by θ 1 (t) = (1 + c) −1 c and θ 0 (t) = (1 + c) −1 , assuming that the system was initially in equilibrium |p 0 = |R 0 .
If the receptor estimates the concentration c by reading out and averaging the fraction of time the ligand is bound, θ 1 , over an interval of duration t, the precision of the estimate is bounded from above by the variance of the local time fraction given by Eq. (47) and reads explicitly Typically one assumes that the measurement t is longer than any correlation time [28,30,81,82], which in the present setting implies t 1/[k(1 + c)], i.e. much longer than the correlation time of two-state Markov switching noise, τ c ≡ λ −1 1 = 1/(k + kc) [30]. In this regime the averaging noise corresponds to shot-noise such that the variance decreases as with the number of statistically independent receptor measurements # t [28,30,81,82], where # t ∼ t/τ c is the number of statistically independent realizations of the two-state process and therefore σ 2 θ1 (t) ∝ 1/# t = τ c /t, according to the central limit theorem.
Based on the bound derived in Eq. (56) the shotnoise limit is in fact an upper bound to fluctuations of receptor occupancy at any duration of measurement, and saturates only in the limit t τ c . Namely, a direct application of the bound (56) indeed yields, using G eq t (1, 1) = P t (1|1)P ∞ (1|1), implying, according to Eq. (56) Therefore, for short, and particularly finite measurements the shot-noise limit of fluctuations for long receptor read-out [28,30,81,82] gives only an upper bound to the uncertainty of the estimate, whereas the inequality becomes sharp at long times. Fundamental bounds on the precision of inferring c from θ 1 (t) can be found in [28]. Using a maximum likelihood estimate of c given the entire time trace it has been found that the of c from a long single measurement can reduce the error of the estimate of c by a factor of two [31,32]. A detailed discussion of the precision of inferring kinetic parameters by means of non-local functionals can be found in [37].

VI. CONCLUDING PERSPECTIVE
We developed a general spectral-theoretic approach to time-average statistical mechanics, i.e. to the statistics of bounded, local additive functionals of (normal) ergodic stochastic processes with continuous and discrete statespaces. In particular, we have shown how to obtain exactly the mean, variance and correlations of time-average observables from the eigenspectrum of the underlying forward or backward generator. We re-derived the famous Feynman-Kac formulas using Itô calculus and included a brief derivation for Markov-jump processes. We combined Feynman-Kac formulas with non-Hermitian perturbation theory to derive an exact spectral representation of the results. We demonstrated explicitly, and quantitatively, the emergence of the universal central limit law in a spectral representation on large deviation time-scales. For the special case of equilibrated initial conditions and dynamics obeying detailed balance we derived a general upper bound on fluctuations of timeaverage observables. We discussed our theoretical results from a physical perspective and provided simple but instructive practical examples to demonstrate how the theory is to be applied. Our work is applicable to continuous as well as discrete state-space processes, reversible as well as irreversible, encompassing a wide and diverse range of phenomena involving time-average observables and additive functionals in physical, chemical, biological systems as well as financial mathematics and econophysics.

Mean, fluctuations and correlations
We now derive ψ t , σ 2 ψ (t) and C ψ 1 ψ 2 (t) presented in Eqs. (41)-(49) using the Dyson series. Starting form Eq. (B1) we can carry out all integrations analytically for arbitrary initial conditions p 0 (x). To first order in u we obtain where for p 0 (x) = P inv (x) we have L k |p 0 = δ k0 leading to ψ inv = V 00 = P inv (x).
To second order in u we find for an arbitrary p 0 (x) when k = l = 0 only V 2 00 /2 survives, while for k = 0 and l = 0 we find Conversely, when k = 0 and l = 0 we end up with while for k = l and k = 0 and l = 0 we have Finally, when k = l = 0 we obtain 1 t 2 k>0 l>0,l =k V k0 V lk λ k − λ l L l |p 0 1 − e −λ l t λ l − 1 − e −λ k t λ k (B7) The sum of these terms yields the result sought for. The corresponding result for stationary initial conditions, p 0 (x) = P inv (x), is obtained using L l |P inv = δ l0 , which leads to Eq. (47) and (49). When considering correlations we make the replacement uV → uV 1 + vV 2 and replace ∂ 2 u → ∂ 2 uv to compute the covariance. The formulas above thereby generalize in a straightforward manner.