# Stationary state of harmonic chains driven by boundary resetting

Ritwick Sarkar and Pritam Roy  S. N. Bose National Centre for Basic Sciences,
Kolkata 700106, India
E-mail<EMAIL_ADDRESS>
###### Abstract

We study the nonequilibrium steady state (NESS) of an ordered harmonic chain
of $N$ oscillators connected to two walls which undergo diffusive motion with
stochastic resetting. The intermittent resetting of the walls effectively
emulates two nonequilibrium reservoirs that exert temporally correlated forces
on the boundary oscillators. These reservoirs are characterized by the
diffusion constants and resetting rates of the walls. We find that, for any
finite $N$, the velocity distribution of the bulk oscillators remains non-
Gaussian, as evidenced by a non-zero bulk kurtosis that decays $\sim N^{-1}$.
We calculate the spatio-temporal correlation of the velocity of the
oscillators $\langle v_{l}(t)v_{l^{\prime}}(t^{\prime})\rangle$ both
analytically as well as using numerical simulation. The signature of the
boundary resetting is present at the bulk in terms of the two-time velocity
correlation of a single oscillator and the equal-time spatial velocity
correlation. For the resetting driven chain, the two-time velocity correlation
decay as $t^{-\frac{1}{2}}$ at the large time, and there exists a non-zero
equal-time spatial velocity correlation $\langle
v_{l}(t)v_{l^{\prime}}(t^{\prime})\rangle$ when $l\neq l^{\prime}$. A non-zero
average energy current will flow through the system when the boundary walls
reset to their initial positions at different rates. This average energy
current can be computed exactly in the thermodynamic limit. Numerically we
show that the distribution of the instantaneous energy current at the boundary
is independent of the system size. However, the distribution of the
instantaneous energy current in the bulk approaches a stationary distribution
in the thermodynamic limit.

##  1 Introduction

It is worthwhile to understand the nonequilibrium steady state (NESS) of an
extended system driven by equilibrium reservoirs to comprehend the transport
properties of that system. A paradigmatic model of energy transport in one-
dimensional systems was studied by Rieder, Lebowitz, and Lieb in 1967, which
consists of a harmonic chain of $N$ oscillators, driven by two thermal
reservoirs at the boundaries [1]. For this model, the thermal driving at the
two boundaries leads to a NESS with a nonzero average energy current flowing
through the system. Several generalizations — e.g.; inclusion of anharmonic
interactions between the oscillators of the chain, pinning potential, and
disorders, were studied in the recent past. These studies conclude that these
generalizations lead to some nontrivial behavior of the steady state including
anomalous transport and nonlinear temperature profile [2, 3, 4, 5, 6, 7, 8, 9,
10].

The action of nonequilibrium reservoirs, which do not satisfy fluctuation-
dissipation theorem (FDT)[11], on single probe particles leads to unusual
features like anomalous relaxation dynamics, negative viscosity, modification
of equipartition theorem, etc [14, 15, 16, 17, 18, 19, 20, 21, 22, 12, 13].
Recently, the NESS of a harmonic chain driven by active baths (nonequilibrium
bath) is studied in Ref. [23, 24]. Here we model a nonequilibrium reservoir
using stochastic resetting. Briefly, in our model, the boundary oscillators of
a harmonic chain are attached to two walls that perform free diffusion with
intermittent resetting. These walls act as nonequilibrium reservoirs at the
two ends and the effect of these reservoirs on the NESS of the harmonic chain
is studied in detail in this work.

Stochastic resetting refers to the process where the dynamics of a system
stochastically stops and restarts from some predefined condition[25]. The
paradigmatic example is where the position of a Brownian particle resets to a
fixed point in space with some specific rate [26]. This leads to some striking
features like a nonequilibrium stationary state, finite mean first-passage
time, and anomalous relaxation behavior. The effect of stochastic resetting on
well-known equilibrium and nonequilibrium processes has been observed
theoretically [27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42,
43] as well as experimentally[44, 45, 46]. Stochastic resetting drives a
system out-of-equilibrium irrespective of the underlying dynamics of the
system. Stochastic resetting has found applications in the context of the
search algorithm [47, 48, 49], generation of controlled NESS [50], and quantum
measurements [51, 52]. Moreover, stochastic resetting has been used to model a
nonequilibrium bath [53] which consists of a collection of over-damped
particles that are reset to a specific position.

In this paper, we introduce a simple model of a nonequilibrium reservoir using
stochastic resetting dynamics. The action of tnonequilibrium reservoir is
realized in terms of the wall that performs Brownian motion and stochastically
resets to its initial position at a specific rate. When the boundary
oscillators of a one-dimensional harmonic chain are attached to two of these
reservoirs, the chain attains NESS with an average nonzero energy current
flowing through the system due to the action of these reservoirs. The average
energy current can be easily calculated from the formalism used in Ref. [23]
and the value of the average energy current depends on the resetting rate of
the boundary walls. We find that the velocity distribution of the bulk
oscillator reaches the Gaussian distribution only at a thermodynamically large
system size. The velocity distributions of the boundary oscillators are found
to be size-independent with exponentially decaying tails. Numerical simulation
results confirm that the kurtosis of the velocity of the bulk oscillators
decays as $N^{-1}$. We also calculate the two-point spatio-temporal
correlation of the velocity of the bulk oscillators. We find that at large
time, the two-time velocity correlation $\langle v_{l}(t)v_{l}(0)\rangle$
decay as $t^{-\frac{1}{2}}$ and there exists a non-zero equal-time spatial
velocity correlation $\langle v_{l}(t)v_{l^{\prime}}(t)\rangle$ in the bulk
when $l\neq l^{\prime}$. Instantaneous energy current distribution at the bulk
reaches a stationary distribution in the thermodynamic limit. The exact
expression of this stationary distribution is calculated in Ref. [24]. In
contrast, the distribution of instantaneous energy current at the boundary
remains size-independent. We also discuss the effective thermal limit which is
achieved at a high resetting rate.

The paper is organized as follows; in Sec. 2, we define the model and
summarize our key findings. In Sec. 3, we discuss the velocity distribution
and kurtosis profile of velocity. The two-point spatio-temporal correlation of
the velocity for the resetting driven chain is discussed in Sec. 4. Sec. 5 and
Sec. 6 are devoted to the distribution of instantaneous energy current and the
effective thermal limit. We conclude in Sec. 7 along with some general
remarks.

##  2 Model and Results

Figure 1:  Schematic figure of a linear chain of oscillators connected to two
walls that perform free diffusion and reset to the initial position with
different rates. The boundary oscillators are also attached to two Langevin
baths of temperature $T_{1}$ and $T_{N}$.

We consider a chain of $N$ oscillators, each of mass $m$, coupled to its
nearest neighbors by springs of stiffness constant $k$. Two boundary
oscillators are in contact with Langevin heat baths at temperatures $T_{1}$
and $T_{N}$. These heat baths exert white noises as well as damping forces,
proportional to the velocities, on the boundary oscillators. The boundary
oscillators are also attached to two walls with springs of the same stiffness
constant. These walls perform free diffusion and undergo intermittent
resetting to their initial positions with different rates $r_{1}$ and $r_{N}$.
We denote the displacement of the $l$-th oscillator from its equilibrium
position as $x_{l}$ and the displacements of the left and right resetting
walls as $X_{1}$ and $X_{N}$ respectively. These moving walls exert stochastic
forces on the boundary oscillators due to the intermittent resetting dynamics.
The equations of motion of the oscillators are,

$m\ddot{x_{l}}=\begin{cases}-k(2x_{1}-x_{2})+f_{1}-\gamma\dot{x}_{1}+\xi_{1},&\text{when
$l=1$, }\\\ -k(2x_{l}-x_{l-1}-x_{l+1}),&\text{$\forall l\in[2,N-1]$, }\\\
-k(2x_{N}-x_{N-1})+f_{N}-\gamma\dot{x}_{N}+\xi_{N}.&\text{when
$l=N$}.\end{cases}$ (1)

$f_{1}=kX_{1}$ and $f_{N}=kX_{N}$ are the stochastic forces exerted by the
boundary walls. $\xi_{1}$ and $\xi_{N}$ are white noise from the Langevin heat
bath acting on the boundary oscillators and satisfy the fluctuation-
dissipation relation [11],

$\langle\xi_{i}(t)\xi_{j}(t^{\prime})\rangle=2\gamma
T_{i}\delta_{ij}\delta(t-t^{\prime})\quad\text{and}~{}i,j=1,N.$ (2)

$T_{1}$ and $T_{N}$ are the temperatures of the thermal reservoir at the two
ends of the chain. $\gamma$ is the damping coefficient of the respective
reservoir. Note that, we assume the motion of the wall is independent of the
boundary oscillators.

The boundary walls perform free diffusion with Poisson resetting [25, 26]. In
a small time interval $\Delta t$, the position of the $i$-th wall, relative to
its equilibrium position, is updated by the following procedure[25],

$X_{i}(t+\Delta t)=\begin{cases}0&\text{with probability $r_{i}\Delta t$, }\\\
X_{i}(t)+\sqrt{2D\Delta t}\,\eta_{i}(t)&\text{with probability $1-r_{i}\Delta
t$, and}~{}i=1,N,\\\ \end{cases}$ (3)

where $D$ is the diffusion constant which we assumed to be the same for both
the boundary walls and $\eta_{i}(t)$ has mean and second moment,
$\langle\eta_{i}(t)\rangle=0$, and
$\langle\eta_{i}(t)\eta_{j}(t)\rangle=\delta_{ij}$. The auto-correlation of
the position $X_{i}(t)$ is given by[54],

$\langle
X_{i}(t)X_{j}(t^{\prime})\rangle=\delta_{ij}\frac{2D}{r_{i}}e^{-r_{i}|t-t^{\prime}|}\Big{(}1-e^{-r_{i}~{}\text{min}(t,t^{\prime})}\Big{)}$
(4)

A stationary state is reached when $t\to\infty,~{}t^{\prime}\to\infty$ and
$|t-t^{\prime}|$ is finite. In the stationary state, the auto-correlation of
the position $X_{i}(t)$ becomes,

$\langle
X_{i}(t)X_{j}(t^{\prime})\rangle=\delta_{ij}\frac{2D}{r_{i}}e^{-r_{i}|t-t^{\prime}|}.$
(5)

