**footnotetext: The speaker has been supported from the EU’s Horizon 2020
research and innovation programme under grant agreement No 952480 (DarkWave),
and from the International Research Agenda Pro- gramme AstroCeNT (MAB/2018/7)
funded by the Foundation for Polish Science from the European Regional
Development Fund.

# Status of direct detection searches of heavy dark matter

Marcin Kuźniak

###### Abstract

Dark matter accounts for 26% of the mass-energy density of the Universe,
however, its nature and origins remain the most important open questions in
physics. The search for Weakly Interacting Massive Particles (WIMPs), one of
the leading dark matter particle candidates, is now in a decisive phase, with
experiments targeting both the high-mass and the low-mass (1–10 GeV/$c^{2}$)
WIMP scenarios. The status of the leading experimental searches is discussed,
together with their prospects and challenges they are facing. Searches of
heavy non-WIMP dark matter candidates are also briefly summarized.

##  1 Introduction

###  1.1 Case for heavy dark matter

The existence of dark matter (DM) is supported by a plethora of observational
evidence. It constitutes some 26% of the mass-energy balance of the Universe,
is in the form of non-baryonic, non-dissipative and, most likely, cold DM, and
yet its exact nature remains one of the main remaining puzzles in physics.
Among the various candidates, Weakly Interacting Massive Particle (WIMP)
produced thermally via freeze-out, called standard or thermal WIMP, is one of
the favoured DM candidates and received much focus over the last decades,
alongside axions. Assuming electroweak strength interactions with the Standard
Model sector and respecting the relic density constraints, implies order of
GeV/$c^{2}$ to TeV/$c^{2}$ WIMP mass, potentially consistent with neutralino,
the lightest stable particle expected from a large class of SUSY-based or
other types of SM extensions (for a recent review, see e.g. [1]). Relaxing
somewhat the assumptions on the interaction strength, the messenger particles
involved and the DM production mechanism (thermal vs. not), shifts the lower
mass constraint down to the eV/$c^{2}$ scale, with multiple possible DM
scenarios beyond the scope of this report. At the same time, super-heavy DM,
with mass from multi-TeV/$c^{2}$ all the way up to the Planck scale (1019
GeV/$c^{2}$), may result from more exotic scenarios, e.g. inflaton decay or
gravitational mechanisms related to inflation (often related to Grand Unified
Theories), millicharged gravitino DM, DM produced by primordial black hole
radiation or extended thermal production in the dark sector [2].

###  1.2 Direct detection

The principle of direct detection of heavy dark matter particles originated 40
years ago from the paper by Drukier and Stodolsky [3], where the coherent
enhancement of the neutrino elastic scattering cross-section on nuclei was
first calculated, resulting in values much larger than previously anticipated,
and within long-term experimental reach. In addition to establishing
foundations for the CE$\upnu$NS measurement (only recently accomplished [4]),
that work was soon generalized by Goodman and Witten [5] to the case of WIMPs,
which led to the first dark matter detector concepts [6], and the first
experimental limits [7].

Since then, more than 7 orders of magnitude of the WIMP-nucleon scattering
cross-section parameter space have been experimentally probed. Ultimately, the
sensitivity of non-directional direct searches is limited by irreducible
backgrounds from solar and atmospheric neutrinos, the so-called neutrino fog
[8].

##  2 Direct detection techniques

The experimental signatures of WIMP dark matter are single-site low-energy
(keV) nuclear recoils, as well as the annual modulation of the dark matter
signal, caused by changes in the velocity of Earth with respect to WIMP wind,
and resulting changes in event rate above the energy threshold.

For the case of more exotic super-heavy dark matter candidates, multi-site
nuclear recoil events or electronic recoil events are also a usable
experimental signature, if sufficiently distinct from the expected
backgrounds.

In essence, dark matter detectors seek for any measurable effects in the
target medium caused by the recoil events depositing their energy. These
typically include scintillation, ionization and heat. Sensitivity to more than
one of these channels, on the one hand, comes at the cost of increased readout
complexity, additional materials and associated radioactive impurities, on the
other hand, provides stronger particle identification and background
discrimination handles. As background events caused by instrumental effects or
natural sources of radioactivity, (e. g. traces of radioactive isotopes in the
detector materials or its surroundings) are the main challenge, designing a
competitive dark matter detector is a non-trivial balancing act.

###  2.1 Crystal detectors

A conceptually simple strategy of using an $\mathcal{O}$(100 kg) matrix of
scintillating NaI(Tl) crystals viewed by photomultiplier tubes (PMT) was
adopted by the DAMA/LIBRA experiment [9]. This is also the case for ANAIS-112
[10], SABRE [11] and COSINE-100 [12], which aim to cross-check the DAMA/LIBRA
signal claim. In this configuration, nuclear recoil events (WIMP-like) are not
distinguishable from electronic recoils events (background). COSINUS [13]
experiment is extending that approach by operating the sodium iodide crystals
cryogenically as bolometers, with added phonon (temperature) readout, which
will enable event-by-event discrimination between e-/$\upgamma$ and nuclear
recoil events.

The strategy of reading out multiple interaction channels has been routinely
used by experiments using other types of crystals. For instance, EDELWEISS
[14] uses cryogenic Ge, and CDMS/SuperCDMS [15] employs Ge and Si crystals,
instrumented with ionization and phonon readout, for e-/$\upgamma$ background
discrimination. A fraction of crystals is operated at higher voltage with
ionization sensed indirectly via Luke-Neganov phonon production, which enables
reaching lower energy thresholds at the expense of background discrimination
power. The CRESST [17] collaboration operates cryogenic scintillating
bolometers based on calcium tungstate or lithium aluminate crystals, which
provide background discrimination capabilities based on the scintillation and
phonon signal ratio.

Typically, large arrays, stacks or "towers" of individually calibrated
crystals, often of mixed types, are operated in a single cryostat, providing
additional capability to veto multi-scatter events. Their common feature is
very low $\mathcal{O}$(10 eV) threshold, but also a difficulty in scaling-up
to large target mass, due to material availability and cost.

###  2.2 CCDs

Silicon based charged coupled devices (CCDs) equipped with skipper readout
have been developed and operated by DAMIC [18] and SENSEI [19] collaborations.
While essentially sensitive to ionization, dense arrays of CCD pixels provide
a powerful way to identify particles and discriminate backgrounds based on the
hit pattern and track topology. Due to low energy threshold and good
background discrimination, this technology is competitive in the low-mass WIMP
regime.

###  2.3 Noble element detectors

Over the past decade, liquid argon or xenon based detectors have dominated the
sensitivity for WIMP masses $>$1 GeV/$c^{2}$. Some of the main reasons include
the ease of purification and scale-up to multi-tonne target masses, high
scintillation and ionization yields as well as powerful self-shielding and
background discrimination capabilities, for Xe and Ar, respectively.
Furthermore, Xe has competitive sensitivity for spin-dependent scattering,
particularly on neutrons.

#### Time projection chambers

TPCs, such as currently running XENONnT [20], LZ [21], PandaX-4T [22] with
LXe, or with LAr DarkSide-50 [23] (already decommissioned), or DarkSide-20k
[24] (currently under construction), register primary scintillation (S1) and
secondary scintillation (S2) from drifted electrons, employing the e- drift
for position reconstruction and fiducialization. The scintillation to
ionization ratio is the main electronic recoil discrimination handle in LXe.
Gaseous TPCs (although not only noble element based) are also the primary
technology of choice for directional detection, providing position resolution
better than the nuclear recoil track length.

#### Calorimeters

Calorimeters such as XMASS [25] with LXe, or currently-running DEAP-3600 [26]
with LAr, rely exclusively on scintillation. A big advantage of this approach
is simplicity and scalability. Uniquely in LAr, pulse-shape discrimination of
S1 signals provides a powerful electronic recoil discrimination power, on the
level of 109 at threshold [27].

###  2.4 Bubble chambers

In bubble chambers process of nuclear recoil induced heat deposition,
nucleation and bubble formation in overheated liquids is employed, with the
thermodynamical conditions of the system optimized to suppress any events from
electronic recoils. Fluorine-rich target liquids used by PICO [28] give it
leading sensitivity for spin-dependent (SD) interactions. Scintillating Bubble
Chamber (SBC) [29] has prototyped a detector combining the features of a LAr
calorimeter and bubble chamber, with a potential to reach low energy
thresholds.

###  2.5 Paleodetectors

Examining ancient minerals for traces of WIMP-nucleus interactions recorded
over timescales as large as 1 Gyr may provide competitive exposure (and at
least in principle, circumvent the limitations of the neutrino fog even
without directional sensitivity). R&D and small-scale demonstrations of this
approach are ongoing [30].

##  3 Status of experimental searches

###  3.1 DAMA/LIBRA detection claim

Very significant efforts invested over the past decade into independent cross-
check of the DAMA/LIBRA signal claim [9] are now beginning to bring fruits.

The claim has been already challenged by modulation searches in other targets.
In LXe it was excluded by XENON100, LUX, and XMASS [25]. Recently, DarkSide-50
[31] reported the first annual modulation search in LAr, with the lowest
achieved energy threshold of 0.04 keV: compatible with no modulation, but
neither confirming nor rejecting DAMA.

Searches using the same target material as DAMA are quickly catching up.
COSINE-100 [12] excluded most of the signal region with the model-dependent
approach and reported a model-independent modulation search with 2.8 yr of
data, consistent with both DAMA and no modulation; results from doubled
exposure at a lower energy threshold are expected soon. ANAIS-112 [10]
provided the first sensitive model-independent cross-check: incompatible with
DAMA at $>$3.9$\upsigma$, consistent with no modulation; it is expecting to
reach 5$\upsigma$ sensitivity by the end of 2025.

###  3.2 Low-mass WIMPs

Recently improved S2-only DarkSide-50 analysis resulted in the leading spin-
independent (SI) sensitivity in the 1.2 to 10 GeV/$c^{2}$ range [23], see Fig.
1 (top), with XENON1T, CDMSLite [15], CRESST-III and DAMIC up to 2 orders of
magnitude worse. Taking advantage of the Migdal effect to reduce the
considered energy threshold extends the leading sensitivity of DarkSide-50
down to 0.04 GeV/$c^{2}$ [32]. Having said that, Migdal effect has not been
experimentally confirmed yet [33] and its presence is, in principle, an open
question investigated experimentally by multiple groups [34]. In the SD sector
CRESST-III [17] holds the most stringent exclusion at 1 GeV/$c^{2}$ mass,
while CDMSLite [16] and PICASSO [35] at 3 GeV/$c^{2}$ for scattering on
neutron and proton, respectively. Beyond 3 GeV/$c^{2}$ XENON1T [36] and
PICO-60 [28] take over.

###  3.3 Standard thermal WIMPs

In the high-mass regime, exposure is the critical factor and noble liquid
detectors, primarily LXe detectors are unchallenged, see Fig. 1. Early results
from the new generation of LXe experiments, XENONnT [20], LZ [21] and
PandaX-4T [22], have already pushed down the exclusion limits by more than a
factor of 2 with respect to the previous generation of detectors. The only
exception is the SD WIMP-proton scattering, where PICO-60 [28] has the leading
edge, due to highly fluorinated target. In case of LAr, the DEAP-3600
experiment is targeting reducing the sensitivity gap, following the ongoing
upgrade addressing the limiting backgrounds [26].

Figure 1: Status of direct detection searches: (top) spin-indepenent, (bottom
left) spin-dependent WIMP-neutron scattering, (bottom right) spin-dependent
WIMP-proton scattering. Generated with DM Limit Plotter

###  3.4 Planck-scale DM

Due to a very distinct expected event signature, with high deposited energy
and a large number of nuclear recoils all along the dark matter particle track
in the detector, this type of searches is essentially background-free, with
the sensitivity driven by the cross-sectional area of a detector. This gives
the leading SI sensitivity at Planck-mass scale to DEAP-3600 [2], while at
lower masses and cross-sections, or overall for the SD case, to XENON1T [39].

##  4 Future prospects

#### Low-mass WIMP region and low-energy excess

