Unfolding Tagged Particle Histories in Single-File Diffusion: Exact Single- and Two-Tag Local Times Beyond Large Deviation Theory

Strong positional correlations between particles render the diffusion of a tracer particle in a single file anomalous and non-Markovian. While ensemble average observables of tracer particles are nowadays well understood, little is known about the statistics of the corresponding functionals, i.e. the time-average observables. It even remains unclear how the non-Markovian nature emerges from correlations between particle trajectories at different times. Here, we first present rigorous results for fluctuations and two-tag correlations of general bounded functionals of ergodic Markov processes with a diagonalizable propagator. They relate the statistics of functionals on arbitrary time-scales to the relaxation eigenspectrum. Then we study tagged particle local times -- the time a tracer particle spends at some predefined location along a single trajectory up to a time t. Exact results are derived for one- and two-tag local times, which reveal how the individual particles' histories become correlated at higher densities because each consecutive displacement along a trajectory requires collective rearrangements. Our results unveil the intricate meaning of projection-induced memory on a trajectory level, invisible to ensemble-average observables, and allow for a detailed analysis of single-file experiments probing tagged particle exploration statistics.

Single-file dynamics refers to the motion of particles in a narrow, effectively one-dimensional channel, which prevents their crossing, and is central to the transport in biological channels [1] the kinetics of transcription regulation [2], transport in zeolites [3] and in superionic conductors [4]. Recent advances in single-particle tracking and nanofluidics enabled experimental studies of single file dynamics in colloidal systems, which directly probe the fundamental physical principles of tagged particle motion to an unprecedented precision [5].
The motion of particles in a single file is strongly correlated, which gives rise to a rich and intricate phenomenology. In a Brownian single file the non-crossing constraint leads to subdiffusion with the ensemble mean squared displacement (MSD) of a tagged particle scaling as [x(t) − x(0)] 2 ∝ √ t [6]. When confined to a finite interval the subdiffusive scaling of the MSD is transient, saturating at an equilibrium variance, with the extent of the subdiffusive regime growing with the particle density (see Fig. 1a and [7]). Concurrently, an effective harmonization emerges at increasing density, with the invariant measure of a tagged particle approaching a Gaussian and a vanishing kurtosis excess γ 2 = x 4 eq / x 2 2 eq − 3 (see inset of Fig. 1a). More generally it holds that the asymptotic MSD of a tagged particle in a single file, and the absolute dispersion of a free particle are related via [x(t) − x(0)] 2 ∝ |x(t)| free [8]. The motion of particles on a many-body level is Markovian, the resulting tagged particle dynamics is, however, highly non-Markovian [7], and displays a staggering dependence on the respective initial conditions [9].
Tremendous effort has been made to study the tagged particle dynamics theoretically [10]. In particular, the tagged particle propagator has been studied using the 'reflection principle' [11], the Jepsen mapping [12], momentum Bethe ansatz [7], harmonization techniques [13], and macroscopic fluctuation [14] and large deviation [15] theory. Notwithstanding, these works, with isolated exceptions [16], focused on ensemble-average properties alone. State-of-the-art experiments, however, probe particle tra- Trajectories of two next-nearest neighbor particles in a single file of 11 particles (red and blue curves) alongside the respective left and right nearest neighbors (gray curves). Overlaid are corresponding local time fractions up to a time t, θ i t in the respective red and blue shaded intervals. The remaining particle trajectories are omitted for convenience.
jectories. Probing functionals of tagged particle trajectories is thus not only feasible but also more natural than studying ensemble-average observables. Moreover, to arrive at a physical understanding of projection-induced memory effects and resulting non-Markovianity, an understanding of the correlations of particle histories and their decorrelation on ergodic time-scales is required.
In particular, we here focus on the trajectory-, or timeaverage analogue of the tagged particle ensemble propagator [7], which in this context takes the form of the tagged particle local time fraction where 1 j y [x(τ )] = 1 if x j ∈ dy centered at y, and zero otherwise [6]. Eq. (1) represents a random quantity denoting the fraction of time the tagged particle j spends in an infinitesimal region around the point y along a trajectory up until time t (see Fig.1b), and where x(t) ≡ (x 1 (t), . . . , x N (t)) T is the many-body trajectory written in vector form. The dynamics of a tagged particle x i (t) irrespective of the other N − 1 is not Markovian, and any two tagged particle trajectories x i (t) and x j (t) are correlated on all but ergodically long times.
A general theory of local times in such correlated non-Markovian dynamics so far remained elusive. And while the statistics of functionals of the form in Eq. (1) in onedimensional stochastic processes have been studied extensively in a variety of fields [18], studies of tagged particle functionals in interacting many-body systems are sparse, and mostly limited to extreme value statistics of vicious walkers (see e.g. [19]).
Here, we present rigorous results for variances and two-tag correlations of bounded linear functionals of projected dynamics on arbitrary time-scales, in terms of the relaxation eigenspectrum of the corresponding Markovian embedding. The theory applies to all ergodic Markovian systems with a diagonalizable embedding propagator. As an example we study tagged particle local times in a single file of Brownian point particles in a box. Diagonalizing the many-body propagator using the coordinate Bethe ansatz, our results uncover striking trajectory-totrajectory fluctuations, and a cross-over from negatively to positively correlated two-tag particle histories upon increasing density, mirroring the emergence of collective fluctuations breaking Markovianity in tagged particle motion. The connection to the relaxation spectrum provides an intuitive understanding of non-Gaussian statistics at sub-ergodic times in a general setting.
Theory.-We consider a trajectory of a general Ndimensional system x(t) evolving according to Fokker-Planck or discrete-state Markovian dynamics. We are interested in ergodic systems with a unique steady-state P (x) and also assume steady state initial conditions, for which θ j t (y) = dx N δ(y − x j )P (x), where we introduced the Dirac delta function δ(x). In the presence of detailed balance (DB) P (x) is the Boltzmann-Gibbs measure P eq (x). We focus on the fluctuations and the correlations of local times where · · · denotes the average over all N -particle trajectories starting from the steady-state (in the case of DB equilibrium) and propagating up to time t. Obtaining Eqs. (2)(3) essentially amounts to computing the probability generating function of the joint local time functional given by the Feynman-Kac path integral A straightforward generalization of the trotterization in Ref. [3] shows that Q u,v (x i , y j |t) is the propagator of a tilted evolution operator [21] whereL andL † denote the 'bare' forward and adjoint (backward) operator [22], and we introduced the 'flat' |-≡ dx|x and steady states |ss = dxP (x)|x in the bra-ket notation, which are the left (right) and right (left) ground eigenstates ofL (L † ), respectively. We obtain exact expressions for the moments in Eqs. (2)(3) by performing a Dyson series-expansion of Eq. (S3) [8], converging for any bounded functional of x(t) (see [21]). Having assumed diagonalizability ofL † (andL), we expand the backward operator in a complete biorthogonal set of left and right eigenstates [24],L † = k λ k |ψ L k ψ R k |, λ k denoting the (possibly degenerate) eigenvalues and ψ L k |ψ R l = δ kl . As the calculation of the moments is tedious we here simply state the result. The details are shown in [21]. Obviously, ψ R 0 |1 i x |ψ L 0 = P (x), since the system is ergodic. The exact results for the variance and correlations are remarkably simple and read where we introduced the auxiliary function The exact large deviation (LD) limits of Eqs. (6-7) readily follow in the limit t λ −1 where denotes asymptotic equality. Analogous formulas for LD limits of local times not connected to a spectral expansion have also been developed (see e.g. [26]). Notably, for systems obeying DB σ 2,LD xi (t) sets a universal upper bound on the variance of θ t (compare Eqs. (6) and (8)). The results in Eqs. (6)(7)(8)(9) readily extend to arbitrary functionals t −1 t 0V [x(τ )]dτ with a bounded and local V , by performing a simple exchange 1 i x →V , modifying only Ω k (x i , y j ) [21]. Our theory applies to all diagonal-izableL, thus including all systems obeying DB.
Eqs. (6)(7)(8)(9) provide an intuitive understanding of local time statistics via a mapping onto relaxation eigenmodes, with fluctuation and correlation amplitudes proportional to the sum of transition amplitudes of excitations from the steady state to excited states and back, Ω k (x i , y j ). On ergodic time scales θ t at different t decorrelate, and hence display features of shot-noise, i.e. σ 2 xi (t) and C ij xy (t) decay inversely proportional to the number of independent observations of each excitation mode, ∼ λ −1 k /t. At finite times t λ −1 k shot-noise statistics are altered due to a finite survival probability of the eigenmodes at given t (see correction terms in brackets of Eqs. (6)(7).) Local times in single-file diffusion.-Consider the dynamics of N identical hard-core interacting Brownian point particles diffusing in the unit interval [0, 1], and set D = 1 without loss of generality. The extension to a finite particle radius follows from a trivial change of coordinates [7] The non-crossing condition on the N individual particle trajectories is imposed via the set of N − 1 reflecting internal boundary conditions lim xi+1→xi (∂ x0,i+1 − ∂ x0,i )P (x 0 , t|x) = 0 ∀i. Confinement into a unit interval is imposed through external reflecting boundary conditions ∂ x1 P (x 0 , t|x)| x0,1=0 = ∂ x N P (x, t|x 0 )| x 0,N =1 = 0. Under these boundary conditions we diagonalizeL † using the coordinate Bethe ansatz [27,28] and obtain the Bethe eigenvalues λ k = π 2 i k 2 i and corresponding left and right eigenvectors where m k is the multiplicity of the Bethe eigenmode |ψ L k [21], and {ki} denotes the sum over all permutations of single-particle eigenvalues with k i ∈ N 0 .
The matrix elements entering Ω k (x i , y j ) follow upon integration over the n l and n r particle coordinates to the left and right, respectively, from the tagged particle i while strictly preserving the particle ordering [7], yielding and This delivers exact results for σ 2 xi (t) and C ij xy (t) in Eqs. (6)(7)(8). An efficient numerical implementation of our analytical results can be made available upon request.
The results for σ 2 xi (t) in Eq. (6) for the central particle in single files with various N are depicted in Fig. 2, and reflect large fluctuations exceeding 200% on time-scales where roughly only 50% of the particles have collided with their neighbors. The fluctuations display a nontrivial dependence on x, which does not follow the shape of P eq (x t ) = N !x n l t (1 − x t ) nr /(n l !n r !), and reveal striking boundary-layer effects. At longer t, where ∼ 50-100 collisions/particle have occured, θ i t at different t become uncorrelated according to the central limit theorem, with σ 2 xc (t) converging to its LD limit (8). On these timescales the ensemble MSD has already saturated (compare Figs. 1a and 2c and d). Notably, LD asymptotics correctly capture only small fluctuations of the order ±10%. As noted above and confirmed by computer simulations, large deviations reflecting Gaussian statistics set an upper bound to the fluctuations of θ i t ( Fig. 2c and d). Single-file diffusion displays no time-scale separation in the relaxation spectrum. As a result, the projection of dynamics onto a tagged particle coordinate induces subdiffusion and strong non-Markovianity on time scales t < λ −1 . The respective onset of the √ t scaling of the tagged particle MSD shifts to shorter t upon increasing N (Fig. 1a). Increasing N in turn leads to a high degeneracy of Bethe eigenmodes, reflecting emerging dynamical symmetries [21]. As a result, fewer Bethe modes are required for a convergence of the sums in Eqs. (6)(7).
To gain deeper insight into the physical origin of the memory on a trajectory level we analyzed twotag correlations between particle histories by means of the reduced covariance of local timesC ij Correlations between the histories of the central particle c and its nearest c + 1 and next-nearest c + 2 neighbors at the midpoint between the maxima of P eq (x c ) and P eq (x c+1,c+2 ) [21] are depicted in Fig. 3. Due to ergodicity, θ i t (x) become very weakly correlated at long t and Gaussian statistics emerge. Consequently, C ij xy (t) vanishes for long times, after 10 2 collisions took place on average. Note that C ij xy (t) measures correlations between particle histories and not particle positions. The latter never decorrelate, i.e. two-tag position correlation functions display an algebraic decay even at equilibrium P ij where n l,r and m l,r are the number of particles to the left/right of the two tagged particles i and j (for details see [21]).
Notably, we observe a transition from negatively to Mean local time (blue lines) and standard deviation (shaded area) for a) the first and third, and b) first and 8-th tagged particle in a single file with N = 3 and N = 10, respectively at three different lengths of trajectories. c) and d) Reduced variance of local time of the central particle σ 2 xc=1/2 (t)/ θ c t (xc = 1/2) 2 for various N . The full lines denote exact results from Eq. (6) and dashed lines large deviation asymptotics Eq. (8). Symbols correspond to Brownian dynamics simulation of an ensemble of 10 6 independent trajectories starting from equilibrium initial conditions.
y) ), reduced two-tag local time correlation functions of the central particle c and its nearest (a) and next-nearest (b) neighbor for different N . Time is expressed in units of the mean collision time. Lines depict the theory in Eq. (7) whereas symbols correspond to Brownian dynamics simulations of 10 6 independent trajectories starting from equilibrium initial conditions. positively correlated tagged particle histories upon increasing density (Fig. 3), mirroring a change in particle dynamics from single-particle to collective fluctuations. The driving force for this transition can be found in an enhanced packing at higher densities resembling a 'crystallization' transition, where invariant tagged particle densities P eq (x) become strongly overlapping, whereas their respective widths shrink only very slowly (see [21]).
The 'critical' density, at which the behavior shifts from negatively to positively correlated histories, depends on the topological separation between the two tagged particles and is shifted to higher values of N for more distant particles (compare a) and b) in Fig. 3). In turn, this reflects a growing dynamical correlation-length with increasing N . Moreover, upon increasing N ,C ii xy (t) of the central particle becomes non-monotonic, with weak anti-correlations at short t turning to weak correlations at large t, before reaching the LD limit of uncorrelated histories, where harmonization [13] ideas apply. The increasingly positive correlations with growing N reflect a persistence and a finite life-time of typical collective fluctuations on a trajectory level, akin to glassy dynamics in kinetically constrained models [29]. Accordingly, positive correlations are are not observed if we tag outer particles at external boundaries (see [21]).
Conclusions.-We established a general method for determining exactly the variance and two-tag correlations of bounded non-negative functionals of stationary ergodic Markov processes with a diagonalizable propagator. The theory relates the statistics of functionals to the relaxation eigenspectrum, and allows for an exact treatment of non-Markovian dynamics from a Markovian embedding. It also holds for diagonalizable irreversible dynamics, where a broken time-reversal symmetry can cause oscillations in higher order terms in Eqs. (6)(7)(8)(9) and/or fluctuations exceeding the large deviation limit in Eq. (8). From the spectrum of the many-body prop-agator obtained via the coordinate Bethe ansatz, we derived exact results for one-and two-tag local times in single file diffusion, which unveiled non-trivial correlations between tagged particle histories and the emergence of collective dynamics at increasing particle densities. Going beyond large deviation time-scales, our results revealed that harmonization concepts, assuming dynamics in-between local equilibria -an assumption that works well for ensemble-average observables [13] -fail on the more fundamental trajectory level. This highlights the intricate physical meaning of projection-induced memory on the level of single trajectories, which is virtually invisible to ensemble-average observables. Our results on local times can readily be tested experimentally (see e.g. [5]), and hopefully our theory will stimulate further research directed towards tagged particle functionals. Particularly interesting would be extensions to tagged particle dynamics in rugged potential landscapes [30].
We thank David Hartich for insightful discussions and critical reading of the manuscript. The financial support from the German Research Foundation (DFG) through the Emmy Noether Program "GO 2762/1-1" (to AG), and an IMPRS fellowship of the Max Planck Society (to AL) are gratefully acknowledged.  (7) Let x(t) be an arbitrary-dimensional ergodic Markov process on a discrete or continuous state-space. The evolution of the probability density function evolves under the corresponding diagonalizable forward operatorL (e.g. Fokker-Planck-or discrete-state master equation-type) with invariant measure P (x) and the adjoint (i.e. backward) operator L † . Let the respective eigenspectra beL = k λ k |ψ R k ψ L k |, λ k andL † = k λ k |ψ L k ψ R k |, λ k denoting the possibly degenerate and in general complex-valued eigenvalues. Note thatL|ψ R k = λ k |ψ R k andL † |ψ L k = λ k |ψ L k , i.e. the left and right eigenstates span a bi-orthogonal eigenspace ψ L k |ψ R l = δ kl [1]. The forward and backward propagators of the process can then be written as [1] The probability density function of a bounded linear functional ϕ t = t 0V [x(τ )]dτ over all paths starting from a (potentially non-equlibrium) steady-state and propagating up to time t, is defined by the path integral with the corresponding stochastic action functional S[x(t)] of the continuous [2,3] or discrete state-space [5] Markov process x(t), and where we introduced the Dirac delta function δ(x). By means of a straightforward vectorial generalization of the trotterization of the the path integral (S2) in Refs. [3,4] (for the backward and forward approach, respectively), one finds that the generating function, corresponding to the Laplace transformF(u|t) = ∞ 0 dϕe −uϕ F(ϕ|t), is the propagator of a tilted operator where we have introduced the 'flat' |-≡ dx|x and steady states |ss = dxP (x)|x , which are the left (right) and right (left) ground eigenstates ofL (L † ), respectively. The last equality follows fromF(u † |t) † =F(u|t). In taking the Laplace transform we assumed that the functional has non-negative support (such as in the case of local times). In case the support extends to negative values one simply needs to take the Fourier transform instead. The moments of F(ϕ|t) at any given t follow immediately from ϕ n t = (−1) n ∂ n uF (u|t)| u=0 , where · · · denotes the average over all trajectories starting from a steady state and propagating up to time t. In case the Fourier transform is used, a corresponding change of the prefactor is required.
For bounded linear functionals of ergodic Markov processes all moments are finite, | ϕ n t | ≤ f (t) lim t→∞ | ϕ n t | < ∞ with a smooth scaling function f (t), which depends on the detailed form ofV [x(t)]. Moreover, F(ϕ|t) obeys a large deviation principle [6,7]. In the specific case of local times, treated in the Letter,V [x(t)] = 1 j y [x(t)] and f (t) ∝ t n for ϕ n t . The finiteness of moments implies thatF(u|t) is an analytic (i.e. holomorphic) function of u at least at and near u = 0 for any t.
To obtain exact results for second moments we simply need to expandF(u|t) in a Dyson series to second order in uV preserving the time-ordering, and afterwards take the second derivative at u = 0. The series is guaranteed to converge, sinceV is bounded. Because trivially -|e −tL |ss = ss|e −tL † |-= 1, the Dyson expansion gives [8] with t ≥ t ≥ t ≥ 0. An equivalent expansion can be obtained forL. The Dyson series (S4) converges for u ∈ C < ∞.
To prove the convergence for any bounded linear operatorB we consider the operator norm. Let Ψ be a complete normed linear space, andB : Ψ → Ψ. The operator norm is then defined as ||B|| = sup ||ψ||=1 ||Bψ|| with ψ ∈ Ψ. The operator norm corresponds to the largest valueB stretches an element of Ψ. SinceB is bounded we have ||B N || ≤ ||B|| N , ∀N ∈ N, which follows simply from ||ÂB|| ≤ ||Â|| ||B||. The operator exponential is defined as the limit eB = lim N →∞ N k=0B k /k! and the convergence is in operator norm, since || N k=0B The series on the right hand side converges absolutely for any number ||B|| ∈ C. Due to the completeness of the space Ψ, eB as well belongs to a complete normed linear space, and moreover ||eB|| ≤ e ||B|| . TakingB = uV with u ∈ C completes the proof of convergence of the series (S4).
Utilizing now the spectral expansionL † = k λ k |ψ L k ψ R k | in Eq. (S4) we obtain for the first order term where we used the fact that ss| and |are the left and right ground states ofL † as well as the bi-orthogonality of the eigenbasis. The second order term follows similarly We can now trivially extend uV → uÂ+vB for u, v ∈ C and any two bounded linear operatorsÂ andB. In the specific case of tagged particle local times studied in the Letter we haveÂ dy centered at y, and zero otherwise [6]. The exact second moments are now obtained from ∂ u ∂ vF (u|t)| u=v=0 and ∂ 2 uF (u|t)| u=0 by considering the corresponding operatorsÂ andB. Finally, since in the Letter we consider the local time fraction and not the total local time, we must take t −2 ∂ u ∂ vF (u|t)| u=v=0 and t −2 ∂ 2 uF (u|t)| u=0 , respectively. This completes the proof of the main general results of the Letter, i.e. Eqs. (6) and (7).

