# Tracer dynamics in the active random average process

Saikat Santra1, Prashant Singh2 and Anupam Kundu1 International Centre for
Theoretical Sciences, Tata Institute of Fundamental Research, Bengaluru
560089, India 1 Niels Bohr International Academy, Niels Bohr Institute,
University of Copenhagen, Blegdamsvej 17, 2100 Copenhagen, Denmark 2
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We investigate the dynamics of tracer particles in the random average process
(RAP), a single-file system in one dimension. In addition to the position,
every particle possesses an internal spin variable $\sigma(t)$ that can
alternate between two values, $\pm 1$, at a constant rate $\gamma$.
Physically, the value of $\sigma(t)$ dictates the direction of motion of the
corresponding particle and for finite $\gamma$, every particle performs a non-
Markovian active dynamics. Herein, we study the effect of this non-
Markovianity in the fluctuations and correlations of the positions of tracer
particles. We analytically show that the variance of the position of a tagged
particle grows sub-diffusively as $\sim\zeta_{\text{q}}\sqrt{t}$ at large
times for the quenched uniform initial condition. While this sub-diffusive
growth is identical to that of the Markovian/non-persistent RAP, the
coefficient $\zeta_{\text{q}}$ is rather different and bears the signature of
the persistent motion of active particles through higher point correlations
(unlike in the Markovian case). Similarly, for the annealed (steady state)
initial condition, we find that the variance scales as
$\sim\zeta_{\text{a}}\sqrt{t}$ at large times with coefficient
$\zeta_{\text{a}}$ once again different from the non-persistent case. Although
$\zeta_{\text{q}}$ and $\zeta_{\text{a}}$ both individually depart from their
Markovian counterparts, their ratio $\zeta_{\text{a}}/\zeta_{\text{q}}$ is
still equal to $\sqrt{2}$, a condition observed for other diffusive single-
file systems. This condition turns out to be true even in the strongly active
regimes as corroborated by extensive simulations and calculations. Finally, we
study the correlation between the positions of two tagged particles in both
quenched uniform and annealed initial conditions. We verify all our analytic
results by extensive numerical simulations.

##  1 Introduction

The dynamics of a tracer particle in a collection of non-overtaking particles
in one dimension is a prototypical example of strongly correlated system in
statistical physics. This non-overtaking constraint, called single-file
constraint, drastically changes the dynamical behaviour of a tracer particle
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]. For example, in single-file diffusion,
the mean squared displacement (MSD) of a tagged particle grows sub-diffusively
as $\sim\sqrt{t}$ at late times in contrast to the linear growth of the MSD of
a free diffusive particle. The coefficient of the sub-diffusive growth depends
on the particle number density and bare diffusion of the particles[2, 5, 4].
For Hamiltonian systems, the single-file constraint, slows down the motion of
a tagged particle [11, 12, 13, 14]. This slowing down of the dynamics is a
common effect in the single-file motion and occurs due to the hindrance in the
motion faced by one particle due to the presence of other particles. It has
been found that the coefficient of the late time growth of the MSD of a tagged
particle crucially depends, in addition to the density of the particles, on
the microscopic dynamics of individual particles, interactions among
themselves and on the statistical properties of the initial state [13, 12, 15,
16, 9, 10]. In this paper, we study how the motion of a tagged particle gets
modified if all particles in the single-file system are active.

Active matter is a class of driven out-of-equilibrium systems where every
individual unit consumes energy from the environment and converts it into a
systematic movement via some internal mechanisms [17, 18]. At the collective
level, these particles exhibit interesting phenomena such as motility-induced
phase separation, absence of equation of state for pressure etc [19, 22, 20,
21]. Non-interacting active particles also show behaviours which are different
than their passive counterpart as exemplified by clustering inside bounded
domain, climbing against potential hill, non-Boltzmann stationary state and
large deviation in position distributions and survival probability which are
different than the thermal particles [23, 24, 28, 29, 25, 26, 27, 30, 31, 32,
33, 34, 35, 36, 37, 38, 40, 41, 39]. At the interacting front, the
distribution of two active particles with mutual exclusion has been studied
and shown to display jamming features [42]. Going beyond two particles, there
have also been attempts to derive fluctuating hydrodynamic descriptions for
active lattice gases that turn out to be useful to study density and current
fluctuations and entropy production [43, 44, 45]. With the growing interest in
active matter in the last decade, people naturally got interested in knowing
how broken detailed balance at the microscopic dynamics manifests itself in
the tracer dynamics in single-file motion. Several numerical studies in this
direction have pointed out that the temporal growth of the MSD of a tagged
particle, at late times, remains same as in the absence of activity. However,
the coefficient associated with this temporal growth gets non-trivially
modified due to the persistent nature of these particles [49, 50, 51, 48, 46,
47, 52]. Attempts to derive this with harmonic chain of active particles
reproduce only the passive result at late times and do not shed light on the
role of activity in the tracer dynamics [52, 53]. It is, however, not
difficult to realize that the presence of activity increases the correlations
among particles and this, in addition to the single-file constraint, should
affect the motion of a tagged particle. Naturally one may ask how such
enhanced correlations affect the motion of a tagged particle? While the above
mentioned studies discuss the overall effect of the presence of activity, the
contribution from the enhanced correlation is not very clear and transparent.
Moreover, how the fluctuations in the initial conditions affect the motion of
tracer particles in active single-file systems is also not explored. In
absence of a general formulation to investigate these questions, it is crucial
to study specific model systems that are amenable to analytical calculations.
In this paper, we consider a version of the random average process (RAP) in
which individual particles are subjected to active noises. For this model, we
provide systematic answers to these questions.

Originally, the RAP model was first studied for non-persistent particles
(without active noises) by Fontes and Ferrari as a generalization of the
smoothing process and the voter model [54]. It has also appeared in several
other physical problems like force fluctuations in bead packs [55], in mass
transport models [56, 57], in models of wealth distribution and traffic [59,
58] and generalized Hammersley process [60]. In this model, motion of the
particles is restricted as each particle can jump, with a fixed rate $p$, on
either side only by a random fraction $\eta$ of the space available till the
next neighbouring particle [54, 61]. The random fraction $\eta$ is chosen from
some distribution $R(\eta)$. Since particles cannot overtake their
neighbouring particles, their initial order remains preserved throughout the
time evolution which gives rise to the single-file motion. Meanwhile, the jump
that a particle makes at a given instant is independent of what it does in the
previous step. Therefore, we will refer to this model as the Markovian RAP
(MRAP). Later, we will contrast this with the active case where every particle
possesses a spin variable $\sigma(t)$ that dictates the direction of its
motion and has non-vanishing correlations at two different times.

Being a paradigmatic model for interacting multi-particle system, the motion
of tagged particles in MRAP was studied previously and several results were
obtained both analytically and numerically. If $x_{i}(t)$ denotes the position
of the $i$-th particle at time $t$, then the MSD $\langle z_{i}^{2}(t)\rangle$
(with $z_{i}(t)=x_{i}(t)-x_{i}(0)$) and the correlation between two tagged
particles $\langle z_{i}(t)z_{j}(t)\rangle$ were computed using microscopic
calculations [61, 57, 62, 63], as well as using hydrodynamic approach [64]. At
late times, these quantities are explicitly given by

$\left.\begin{aligned} &\langle
z_{0}^{2}(t)\rangle\simeq\frac{\rho^{-2}\mu_{2}\sqrt{\mu_{1}p}}{(\mu_{1}-\mu_{2})\sqrt{\pi}}~{}\sqrt{t},\\\
&\langle
z_{0}(t)z_{i}(t)\rangle\simeq\frac{\rho^{-2}\mu_{2}\sqrt{\mu_{1}p}}{(\mu_{1}-\mu_{2})\sqrt{\pi}}~{}\sqrt{t}~{}f\left(\frac{|i|}{\sqrt{4\mu_{1}pt}}\right),\end{aligned}\right\\}\text{
(MRAP)}$ (1)

where $\mu_{i}$ is the $i^{\rm th}$ moment of the jump distribution $R(\eta)$
and $\rho$ is the stationary density of the particles. Explicit expression of
the scaling function $f(y)$ is given by [61, 64]

$\displaystyle f(y)=e^{-y^{2}}-\sqrt{\pi}y~{}\text{erfc}(y).$ (2)

The same scaling function also appears in many single-file systems that
possess diffusive hydrodynamics at the macroscopic scales [64, 65, 3, 66, 67].
It is, however, important to note that for MRAP, the equations for two-point
correlation functions close onto themselves. Therefore, two-point correlations
are enough to decide the pre-factor in the expressions of MSD. Contrarily,
this is found not to be true when the particles are subjected to active
noises. Then, the two-point correlations would depend on three-point
correlations and so on. Consequently, higher point correlations start to
contribute to the growth of the MSD through two-point correlations. For
example, such dependence on higher order correlations was also observed for
gap statistics in hardcore run and tumble particles [68]. Question then
arises: does activity facilitate in the growth of MSD? If so, how? In this
paper, we present an example of a model where these questions can be
thoroughly addressed through microscopic analytic computations aided by
numerical simulations.

Our paper is organised as follows: In Section 2, we introduce the model, fix
notations and also summarize main results of the paper. We then compute the
spin-position correlation function $C_{i}(t)=\langle
z_{i}(t)\sigma_{0}(t)\rangle$ in Section 3. We use this result to study the
properties of tracer particles with fixed initial condition in Section 4. More
specifically, we look at the MSD of the position of a tagged particle in
Section 4.1 and position correlation of two tagged particles in Section 4.2.
Section 5 discusses these quantities in the annealed case with Section 5.1
devoted to the MSD and Section 5.2 to the correlation function. Since most of
our analytic results are derived with mean-field approximation, we discuss the
validity of our results for strongly active regime (small $\gamma$) in Section
6. Finally we conclude in Section 7.

Figure 1: Schematic illustration of the random average process with active
particles. Each particle has an internal spin denoted by the arrows inside the
circles and the direction of the arrow represents the state of the spin. The
spin variable of individual particles changes direction independently with
rate $\gamma$. In a small time interval $dt$ a particle chosen at random
either makes a jump in the direction of the spin with probability $pdt$ or
does not jump with the remaining probability $(1-pdt)$. Every successful jump
of a particle, say $i^{\rm th}$ with $\sigma_{i}=1$, takes place by a random
fraction $\eta_{i}$ of the space available in jump direction i.e. by an amount
$\eta_{i}\left(x_{i+1}(t)-x_{i}(t)\right)$ where $\eta_{i}\in[0,1)$ is a
random variable drawn from the distribution $R(\eta_{i})$.

##  2 Model and summary of our main results

We consider active particles moving in an infinite line distributed with
density $\rho$. We denote the position of the $i$-th particle at time $t$ by
$x_{i}(t)$ where $i\in\mathbb{Z}$ and $x_{i}(t)\in\mathbb{R}$. In addition,
every particle has an internal variable $\sigma_{i}(t)$ (called spin) which
can alternate between $\pm 1$ at a rate $\gamma>0$. The variable
$\sigma_{i}(t)$ represents the usual dichotomous noise widely studied for the
run and tumble particles [23]. Initially, the positions of these particles are
fixed and are kept at a fixed distance $a=1/\rho$ apart. However, the initial
value of $\sigma_{i}(0)$ can be random which, for simplicity, we choose to be
$\pm 1$ with equal probability $1/2$. Thus for all $i\in\mathbb{Z}$, we have

$\displaystyle x_{i}(0)$ $\displaystyle=ia=i/\rho,$ (3)
$\displaystyle\sigma_{i}(0)$
$\displaystyle=\begin{cases}~{}~{}1~{}~{}\text{with probability }1/2,\\\
-1~{}~{}\text{with probability }1/2.\end{cases}$ (4)

At a small time interval $[t,t+dt]$, the direction of motion of the $i$-th
particle depends on its spin variable $\sigma_{i}(t)$. If $\sigma_{i}(t)=1$,
then the particle jumps to the right with probability $pdt$ and does not jump
with probability $(1-pdt)$. On the other hand, if $\sigma_{i}(t)=-1$, then the
particle jumps to its left with probability $pdt$ and with probability
$(1-pdt)$, stays at $x_{i}(t)$. The jump, either to the left or to the right,
is by a random fraction $\eta_{i}$ of the space available between the particle
and its neighbour. This means, when the particle jumps to the right, it will
jump by an amount $\eta_{i}\left[x_{i+1}(t)-x_{i}(t)\right]$ whereas jump to
the left takes place by an amount $\eta_{i}\left[x_{i-1}(t)-x_{i}(t)\right]$.
The jump fraction $\eta_{i}\in[0,1)$ is a random variable drawn independently
from the distribution $R(\eta)$ and characterized by the moments
$\mu_{k}=\langle\eta^{k}\rangle$. In contrary to the original MRAP, we see
that the motion of a particle at time interval $[t,t+dt]$ depends on the value
of $\sigma(t)$ which itself depends on its previous history. Thus, every
particle performs a non-Markovian active dynamics. A schematic illustration of
this model is shown in Figure (1). We refer to this model as the active random
average process (ARAP).

The time evolution equation for the position $x_{i}(t)$ and spin
$\sigma_{i}(t)$ can be written as

$\displaystyle x_{i}(t+dt)=x_{i}(t)+\Gamma_{i}(t),$ (5)
$\displaystyle\sigma_{i}(t+dt)=\begin{cases}-\sigma_{i}(t),~{}~{}~{}\text{with
probability}~{}~{}~{}~{}~{}\gamma dt,\\\
~{}~{}\sigma_{i}(t),~{}~{}~{}\text{with probability}~{}~{}(1-\gamma
dt),\end{cases}$ (6)

where the increment $\Gamma_{i}(t)$ reads

$\displaystyle\Gamma_{i}(t)=\begin{cases}\eta_{i}\left[x_{i+1}(t)-x_{i}(t)\right],~{}~{}~{}\text{with
probability}~{}~{}~{}~{}\left(\frac{1+\sigma_{i}(t)}{2}\right)pdt,\\\
\eta_{i}\left[x_{i-1}(t)-x_{i}(t)\right],~{}~{}~{}\text{with
probability}~{}~{}~{}~{}\left(\frac{1-\sigma_{i}(t)}{2}\right)pdt,\\\
~{}~{}~{}~{}~{}~{}~{}0,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\text{with
probability}~{}~{}~{}~{}~{}~{}~{}(1-pdt).\end{cases}$ (7)

Since we are interested in studying the variance and correlation of the
displacement of the tagged particles, it seems convenient to work in terms of
the displacement variable $z_{i}(t)=x_{i}(t)-x_{i}(0)$. Also, the rate $p$,
essentially fixes a time scale in the model and can be scaled out by
redefining $\gamma\to\gamma/p$. Hence, from now on, we choose $p=1$ without
any loss of generality. The aforementioned update rules then become

$\displaystyle z_{i}(t+dt)=z_{i}(t)+\Gamma_{i}(t),~{}~{}~{}\text{with }$ (8)
$\displaystyle\Gamma_{i}(t)=\begin{cases}\eta_{i}\left[z_{i+1}(t)-z_{i}(t)+a\right],~{}~{}~{}\text{with
probability}~{}~{}~{}~{}\left(\frac{1+\sigma_{i}(t)}{2}\right)dt,\\\
\eta_{i}\left[z_{i-1}(t)-z_{i}(t)-a\right],~{}~{}~{}\text{with
probability}~{}~{}~{}~{}\left(\frac{1-\sigma_{i}(t)}{2}\right)dt,\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}0,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\text{with
probability}~{}~{}~{}~{}~{}~{}~{}(1-dt).\end{cases}$ (9)

For the Markov case $(\gamma\to\infty)$, we saw in Eq. (1) that the variance
and the equal time correlation function exhibit interesting scaling behaviours
with time. Here, we derive them for finite $\gamma$ illustrating the interplay
of single-file constraint and non-Markovianity in the dynamics. Before,
delving into the main calculation, we present a brief summary of all our
results in the paper.

###  2.1 Summary of the main results

  * •

We first study the two-point correlation $C_{i}(t)\equiv\langle
z_{0}(t)\sigma_{i}(t)\rangle$ between the displacement of one particle and the
spin of another particle. This basic quantity turns out to be instrumental in
determining the two-point position correlations for the active RAP. However,
it is absent in systems without persistent noise $\sigma_{i}(t)$ i.e. in the
Markovian RAP system [61]. In Section 3, we obtain an explicit expression of
$C_{i}(t)$ based on mean-field (MF) approximation (see Eq. (25)). We observe
that this correlation, starting from zero, increases with increasing $t$ (see
Figure (2)) and saturates to a constant value as $t\to\infty$. This value is
given by