Cryogenic bolometers and LAr TPCs are currently best positioned to push
towards the neutrino fog for WIMP masses $<$10 GeV/$c^{2}$. For bolometers, a
limiting factor currently is a low-energy event excess observed by multiple
technologically different experiments (including DAMIC, SuperCDMS, CRESST)
[40]. The community has joined forces to address and resolve the issue in a
collaborative manner. The mitigation strategy for such events and the exact
origin of the excess remain open issues, with hints that energy release from
long-lived metastable states is at play. A concept of DarkSide-LM, a tonne-
scale LAr TPC optimized for S2-only operation and low mass sensitivity, has
been proposed [41], inspired by DarkSide-50. Also in this case, a low-energy
excess was observed in DarkSide-50, and at least in part attributed to
impurities, capturing electrons and releasing them with a delay. Success of
this project depends on securing funding, and on the efficiency of the
impurity removal and spurious electron event mitigation.

#### Directional searches

Directional sensitivity allows to distinguish the expected WIMP signal from
neutrino backgrounds, although so far it was only demonstrated in gaseous
detectors with low target mass, and far from competitive exposure. Possibility
to employ the columnarity effect in LAr to, at least partly, gain directional
sensitivity in a scalable medium has been recently excluded [42]. However, the
community has converged into the CYGNUS collaboration, with the goal of scale-
up. Particularly, in the very-low WIMP mass regime, for WIMP-electron
scattering, or for spin-dependent scenarios interesting sensitivities could be
reached [43].

#### "Ultimate" large detectors

Above 10 GeV/$c^{2}$, exposure of 200 tonne-year in LXe or 3000 tonne-year in
LAr would be needed to reach the neutrino fog. Both communities have
consolidated efforts, with LZ, XENON and DARWIN, merging into XLZD [44] with
targeted 60 tonne LXe mass, and independently PandaX-xT [45] with $>$30 tonne
mass. The LAr community has merged into the Global Argon Dark Matter
Collaboration, aiming to build a 300 tonne fiducial mass detector, ARGO [1],
following the commissioning of DarkSide-20k. While technological challenges to
resolve remain for both targets (and the community, guided by ECFA, has
organized to address them [46]), due to excellent pulse shape discrimination,
in particular for events from elastic solar neutrino-electron scattering, LAr
is in principle capable of achieving background-free exposure, which would
result in a higher 5$\upsigma$ discovery sensitivity [1].

##  5 Summary

Despite the lack of uncontested WIMP signal, the field has made a tremendous
progress over the past few years, with a new set of exclusion limits from the
current generation of detectors recently published. In addition to this, many
important efforts are progressing, with the columnarity in LAr recently
excluded, and the conclusive verification of the Migdal effect and a cross-
check of the DAMA/LIBRA claim already on the horizon. While the roadmap
towards the neutrino floor remains sound, significant challenges remain ahead,
particularly in the context of understanding and mitigating the low energy
backgrounds. Growing consolidation of efforts in the community gives many
reasons for optimism and hope that 40 years into direct detection searches we
now hawe most of the distance to the neutrino fog already behind us.

## References

  * [1] J. Billard et al., Rep. Prog. Phys. 85 (2022) 056201
  * [2] P. Adhikari et al, Phys. Rev. Lett. 128 (2022) 011801
  * [3] A. Drukier, L. Stodolsky, Phys. Rev. D 30 (1984) 2295
  * [4] D. Akimov et al., Science 357 (2017) 1123
  * [5] M. Goodman, E. Witten, Phys. Rev. D 31 (1985) 3059
  * [6] A. K. Drukier et al., Phys. Rev. D 33 (1986) 3495
  * [7] S. P. Ahlen et al., Phys. Lett. B 195 (1987) 603
  * [8] C. O’Hare, Phys. Rev. Lett. 127 (2021) 251802
  * [9] R. Bernabei et al., Nucl. Part. Phys. Proc 303–305 (2018) 74
  * [10] I. Coarasa et al., arXiv:2311.03392 (2023)
  * [11] E. Barberio et al., PoS(ICRC2023) 1370
  * [12] G. Adhikari et al, Phys. Rev. D106 (2022) 052005
  * [13] G. Angloher et al., arXiv:2307.11139 (2023)
  * [14] J. Gascon, J. Phys.: Conf. Ser. 1468 (2020) 012018
  * [15] R. Agnese et al., Phys. Rev. D 99 (2019) 062001
  * [16] R. Agnese et al., Phys. Rev. D 97 (2018) 022002
  * [17] G. Angloher et al., Phys. Rev. D 106 (2022) 092008
  * [18] A. Aguilar-Arevalo et al., Phys. Rev. Lett. 125 (2020) 241803
  * [19] L. Barak et al., Phys. Rev. Lett. 125 (2020) 171802
  * [20] E. Aprile et al., Phys. Rev. Lett. 131 (2023) 041003
  * [21] J. Aalbers et al., Phys. Rev. Lett.131(2023) 041002
  * [22] Y. Meng et al., Phys. Rev. Lett. 127 (2021) 261802
  * [23] P. Agnes et al., Phys. Rev. D 107 (2023) 063001
  * [24] P. Agnes, EPJ Web of Conf. 280 (2023) 06003
  * [25] K. Abe et al., Phys. Rev. D 108 (2023) 083022
  * [26] M. Lai, JINST 18 (2023) C02046
  * [27] P. Adhikari et al., Eur. Phys. J. C 81 (2021) 823
  * [28] C. Amole et al., Phys. Rev. D 100 (2019) 022001
  * [29] B. Broermann, SciPost Phys. Proc.12 (2023) 023
  * [30] S. Baum et al., Instruments 5 (2021) 21
  * [31] P. Agnes et al., arXiv:2307.07249 (2023)
  * [32] P. Agnes et al., Phys. Rev. Lett. 130 (2023) 101001
  * [33] J. Xu et al., arXiv:2307.12952 (2023)
  * [34] H. M. Araújo, Astropart. Phys. 151 (2023) 102853
  * [35] E. Behnke et al., Astropart. Phys. 90 (2017) 85
  * [36] E. Aprile et al., Phys. Rev. Lett. 123 (2019) 251801
  * [37] W. Ma et al., Phys. Rev. Lett. 130 (2023) 021802
  * [38] Y. Meng et al., Phys. Rev. Lett. 130 (2023) 261001
  * [39] E. Aprile et al., Phys. Rev. Lett. 130 (2023) 261002
  * [40] P. Adari, SciPost Phys. Proc. 9 (2022) 001
  * [41] P. Agnes et al., Phys. Rev. D 107 (2023) 112006
  * [42] P. Agnes et al., arXiv:2307.15454 (2023)
  * [43] S. Vahsen et al., Annu. Rev. Nucl. Part. Sci. 71 (2021) 189
  * [44] J. Aalbers et al., J. Phys. G: Nucl. Part. Phys. 50 (2023) 013001
  * [45] X. Wang et al., JINST 18 (2023) P05028
  * [46] ECFA Detector R&D Roadmap Task Force 2: https://indico.cern.ch/event/1214404


Activity driven transport in harmonic chains

Ion Santra1, Urna Basu1,2*,

1 Raman Research Institute, C. V. Raman Avenue, Bengaluru 560080, India

2 S. N. Bose National Centre for Basic Sciences, JD Block, Saltlake, Kolkata
700106, India

<EMAIL_ADDRESS>
## Abstract

The transport properties of an extended system driven by active reservoirs is
an issue of paramount importance, which remains virtually unexplored. Here we
address this issue, for the first time, in the context of energy transport
between two active reservoirs connected by a chain of harmonic oscillators.
The couplings to the active reservoirs, which exert correlated stochastic
forces on the boundary oscillators, lead to fascinating behavior of the energy
current and kinetic temperature profile even for this linear system. We
analytically show that the stationary active current (i) changes non-
monotonically as the activity of the reservoirs are changed, leading to a
negative differential conductivity (NDC), and (ii) exhibits an unexpected
direction reversal at some finite value of the activity drive. The origin of
this NDC is traced back to the Lorentzian frequency spectrum of the active
reservoirs. We provide another physical insight to the NDC using
nonequilibrium linear response formalism for the example of a dichotomous
active force. We also show that despite an apparent similarity of the kinetic
temperature profile to the thermally driven scenario, no effective thermal
picture can be consistently built in general. However, such a picture emerges
in the small activity limit, where many of the well-known results are
recovered.

###### Contents

  1. I Introduction
  2. II Model
  3. III Results
    1. III.1 Stationary energy current
      1. III.1.1 Negative differential conductivity
      2. III.1.2 Current reversal
    2. III.2  Kinetic temperature
  4. IV Conclusions
  5. A Stationary state Current
  6. B Kinetic temperature profile
    1. B.1 $\hat{T}_{l}$ near left boundary
    2. B.2 $\hat{T}_{l}$ near right boundary
  7. C Linear response: Differential conductivity

##  I Introduction

Understanding energy transport properties of driven systems is a central issue
of nonequilibrium statistical physics. Theoretical attempts in this regard
often rely on the study of simple, yet analytically tractable model systems
DharReview2008 ; Transportbook . A paradigmatic example is a chain of harmonic
oscillators connected to thermal reservoirs of different temperatures at the
two ends, first studied by Rieder, Lieb, and Lebowitz (RLL) in a seminal work
RLL . They showed that this system reaches a nonequilibrium stationary state
carrying a thermal current, which survives in the thermodynamic limit. Several
generalizations of this simple model have been studied by introducing
disorders, anharmonic interactions, pinning potentials and activity in the
bulk Nakazawa ; Dhar2001 ; RoyDhar2008 ; FPUT ; FPUT_alternatingmass ; ldf ;
Kannan2012 ; reservoirchain3 . In almost all of these studies, however, the
reservoirs attached to the system are taken to be equilibrium ones — the
random and dissipative forces exerted by each reservoir on the boundary
oscillators satisfy the Fluctuation-Dissipation theorem (FDT) Kubo .

Nonequilibrium reservoirs, on the other hand, do not respect any such FDT,
giving rise to a wide range of new possibilities maes2013 ; maes2014 ;
maes2015 ; vandebroek ; Sabhapandit2017 . For example, energy transport in
systems connected to nonequilibrium reservoirs show non-monotonic kinetic
temperature profile, negative differential thermal conductivity and non-
reciprocal heat transport Iacobucci2011 ; prosen2011 ; Bagchi2013 ;
Hayakawa2013 . Active reservoirs refer to a special class of nonequilibrium
reservoirs, consisting of self-propelled particles like bacteria or Janus
beads, which are inherently out of equilibrium by consuming energy from the
environment at an individual level Libchaber2000 ; Soni2003 ; SoodNature2016 .
Recent studies, both theoretical and experimental, show that individual probe
particles immersed in such active reservoirs exhibit many unusual features
including emergence of negative friction, modification of equipartition
theorem and anomalous relaxation dynamics bacterialbath2011 ; maggi2014 ;
gopal2021 ; maes2020 ; kafri2021 ; ahmed2021 ; ABP_polymer ; work_activebath ;
Fodor2020 ; Chakrabarti2019 . A natural question is how the transport
properties of an extended system are affected when connected to active
reservoirs at the boundaries. To the best of our knowledge, this has not been
studied so far.

In this article, we ask this question in a simple setting similar to RLL
model—an ordered chain of harmonic oscillators connected to two active
reservoirs at the two ends. The active reservoirs exert stochastic forces on
the boundary oscillators, which do not satisfy FDT. As a simple model, we
consider that this stochastic force has an exponentially decaying
autocorrelation, which is a common feature of active dynamics, the
autocorrelation time-scale being a measure of the activity of the reservoirs.
In the long-time limit the system reaches a nonequilibrium stationary state
(NESS) carrying an energy current which we compute exactly. We find that this
current shows two remarkable features, namely, an unexpected direction
reversal and a negative differential conductivity (NDC) whose origin lies with
the Lorentzian frequency spectra of the active reservoirs. The emergence of
the NDC and current reversal in a linear system without any kinetic
constraints sets it apart from the few similar phenomena observed
previouslyIacobucci2011 ; Casati2006 ; Franosch2013 ; ndc_2013 ;
chatterjee2018_ndc ; chatterjee2018_aem . For a specific model of a
dichotomous active force, we illustrate that the NDC can also be viewed as a
result of a positive correlation of the current and the number of directional
flips of the force. We also show that the kinetic temperature profile retains
strong signatures of activity despite attaining a uniform value at the bulk.
In the limit of small activity, the reservoirs behave somewhat similar to
thermal ones and the well-known properties of RLL-model are recovered.