The forces experienced by the boundary oscillators at the two boundaries are
$f_{1}=kX_{1}$ and $f_{N}=kX_{N}$ which are stochastic hin nature. These
stochastic forces are correlated as,

$\langle
f_{i}(t)f_{j}(t^{\prime})\rangle=\delta_{ij}\,a_{i}^{2}\,e^{-r_{i}|t-t^{\prime}|}\quad\text{with}\quad
a_{i}^{2}=\frac{2Dk^{2}}{r_{i}}\quad i=1,N.$ (6)

Here $a_{i}$ is the strength of the noise coming from the $i$-th resetting
wall. The stochastic force $f_{i}(t)$ acting on the boundary oscillator
violates the fluctuation-dissipation relation. For a resetting driven linear
chain, the equations of motion Eq. (1) can be solved in the frequency domain,

$x_{l}(t)=\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}e^{-i\omega
t}[G_{l1}(\omega)\tilde{f}_{1}(\omega)+G_{lN}(\omega)\tilde{f}_{N}(\omega)].$
(7)

Here $G(\omega)$ denotes the Green’s function matrix [See A for a detailed
derivation] and $\tilde{f}_{i}(\omega)$ is the Fourier transform of $f_{i}(t)$
with respect to $t$. The two-point correlation of the stochastic force
$\tilde{f}_{i}(t)$ in the Fourier domain and can be written using Eq. (6) as,

$\langle\tilde{f}_{i}(\omega)\tilde{f}_{j}(\omega^{\prime})\rangle=2\pi\delta_{ij}\tilde{g}(\omega,r_{i})\delta(\omega+\omega^{\prime})\quad\text{and}\quad\tilde{g}(\omega,r_{i})=\frac{2a_{i}^{2}r_{i}}{r_{i}^{2}+\omega^{2}},$
(8)

where $\tilde{g}(\omega,r_{i})$ is the frequency spectrum of the resetting
force.

It can be easily shown that for the harmonic chain, the average energy current
has two independent components— a thermal component, $J_{\text{therm}}$ and a
resetting component $J_{\text{r}}$. $J_{\text{therm}}$ is proportional to the
temperature difference at the two ends of the chain, $T_{1}-T_{N}$, and
$J_{\text{r}}$ will depend on the resetting rates of the boundary walls. The
main objective of this paper is to characterize the NESS driven by the
boundary resetting and we choose $T_{1}=T_{N}=0$ for the rest of the
calculation. In this paper the numerical simulations are done using the scheme
of [55] which are accurate to order $(\Delta t)^{2}$.

Before going into the detailed discussion, a summary of the result is
presented.

  * •

Velocity fluctuation: We measure the velocity distribution of the oscillators
$P(v_{l})$ in the NESS. We find that, $P(v_{l})$ of the bulk oscillators
approach Gaussian distribution in the thermodynamic limit, $N\to\infty$ with
width $\hat{T}_{\text{bulk}}=m\langle\dot{x}_{l}^{2}\rangle$, which is a
measure of the kinetic temperature. To characterize the finite-size
dependence, we measure the kurtosis profile of velocity, $\kappa_{l}=\langle
v_{l}^{4}\rangle/\langle v_{l}^{2}\rangle^{2}-3$. We find that for the bulk
oscillators, $\kappa_{l}$ decays as $N^{-1}$. The velocity distributions of
the boundary oscillators [$P(v_{1})$ or $P(v_{N})$] are found to be size-
independent as well as non-Gaussian. The NESS of the thermally driven chain is
Gaussian, independent of the size of the harmonic chain and this is a sharp
difference between the resetting driven and thermally driven case.

  * •

Spatio-temporal correlation of velocity: Another physical observable of
interest is the spatio-temporal velocity correlation of the oscillators
$\langle v_{l}(t)v_{l^{\prime}}(t^{\prime})\rangle$, in the steady state. We
calculate two-time velocity correlation of single oscillator ($\langle
v_{l}(t)v_{l}(t^{\prime})\rangle$) and equal-time spatial velocity correlation
($\langle v_{l}(t)v_{l^{\prime}}(t)\rangle$) for resetting driven chain
explicitly. We show that in the thermodynamic limit, the two-time velocity
correlation of a single oscillator in the bulk is,

$\langle v_{l}(t)v_{l}(t^{\prime})\rangle=\frac{1}{\gamma
m}\int_{0}^{\pi}\frac{dq}{2\pi}\cos\Big{[}\omega_{c}\sin{\Big{(}\frac{q}{2}\Big{)}}(t-t^{\prime})\Big{]}\Bigg{[}\frac{a_{1}^{2}r_{1}}{r_{1}^{2}+\omega_{c}^{2}\sin^{2}{\frac{q}{2}}}+\frac{a_{N}^{2}r_{N}}{r_{N}^{2}+\omega_{c}^{2}\sin^{2}{\frac{q}{2}}}\Bigg{]}.$
(9)

where $\omega_{c}=2\sqrt{k/m}$. When $t-t^{\prime}\gg\omega_{c}^{-1}$,
$\langle v_{l}(t)v_{l}(t^{\prime})\rangle\propto
J_{0}\big{(}\omega_{c}(t-t^{\prime})\big{)}$, here $J_{0}(z)$ is the Bessel
function of first kind [56]. The equal-time spatial velocity correlation in
thermodynamic limit is given by,

$\langle
v_{l}(t)v_{l^{\prime}}(t)\rangle=\frac{1}{2\gamma}\Big{[}a_{1}^{2}\theta(r_{1},l,l^{\prime})+a_{N}^{2}\theta(r_{N},l,l^{\prime})\Big{]}.$
(10)

where,

$\theta(r_{i},l,l^{\prime})=\frac{r_{i}}{mr_{i}^{2}+4k}\,{}_{3}\tilde{F}_{2}\,\Big{[}1/2,1,1;1-l^{\prime}+l,1+l^{\prime}-l;\frac{4k}{4k+mr_{i}^{2}}\Big{]}.$
(11)

Here $l^{\prime}-l$ is an integer and ${}_{p}\tilde{F}_{q}$ is a generalized
regularized hypergeometric function [56].

  * •

Fluctuation of energy current: We also measure the distribution of
instantaneous energy current flowing through the system in the steady state.
For the resetting driven oscillator chain, the distribution of instantaneous
energy current at the bulk, $\mathcal{J}_{l}$ is found to approach

$P(\mathcal{J}_{l})=\frac{1}{\pi\sqrt{g_{l}}}e^{\frac{J_{\text{r}}}{g_{l}}\mathcal{J}_{l}}K_{0}\Big{(}\frac{u_{l}}{g_{l}}|\mathcal{J}_{l}|\Big{)}.$
(12)

in the thermodynamic limit. Here $K_{0}(z)$ is the zeroth-order modified
Bessel function of the second kind and the definition of $J_{\text{r}}$,
$g_{l}$ and $u_{l}$ is given in Eq. (34) and Eq. (37). The numerical
simulation also confirms that the distribution of instantaneous energy current
at the boundary $\mathcal{J}_{1}$, is size-independent.

  * •

Effective thermal limit: At high resetting rate, the $i$-th boundary wall acts
as effective Langevin bath with effective temperature $T_{i}^{\text{eff}}$. In
this limit, the known result of the average energy current for a thermally
driven oscillator chain is recovered [1]. The equal-time spatial velocity
correlation also becomes uncorrelated for $l\neq l^{\prime}$ in this effective
thermal limit.

##  3 Velocity distribution and kurtosis profile of the velocity

The probability distributions of the velocity of the constituent oscillators
play an important role in characterizing the NESS of the driven oscillator
chain. For thermally or activity-driven oscillator chains, the velocity
distributions of the oscillators at the bulk are found to be Gaussian with its
width given by local kinetic temperature $\hat{T}_{\text{bulk}}=m\langle
v_{l}^{2}\rangle$ [1, 24]. However, the velocity distribution of boundary
oscillators depends on the specific dynamics of driving. For resetting driven
harmonic chain also, we first measure the velocity distributions of bulk and
boundary oscillators for different system sizes and different resetting rates
at the two boundaries.

###  3.1 Distribution of velocity

The velocity of the $l$-th oscillator for resetting driven chain is [from Eq.
(7)],

$v_{l}(t)=\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}e^{-i\omega
t}(-i\omega)\big{[}G_{l1}(\omega)\tilde{f}_{1}(\omega)+G_{lN}(\omega)\tilde{f}_{N}(\omega)\big{]},$
(13)

and we numerically measure the velocity distribution of the bulk oscillators
$P(v_{l})$ for different system sizes. In Fig. 2(a), the numerically measured
$P(v_{l})$ is plotted for different system sizes. For any finite $N$,
$P(v_{l})$ shows significant deviation at the tails from the Gaussian
distribution with standard deviation $\hat{T}_{\text{bulk}}$. However, the
deviation tends to decrease with an increase in system size. Therefore we
conclude that,

$\text{For}~{}N\to\infty,~{}~{}P(v_{l})=\frac{1}{\sqrt{2\pi\hat{T}_{bulk}}}\exp{\Bigg{(}\frac{-mv_{l}^{2}}{2\hat{T}_{\text{bulk}}}\Bigg{)}},$
(14)

where the value of kinetic temperature at the bulk is [see B.1 for the
derivation],

$\hat{T}_{\text{bulk}}=\frac{a_{1}^{2}}{2\gamma\sqrt{r_{1}^{2}+\frac{4k}{m}}}+\frac{a_{N}^{2}}{2\gamma\sqrt{r_{N}^{2}+\frac{4k}{m}}}.$
(15)

This finite-size dependence of velocity distribution makes the resetting
driven harmonic chain different than the thermally driven harmonic chain,
where the steady state is Gaussian irrespective of the size. Numerical
simulations with different system sizes suggest, unlike velocity distributions
of the bulk oscillators, $P(v_{l})$, the velocity distribution of the boundary
oscillator $P(v_{1})$ [or $P(v_{N})$] is size-independent, [see Fig. 3].