EXTENDED PHASE-SPACE INTEGRATION IN SINGLE-FILE DIFFUSION
The integrals involved in the evaluation of invariant measures and matrix elements in single-file diffusion involve nesting, i.e. the ordering of particles is strictly preserved This imposes non-trivial topology of the phase space of the system. A tremendous simplification is achieved through the so-called 'Extended Phase-Space Integration' developed by Lizana and Ambjörnsson, which exactly reduces the nested high-dimensional integrals to scaled single particle integrals, e.g. [9]: where n l and n r are the number of particles (integrals) to the left and right of the tagged particle m, respectively. The extended phase-space integration in Eq. (S8) applies to all functions f (x), which are invariant under the exchange x i ↔ x i+1 . Throughout our work all nested integrals included in the bra-s ψ k | (scalar products, matrix elements etc.) are evaluated using the extended phase-space integration.

EIGENMODE MULTIPLICITY AND EIGENVALUE DEGENERACY
As described in the Letter we diagonalize the many-body Fokker-Planck operator using the coordinate Bethe ansatz method. Each Bethe eigenstate of a Single-File of N particles is uniquely defined by a tuple k = (k 1 , k 2 , . . . , k N ). To each tuple corresponds one eigenvalue through the relation: since more than one tuple may correspond to the same eigenvalue, these are degenerate. To each tuple k it is possible to associate a set K containing the elements of k counted once. Defining n Ki as the number of times the element K i appears in the tuple k, we define the multiplicity of the eigenvectors associated to k as TAGGED PARTICLE EQUILIBRIUM PROBABILITY DENSITIES The exact tagged particle equilibrium probability density function of the tagged particle i is obtained by a nested integration of all other particle positions where n l and n r are, respectively, the number of particles to the left and to the right of the tagged particle i. Fig. 4 depicts results for P eq i (x) for the central particle c and the two nearest neighbors to the right, c + 1 and c + 2, respectively. The probability density of the central particle approaches a Gaussian shape as the number density N increases. For large enough N the width of P eq i (x) stops decreasing appreciably, while the probability densities of neighboring particles begin to overlap strongly. This has important physical consequences for correlations of particle histories, as we explained in the discussion of Fig. 3 in the Letter.  In order to allow for a meaningful comparison of results for different particle numbers N we need to choose appropriate reference conditions. To do so, we focus only on odd particle numbers, for which the system is symmetric with respect to the peak of the invariant measure of the central particle P eq (x c ). This way a comparison of correlations with nearest c + 1 and next-nearest c + 2 neighbors at different densities is indeed consistent. Moreover, in order to compare equilibrium and near-equilibrium tagged particle excursions with far-from equilibrium fluctuations we choose the following reference points with respect to P eq (x i ): The point x 50 , in which x m N = 31 P eq x 16th particle 17th particle FIG. 5. Invariant measures for the central particle and its nearest neighbor, denoting the different kinds of reference points.
particle histories for two particles i and j we focus on the mid-point x m (i, j) = (x 50,i + x 50,j )/2.