##  II Model

We consider a one-dimensional chain with $N$ particles, each with mass $m$,
connected by harmonic springs of stiffness $k$, attached to two different
active reservoirs at the boundaries [see Fig. 1]. The coupling to the active
reservoir is modeled by including a stochastic force on the boundary particle,
in addition to the usual dissipative and white-noise forces coming from an
equilibrium thermal reservoir. The equations of motion for $x_{l}$, the
displacement of the $l$-th particle from its equilibrium position, read,

$\displaystyle m\ddot{x}_{1}$ $\displaystyle=$
$\displaystyle-k(2x_{1}-x_{2})-\gamma\,\dot{x_{1}}+\xi_{1}(t)+f_{1}(t),$ (1a)
$\displaystyle m\ddot{x}_{l}$ $\displaystyle=$
$\displaystyle-k(2x_{l}-x_{l-1}+x_{l+1}),~{}~{}~{}\forall\,l\in[2,N-1],$ (1b)
$\displaystyle m\ddot{x}_{N}$ $\displaystyle=$
$\displaystyle-k(2x_{N}-x_{N-1})-\gamma\,\dot{x_{N}}+\xi_{N}(t)+f_{N}(t),~{}~{}~{}$
(1c)

where we have used fixed boundary conditions $x_{0}=0=x_{N+1}.$ We assume that
the thermal components of the reservoirs are at temperatures $T_{1}$ and
$T_{N}$, so that the white noises $\xi_{1,N}(t)$ acting on the boundary
particles $l=1,N$ are related to the dissipation $\gamma$ through FDT,

$\langle\xi_{l}(t)\xi_{j}(t^{\prime})\rangle=2\gamma
T_{i}\delta_{l,j}\delta(t-t^{\prime}).$ (2)

The FDT is violated by the presence of the active forces $f_{1,N}(t)$ which
are assumed to be independent stationary colored noises. Most commonly, such
active noises have an exponentially decaying correlation, $\langle
f_{j}(t)f_{l}(t^{\prime})\rangle=\delta_{jl}\,a_{j}^{2}\exp(-|t-t^{\prime}|/\tau_{j})$,
where $a_{j}$ denotes the strength of the noise and the correlation-time
$\tau_{j}$ is a measure of the activity. As a specific example, we consider
the dichotomous noise

$\displaystyle f_{j}(t)=a_{j}\sigma_{j}(t),$ (3)

where $\sigma_{j}$ alternates between $\pm 1$ at a constant rate $\alpha_{j}$,
giving rise to an exponential correlation with $\tau_{j}=1/(2\alpha_{j}).$
However, our main results remain quite robust for general active driving, as
such exponential correlations generically appear in active processes including
run-and-tumble motion, active Brownian motion and direction reversing active
Brownian motion rtp ; abp ; drabp .

Figure 1: Schematic representation of a chain of oscillators connected to two
nonequilibrium reservoirs which exert active forces $f_{1,N}(t)$ on the
boundary oscillators.

##  III Results

We first present a brief summary of our main results. The primary observables
of interest here are the energy current and the kinetic temperature profile,
both of which we compute exactly. The energy current flowing through the
system is most conveniently expressed as DharReview2008 ,

$\displaystyle
J\equiv\langle\mathscr{J}(t)\rangle=\big{\langle}\dot{x}_{1}[-\gamma\dot{x}_{1}+\xi_{1}(t)+f_{1}(t)]\big{\rangle},$
(4)

where the average is over the NESS. Because of the linear nature of the
equations of motion the stationary current naturally separates into two
components, an _active_ one induced by the activity driving and a thermal one
proportional to the temperature difference of the two reservoirs (same as in
the usual RLL setup RLL ). We show that the active current in the
thermodynamic limit is given by,

$\displaystyle J_{\text{act}}$ $\displaystyle=$
$\displaystyle\frac{m}{2\gamma^{2}}(a_{1}^{2}\,\mathscr{E}_{1}-a_{N}^{2}\,\mathscr{E}_{N})\quad\text{with},$
(5) $\displaystyle\mathscr{E}_{j}$ $\displaystyle=$
$\displaystyle\frac{\tau_{j}^{2}k^{2}\left[\sqrt{1+\frac{4\gamma^{2}}{mk}}-1\right]+\gamma^{2}\left[1-\sqrt{1+\frac{4k\tau_{j}^{2}}{m}}\right]}{2\tau_{j}(\tau_{j}^{2}k^{2}-\gamma^{2})}.$
(6)

There are a number of striking features of this active current which
distinguishes it from the usual thermal current. First, $J_{\text{act}}$
exhibits a non-monotonic behavior as the activity of either of the reservoirs
is changed, giving rise to a negative differential conductivity [see Fig. 2].
More surprisingly, the current reverses its direction as the activity of one
of the reservoirs, say $\tau_{1}$, is changed at a non-trivial value
$\tau_{1}^{*}\neq\tau_{N}$ [see the phase diagram in Fig. 4].

We also show that the stationary kinetic temperature profile
$\hat{T}_{l}=m\langle\dot{x}_{l}^{2}\rangle$ attains a constant value in the
bulk $1\ll l\ll N$ with an exponentially decaying boundary layer.
Surprisingly, we find that, the bulk temperature can be expressed in a form
similar to the famous RLL result RLL ,

$\displaystyle\hat{T}_{bulk}=\frac{1}{2}(\mathscr{T}_{1}+\mathscr{T}_{N}),\quad\text{with}~{}~{}\mathscr{T}_{j}=\frac{a_{j}^{2}\tau_{j}}{\gamma\sqrt{1+4\tau_{j}^{2}k/m}}.$
(7)

This suggests the possibility of interpreting $\mathscr{T}_{1,N}$ as
‘effective temperatures’ associated to the two active reservoirs. However, we
show that such an interpretation is not acceptable and the active reservoirs
remain essentially different than thermal ones.

In the limit of small activity $\tau_{1},\tau_{N}\ll\sqrt{m/k}$, however, an
effective thermal picture emerges. In this case, we show that, the active
forces behave somewhat similar to white noises and the energy current and bulk
kinetic temperature are consistent with the system being connected to thermal
reservoirs with effective temepartures
$T_{j}^{\text{eff}}=T_{j}+a_{j}^{2}\tau_{j}/\gamma$. However, the signatures
of activity still remain in some atypical features, like the presence of a
non-trivial boundary layer even when $T_{1}^{\text{eff}}=T_{N}^{\text{eff}}$.

We start by rewriting Eqs. (1) as,

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

where $X(t)=\\{x_{l}(t);l=1,\dotsc,N\\}$ is a vector and $M$ is an
$N$-dimensional diagonal matrix with $M_{lj}=m\delta_{l,j}$. Moreover,
$\Gamma$ and $\Phi$ are $N$-dimensional matrices given by

$\displaystyle\Gamma_{lj}$ $\displaystyle=$
$\displaystyle\gamma(\delta_{l,1}\delta_{j,1}+\delta_{l,N}\delta_{j,N}),\Phi_{lj}=k\left(2\delta_{l,j}-\delta_{l,j-1}-\delta_{l,j+1}\right).$

Finally, $\Xi_{j}=\xi_{1}(t)\delta_{j1}+\xi_{N}(t)\delta_{jN}$ and
$F_{j}=f_{1}(t)\delta_{j1}+f_{N}(t)\delta_{jN}$ are vectors.

Figure 2: Activity induced current $J_{\text{act}}$ vs $\tau_{1}$ for
different values of $\tau_{N}$ and $\gamma=1=m=a_{1}=a_{N}$ and $k=2$; symbols
indicate the data obtained from numerical simulations with $N=64$ oscillators
and the solid black lines indicate the analytical prediction Eq. (6). The red
crosses mark the non-trivial current reversal points.

We are interested in the solution of Eq. (8) in the stationary state, which is
most conveniently obtained by taking a Fourier transform with respect to time,
$\tilde{X}(\omega)=\int_{-\infty}^{\infty}dt\,e^{i\omega t}X(t)$. This leads
to,

$\displaystyle\tilde{X}(\omega)$ $\displaystyle=$ $\displaystyle
G(\omega)\left[\tilde{\Xi}(\omega)+\tilde{F}(\omega)\right],$ (9)

where $G(\omega)=[-M\omega^{2}+\Phi-i\omega(\Gamma_{L}+\Gamma_{R})]^{-1}$.
Here, $\tilde{\Xi}(\omega)$ and $\tilde{F}(\omega)$ are the Fourier transforms
of $\Xi(t)$ and $F(t)$ respectively. Inverting the transform, we get from Eq.
(9),

$\displaystyle X(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}d\omega\,e^{-i\omega
t}G(\omega)\big{[}\tilde{\Xi}+\tilde{F}\big{]}.$ (10)

To compute the steady state energy current $J$ defined in Eq. (4), we need the
autocorrelation of the stochastic forces $\xi_{j}(t)$ and $f_{j}(t)$ in the
Fourier-space,

$\displaystyle\langle\tilde{\xi}_{j}(\omega)\tilde{\xi}_{l}(\omega^{\prime})\rangle$
$\displaystyle=$ $\displaystyle 4\pi\gamma
T_{j}\delta_{jl}\delta(\omega+\omega^{\prime}),$ (11a)
$\displaystyle\langle\tilde{f}_{j}(\omega)\tilde{f}_{l}(\omega^{\prime})\rangle$
$\displaystyle=$ $\displaystyle
2\pi\delta_{jl}\delta(\omega+\omega^{\prime})\tilde{g}(\tau_{j},\omega).$
(11b)

Here
$\tilde{g}(\tau_{j},\omega)=\frac{2a_{j}^{2}\tau_{j}}{1+\omega^{2}\tau_{j}^{2}}$
denotes the spectral density of the active force from the $j$th reservoir,
which clearly is a Lorentzian with corner frequency $\tau_{j}^{-1}.$

###  III.1 Stationary energy current

The independence of the thermal and active noises along with the linear nature
of the couplings lead to the current in Eq. (4) to separate into two
components $J=J_{\text{th}}+J_{\text{act}}$; see Appendix. A for details. The
thermal current, generated due to the temperature gradient,

$\displaystyle J_{\text{th}}$ $\displaystyle=$
$\displaystyle\gamma^{2}(T_{1}-T_{N})\int_{0}^{\infty}\frac{d\omega}{\pi}\,\omega^{2}|G_{1N}(\omega)|^{2},$
(12)

remains same as in the case of equilibrium reservoirs and can be computed
explicitly RLL ; Dhar2001 . The active nature of the reservoirs gives rise to
the additional current,

$\displaystyle J_{\text{act}}$ $\displaystyle=$
$\displaystyle\gamma\int_{0}^{\infty}\frac{d\omega}{\pi}\,\omega^{2}|G_{1N}(\omega)|^{2}\Big{[}\tilde{g}(\tau_{1},\omega)-\tilde{g}(\tau_{N},\omega)\Big{]},~{}~{}~{}$
(13)

where $\tilde{g}(\tau_{j},\omega)$ contains information about the reservoir
activity. Equation (13) is a Landauer-like formula, where the transmission
coefficient depends explicitly on the Lorentzian reservoir spectra
$\tilde{g}(\tau_{j},\omega)$ along with the system phonon spectrum
$\omega^{2}G_{1N}(\omega)$.

To compute the currents explicitly we need $G_{1N}(\omega)$, which is obtained
exploiting the tridiagonal structure of $G^{-1}(\omega)$ Dhar2001 ;
RoyDhar2008 ; Kannan2012 ; Usmani . We are particularly interested in the
thermodynamic limit $N\to\infty$, where $G_{1N}(\omega)$ vanishes
exponentially outside the phonon band $|\omega|>\omega_{c}=2\sqrt{k/m}$
RoyDhar2008 . In that limit, we show that, the contribution from the $j$-th
reservoir is given by [see Appendix. A for details],

$\displaystyle\gamma\int_{0}^{\infty}\frac{d\omega}{\pi}\,\omega^{2}|G_{1N}(\omega)|^{2}\tilde{g}(\tau_{j},\omega)$
$\displaystyle=$
$\displaystyle\int_{0}^{\pi}\frac{dq}{\pi}\frac{mka_{j}^{2}\tau_{j}\sin^{2}q}{[mk+2\gamma^{2}(1-\cos
q)][m+2k\tau_{j}^{2}(1-\cos q)]},\quad$ (14)