$\displaystyle
C_{i}(t\to\infty)=\frac{\mu_{1}a}{2\sqrt{\gamma^{2}+\mu_{1}\gamma}}~{}\text{exp}{[-|i|/\xi]},~{}~{}~{}~{}~{}\text{with
}$ (10)
$\displaystyle\xi^{-1}=\log\left[\frac{2\gamma+\mu_{1}+2\sqrt{\gamma^{2}+\mu_{1}\gamma}}{\mu_{1}}\right].$
(11)

We find in Figure (3) that this analytical form shows a good agreement with
numerical simulations except for small $\gamma$.

  * •

We next study the temporal behaviour of the MSD $\langle z_{0}^{2}(t)\rangle$
of a tagged particle. Similar to MRAP, we find a crossover of $\langle
z_{0}^{2}(t)\rangle$ from a linear growth at small times to a sub-diffusive
growth $(\sim\sqrt{t})$ at large times as

$\langle
z_{0}^{2}(t)\rangle\simeq\begin{cases}\mu_{2}a^{2}t,~{}~{}~{}~{}~{}~{}~{}\text{for
small}~{}t\\\ \zeta_{1}\sqrt{t}+\zeta_{2},~{}~{}\text{for
large}~{}t.\end{cases}$ (12)

However, in contrast to the Markov case, constants $\zeta_{1}$ and $\zeta_{2}$
are found to depend on the higher-order correlation functions for the ARAP as

$\displaystyle\zeta_{1}=\sqrt{\frac{\mu_{1}}{\pi}}\left[2aC_{I}+T_{I}+\frac{a\mu_{2}}{\mu_{1}-\mu_{2}}\left\\{a+2C_{1}(t\to\infty)\right\\}\right],$
(13)
$\displaystyle\text{where}~{}~{}C_{I}=\sum_{i=-\infty}^{\infty}C_{i}(t\to\infty),~{}~{}~{}T_{I}=\sum_{i=-\infty}^{\infty}T_{i}(t\to\infty),$
(14)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}\text{and}~{}~{}\zeta_{2}=-\frac{\mu_{2}\zeta_{1}}{4(\mu_{1}-\mu_{2})}~{}\sqrt{\frac{\pi}{\mu_{1}}}.$
(15)

Here $C_{I}$ and $T_{I}$ are constants that depend on the large time
saturation values of the two-point correlation $C_{i}(t)$ and the three-point
correlation $T_{i}(t)$ defined in Eq.(29). In the Markovian limit
$(\gamma\to\infty)$, both of them vanish and we recover Eq. (1) for $\langle
z_{0}^{2}(t)\rangle$. However, for finite $\gamma$, these higher order
correlations are non-vanishing and we observe a substantial enhancement of the
MSD in comparison to the non-persistent case (see Figure (6)).

  * •

We mentioned earlier that the hierarchy of the correlation functions does not
close. This is also seen in the time evolution equation of the three-point
correlations $T_{i}(t)$ which reveals that they depend on the four-point
correlations which in turn, depend on the higher-point correlations. We
numerically demonstrate that any decoupling approximation to break this
hierarchy such as decomposing the four-point functions into lower point
correlation functions does not provide a good approximation (see A).

  * •

We also compute the position correlation $g_{i}(t)=\langle
z_{0}(t)z_{i}(t)\rangle$ and find the following scaling behaviour at large
times:

$\displaystyle
g_{i}(t)\simeq\zeta_{1}\sqrt{t}~{}f\left(\frac{|i|}{\sqrt{4\mu_{1}t}}\right),$
(16)

where $f(y)$ is the same scaling function given in Eq. (2). Once again we
notice that while the scaling function $f(y)$ is same as the Markov case [61],
the pre-factor $\zeta_{1}$ in Eq. (16) is different and therefore, carries the
effect of the persistent dynamics of the active particles.

  * •

The previous results are derived for the quenched uniform initial condition
where the initial positions of the particles remain fixed for different
realisations. We also investigate the variance and the correlations with the
steady state initial condition. To achieve this, we first evolve the system
till time $t_{0}$ and then start measuring the position till further time
$(t_{0}+t)$. Observe that the position of the particle at the onset of the
measurement is different for different realisations. Taking $t_{0}\to\infty$,
we obtain the MSD and correlation in the steady state to be

$\displaystyle l_{0}(t)$
$\displaystyle=\lim_{t_{0}\to\infty}\langle\left[x_{0}(t_{0}+t)-x_{0}(t_{0})\right]^{2}\rangle\simeq\zeta_{1}\sqrt{2t},$
(17) $\displaystyle l_{i}(t)$
$\displaystyle=\lim_{t_{0}\to\infty}\langle\left[x_{i}(t_{0}+t)-x_{i}(t_{0})\right]\left[x_{0}(t_{0}+t)-x_{0}(t_{0})\right]\rangle,$
$\displaystyle\simeq\zeta_{1}\sqrt{2t}~{}f\left(\frac{|i|}{\sqrt{2\mu_{1}t}}\right).$
(18)

Both these results are valid only at large times. Once again, compared to the
MRAP, the persistent nature only affects the coefficient $\zeta_{1}$ but does
not change the sub-diffusive exponent. Another interesting observation is that
for ARAP also, the ratio of the MSDs in the annealed and quenched initial
settings $l_{0}(t)/g_{0}(t)$ is equal to $\sqrt{2}$ at large times, a
condition valid for the Markov case [61]. This means while both $l_{0}(t)$ and
$g_{0}(t)$ individually change due to the persistent dynamics, their ratio is
still fixed to the value $\sqrt{2}$, same as the Markov case, even at finite
$\gamma$.

##  3 Correlation $C_{i}(t)=\langle z_{i}(t)\sigma_{0}(t)\rangle$

Let us begin by computing the spin-position correlation function
$C_{i}(t)=\langle z_{i}(t)\sigma_{0}(t)\rangle$ which will be useful in
calculating the correlations and fluctuations of the tagged particles later.
First note that due to the symmetry of our model, one also has
$C_{i}(t)=\langle z_{0}(t)\sigma_{i}(t)\rangle$. In order to evaluate the
evolution of $C_{i}(t)$ in a small time interval $dt$, let us look at the
different contributions to $C_{i}(t+dt)=\langle
z_{i}(t+dt)\sigma_{0}(t+dt)\rangle$. From this expression, it is clear that
the change in $C_{i}(t+dt)$ can occur either due to the change in
$z_{i}(t+dt)$ or in $\sigma_{0}(t+dt)$. Following Eq. (6), we know that for a
small $dt$, $\sigma_{0}(t+dt)=-\sigma_{0}(t)$ with probability $\gamma dt$
while $\sigma_{0}(t+dt)=\sigma_{0}(t)$ with the complementary probability
$(1-\gamma dt)$. This gives rise to the following change in $C_{i}(t+dt)$:

$\displaystyle C_{i}(t+dt)$ $\displaystyle=-\gamma dt~{}\langle
z_{i}(t)\sigma_{0}(t)\rangle+(1-\gamma dt)~{}\langle
z_{i}(t+dt)\sigma_{0}(t)\rangle.$ (19)

Next, we plug the expression of $z_{i}(t+dt)$ from Eq. (8) and retain all
terms up to linear order in $dt$. It turns out that one then gets different
contributions in $C_{i}(t+dt)$ depending on whether $i=0$ or $i\neq 0$. In
particular, for the latter case, we obtain

$\displaystyle C_{i}(t+dt)$ $\displaystyle~{}=~{}C_{i}(t)-2\gamma
C_{i}(t)dt+\frac{\mu_{1}dt}{2}\left[C_{i+1}(t)+C_{i-1}(t)-2C_{i}(t)\right]$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\frac{\mu_{1}dt}{2}\left[\langle\sigma_{0}(t)\sigma_{i}(t)z_{i+1}(t)\rangle-\langle\sigma_{0}(t)\sigma_{i}(t)z_{i-1}(t)\rangle\right],$
(20)

where $\mu_{k}=\langle\eta^{k}\rangle=\int_{0}^{1}\eta^{k}R(\eta)d\eta$. On
the other hand for $i=0$, similar treatment yields

$\displaystyle C_{0}(t+dt)$ $\displaystyle~{}=~{}C_{0}(t)-2\gamma
C_{0}(t)dt+\frac{\mu_{1}dt}{2}\left[C_{1}(t)+C_{-1}(t)-2C_{0}(t)+2a\right].$
(21)

Now by combining Eqs. (20) and (21), one can appropriately write $C_{i}(t+dt)$
for any value of $i$ as

$\displaystyle C_{i}(t+dt)$ $\displaystyle~{}=~{}C_{i}(t)-2\gamma
C_{i}(t)dt+\frac{\mu_{1}dt}{2}\left[C_{i+1}(t)+C_{i-1}(t)-2C_{i}(t)+2a\delta_{i,0}\right]$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\frac{\mu_{1}dt}{2}\left[\langle\sigma_{0}(t)\sigma_{i}(t)z_{i+1}(t)\rangle-\langle\sigma_{0}(t)\sigma_{i}(t)z_{i-1}(t)\rangle\right].$
(22)

Taking $dt\to 0$ limit, one arrives at the following equation for $C_{i}(t)$:

$\displaystyle\frac{dC_{i}(t)}{dt}$ $\displaystyle=-2\gamma
C_{i}(t)+\frac{\mu_{1}}{2}\left[C_{i+1}(t)+C_{i-1}(t)-2C_{i}(t)+2a\delta_{i,0}\right]$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}+\frac{\mu_{1}}{2}\left[\langle\sigma_{0}(t)\sigma_{i}(t)z_{i+1}(t)\rangle-\langle\sigma_{0}(t)\sigma_{i}(t)z_{i-1}(t)\rangle\right].$
(23)

Figure 2: Comparison of the analytical (mean-field) expression of the
correlation $C_{i}(t)=\langle z_{i}(t)\sigma_{0}(t)\rangle$ in Eq. (25) with
the numerical simulation for (a) $i=0,~{}\gamma=1$ (left panel) and (b)
$i=0,~{}\gamma=2$ (right panel). In both panels, insets show the same
comparison for $i=1$ with the same set of parameters. Simulation is conducted
with $N=201$ particles and the averaging is done over $10^{7}$ realisations.

While this is an exact time evolution equation, it is not closed due to the
presence of higher point correlations. In fact, this turns out to be a general
property of the persistent case that the dynamics of any correlation function
requires knowledge of higher order correlation functions. This makes the
problem analytically challenging. However, as often done, one can make
progress by performing the mean field approximation under which we break the
three point correlation in Eq. (23) into a product of lower order
correlations. The validity of this approximation will be discussed later.
Proceeding ahead, Eq. (23) then simplifies to

$\frac{dC_{i}(t)}{dt}=-2\gamma
C_{i}(t)+\frac{\mu_{1}}{2}\left[C_{i+1}(t)+C_{i-1}(t)-2C_{i}(t)+2a\delta_{i,0}\right].$
(24)

One needs to solve this equation with the initial condition $C_{i}(0)=0$
because $z_{i}(0)=0$ by definition. In B, we have explicitly solved it and
obtained the expression of $C_{i}(t)$ as

$C_{i}(t)=\mu_{1}a\int_{0}^{t}e^{-(2\gamma+\mu_{1})\tau}~{}I_{|i|}(\mu_{1}\tau)~{}d\tau.$
(25)

where $I_{i}(y)$ denotes the modified Bessel function. We also see
$C_{i}(t)=C_{-i}(t)$ since the dynamics of the $0^{\rm th}$ particle
experiences (statistically) same contributions from particles on its either
sides. It turns out that for later calculations, one needs to specify
$C_{i}(t\to\infty)$ which can be easily computed from Eq. (25). Once again we
refer to B for details on this calculation and quote only the final result
here as

$\displaystyle
C_{i}(t\to\infty)=\frac{\mu_{1}a}{2\sqrt{\gamma^{2}+\mu_{1}\gamma}}~{}\text{exp}{[-|i|/\xi]},~{}~{}~{}~{}~{}\text{with
}$ (26)
$\displaystyle\xi^{-1}=\log\left[\frac{2\gamma+\mu_{1}+2\sqrt{\gamma^{2}+\mu_{1}\gamma}}{\mu_{1}}\right].$
(27)

Figure 3: Numerical verification of the correlation $C_{i}(t)=\langle
z_{i}(t)\sigma_{0}(t)\rangle$ at $t\to\infty$ for different values of
$\gamma$. Solid lines represent the analytical (mean-field) expression in Eq.
(26) and symbols are the simulation data. For very small $\gamma$, we observe
deviation of the numerical data from the theoretical expression because the
mean-field approximation (used in the analytical calculation) becomes less and
less valid as $\gamma$ becomes small. The numerical data is obtained using
$10^{7}$ realisations.

In Figures (2) and (3), we have compared our analytical results based on the
mean field approximations with numerical simulations for different values of
$\gamma$. From this comparison, we find that Eq. (26) matches with the
numerics only for moderate and large values of $\gamma$ [see Figure (3)].
However, for small $\gamma$, our results deviate significantly as seen for
$\gamma=0.1$. This is because, at smaller values of $\gamma$, the effect of
activity is so strong that the mean field (decoupling) approximation fails and
one cannot really neglect the three-point (connected) correlation. We will
delve more into this in Section 6. Moreover, observe that the correlation
$C_{i}(t\to\infty)$ decays to zero in the Markov $(\gamma\to\infty)$ limit.
However, for any finite $\gamma$, it possesses a substantial non-zero value.
In what follows, we show that the knowledge of the spin-position correlation
$C_{i}(t)$ is essential to compute the fluctuations and correlations of the
displacements of the tagged particles for active random average process.

##  4 Mean squared displacement and correlations in the quenched initial
condition

We now look at the mean squared displacement and the equal time correlations
of the positions of the tagged particles when their initial positions are
fixed as given in Eq. (3). However, the initial spin $\sigma_{i}(0)$ can still
fluctuate and take values $\sigma_{i}(0)=\pm 1$ with equal probability $1/2$
independently for individual particles. First notice that due to the
translational symmetry in our model, the correlation $\langle
z_{i}(t)z_{j}(t)\rangle$ will depend only on the separation $|i-j|$.
Therefore, without any loss of generality, we put $j=0$ and denote the
correlation $\langle z_{0}(t)z_{i}(t)\rangle$ by $g_{i}(t)$.

To derive the time evolution equation for $g_{i}(t)$, we follow the same
procedure as done for $C_{i}(t)$ in Section 3. At a small time interval
$[t,t+dt]$, we evaluate $g_{i}(t+dt)=\langle z_{0}(t+dt)z_{i}(t+dt)\rangle$
using the update rule in Eq.(8). Keeping all terms up to linear order in $dt$,
one finds the following evolution equation for $i\neq 0$

$\displaystyle
g_{i}(t+dt)=g_{i}(t)+\mu_{1}dt\Big{[}g_{i+1}(t)+g_{i-1}(t)-2g_{i}(t)+2aC_{i}(t)+T_{i}(t)\Big{]},$
(28)

where $C_{i}(t)=\langle z_{i}(t)\sigma_{0}(t)\rangle=\langle
z_{0}(t)\sigma_{i}(t)\rangle$, and $T_{i}(t)$ is a combination of the three-
point correlations defined as

$\displaystyle T_{i}(t)=\frac{1}{2}\Big{[}\langle
z_{0}(t)z_{i+1}(t)\sigma_{i}(t)\rangle-\langle
z_{-1}(t)z_{i}(t)\sigma_{0}(t)\rangle+\langle
z_{1}(t)z_{i}(t)\sigma_{0}(t)\rangle-\langle
z_{0}(t)z_{i-1}(t)\sigma_{i}(t)\rangle\Big{]}.$ (29)

On the other hand, the same procedure for $i=0$ yields

$\displaystyle g_{0}(t+dt)=$
$\displaystyle~{}g_{0}(t)+\mu_{1}dt\Big{[}g_{1}(t)+g_{-1}(t)-2g_{0}(t)+2aC_{0}(t)+T_{0}(t)\Big{]}$
$\displaystyle+\mu_{2}dt\left[a^{2}+2g_{0}(t)-2g_{1}(t)+2aC_{1}(t)-2aC_{0}(t)-T_{0}(t)\right].$
(30)

Combining Eqs. (28) and (30) and taking $dt\to 0$ limit, we find