CONVERGENCE RATES OF SERIES AND EIGENVALUE DEGENERACY
The exact expressions for variance and covariance of local time of a tagged particles in Eqs. (6) and (7) in the Letter involve an infinite series, whose rate of convergence is difficult to predict, as it strongly depends on the particular position x of the tagged particle under inspection, as well as on the number of particles N and D(k), the degeneracy of Bethe eigenvalue λ k . To inspect the rate of convergence of the series we compute the relative deviation of the results for the variance of local time of the central particle truncated at the kth Bethe eigenvalue, |σ 2 xc (t) − σ 2 k (t)|/σ 2 (t) as a function of k at different positions x and at different lengths of trajectories t. Fig. 6 depicts how fast the series for the variance of local time of the central particle (Eq. (6) in the Letter) truncated at the kth term converges to the exact value k → ∞. n order to compare systems with different N we focused on points x 50 and x 75 of the central particle, with the specific values given in Tab. I).
Intuitively, the convergence rate increases with increasing length of the observation t, since faster modes must become less and less important. The convergence rate also increases with increasing N , which is due to an increasing degeneracy of lower-lying eigenvalues at larger N . Degenerate low-lying eigenvalues allow for a mixing of different collective slow modes, which become dominant. Finally, by comparing the columns of Fig. 6 we notice that the rate of convergence also depends on the tagging position, which in turn depends on the curvature of the modes at different N . Analytical results for |σ 2 xc (t) − σ 2 k (t)|/σ 2 (t) of the central particle for different particle numbers (N ) as a function of k at different tagging positions. The black symbols depict the eigenvalue degeneracy D(k).