where $\omega$ and $q$ are related by $m\omega^{2}=2k(1-\cos q)$. Computing
the $q$-integral and combining the contributions from both the reservoirs, we
get the active current flowing through the system in the thermodynamic limit
which is quoted in Eq. (6).

Figure 2 shows a plot of the predicted $J_{\text{act}}$ as a function of the
left reservoir activity $\tau_{1}$ for a set of different values of
$\tau_{N}$. This shows an excellent match with the current measured from
numerical simulations with a chain of oscillators driven by the dichotomous
noise given in Eq. (3). The figure illustrates some remarkable features of the
active current which we discuss below.

Figure 3: (a) Plot of the system phonon band $\omega^{2}|G_{1N}(\omega)|^{2}$
in the $N\to\infty$ limit [see Eq. (50)] and the reservoir spectrum
$\tilde{g}(\tau,\omega)$ [see Eq. (30)] for $\tau=0.5$ as functions of
$\omega$. (b) Plot of the single reservoir transmission coefficient
$\omega^{2}|G_{1N}(\omega)|^{2}\tilde{g}(\tau,\omega)$ vs $\omega$ for
different values of $\tau$. Here we have taken $m=1,k=2$ and $\gamma=1.$

####  III.1.1 Negative differential conductivity

The active current shows a non-monotonic behavior—as $\tau_{1}$ is increased,
$J_{\text{act}}$ initially increases until reaching a maximum value after
which it starts to decrease. It is clear from Eq. (6) that this non-monotonic
behavior is inherent to the individual contributions from both the reservoirs
— if $\tau_{N}$ is increased, keeping $\tau_{1}$ fixed, a similar behavior is
seen where the current first decreases and then starts to increase. The
existence of this non-monotonic behavior becomes qualitatively clear by
looking at the frequency spectrum of the reservoir $\tilde{g}(\tau,\omega)$.
From Eq. (11b), it is clear that $\tilde{g}(\tau,\omega)$ is a Lorentzian,
peaked around $\omega=0$ with width $\sim\tau^{-1}$. On the other hand, the
system phonon band is peaked around the characteristic frequency
$\omega_{c}=2\sqrt{k/m}$, with a minimum at $\omega=0$ [see Fig. 3(a)].
Consequently, the overlap of the system and reservoir spectra changes non-
monotonically as $\tau$ is changed, reaching a maximum at some intermediate
value of $\tau^{-1}\in[0,\omega_{c}]$ [see Fig. 3(b)]. This, in turn, gives
rise to the non-monotonic behavior of $J_{\text{act}}$, which shows a maximum
(minimum) as $\tau_{1}$ ($\tau_{N}$) is varied. In fact, it can be easily seen
from Eq. (6) that for large $k$, the current is maximum at a value of
$\tau_{1}\propto\omega_{c}^{-1}$.

The non-monotonic behavior implies that the differential conductivity
$\chi_{j}=\frac{dJ_{\text{act}}}{d\tau_{j}}$, which is nothing but the linear
response of the current to a small change in the activity of the $j$-th
reservoir, becomes negative in some parameter regimes. Noneqilibrium response
theory provides a way to express this coefficient in terms of correlations of
some physical observables MaesPRL2009 ; Maes2020Front . For the simple
dynamics (3), using a trajectory based approach, we find [see Appendix C],

$\displaystyle\chi_{j}=\lim_{t\to\infty}-\frac{1}{\tau_{j}}[\langle
n_{j}(t)\mathscr{J}(t)\rangle-\langle
n_{j}(t)\rangle\langle\mathscr{J}(t)\rangle]$ (15)

where $n_{j}(t)$ denotes the total number of flips of $\sigma_{j}$ during a
time interval $[0,t]$, $\mathscr{J}(t)$ is the instantaneous current and the
average is computed in the unperturbed system. The above equation implies that
when the number of flips $n_{j}(t)$ is positively correlated with the current
an NDC emerges.

####  III.1.2 Current reversal

There is another, more striking, behavior induced by the presence of the
active driving, namely, reversal of the direction of the current. We see from
Fig. 2, that for any given $\tau_{1}$, $J_{\text{act}}$ reverses its direction
twice—once (trivially) at $\tau_{1}=\tau_{N}$ and again at another value
$\tau_{1}=\tau_{1}^{*}$ which depends non-trivially on $\tau_{N}$. For a fixed
$\tau_{N}$, $J_{\text{act}}$ begins with a negative value (energy flowing from
right to left reservoir) for $\tau_{1}=0$, which becomes positive (energy
flowing from left to right reservoir) with increase in $\tau_{1}$. However, on
increasing $\tau_{1}$ further, the current again reverses its direction and
becomes negative. Mathematically, this additional reversal can be understood
from the observation that for a fixed value of $\tau_{N}$, $\mathscr{E}_{1}\to
0$ for both $\tau_{1}\to 0$ and $\tau_{1}\to\infty$ [see Eq. (6)], and
consequently $J_{\text{act}}$ has the same negative value at these two limits.
Now, since $J_{\text{act}}$ must reverse sign at $\tau_{1}=\tau_{N}$, an
additional reversal is required to reach the limiting negative values. A
similar scenario is observed when $\tau_{N}$ is changed keeping $\tau_{1}$
fixed, as expected from the symmetry of the system.

This behavior is illustrated in Fig. 4; panel (a) shows a three-dimensional
plot of $J_{\text{act}}$ on the $(\tau_{1},\tau_{N})$ plane, while Fig. 4(b)
shows the two-dimensional projection of (a) indicating the regions
$J_{\text{act}}>0$ and $J_{\text{act}}<0$. For any given $\tau_{N}$, the
current reverses its direction at $\tau_{1}=\tau_{N}$ and another non-trivial
point $\tau_{1}=\tau_{1}^{*}(\tau_{N})$. The latter is given by the non-
trivial solution of
$a_{1}^{2}\mathscr{E}_{1}(\tau_{1})=a_{N}^{2}\mathscr{E}_{N}(\tau_{N})$.
Similarly, for any given $\tau_{1}$, the current reversal occurs at
$\tau_{N}=\tau_{1}$ and $\tau_{N}^{*}(\tau_{1})$ [indicated by the solid red
curve in 4(b)]. Interestingly, the intersection of the curves
$\tau_{1}=\tau_{N}$ and $\tau_{1}=\tau_{1}^{*}(\tau_{N})$ denoted by
$(\bar{\tau},\bar{\tau})$ is a saddle point, as can be seen from Fig. 4(a).
The current does not change direction when one passes through the saddle
point—for $\tau_{N}=\bar{\tau}$, the current remains negative for all values
of $\tau_{1}\neq\tau_{N}$, while for $\tau_{1}=\bar{\tau}$, the current
remains positive for all values of $\tau_{N}\neq\tau_{1}$.

Figure 4: (a) Three dimensional plot of $J_{\text{act}}$ on the
$(\tau_{1},\tau_{N})$ plane: the meshed orange surface denotes
$J_{\text{act}}$ given by Eq. (5) in the main text, while the un-meshed, semi-
transparent blue surface corresponds to $J_{\text{act}}=0$. The curves formed
by intersection of these two surfaces give the locus of the zeros of the
active current, the intersection of these two curves are denoted by
$(\bar{\tau},\bar{\tau})$, which clearly is a saddle point of
$J_{\text{act}}$. $J_{\text{act}}>0$ in the region above the blue surface,
while $J_{\text{act}}<0$ in the region below the blue surface. (b) Two
dimensional projection of (a) showing phase diagram of $J_{\text{act}}$ in the
$(\tau_{1},\tau_{N})$ plane—: The light blue (green) shade indicates the
region where the active current is negative (positive). The continuous red
curve shows $\tau_{N}^{*}$ as a function of $\tau_{1}$ whereas the dashed
curve indicates the line $\tau_{N}=\tau_{1}$.

NDC and current reversal have been observed in certain nonequilibrium systems
with non-linearity, presence of obstacles or kinetic constraints Iacobucci2011
; Casati2006 ; Franosch2013 ; ndc_2013 ; chatterjee2018_ndc ;
chatterjee2018_aem . Surprisingly, the dynamical active driving here gives
rise to both features even in a linear chain.

###  III.2  Kinetic temperature

The average kinetic energy of the oscillators provides a way to define a local
‘temperature’ for driven oscillator chains RLL ; Transportbook . For purely
thermal drive, this kinetic temperature is uniform in the bulk of the system
and is given simply by $(T_{1}+T_{N})/2$ in the $N\to\infty$ limit. Here we
are interested in the effect of the active drive on the the kinetic
temperature $\hat{T}_{l}=m\langle\dot{x}_{l}^{2}(t)\rangle$ and thus consider
$T_{1}=T_{N}=0.$ In this case, using Eq. (10), we get,

$\displaystyle\hat{T}_{l}=m\int\frac{d\omega}{2\pi}\omega^{2}\Big{[}|G_{l1}(\omega)|^{2}\tilde{g}(\tau_{1},\omega)+|G_{lN}(\omega)|^{2}\tilde{g}(\tau_{N},\omega)\Big{]}.~{}~{}~{}~{}~{}$
(16)

The matrix elements can again be computed exploiting the tridiagonal structure
of $G^{-1}(\omega)$. Performing a similar calculation as before [see Appendix.
B for details], we find that, in the thermodynamic limit $N\to\infty$, the
steady state temperature profile is flat at the bulk, accompanied by
exponentially decaying boundary layers. The bulk temperature
$\hat{T}_{\text{bulk}}$ can be obtained explicitly and is quoted in Eq. (7).
The predicted value of bulk temperatures for a fixed $\tau_{1}$ and different
values of $\tau_{N}$, compared with numerical simulations performed with Eq.
(3) in Fig. 5(a) shows an excellent agreement. Interestingly, boundary kinks
in the $\hat{T}_{l}$ profile, which are generically present for coupling with
thermal reservoirs LepriLiviPoliti , are absent here.

Figure 5: (a) The kinectic temperature profile $\hat{T}_{l}$ for $\tau_{1}=1$
and different values of $\tau_{N}$ measured from simulations with a chain of
$N=64$ oscillators driven by the active noises (3). The dashed black lines
show the predicted bulk temperature (7). (b) Comparison of $J_{\text{act}}$
(solid lines) with the expected current for ‘effective’ temperature gradient
$\mathscr{T}_{1}-\mathscr{T}_{N}$ (dashed lines) for two different values of
$\tau_{N}$. Here $m=1,k=1,\,\gamma=1$ and $a_{1}=a_{N}=1$.

The form of Eq. (7) raises a possibility of associating an effective
temperature $\mathscr{T}_{j}$ to the $j$-th active reservoir. At first glance,
this identification also appears to be consistent with a ‘zeroth law’ — when
$\tau_{1}=\tau_{N}$, i.e., $\mathscr{T}_{1}=\mathscr{T}_{N}$, the bulk of the
system is at the same ‘temperature’ as the reservoirs. However, such an
interpretation is not acceptable for several reasons. First, note that the
kinetic temperature of the boundary sites $\hat{T}_{1,N}$ remain different
than $\hat{T}_{\text{bulk}}$ giving rise to a boundary layer even when
$\tau_{1}=\tau_{N}$ [see Fig. 5(a)] which is absent for ordinary equilibrium
reservoirs. Moreover, the stationary active current (6) is very different than
the energy current which would have been generated if the system were
connected to thermal reservoirs of temperatures $\mathscr{T}_{1}$ and
$\mathscr{T}_{N}$ at the two ends. This is illustrated in Fig. 5(b) which
shows neither current reversal nor any NDC in the ‘effective’ thermal
scenario. However, the effective temperature picture becomes viable in the
limit of small activity, which we discuss next.