$\displaystyle\frac{dg_{i}(t)}{dt}=$
$\displaystyle\mu_{1}\Big{[}g_{i+1}(t)+g_{i-1}(t)-2g_{i}(t)+2aC_{i}(t)+T_{i}(t)\Big{]}$
$\displaystyle+\delta_{i,0}~{}\mu_{2}\left[a^{2}+2g_{0}(t)-2g_{1}(t)+2aC_{1}(t)-2aC_{0}(t)-T_{0}(t)\right],$
(31)

for all $i=-\infty,...,-1,0,1,...\infty$. For the original MRAP, the
corresponding equation was derived in [61]. Unlike in the Markov case, once
again we find that Eq. (31) does not satisfy the closure property and involves
higher order correlations in the form of $T_{i}(t)$. Also the function
$C_{i}(t)$ needs the knowledge of higher order correlations as illustrated in
Eq. (23). Overall this makes the computation of $g_{i}(t)$ for the persistent
case rather challenging. To proceed, we perform the joint Fourier-Laplace
transformation

$\displaystyle\bar{g}(q,t)$ $\displaystyle=\sum_{i=-\infty}^{\infty}e^{\iota
iq}~{}g_{i}(t),~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\mathcal{G}(q,s)=\int_{0}^{\infty}dt~{}e^{-st}~{}\bar{g}(q,t),$
(32) $\displaystyle\bar{C}(q,t)$
$\displaystyle=\sum_{i=-\infty}^{\infty}~{}e^{\iota
iq}C_{i}(t),~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\mathcal{C}(q,s)=\int_{0}^{\infty}dt~{}e^{-st}~{}\bar{c}(q,t),$
(33) $\displaystyle\bar{T}(q,t)$
$\displaystyle=\sum_{i=-\infty}^{\infty}e^{\iota
iq}~{}T_{i}(t),~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\mathcal{T}(q,s)=\int_{0}^{\infty}dt~{}e^{-st}~{}\bar{T}(q,t),$
(34)

where $\iota^{2}=-1$ and plug them in Eq. (31) to obtain

$\displaystyle\mathcal{G}(q,s)$
$\displaystyle=\frac{\mu_{2}}{s+\beta(q)}\left[\frac{a^{2}}{s}+2\left(\tilde{g}_{0}(s)-\tilde{g}_{1}(s)\right)+2a\left(\tilde{C}_{1}(s)-\tilde{C}_{0}(s)\right)-\tilde{T}_{0}(s)\right]$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\frac{\mu_{1}}{s+\beta(q)}\left[2a~{}\mathcal{C}(q,s)+\mathcal{T}(q,s)\right],$
(35)

where $\beta(q)=2\mu_{1}(1-\cos(q))$ and the functions $\tilde{g}_{i}(s)$,
$\tilde{C}_{i}(s)$ and $\tilde{T}_{i}(s)$ denote the Laplace transformations
of $g_{i}(t)$, $C_{i}(t)$ and $T_{i}(t)$ respectively. For given
$\tilde{C}_{i}(s)$ and $\tilde{T}_{i}(s)$, Eq. (35) has two unknowns, namely
$\tilde{g}_{0}(s)$ and $\tilde{g}_{1}(s)$. To get rid of $\tilde{g}_{1}(s)$,
we take the Laplace transformation of Eq. (31) for $i=0$ and get

$\tilde{g}_{0}(s)-\tilde{g}_{1}(s)=\frac{a\mu_{2}\tilde{C}_{1}(s)}{\mu_{1}-\mu_{2}}+a\tilde{C}_{0}(s)+\frac{\tilde{T}_{0}(s)}{2}+\frac{\mu_{2}a^{2}}{2s(\mu_{1}-\mu_{2})}-\frac{s~{}\tilde{g}_{0}(s)}{2(\mu_{1}-\mu_{2})}.$
(36)

Inserting this in Eq. (35), we obtain

$\displaystyle\mathcal{G}(q,s)$
$\displaystyle=\frac{\mu_{2}(\mu_{1}-\mu_{2})^{-1}}{s+\beta(q)}\left[\frac{\mu_{1}a^{2}}{s}+2a\mu_{1}\tilde{C}_{1}(s)-s\tilde{g}_{0}(s)\right]+\frac{\mu_{1}}{s+\beta(q)}\left[2a\mathcal{C}(q,s)+\mathcal{T}(q,s)\right].$
(37)

We now have to specify only $\tilde{g}_{0}(s)$ to calculate
$\mathcal{G}(q,s)$. However, the Laplace transform $\tilde{g}_{0}(s)$ can be
obtained self-consistently from Eq. (37) by integrating $\mathcal{G}(q,s)$
with respect to $q$. Then plugging back $\tilde{g}_{0}(s)$ in Eq. (37) gives
the function $\mathcal{G}(q,s)$ exactly. Below, we analyse this equation first
to calculate the mean squared displacement $g_{0}(t)$ and then the equal time
correlation.

###  4.1 Mean-squared displacement $g_{0}(t)=\langle z_{0}^{2}(t)\rangle$

The MSD can be obtained by taking the inverse Fourier transform
$\tilde{g}_{0}(s)=\int_{-\pi}^{\pi}\frac{dq}{2\pi}~{}\mathcal{G}(q,s)$ with
the expression of $\mathcal{G}(q,s)$ given in Eq. (37). We then obtain

$\displaystyle\tilde{g}_{0}(s)\left[1+\frac{\mu_{2}sY(s)}{\mu_{1}-\mu_{2}}\right]=$
$\displaystyle\mu_{1}W(s)+\frac{\mu_{1}\mu_{2}~{}Y(s)}{\mu_{1}-\mu_{2}}\left[\frac{a^{2}}{s}+2a\tilde{C}_{1}(s)\right],$
(38) $\displaystyle\text{where}~{}~{}~{}~{}~{}Y(s)=$
$\displaystyle~{}\frac{1}{2\pi}\int_{-\pi}^{\pi}\frac{dq}{s+\beta(q)}=\frac{1}{\sqrt{s^{2}+4\mu_{1}s}},$
(39) $\displaystyle\text{and}~{}~{}~{}~{}~{}W(s)=$
$\displaystyle~{}\frac{1}{2\pi}\int_{-\pi}^{\pi}dq~{}\frac{2a\mathcal{C}(q,s)+\mathcal{T}(q,s)}{s+\beta(q)}.$
(40)

This equation formally gives the exact MSD of the position of the tagged
particle in the Laplace domain given the two and three point correlations
$C_{i}(t)$ and $T_{i}(t)$ are known.

Figure 4: (a) Comparison of the mean squared displacement $\langle
z_{0}^{2}(t)\rangle$ with the numerical simulation for small values of $t$.
Solid lines represent the analytical expression in Eq. (41) and symbols are
from simulation. (b) Numerical verification of the crossover behaviour of
$\langle z_{0}^{2}(t)\rangle$ from linear growth at small times to sub-
diffusive $\sim\sqrt{t}$ growth at late times. The corresponding expressions
are given in Eq. (41) for small $t$ and in Eq. (48) for large $t$. For panels
(a) and (b), the MSD of a tagged particle is computed using $10^{7}$ and
$10^{6}$ realizations respectively.

Though the functions $\mathcal{C}(q,s)$ and $\mathcal{T}(q,s)$ are not known
exactly, one can still derive some scaling behaviors of $g_{0}(t)$ for
different values of $t$. For small $t$, one has $C_{i}(t\to 0)=0$ and
$T_{i}(t\to 0)=0$ because at very small times, the displacement is negligibly
small. In the Laplace domain, this implies that both
$\left[s\tilde{C}_{i}(s)\right]$ and $\left[s\tilde{T}_{i}(s)\right]$ converge
to zero for large values of $s$. Hence, the function $W(s)$ in Eq. (40) decays
faster than $\sim 1/s$ as $s$ becomes large. On the other hand, from Eq. (39),
we see $Y(s)\simeq 1/s$ for large $s$. Using these approximations in Eq. (38),
we obtain for large $s$,
$\tilde{g}_{0}\left(s\right)~{}\simeq~{}\mu_{2}a^{2}/s^{2}$ which in the time
domain gives

$g_{0}(t)=\langle
z_{0}^{2}(t)\rangle~{}\simeq~{}\mu_{2}a^{2}t,~{}~{}~{}~{}\text{as }t\ll 1.$
(41)

This linear growth of the MSD at small times has been numerically verified in
Figure (4) (left panel). It is easy to understand the small $t$ asymptotic
from the following physical reasoning. At small times, the leading order
contribution to the MSD comes from those realisations where the tagged
particle has made one jump while the other particles have not moved at all.
The probability of observing such an event is $t$ ($pt$ in term of the
unscaled time). Now the particle jumps by a random amount $\pm\eta_{0}a$
depending on its spin $\sigma_{0}(0)$. Since $\sigma_{0}(0)=\pm 1$ with equal
probability $1/2$, we obtain the MSD $\langle
z_{0}^{2}(t)\rangle=\mu_{2}a^{2}\times t=\mu_{2}a^{2}t$ as given in Eq. (41).

Figure 5: Numerical plot of $\sum_{i}T_{i}(t)$ for $\gamma=1$, where
$T_{i}(t)$ is defined in Eq. (29). At late times, we find that this sum
saturates to a non-zero value $\sum_{i}T_{i}(t)\simeq 0.438$ which is used in
Eq. (13) while evaluating $g_{0}(t)$. We have used $10^{7}$ realisations to
obtain the numerical data.

Next we focus on the large-$t$ behaviour of the MSD from Eq. (38). Recall that
for the MRAP, the MSD scales sub-diffusively as $\sim\sqrt{t}$ with a
prefactor that depends on the model parameters [see Eq. (1)]. To see this for
the non-Markov active case, we perform the small-$s$ expansion of
$\tilde{g}_{0}(s)$ in Eq. (38) for which we need to specify $Y(s)$,
$\tilde{C}_{1}(s)$ and $W(s)$ for smaller values of $s$. Expression of $Y(s)$
for $s\ll 1$ follows easily from Eq. (39) as $Y(s)\simeq 1/\sqrt{4\mu_{1}s}$.
As mentioned earlier, due to the hierarchical dependence of the correlations
it is difficult to find the small-$s$ behaviours of $\tilde{C}_{1}(s)$ and
$W(s)$. We, however, numerically observe that the correlations $C_{i}(t)$ and
$\sum_{i}T_{i}(t)$ at late times saturate to their $i$-dependent values and
become time independent. This is numerically illustrated in Figures (3) and
(5). Hence for small $s$, one gets $s\tilde{C}_{1}(s)\simeq C_{1}(t\to\infty)$
and

$\displaystyle
W(s)\simeq\frac{Y(s)}{s}~{}\left[2aC_{I}+T_{I}\right],~{}~{}(\text{for }s\ll
1)$ (42)

where we identify

$\displaystyle\begin{split}\bar{C}(q\to
0,t\to\infty)=\sum_{i}C_{i}(t\to\infty)=C_{I},\\\ \bar{T}(q\to
0,t\to\infty)=\sum_{i}T_{i}(t\to\infty)=T_{I}\end{split}$ (43)

as defined in Eq. (14). These constants can be obtained from the saturation
values of $C_{i}(t)$ and $\sum_{i}T_{i}(t)$ at large $t$ [see Figures. (3) and
(5)]. Furthermore from Eq. (39), it is easy to see that
$Y(s)\simeq\frac{1}{\sqrt{4\mu_{1}s}}$ for small $s$. Plugging these small
$s$-asymptotics in Eq. (38), we find that the Laplace transform
$\tilde{g}_{0}(s)$ reads

$\displaystyle\tilde{g}_{0}(s)\simeq\frac{\sqrt{\pi}\zeta_{1}}{2s^{3/2}},~{}~{}~{}~{}\text{as
}s\to 0,$ (44)

with the constant $\zeta_{1}$ given explicitly in Eq. (13). We emphasize that
this expression is exact at late times and does not involve any approximation.
However, it can be simplified further by using the approximate expressions of
$C_{i}(t)$ at large $t$ given in Eq. (26). Under this approximation, it is
easy to compute $C_{I}=\bar{C}(q\to 0,t\to\infty)$ by taking $t\to\infty$
limit of Eq. (78) at $q=0$ and one finds

$C_{I}\simeq\frac{\mu_{1}a}{2\gamma},$ (45)

plugging which in Eq. (13), one finds the following simpler expression for
$\zeta_{1}$

$\zeta_{1}=\simeq\sqrt{\frac{\mu_{1}}{\pi}}\left[\frac{a^{2}\mu_{1}}{\gamma}+T_{I}+\frac{a^{2}\mu_{2}}{\mu_{1}-\mu_{2}}\left\\{1+\frac{\mu_{1}}{\sqrt{\gamma^{2}+\mu_{1}\gamma}}~{}\text{exp}{[-1/\xi]}\right\\}\right],$
(46)

where $\xi^{-1}$ is given in Eq.(11). It is now straightforward to perform the
inverse Laplace transformation of Eq. (44) and obtain

$\displaystyle g_{0}(t)=\langle
z_{0}^{2}(t)\rangle\simeq\zeta_{1}\sqrt{t}~{}~{}~{}~{}(\text{for large }t).$
(47)

Figure 6: Comparison of analytical (mean-field) expression of the mean squared
displacement $g_{0}(t)=\langle z_{0}^{2}(t)\rangle$ given in Eq. (48) with the
same obtained from numerical simulation for $\gamma=1$ (left panel) and
$\gamma=2$ (right panel). To contrast our result, we have also plotted
$\langle z_{0}^{2}(t)\rangle$ for the Markov case whose expression is given in
Eq. (1) [solid black line labeled by MRAP]. For comparison, we have taken
$T_{I}$ is $0.438$ for left panel and $0.216$ for right panel. The steady-
state value of the correlator $C_{i}(t\to\infty)$ is approximated by the
formula as given in Eq. (26). For both plots, simulation is performed with
$N=501$ particles.

This gives the leading order contribution to the MSD at large times. One can
also obtain the sub-leading term which just turns out to be a constant. To
maintain continuity of our presentation, we relegate this calculation to C and
present only the final result as

$\displaystyle g_{0}(t)=\langle
z_{0}^{2}(t)\rangle\simeq\zeta_{1}\sqrt{t}+\zeta_{2}~{}~{}~{}~{}(\text{for
large }t),$ (48)

with $\zeta_{1}$ and $\zeta_{2}$ given in Eqs. (13) and (15), respectively. In
conjunction to the Markov case, we find that the MSD for the persistent active
case also scales sub-diffusively as $\sim\sqrt{t}$ at large times. Similar
temporal growth patterns have been consistently observed in other active
particle systems. For instance, the large-time growth of current fluctuations
in run-and-tumble particles is the same as that of Brownian particles[43, 69].
However the associated coefficients $\zeta_{1}$ are different for two cases.
While for $\gamma\to\infty$, the two coefficients converge, we see a clear
difference between them for finite $\gamma$. This difference is also
illustrated in Figure (6) where we have also compared with the numerical
simulations. This implies that at large times, the tracer particle performs
sub-diffusion with exponent $1/2$ for both Markov as well as the active RAP.
But the persistent nature of the active particles enhances the coefficient of
the MSD. For the dynamics of a tracer particle in an active single-file
system, such a difference in the MSD with respect to the Markov case was also
numerically observed recently in [47]. Herein, we are able to establish this
analytically based on the mean field approximations. Expectedly, this
approximation breaks down for small values of $\gamma$ and one then needs to
consider the exact form of $\zeta_{1}$. Later, we show that Eq. (44) still
remains valid for small $\gamma$ and a number of results can still be derived.

To summarize, we have shown that the MSD $\langle z_{0}^{2}(t)\rangle$
exhibits a crossover from the diffusive scaling at small times to the sub-
diffusive $(\simeq\zeta_{1}\sqrt{t})$ scaling at large times with the
coefficient $\zeta_{1}$ containing the effect of the persistent motion of
active particles. This crossover behaviour has been shown in Figure (4) (right
panel). Before we end this section, we remark that we have used the random
sequential update scheme to perform the numerical simulation. For
completeness, we have provided the details of this scheme in D. We use this
numerical strategy to verify other analytical results throughout the paper.

###  4.2 Correlation $g_{i}(t)=\langle z_{0}(t)z_{i}(t)\rangle$

We now look at the equal time correlation of the positions of two tagged
particles. For this, we have to perform the inverse Fourier transformation
$\tilde{g}_{i}(s)=\int_{-\pi}^{\pi}\frac{dq}{2\pi}~{}e^{-i\iota
q}~{}\mathcal{G}(q,s)$ for arbitrary $i$ and use $\mathcal{G}(q,s)$ from Eq.
(37). We then obtain