The scaled velocity distribution of the oscillator at the left boundary is
shown in Fig. 3. Here the numerically measured velocity distribution
$P(v_{1})$ is scaled with numerically measured $\sqrt{\langle
v_{1}^{2}\rangle}$, which shows a scaling collapse. The numerical simulation
also indicates that the velocity distribution of the boundary oscillator has
an exponentially decaying tail [which is shown with the red dashed line in
Fig. 3].

Figure 2:  (a) Velocity distribution of the bulk oscillator, $l=N/2$ measured
from simulation with chains of different sizes for the same resetting rate,
$r_{1}=r_{N}=2.5$. The black solid line corresponds to the Gaussian
distribution for thermodynamically large system size ($N\to\infty$) given by
Eq. (14). (b) Velocity distribution of the boundary oscillator, $l=1$ measured
from simulation with chains of different sizes for the same resetting rate,
$r_{1}=r_{N}=2.5$. Other parameters are $D=m=\gamma=k=1$. Figure 3:  The plot
of the scaled velocity distribution of the oscillator at the left boundary for
fixed $r_{N}=2.5$ and three different values of $r_{1}$. The symbols denote
the numerical simulation with a chain of $N=1024$. Other parameters are
$D=m=\gamma=k=1$.

###  3.2 Kurtosis profile of velocity

To characterize the finite-size dependence of the velocity distribution of the
bulk oscillators, we measure the kurtosis profile of velocity. The definition
of the kurtosis of the velocity of the $l$-th oscillator is,

$\kappa_{l}=\frac{\langle v_{l}^{4}\rangle}{\langle v_{l}^{2}\rangle^{2}}-3.$
(16)

If the velocity distribution of the $l$-th oscillator is Gaussian, then the
$\kappa_{l}=0$. Typically the distribution will have a fatter tail when
$\kappa_{l}>0$. We compute the velocity kurtosis profile for different system
sizes using numerical simulation. In Fig. 4(a), $\kappa_{l}$ for different
system sizes is plotted for $r_{1}=r_{N}$ and from this plot, it is clear that
the value of the kurtosis at the bulk decreases as the system size increases.
In Fig. 4(b) $\kappa_{l}$ for a chain of length $N=1024$ is plotted for a
fixed value of $r_{1}$ and three different values of $r_{N}$. Note that, the
velocity kurtosis profile is symmetric for $r_{1}=r_{N}$ and the velocity
kurtosis profile is symmetric in $|r_{1}-r_{N}|$. It is clear from Fig. 4(b)
that the velocity distribution of the bulk oscillators will reach the Gaussian
distribution if the resetting rate is high. However, the boundary oscillators
will have a fatter tail even at a higher resetting rate.

Figure 4: (a) Numerically measured velocity kurtosis profile for the same
resetting rate at the two ends, $r_{1}=r_{N}=2.5$ with chains of different
sizes. (b) Velocity kurtosis profile for a chain of $N=1024$ for fixed
$r_{1}=2.5$ and different $r_{N}$, other parameters are $D=m=\gamma=k=1$.
Figure 5: Numerically measured velocity kurtosis of the $l$-th oscillator at
the bulk (where $l=N/2$) as a function of system size $N$ for fixed
$r_{1}=2.5$ and different value of $r_{N}$. Other parameters are
$D=m=\gamma=k=1$ and the red dashed line shows $\sim N^{-1}$ decay.

It is intriguing to investigate the system size dependence of the velocity
kurtosis profile, therefore we numerically estimate $\kappa_{\text{bulk}}$
[$\kappa_{l}$ at $l=N/2$] with increasing system size $N$ for fixed $r_{1}$
and three different values of $r_{N}$. The result is shown in Fig. 5. From the
result of numerical simulation, we conclude that,

$\kappa_{\text{bulk}}\sim N^{-1}.$ (17)

Therefore only in the thermodynamic limit, the velocity distributions of the
bulk oscillators become Gaussian. Size-dependent NESS is the signature of the
resetting driven harmonic chain because such behavior is absent for the
thermally driven case.

##  4 Spatio-temporal correlation of velocity

Correlation functions play an important role in characterizing transport
processes. For example, the integral of equilibrium velocity-velocity
correlation $\langle v(t)v(t^{\prime})\rangle$ is related to the diffusion
coefficient of a Brownian particle. For thermally driven transport in the
mass-disordered harmonic chain, the two-time correlation of the current
determines the asymptotic size dependence of the current fluctuation in the
steady state [57]. In this section, we discuss the two-point spatio-temporal
correlation of velocity $\langle v_{l}(t)v_{l^{\prime}}(t^{\prime})\rangle$ in
the bulk.

Using Eq. (13), we can write the two-point spatio-temporal correlation of
velocity as,

$\langle
v_{l}(t)v_{l^{\prime}}(t^{\prime})\rangle=\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\omega^{2}e^{-i\omega(t-t^{\prime})}\big{[}G_{l1}(\omega)G_{l^{\prime}1}^{*}(\omega)\tilde{g}(\omega,r_{1})+G_{lN}(\omega)G_{l^{\prime}N}^{*}(\omega)\tilde{g}(\omega,r_{N})\big{]}.$
(18)

In the following, we calculate the two-time velocity correlation of a single
oscillator $\langle v_{l}(t)v_{l}(t^{\prime})\rangle$ and equal-time spatial
velocity correlation $\langle v_{l}(t)v_{l^{\prime}}(t)\rangle$ in the
stationary state explicitly for the resetting driven chain.

###  4.1 Two-time velocity correlation of single oscillator

In the stationary state, the two-time velocity correlation of a single
oscillator, $\langle v_{l}(t)v_{l}(t^{\prime})\rangle$ is a function of
$t-t^{\prime}$ [see Eq. (18)]. To simplify the calculation, we consider
$\langle v_{l}(t)v_{l}(0)\rangle$ which can be written as,

$\langle
v_{l}(t)v_{l}(0)\rangle=\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\omega^{2}e^{-i\omega
t}\big{[}|G_{l1}(\omega)|^{2}\tilde{g}(\omega,r_{1})+|G_{lN}(\omega)|^{2}\tilde{g}(\omega,r_{N})\big{]}.$
(19)

Figure 6: (a) Plot of $\langle v_{l}(t)v_{l}(0)\rangle$ as a function of time
for fixed $r_{1}=1.0$ and three different value of $r_{N}$. Symbols correspond
to the data obtained from numerical simulation with a chain of $N=256$. The
black solid line corresponds to Eq. (22). (b) Plot of $\langle
v_{l}(t)v_{l}(0)\rangle$ in the large time limit is shown along with the red
dashed line which corresponds to $1/\sqrt{t}$ behavior. The other parameters
are $D=m=\gamma=k=1$.

The above equation has two separate contributions coming from two resetting
walls at the two boundaries and we can write,

$\langle v_{l}(t)v_{l}(0)\rangle=a_{1}^{2}{C}(r_{1},t)+a_{N}^{2}{C}(r_{N},t).$
(20)

where ${C}(r_{i},t)$ is the contribution from the $i$-th resetting wall. For
the oscillators in the bulk [$1\ll l\ll N$], $C(r_{i},t)$ can be computed
explicitly and takes a simple form in the thermodynamic limit,

${C}(r_{i},t)=\frac{r_{i}}{\gamma
m}\int_{0}^{\pi}\frac{dq}{2\pi}\frac{\cos\Big{[}\omega_{c}t\sin{\Big{(}\frac{q}{2}\Big{)}}\Big{]}}{r_{i}^{2}+\omega_{c}^{2}\sin^{2}{\frac{q}{2}}},$
(21)

see B for the detailed derivation. Combining the contributions from both the
reservoirs, we can write $\langle v_{l}(t)v_{l}(0)\rangle$ in the
thermodynamic limit as,

$\langle v_{l}(t)v_{l}(0)\rangle=\frac{1}{\gamma
m}\int_{0}^{\pi}\frac{dq}{2\pi}\cos\Big{[}\omega_{c}t\sin{\Big{(}\frac{q}{2}\Big{)}}\Big{]}\Bigg{[}\frac{a_{1}^{2}r_{1}}{r_{1}^{2}+\omega_{c}^{2}\sin^{2}{\frac{q}{2}}}+\frac{a_{N}^{2}r_{N}}{r_{N}^{2}+\omega_{c}^{2}\sin^{2}{\frac{q}{2}}}\Bigg{]}.$
(22)

The above integral can be evaluated numerically for arbitrary values of $t$.
This is illustrated in Fig 6(a) where the symbols correspond to the data
obtained from numerical simulation for a fixed $r_{1}$ and different values of
$r_{N}$. The black solid lines in Fig 6(a) correspond to the numerical
integration of Eq. (22).

It is difficult to get the closed form of Eq. (22). However, it is possible to
predict the small and large time behavior of $\langle
v_{l}(t)v_{l}(0)\rangle$. For very small value of $t$ or
$t\ll\omega^{-1}_{c}$, $\cos[\omega_{c}t\sin{(q/2)}]$ can be expanded up to
second order in time $t$. Therefore, for small value of $t$, $\langle
v_{l}(t)v_{l}(0)$ is proportional to $t^{2}$. Similarly, when
$t\gg\omega_{c}^{-1}$, the integral can be performed exactly [see B.1 for the
detailed derivation]. In this large time limit,

$\langle
v_{l}(t)v_{l}(0)\rangle\simeq\Big{[}\frac{a_{1}^{2}r_{1}}{2\gamma(mr_{1}^{2}+4k)}+\frac{a_{N}^{2}r_{N}}{2\gamma(mr_{N}^{2}+4k)}\Big{]}J_{0}(\omega_{c}t),$
(23)

where $\omega_{c}=2\sqrt{k/m}$, and $J_{0}(z)$ is the Bessel function of the
first kind. In Fig 6(b) numerically measured $\langle v_{l}(t)v_{l}(0)\rangle$
($l=N/2$) is plotted for large $t$, which is an oscillatory function of $t$
and the envelope of $\langle v_{l}(t)v_{l}(0)\rangle$ seemed to decay as
$t^{-\frac{1}{2}}$, this is similar to the behavior of $J_{0}(\omega_{c}t)$
when $t\gg\omega_{c}^{-1}$.

###  4.2 Equal-time spatial velocity correlation