Passive limit- It is well known that active systems show an effective passive
behavior in the limit of vanishing correlation time rtp ; abp ; drabp .
Similarly, in our case, when $\tau_{j}\to 0$, the active force $f_{j}(t)$
resembles a white noise with effective correlation $\langle
f_{j}(t)f_{j}(t^{\prime})\rangle\to a_{j}^{2}\tau_{j}\delta(t-t^{\prime})$. In
this limit, the active forces in Langevin Eqs. (1) can be thought of
representing thermal reservoirs with effective temperatures
$a_{j}^{2}\tau_{j}/\gamma$ and satisfying FDT. The well known results of the
RLL model are expected to be recovered in this ‘thermal’ limit. Indeed we see
from Eq. (7) that when the active time-scales are much smaller than the
coupling time-scale, i.e., $\tau_{1},\tau_{N}\ll\sqrt{m/k}$, the kinetic
temperature associated with the reservoirs $T^{\text{eff}}_{j}\simeq
a_{j}^{2}\tau_{j}/\gamma$ are consistent with the thermal picture. Moreover,
in this limit, it can be easily seen from Eq. (6) that,

$\displaystyle
J_{\text{act}}=\frac{k(T^{\text{eff}}_{1}-T^{\text{eff}}_{N})}{2\gamma}\left[1+\frac{mk}{2\gamma^{2}}-\frac{mk}{2\gamma^{2}}\sqrt{1+\frac{4\gamma^{2}}{mk}}\right]+O(\tau_{j}^{2}),$

which is same as the well-known form of the thermal current RLL ; RoyDhar2008
to leading order in $\tau_{1},\,\tau_{N}$. This can also be seen from Fig.
5(b) where $J_{\text{act}}$ converges to the effective thermal current for
$\tau_{1},\tau_{N}\ll\sqrt{m/k}$.

##  IV Conclusions

In summary, we have analytically studied the transport properties of a
harmonic chain coupled to two active reservoirs which exert exponentially
correlated stochastic forces on the boundary oscillators. We find that this
active drive leads to an NESS carrying an energy current, which exhibits
intriguing features like NDC and current reversal. For a simple model of
dichotomous active force, we show that the negative differential conductivity
results from a positive correlation of the energy current and number of flips
of the active force. The kinetic temperature profile, which, similar to the
thermally driven scenario, remains uniform at the bulk of the system, also
carries strong signatures of activity and an effective temperature picture
cannot be consistently built.

Our work is the first to study the effect of active reservoirs on transport
properties of extended systems. The results presented here are quite robust as
the exponential correlation is a generic feature of active dynamics. However,
signatures of specific dynamics are expected to be seen in the fluctuations of
the current. It would be interesting to see if our results can be
qualitatively verified in experiments with active reservoirs, say a collection
of active Brownian particles ABP_polymer , connected by passive polymers. Some
other interesting questions are: What are the effects of disorder,
anharmonicity and pinning in the presence of active driving? How do our
results change, if the nonequilibrium reservoir is modeled by a chain of
active particles in the spirit of Caprini2020 ; reservoirchain2 ;
reservoirchain3 ?

## Acknowledgements

The authors would like to thank Abhishek Dhar, Christian Maes and P. K.
Mohanty for useful discussions.

##### Funding information

U. B. acknowledges support from the Science and Engineering Research Board
(SERB), India, under a Ramanujan Fellowship (Grant No. SB/S2/RJN-077/2018).

##  Appendix A Stationary state Current

In this section, we sketch the main steps of the computation of the current
starting from Eq. (4) in the main text. For the sake of completeness we first
rewrite the Langevin equations Eqs. (1),

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

where, $X(t)=\\{x_{l}(t);l=1,\dotsc,N\\}$ is a vector, $M$ is an
$N$-dimensional diagonal matrix with $M_{lj}=m\delta_{l,j}$; $\Phi$ and
$\Gamma$ are $N$-dimensional matrices given by

$\displaystyle\Phi_{jl}$ $\displaystyle=$ $\displaystyle
k\left(2\delta_{j,l}-\delta_{j,l-1}-\delta_{j,l+1}\right),$ (18)
$\displaystyle\Gamma$ $\displaystyle=$
$\displaystyle\Gamma_{L}+\Gamma_{R}\quad~{}\text{with}\qquad\qquad(\Gamma_{L})_{jl}=\gamma\delta_{j,1}\delta_{l,1},\quad(\Gamma_{R})_{jl}=\gamma\delta_{j,N}\delta_{l,N}.$
(19)

Moreover, the vectors $\Xi(t)$ and $F(t)$ represent the thermal and active
forces exerted by the reservoirs on the boundary oscillators,

$\displaystyle\Xi(t)=\Xi_{L}(t)+\Xi_{R}(t)\quad$ with
$\displaystyle\quad(\Xi_{L})_{j}(t)=\xi_{1}(t)\delta_{j1}\quad\text{and}\quad(\Xi_{R})_{j}(t)=\xi_{N}(t)\delta_{jN},$
(20a) $\displaystyle F(t)=F_{L}(t)+F_{R}(t)\quad$ with
$\displaystyle\quad(F_{L})_{j}(t)=f_{1}(t)\delta_{j1}\quad\text{and}\quad(F_{R})_{j}(t)=f_{N}(t)\delta_{jN}.$
(20b)

Here, $\xi_{1,N}(t)$ are delta correlated white-noises, while the active
noises $f_{1,N}(t)$ have an exponentially decaying auto-correlation,

$\displaystyle\langle\xi_{j}(t)\xi_{l}(t^{\prime})\rangle$ $\displaystyle=$
$\displaystyle\delta_{jl}\,2\gamma
T_{j}\delta(t-t^{\prime}),\quad\text{and},\quad\langle
f_{j}(t)f_{l}(t^{\prime})\rangle=\delta_{jl}\,a_{j}^{2}e^{-|t-t^{\prime}|/\tau_{j}}.$
(21)

Note that, even though Eq. (17) formally appears to be a limiting case of
reservoirchain3  with vanishing bulk activity, the two scenarios differ by
their physical nature as well as emergent phenomena, as we will see below.

The stationary energy current flowing through the system can be expressed as
$J=\langle\mathscr{J}(t)\rangle$ where

$\displaystyle\mathscr{J}(t)=\dot{x}_{1}[-\gamma\dot{x}_{1}+\xi_{1}(t)+f_{1}(t)]$
(22)

denotes the instantaneous work done by, the left reservoir on the left
boundary oscillator and the statistical averaging is done over the stationary
state. It is convenient to recast this energy current using the above matrix
notation and separate it into two terms,

$\displaystyle J$ $\displaystyle=$ $\displaystyle
J_{1}+J_{2}\quad\text{with},\quad
J_{1}=-\text{Tr}\left[\langle\dot{X}(t)\dot{X}^{T}(t)\Gamma_{L}\rangle\right]\quad\text{and}\quad
J_{2}=\text{Tr}\left[\langle(\Xi_{L}+F_{L})\dot{X}^{T}(t)\rangle\right],$ (23)

where $\dot{X}^{T}$ denotes the transpose of the vector $\dot{X}$. In the
following we compute $J_{1}$ and $J_{2}$ separately using the solution of Eq.
(17),

$\displaystyle X(t)=\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\,e^{-i\omega
t}G(\omega)[\Xi(\omega)+F(\omega)].$ (24)

where $G(\omega)=[-M\omega^{2}+\Phi-i\omega(\Gamma_{L}+\Gamma_{R})]^{-1}$ [see
Eq. (10)]. Let us first consider,

$\displaystyle J_{1}$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\int_{-\infty}^{\infty}\frac{d\omega^{\prime}}{2\pi}\,\omega\,\omega^{\prime}\,e^{-i(\omega+\omega^{\prime})t}\,\text{Tr}\left[\langle\tilde{X}(\omega)\tilde{X}^{T}(\omega^{\prime})\rangle\Gamma_{L}\right]$
(25) $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\int_{-\infty}^{\infty}\frac{d\omega^{\prime}}{2\pi}\,\omega\,\omega^{\prime}\,e^{-i(\omega+\omega^{\prime})t}\,\text{Tr}\left[G(\omega)\left\langle[\Xi(\omega)+F(\omega)][\Xi(\omega^{\prime})+F(\omega^{\prime})]\right\rangle
G(\omega^{\prime})\Gamma_{L}\right],$ (26)

where we have used the fact that $G^{T}(\omega^{\prime})=G(\omega^{\prime})$
as $G$ is a symmetric matrix. The noise correlations appearing in the above
equation can be evaluated in a straightforward manner using Eqs. (20)-(21).
Since the noises from the two reservoirs are independent, it is natural to
separate the corresponding contributions and write,

$\displaystyle\left\langle\big{[}\Xi(\omega)+F(\omega)\big{]}\big{[}\Xi(\omega^{\prime})+F(\omega^{\prime})\big{]}\right\rangle=2\pi\delta(\omega+\omega^{\prime})\big{[}S_{L}(\omega)+S_{R}(\omega)\big{]},$
(28)

where the matrix elements of $S_{L,R}(\omega)$ are given by,

$\displaystyle\big{(}S_{L}(\omega)\big{)}_{jl}$ $\displaystyle=$
$\displaystyle[2\gamma
T_{1}+\tilde{g}(\tau_{1},\omega)]\delta_{j,1}\delta_{l,1},\quad\text{and}\quad(H_{R})_{jl}=[2\gamma
T_{N}+\tilde{g}(\tau_{N},\omega)]\delta_{j,N}\delta_{l,N}.$ (29)

Here $\tilde{g}(\tau_{j},\omega)$ denotes the Fourier transform of the active
force auto-correlation

$\displaystyle\tilde{g}(\tau_{j},\omega)=a_{j}^{2}\int_{-\infty}^{\infty}ds~{}e^{i\omega
s}e^{-|s|/\tau_{j}}=\frac{2a_{j}^{2}\tau_{j}}{(1+\omega^{2}\tau_{j}^{2})}.$
(30)

Using Eqs. (28) and (29) in Eq. (LABEL:A:j1_1), we get,

$\displaystyle
J_{1}=-\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\,\omega^{2}\text{Tr}\Big{[}G(\omega)\Big{(}S_{L}(\omega)+S_{R}(\omega)\Big{)}G^{*}(\omega)\Gamma_{L}\Big{]},$
(31)

where $G^{*}(\omega)=G(-\omega)$ denotes the complex conjugate of $G(\omega)$.
Proceeding similarly for $J_{2}$, we have from Eq. (23) and Eq. (28),

$\displaystyle
J_{2}=i\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\,\omega\,\text{Tr}[G^{*}(\omega)S_{L}(\omega)].$
(32)

Combining Eqs. (31) and (32) and rearranging the terms, we have,

$\displaystyle
J=\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\,\omega\,\text{Tr}\left[\bigg{(}iG^{*}(\omega)-\omega
G^{*}(\omega)\Gamma_{L}G(\omega)\bigg{)}S_{L}(\omega)\right]-\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\,\omega^{2}\,\text{Tr}\Big{[}G^{*}(\omega)\Gamma_{L}S_{R}(\omega)\Big{]}.$
(33)

Now, remembering the definition of
$G(\omega)=[-M\omega^{2}+\Phi-i\omega(\Gamma_{L}+\Gamma_{R})]^{-1}$, it can be
easily shown that,

$\displaystyle
G^{*}(\omega)\Gamma_{L}G(\omega)=\frac{G(\omega)-G^{*}(\omega)}{2i\omega}-G^{*}(\omega)\Gamma_{R}G(\omega).$
(34)

Using the above relation the first term of Eq. (33) can be further simplified,

$\displaystyle\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\,\omega\,\text{Tr}\left[\bigg{(}iG^{*}(\omega)-\omega
G^{*}(\omega)\Gamma_{L}G(\omega)\bigg{)}S_{L}(\omega)\right]$ (35)
$\displaystyle=$
$\displaystyle\frac{i}{2}\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\,\omega\,\text{Tr}\Big{[}\Big{(}G(\omega)+G^{*}(\omega)\Big{)}S_{L}(\omega)\Big{]}+\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\,\omega^{2}\,\text{Tr}\big{[}G^{*}(\omega)\Gamma_{R}G(\omega)S_{L}(\omega)\big{]}.$
(36)

The first term on the second line vanishes as
$\omega(G(\omega)+G^{*}(\omega))S_{L}(\omega)$ is an odd function of $\omega,$
and we finally have, from Eqs. (33) and (36),

$\displaystyle
J=\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\,\omega^{2}\,\text{Tr}\bigg{[}G(\omega)\,\Gamma_{R}\,G^{*}(\omega)\,S_{L}(\omega)-G(\omega)\,S_{R}(\omega)\,G^{*}(\omega)\,\Gamma_{L}\bigg{]}.$
(37)