$\displaystyle\tilde{g}_{i}(s)=$
$\displaystyle\mu_{1}~{}W_{i}(s)+\frac{\mu_{2}~{}Y_{i}(s)}{\mu_{1}-\mu_{2}}\left[\frac{\mu_{1}a^{2}}{s}+2a\mu_{1}\tilde{C}_{1}(s)-s\tilde{g}_{0}(s)\right],$
(49) $\displaystyle\text{where}~{}~{}~{}~{}~{}Y_{i}(s)=$
$\displaystyle\frac{1}{2\pi}\int_{-\pi}^{\pi}dq~{}\frac{e^{-\iota
iq}}{s+\beta(q)}$ (50) $\displaystyle\text{and}~{}~{}~{}~{}~{}~{}~{}W_{i}(s)=$
$\displaystyle\frac{1}{2\pi}\int_{-\pi}^{\pi}dq~{}e^{-\iota
iq}~{}\left[\frac{2a\mathcal{C}(q,s)+\mathcal{T}(q,s)}{s+\beta(q)}\right],$
(51)

with $\mathcal{C}(q,s)$ and $\mathcal{T}(q,s)$ defined in Eq. (33) and Eq.
(34), respectively. As discussed for the MSD, carrying out the integration
over $q$ in these expressions turns out to be difficult. However, for small
$s$ (which corresponds to large $t$ in the time domain), one can still perform
this integration approximately which then substantially simplifies the
expression of $\tilde{g}_{i}(s)$ in Eq. (49). In E, we have shown that the
functions $Y_{i}(s)$ and $W_{i}(s)$ behave as

$\displaystyle
Y_{i}(s)\simeq\frac{1}{\sqrt{4\mu_{1}s}}~{}\exp\left[-|i|\sqrt{\frac{s}{\mu_{1}}}\right],$
(52) $\displaystyle
W_{i}(s)\simeq\frac{Y_{i}(s)}{s}\left[2a~{}C_{I}+T_{I}\right],$ (53)

for small $s$. In addition to these quantities, we also need
$\tilde{g}_{0}(s)$ and $\tilde{C}_{1}(s)$ to evaluate $\tilde{g}_{i}(s)$ in
Eq. (49). For this, we use Eq. (44) to get $\tilde{g}_{0}(s)\sim s^{-3/2}$ and
Eq. (26) to get $\tilde{C}_{1}(s)\simeq C_{1}(t\to\infty)/s$ for smaller
values of $s$. Using these approximations in Eq. (49), we get the leading
order behaviour of $\tilde{g}_{i}(s)$ as

$\displaystyle\tilde{g}_{i}(s)\simeq\frac{\sqrt{\pi}\zeta_{1}}{2s^{3/2}}~{}\text{exp}\left[-|i|\sqrt{\frac{s}{\mu_{1}}}\right],~{}~{}~{}~{}~{}\text{for
}s\ll 1.$ (54)

To get the correlation $g_{i}(t)$, we now use the following standard Laplace
transformation [61]:

$\displaystyle\int_{0}^{\infty}dt~{}e^{-st}\sqrt{\frac{t}{\pi}}\left[e^{-\frac{b^{2}}{4t}}-\frac{b\sqrt{\pi}}{2\sqrt{t}}~{}\text{erfc}\left(\frac{b}{2\sqrt{t}}\right)\right]=\frac{e^{-b\sqrt{s}}}{2s^{3/2}},~{}~{}~{}~{}~{}\text{with
}b\geq 0.$ (55)

Plugging this in Eq. (54), we find that $g_{i}(t)$ satisfies the scaling
relation

$\displaystyle g_{i}(t)=\langle
z_{0}(t)z_{i}(t)\rangle\simeq~{}\zeta_{1}\sqrt{t}~{}f\left(\frac{|i|}{\sqrt{4\mu_{1}t}}\right),$
(56)

with the scaling function $f(y)$ given in Eq. (2). Once again, we observe that
the correlation $g_{i}(t)$ is characterised by the same scaling function
$f(y)$ as the MRAP in Eq. (1). However, in conjunction to the MSD, here also
the signature of activity is found in the coefficient $\zeta_{1}$ that appears
in the scaling relation. This means while the scaling function $f(y)$ is same
for the two cases, the overall scaling form is slightly different for any
finite $\gamma$. In Figure (7), we have numerically illustrated this scaling
behaviour for three different values of $t$ and for two different values of
$\gamma$. For all cases, the simulation data converge to Eq. (2) under
appropriate scaling.

Figure 7: Illustration of the scaling behaviour of the correlation
$g_{i}(t)=\langle z_{0}(t)z_{i}(t)\rangle$ in Eq. (56) for $\gamma=1$ and
$\gamma=2$. Symbols are the simulation data for three different times which
converge with the theoretical scaling function $f(y)$ in Eq. (2). For both
plots, simulation has been done with $N=501$ particles.

##  5 Mean squared displacement and correlations in the annealed initial
condition

In the previous sections, we calculated the MSD and the correlations of the
positions of tagged particles in the quenched case during which their initial
positions are fixed to Eq. (3) for all realisations. For this case, we saw
that the persistence nature of the active particles has effects on the
dynamics of a tagged particle even at large times. In this section, we are
interested in carrying out this analysis for the annealed case where we assume
the initial positions are chosen from the stationary state of the system.
Consequently, the initial positions fluctuate across different realisations.
Various studies for the single-file motion of passive particles have shown
that the fluctuations in initial positions have a long term effect on the
dynamics of a tracer particle [10, 9]. In particular, the MSD in the annealed
initial condition at late times is $\sqrt{2}$ times that in the quenched
initial condition. In the remaining of our paper, we address two main
questions: (a) How does MSD for the persistent particles behave in the
annealed case? (b) Are the MSDs for two cases related as for the non-
persistent single-file systems?

In order to study the annealed case, we follow the standard technique which we
briefly discuss here [9]. Starting from the quenched initial condition, we
evolve the system up to a time $t_{0}$ and start measuring the position till
further time $(t_{0}+t)$. We then define the following two-time correlation
function of two tagged particles:

$\displaystyle h_{i}\left(t_{0},t_{0}+t\right)$
$\displaystyle=\langle\left[x_{i}(t_{0}+t)-x_{i}(t_{0})\right]\left[x_{0}(t_{0}+t)-x_{0}(t_{0})\right]\rangle,$
(57)
$\displaystyle=\langle\left[z_{i}(t_{0}+t)-z_{i}(t_{0})\right]\left[z_{0}(t_{0}+t)-z_{0}(t_{0})\right]\rangle,$
(58) $\displaystyle=g_{i}(t_{0}+t)+g_{i}(t_{0})-2S_{i}(t_{0},t_{0}+t),$ (59)

where $S_{i}(t_{0},t_{0}+t)=\langle z_{0}(t_{0})z_{i}(t_{0}+t)\rangle=\langle
z_{i}(t_{0})z_{0}(t_{0}+t)\rangle$. The later equality can be easily proved by
appropriately translating and reflecting the index $i$. In the limit
$t_{0}\to\infty$, the two-time correlation function
$h_{i}\left(t_{0},t_{0}+t\right)$ reduces to the position correlation
$l_{i}(t)$ in the annealed initial condition, i.e.

$\displaystyle l_{i}(t)$
$\displaystyle=\lim_{t_{0}\to\infty}h_{i}\left(t_{0},t_{0}+t\right),$ (60)
$\displaystyle=\lim_{t_{0}\to\infty}\left[g_{i}(t_{0}+t)+g_{i}(t_{0})-2S_{i}(t_{0},t_{0}+t)\right].$
(61)

This means that we do not need to perform the averaging over the initial
positions and simply use the results for the quenched case. Since by
performing the time shift by $t_{0}$ and taking $t_{0}\to\infty$, we are
effectively putting the system to its steady state, we anticipate the two
methods to be equivalent. In what follows, we use the relation (61) and
compute the correlation and fluctuation in the annealed setting. As clear from
this relation, this reduces to calculating the correlation
$S_{i}(t_{0},t_{0}+t)$ which we carry out below.

For this, we again start with the evolution of
$S_{i}(t_{0},t_{0}+t+dt)=\langle z_{0}(t_{0})z_{i}(t_{0}+t+dt)\rangle$ in a
small time interval $[t,t+dt]$ and use Eq. (8) to plug $z_{i}(t_{0}+t+dt)$.
Following same steps as before, the dynamics of $S_{i}(t_{0},t_{0}+t)$ can be
shown to be

$\displaystyle\frac{\partial S_{i}(t_{0},t_{0}+t)}{\partial t}$
$\displaystyle=\frac{\mu_{1}}{2}\left[S_{i+1}(t_{0},t_{0}+t)+S_{i-1}(t_{0},t_{0}+t)-2S_{i}(t_{0},t_{0}+t)+2aC_{i}(t_{0},t_{0}+t)\right]$
$\displaystyle+\frac{\mu_{1}}{2}\left[\langle
z_{0}(t_{0})\sigma_{i}(t_{0}+t)z_{i+1}(t_{0}+t)\rangle-\langle
z_{0}(t_{0})\sigma_{i}(t_{0}+t)z_{i-1}(t_{0}+t)\rangle\right].$ (62)

where we have defined $C_{i}(t_{0},t_{0}+t)=\langle
z_{0}(t_{0})\sigma_{i}(t_{0}+t)\rangle$. Again, we see that this equation is
not closed and involves higher order correlations. Under the mean field
approximation, we break these higher correlations as a product of lower order
correlations. Unlike in the equal time case, this approximation turns out to
be quite good here since at large $t$, these three point correlations rapidly
decay to zero. Thus within this approximation, Eq. (62) can be simplified as

$\displaystyle\frac{\partial S_{i}(t_{0},t_{0}+t)}{\partial t}$
$\displaystyle\simeq\frac{\mu_{1}}{2}\left[S_{i+1}(t_{0},t_{0}+t)+S_{i-1}(t_{0},t_{0}+t)-2S_{i}(t_{0},t_{0}+t)+2aC_{i}(t_{0},t_{0}+t)\right].$
(63)

Figure 8: Numerical verification of the MSD $l_{0}(t)$ with the steady state
initial condition for two values of $\gamma$. The theoretical expression is
given in Eq. (67). For simulation, we have first allowed the system to evolve
till time $t_{0}=3000$ and then start measuring the position with $N=501$
particles.

Now to solve this equation, we have to find the source term
$C_{i}(t_{0},t_{0}+t)$ on the right hand side. Once again this term can be
easily calculated using the update rule for $\sigma_{i}(t)$ in Eq. (6). Since
this derivation is exactly same as the previous ones, we have presented it in
F and quote only the final result as

$\displaystyle C_{i}(t_{0},t+t_{0})=C_{i}(t_{0})~{}e^{-2\gamma t},$ (64)

where $C_{i}(t_{0})=\langle z_{0}(t_{0})\sigma_{i}(t_{0})\rangle$ is given in
Eq. (25). We now have all terms in the right hand side of Eq. (63) which can
now be straightforwardly solved by taking the joint Fourier-Laplace
transformation. As shown in G, we obtain the expression of
$S_{i}(t_{0},t_{0}+t)$ as

$\displaystyle
S_{i}(t_{0},t_{0}+t)\simeq\zeta_{1}\sqrt{t_{0}}-\frac{\zeta_{1}\sqrt{t}}{\sqrt{2}}~{}\mathcal{W}\left(\frac{|i|}{\sqrt{2\mu_{1}t}}\right),$
(65)

where the constant $\zeta_{1}$ is given in Eq. (13) and the function
$\mathcal{W}(y)$ is defined as

$\displaystyle\mathcal{W}(y)=e^{-y^{2}}+\sqrt{\pi}y~{}\rm{erf}(y).$ (66)

We emphasize that the expression of $S_{i}(t_{0},t_{0}+t)$ in Eq. (65) is
valid only for large $t_{0}$ and large $t$ but with their ratio $t/t_{0}$
fixed to a value much smaller than $1$. Below we use this form in Eq. (61) to
compute the asymptotic behaviours of the MSD and correlation in the annealed
initial setting.

###  5.1 Mean squared displacement $l_{0}(t)$

Using the result in Eq. (65), we get
$S_{0}(t_{0},t_{0}+t)\simeq\zeta_{1}\sqrt{t_{0}}-\zeta_{1}\sqrt{t}/\sqrt{2}$.
Plugging this in Eq. (61) along with $g_{0}(t_{0})$ from Eq. (47) for large
$t_{0}$, we obtain the MSD $l_{0}(t)$ as

$\displaystyle
l_{0}(t)\simeq\zeta_{1}\sqrt{2t},~{}~{}~{}~{}~{}~{}\text{for~{}~{}~{}}t\gg 1.$
(67)

We have verified this analytic expression against simulations in Figure (8).
Few remarks are in order. First, the MSD of a tracer in ARAP again scales sub-
diffusively as $\sim\sqrt{t}$ at large times (reminiscent of the Markov case).
However, due to the persistent motion of active particles, once again we see
that the coefficient accompanying this sub-diffusive growth is different than
the Markov case. Only for $\gamma\to\infty$, the coefficient converges to the
non-persistent value. Second interesting observation is that while MSDs in the
annealed and quenched initial settings change from their Markov counterparts,
their ratio is still equal to $\sqrt{2}$, a result seen for many single-file
systems [10, 9, 61]. This implies that even though $l_{0}(t)$ and $g_{0}(t)$
individually carry the signature of activity, their ratio however is still
fixed to the value $\sqrt{2}$ even for finite $\gamma$.

Figure 9: Numerical verification of the scaling function $f(y)$ in Eq. (68)
for the correlation $l_{i}(t)$ in the steady state. In both panels, we have
performed the comparison for three different times shown in different colours.
We observe that the data for different times (symbols) converge to the
theoretical curve (solid line) under scaling with respect to time. For
simulation, we have first evolved the system till time $t_{0}=3000$ with
$N=501$ and then start measuring the position.

###  5.2 Correlation $l_{i}(t)$

We next look at the expression of $l_{i}(t)$ in Eq. (61) for general $i$ and
insert $S_{i}(t_{0},t_{0}+t)$ from Eq. (65) and $g_{i}(t_{0})$ from Eq. (56).
The correlation then turns out to be

$\displaystyle l_{i}(t)$ $\displaystyle\simeq
2\zeta_{1}\sqrt{t_{0}}~{}f\left(\frac{|i|}{\sqrt{4\mu_{1}t_{0}}}\right)-2\zeta_{1}\sqrt{t_{0}}+\zeta_{1}\sqrt{2t}~{}\mathcal{W}\left(\frac{|i|}{\sqrt{2\mu_{1}t}}\right),$
$\displaystyle\simeq-\frac{\sqrt{\pi}\zeta_{1}|i|}{\sqrt{\mu_{1}}}~{}\text{erfc}\left(\frac{|i|}{\sqrt{4\mu_{1}t_{0}}}\right)+\zeta_{1}\sqrt{2t}~{}\mathcal{W}\left(\frac{|i|}{\sqrt{2\mu_{1}t}}\right),$
$\displaystyle\simeq-\frac{\sqrt{\pi}\zeta_{1}|i|}{\sqrt{\mu_{1}}}+\zeta_{1}\sqrt{2t}~{}\mathcal{W}\left(\frac{|i|}{\sqrt{2\mu_{1}t}}\right),$
$\displaystyle\simeq\zeta_{1}\sqrt{2t}~{}f\left(\frac{|i|}{\sqrt{2\mu_{1}t}}\right),$
(68)

where the scaling function $f(y)$ is given in Eq. (2). Also, in going from
second line to the third line, we have used the asymptotic behaviour of
complementary error function as
$\text{erfc}\left(\frac{|i|}{\sqrt{4\mu_{1}t_{0}}}\right)\to 1$ as
$t_{0}\to\infty$ for finite $i$. In Figure (9), we have compared the scaling
behaviour of $l_{i}(t)$ with the numerical simulation for $\gamma=1$ in left
panel and $\gamma=2$ in right panel. For both panels, we have carried out the
comparison for three different values of $t$. We observe excellent match of
our analytical results with the simulation for all cases. Compared to the
Markov case, once again we see that the persistence only changes the prefactor
in the scaling relation (68) but keeps the form of the scaling function the
same.

##  6 Effect of small $\gamma$ on the MSD