Figure 7: Plot of $\langle v_{N/2}(t)v_{l^{\prime}}(t)\rangle$ as a function
of $N/2-l^{\prime}$ for fixed $r_{1}=2.0$ and three different value of
$r_{N}$. Symbols correspond to the data obtained from numerical simulation
with a chain of size $N=64$. The black dashed lines correspond to Eq. (28)
with $l-l^{\prime}$ as an integer. The other parameters are the same as the
previously mentioned simulations.

Equal-time spatial velocity correlation, $\langle
v_{l}(t)v_{l^{\prime}}(t)\rangle$, can be calculated by taking $t=t^{\prime}$
in Eq. (18).

$\langle
v_{l}(t)v_{l^{\prime}}(t)\rangle=\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\omega^{2}[G_{l1}(\omega)G_{l^{\prime}1}^{*}(\omega)\tilde{g}(\omega,r_{1})+G_{lN}(\omega)G_{l^{\prime}N}^{*}(\omega)\tilde{g}(\omega,r_{N})].$
(24)

In the stationary state $\langle v_{l}(t)v_{l^{\prime}}(t)\rangle$ does not
depend on the time and has two separate contributions coming from the two
walls performing intermittent resetting. Therefore $\langle
v_{l}(t)v_{l^{\prime}}(t)\rangle$ can be written as,

$\langle
v_{l}(t)v_{l^{\prime}}(t)\rangle=a_{1}^{2}\theta(r_{1},l,l^{\prime})+a_{N}^{2}\theta(r_{N},l,l^{\prime}).$
(25)

In the bulk, i.e. $1\ll l,l^{\prime}\ll N$, the contribution from the $i$-th
wall in the limit $N\rightarrow\infty$ is [see B for detailed derivation],

$\theta(r_{i},l,l^{\prime})=\frac{r_{i}^{2}}{2\pi\gamma}\int_{0}^{\pi}dq\frac{\cos{(l^{\prime}q-lq)}}{mr_{i}^{2}+4k\sin^{2}\frac{q}{2}}$
(26)

where,

$\theta(r_{i},l,l^{\prime})=\frac{r_{i}}{mr_{i}^{2}+4k}\,{}_{3}\tilde{F}_{2}\,\Big{[}1/2,1,1;1-l^{\prime}+l,1+l^{\prime}-l;\frac{4k}{4k+mr_{i}^{2}}\Big{]}~{}\text{and}~{}l^{\prime}-l\in\mathbb{I}.$
(27)

Here ${}_{p}\tilde{F}_{q}$ is a generalized regularized hypergeometric
function [56]. Combining the contribution of both reservoirs we arrive at the,

$\langle
v_{l}(t)v_{l^{\prime}}(t)\rangle=\frac{a_{1}^{2}}{2\gamma}\theta(r_{1},l,l^{\prime})+\frac{a_{N}^{2}}{2\gamma}\theta(r_{N},l,l^{\prime}).$
(28)

In Fig 7 we plotted numerically measured $\langle
v_{l}(t)v_{l^{\prime}}(t)\rangle$ for a fixed value of $r_{1}$ and different
values of $r_{N}$ which matches with the analytic prediction, Eq. (28) well.

The two time velocity correlation of single oscillator $\langle
v_{l}(t)v_{l}(0)\rangle$ and the equal-time spatial velocity correlation
$\langle v_{l}(t)v_{l^{\prime}}(t)\rangle$ in the stationary state depends on
the frequency spectrum of the nonequilibrium reservoir
$\tilde{g}(\omega,r_{i})$ [Eq. (18)]. Therefore, these two observables can be
evaluated for the nonequilibrium driving described in Ref. [23, 24].

##  5 Fluctuation and second moment of instantaneous energy current

The definition of instantaneous energy current from the left reservoir to the
system is,

$\mathcal{J}_{1}=(-\gamma v_{1}+f_{1})v_{1}.$ (29)

Figure 8: Plot of average energy current $J_{\text{r}}$ as functions of
resetting rate $r_{1}$ for different values of $r_{N}$. Symbols obtained from
the numerical simulation with $N=1024$ oscillators and other parameters are
$D=m=\gamma=k=1$. Black solid lines correspond to Eq. (34).

Similarly instantaneous energy current from the right reservoir to the system
is

$\mathcal{J}_{N+1}=(-\gamma v_{N}+f_{N})v_{N}.$ (30)

In the bulk, the expression for the instantaneous energy current between
$l-1$th and $l$-th oscillator,

$\mathcal{J}_{l}=\frac{k}{2}\left(v_{l-1}+v_{l}\right)\left(x_{l-1}-x_{l}\right).$
(31)

The steady state of the resetting-driven harmonic chain is characterized by
the existence of non-zero average energy flux through the chain. In NESS, the
average energy current flowing through the chain is,

$J_{\text{r}}=\langle\mathcal{J}_{1}\rangle=\langle\mathcal{J}_{2}\rangle=\cdots=\langle\mathcal{J}_{l}\rangle=\cdots=-\langle\mathcal{J}_{N+1}\rangle.$
(32)

The analytical expression of the average energy current $J_{\text{r}}$ can be
calculated easily following the method described in Ref. [23]. For
completeness, we quote the result here.

$\displaystyle J_{\text{r}}$ $\displaystyle=$
$\displaystyle\frac{m}{4\gamma^{2}}\Big{[}a_{1}^{2}{\zeta}(r_{1})-a^{2}_{N}{\zeta}(r_{N})\Big{]}\quad\text{where}$
(33) $\displaystyle{\zeta}(r_{j})$ $\displaystyle=$
$\displaystyle\frac{r_{j}k^{2}\left(\sqrt{1+\frac{4\gamma^{2}}{mk}}-1\right)+r_{j}^{3}\gamma^{2}\left(1-\sqrt{1+\frac{4k}{mr_{j}^{2}}}\right)}{(k^{2}-r_{j}^{2}\gamma^{2})}.$
(34)

In Fig.8, we have shown the plot of average energy current $J_{\text{r}}$ for
fixed $r_{1}$ and different $r_{N}$ along with the analytical prediction Eq.
(34). In the next subsection, we will discuss the distribution of the
instantaneous energy current at the bulk and the boundary.

###  5.1 Distribution of instantaneous energy current at the bulk

In this section, we will concentrate on the stationary state distribution of
instantaneous current at the bulk $P(\mathcal{J}_{l})$. It is clear from the
definition of $\mathcal{J}_{l}$ mentioned in Eq. (31) that the
$P(\mathcal{J}_{l})$ will depend on the joint distribution
$\mathcal{P}(\\{x_{l},v_{l}\\})$. In Sec. 3, we have shown that the velocity
distribution of the bulk oscillators of the resetting driven harmonic chain
approaches the Gaussian distribution Eq. (14) for thermodynamically large
system size. A similar result is also expected for the position variables
$x_{l}$ in the bulk. Therefore, we assume that the joint probability
distribution of $\\{x_{l},v_{l}\\}$ of the bulk oscillators,
$\mathcal{P}(\\{x_{l},v_{l}\\})$ is a multivariate Gaussian in the
thermodynamic limit. Therefore, we can write the joint distribution of
$\\{v_{l-1},v_{l},x_{l-1},x_{l}\\}$ in that limit as,

$\mathcal{P}(v_{l-1},v_{l},x_{l-1},x_{l})=\frac{\exp\left[-\frac{1}{2}{W}_{l}^{T}\mathcal{E}^{-1}_{l}{W}_{l}\right]}{\sqrt{(2\pi)^{4}\text{det}(\mathcal{E}_{l})}},$
(35)

where, $W_{l}^{T}=(v_{l-1}~{}v_{l}~{}x_{l-1}~{}x_{l})$ and
$\mathcal{E}_{l}=\langle W_{l}W_{l}^{T}\rangle$, the corresponding $4\times 4$
correlation matrix. For Gaussian NESS, the instantaneous current distribution
$P(\mathcal{J}_{l})$ at the bulk is derived in Ref. [24]. Here we quote the
final result,

$P(\mathcal{J}_{l})=\frac{1}{\pi\sqrt{g_{l}}}\exp\Big{(}{\frac{J_{\text{r}}}{g_{l}}\mathcal{J}_{l}\Big{)}}K_{0}\Big{(}\frac{u_{l}}{g_{l}}|\mathcal{J}_{l}|\Big{)}.$
(36)

Here $K_{0}(z)$ is the zeroth-order modified Bessel function of the second
kind, and $J_{\text{r}}$ is the average energy current flowing through the
resetting driven harmonic chain. Note that, the specific form of distribution
function mentioned in Eq. (36) also appeared in the context of time-integrated
heat current fluctuations of Brownian particles in an active environment [58].
For resetting driven chain, the $g_{l}$ and $u_{l}$ appeared in Eq. (36) are
given by,

$g_{l}=\frac{k\hat{T}_{\text{bulk}}}{2}[\hat{T}_{\text{bulk}}+\langle
v_{l-1}v_{l}\rangle]-J_{\text{r}}^{2},~{}~{}\text{and}~{}~{}u_{l}=\sqrt{g_{l}+J_{\text{r}}^{2}}.$
(37)

Here $\langle v_{l-1}v_{l}\rangle$ is equal time spatial velocity correlation
for a driven harmonic chain.

Figure 9: (a) Numerically measured distribution of instantaneous current
distribution at bulk $P(\mathcal{J}_{l})$ for different resetting rate at
boundary, $r_{1}=2.5$ and $r_{N}=3.5$. In (b), the behavior of
$P(\mathcal{J}_{l})$ near the origin is shown. The black solid line in both
(a) and (b) corresponds to Eq. (36) and Eq.(38) respectively. Other parameters
are set at $D=m=\gamma=k=1$.

In Fig. 9(a), numerically measured $P(\mathcal{J}_{l})$ of different system
sizes are plotted for different values of $r_{1}$ and $r_{N}$ and for a large
system size (e.g. $N=1024$), there is a good agreement with Eq. (36), shown by
the black solid line. This serves as a cross-check of our assumption that in
the thermodynamic limit, the joint distribution
$\mathcal{P}(\\{v_{l},x_{l}\\})$ for the resetting driven harmonic chain, is a
multivariate Gaussian, where $1\ll l\ll N$. It is noteworthy that such strong
finite-size dependence of $P(\mathcal{J}_{l})$ is not reported for activity-
driven harmonic chain [24]. This finite-size dependence is absent for
thermally driven harmonic chains as the NESS is Gaussian irrespective of the
size of the system. For small system sizes, deviations from the analytic curve
are visible on the tails of the distributions which decreases with the
increase of the system size.