From the expressions of $S_{L}(\omega)$ and $S_{R}(\omega)$ given in Eq. (29)
it is immediately clear that $J$ separates into two parts —
$J=J_{\text{th}}+J_{\text{act}}$, where,

$\displaystyle J_{\text{th}}$ $\displaystyle=$
$\displaystyle\gamma^{2}(T_{1}-T_{N})\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\,\omega^{2}|G_{1N}(\omega)|^{2},~{}~{}\text{and}$
(38a) $\displaystyle J_{\text{act}}$ $\displaystyle=$ $\displaystyle
J_{\text{act}}^{1}-J_{\text{act}}^{N}\quad\text{with}\quad
J_{\text{act}}^{j}=\gamma\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\
\omega^{2}|G_{1N}(\omega)|^{2}\tilde{g}(\tau_{j},\omega).$ (38b)

The thermal current $J_{\text{th}}$ is well known in the literature RLL ;
DharReview2008  and is given by,

$\displaystyle
J_{\text{th}}=\frac{k(T_{1}-T_{N})}{2\gamma}\Bigg{[}1+\frac{mk}{2\gamma^{2}}-\frac{mk}{2\gamma^{2}}\sqrt{1+\frac{4\gamma^{2}}{mk}}\Bigg{]}.$
(39)

In the following we compute the active current $J_{\text{act}}$ exactly. To
this end, we first need the explicit form for the matrix element
$G_{1N}(\omega)$. This has been calculated in the context of thermal transport
DharReview2008 , we revisit the calculation here for the sake of completeness.

By definition, $G(\omega)$ is the inverse of a tri-diagonal matrix,

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

The matrix elements $G_{ij}(\omega)$ can be computed explicitly exploiting
this tridiagonal structure of $G^{-1}(\omega)$ Usmani . In particular, we will
need the following elements,

$\displaystyle G_{l1}(\omega)$ $\displaystyle=$
$\displaystyle(-k)^{l-1}\frac{\theta_{N-l}}{\theta_{N}},~{}~{}\text{and}$
(41a) $\displaystyle G_{lN}(\omega)$ $\displaystyle=$
$\displaystyle(-k)^{N-l}\frac{\theta_{l-1}}{\theta_{N}}$ (41b)

where $\theta_{l}$ satisfies the recursion relation,

$\displaystyle\theta_{l}$ $\displaystyle=$
$\displaystyle(-m\omega^{2}+2k)\,\theta_{l-1}-k^{2}\,\theta_{l-2}\qquad\quad\text{for
}l=2,3,\dotsc,N-1,$ (42a) $\displaystyle\text{and}\quad\theta_{N}$
$\displaystyle=$
$\displaystyle(-m\omega^{2}+2k-i\omega\gamma)\,\theta_{N-1}-k^{2}\,\theta_{N-2}.$
(42b)

Using the boundary conditions $\theta_{0}=1$ and
$\theta_{1}=(-m\omega^{2}+2k-i\omega\gamma)$ Usmani , the recursion relation
(42a) can be solved in a straightforward manner. It is convenient to express
the solution as,

$\displaystyle\theta_{l}$ $\displaystyle=$
$\displaystyle\frac{(-k)^{l-1}}{\sin(q)}\left[k\sin\left((l+1)q\right)-i\omega\gamma\sin(lq)\right]\quad\text{for
}l=2,3,\dotsc,N-1,$ (43)

where $q$ and $\omega$ are related by,

$\displaystyle\cos
q=\left(1-\frac{m\omega^{2}}{2k}\right)\quad\Rightarrow\quad\omega=\omega_{c}\,\sin\frac{q}{2},$
(44)

where $\omega_{c}=2\sqrt{k/m}.$ Using Eq. (43) in Eq. (42b) we then have,

$\displaystyle\theta_{N}=\frac{(-k)^{N}}{\sin
q}\,\Big{[}a(q)\sin(Nq)+b(q)\cos(Nq)\Big{]},$ (45)

where,

$\displaystyle a(q)$ $\displaystyle=$
$\displaystyle-\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).$ (46)

Now we can proceed to compute the active current. given by . Using Eq. (45)
and Eq. (41b) for $l=1$ in Eq. (38b), we get,

$\displaystyle J_{\text{act}}^{1}=\frac{\gamma}{\pi
k^{2}}\int_{0}^{\infty}d\omega\,\omega^{2}\frac{\sin^{2}q}{\left|a(q)\sin(Nq)+b(q)\cos(Nq)\right|^{2}}\tilde{g}(\tau_{1},\omega).$
(47)

At this point, it is important to note that, for $\omega>\omega_{c}$, $q$
becomes complex. Thus, for large $N$, in the region $\omega>\omega_{c}$, the
integrand vanishes exponentially as $\exp(-2N\bar{q})$, where $\bar{q}$ is
real. Thus, to compute the current for thermodynamically large systems, we can
limit the range of integration in Eq. (47) to be $0\leq\omega\leq\omega_{c}$
or equivalently, $0\leq q\leq\pi$. Moreover, the functions $\sin(Nq)$ and
$\cos(Nq)$ are highly oscillatory for large $N$ and in the $N\to\infty$ limit,
we can average over $x=Nq$ and write Kannan2012 ,

$\displaystyle J_{\text{act}}^{1}=\frac{\gamma}{\pi
k^{2}}\int_{0}^{\omega_{c}}d\omega\,\omega^{2}\sin^{2}q~{}\,\tilde{g}(\tau_{1},\omega)\int_{0}^{2\pi}\frac{dx}{2\pi}\frac{1}{\left|a(q)\sin
x+b(q)\cos x\right|^{2}}.$ (48)

The $x$-integral has a simple form and can be evaluated exactly,

$\displaystyle\int_{0}^{2\pi}\frac{dx}{2\pi}\frac{1}{(c_{1}\sin x+d\cos
x)^{2}+c_{2}^{2}\sin^{2}x}=-\frac{1}{dc_{2}},$ (49)

where we have denoted $c_{1}=\text{Re}[a(q)]$, $c_{2}=\text{Im}[a(q)]$ and
$d=b(q)$ for notational simplicity. Substituting Eq. (49) in Eq. (48), we get,

$\displaystyle
J_{\text{act}}^{1}=\frac{k}{2\pi}\int_{0}^{\omega_{c}}d\omega\,~{}\,\frac{\omega\sin
q}{k^{2}+\gamma^{2}\omega^{2}}\tilde{g}(\tau_{1},\omega)=\frac{k}{2}\int_{0}^{\pi}\frac{dq}{\pi}\Big{|}\frac{d\omega}{dq}\Big{|}~{}\,\frac{\omega\sin
q}{k^{2}+\gamma^{2}\omega^{2}}\tilde{g}(\tau_{1},\omega).$ (50)

Thereafter, using the Jacobian $|\frac{d\omega}{dq}|=\frac{k\sin q}{m\omega}$,
we arrive at,

$\displaystyle J_{\text{act}}^{1}$ $\displaystyle=$
$\displaystyle\int_{0}^{\pi}\frac{dq}{\pi}\frac{mk\tau_{1}a_{1}^{2}\sin^{2}q}{[mk+2\gamma^{2}(1-\cos
q)][m+2k\tau_{1}^{2}(1-\cos q)]},$ (51)

where we have also expressed $\tilde{g}(\tau_{1},\omega)$ as a function of
$q$. This integral can be evaluated exactly and leads to,

$\displaystyle J_{\text{act}}^{1}$ $\displaystyle=$
$\displaystyle\frac{m}{2\gamma^{2}}\,a_{1}^{2}\mathscr{E}_{1}\quad\text{with}\quad\mathscr{E}_{1}=\frac{\tau_{1}^{2}k^{2}\left[\sqrt{1+\frac{4\gamma^{2}}{mk}}-1\right]+\gamma^{2}\left[1-\sqrt{1+\frac{4k\tau_{1}^{2}}{m}}\right]}{2\tau_{1}(\tau_{1}^{2}k^{2}-\gamma^{2})}.$
(52)

One can similarly obtain
$J_{\text{act}}^{N}=\frac{m}{2\gamma^{2}}\,a_{N}^{2}\mathscr{E}_{N}$, where,

$\displaystyle\mathscr{E}_{N}=\frac{\tau_{N}^{2}k^{2}\left[\sqrt{1+\frac{4\gamma^{2}}{mk}}-1\right]+\gamma^{2}\left[1-\sqrt{1+\frac{4k\tau_{N}^{2}}{m}}\right]}{2\tau_{N}(\tau_{N}^{2}k^{2}-\gamma^{2})}.$
(53)

The total active current is obtained by combining Eq. (53) and (52), which is
quoted in Eq. (5).

##  Appendix B Kinetic temperature profile

The kinetic temperature of the $l^{th}$ oscillator as defined in the main text
is given by,

$\displaystyle\hat{T}_{l}=m\langle\dot{x}_{l}^{2}(t)\rangle.$ (54)

Since we are primarily interested in the effect of the active driving, we put
$T_{1}=T_{N}=0.$ Then, from Eq. (24), we get,

$\displaystyle\hat{T}_{l}=m\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\omega^{2}\Big{[}|G_{l1}(\omega)|^{2}\tilde{g}(\tau_{1},\omega)+|G_{lN}(\omega)|^{2}\tilde{g}(\tau_{N},\omega)\Big{]}.$
(55)

From Eqs. (41) we have,

$\displaystyle|G_{l1}(\omega)|^{2}$ $\displaystyle=$
$\displaystyle\frac{k^{2N-4}}{\sin^{2}q|\theta_{N}|^{2}}\Big{|}k\,\sin(N-l+1)q-i\omega\gamma\,\sin(N-l)q\,\Big{|}^{2},$
(56) $\displaystyle|G_{lN}(\omega)|^{2}$ $\displaystyle=$
$\displaystyle\frac{k^{2N-4}}{2\sin^{2}q|\theta_{N}|^{2}}\Big{|}k\,\sin(lq)-i\omega\gamma\,\sin(l-1)q\,\Big{|}^{2}.$
(57)

We are particularly interested in the behavior of the kinetic temperature at
the bulk in the thermodynamic limit $N\to\infty$. For this purpose we evaluate
$\hat{T}_{l}$ for $l=N/2+\ell$ where $\ell\ll N$. Let us first consider the
contribution from the left reservoir, i.e., the first term in Eq. (55). Once
again, the integrand vanishes exponentially for $\omega>\omega_{c}$ in the
large $N$ limit, and we can write,

$\displaystyle I_{1}$ $\displaystyle\equiv$
$\displaystyle\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\omega^{2}|G_{l1}(\omega)|^{2}\tilde{g}(\tau_{1},\omega)$
(58) $\displaystyle=$
$\displaystyle\frac{1}{k^{4}}\int_{0}^{\pi}\frac{dq}{2\pi}\left|\frac{d\omega}{dq}\right|\omega^{2}\,\frac{k^{2}\,(1-\cos(N-2\ell+2)q)+\omega^{2}\gamma^{2}\,(1-\cos(N-2\ell)q)}{|a(q)\sin(Nq)+b(q)\sin(Nq)|^{2}}\tilde{g}(\tau_{1},\omega).$
(59)

As before, in the $N\to\infty$ limit, we can average over the fast
oscillations in $x=Nq$. For this purpose, let us note,

$\displaystyle\int_{0}^{2\pi}\frac{dx}{2\pi}\frac{\sin x}{(c_{1}\sin x+d\cos
x)^{2}+c_{2}^{2}\sin^{2}x}=\int_{0}^{2\pi}\frac{dx}{2\pi}\frac{\cos
x}{(c_{1}\sin x+d\cos x)^{2}+c_{2}^{2}\sin^{2}x}=0.$ (60)

Using these identities and Eq. (49), Eq. (59) reduces to,

$\displaystyle I_{1}$ $\displaystyle=$ $\displaystyle\frac{1}{2\gamma
k}\int_{0}^{\pi}\frac{dq}{2\pi}\Big{|}\frac{d\omega}{dq}\Big{|}\,\frac{\omega}{\sin
q}\tilde{g}(\tau_{1},\omega)=\frac{1}{4\pi\gamma
m}\int_{0}^{\pi}dq~{}\tilde{g}(\tau_{1},\omega(q))$ (61)