In the previous sections, we looked at the fluctuations and correlations of
the positions of tagged particles and studied their dependence on the initial
condition. Moreover, based on a mean field approximation, we provided semi-
analytic expressions of these quantities. However, it turns out that this
approximation is valid only for large and intermediate values of the flipping
rate $\gamma$ and breaks down for its smaller values. To illustrate this, we
have plotted the simulation data of the MSD $g_{0}(t)/\sqrt{t}$ for
$\gamma=0.5$ and $\gamma=0.25$ in Figure (10). Furthermore, we compare our
numerical result with its analytic form in Eq. (13) by computing the prefactor
$\zeta_{1}$ in two ways: (i) first we obtain $C_{i}(t)$ from numerics and plug
it in Eq. (13) to get a complete numerical estimate of $\zeta_{1}$, (ii)
second we use the approximated theoretical form of $\zeta_{1}$ in Eq. (46). As
seen in Figure (10), while $\zeta_{1}$ for case (i) matches with the
simulation data, there is a clear departure of $\zeta_{1}$ obtained for case
(ii). This departure becomes more pertinent as we go to smaller and smaller
values of $\gamma$. Hence, we still find that, $g_{0}(t)$ scales sub-
diffusively as $\sim\sqrt{t}$ at late times even for small $\gamma$ with
$\zeta_{1}$ given by its exact form in Eq. (13).

Figure 10: Numerical verification of the prefactor $\zeta_{1}$ associated with
the variance $g_{0}(t)\simeq\zeta_{1}\sqrt{t}$. We have also performed a
comparison with the theoretical mean-field expression in Eq. (46) (shown in
blue) and the exact expression in Eq. (13) (shown in green) with $C_{i}(t)$
obtained from simulations. Based on the simulation, we have obtained (a)
$C_{I}=0.485$ and $T_{I}=0.91$ (for the left panel) and (b) $C_{I}=0.91$ and
$T_{I}=2.08$ (for the right panel). Simulation for both plots has been
conducted with $N=1001$ particles and the averaging is done over $5\times
10^{5}$ realizations.

We next look at the MSD $l_{0}(t)$ in the annealed case where the initial
positions are drawn from the steady state. For large and intermediate
$\gamma$, we proved before that the ratio $l_{0}(t)/g_{0}(t)$, at late times,
is still equal to $\sqrt{2}$, a relation true for many single-file systems. In
the remaining part of this section, we analyze how this ratio changes for
smaller values of $\gamma$. Extensive numerical simulations suggest that this
ratio is still equal to the factor $\sqrt{2}$ even for small $\gamma$. For
example, in Figure (11), we have shown the simulation data for the ratio
$l_{0}(t)/g_{0}(t)$ for $\gamma=0.5$ and $\gamma=0.25$. For both cases, we
find that the ratio approaches the value $\sqrt{2}$ at late times. This means
while both $l_{0}(t)$ and $g_{0}(t)$ deviate from their mean field forms,
their ratio is still fixed to $\sqrt{2}$. To prove this, we first have to
evaluate the behavior of $S_{0}(t_{0},t_{0}+t)$ [see Eq. (61)]. Rewriting the
time evolution equation for $S_{i}(t_{0},t_{0}+t)$ in Eq. (62), we get

$\displaystyle\frac{\partial S_{i}(t_{0},t_{0}+t)}{\partial t}$
$\displaystyle=\frac{\mu_{1}}{2}\left[S_{i+1}(t_{0},t_{0}+t)+S_{i-1}(t_{0},t_{0}+t)-2S_{i}(t_{0},t_{0}+t)+2aU_{i}(t_{0},t_{0}+t)\right],$
(69)

where the function $U_{i}(t_{0},t_{0}+t)$ denotes

$\displaystyle U_{i}(t_{0},t_{0}+t)=$ $\displaystyle\frac{1}{2a}\left[\langle
z_{0}(t_{0})\sigma_{i}(t_{0}+t)z_{i+1}(t_{0}+t)\rangle-\langle
z_{0}(t_{0})\sigma_{i}(t_{0}+t)z_{i-1}(t_{0}+t)\rangle\right]$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+C_{i}(t_{0},t_{0}+t).$
(70)

Figure 11: Numerical comparison of the ratio of the MSD $l_{0}(t)$ measured in
the steady state and the MSD $g_{0}(t)$ with quenched uniform initial state
for two small values of $\gamma$. We observe that even for small $\gamma$,
where the predictions from mean field approximation do not hold, the ratio
$l_{0}(t)/g_{0}(t)$, at large $t$ approaches, $\sqrt{2}$ implying that the
relation in Eq. (73) to be valid for all non-zero $\gamma$. For the annealed
case, we start measuring position only after time $t_{0}=10000$ in this
simulation. We have performed averaging over $5\times 10^{5}$ realisations.

Due to the presence of $\sigma$-variable, we anticipate $U_{i}(t_{0},t_{0}+t)$
to decay, at late times, as

$\displaystyle
U_{i}(t_{0},t_{0}+t)~{}~{}\overset{\mathrm{t_{0}\to\infty}}{\simeq}~{}~{}\psi_{i}~{}e^{-\chi_{i}t},~{}~{}~{}~{}\text{for
}t\gg 1,$ (71)

where both $\psi_{i}$ and $\chi_{i}$ are time-independent functions of the
index $i$. Under mean-field approximation, we explicitly showed that
$\psi_{i}=C_{i}(t_{0}\to\infty)$ and $\chi_{i}=2\gamma$ (see Eq. (64)). Here,
we have carried out extensive numerical simulations to find out that
$\chi_{0}=2\gamma$ and $\chi_{1}=\gamma/4$. We have illustrated this in Figure
12 for $\gamma=0.5$ and $\gamma=0.25$. Proceeding with this form, we have
shown in H that even for small $\gamma$ we obtain

$\displaystyle
S_{0}\left(t_{0},t_{0}+t\right)\simeq\zeta_{1}\sqrt{t_{0}}-\zeta_{1}\sqrt{\frac{t}{2}},$
(72)

for both $t_{0}$ and $t$ large. We emphasize that this expression is valid for
all non-zero values of $\gamma$ with $\zeta_{1}$ given exactly in Eq. (13).
Plugging this in Eq. (61) yields

$\displaystyle l_{0}(t)\simeq\sqrt{2}~{}g_{0}(t),$ (73)

with $g_{0}(t)=\zeta_{1}\sqrt{t}$. This means that while both $l_{0}(t)$ and
$g_{0}(t)$ depart individually from their mean field forms, the ratio is still
given by $\sqrt{2}$. We have numerically verified this result in Figure (11)
for $\gamma=0.5$ and $\gamma=0.25$. Demonstrating this result numerically for
very small values of $\gamma$ turns out to be computationally expensive.
Recall from Eq. (61) that one needs to go to very large $t_{0}$ to measure
$l_{0}(t)$. Numerically, we see that smaller the value of $\gamma$, larger is
the value of $t_{0}$ that one has to consider. On the other hand at very large
time, the boundary effects start to become important which alter the MSD. In
order to avoid boundary effects, we have to increase the number of particles
in the simulation, which makes the computation highly expensive. In our study,
we have fixed the smallest value as $\gamma=0.25$. Already for this value, we
observe departure of the simulation data from mean field results.

Figure 12: Numerical plot of the correlation $U_{i}(t_{0},t_{0}+t)$ in Eq.
(70) for $i=0$ as a function of time for (a) $\gamma=0.5$ and (b)
$\gamma=0.25$ with $t_{0}=10000$, $N=1001$ and $5\times 10^{5}$ realisations.
For large $t$, we find that the simulation data (shown in red) can be fitted
(shown in black) by
$\lim_{t_{0}\to\infty}U_{0}(t_{0},t_{0}+t)\simeq\psi_{0}~{}e^{-2\gamma t}$
with $\psi_{0}=0.95$ for panel (a) and $\psi_{0}=1.8$ for panel (b). In the
respective insets, we have shown the corresponding plots of
$U_{1}(t_{0},t_{0}+t)$ along with the fits
$\lim_{t_{0}\to\infty}U_{1}(t_{0},t_{0}+t)\simeq\psi_{1}~{}e^{-\gamma t/4}$
with $\psi_{1}=0.009$ for $\gamma=0.5$ and $\psi_{1}=0.04$ for $\gamma=0.25$.
For this case, we have chosen $t_{0}=200$, $N=101$ and $10^{9}$ realisations.

##  7 Conclusion

In this paper, we have investigated the motion of tracer particles in the
random average process of persistent active particles in an infinite line.
Using mean field approximation, we calculated the mean squared displacement
and correlation of the positions of tracer particles both in the quenched
initial condition and in the steady state. In particular, for the quenched
case, we showed that the MSD exhibits a crossover from a diffusive scaling at
small times to a sub-diffusive $(\sim\sqrt{t})$ scaling at late times.
Interestingly we find that the coefficient associated with this sub-diffusive
growth is different from the corresponding non-persistent case and the two
converge only in the limit $\gamma\to\infty$. For finite $\gamma$, we see a
clear difference between them as illustrated in Figure (6). Similarly, for the
position correlation, we find slight difference in Eq. (56) compared to the
Markov case. While the overall scaling function $f(y)$ in Eq. (56) is same as
the MRAP, the prefactor $\zeta_{1}$ is different and therefore carries the
effect of the activity.

Next, we studied these quantities in the steady state where we first evolve
the system till time $t_{0}\to\infty$ and then start measuring the positions.
Unlike in the quenched case, here the positions at the onset of the
measurement fluctuate for different realisations. For this case, we
analytically showed that the MSD at late times grows sub-diffusively as
$\sim\sqrt{t}$ with the associated coefficient once again different than the
Markov case. Only for $\gamma\to\infty$, the two become equal. Quite
remarkably while both MSDs in the quenched initial condition and in the steady
state individually change due to the persistent nature of the particles, their
ratio is still equal to $\sqrt{2}$ at large times. This is a common result
known to be true for many single-file systems [61, 9, 10]. Our study here
reveals it to be valid even for the ARAP for all non-zero values of $\gamma$.
Finally, we calculated the correlation between positions of two tagged
particles at steady state in Eq. (68).

Most of our derivations in this paper relied on the mean-field approximation.
Due to the non-Markovian nature of the noise, the equations for correlation
functions do not satisfy closure property and involve higher order
correlations. Moreover, we observe that these higher order correlations
possess non-negligible values. To compute these in Eq. (23) for
$C_{i}(t)=\langle z_{i}(t)\sigma_{0}(t)\rangle$, we have used mean-field
approximation. For this case, we found that the cross-correlation involved in
this approximation is negligible for moderate and large tumbling rate
$\gamma$. However, for the three-point correlation $T_{i}(t)$ in A, such
approximation fails even for moderate $\gamma$. Only in the Markovian limit
($\gamma\to\infty$), this can be neglected.

Solving single-file motion for active particles is a notoriously challenging
problem. Here, we showcased an example of active single-file motion for which
we could derive many results analytically. Carrying out this study for active
particles with hardcore exclusions is an interesting and important direction
to explore. Recent numerical studies in this direction have pointed out at
some interesting qualitative differences than the usual single-file diffusion
[47]. Proving this analytically is still an open problem. For systems obeying
diffusive hydrodynamics, the coefficient of the sub-diffusive growth of the
MSD of a tracer particles is specified by the diffusivity and mobility of the
system [9]. In short range interacting systems such transport coefficients are
usually determined by the two point correlations [65] as in the Markov RAP
case where only $\mu_{1}$ and $\mu_{2}$ appears in the expression of the MSD
[see Eq. (1)]. In contrast, for our active RAP system we observe that the MSD
gets contribution from higher point correlations also. It would be interesting
to investigate if it is possible to derive the same MSD form as in Eq. (12)
from hydrodynamic evolution for the density [70], more precisely for the
inter-particle separation field similar to the Markov RAP case [64]. Also, in
this work, we have only looked at the MSD and the two-point correlation
functions. Obtaining higher moments and distribution of the position of a
tagged particle are interesting problems even for the Markov RAP. Finally it
would be interesting to study the effect of local biases in the single-file
model of active particles in the same spirit as in [14, 64].

## Acknowledgements

The authors thank Tirthankar Banerjee for stimulating discussions on the
paper. A. K. acknowledges the support of the core research grant no.
$\text{CRG}/2021/002455$ and MATRICS grant $\text{MTR}/2021/000350$ from the
Science and Engineering Research Board (SERB), Department of Science and
Technology, Government of India. P. S. acknowledges the support of Novo
Nordisk Foundation grant $\text{NNF}21\text{OC}0071284$. The authors also
acknowledge the support from the Department of Atomic Energy, Government of
India, under Project No. RTI4001.

##  Appendix A Failure of mean field approximation to specify $T_{i}(t)$ in
Eq. (29)

In this appendix, we will demonstrate why the mean field approximation does
not correctly characterize the correlation $T_{i}(t)$ in Eq. (29). Looking at
this expression, we see that we have to calculate three-point correlation
functions like $\langle z_{0}(t)z_{i+1}(t)\sigma_{i}(t)\rangle$ and $\langle
z_{0}(t)z_{i-1}(t)\sigma_{i}(t)\rangle$. To see if the mean field approach
works for this case, we define a general correlation
$\mathcal{T}_{ij}(t)\equiv\langle z_{0}(t)z_{i+1}(t)\sigma_{j}(t)\rangle$ and
write its time evolution equations as