EQUILIBRIUM POSITION CORRELATION FUNCTION
In the Letter we focused on the covariance of tagged particle local timesC ij reflecting the correlations between particle histories. We found that histories decorrelate at long times as a consequence of the central limit theorem. Conversely, the particle positions on the ensemble average level do not decorrelate, not even in equilibrium. To demonstrate this we compute exactly the pair correlation function where n l,r and m l,r are the number of particles to the left/right of the two tagged particles i and j. Here we want to focus on the ensemble-average reduced pair correlation function P ij eq (x i , x j )/(P eq (x i )P eq (x j )) depicted in Fig. 7. The latter is in general different from 0, while the former goes to 0 in the limit t → ∞. P c,c+1 eq (x, y) depends on N and for large N changes monotonically from positive to negative correlations as a function of particle separation. At small N the correlations becomes weaker with increasing separation, but remains positive. Intuitively, as one particle lies in-between, the dependence of P c,c+2 eq on the interparticle separation is for large N non-monotonic, going from perfectly anticorrelated to correlated and back to anti-correlation. At low N P c,c+2 eq depends non-trivially on the interparticle separation, such that the aforementioned terminal anticorrelation disappears for small enough N . Notably, for x c = x 90 the correlations are much stronger, which suggest that more extensive excursions are entropically penalized as they demand collective fluctuations.