The $q$-integral can be evaluated exactly, and yields,

$\displaystyle I_{1}=\frac{1}{2\gamma
m}\frac{a_{1}^{2}\tau_{1}}{\sqrt{1+4\tau_{1}^{2}k/m}}.$ (62)

The integral involving $G_{lN}$ can also be performed following the same
procedure and results in,

$\displaystyle
I_{2}\equiv\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\omega^{2}|G_{lN}(\omega)|^{2}\tilde{g}(\tau_{N},\omega)=\frac{1}{2\gamma
m}\frac{a_{N}^{2}\tau_{N}}{\sqrt{1+4\tau_{N}^{2}k/m}}.$ (63)

Combining these results, we see that the kinetic temperature remains uniform
at the bulk and is given by,

$\displaystyle\hat{T}_{bulk}$ $\displaystyle=$
$\displaystyle\frac{1}{2\gamma}\left(\frac{a_{1}^{2}\tau_{1}}{\sqrt{1+4\tau_{1}^{2}k/m}}+\frac{a_{N}^{2}\tau_{N}}{\sqrt{1+4\tau_{N}^{2}k/m}}\right).$
(64)

This is the result presented in Eq. (6).

Figure 6: Boundary layer properties of the kinetic temperature profile: (a)
shows $\hat{T}_{l}$ profile near the left boundary. The main plot compares the
contributions from Eq. (68), (73) and (77) (in solid black lines) with
numerical simulations (in colored symbols). The inset plot shows the
exponential decay of $\hat{T}_{l}$ from $\hat{T}_{\text{bulk}}$ at the left
boundary for a system of $N=64$ oscillators with $m=1,\,k=1,\,\gamma=1$ and
$\tau_{1}=1$. (b) shows the absence (presence) of the boundary kinks when
$\tau_{1,N}$ is much larger (smaller) than $\omega_{c}=2\sqrt{k/m}$. Here
$N=32$ with $m=1,\,k=0.5,\,\gamma=1$ and $\tau_{1}=1$.

For a finite chain the kinetic temperature deviates from $\hat{T}_{bulk}$ near
the boundaries giving rise to exponentially decaying boundary layers; see Fig.
6(a). To obtain the behavior of the boundary layers, we need to evaluate Eq.
(55) in the limits $l\ll N$ and $l\sim N$.

###  B.1 $\hat{T}_{l}$ near left boundary

Let us first concentrate near the left boundary, where $l=1,2,3\dotsc\ll N$.
For convenience, we rewrite Eq. (55) as,

$\displaystyle\hat{T}_{l}=m[L_{1}(l,\tau_{1})+L_{N}(l,\tau_{N})],$ (65)

where $L_{1}(l,\tau_{1})$ and $L_{N}(l,\tau_{N})$ denote the contributions
from the left and right reservoirs respectively,

$\displaystyle L_{1}(l,\tau)$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\omega^{2}|G_{l1}(\omega)|^{2}\tilde{g}(\tau,\omega)$
(66a) $\displaystyle L_{N}(l,\tau)$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\omega^{2}|G_{lN}(\omega)|^{2}\tilde{g}(\tau,\omega)$
(66b)

We first evaluate the contribution from the right reservoir $L_{N}(l,\tau)$.
In this case, once again, the contribution coming from $|\omega|>\omega_{c}$
vanishes exponentially for large $N$ and in the thermodynamic limit Eq. (66b)
reduces to,

$\displaystyle
L_{N}(l,\tau)=\frac{1}{k^{4}}\int_{0}^{\omega_{c}}\frac{d\omega}{\pi}\omega^{2}\,\frac{k^{2}\sin^{2}(lq)+\omega^{2}\gamma^{2}\sin^{2}(l-1)q}{|a(q)\sin(Nq)+b(q)\sin(Nq)|^{2}}\,\tilde{g}(\tau,\omega)$
(67)

Averaging over the fast oscillations in the $N\to\infty$ limit and using Eq.
(49), we get,

$\displaystyle L_{N}(l,\tau)=\frac{1}{2\gamma
m}\int_{0}^{\pi}\frac{dq}{\pi}\,\frac{k^{2}\sin^{2}(lq)+\omega^{2}\gamma^{2}\sin^{2}(l-1)q}{(k^{2}+\gamma^{2}\omega^{2})}\,\,\tilde{g}(\tau,\omega).$
(68)

Though this integral does not yield any closed form expression, it can be
evaluated numerically for arbitrary $l$ and $\tau$.

Next, we consider the contribution from the left reservoir $L_{1}(l,\tau)$. It
turns out that $L_{1}(l,\tau)$ (Eq. (66a)) has non-vanishing contribution from
both $|\omega|<\omega_{c}$ and $|\omega|>\omega_{c}$. Thus, it is convenient
to rewrite Eq. (66a) as,

$\displaystyle L_{1}(l,\tau)=L_{1}^{b}(l,\tau)+L_{1}^{o}(l,\tau),$ (69)

where $L_{1}^{b}(l,\tau)$ and $L_{1}^{o}(l,\tau)$ denote the contributions
from $|\omega|<\omega_{c}$ and $|\omega|>\omega_{c}$ respectively. For
$|\omega|>\omega_{c}$, Eq. (44) implies that $q=\pi-i\bar{q}$, where $\bar{q}$
is real. We first evaluate the contribution from this region,

$\displaystyle
L_{1}^{o}(\tau)=\int_{\omega_{c}}^{\infty}\frac{d\omega}{\pi}\omega^{2}|G_{l1}|^{2}\,\tilde{g}(\tau,\omega)=\frac{1}{k^{4}}\int_{\omega_{c}}^{\infty}\frac{d\omega}{\pi}\omega^{2}\frac{|ik\sinh(N-l+1)\bar{q}-\omega\gamma\sinh(N-l)\bar{q}|^{2}}{|ia(\bar{q})\sinh(Nq)-b(\bar{q})\cosh(N\bar{q})|^{2}}\,\tilde{g}(\tau,\omega),$
(70)

where we have used the identities,

$\displaystyle\sin(nq)=(-1)^{n+1}~{}i~{}\sinh(n\bar{q}),\quad\text{and}\quad\cos(nq)=(-1)^{n}\cosh(n\bar{q}),\qquad
n=0,1,2,\dotsc.$ (71)

In the $N\to\infty$ limit, Eq. (70) reduces to

$\displaystyle
L_{1}^{o}(\tau)=\frac{1}{k^{4}}\int_{\omega_{c}}^{\infty}\frac{d\omega}{\pi}\omega^{2}e^{-2l\bar{q}}\frac{k^{2}e^{2\bar{q}}+\omega^{2}\gamma^{2}}{|ia(\bar{q})-b(\bar{q})|^{2}}\tilde{g}(\tau,\omega)$
(72)

The integral over $\omega\in[\omega_{c},\infty]$ can be converted to an
integral over $\bar{q}\in[0,\infty]$ using the relation
$\omega=\omega_{c}\cosh(\bar{q}/2)$ [see Eq. (44)], to get,

$\displaystyle
L_{1}^{o}(\tau)=\frac{\omega_{c}^{3}}{2}\int_{0}^{\infty}\frac{d\bar{q}}{\pi}\,e^{-2l\bar{q}}\,\frac{\sinh(\bar{q}/2)\cosh^{2}(\bar{q}/2)}{k^{2}+\gamma^{2}\omega_{c}^{2}\cosh^{2}(\bar{q}/2)e^{-2\bar{q}}}\tilde{g}\Big{(}\tau,\omega_{c}\cosh\frac{\bar{q}}{2}\Big{)}.$
(73)

This, again, can be evaluated numerically for arbitrary $l$. For
$|\omega|<\omega_{c}$, $q$ is real and the contribution to Eq. (66a) is given
by,

$\displaystyle
L_{1}^{b}(\tau)=\int_{0}^{\omega_{c}}\frac{d\omega}{\pi}\omega^{2}|G_{l1}|^{2}\,\tilde{g}(\tau,\omega)=\frac{1}{k^{4}}\int_{0}^{\omega_{c}}\frac{d\omega}{\pi}\omega^{2}\,\frac{k^{2}\sin^{2}(N-l+1)q+\omega^{2}\gamma^{2}\sin^{2}(N-l)q}{|a(q)\sin(Nq)+b(q)\sin(Nq)|^{2}}\,\tilde{g}(\tau,\omega)$
(74)

For $l\ll N$ and $N\to\infty$ limit, averaging over the fast oscillations
$x=Nq$ involves integrals of the form,

$\displaystyle
Q(v)=\int_{0}^{2\pi}\frac{dx}{2\pi}\frac{\sin^{2}(x-v)}{(c_{1}\sin x+d\cos
x)^{2}+c_{2}^{2}\sin^{2}x}=\frac{(c_{1}^{2}+c_{2}^{2}-d^{2})\cos(2v)-(c_{1}^{2}+(c_{2}-d)^{2}+2c_{1}d\sin(2v))}{2dc_{2}(c_{1}^{2}+(c_{2}-d)^{2})},$
(75)

where $v>0$ is arbitrary and $c_{1}=\text{Re}[a(q)]$, $c_{2}=\text{Im}[a(q)]$
and $d=b(q)$ as before. Using the above result in Eq. (74) with appropriate
values of $v$,

$\displaystyle L_{1}^{b}(\tau)$ $\displaystyle=$
$\displaystyle\frac{1}{k^{4}}\int_{0}^{\omega_{c}}\frac{d\omega}{\pi}\omega^{2}\Big{[}k^{2}\,Q((l-1)q)+\omega^{2}\gamma^{2}\,Q(lq)\Big{]}\tilde{g}(\tau,\omega)$
(76) $\displaystyle=$
$\displaystyle\frac{\omega_{c}^{3}}{2k^{4}}\int_{0}^{\infty}\frac{d\bar{q}}{\pi}\sin^{2}\frac{q}{2}\cos\frac{q}{2}\Big{[}k^{2}\,Q((l-1)q)+\omega_{c}^{2}\sin^{2}\frac{q}{2}\gamma^{2}\,Q(lq)\Big{]}\tilde{g}\Big{(}\tau,\omega_{c}\sin\frac{q}{2}\Big{)}$
(77)

Adding the contributions given by Eq. (68), (73) and (77), we can evaluate the
kinetic temperature profile near the left boundary, which is shown in Fig. 6.

###  B.2 $\hat{T}_{l}$ near right boundary

The behavior near the right boundary can be obtained in a similar manner. For
this purpose, it is convenient to define, $\ell=N-l+1$, such that
$\ell=1,2,3,\dotsc\ll N$ corresponds to the oscillators near the right
boundary. Next, we note that, from Eqs. (56) and (57),

$\displaystyle|G_{l1}(\omega)|^{2}=|G_{N-l+1,N}(\omega)|^{2}$ (78)
$\displaystyle|G_{lN}(\omega)|^{2}=|G_{N-l+1,1}(\omega)|^{2}.$ (79)

Then the $\hat{T}_{l}$ profile near the right boundary ($l\sim N$) is given
by,

$\displaystyle\hat{T}_{N-\ell+1}$ $\displaystyle=$ $\displaystyle
m\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\omega^{2}\Big{[}|G_{{N-\ell+1},1}(\omega)|^{2}\tilde{g}(\tau_{1},\omega)+|G_{{N-\ell+1,N}}(\omega)|^{2}\tilde{g}(\tau_{N},\omega)\Big{]}$
(80) $\displaystyle=$ $\displaystyle
m\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\omega^{2}\Big{[}|G_{\ell,N}(\omega)|^{2}\tilde{g}(\tau_{1},\omega)+|G_{\ell,1}(\omega)|^{2}\tilde{g}(\tau_{N},\omega)\Big{]}$
(81) $\displaystyle=$ $\displaystyle
m\Big{[}L_{N}(\ell,\tau_{1})+L_{1}(\ell,\tau_{N})\Big{]},$ (82)

where $L_{1}(\ell,\tau)$ and $L_{N}(\ell,\tau)$ are obtained from Eqs. (68),
(73) and (77).

Interestingly, boundary kinks, which are absent in the active regime appear in
the passive limit, similar to the thermal scenario. This is shown in Fig. 6(b)
where kinks are visible near the right boundary as the activity of the
corresponding reservoirs is small, whereas no kinks are visible near the left
reservoir, which remains in the strongly active regime.