$\displaystyle\frac{d\mathcal{T}_{ij}(t)}{dt}=$
$\displaystyle-2\gamma\mathcal{T}_{ij}(t)+\frac{\mu_{1}}{2}\Big{[}\mathcal{T}_{i+1,j}(t)+\mathcal{T}_{i-1,j}(t)+\mathcal{T}_{i+1,j+1}(t)+\mathcal{T}_{i-1,j-1}(t)-4\mathcal{T}_{i,j}(t)\Big{.}$
$\displaystyle\Big{.}+2a\\{\langle\sigma_{0}z_{i+1}(t)\sigma_{j}(t)\rangle+\langle\sigma_{i+1}(t)z_{0}(t)\sigma_{j}(t)\rangle\Big{]}+\frac{\mu_{1}}{2}\Big{[}\langle
z_{0}(t)\sigma_{i+1}(t)z_{i+2}(t)\sigma_{j}(t)\rangle\Big{.}$
$\displaystyle\Big{.}-\langle
z_{0}(t)\sigma_{i+1}(t)z_{i}(t)\sigma_{j}(t)\rangle+\langle
z_{1}(t)\sigma_{0}(t)z_{i+1}(t)\sigma_{j}(t)\rangle-\langle
z_{-1}(t)\sigma_{0}(t)z_{i+1}(t)\sigma_{j}(t)\rangle\Big{]}$
$\displaystyle+\mu_{2}\Big{[}2\langle
z_{0}(t)z_{0}(t)\sigma_{j}(t)\rangle-\langle
z_{1}(t)z_{1}(t)\sigma_{j}(t)\rangle-\langle
z_{-1}(t)z_{-1}(t)\sigma_{j}(t)\rangle\Big{.}$ $\displaystyle\Big{.}-\langle
z_{0}(t)z_{1}(t)\sigma_{0}(t)\sigma_{j}(t)\rangle+\langle
z_{0}(t)z_{-1}(t)\sigma_{0}(t)\sigma_{j}(t)\rangle\Big{]}.$ (74)

Note that we are interested in calculating $T_{i}(t)$ which is obtained by
putting $i=j$ in $\mathcal{T}_{ij}(t)$. Coming to Eq. (74), we observe that it
does not satisfy the closure property as it contains four-point correlation
functions. Once again, we use mean field approximations to break the four-
point correlation in terms of lower point correlations as

$\displaystyle\langle
z_{0}(t)\sigma_{i+1}(t)z_{i+2}(t)\sigma_{j}(t)\rangle\simeq$
$\displaystyle\langle z_{0}(t)\sigma_{i+1}(t)\rangle\langle
z_{i+2}(t)\sigma_{j}(t)\rangle+\langle z_{0}(t)\sigma_{j}(t)\rangle\langle
z_{i+2}(t)\sigma_{i+1}(t)\rangle$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\langle
z_{0}(t)z_{i+2}(t)\rangle\langle\sigma_{i+1}(t)\sigma_{j}(t)\rangle.$ (75)

Figure 13: Comparison of the numerically obtained four-point correlation
function $P_{0}(t)$ as defined in Eq. (76) with its approximated value
$C_{1}(t)(C_{2}(t)-C_{0}(t))$ using mean field. For both $\gamma=1$ and
$\gamma=2$, the two deviate substantially indicating that the mean field
approximation does not work for $P_{0}(t)$. For this simulation, we have used
$4\times 10^{6}$ realisations.

With this approximation, the following four-point point correlation function
appearing in Eq. (74) becomes

$\displaystyle P_{0}(t)$ $\displaystyle=\Big{[}\langle
z_{0}(t)\sigma_{1}(t)z_{2}(t)\sigma_{0}(t)\rangle-\langle
z_{0}(t)\sigma_{1}(t)z_{0}(t)\sigma_{0}(t)\rangle+\langle
z_{1}(t)\sigma_{0}(t)z_{1}(t)\sigma_{0}(t)\rangle\Big{.}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\Big{.}-\langle
z_{-1}(t)\sigma_{0}(t)z_{1}(t)\sigma_{0}(t)\rangle\Big{]},$
$\displaystyle\simeq C_{1}(t)\big{[}C_{2}(t)-C_{0}(t)\big{]}.$ (76)

We now test the validity of this approximation. For this, we measure both
$P_{0}(t)$ and $C_{1}(t)\big{[}C_{2}(t)-C_{0}(t)\big{]}$ from the numerical
simulations and compare them in Figure (13) for $\gamma=1$ and $\gamma=2$. For
both cases, we observe that $P_{0}(t)$ has a large positive value whereas
$C_{1}(t)\big{[}C_{2}(t)-C_{0}(t)\big{]}$ has a small negative value. Clearly,
this implies $P_{0}(t)\neq C_{1}(t)\big{[}C_{2}(t)-C_{0}(t)\big{]}$. Hence we
numerically find that the decoupling approximation of breaking four-point
correlation function in terms of two-point correlation functions in Eq. (75)
is not valid and thus the analytical calculation of obtaining three point
correlation functions seems difficult.

##  Appendix B Expression of $C_{i}(t)$ in Eq. (25)

In this appendix, we will derive the expression of $C_{i}(t)$ quoted in Eq.
(25) of the main text. For this, we take the Fourier transform
$\bar{C}(q,t)=\Sigma_{i=-\infty}^{\infty}e^{\iota iq}C_{i}(t)$ (where
${\iota}^{2}=-1$) and insert this in Eq. (24) to obtain

$\displaystyle\frac{d\bar{C}(q,t)}{dt}=-\alpha(q)\tilde{C}(q,t)+\mu_{1}a,$
(77)

where $\alpha(q)=\mu_{1}(1-\cos(q))+2\gamma$. Solving this equation, we get

$\bar{C}(q,t)=\mu_{1}a\left(\frac{1-e^{-\alpha(q)t}}{\alpha(q)}\right).$ (78)

Performing the inverse Fourier transformation yields

$C_{i}(t)=\frac{\mu_{1}a}{2\pi}\int_{-\pi}^{\pi}e^{-\iota
iq}\left(\frac{1-e^{-\alpha(q)t}}{\alpha(q)}\right)dq.$ (79)

Performing the integration over $q$ analytically in this equation is
difficult. However, one can get a simplified expression by proceeding as
follows. Differentiating on both sides of Eq. (79), we get

$\displaystyle\frac{dC_{i}(t)}{dt}$
$\displaystyle=\frac{\mu_{1}a}{2\pi}\int_{-\pi}^{\pi}e^{-\iota
iq}e^{-\alpha(q)t}dq,$
$\displaystyle=\frac{\mu_{1}a}{2\pi}e^{-(2\gamma+\mu_{1})t}\int_{-\pi}^{\pi}\cos(iq)e^{\mu_{1}t\cos(q)}dq,$
$\displaystyle=\mu_{1}ae^{-(2\gamma+\mu_{1})t}I_{|i|}(\mu_{1}t).$ (80)

Next, we integrate both sides with respect to $t$ and use the initial
condition $C_{i}(0)=0$ to obtain

$C_{i}(t)=\mu_{1}a\int_{0}^{t}e^{-(2\gamma+\mu_{1})\tau}I_{|i|}(\mu_{1}\tau)~{}d\tau.$
(81)

In the limit $t\to\infty$, one can exactly carry out the integration over
$\tau$ to get

$\displaystyle C_{i}(t\to\infty)$
$\displaystyle=\mu_{1}a\int_{0}^{\infty}e^{-2\gamma\tau-\mu_{1}\tau}I_{|i|}(\mu_{1}\tau)~{}d\tau,$
$\displaystyle=a\int_{0}^{\infty}e^{-\left(\frac{2\gamma+\mu_{1}}{\mu_{1}}\right)w}I_{|i|}(w)~{}dw,$
$\displaystyle=\frac{\mu_{1}a}{\sqrt{4\gamma^{2}+4\mu_{1}\gamma}}\left(\frac{2\gamma+\mu_{1}+\sqrt{4\gamma^{2}+4\mu_{1}\gamma}}{\mu_{1}}\right)^{-|i|}.$
(82)

This result has been used in Eq. (26) in the main text.

##  Appendix C Sub-leading term in $g_{0}(t)$ in Eq. (48)

In this appendix, we derive the expression of the sub-leading term in the
variance $g_{0}(t)$ for large $t$. As seen in Eq. (48), for large $t$, the
variance scales sub-diffusively as $g_{0}(t)\simeq\zeta_{1}\sqrt{t}$ with
prefactor $\zeta_{1}$ given in Eq. (13). Here, we are interested in
calculating the next order correction which turns out to be a constant. To
derive this, let us quote here the Laplace transform $\tilde{g}_{0}(s)$ from
Eq. (38)

$\displaystyle\tilde{g}_{0}(s)\left[1+\frac{\mu_{2}sY(s)}{\mu_{1}-\mu_{2}}\right]=$
$\displaystyle\mu_{1}W(s)+\frac{\mu_{1}\mu_{2}~{}Y(s)}{\mu_{1}-\mu_{2}}\left[\frac{a^{2}}{s}+2a\tilde{C}_{1}(s)\right],$
(83)

where $Y(s)$ and $W(s)$ are defined in Eqs. (39) and (40) respectively. Note
that $Y(s)\simeq 1/\sqrt{4\mu_{1}s}$ for small $s$, inserting which in Eq.
(83) gives

$\displaystyle\tilde{g}_{0}(s)\simeq\underbrace{\frac{\mu_{1}W(s)}{1+\phi\sqrt{s}}}_{\text{first
term}}+\underbrace{\frac{\mu_{2}\sqrt{\mu_{1}}}{2\sqrt{s}~{}\left(\mu_{1}-\mu_{2}\right)\left(1+\phi\sqrt{s}\right)}~{}\left[\frac{a^{2}}{s}+2a\tilde{C}_{1}(s)\right]}_{\text{second
term}},$ (84)

where $\phi=\frac{\mu_{2}}{\sqrt{4\mu_{1}}~{}(\mu_{1}-\mu_{2})}$. For
computational ease, we have written the two terms separately in the right hand
side. As evident, for the first term, we have to specify the function $W(s)$.
For small $s$, the integrand in Eq. (40) is dominated by smaller values of
$q$. We therefore anticipate the major contribution to the integration to come
from smaller values of $q$. With this approximation, the expression of $W(s)$
reduces to

$\displaystyle
W(s)\simeq~{}\frac{1}{\sqrt{4\mu_{1}s}}\left[2a\mathcal{C}\left(q\to
0,s\right)+\mathcal{T}\left(q\to 0,s\right)\right],$ (85)

and the first term in Eq. (84) for $s\to 0$ becomes

$\displaystyle\text{first
term}\simeq\sqrt{\frac{\mu_{1}}{4}}~{}\left(\frac{1}{\sqrt{s}}-\phi\right)~{}\left[2a\mathcal{C}\left(q\to
0,s\right)+\mathcal{T}\left(q\to 0,s\right)\right].$ (86)

We now have to compute the asymptotic forms of $\mathcal{C}\left(q\to
0,s\right)$ and $\mathcal{T}\left(q\to 0,s\right)$ as $s\to 0$. To do this, we
first recall that both $C_{i}(t)$ and $T_{i}(t)$ involve the
$\sigma$-variables in their definitions. Due to this, both of these functions
$\mathcal{C}\left(q\to 0,s\right)$ and $\mathcal{T}\left(q\to 0,s\right)$
relax exponentially to their steady values with relaxation time scale $\sim
1/\gamma$. For $\mathcal{C}\left(q\to 0,s\right)$, one can show this from Eq.
(78) whereas for $\mathcal{T}\left(q\to 0,s\right)$, one can see this
numerically. In the Laplace domain, this implies

$\displaystyle\mathcal{C}(q\to 0,s)$ $\displaystyle\simeq~{}\frac{\bar{C}(q\to
0,t\to\infty)}{s}+\text{constant term},$ (87) $\displaystyle\mathcal{T}(q\to
0,s)$ $\displaystyle\simeq~{}\frac{\bar{T}(q\to
0,t\to\infty)}{s}+\text{constant term},$ (88)

for small $s$. Here the constant terms do not involve $s$. Inserting these
forms in Eq. (86) and using the identification in Eq. (43), we obtain

first term
$\displaystyle\simeq\frac{\sqrt{\mu_{1}}}{2}~{}\left[2aC_{I}+T_{I}\right]~{}\left(\frac{1}{s^{3/2}}-\frac{\phi}{s}\right),$
(89)

which can be simplified further as

first term
$\displaystyle\simeq\frac{1}{2}\left[\frac{a^{2}(\mu_{1})^{3/2}}{\gamma}+\sqrt{\mu_{1}}~{}T_{I}\right]~{}\left(\frac{1}{s^{3/2}}-\frac{\phi}{s}\right),$
(90)

using the approximation $C_{I}=\frac{\mu_{1}a}{2\gamma}$ as derived in Eq.
(45). So far, we have obtained the approximate form of the first term in Eq.
(84) for small $s$. Next, we carry out the same analysis for the second term.
Looking at Eq. (84), we observe that one needs to calculate the Laplace
transform $\tilde{C}_{1}(s)$ for small $s$. To obtain this, we consider the
expression of $C_{1}(t)$ in Eq. (25) and rewrite it as

$C_{1}(t)=C_{1}(t\to\infty)-\mu_{1}a\int_{t}^{\infty}e^{-(2\gamma+\mu_{1})\tau}~{}I_{1}(\mu_{1}\tau)~{}d\tau.$
(91)

For large $t$, the argument of the Bessel function is also large which enables
us to use the approximation $I_{1}(\mu_{1}\tau)\simeq
e^{\mu_{1}\tau}/\sqrt{2\pi\mu_{1}\tau}$ for $\mu_{1}\tau\gg 1$. With this
approximation, the integration over $\tau$ can be easily carried out and we
get

$C_{1}(t)\simeq
C_{1}(t\to\infty)-a\sqrt{\frac{\mu_{1}}{4\gamma}}~{}\text{erfc}\left(\sqrt{2\gamma
t}\right).$ (92)

Taking the Laplace transformation of this equation gives

$\displaystyle\tilde{C}_{1}(s)$
$\displaystyle\simeq\frac{C_{1}(t\to\infty)}{s}-a\sqrt{\frac{\mu_{1}}{4\gamma}}~{}\frac{1}{s+2\gamma+\sqrt{2\gamma
s+4\gamma^{2}}},$ (93)
$\displaystyle\simeq\frac{C_{1}(t\to\infty)}{s}-\frac{a}{4\gamma}\sqrt{\frac{\mu_{1}}{4\gamma}},~{}~{}~{}~{}\text{as
}s\to 0.$ (94)

Plugging this approximate form in Eq. (84), we obtain the second term as

$\displaystyle\text{second
term}=\frac{a\mu_{2}\sqrt{\mu_{1}}}{2(\mu_{1}-\mu_{2})}\left(a+2C_{1}(t\to\infty)\right)~{}\left(\frac{1}{s^{3/2}}-\frac{\phi}{s}\right).$
(95)

Therefore, we have computed the approximate forms of both terms in Eq. (84)
for small $s$. Using these forms, the expression of the Laplace transform
$\tilde{g}_{0}(s)$ simplifies to

$\displaystyle\tilde{g}_{0}(s)\simeq\frac{\sqrt{\pi}}{2}~{}\zeta_{1}\left(\frac{1}{s^{3/2}}-\frac{\phi}{s}\right).$
(96)

Finally performing the inverse Laplace transformation, we obtain the result
written in Eq. (48) in the main text.

##  Appendix D Details about numerical simulations

This appendix provides details about the numerical simulation adopted to
verify various analytical results in the paper. To begin with, we have
$N(=2n+1)$ number of particles initially placed at a uniform distance apart as

$\displaystyle x_{i}(0)=ia,$ (97)

where $i$ is an integer that lies between $-n\leq i\leq n$. On the other hand,
we choose initial $\sigma_{i}(0)$ from $\pm 1$ with equal probability $1/2$
independently for all particles. This means while the initial positions are
fixed for different realisations, the initial $\sigma_{i}(0)$ still fluctuate.
For a given realisation, we then implement random sequential rule to update
the position of the particles. During a small time interval $[t,t+dt]$, we
choose a random integer $m$ uniformly between $[-n,n]$ and update the position
$x_{m}(t)$ and spin variable $\sigma_{m}(t)$ according to Eqs. (5) and (6)
using the Gillespie algorithm [71]. We then repeat this step for $N$ times.
This implies that in the time interval $dt$, we perform the update rule
randomly for $N$ times. Finally, we iterate this process till the observation
time $t$ is reached. For all figures, we chose $dt=0.01$ except for the
Figures (2) and (3) for which we chose $dt=0.002$. The random variables
$\eta_{i}$ in Eq. (7) are chosen uniformly from $[0,1)$, hence $R(\eta)=1$ for
$0\leq\eta<1$ and zero otherwise. To obtain the numerical data for Figs. 6-9,
we have used $10^{6}$ realisations.

##  Appendix E Expressions of $Y_{i}(s)$ and $W_{i}(s)$ as $s\to 0$

Here, we derive the approximate expressions of $Y_{i}(s)$ and $W_{i}(s)$ in
Eqs. (50) and (51) for small values of $s$. Let us first present the
calculation for $Y_{i}(s)$ which reads

$\displaystyle Y_{i}(s)=$
$\displaystyle\frac{1}{2\pi}\int_{-\pi}^{\pi}dq~{}\frac{e^{-\iota
iq}}{s+\beta(q)},$ (98)

where $\beta(q)=2\mu_{1}(1-\cos(q))$. For $s\to 0$, the integrand in Eq. (98)
diverges as $q\to 0$. Therefore, we expect the major contribution to the
integration should come from the small values of $q$. Taking the $q\to 0$
limit, we get $\beta(q)\simeq\mu_{1}q^{2}$ plugging which in Eq. (98), we get

$\displaystyle Y_{i}(s)\simeq$
$\displaystyle\frac{1}{2\pi}\int_{-\pi}^{\pi}dq~{}\frac{e^{-\iota
iq}}{s+\mu_{1}q^{2}}.$ (99)

Changing the variable $q=\sqrt{s/\mu_{1}}~{}w$ and taking $s\to 0$, we get

$\displaystyle Y_{i}(s)$
$\displaystyle\simeq\frac{1}{2\pi\sqrt{s\mu_{1}}}~{}\int_{-\pi}^{\pi}\frac{dw}{1+w^{2}}~{}\text{exp}\left[-\iota
iw\sqrt{\frac{s}{\mu_{1}}}\right],$ (100)
$\displaystyle\simeq\frac{1}{\sqrt{4\mu_{1}s}}~{}\text{exp}\left[-|i|\sqrt{\frac{s}{\mu_{1}}}\right].$
(101)

This result has been quoted in Eq. (52) which was instrumental in obtaining
the asymptotic behaviour of the MSD and correlation for the positions of the
particles.

It turns out that for $W_{i}(s)$ also, one can proceed similarly to get its
small $s$ behaviour. To see this, let us first rewrite its expression from Eq.
(51)

$\displaystyle W_{i}(s)=$
$\displaystyle\frac{1}{2\pi}\int_{-\pi}^{\pi}dq~{}e^{-\iota
iq}~{}\left[\frac{2a\mathcal{C}(q,s)+\mathcal{T}(q,s)}{s+\beta(q)}\right],$
(102)

where $\mathcal{C}(q,s)$ and $\mathcal{T}(q,s)$ are the joint Fourier-Laplace
transforms given in Eqs. (33) and (34) respectively. As discussed before for
small $s$, the integration in Eq. (102) will be dominated by small values of
$q$ which gives

$\displaystyle W_{i}(s)\simeq$
$\displaystyle\frac{1}{2\pi}\int_{-\pi}^{\pi}dq~{}e^{-\iota
iq}~{}\left[\frac{2a\mathcal{C}(q\to 0,s\to 0)+\mathcal{T}(q\to 0,s\to
0)}{s+\mu_{1}q^{2}}\right].$ (103)

Numerically, we see that both $\bar{C}(q,t)$ and $\bar{T}(q,t)$ in Eqs. (33)
and (34) attain stationary values as $t\to\infty$. This implies, as mentioned
previously, that in the Laplace domain, one gets

$\displaystyle\mathcal{C}(q\to 0,s\to 0)$
$\displaystyle\simeq~{}\frac{\bar{C}(q\to 0,t\to\infty)}{s},$ (104)
$\displaystyle\mathcal{T}(q\to 0,s\to 0)$
$\displaystyle\simeq~{}\frac{\bar{T}(q\to 0,t\to\infty)}{s},$ (105)

Plugging these forms in Eq. (103) gives

$\displaystyle W_{i}(s)\simeq$
$\displaystyle~{}\frac{Y_{i}(s)}{s}~{}\left[2aC_{I}+T_{I}\right].$ (106)

where we have used Eq. (43) and $Y_{i}(s)$ is given in Eq. (99). Inserting the
form of $Y_{i}(s)$ for small $s$ from Eq. (101), we finally get

$\displaystyle
W_{i}(s)\simeq\frac{Y(s)}{s}~{}\left[2aC_{I}+T_{I}\right]~{}~{}\text{exp}\left[-|i|\sqrt{\frac{s}{\mu_{1}}}\right].$
(107)

To summarize, in this appendix, we have derived the forms of $Y_{i}(s)$ and
$W_{i}(s)$ as $s\to 0$. Observe that for $i=0$, the function $Y_{i}(s)$ from
Eq. (98) reduces to $Y(s)$ in Eq. (39). Consequently, the Eq. (107) with $i=0$
provides the small $s$ behaviour of $W(s)$ as quoted in Eq. (42).

##  Appendix F Computation of $C_{i}(t_{0},t_{0}+t)=\langle
z_{0}(t_{0})\sigma_{i}(t_{0}+t)\rangle$

This appendix presents a derivation of the expression of
$C_{i}(t_{0},t_{0}+t)$ quoted in Eq. (64). Let us begin by writing the
dynamics of $C_{i}(t_{0},t_{0}+t+dt)$ in small time interval $dt$. Using the
time evolution of $\sigma_{i}(t)$ in Eq. (6), we have

$\displaystyle C_{i}(t_{0},t_{0}+t+dt)$ $\displaystyle=\langle
z_{0}(t_{0})\sigma_{i}(t_{0}+t+dt)\rangle,$ (108) $\displaystyle=-\gamma
dt\langle z_{0}(t_{0})\sigma_{i}(t_{0}+t)\rangle+(1-\gamma dt)\langle
z_{0}(t_{0})\sigma_{i}(t_{0}+t)\rangle,$ (109)
$\displaystyle=C_{i}(t_{0},t_{0}+t)-2\gamma dt~{}C_{i}(t_{0},t_{0}+t).$ (110)

Taking the $dt\to 0$ limit, we get the dynamics of $C_{i}(t_{0},t_{0}+t)$ as

$\displaystyle\frac{\partial C_{i}(t_{0},t_{0}+t)}{\partial t}=-2\gamma
C_{i}(t_{0},t_{0}+t).$ (111)

In order to solve this equation, we need to specify appropriate initial
condition. Recall that as $t\to 0$, one has $C_{i}(t_{0},t_{0})=\langle
z_{0}(t_{0})\sigma_{i}(t_{0})\rangle$ which is simply $C_{i}(t_{0})$ defined
in Section 3. The solution of $C_{i}(t_{0})$ has been obtained in Eq. (25).
Solving Eq. (111) with this initial condition, we obtain

$\displaystyle C_{i}(t_{0},t+t_{0})=C_{i}(t_{0})~{}e^{-2\gamma t},$ (112)

which has also been quoted in Eq. (64).

##  Appendix G Computation of $S_{i}(t_{0},t_{0}+t)$ in Eq. (65)

In this appendix, we derive the expression of $S_{i}(t_{0},t_{0}+t)$ in Eq.
(65) which was used to calculate the MSD and correlation for the active RAP
with annealed initial condition. To this aim, we take the joint Fourier-
Laplace transformation of $S_{i}(t_{0},t_{0}+t)$ as

$\displaystyle\mathcal{S}(q,s,t)=\sum_{i=-\infty}^{\infty}e^{\iota
iq}~{}\tilde{\mathcal{S}}_{i}(s,t),~{}~{}~{}\text{with
}\tilde{\mathcal{S}}_{i}(s,t)=\int_{0}^{\infty}dt_{0}~{}e^{-st_{0}}~{}S_{i}(t_{0},t_{0}+t),$
(113)

and insert it in Eq. (63) to obtain

$\displaystyle\mathcal{S}(q,s,t)=\mathcal{G}(q,s)e^{-\frac{\beta(q)t}{2}}+\frac{\mu_{1}a~{}\mathcal{C}(q,s)}{\left(\frac{\beta(q)}{2}-2\gamma\right)}\left(e^{-2\gamma
t}-e^{-\frac{\beta(q)}{2}t}\right),$ (114)

where $\mathcal{G}(q,s)$ and $\mathcal{C}(q,s)$ denote, respectively, the
joint Fourier-Laplace transformations of $g_{i}(t_{0})$ and $C_{i}(t_{0})$
given in Eqs. (37) and (78). Also, we have defined
$\beta(q)=2\mu_{1}(1-\cos(q))$. Now to get $S_{i}(t_{0},t_{0}+t)$ from Eq.
(114), one needs to perform two inversions: one is the inverse Fourier
transformation with respect to $q$ and the other is the inverse Laplace
transformation with $s$. Let us first write the inversion with respect to $q$
as

$\displaystyle\tilde{\mathcal{S}}_{i}(s,t)=\int_{-\pi}^{\pi}\frac{dq}{2\pi}~{}e^{-\iota
iq}~{}\mathcal{S}(q,s,t)=\tilde{\mathcal{S}}^{(1)}_{i}(s,t)+\tilde{\mathcal{S}}^{(2)}_{i}(s,t),$
(115)

where the two terms are

$\displaystyle\left.\begin{aligned}
\tilde{\mathcal{S}}^{(1)}_{i}(s,t)&=\int_{-\pi}^{\pi}\frac{dq}{2\pi}~{}e^{-\iota
iq}~{}e^{-\frac{\beta(q)t}{2}}~{}\tilde{\mathcal{G}}(q,s)\\\
\tilde{\mathcal{S}}^{(2)}_{i}(s,t)&=\frac{\mu_{1}a~{}}{2\pi}\int_{-\pi}^{\pi}\frac{dq~{}e^{-\iota
iq}}{\left[\frac{\beta(q)}{2}-2\gamma\right]}~{}~{}\mathcal{C}(q,s)\left(e^{-2\gamma
t}-e^{-\frac{\beta(q)}{2}t}\right).\end{aligned}\right\\}$ (116)