Using the asymptotic behavior of $K_{0}(z)$ for $z\to 0$, we get,

$P(\mathcal{J}_{l})=-\frac{1}{\pi\sqrt{g_{l}}}\left(\ln{\frac{u_{l}}{2g_{l}}|\mathcal{J}_{l}|}+E_{\gamma}\right)+O(\mathcal{J}_{l}),$
(38)

near $\mathcal{J}_{l}=0$. Here $E_{\gamma}\simeq 0.577216$ is the Euler’s
constant [56]. From this figure, it is clear that $P(\mathcal{J}_{l})$ has a
logarithmic divergence even for small system sizes and approaches Eq. (38)
with increasing system size.

###  5.2 Distribution of the instantaneous energy current at the boundary

The definition of instantaneous current from the left reservoir to the system
is,

$\mathcal{J}_{1}=(-\gamma v_{1}+f_{1})v_{1},$ (39)

and corresponding distribution can be written as,

$P(\mathcal{J}_{1})=\int dv_{1}df_{1}\,\delta[\mathcal{J}_{1}-(-\gamma
v_{1}+f_{1})v_{1}]\,\mathcal{P}(v_{1},f_{1}).$ (40)

We mentioned in Sec. 3.1 that the velocity distributions of the boundary
oscillators [$P(v_{1})$ or $P(v_{N})$] are independent of the size of the
system. It is important to note that, in the stationary state, the position
distribution $P(X_{i})$ of the boundary wall [whose dynamics is given in Eq.
(3)] is given by [26],

$P(X_{i})=\sqrt{\frac{r_{i}}{4D}}\exp{\Big{(}-\sqrt{\frac{r_{i}}{D}}|X_{i}|\Big{)}}.$
(41)

Therefore the force $f_{1}=kX_{1}$ and $f_{N}=kX_{N}$ exerted on the boundary
oscillators are of the same form and independent of the size of the oscillator
chain.

Figure 10: (a) Numerically measured distribution of instantaneous current
distribution at left boundary $P(\mathcal{J}_{1})$ for $r_{1}=2.5$ and
$r_{N}=3.5$. (b) Behavior of $P(\mathcal{J}_{1})$ near the origin. The other
parameters are set at $D=m=\gamma=k=1$.

As the marginal distributions $P(v_{1})$ and $P(f_{1})$ are both size-
independent, it can be inferred that the joint distribution
$\mathcal{P}(v_{1},f_{1})$ is size-independent. Therefore, the distribution of
instantaneous current at the boundary will also be size-independent [see Eq.
(40)]. The analytical form of $P(\mathcal{J}_{1})$ is hard to calculate as the
analytical form of $P(v_{1})$ is unknown. However, in Fig. 10(a), the
numerical simulation result for $P(\mathcal{J}_{1})$ with different system
size is shown, which confirms that the $P(\mathcal{J}_{1})$ is size-
independent. Behavior of $P(\mathcal{J}_{1})$ near origin is shown in Fig.
10(b). Similar to $P(\mathcal{J}_{l})$, $P(\mathcal{J}_{1})$ also has
logarithmic divergence near the origin.

###  5.3 Second moment of the energy current

The second moment of the instantaneous energy current at bulk
$\mathcal{J}_{l}$ and boundary $\mathcal{J}_{1}$ can be calculated from the
moment generating functions given by,

$\langle\mathcal{J}_{l}^{2}\rangle=-\frac{d^{2}}{d\mu^{2}}\langle
e^{i\mu\mathcal{J}_{l}}\rangle\Big{|}_{\mu=0},$ (42)
$\langle\mathcal{J}_{1}^{2}\rangle=-\frac{d^{2}}{d\mu^{2}}\langle
e^{i\mu\mathcal{J}_{1}}\rangle\Big{|}_{\mu=0}.$ (43)

Figure 11: (a) Plot of second moment of $\mathcal{J}_{l}$ as functions of
resetting rate $r_{1}$. The black solid line corresponds to the Eq. (44). (b)
The plot of the second moment of $\mathcal{J}_{1}$ as functions of resetting
rate $r_{1}$. Symbols obtained from the numerical simulation with $N=1024$
oscillator and for all numerical simulations other parameters are
$D=m=\gamma=k=1$.

In Sec. 5.1, it has been shown that the current distribution
$P(\mathcal{J}_{l})$ approaches Eq. (36) for thermodynamically large system
and the second moment of $\mathcal{J}_{l}$ is given by [24],

$\langle\mathcal{J}_{l}^{2}\rangle=2J_{\text{r}}+u_{l}^{2}.$ (44)

In Fig. 11(a), numerically measured $\langle\mathcal{J}_{l}^{2}\rangle$ is
shown as a function of $r_{1}$ for different values of $r_{N}$ along with the
analytical plot of Eq. (44), and they are in a good agreement.

We do not have the analytical expression of the second moment of
$\mathcal{J}_{1}$. However, $\langle\mathcal{J}_{1}^{2}\rangle$ is size-
independent as the probability distribution $P(\mathcal{J}_{1})$ is size-
independent. In Fig. 11(b) numerically measured
$\langle\mathcal{J}_{1}^{2}\rangle$ is shown as a function of $r_{1}$ for
different values of $r_{N}$.

##  6 Effective thermal limit

In the limit of a high resetting rate, the colored noise coming from the
resetting wall at the boundary becomes an effective white noise. Therefore an
effective fluctuation-dissipation relation can be written in terms of an
effective temperature, $T_{i}^{\text{eff}}$.

$\langle
f_{i}(t)f_{j}(t^{\prime})\rangle=\delta_{ij}a_{i}^{2}e^{-r_{i}|t-t^{\prime}|}\xrightarrow{r_{i}\to\infty}2\gamma\delta_{ij}T_{i}^{\text{eff}}\delta(t-t^{\prime})~{}\text{where,}~{}T_{i}^{\text{eff}}=\frac{a_{i}^{2}}{\gamma
r_{i}}.$ (45)

Therefore the resetting wall at the boundary acts as an effective Langevin
bath with effective temperature $T_{i}^{\text{eff}}$. In this high resetting
rate, the average energy flux in the NESS Eq. (34) becomes identical to the
thermally driven harmonic chain [1] with higher order correction,

$J_{\text{r}}=\frac{k(T_{1}^{\text{eff}}-T_{N}^{\text{eff}})}{2\gamma}\Bigg{[}1+\frac{mk}{2\gamma^{2}}-\frac{mk}{2\gamma^{2}}\sqrt{1+\frac{4\gamma^{2}}{mk}}\Bigg{]}+O\Big{(}\frac{1}{r_{j}^{2}}\Big{)}.$
(46)

Similarly, we can evaluate the two-time velocity correlation of a single
oscillator and equal-time spatial velocity correlation in the limit of a high
resetting rate. For large $r_{i}$, Eq. (22) can be approximated as,

$\langle
v_{l}(t)v_{l}(0)\rangle\simeq\frac{T_{1}^{\text{eff}}+T_{N}^{\text{eff}}}{2m}J_{0}(\omega_{c}t).$
(47)

Here $J_{0}(z)$ is the Bessel function of the first kind. Equal-time spatial
velocity correlation in the limit of high resetting rate can be evaluated
using Eqs.(25) and (26) and the result is,

$\langle
v_{l}(t)v_{l^{\prime}}(t)\rangle\simeq\frac{T_{1}^{\text{eff}}+T_{N}^{\text{eff}}}{2m}\delta_{ll^{\prime}}.$
(48)

In Sec. 4 we discuss the two-time velocity correlation of single oscillator
and equal-time spatial velocity correlation for resetting driven oscillator
chain. For the resetting driven chain, the nonzero value of $\langle
v_{l}(t)v_{l^{\prime}}(t)\rangle$ for $l\neq l^{\prime}$ in the NESS implies
that there exists a finite correlation between the velocity of the oscillators
in the bulk. However, in the effective thermal limit, this correlation
vanishes. Similarly, in the effective thermal limit, the two-time velocity
correlation of single oscillator $\langle v_{l}(t)v_{l}(0)\rangle$ is given in
Eq. (47). The small time behavior ($t\ll\omega^{-1}_{c}$) of two-time velocity
correlation $\langle v_{l}(t)v_{l}(0)\rangle$ is different for the resetting
driven chain [Fig. 6(a)] but the large time behavior is similar for both the
cases.

##  7 Conclusion

In this work, we study the NESS of a harmonic chain driven by resetting
dynamics at the boundary. The harmonic chain is attached to two walls at the
boundaries which undergo one-dimensional Brownian motion and stochastically
reset to their initial position at different rates. As a result, the boundary
oscillators experience exponentially correlated forces that lead to a NESS.
For the resetting driven harmonic chain, we conclude from the numerical
simulation that the velocity distributions of the bulk oscillators reach
Gaussian distribution when the system size is thermodynamically large. The
velocity distributions of the boundary oscillators are non-Gaussian and size-
independent. We have numerically computed the velocity kurtosis profile and
found that velocity kurtosis of the bulk oscillators decay as $N^{-1}$, where
$N$ is the system size. We also show that the two-time velocity correlations
of the bulk oscillators $\langle v_{l}(t)v_{l}(0)\rangle$ decay as
$t^{-\frac{1}{2}}$ when $t\gg\omega_{c}^{-1}$. In the bulk, there exists a
nonzero equal-time spatial velocity correlation $\langle
v_{l}(t)v_{l^{\prime}}(t)\rangle$ when $l\neq l^{\prime}$. The distribution of
the instantaneous energy current at the bulk is found to be size-dependent and
reaches a stationary distribution only when $N\to\infty$. However, The
distribution of instantaneous current at the boundary is size-independent and
depends only on the resetting rate of the boundary walls. It is important to
note that, at a high resetting rate, an effective thermal picture emerges
where we recover expression for the average energy current of a thermally
driven harmonic chain. In this effective thermal limit, equal-time spatial
velocity correlation $\langle v_{l}(t)v_{l^{\prime}}(t)\rangle$ vanishes when
$l\neq l^{\prime}$.