TWO-TAG CORRELATION FUNCTION OF LOCAL TIMES
In the Letter (in particular in Fig. 3) we analyzed one-point two-particle historiesC c,c+1,2 xm,xm (t) at the respective mid-point positions x m listed in Table III. . P ij eq (xi, xj) between the central particle c and the right nearest and next-nearest neighbor, c + 1 and c + 2 respectively, when the central particle is tagged at xc = x50, x75, x90 and when xc = xm (see Tab. II for the exact positions). The vertical dashed lines are drawn at the position of the central particle (see Tab. II) to denote the non-crossing boundary condition.
To gain further insight we also compute the first and central-particle two-point reduced correlation functions of local timesC 1,1 xy (t) andC c,c xy (t) with x and y given in Table IV. Fig. 8). In the left plot we show the self-reduced reveals weak anti-correlations at short times t turning to weak correlations at longer t, before reaching the large deviation limit of uncorrelated histories. As already mentioned in the Letter, the correlations for the outer particles are weak and become weaker with increasing number of particles N , since the outer particles are constrained between the reflecting wall and the right nearest neighbor. In the case of the central particle we observe further evidence of the emergence of persistent collective fluctuations at higher densities, observed and described in the Letter. Notably here even near-equilibrium fluctuations reveal signatures of collective behavior in the form of persistent histories (note that we are tagging at x 50 and x 75 ), i.e.C c,c x50,x75 (t) turns from purely negative to weakly positive correlations.  x 50 y 75 (t) of the first particle and b)C c,c x 50 y 75 (t) of the central particle.