We now proceed to evaluate these two terms separately.

###  G.1 Calculation of $\tilde{\mathcal{S}}^{(1)}_{i}(s,t)$

Since we are interested in the $t_{0}\to\infty$ limit, we evaluate these terms
for small values of $s$. Using the approximation
$\tilde{\mathcal{G}}(q,s)\simeq\zeta_{1}\sqrt{\mu_{1}\pi}/s(s+\mu_{1}q^{2})$
for $s\to 0$ from Eq. (54), the first term in Eq. (116) becomes

$\displaystyle\tilde{\mathcal{S}}^{(1)}_{i}(s,t)$
$\displaystyle\simeq\frac{\sqrt{\mu_{1}}~{}\zeta_{1}}{2\sqrt{\pi}s}\int_{-\pi}^{\pi}\frac{dq}{\left(s+\mu_{1}q^{2}\right)}~{}~{}\text{exp}{\Big{[}-\mu_{1}t(1-\cos(q))-i\iota
q\Big{]}}.$ (117)

Since the integrand is exponentially decaying in $t$, the leading contribution
to the integral comes from the smaller values of $q$ which enables us to
approximate it as

$\tilde{\mathcal{S}}^{(1)}_{i}(s,t)\simeq\frac{\sqrt{\mu_{1}}~{}\zeta_{1}}{2\sqrt{\pi}s}\int_{-\pi}^{\pi}\frac{dq}{\left(s+\mu_{1}q^{2}\right)}~{}~{}\text{exp}{\Big{[}-\frac{\mu_{1}t}{2}q^{2}-\iota
iq\Big{]}}.$ (118)

Under the variable transformation $q=\sqrt{\frac{s}{\mu_{1}}}w$, it simplifies
to

$\displaystyle\tilde{\mathcal{S}}^{(1)}_{i}(s,t)$
$\displaystyle\simeq\frac{\zeta_{1}}{2\sqrt{\pi}s^{3/2}}\int_{-\infty}^{\infty}dw~{}~{}\frac{\text{exp}{\Big{[}-\frac{stw^{2}}{2}-\iota
i\sqrt{\frac{s}{\mu_{1}}}w}\Big{]}}{1+w^{2}},$ (119)
$\displaystyle\simeq\frac{\sqrt{\pi}\zeta_{1}}{4s^{3/2}}\left[e^{\frac{st}{2}-|i|\sqrt{\frac{s}{\mu_{1}}}}~{}\text{erfc}\left(\sqrt{\frac{st}{2}}-\frac{|i|}{\sqrt{2\mu_{1}t}}\right)+e^{\frac{st}{2}+|i|\sqrt{\frac{s}{\mu_{1}}}}\right.$
$\displaystyle\left.~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\times~{}\text{erfc}\left(\sqrt{\frac{st}{2}}+\frac{|i|}{\sqrt{2\mu_{1}t}}\right)\right].$
(120)

Once again we take the $s\to 0$ limit to recast this equation as

$\tilde{\mathcal{S}}^{(1)}_{i}(s,t)\simeq\frac{\sqrt{\pi}\zeta_{1}}{2}\Big{[}\frac{1}{s^{3/2}}-\frac{1}{s}\Big{\\{}\frac{\sqrt{2t}}{\sqrt{\pi}}e^{-\frac{i^{2}}{2\mu_{1}t}}+\frac{i}{\sqrt{\mu_{1}}}{\rm{erf}}\left(\frac{i}{\sqrt{2\mu_{1}t}}\right)\Big{\\}}\Big{]}.$
(121)

Now it is straightforward to perform the inverse Laplace transformation with
respect to $s$. However, before that, we calculate the second term in Eq.
(116).

###  G.2 Calculation of $\tilde{\mathcal{S}}^{(2)}_{i}(s,t)$

Observe that the second term in Eq. (116) depends on the function
$\mathcal{C}(q,s)$ whose expression is given in Eq. (78) in the time domain.
For small $s$, this function simplifies to

$\displaystyle\mathcal{C}(q,s)\simeq\frac{\mu_{1}a}{s\left[\mu_{1}(1-\cos(q))+2\gamma\right]}$
(122)

Plugging this in Eq. (116) gives

$\displaystyle\tilde{\mathcal{S}}^{(2)}_{i}(s,t)\simeq\frac{(\mu_{1}a)^{2}}{2\pi
s}\int_{-\pi}^{\pi}\frac{dq~{}e^{-\iota
iq}}{\left[\frac{\beta(q)}{2}-2\gamma\right]\left[\frac{\beta(q)}{2}+2\gamma\right]}~{}\left(e^{-2\gamma
t}-e^{-\frac{\beta(q)}{2}t}\right).$ (123)

Performing integration over $q$ at this stage turns out to be difficult.
However, we can carry out this by taking an additional Laplace transform with
respect to $t\to\lambda$ under which the above expression becomes

$\displaystyle\mathbb{S}_{i}^{(2)}(s,\lambda)$
$\displaystyle\simeq\frac{(\mu_{1}a)^{2}}{2\pi
s(\lambda+2\gamma)}\int_{-\pi}^{\pi}\frac{dq~{}e^{-\iota
iq}}{\left[2\gamma+\mu_{1}(1-\cos(q))\right]\left[\lambda+\mu_{1}(1-\cos(q))\right]},$
(124) $\displaystyle\simeq\frac{(\mu_{1}a)^{2}}{8\pi
s\gamma^{2}}~{}\int_{-\pi}^{\pi}~{}dq~{}\frac{e^{-\iota
iq}}{\lambda+\mu_{1}q^{2}/2},$ (125)

where we have used the notation $\mathbb{S}_{i}^{(2)}(s,\lambda)$ to denote
the Laplace transform of $\tilde{\mathcal{S}}^{(2)}_{i}(s,t)$. Moreover, in
going to the second line, we have used the approximation that for small
$\lambda$ (which corresponds to large $t$), the integral is dominated by small
values of $q$ which then allows us to write $1-\cos(q)\simeq q^{2}/2$. Next,
we change the variable $q=\sqrt{\frac{2\lambda}{\mu_{1}}}~{}w$ and take
$\lambda\to 0$ to rewrite Eq. (125) as

$\displaystyle\mathbb{S}_{i}^{(2)}(s,\lambda)$
$\displaystyle\simeq\frac{a^{2}~{}\mu_{1}^{3/2}}{4\pi
s\gamma^{2}\sqrt{2\lambda}}~{}\int_{-\infty}^{\infty}dw~{}\frac{e^{-\iota
iw\sqrt{\frac{2\lambda}{\mu_{1}}}}}{1+w^{2}},$ (126)
$\displaystyle\simeq\frac{a^{2}(\mu_{1})^{3/2}}{4s\gamma^{2}\sqrt{2\lambda}}e^{-|i|\sqrt{\frac{2\lambda}{\mu_{1}}}}.$
(127)

Finally performing the inverse Laplace transformation from $\lambda\to t$
yields

$\displaystyle\tilde{\mathcal{S}}^{(2)}_{i}(s,t)$
$\displaystyle\simeq\frac{a^{2}~{}\mu_{1}^{3/2}}{4s\gamma^{2}\sqrt{2\pi
t}}~{}\text{exp}\left(-|i|^{2}/2\mu_{1}t\right).$ (128)

###  G.3 $\tilde{\mathcal{S}}_{i}(s,t)$ in Eq. (115)

Comparing the two terms in Eqs. (121) and (128) respectively, we find that the
leading order contribution to $\tilde{\mathcal{S}}_{i}(s,t)$ in Eq. (115) at
large $t$ comes only from the first term. This allows us to write

$\displaystyle\tilde{\mathcal{S}}_{i}(s,t)\simeq\frac{\sqrt{\pi}\zeta_{1}}{2}\Big{[}\frac{1}{s^{3/2}}-\frac{1}{s}\Big{\\{}\frac{\sqrt{2t}}{\sqrt{\pi}}e^{-\frac{i^{2}}{2\mu_{1}t}}+\frac{i}{\sqrt{\mu_{1}}}\rm{erf}\left(\frac{i}{\sqrt{2\mu_{1}t}}\right)\Big{\\}}\Big{]}.$
(129)

Notice that our result is valid only for large $t_{0}$ and large $t$ but with
the ratio $t/t_{0}\ll 1$. Finally taking the inverse Laplace transformation,
we obtain the expression of $S_{i}(t_{0},t_{0}+t)$ presented in Eq. (65).

##  Appendix H Derivation of $S_{0}\left(t_{0},t_{0}+t\right)$ in Eq. (72)

In this appendix, we derive the expression of
$S_{0}\left(t_{0},t_{0}+t\right)$ in Eq. (72) which is instrumental in
deriving the form of $l_{0}(t)$ in Eq. (73) for small $\gamma$. To this end,
we take the joint Fourier-Laplace transformation of Eq. (69) with respect to
$i~{}(\to q)$ and $t_{0}~{}(\to s)$ to obtain

$\frac{\partial\tilde{\mathcal{S}}(q,s,t)}{\partial
t}=-\frac{\beta(q)}{2}\tilde{\mathcal{S}}(q,s,t)+\mu_{1}a~{}\tilde{\mathcal{U}}(q,s,t),$
(130)

where $\tilde{\mathcal{S}}(q,s,t)$ and $\tilde{\mathcal{U}}(q,s,t)$ are the
joint Fourier-Laplace transformations of $S_{i}(t_{0},t_{0}+t)$ and
$U_{i}(t_{0},t_{0}+t)$ respectively and $\beta(q)=2\mu_{1}(1-\cos q)$. Solving
this equation with the initial condition
$\mathcal{S}(q,s,t=0)=\mathcal{G}(q,s)$, we obtain

$\displaystyle\tilde{\mathcal{S}}(q,s,t)=\mathcal{G}(q,s)e^{-\frac{\beta(q)t}{2}}+\mu_{1}a\int_{0}^{t}d\tau~{}\exp\left[-\frac{\beta(q)}{2}(t-\tau)\right]~{}\tilde{\mathcal{U}}(q,s,\tau).$
(131)