In this paper, the boundary walls perform free diffusion with Poisson
resetting where the waiting time between the two resetting events is drawn
from an exponential distribution. However, the dynamics of the boundary walls
can be generalized where the waiting time between two consecutive resetting
events follow a power-law distribution [34, 59, 60] or the resetting mechanism
becomes non-instantaneous [61, 62, 63, 64]. It is expected that such changes
in the dynamics of the boundary walls will affect the NESS of the ordered
harmonic chain and it is intriguing to explore the transport properties of
such NESS. Another interesting question is how our results change when instead
of our simple model of nonequilibrium bath, the boundary oscillators of the
harmonic chain are attached to explicit nonequilibrium baths that are modeled
using the stochastic resetting e.g. the model of nonequilibrium bath described
in Ref. [53]. It is also interesting to study how the energy transport in
resetting driven oscillator chains gets affected if nonlinearity, mass
disorder, and momentum non-conserving dynamics are introduced.

## Acknowledgments

The authors would like to thank Urna Basu and Ion Santra for useful
discussions. R.S. acknowledges support from the Council of Scientific and
Industrial Research, India [Grant No. 09/0575(11358)/2021-EMR-I].

##  Appendix A Matrix formulation

In this appendix, we will briefly discuss the correlation calculations. To
begin with, let us start with the matrix form the Langevin equations Eq. (1),

$M\ddot{X}=-\Phi X(t)-\Gamma\dot{X}(t)+F(t),$ (49)

where $X(t)$ is a state vector of dimension $N\times 1$, and the $l$-th
component is $x_{l}(t)$, displacement of the $l$-th oscillator from the
equilibrium position. $\Phi$ and $\Gamma$ are $N\times N$ matrix and the
definitions of $\Phi$ and $\Gamma$ are,

$\Phi_{ij}=2k\delta_{ij}-k\delta_{i,j-1}-k\delta_{i,j+1}\quad\text{and}\quad\Gamma_{ij}=\gamma\delta_{i1}\delta_{j1}+\gamma\delta_{iN}\delta_{jN}.$
(50)

$F(t)$ is the noise coming from the resetting wall. We can write these as,

$F_{j}(t)=f_{1}(t)\delta_{j1}+f_{N}(t)\delta_{jN}.$ (51)

The correlation of this noise in the frequency domain is,

$\langle\tilde{F}(\omega)\tilde{F}^{T}(\omega^{\prime})\rangle_{ij}=2\pi\delta(\omega+\omega^{\prime})\left[\tilde{g}(\omega,r_{1})\delta_{i1}\delta_{j1}+\tilde{g}(\omega,r_{N})\delta_{iN}\delta_{jN}\right],$
(52)

where $\tilde{g}(\omega,r_{j})=\frac{2a_{j}^{2}r_{j}}{r_{j}^{2}+\omega^{2}}$.
With the help of the Fourier transform, we can get rid of the time
derivatives. The Fourier transform of the displacement matrix,

$\tilde{X}(\omega)=\int^{\infty}_{-\infty}dte^{i\omega t}X(t).$ (53)

Therefore the Langevin equation Eq.(49) takes the form,

$\tilde{X}(\omega)=G(\omega)\tilde{F}(\omega)\quad\text{where}\quad
G(\omega)=[-M\omega^{2}+\Phi-i\omega\Gamma]^{-1}.$ (54)

$G(\omega)$ is the inverse of the tridiagonal matrix. The explicit form of
this matrix is,

$G(\omega)=\begin{bmatrix}-m\omega^{2}+2k-i\omega\gamma&-k&\cdots\\\
-k&-m\omega^{2}+2k&\cdots\\\ \vdots&\ddots&\cdots\\\
0&\cdots&-m\omega^{2}+2k-i\omega\gamma\\\ \end{bmatrix}^{-1}.$ (55)

$G$ is the inverse of a symmetric matrix, therefore $G$ it self is a symmetric
matrix and $G^{*}(\omega)=G(-\omega)$. We can write the components of
$G(\omega)$ [65],

$G_{l1}=(-k)^{l-1}\frac{\theta_{N-l}}{\theta_{N}}\quad\text{and}\quad
G_{lN}=(-k)^{N-l}\frac{\theta_{l-1}}{\theta_{N}}.$ (56)

where $\theta_{l}$ satisfy recursion relations,

$\theta_{l}=(-m\omega^{2}+2k)\theta_{l-1}-k^{2}\theta_{l-2}\quad\forall
l=2,3,\cdots N-1,$ (57)

and,

$\theta_{N}=(-m\omega^{2}+2k-i\omega\gamma)\theta_{N-1}-k^{2}\theta_{N-2},$
(58)

with the boundary condition,

$\theta_{0}=1,\qquad\theta_{1}=-m\omega^{2}+2k-i\omega\gamma.$ (59)

Solving the recursion relations of Eq. (57) we get,

$\theta_{l}=\frac{(-k)^{l-1}}{\sin{q}}[k\sin{(l+1)q}-i\omega\gamma\sin{lq}]\quad\text{and}\quad\theta_{N}=\frac{(-k)^{N}}{\sin{q}}[a(q)\sin{Nq}+b(q)\cos{Nq}].$
(60)

Here,

$a(q)=-\frac{2i\gamma\omega}{k}+\cos{q}\left(1-\frac{\gamma^{2}\omega^{2}}{k^{2}}\right)\quad\text{and}\quad
b(q)=\sin{q}\left(1-\frac{\gamma^{2}\omega^{2}}{k^{2}}\right),$ (61)

and $\omega$ and $q$ are related through,

$\cos{q}=\left(1-\frac{m\omega^{2}}{2k}\right),\quad\text{and}\quad\omega=\omega_{c}\sin{\frac{q}{2}}.$
(62)

with $\omega_{c}=2\sqrt{k/m}$.

##  Appendix B Calculation for spatio-temporal two-point correlation of
velocity

In this section, we will calculate the spatio-temporal two-point correlation
of the velocity of $l$-th oscillator $\langle
v_{l}(t)v_{l^{\prime}}(t^{\prime})\rangle$ in the steady state for resetting
driven oscillator chain. The velocity of the $l$-th oscillator of the
resetting driven chain can be written using Eqs. (53) and (54),

$v_{l}(t)=\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}e^{-i\omega
t}(-i\omega)\big{[}G_{l1}(\omega)\tilde{f}_{1}(\omega)+G_{lN}(\omega)\tilde{f}_{N}(\omega)\big{]}.$
(63)

Here $\tilde{f}_{j}(\omega)$ is the Fourier transform of the resetting force,
$f_{j}(t)$. Using Eq. (52), we can write the spatio-temporal correlation of
the velocity as,

$\langle
v_{l}(t)v_{l^{\prime}}(t^{\prime})\rangle=\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\omega^{2}e^{-i\omega(t-t^{\prime})}\big{[}G_{l1}(\omega)G_{l^{\prime}1}^{*}(\omega)\tilde{g}(\omega,r_{1})+G_{lN}(\omega)G_{l^{\prime}N}^{*}(\omega)\tilde{g}(\omega,r_{N})\big{]}.$
(64)

###  B.1 Two-time velocity correlation of single oscillator:

For $l=l^{\prime}$ and $t>t^{\prime}$,

$\langle
v_{l}(t)v_{l}(t^{\prime})\rangle=\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\omega^{2}e^{-i\omega(t-t^{\prime})}\big{[}|G_{l1}(\omega)|^{2}\tilde{g}(\omega,r_{1})+|G_{lN}(\omega)|^{2}\tilde{g}(\omega,r_{N})\big{]}.$
(65)

From the above equation, it is clear that the two-time velocity correlation of
single oscillator $\langle v_{l}(t)v_{l}(t^{\prime})\rangle$ for resetting
driven chain, can be written as a sum of two separate contributions coming
from two resetting walls at the two boundaries, i.e.,

$\langle
v_{l}(t)v_{l}(t^{\prime})\rangle=\bar{C}(r_{1},t,t^{\prime})+\bar{C}(r_{N},t,t^{\prime}).$
(66)

The contribution from the left resetting wall,

$\bar{C}(r_{1},t,t^{\prime})=\frac{1}{2\pi}\int_{-\infty}^{\infty}d\omega\omega^{2}e^{-i\omega(t-t^{\prime})}|G_{l1}(\omega)|^{2}\tilde{g}(\omega,r_{1}).$
(67)

$\bar{C}(r_{1},t,t^{\prime})$ will have a non-zero contribution from the even
terms with respect to $\omega$. Using Eq. (56) and keeping only the terms that
are even in $\omega$, we get,

$\displaystyle\bar{C}(r_{1},t,t^{\prime})$ $\displaystyle=$
$\displaystyle\frac{1}{k^{4}}\int_{0}^{\infty}\frac{d\omega}{\pi}\omega^{2}\cos{[\omega(t-t^{\prime})]}\Bigg{(}\frac{k^{2}\sin^{2}{(Nq-
lq+q)}}{|a(q)\sin{Nq}+b(q)\cos{Nq}|^{2}}$ (68) $\displaystyle+$
$\displaystyle\frac{\omega^{2}\gamma^{2}\sin^{2}{(Nq-
lq)}}{|a(q)\sin{Nq}+b(q)\cos{Nq}|^{2}}\Bigg{)}\tilde{g}(\omega,r_{1})$ (69)

We are interested in the time correlation of velocity in the bulk of the
chain, therefore we take $l=N/2+\epsilon$, and take the limit $\epsilon\ll N$,
which leads to,

$\displaystyle\bar{C}(r_{1},t,t^{\prime})$ $\displaystyle=$
$\displaystyle\frac{1}{k^{4}}\int_{0}^{\infty}\frac{d\omega}{2\pi}\omega^{2}\cos{[\omega(t-t^{\prime})]}\Bigg{(}\frac{(k^{2}+\omega^{2}\gamma^{2})}{|a(q)\sin{Nq}+b(q)\cos{Nq}|^{2}}$
(70) $\displaystyle-$ $\displaystyle\frac{k^{2}\cos{(Nq-2\epsilon
q+2q)}+\omega^{2}\gamma^{2}\cos{(Nq-2\epsilon
q)}}{|a(q)\sin{Nq}+b(q)\cos{Nq}|^{2}}\Bigg{)}\tilde{g}(\omega,r_{1}).$ (71)