##  Appendix C Linear response: Differential conductivity

In this section we derive an expression for the differential conductivity for
the energy current using nonequilibrium response theory. Linear response
relations in nonequilibrium are most conveniently derived using a trajectory
based approach Maes2020Front , and we take the the specific example of
dichotomous active forces which flips sign with rates $\alpha_{1}$ and
$\alpha_{N}$ at the two reservoirs [see Eq. (3)]. In this case, it is most
natural to consider a perturbation $\alpha_{j}\to\alpha_{j}+d\alpha_{j}$ and
express the differential conductivity as,

$\displaystyle\chi_{j}\equiv\frac{dJ_{\text{act}}}{d\tau_{j}}=-\frac{1}{2\tau_{j}^{2}}\frac{dJ_{\text{act}}}{d\alpha_{j}},$
(83)

where we have used the fact that $\tau_{j}=1/(2\alpha_{j})$ in this scenario.
Let $\omega=\\{x_{i}(s),\sigma_{1}(s),\sigma_{N}(s);0\leq s\leq t\\}$ denote a
trajectory of the system during the interval $[0,t]$ and
$P_{\alpha_{1},\alpha_{N}}(\omega)$ denote the corresponding probability. Of
course, the trajectory probability depends on the various system parameters,
but since we are interested in the response to a change in the flip rate, it
suffices to consider the $\alpha_{j}$ dependence. Hence, we can write,

$\displaystyle
P_{\alpha_{1},\alpha_{N}}(\omega)=e^{-(\alpha_{1}+\alpha_{N})t}\alpha_{1}^{n_{1}}\alpha_{N}^{n_{N}}{\cal
U}(\omega)$ (84)

where $n_{1}$ and $n_{N}$ denote the number of flips of the active force at
the left and right boundaries during time $[0,t]$ and ${\cal U}(\omega)$
contains the $x_{j}$-dependent components. The weight of the same trajectory
$\omega$ changes upon adding the perturbation i.e., say, changing
$\alpha_{1}\to\alpha_{1}+d\alpha_{1}$. Note that this change does not affect
${\cal U}(\omega)$. The linear response of the expectation value of any
observable $\langle O\rangle$ to this change can be expressed as a connected
correlation in the unperturbed state Maes2020Front ,

$\displaystyle\frac{d\langle O\rangle}{d\alpha_{1}}=-\langle{\cal
A}(\omega);O(\omega)\rangle$ (85)

where ${\cal A}(\omega)=-\frac{d}{d\alpha_{1}}\log
P_{\alpha_{1},\alpha_{N}}(\omega)$ is the excess action associated to the
trajectory $\omega$ due to the perturbation.

Figure 7: Plot of $\frac{dJ_{\text{act}}}{d\alpha_{1}}$ versus $\alpha_{1}$
for the dichotomous active noise for two different values of $k$: the solid
lines show the exact expression obtained from Eq. (52) while the symbols show
the predicted correlation [see Eq. (86)] measured from numerical simulations.

Then, from Eq. (84), we have, for the active current
$J_{\text{act}}=\langle\mathscr{J}(t)\rangle$,

$\displaystyle\frac{dJ_{\text{act}}}{d\alpha_{1}}=\frac{1}{\alpha_{1}}[\langle
n_{1}\mathscr{J}(t)\rangle-\langle
n_{1}\rangle\langle\mathscr{J}(t)\rangle]].$ (86)

The response to a change in the right reservoir is also given by a similar
expression. The stationary response is obtained by taking the $t\to\infty$
limit which is quoted in Eq. (15), in terms of the activity parameter
$\tau_{j}=1/(2\alpha_{j})$. Figure 7 compares the exact analytical response
$\frac{dJ_{\text{act}}}{d\alpha_{1}}$ obtained from Eqs. (52)-(53) using
$\alpha_{1}=1/(2\tau_{1})$ with the prediction (86), measured from numerical
simulations, which shows an excellent match.

We close this discussion with a final remark. The nonequilibrium linear
response relation obtained in Eq. (86) is purely frenetic if the active noise
$f_{j}=a_{j}\sigma_{j}$ is considered as a force, and hence symmetric under
time-reversal. Frenetic contribution to the linear response is known to result
in negative differential response in various contexts ndc_2013 . The activity
driven harmonic chain provides another example where the same mechanism works,
although the absence of any equilibrium limit and the nature of the
perturbation here means that there is no traditional regime where one recovers
a Kubo-like formula.

## References

  * (1) Thermal Transport in Low Dimensions, Ed. Stefano Lepri, Springer Heidelberg (2016).
  * (2) A. Dhar, Adv. in Phys., 57, 457 (2008). 10.1080/00018730802538522
  * (3) Z. Rieder, J. L. Lebowitz, and E. Lieb, J. Math. Phys. 8, 1073 (1967). 10.1007/978-3-662-10018-9_21
  * (4) H. Nakazawa, Prog. Theor. Phys. Suppl. 45, 231 (1970). 10.1143/PTPS.45.231
  * (5) D. Roy and A. Dhar, J. Stat. Phys. 131, 535 (2008). 10.1007/s10955-008-9487-1
  * (6) A. Dhar, Phys. Rev. Lett. 86, 5882 (2001). 10.1103/PhysRevLett.86.5882
  * (7) S. Lepri, R. Livi, and A. Politi, Chaos 15, 015118 (2005). 10.1063/1.1854281
  * (8) T. Mai, A. Dhar, and O. Narayan, Phys. Rev. Lett. 98, 184301 (2007). 10.1103/PhysRevLett.98.184301
  * (9) A. Kundu, S. Sabhapandit, and A. Dhar, J. Stat. Mech. 2011, P03007 (2011). 10.1088/1742-5468/2011/03/P03007
  * (10) V. Kannan, A. Dhar, and J. L. Lebowitz, Phys. Rev. E 85, 041118 (2012). 10.1103/PhysRevE.85.041118
  * (11) D. Gupta and D. A. Sivak, Phys. Rev. E 104, 024605 (2021). 10.1103/PhysRevE.104.024605
  * (12) R. Kubo, Rep. Prog. Phys. 29, 255 (1966). 10.1088/0034-4885/29/1/306
  * (13) C. Maes and S. R. Thomas, Phys. Rev. E 87, 022145 (2013). 10.1103/PhysRevE.87.022145
  * (14) C. Maes, J. Stat. Phys. 154, 705 (2014). 10.1007/s10955-013-0904-8
  * (15) C. Maes and S. Steffenoni, Phys. Rev. E 91, 022128 (2015). 10.1103/PhysRevE.91.022128
  * (16) H. Vandebroek and C. Vanderzande, J. Stat. Phys. 167, 14, (2017). 10.1007/s10955-017-1734-x
  * (17) D. Gupta and S. Sabhapandit, Phys. Rev. E 96, 042130 (2017). 10.1103/PhysRevE.96.042130
  * (18) A. Iacobucci, F. Legoll, S. Olla, and G. Stoltz, Phys. Rev. E 84, 061108 (2011). 10.1103/PhysRevE.84.061108
  * (19) M. C. Zheng, F. M. Ellis, T. Kottos, R. Fleischmann, T. Geisel, and T. Prosen, Phys. Rev. E 84, 021119 (2011). 10.1103/PhysRevE.84.021119
  * (20) D. Bagchi, J. Stat. Mech. 2013, P12005 (2013). 10.1088/1742-5468/2013/12/P12005
  * (21) K. Kanazawa, T. Sagawa, and H. Hayakawa, Phys. Rev. E 87, 052124, (2013). 10.1103/PhysRevE.87.052124
  * (22) X. L. Wu and A. Libchaber, Phys. Rev. Lett. 84, 3017 (2000). 10.1103/PhysRevLett.84.3017
  * (23) G. V. Soni, B. M. J. Ali, Y. Hatwalne, and G. Shivshankar, Biophys. J. 84, 2634 (2003). 10.1016/s0006-3495(03)75068-1
  * (24) S. Krishnamurthy, S. Ghosh, D. Chatterji, R. Ganapathy, and A. K. Sood, Nature Phys 12, 1134 (2016). 10.1038/nphys3870
  * (25) C. Valeriani, M. Li, J. Novosel, J. Arlta and D. Marenduzzoa, Soft Matter, 7, 5228 (2011). 10.1039/C1SM05260H
  * (26) C. Maggi, M. Paoluzzi, N. Pellicciotta, A. Lepore, L. Angelani, and R. Di Leonardo, Phys. Rev. Lett. 113, 238303 (2014). 10.1103/PhysRevLett.113.238303
  * (27) A. Gopal, É. Roldán, and S. Ruffo, J. Phys. A: Math. Theor. 54, 164001 (2021). 10.1088/1751-8121/abe5cb
  * (28) C. Maes, Phys. Rev. Lett. 125, 208001 (2020). 10.1103/PhysRevLett.125.208001
  * (29) O. Granek, Y. Kafri and J. Tailleur, arXiv:2108.11970. 10.48550/arXiv.2108.11970
  * (30) H. Seyforth, M. Gomez, W. B. Rogers, J. L. Ross, W. W. Ahmed, Phys. Rev. Research 4, 023043 (2022). 10.1103/PhysRevResearch.4.023043
  * (31) S. M. Mousavi, G. Gompper, and R. G. Winkler, J. Chem. Phys. 155, 044902 (2021). 10.1063/5.0058150
  * (32) A. Pal and S. Sabhapandit, Phys. Rev. E 90, 052116 (2014). doi.org/10.1103/PhysRevE.90.052116
  * (33) É. Fodor, T. Nemoto, and S. Vaikuntanathan, New J. Phys. 22, 013052 (2020). 10.1088/1367-2630/ab6353
  * (34) S. Chaki and R. Chakrabarti, Physica A, 530, 121574 (2019). 10.1016/j.physa.2019.121574
  * (35) P. Baerts, U. Basu, C. Maes, S. Safaverdi, Phys. Rev. E 88, 052109 (2013). 10.1103/PhysRevE.88.052109
  * (36) B. Li, L. Wang, and G. Casati, Appl. Phys. Lett. 88, 143501 (2006). 10.1063/1.2191730
  * (37) S. Leitmann and T. Franosch, Phys. Rev. Lett. 111, 190603 (2013). 10.1103/PhysRevLett.111.190603
  * (38) A. K. Chatterjee, U. Basu and P. K. Mohanty, Phys. Rev. E 97, 052137 (2018). 10.1103/PhysRevE.97.052137
  * (39) A. K. Chatterjee and P. K. Mohanty, Phys. Rev. E 98, 062134 (2018). 10.1103/PhysRevE.98.062134
  * (40) D. Martin, J. O’Byrne, M. E. Cates, É. Fodor, C. Nardini, J. Tailleur, F. van Wijland, Phys. Rev. E 103, 032607 (2021). 10.1103/PhysRevE.103.032607
  * (41) I. Santra, U. Basu, S. Sabhapandit, Phys. Rev. E 101, 062120 (2020). 10.1103/PhysRevE.101.062120
  * (42) U. Basu, S. N. Majumdar, A. Rosso, G. Schehr, Phys. Rev. E 98, 062121 (2018). 10.1103/PhysRevE.98.062121
  * (43) I. Santra, U. Basu, S. Sabhapandit, Phys. Rev. E 104, L012601 (2021). 10.1103/PhysRevE.104.L012601
  * (44) R. Usmani, Comput. Math. Appl., 27, 59 (1994). 10.1016/0898-1221(94)90066-3
  * (45) M. Baiesi, C. Maes, and B. Wynants, Phys. Rev. Lett. 103, 010602 (2009). 10.1103/PhysRevLett.103.010602
  * (46) C. Maes, Front. Phys. 8, 229, (2020). 10.3389/fphy.2020.00229 Focus to learn more
  * (47) S. Lepri, R. Livi, and A. Politi, Phys. Rep. 377, 1 (2003). 10.1016/S0370-1573Focus to learn more
  * (48) L. Caprini and U. M. B. Marconi, Phys. Rev. Research 2, 033518 (2020). 10.1103/PhysRevResearch.2.033518
  * (49) P. Singh and A. Kundu, J. Phys. A: Math. Theor. 54, 305001 (2021). 10.1088/1751-8121/ac0a9f