We now carry out the inversion of this expression with respect to $q$ to
obtain

$\displaystyle\mathcal{S}_{0}(s,t)=\underbrace{\int_{-\pi}^{\pi}\frac{dq}{2\pi}~{}e^{-\frac{\beta(q)t}{2}}~{}\mathcal{G}(q,s)}_{\text{first
term}}+\underbrace{\mu_{1}a\int_{-\pi}^{\pi}\frac{dq}{2\pi}~{}\int_{0}^{t}d\tau~{}\exp\left[-\frac{\beta(q)}{2}(t-\tau)\right]~{}\tilde{\mathcal{U}}(q,s,\tau)}_{\text{second
term}}.$ (132)

Now the first term is exactly same as Eq. (116) and has been explicitly
calculated in G to be

$\displaystyle\text{first
term}=\frac{\sqrt{\pi}\zeta_{1}}{2}\left(\frac{1}{s^{3/2}}-\frac{1}{s}\sqrt{\frac{2t}{\pi}}\right).$
(133)

On the other hand, to calculate second term, we further take a Laplace
transformation with respect to $t~{}(\to\lambda)$ to yield

$\displaystyle\text{second
term}=\mu_{1}a\int_{-\pi}^{\pi}\frac{dq}{2\pi}~{}\frac{\tilde{\mathcal{U}}(q,s,\lambda)}{\lambda+\mu_{1}(1-\cos
q)}.$ (134)

Once again, performing integration over $q$ turns out to be difficult since we
do not know the exact form of $\tilde{\mathcal{U}}(q,s,\lambda)$. However, for
small $\lambda$, the integral over $q$ will be dominated by small values of
$q$ and one can then use the asymptotic form of of
$\tilde{\mathcal{U}}(q,s,\lambda)$ in Eq. (71). To see this, we approximate
$1-\cos q\simeq q^{2}/2$ and rewrite (134) for small $\lambda$ as

second term $\displaystyle\simeq\mu_{1}a~{}\tilde{\mathcal{U}}(q\to 0,s\to
0,\lambda\to
0)\int_{-\pi}^{\pi}\frac{dq}{2\pi}~{}\frac{1}{\lambda+\mu_{1}q^{2}/2},$ (135)
$\displaystyle\simeq
a\sqrt{\frac{\mu_{1}}{2\lambda}}~{}\tilde{\mathcal{U}}(q\to 0,s\to
0,\lambda\to 0).$ (136)

Using the ansatz in Eq. (71), one can show that $\tilde{\mathcal{U}}(q\to
0,s\to 0,\lambda\to 0)\sim~{}1/s$ where the prefactor does not depend on
$\lambda$. This means that the second term diverges as $\sim 1/\sqrt{\lambda}$
as $\lambda\to 0$ which in the time domain implies a decay of the form $\sim
t^{-1/2}$. Therefore, for large $t$, we obtain sub-leading contribution of the
second term compared to the first term in Eq. (132) and $\mathcal{S}_{0}(s,t)$
becomes

$\displaystyle\mathcal{S}_{0}(s,t)\simeq\frac{\sqrt{\pi}\zeta_{1}}{2}\left(\frac{1}{s^{3/2}}-\frac{1}{s}\sqrt{\frac{2t}{\pi}}\right).$
(137)

Taking the inverse Laplace transformation, we obtain the expression of
$S_{0}\left(t_{0},t_{0}+t\right)$ which has been quoted in Eq. (72) in the
main text.

## References

## References

  * [1] T. E. Harris, Diffusion with “collisions” between particles, J. Appl. Probab. 2(2), 323–338 (1965).
  * [2] D. W. Jepsen, Dynamics of a Simple Many-Body System of Hard Rods, J. Math. Phys. 6(3), 405 (1965).
  * [3] O. Bénichou, P. Illien, G. Oshanin, A. Sarracino and R. Voituriez, J. Phys.: Condens. Matter 30(44), 443001 (2018).
  * [4] R. Arratia, The motion of a tagged particle in the simple symmetric exclusion system on Z, Ann. Prob. 11(2), 362 (1983).
  * [5] S. Alexander and P. Pincus, Diffusion of labeled particles on one-dimensional chains, Phys. Rev. B. 18, 2011 (1978).
  * [6] C. R$\ddot{\rm o}$odenbeck , J. K$\ddot{\rm a}$arger and K. Hahn, Calculating exact propagators in single-file systems via the reflection principle, Phys. Rev. E. 57, 4382 (1998).
  * [7] E. Barkai and R. Silbey, Theory of Single File Diffusion in a Force Field, Phys. Rev. Lett. 102, 050602 (2009).
  * [8] P. L. Krapivsky, K. Mallick and T. Sadhu, Large Deviations in Single-File Diffusion, Phys. Rev. Lett. 113, 078101 (2014).
  * [9] P. L. Krapivsky, K. Mallick and T. Sadhu, Dynamical properties of single-file diffusion, J. Stat. Mech. P09007 (2015).
  * [10] T. Banerjee, R. L. Jack and M. E. Cates, Role of initial conditions in one-dimensional diffusive systems: Compressibility, hyperuniformity, and long-term memory, Phys. Rev. E 106, L062101 (2022).
  * [11] C. Hegde, S. Sabhapandit and A. Dhar, Universal Large Deviations for the Tagged Particle in Single File Motion, Phys. Rev. Lett. 113, 120601 (2014).
  * [12] A. Roy, O. Narayan, A. Dhar and S. Sabhapandit, Tagged Particle Diffusion in One-Dimensional Gas with Hamiltonian Dynamics, J. Stat. Phys. 150, 851-866 (2013).
  * [13] M. Kollmann, Single-file Diffusion of Atomic and Colloidal Systems: Asymptotic Laws, Phys. Rev. Lett. 90, 180602 (2003).
  * [14] J. Cividini, A. Kundu, S. N. Majumdar and D. Mukamel, Correlation and fluctuation in Random Average Process on an infinite line with a driven tracer, J. Stat. Mech., 053212 (2016).
  * [15] J. Cividini and A. Kundu, Tagged particle in single-file diffusion with arbitrary initial conditions, J. Stat. Mech., 083203 (2017).
  * [16] N. Leibovich and E. Barkai, Everlasting effect of initial conditions on single-file diffusion, Phys. Rev. E. 88, 032107 (2013).
  * [17] S. Ramaswamy, The mechanics and statistics of active matter, Annu. Rev. Condens. Matter Phys. 1, 323–45 (2010).
  * [18] M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao and R. A. Simha, Hydrodynamics of soft active matter, Rev. Mod. Phys. 85, 1143 (2013).
  * [19] M. E. Cates and J. Tailleur, Motility-induced phase separation, Annu. Rev. Condens. Matter Phys. 6, 219–44 (2015).
  * [20] R. Soto and R. Golestanian, Run-and-tumble dynamics in a crowded environment: Persistent exclusion process for swimmers, Phys. Rev. E 89, 012706 (2014).
  * [21] N. Sepúlveda and R. Soto, Coarsening and clustering in run-and-tumble dynamics with short-range exclusion, Phys. Rev. E 94, 022603 (2016).
  * [22] A. P. Solon, Y. Fily, A. Baskaran, M. E. Cates, Y. Kafri, M. Kardar and J. Tailleur, Pressure is not a state function for generic active fluids, Nat. Phys. 11, 673 (2015).
  * [23] K. Malakar, V. Jemseena, A. Kundu, K. V. Kumar, S. Sabhapandit, S. N. Majumdar, S. Redner and A. Dhar, Steady state, relaxation and first-passage properties of a run-and-tumble particle in one-dimension, J. Stat. Mech., 043215 (2018).
  * [24] U. Basu, S. N. Majumdar, A. Rosso and G. Schehr, Active Brownian motion in two dimensions, Phys. Rev. E. 98, 062121 (2018).
  * [25] T. Demaerel and C. Maes, Active processes in one dimension, Phys. Rev. E. 97, 032604 (2018).
  * [26] I. Santra, U. Basu, S. Sabhapandit, Run-and-tumble particles in two dimensions: Marginal position distributions, Phys. Rev. E. 101, 062120 (2020).
  * [27] D. S. Dean, S. N. Majumdar and H. Schawe, Position distribution in a generalized run-and-tumble process, Phys. Rev. E 103, 012130 (2021).
  * [28] A. Dhar, A. Kundu, S. N. Majumdar, S. Sabhapandit and G. Schehr, Run-and-tumble particle in one-dimensional confining potentials: Steady-state, relaxation, and first-passage properties, Phys. Rev. E. 99, 032132 (2019).
  * [29] U. Basu and S. N. Majumdar, A. Rosso, S. Sabhapandit and G. Schehr, Exact stationary state of a run-and-tumble particle with three internal states in a harmonic trap, J. Phys. A: Math. Theor. 53, 09LT01 (2020).
  * [30] L. Angelani, R. D. Leonardo and M. Paoluzzi, First-passage time of run-and-tumble particles, Eur. Phys. J. E, 37:59 (2014).
  * [31] F. Mori, P. L. Doussal, S. N. Majumdar and G. Schehr, Universal Survival Probability for a $d$-Dimensional Run-and-Tumble Particle Phys. Rev. Lett. 124, 090603 (2020).
  * [32] B. D. Bruyne, S. N. Majumdar and G. Schehr, Survival probability of a run-and-tumble particle in the presence of a drift, J. Stat. Mech. 043211 (2021).
  * [33] P. Singh and A. Kundu, Generalised ‘Arcsine’ laws for run-and-tumble particle in one dimension, J. Stat. Mech., 83205 (2019).
  * [34] P. Singh, S. Sabhapandit and A. Kundu, Run-and-tumble particle in inhomogeneous media in one dimension, J. Stat. Mech. 083207 (2020).
  * [35] L. Angelani, Run-and-tumble particles, Telegrapher’s equation and absorption problems with partially reflecting boundaries, J. Phys. A: Math. Theor. 48, 495003 (2015).
  * [36] P. Singh and A. Kundu, Local time for run and tumble particle, Phys. Rev. E. 103, 042119 (2021).
  * [37] P. Singh, A. Kundu, S. N. Majumdar, H. Schawe, Mean area of the convex hull of a run and tumble particle in two dimensions, J. Phys. A: Math. Theor. 55, 225001 (2022).
  * [38] A. K. Hartmann, S. N. Majumdar and H. Schawe and G. Schehr, The convex hull of the run-and-tumble particle in a plane, J. Stat. Mech., 053401 (2020).
  * [39] T. Banerjee, S. N. Majumdar, A. Rosso and G. Schehr, Current fluctuations in noninteracting run-and-tumble particles in one dimension, Phys. Rev. E. 101, 052101 (2020).
  * [40] E. Woillez, Y. Zhao, Y. Kafri, V. Lecomte and J. Tailleur, Activated Escape of a Self-Propelled Particle from a Metastable State, Phys. Rev. Lett. 122, 258001 (2019).
  * [41] P. Singh, S. Santra and A. Kundu, Extremal statistics of a one-dimensional run and tumble particle with an absorbing wall, J. Phys. A: Math. Theor. 55, 465004 (2022).
  * [42] A. B. Slowman, M. R. Evans and R. A. Blythe, Jamming and attraction of interacting run-andtumble random walkers, Phys. Rev. Lett. 116, 218101 (2016).
  * [43] S. Jose, R. Dandekar and K. Ramola, Current fluctuations in an interacting active lattice gas, J. Stat. Mech. 083208 (2023).
  * [44] T. Agranov, S. Ro, Y. Kafri and V. Lecomte, Macroscopic fluctuation theory and current fluctuations in active lattice gases, SciPost Phys. 14, 045 (2023).
  * [45] T. Agranov, M. E. Cates and R. L. Jack, Entropy production and its large deviations in an active lattice gas, J. Stat. Mech. 123201 (2022).
  * [46] P. Dolai, A. Das, A. Kundu, C. Dasgupta, A. Dhar and K. Vijay Kumar, Universal scaling in active single-file dynamics, SOFT MATTER 16, 7077-7087 (2020).
  * [47] T. Banerjee, R. L. Jack and M. E. Cates, Tracer dynamics in one dimensional gases of active or passive particles, J. Stat. Mech., 013209 (2022).
  * [48] C. Merrigan, K. Ramola, R. Chatterjee, N. Segall, Y. Shokef and B. Chakraborty, Arrested states in persistent active matter: Gelation without attraction, Phys. Rev. R. 2, 013260 (2020).
  * [49] E. Teomy and R. Metzler, Transport in exclusion processes with one-step memory: density dependence and optimal acceleration, J. Phys. A: Math. Theor. 52, 385001 (2019).
  * [50] E. Teomy and R. Metzler, Correlations and transport in exclusion processes with general finite memory, J. Stat. Mech., 103211 (2019).
  * [51] M. Galanti, D. Fanelli and F. Piazza, Persistent random walk with exclusion, Eur. Phys. J. B 86, 456 (2013).
  * [52] S. Put, J. Berx and C. Vanderzande, Non-Gaussian anomalous dynamics in systems of interacting run-and-tumble particles, J. Stat. Mech., 123205 (2019).
  * [53] P. Singh and A. Kundu, Crossover behaviours exhibited by fluctuations and correlations in a chain of active particles, J. Phys. A: Mathematical and Theoretical, 54, 305001 (2021).
  * [54] P. Ferrari and L. Fontes, Fluctuations of a surface submitted to a random average process, Electronic J. Probab. 3, 1 (1998).
  * [55] S. N. Coppersmith, C. h. Liu, S. Majumdar, O. Narayan and T. A. Witten, Model for force fluctuations in bead packs, Phys. Rev. E., 53, 46734685, (1996).
  * [56] R. Rajesh and S. N. Majumdar, Conserved mass models and particle systems in one dimension, J. Stat. Phys. 99, 943 (2000).
  * [57] J. Krug and J. Garcia, Asymmetric particle systems on R, J. Stat. Phys. 99, 31 (2000).
  * [58] Z. A. Melzak, Mathematical ideas, modeling and applications, John Wiley & Sons Inc., New Jersey, U.S., ISBN-13:978-0471593416 2 (1976).
  * [59] S. Ispolatov, P. L. Krapivsky and S. Redner, Wealth distributions in asset exchange models, Eur. Phys. J. B 2, 267 (1998).
  * [60] D. Aldous and P. Diaconis, Hammersley’s interacting particle process and longest increasing subsequences, Probab. Theory Relat. Fields, 103, 199213 (1995).
  * [61] R. Rajesh and S. N. Majumdar, Exact tagged particle correlations in the random average process, Phys. Rev. E. 64, 036103 (2001).
  * [62] G. M. Sch$\ddot{\rm u}$tz, Exact tracer diffusion coefficient in the asymmetric random average pro- cess, J. Stat. Phys. 99(3), 1045 (2000).
  * [63] J. Cividini, A. Kundu, S. N. Majumdar and D. Mukamel, Exact gap statistics for the random average process on a ring with a tracer, J. Phys. A: Math. Theor. 49, 085002 (2016).
  * [64] A. Kundu and J. Cividini, Exact correlations in a single-file system with a driven tracer, EPL 115, 54003 (2016).
  * [65] A. Das, A. Kundu and P. Pradhan, Einstein relation and hydrodynamics of nonequilibrium mass transport processes, Phys. Rev. E. 95, 062128 (2017).
  * [66] A. Poncet, A. Grabsch, P. Illien, and O. Bénichou, Generalized Correlation Profiles in Single-File Systems Phys. Rev. Lett. 127, 220601 (2021).
  * [67] A. Grabsch, P. Rizkallah, A. Poncet, P. Illien, and O. Bénichou, Exact spatial correlations in single-file diffusion, Phys. Rev. E. 107, 044131 (2023).
  * [68] R. Dandekar, S. Chakraborti, and R. Rajesh, Hard core run and tumble particles on a one-dimensional lattice, Phys. Rev. E 102, 062111 (2020).
  * [69] T. Chakraborty and P. Pradhan, Time-dependent properties of run-and-tumble particles. II.: Current fluctuations, Phys. Rev. E 109, 044135 (2024).
  * [70] T. Chakraborty and P. Pradhan, Time-dependent properties of run-and-tumble particles: Density relaxation, Phys. Rev. E 109, 024124 (2024).
  * [71] D. T. Gillespie, A General Method for Numerically Simulating the Stochastic Time Evolution of Coupled Chemical Reactions, J. Comp. Phys. 2, 403-434 (1976).