For $\omega>\omega_{c}$, $q$ becomes complex. In the thermodynamically large
system size, the integrand vanishes exponentially as $e^{-2N\bar{q}}$ in the
region $\omega>\omega_{c}$, here $\bar{q}$ is real. The range of the
integration reduces to $0\leq\omega\leq\omega_{c}$, or $0\leq q\leq\pi$. In
the large-$N$ limit, $\cos{Nq}$ is a highly oscillatory function and we can
evaluate these integrations by averaging over fast oscillation in $x=Nq$,
using the following identities (see 2.558 in [66]),

$\displaystyle\frac{1}{2\pi}\int^{2\pi}_{0}\frac{dx}{(c_{1}\sin{x}+d\cos{x})^{2}+c_{2}^{2}\sin^{2}{x}}$
$\displaystyle=$ $\displaystyle-\frac{1}{c_{2}d},~{}\text{for}~{}c_{2}<0,$
$\displaystyle\frac{1}{2\pi}\int^{2\pi}_{0}\frac{dx\cos{x}}{(c_{1}\sin{x}+d\cos{x})^{2}+c_{2}^{2}\sin^{2}{x}}$
$\displaystyle=$ $\displaystyle 0.$ (72)

with $c_{2}=\text{Im}[a(q)]=-\frac{2\gamma\omega}{k}$ and
$d=\text{Re}[b(q)]=\sin{q}\left(1+\frac{\gamma^{2}\omega^{2}}{k^{2}}\right)$.
Therefore Eq. (71) can be written as,

$\displaystyle\bar{C}(r_{1},t,t^{\prime})=\frac{1}{k^{4}}\int_{0}^{\infty}\frac{d\omega}{2\pi}\omega^{2}\cos{[\omega(t-t^{\prime})]}\frac{k^{2}+\gamma^{2}\omega^{2}}{-dc_{2}}\tilde{g}(\omega,r_{1}).$
(73)

Finally, we arrive at,

$\displaystyle\bar{C}(r_{1},t,t^{\prime})$ $\displaystyle=$
$\displaystyle\frac{1}{k}\int_{0}^{\pi}\frac{dq}{\pi}\Big{|}\frac{d\omega}{dq}\Big{|}\frac{\omega\cos{[\omega(t-t^{\prime})]}}{4\gamma\sin{q}}\tilde{g}(\omega,r_{1})$
(74) $\displaystyle=$ $\displaystyle\frac{a_{1}^{2}r_{1}}{\gamma
m}\int_{0}^{\pi}\frac{dq}{2\pi}\frac{\cos\Big{[}\omega_{c}\sin{\Big{(}\frac{q}{2}\Big{)}}(t-t^{\prime})\Big{]}}{r_{1}^{2}+\omega_{c}^{2}\sin^{2}{\frac{q}{2}}}=a_{1}^{2}C(r_{1},t,t^{\prime}),$
(75)

where we have used the explicit form of $\tilde{g}(\omega,r_{1})$ and
$\omega(q)$ from Eq. (52) and (62). Similarly, we can evaluate the
contribution from the right resetting wall, and by combining these two results
we arrive at the,

$\langle v_{l}(t)v_{l}(t^{\prime})\rangle=\frac{1}{\gamma
m}\int_{0}^{\pi}\frac{dq}{2\pi}\cos\Big{[}\omega_{c}\sin{\Big{(}\frac{q}{2}\Big{)}}(t-t^{\prime})\Big{]}\Bigg{[}\frac{a_{1}^{2}r_{1}}{r_{1}^{2}+\omega_{c}^{2}\sin^{2}{\frac{q}{2}}}+\frac{a_{N}^{2}r_{N}}{r_{N}^{2}+\omega_{c}^{2}\sin^{2}{\frac{q}{2}}}\Bigg{]}.$
(76)

This integral can be evaluated numerically.

Large time behavior ($t-t^{\prime}\gg\omega_{c}^{-1}$): To calculate the large
time behavior of $\langle v_{l}(t)v_{l}(t^{\prime})\rangle$, we first
substitute $z=\omega_{c}\sin\Big{(}\frac{q}{2}\Big{)}$ in Eq. (75) and arrive
at,

$C(r_{1},t,t^{\prime})=\frac{r_{1}}{\gamma
m}\int_{0}^{\omega_{c}}\frac{dz}{\pi}\frac{\cos\big{(}z(t-t^{\prime})\big{)}}{\sqrt{\omega_{c}^{2}-z^{2}}(r_{1}^{2}+z^{2})}.$
(77)

For large $t-t^{\prime}$, $\cos\big{(}z(t-t^{\prime})\big{)}$ is a bounded and
fast oscillatory function. Therefore, the dominating contribution to the
integral is coming from the region $z=\omega_{c}$ and contributions from small
$z$ are negligible. We can now approximate $r_{1}^{2}+z^{2}$ as
$r_{1}^{2}+\omega_{c}^{2}$. Therefore the last integral becomes,

$C(r_{1},t,t^{\prime})\simeq\frac{r_{1}}{\gamma
m}\int_{0}^{\omega_{c}}\frac{dz}{\pi}\frac{\cos\big{(}z(t-t^{\prime})\big{)}}{\sqrt{\omega_{c}^{2}-z^{2}}(r_{1}^{2}+\omega_{c}^{2})}=\frac{r_{1}}{2\gamma}\frac{J_{0}\big{(}\omega_{c}(t-t^{\prime})\big{)}}{(mr_{1}^{2}+4k)}.$
(78)

Here $J_{0}(z)$ is the Bessel function of the first kind. Therefore when
$t-t^{\prime}$ is large, using Eqs. (75) and (78) we get

$\langle
v_{l}(t)v_{l}(t^{\prime})\rangle\simeq\Bigg{[}\frac{a_{1}^{2}r_{1}}{2\gamma(mr_{1}^{2}+4k)}+\frac{a_{N}^{2}r_{N}}{2\gamma(mr_{N}^{2}+4k)}\Bigg{]}J_{0}\big{(}\omega_{c}(t-t^{\prime})\big{)}.$
(79)

Kinetic temperature profile at the bulk: The bulk value of the kinetic
temperature profile can be evaluated using Eq. (76) by taking $t=t^{\prime}$.
Therefore, in the stationary state,

$\langle
v_{l}(t)^{2}\rangle=\frac{1}{\gamma}\int_{0}^{\pi}\frac{dq}{2\pi}\Bigg{[}\frac{a_{1}^{2}r_{1}}{mr_{1}^{2}+4k\sin^{2}{\frac{q}{2}}}+\frac{a_{N}^{2}r_{N}}{mr_{N}^{2}+4k\sin^{2}{\frac{q}{2}}}\Bigg{]}.$
(80)

The last integral can be evaluated exactly, which leads to the following
expression of the bulk value of the kinetic temperature profile,

$\hat{T}_{\text{bulk}}=m\langle
v_{l}^{2}(t)\rangle=\frac{a_{1}^{2}}{2\gamma\sqrt{r_{1}^{2}+\frac{4k}{m}}}+\frac{a_{N}^{2}}{2\gamma\sqrt{r_{N}^{2}+\frac{4k}{m}}}.$
(81)

###  B.2 Equal-time spatial velocity correlation:

For $l\neq l^{\prime}$ and $t=t^{\prime}$, using Eq. (64),

$\langle
v_{l}(t)v_{l^{\prime}}(t)\rangle=\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\omega^{2}\big{[}G_{l1}(\omega)G_{l^{\prime}1}^{*}(\omega)\tilde{g}(\omega,r_{1})+G_{lN}(\omega)G_{l^{\prime}N}^{*}(\omega)\tilde{g}(\omega,r_{N})\big{]}.$
(82)

Similar to the previous calculation, $\langle
v_{l}(t)v_{l^{\prime}}(t)\rangle$ can be written as a sum of two separate
contributions coming from resetting walls at the boundaries as the following,

$\langle
v_{l}(t)v_{l^{\prime}}(t)\rangle=\bar{\theta}(r_{1},l,l^{\prime})+\bar{\theta}(r_{N},l,l^{\prime}).$
(83)

The contribution from the left resetting wall,

$\bar{\theta}(r_{1},l,l^{\prime})=\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\omega^{2}G_{l1}(\omega)G_{l^{\prime}1}^{*}(\omega)\tilde{g}(\omega,r_{1}).$
(84)

Using Eq. (56) and keeping only the terms that are even in $\omega$, we get,

$\displaystyle\bar{\theta}(r_{1},l,l^{\prime})$ $\displaystyle=$
$\displaystyle\frac{1}{k^{4}}\int_{0}^{\infty}\frac{d\omega}{2\pi}\omega^{2}\Bigg{(}\frac{(k^{2}+\gamma^{2}\omega^{2})\cos{(l^{\prime}q-lq)}}{|a(q)\sin{Nq}+b(q)\cos{Nq}|^{2}}$
(85) $\displaystyle-$ $\displaystyle\frac{k^{2}\cos{(Nq-\epsilon
q+2q)+\omega^{2}\gamma^{2}\cos{(Nq-\epsilon
q)}}}{|a(q)\sin{Nq}+b(q)\cos{Nq}|^{2}}\Bigg{)}\tilde{g}(\omega,r_{1}).$ (86)

Now we are interested in $l=\frac{N}{2}$, $l^{\prime}=\frac{N}{2}+\alpha$
where $\alpha$ is an integer. For $\omega>\omega_{c}$, the integrand vanishes
[see the discussion after Eq. (71)]. The final expression after averaging over
fast oscillation using Eq. (72) is,

$\bar{\theta}(r_{1},l,l^{\prime})=\frac{1}{k^{4}}\int_{0}^{\pi}\frac{dq}{2\pi}\Big{|}\frac{d\omega}{dq}\Big{|}\omega^{2}\frac{(k^{2}+\gamma^{2}\omega^{2})\cos{(l^{\prime}-l)q}}{-dc_{2}}\tilde{g}(\omega,r_{1}).$
(87)

Using $c_{2}=\text{Im}[a(q)]=-\frac{2\gamma\omega}{k}$ and
$d=\text{Re}[b(q)]=\sin{q}\left(1+\frac{\gamma^{2}\omega^{2}}{k^{2}}\right)$,
and explicit form of $\tilde{g}(\omega,r_{1})$, we arrive at,

$\bar{\theta}(r_{1},l,l^{\prime})=\frac{a_{1}^{2}r_{1}^{2}}{\gamma}\int_{0}^{\pi}\frac{dq}{2\pi}\frac{\cos{(l^{\prime}q-lq)}}{mr_{1}^{2}+4k\sin^{2}\frac{q}{2}}=\frac{a_{1}^{2}}{2\gamma}\theta(r_{1},l,l^{\prime}),$
(88)

where,

$\theta(r_{i},l,l^{\prime})=\frac{r_{i}}{mr_{i}^{2}+4k}\,{}_{3}\tilde{F}_{2}\,\Bigg{[}1/2,1,1;1-l^{\prime}+l,1+l^{\prime}-l;\frac{4k}{4k+mr_{i}^{2}}\Bigg{]}.$
(89)

Here $l^{\prime}-l$ is an integer and ${}_{p}\tilde{F}_{q}$ is a generalized
regularized hypergeometric function [56]. The definition of generalized
regularized hypergeometric function ${}_{p}\tilde{F}_{q}$ is,

${}_{p}\tilde{F}_{q}[x_{1},x_{2}\cdots x_{p};y_{1},y_{2},\cdots
y_{q};z]=\frac{{}_{p}{F}_{q}[x_{1},x_{2}\cdots x_{p};y_{1},y_{2},\cdots
y_{q};z]}{\Gamma(y_{1})\Gamma(y_{2})\cdots\Gamma(y_{q})},$ (90)

where $\Gamma(z)$ is the gamma function and ${}_{p}{F}_{q}$ is the generalized
hypergeometric function. The definition of generalized hypergeometric function
${}_{p}{F}_{q}$ is,

${}_{p}{F}_{q}[x_{1},x_{2}\cdots x_{p};y_{1},y_{2},\cdots
y_{q};z]=\sum_{j=0}^{\infty}\frac{(x_{1})_{j}(x_{2})_{j}\cdots(x_{p})_{j}}{(y_{1})_{j}(y_{2})_{j}\cdots(y_{q})_{j}}\Big{(}\frac{z^{j}}{j!}\Big{)}.$
(91)

$(x)_{j}$ is the Pochhammer symbol defined as
$(x)_{j}\equiv\Gamma(x+j)/\Gamma(x)$ [56]. We can calculate the contribution
coming from the right resetting wall $\bar{\theta}(r_{N},l,l^{\prime})$ in a
similar manner. Combining both contributions we arrive at the,

$\langle
v_{l}(t)v_{l^{\prime}}(t)\rangle=\frac{1}{2\gamma}\Big{[}a_{1}^{2}\theta(r_{1},l,l^{\prime})+a_{N}^{2}\theta(r_{N},l,l^{\prime})\Big{]}.$
(92)

References

## References

  * [1] Rieder Z, Lebowitz J L and Lieb E 1967 J. Math. Phys. 8 1073
  * [2] Dhar A 2008 Adv. in Phys. 57 457
  * [3] Thermal Transport in Low Dimensions 2016, Ed. Stefano Lepri, Springer Heidelberg
  * [4] Nakazawa H 1970 Prog. Theor. Phys. Suppl. 45 231
  * [5] Roy D and Dhar A 2008 J. Stat. Phys. 131 535
  * [6] Dhar A 2001 Phys. Rev. Lett. 86 5882
  * [7] Lepri S, Livi R and Politi A 2005 Chaos 15 015118
  * [8] Mai T, Dhar A and Narayan O 2007 Phys. Rev. Lett. 98 184301
  * [9] Kundu A, Sabhapandit S and Dhar A 2011 J. Stat. Mech. P03007
  * [10] Kannan V, Dhar A and Lebowitz J L 2012 Phys. Rev. E 85 041118
  * [11] Kubo R 1966 Rep. Prog. Phys. 29 255
  * [12] Valeriani C, Li M, Novosel J, Arlta J and Marenduzzoa D 2011 Soft Matter 7 5228
  * [13] Gopal A, Roldán É and Ruffo S 2021 J. Phys. A: Math. Theor. 54 164001
  * [14] Granek O, Kafri Y and Tailleur J 2022 Phys. Rev. Lett. 129 038001
  * [15] Maggi C, Paoluzzi M, Pellicciotta N, Lepore A, Angelani L and Di Leonardo R 2014 Phys. Rev. Lett. 113 238303
  * [16] Maes C 2020 Phys. Rev. Lett. 125 208001
  * [17] Seyforth H, Gomez M, Rogers W. B., Ross J L and Ahmed W W 2022 Phys. Rev. Research 4 023043
  * [18] Mousavi S M, Gompper G and Winkler R G 2021 J. Chem. Phys. 155 044902
  * [19] Pal A and Sabhapandit S 2014 Phys. Rev. E 90 052116
  * [20] Fodor É, Nemoto T and Vaikuntanathan S 2020 New J. Phys. 22 013052
  * [21] Chaki S and Chakrabarti R 2019 Physica A: Stat. Mech. Appl. 530 121574
  * [22] Santra I 2023 J. Phys. Complex. 4, 015013
  * [23] Santra I and Basu U 2022 Scipost Phys. 13 041
  * [24] Sarkar R, Santra I and Basu U 2023 Phys. Rev. E 107 014123
  * [25] Evans M R, Majumdar S N and Schehr G 2020 J. Phys. A: Math. Theor. 53 193001
  * [26] Evans M R and Majumdar S N 2011 Phys. Rev. Lett. 106 160601
  * [27] Magoni M, Majumdar S. N., and Schehr G 2020 Phys. Rev. Research 2, 033182
  * [28] Ray S and Reuveni S 2020 J. Chem. Phys 152, 234110
  * [29] Masoliver J 2019 Phys. Rev. E 99, 012121
  * [30] Kuśmierz L, Majumdar S N, Sabhapandit S, and Schehr G 2014 Phys. Rev. Lett. 113, 220602
  * [31] Majumdar S N, Mounaix P, Sabhapandit S and Schehr G 2022 J. Phys. A: Math. Theor. 55, 034002
  * [32] Durang X, Henkel M, and Park H 2014 J. Phys. A: Math. Theor. 47, 045002
  * [33] Basu U, Kundu A and Pal A 2019 Phys. Rev. E 100, 032136
  * [34] Mishra S and Basu U 2023 J. Stat. Mech. 053202
  * [35] Jain S, Boyer D, Pal A, and Dagdug L 2023 J. Chem. Phys. 158, 054113
  * [36] Di Bello C, Hartmann A K, Majumdar SN, Mori F, Rosso A and Schehr G 2023 Phys. Rev. E 108, 014112
  * [37] Evans M R, Majumdar S N 2018 J. Phys. A: Math. Theor. 51, 475003
  * [38] Santra I, Basu U, and Sabhapandit S 2020 J. Stat. Mech. 113206
  * [39] Stojkoski V, Sandev T, Kocarev L, and Pal A 2021 Phys. Rev. E 104, 014121
  * [40] Singh P, and Pal A 2021 Phys. Rev. E 103, 052119
  * [41] Gupta D, Plata C A, and Pal A 2020 Phys. Rev. E 124, 110608
  * [42] Gupta S, Majumdar S.N., and Schehr G 2014 Phys. Rev. Lett. 112, 220601
  * [43] Gupta S and Nagar A 2016 J. Phys. A: Math. Theor. 49, 445001
  * [44] Tal-Friedman O, Pal A, Sekhon A, Reuveni S and Roichman Y 2020 J. Phys. Chem. Lett. 11, 7350-7355
  * [45] Besga B, Bovon A, Petrosyan A, Majumdar S N and Ciliberto S 2020 Phys. Rev. Res. 2, 032029
  * [46] Faisant F, Besga B, Petrosyan A, Ciliberto S and Majumdar S N 2021 J. Stat. Mech. 113203
  * [47] Montanari A and Zecchina R 2002 Phys. Rev. Lett. 88, 178701
  * [48] Tong H, Faloutsos C, and Pan J Y 2008 Knowl. Inf. Syst. 14, 327
  * [49] Avrachenkov K, Piunovskiy A, and Zhang Y, 2013 J. Appl. Prob. 50, 960
  * [50] Perfetto G, Carollo F, Magoni M, and Lesanovsky I 2021 Phys. Rev. B 104 L180302
  * [51] Magoni M, Carollo F, Perfetto G, and Lesanovsky I 2022 Phys. Rev. A 106, 052210
  * [52] Dubey V, Chetrite R and Dhar A 2023 J. Phys. A: Math. Theor. 56 154001
  * [53] Maes C and Thiery T 2017 J. Phys. A: Math. Theor. 50, 415001
  * [54] Majumdar S N and Oshanin G 2018 J. Phys. A: Math. Theor. 51 435001
  * [55] Vanden-Eijndena E and Ciccotti G 2006 Chem. Phys. Lett. 429, 310
  * [56] NIST Digital Library of Mathematical Functions, Olver F W J, Olde Daalhuis A B, Lozier D W, Schneider B I, Boisvert R F, Clark C W, Miller B R, Saunders B V, Cohl H S, and McClain M A, eds. http://dlmf.nist.gov/, Release 1.1.6 of 2022-06-30
  * [57] Kundu A 2010 Phys. Rev. E 82 031131
  * [58] Goswami K 2019 Phys. Rev. E 99 012112
  * [59] Nagar A and Gupta S 2016 Phys. Rev. E, 93, 060102
  * [60] Nagar A and Gupta S 2023 J. Phys. A: Math. Theor., 56, 283001
  * [61] Besga B., Bovon A, Petrosyan A, Majumdar S. N. and Ciliberto S. 2020 Phys. Rev. Res., 2, 032029
  * [62] Gupta D, Plata C. A., Kundu A, and Pal A 2020. J. Phys. A: Math. Theor., 54, 025003.
  * [63] Santra I, Das S and Nath S. K. 2021 J. Phys. A: Math. Theor., 54, 334001
  * [64] Gupta D, Pal A and Kundu A 2021 J. Stat. Mech., 043202
  * [65] Usmani R. A 1994 Comput. Math. Appl. 27, 59
  * [66] Gradshteyn L. S. and Ryzhik I. M. 1943 Academic Press, New York, US, Table of integrals, series, and products, ISBN 9780123736376