# Microscopic modifications of equilibrium probabilities due to non-
conservative perturbations with applications to anharmonic systems

Dino Osmanović  Department of Mechanical and Aeronautical Engineering,
University of California, Los Angeles, Los Angeles, CA, USA
<EMAIL_ADDRESS>
###### Abstract

The standard relationships of statistical mechanics are upended my the
presence of active forces. In particular, it is no longer usually possible to
simply write down what the stationary probability of a state of such a system
will be, as can be done in ordinary statistical mechanics. While exact
expressions are possible for harmonic systems, anharmonic systems tend to be
much more difficult to treat. In this manuscript, we investigate how the
microscopic probability for anharmonic systems is modified in the presence of
non-conservative forces. We recount how non-conservative corrections to
microscopic probabilities in generic systems can be represented as integrals
over Green’s function kernels, and that these Green’s functions take the form
of path integrals. We show how using analytically tractable form of these
functions allows us to calculate corrections to the microscopic probability
for a generic anharmonic system. These results are compared to active Brownian
particles inside a harmonic or anharmonic well. From this it can be shown that
accumulation of probability away from its equilibrium minima is a natural
effect arising in anharmonic but not harmonic systems. Finally, we extend the
microscopic results to study small active polymers with anharmonic backbone
potentials. The interplay of the active driving with the anharmonic backbone
potential can be used to understand how the end-to-end distance and mode space
distributions of anharmonic active polymers scales with active driving. In
particular, the polymer modes with the smallest eigenvalues will be affected
by non-conservative forces the most.

##  I Introduction

Active systems constitute a class of physical models whose continual
consumption of energy requires the use of new theoretical methods and ideas to
describe their behavior[1]. The presence of energy-consuming processes leads
to violations of the conditions necessary for the establishment of equilibrium
states. Such systems are of interest as many phenomena belong to this domain;
most importantly, this is a defining characteristic of biological systems[2].

Analysis of problems in these areas would be deemed to be related to non-
equilibrium physics. However, the analysis of non-equilibrium systems is
complicated by the fact that ”non-equilibrium” denotes a variety of different
systems, each of which has its own peculiarities. For example, glasses are not
in equilibrium[3], though intuitively appear as a different “kind” of non-
equilibrium to a living cell. Other examples of non-equilibrium could include
maintaining two different parts of a system at two different temperatures,
leading to continuous heat transfer from the hotter part of the system to the
cooler one[4]. Lastly, “classical” active matter (as described above) is often
thought of as using consumption of energy for directed motion[5], though its
definition can also include consumption of energy to drive non-directed
processes[6].

Focusing on problems which have relevance for biology, the lack of any
unifying description to such systems has necessitated a variety of methods for
describing the non-equilibrium stationary states in active systems. A variety
of active problems are well described through field theoretic methods, where
analysis of the symmetries of the problem allow construction of functionals
that can give mean-field behavior [7]. Alternatively, one can describe all the
microscopic degrees of freedom, however, such a level of description usually
requires coarse grained simulations, given the lack of useful expressions for
the microscopic distributions over states at non-equilibrium, nevertheless, a
variety of different systems can be studied in such a manner, including
flocking behaviors, active Brownian motion and active chemistry[8, 9, 10].
More abstractly, stochastic thermodynamics generalizes standard thermodynamic
notions towards the analysis of non-equilibrium steady states[11], which can
be extended towards the study of active matter [12].

Common to all the above is that the standard Boltzmann distribution over
states is no longer sufficient to describe these scenarios[13]. Notably,
however, occasionally such scenarios can be described rather simply through
renormalization of temperature (so called active (or effective)
temperature)[14, 15, 16, 17, 18]. One is left with the question of whether any
generalized microscopic description can tell us something about the properties
of non-equilibrium steady states.

One way to explore this issue is by treating excess forces in these systems as
a perturbation[19, 20, 21, 22, 23, 24]. The perturbations can be represented
as integrals over Green’s function kernels[25], however, the precise form of
these kernels is difficult to give in analytic form, other than as the
solution to a complicated partial differential equation. This Green’s function
governs the response of the microstate probability distribution to additional
(not necessarily arising from a potential) forces in the system, leading to a
new probability of each microstate. As deviations from Boltzmann distributions
can be characterized in terms of integrals over this Green’s function kernel,
which serves then as a microscopic response function, what can we say about
this function that governs microscopic responses to non-equilibrium
perturbations? Response functions are a canonical way of treating non-
equilibrium states[26, 27], however they are usually thought of in terms of
macroscopic rather than microscopic quantities. We note here the caveat that
we use the term “response” rather generally and in a way that is not
necessarily related to its classical usage in fluctuation-dissipation
theories. We are interested in the function that tells us how the microscopic
distribution of states is modified through the presence of non-equilibrium
perturbations, as opposed to the dynamical features of the system in the
presence of an additional time dependent force.

Prior works have demonstrated interpretable analytic expressions for harmonic
systems[29] (systems where the forces are linear in the degrees of freedom)
subject to dissipative driving. Anharmonic systems have attracted some
analytical attention[30, 31, 32], though less than harmonic systems, possibly
due to the general difficulty involved in treating anharmonic forces.

In this manuscript, we set out to explore how this microscopic distribution
over states for non-conservative perturbation can be calculated, for a generic
system including anharmonic conservative forces. To this end, we first
demonstrate that the Green’s function can be written as a path integral over
all paths in the equilibrium system and give an interpretation of the
expression that arises. The similarity of this path integral to those that
arise in quantum theory allows us to derive analytical results for anharmonic
systems. We derive these corrections for an active Brownian particle inside an
anharmonic trap and compare these results against simulations of the same. We
demonstrate that approximate forms of the Green’s function kernel can generate
correct results even up to quite large values of the activity. We then apply
the microscopic expressions that arise from the study of the anharmonic trap
to an active polymer with an anharmonic backbone potential.

##  II Results

###  II.1 Microscopic stationary Green’s functions

We imagine an equilibrium system by positing a Hamiltonian $H(\mathbf{x})$
describing $N$ particles which are also coupled to an external thermal bath at
an inverse temperature $\beta=1/k_{b}T$. At this juncture we introduce a point
of nomenclature to help ease confusion. While Equilibrium can be said to
correspond to a set of states of a system with Hamiltonian $H$, such a system
could be instantiated in a particular region of phase space that is not in
this “equilibrium” subset. However, it is (given sufficient time) still
possible for such a system to reach equilibrium states. This differs in kind
to a system where not all the forces arise from the Hamiltonian, which we
refer to as “non-equilibrium”.

For simplicity, we consider the case where each of these particles are
overdamped (i.e. a colloidal system), such that a state contains only the
position degrees of freedom of all the particles. In standard statistical
mechanics, the probability of a state $\mathbf{x}$ is proportional to the
Boltzmann factor $P(\mathbf{x})\sim\exp(-\beta H(\mathbf{x}))$. However, in
the presence of additional forces which cannot be described through potentials
$\mathbf{F}_{nc}(\mathbf{x})$, this probability doesn’t describe the true
probability of a state probability of the system. In a previous
manuscript[25], we demonstrated through the Smoluchowski equation that the
effect of a non-conservative force (to first order) on the microscopic
probability appears as:

$\log(P(\mathbf{x}))\approx-\beta
H(\mathbf{x})+\int\mathrm{d}\mathbf{x}^{\prime}G(\mathbf{x},\mathbf{x}^{\prime})\rho(\mathbf{F}(\mathbf{x}^{\prime}))$
(1)

where we have suppressed additional constants that would appear on the right
hand side. One could also assert (1) quite generally from linear response
theory. The Green’s function $G(\mathbf{x}^{\prime},\mathbf{x})$ quantifies
how a perturbation at a point $\mathbf{x}^{\prime}$ affects some quantity at
point $\mathbf{x}$. The “source function” $\rho(\mathbf{F}(\mathbf{x}))$ is a
function of the non-conservative perturbation on the system, but we shall see
that it’s generally related to the total power dissipated by the process that
is being introduced. The precise forms of these functions can be calculated
for whatever form of non-conservative perturbation is present. For example, an
equilibrium system subject to some divergence free non-conservative force
$\mathbf{F}_{\text{nc}}(\mathbf{x})$ has the following source function:

$\rho(\mathbf{x})=\beta^{2}\nabla
H(\mathbf{x}).\mathbf{F}_{\text{nc}}(\mathbf{x})$ (2)

Higher order perturbations (see below) are also integrated over the kernel
$G(\mathbf{x},\mathbf{x}^{\prime})$, reinforcing its importance in governing
how the probability of a state depends on non-conservative forces. The Green’s
function is also sometimes referred to as a response function. We again note
the caveat stated in the introduction over the term “response”. This function
has no time dependence, but tell us how the probability at a point
$\mathbf{x}$ “responds” to a delta function impulse at another point
$\mathbf{x}^{\prime}$. We can make this clearer by referring to the function
as the “stationary response function”.

It can be demonstrated via analysis of the Smoluchowski equation[25] that the
Green’s function for perturbation by a non-conservative force is given by the
solution to the following equation:

$\left(\nabla^{2}-\beta\nabla
H(\mathbf{x}).\nabla\right)G(\mathbf{x},\mathbf{x}^{\prime})=\delta(\mathbf{x}-\mathbf{x}^{\prime})$
(3)

which applies regardless of the type of perturbation involves coupling to
different temperature baths or action by deterministic forces. (See appendix A
for a derivation of this equation for a multi-bath system)

As was mentioned previously, the Green’s function characterizes how the
microscopic probability is modified through the presence of non-potential
forces. However, equation 3 lacks a very clear intuitive physical
interpretation of the important elements of this function. Remarkably,
however, equation 3 can be solved by a path integral, a derivation which arose
in a rather different context (electrostatics)[33] but appears nevertheless to
apply to non-equilibrium systems.

The function $G(\mathbf{x},\mathbf{x}^{\prime})$ can be written as an integral
over all paths $\mathbf{r}(s)$:

$\displaystyle G(\mathbf{x},\mathbf{x}^{\prime})$
$\displaystyle=\int_{0}^{\infty}\\!\\!\\!\\!\\!\\!\mathrm{d}s\int_{\mathbf{r}(0)=\mathbf{x}^{\prime}}^{\mathbf{r}(s)=\mathbf{x}}\\!\\!\\!\\!\\!\\!\mathcal{D}\left[\mathbf{r}(s)\right]e^{-\int_{0}^{s}\mathrm{d}s^{\prime}\frac{1}{4}\left(\frac{d\mathbf{r}(s^{\prime})}{ds^{\prime}}\right)^{2}-\frac{1}{2}\beta\nabla
H(\mathbf{r}(s^{\prime})).\frac{d\mathbf{r}(s^{\prime})}{ds^{\prime}}+\frac{1}{4}\beta^{2}\nabla
H(\mathbf{r}(s^{\prime})).\nabla
H(\mathbf{r}(s^{\prime}))-\frac{1}{2}\beta\nabla^{2}H(\mathbf{r}(s^{\prime}))}$
(4)
$\displaystyle=\int_{0}^{\infty}\\!\\!\\!\\!\\!\\!\mathrm{d}s\int_{\mathbf{r}(0)=\mathbf{x}^{\prime}}^{\mathbf{r}(s)=\mathbf{x}}\\!\\!\\!\\!\\!\\!\mathcal{D}\left[\mathbf{r}(s)\right]e^{-\int_{0}^{s}\mathrm{d}s^{\prime}L\left(\frac{d\mathbf{r}(s^{\prime})}{ds^{\prime}},\mathbf{r}(s^{\prime}),s^{\prime}\right)}$
(5)

Which can be rewritten as:

$\displaystyle G(\mathbf{x},\mathbf{x}^{\prime})$
$\displaystyle=e^{\frac{\beta}{2}\left(H(x)-H(x^{\prime})\right)}K(\mathbf{x},\mathbf{x}^{\prime})$
(6) $\displaystyle K(\mathbf{x},\mathbf{x}^{\prime})$
$\displaystyle=\int_{0}^{\infty}\\!\\!\\!\\!\\!\\!\mathrm{d}s\int_{\mathbf{r}(0)=\mathbf{x}^{\prime}}^{\mathbf{r}(s)=\mathbf{x}}\\!\\!\\!\\!\\!\\!\mathcal{D}\left[\mathbf{r}(s)\right]e^{-\int_{0}^{s}\mathrm{d}s^{\prime}\frac{1}{4}\left(\frac{d\mathbf{r}(s^{\prime})}{ds^{\prime}}\right)^{2}+\frac{1}{4}\beta^{2}\nabla
H(\mathbf{r}(s^{\prime})).\nabla
H(\mathbf{r}(s^{\prime}))-\frac{1}{2}\beta\nabla^{2}H(\mathbf{r}(s^{\prime}))}$
(7)

Upon the susbtitution $V(\mathbf{r})=\frac{1}{4}\beta^{2}\nabla
H(\mathbf{r}(s^{\prime})).\nabla
H(\mathbf{r}(s^{\prime}))-\frac{1}{2}\beta\nabla^{2}H(\mathbf{r}(s^{\prime}))$
it can be observed that $K(\mathbf{x},\mathbf{x}^{\prime})$ is the path
integral solution to the wick-rotated Schrödinger equation, integrated over
all “times”, but where time in this system is equal to spatial path length
$s$. We have introduced the variable $s$ here, which has units of $[L]^{2}$.
This variable characterizes distance along a path. We make the analogy more
complete by characterizing the term in the exponential as a “Lagrangian” (as
defined in equation (5)). We note that the origin of expression (4) needs to
take care of Weyl ordering of operators, for a full discussion, see [33].

This result isn’t novel in the sense that it can be compared to the standard
results obtained in linear response theory, however in this case we choose to
represent the Green’s function in path integral form, which is not standard in
previous literature, for reasons of notation and interpretation. As shall be
described immediately below, each of the terms that appear in the expression 4
can be associated with a physical meaning that can make it more apparent how
the presence of dissipative forces in a complicated system can affect the
statistical weight of a particular state through understanding of the
equilibrium landscape of $H$. This also allows us to bring certain physically
justified approximations to bear in the evaluation of microstate probability.
Moreover, and more practically, the usage of path integrals in quantum field
theory has established methods for evaluation in different contexts, including
anharmonic corrections, which we shall occupy ourselves with later in the
manuscript, and poses quite general challenges for out of equilibrium systems.

The derivation of this expression (4) results from mathematical
considerations, interestingly however, we observe in the final expression (4)
that the terms in the exponential have physical interpretation. The effective
“Lagrangian” (weight of the exponential) has three terms, which can be
summarized as:

$L=\text{total length of path}+\text{total power dissipated along
path}+\text{total work done along path}$ (8)

This is using the fact that the total power dissipated along the path has been
associated with the expression $\nabla H.\nabla H-\nabla^{2}H$ in prior works
[34, 35]. The total path length is given by
$\big{|}\frac{\mathrm{d}\mathbf{r}(s)}{\mathrm{d}s}\big{|}^{2}$ and the work
along the path is given by $\nabla
H.\frac{\mathrm{d}\mathbf{r}(s)}{\mathrm{d}s}$. As the equilibrium system is
conservative, the final term can be integrated out explicitly. This suggests
that even in equilibrium, every two points in the state space are connected by
the integral over all paths weighted by this quantity. Therefore, overly long
paths have little effect, but so do paths where a large amount of work is done
by the system. Interestingly however, paths which also dissipate a lot of
power also have minimal effects. This is a purely equilibrium function, and it
could therefore be a source of confusion to speak of “dissipation” at
equilibrium. While the average of $\nabla H.\nabla H-\nabla^{2}H$ can be shown
to be exactly zero at equilibrium (meaning on average the system doesn’t
dissipate power), the same is not true for the evaluation of this quantity
along a path connecting any two points.

The presence of this path dependent, non-local term makes this a difficult
path integral to evaluate in general. However, a simplification arises in the
cases that the Hamiltonian is slowly varying between points $\mathbf{x}$ and
$\mathbf{x}^{\prime}$, in this case $\frac{1}{4}\beta^{2}\nabla
H(\mathbf{r}(s^{\prime})).\nabla
H(\mathbf{r}(s^{\prime}))-\frac{1}{2}\beta\nabla^{2}H(\mathbf{r}(s^{\prime}))\approx
0$ we can ignore the “difficult” parts of the path integral, leading to a
form:

$G(\mathbf{x},\mathbf{x}^{\prime})\approx\frac{e^{\frac{\beta}{2}\left(H(\mathbf{x})-H(\mathbf{x}^{\prime})\right)}}{|\mathbf{x}-\mathbf{x}^{\prime}|^{3N-2}}$
(9)

as was seen previously[25]. In most situations, as the effect of a
perturbation at a point $\mathbf{x}$ declines very quickly in the total
distance between two microstates $|\mathbf{x}-\mathbf{x}^{\prime}|$, the
effect of the dissipated power term will be rather small, as the variation in
the Hamiltonian doesn’t have an opportunity to be important before the other
terms come to dominate the Green’s function and therefore set $G$ to a
infinitesimally small value.

Therefore, one way of understanding the microscopic probability in a non-
equilibrium situation involves summing all the paths (weighted by the
Lagrangian given above) between those two points in equilibrium. This tells us
the effect of a impulse of “dissipated power” at a point $\mathbf{x}^{\prime}$
on the microscopic probability at any point $\mathbf{x}$.

For completeness, we note that the same function $G$ given above appears in
the second order perturbation. The second order perturbation, which goes
beyond linear changes in the nonconservative force $\mathbf{F}_{nc}$ is given
by:

$d\log
P_{2}(x)=\int\mathrm{d}\mathbf{x}^{\prime}\int\mathrm{d}\mathbf{x}^{\prime\prime}G(\mathbf{x},\mathbf{x}^{\prime})\left(\beta\mathbf{F}_{nc}(\mathbf{x}^{\prime}).\nabla^{\prime}-\rho(\mathbf{x}^{\prime})\right)G(\mathbf{x}^{\prime},\mathbf{x}^{\prime\prime})\rho(\mathbf{x}^{\prime\prime})-\frac{1}{2}\left(\int\mathrm{d}\mathbf{x}^{\prime}G(\mathbf{x},\mathbf{x}^{\prime})\rho(\mathbf{F}(\mathbf{x}^{\prime}))\right)^{2}$
(10)

In fact, using the operator given above the full microscopic probability for a
system subject to a non-conservative force $\mathbf{F}_{nc}$ can be given in
complete form as:

$P(\mathbf{x})=\frac{1}{Z_{nc}}\exp\left(-\beta
H(\mathbf{x})+\log\left(\left[1-\int\mathrm{d}\mathbf{x}^{\prime}G(\mathbf{x},\mathbf{x}^{\prime})\left(\beta\mathbf{F}_{nc}(\mathbf{x}^{\prime}).\nabla^{\prime}-\beta^{2}\nabla
H(\mathbf{x}^{\prime}).\mathbf{F}_{\text{nc}}(\mathbf{x}^{\prime})\right)\right]^{-1}\mu_{0}\right)\right)$
(11)

where $\mu_{0}$ is any constant and $Z_{nc}$ is a normalization constant. The
term in the square brackets should be understood as an operator acting on
$\mu_{0}$. This is assuming that the Neumann series given here converges.
Therefore, given a closed form of $G$, we should be able to calculate
arbitrary microscopic probability distributions under non-conservative forces.
The next sections of the paper will be occupied with calculating integrals
over $G$, and in which limits the simpler approximate forms presented here can
be utilized.

##  III Evaluation of Green’s function integrals

In the previous section, we stated that the probability of a state of any
system perturbed with a non-conservative force is mediated through its Green
function, which is an equilibrium property of the system. A useful example to
explore our expressions occurs when the Hamiltonian is harmonic:

$H_{h}(\mathbf{x})=\frac{1}{2}\mathbf{x}.\underline{\kappa}.\mathbf{x}$ (12)

Now we imagine a system with such a Hamiltonian being perturbed with a source
term $\rho(\mathbf{x})=\mathbf{x}.\underline{B}.\mathbf{x}$

A long calculation can follow where we can calculate what the first order
effect of such a non-equilibrium perturbation would be. The details of this
calculation are rather cumbersome, but are reproduced in the appendix B. The
final result for the appearance of the first order log microscopic probability
is given by:

$\log(P(\mathbf{x}))\approx-\frac{1}{2}\beta\mathbf{x}.\underline{\kappa}.\mathbf{x}-\frac{1}{2}\beta^{2}\mathbf{x}.\underline{C}.\mathbf{x}$
(13)

We have also introduced a new matrix $\underline{C}$, which results from
integration over the Green’s function, it is given by:

$\displaystyle\int\mathrm{d}\mathbf{x}^{\prime}G(\mathbf{x},\mathbf{x}^{\prime})\left(\mathbf{x}^{\prime}.\underline{B}.\mathbf{x}^{\prime}\right)=\mathbf{x}.\underline{C}.\mathbf{x}+\text{Tr}(\underline{\kappa}^{-1}\underline{B})$
(14)
$\displaystyle\underline{C}=\int_{0}^{\infty}\mathrm{d}s\exp(-\beta\underline{\kappa}s).\underline{B}.\exp(-\beta\underline{\kappa}s)$
(15)

We can see that the Green’s function integral over the term
$\mathbf{x}^{\prime}.\underline{B}.\mathbf{x}^{\prime}$ reproduces a harmonic
term, with an additional constant given by a matrix trace. In other words, a
harmonic perturbation has a harmonic form after integration, plus some
constant. This constant is irrelevant to the relative probability of any two
states, however, and for divergence free non-conservative force the trace of
these matrices will be equal to zero anyway.

Of interest is how the integral in (14) applies to source functions that go
as, for example $\rho(\mathbf{x})=\Theta_{ijkl}x_{i}x_{j}x_{k}x_{l}$ where
$\Theta$ is an object which gives coupling constants for entries $i,j,k,l$ of
$\mathbf{x}$ (where repeated indices are summed, Einstein convention is
employed, each index corresponds to one entry of the vector $\mathbf{x}$).
This would correspond to, for example a system with an anharmonic Hamiltonian
being driven by a non-conservative force. A similar calculation can be
performed in this case (see appendix B), where the following result is
obtained:

$\displaystyle\int\mathrm{d}\mathbf{x}^{\prime}G(\mathbf{x},\mathbf{x}^{\prime})\left(\Theta_{ijkl}(x^{\prime})_{i}(x^{\prime})_{j}(x^{\prime})_{k}(x^{\prime})_{l}\right)=$
$\displaystyle\int_{0}^{\infty}\mathrm{d}s\Theta_{ijkl}(\chi(s))^{i}(\chi(s))^{j}(\chi(s))^{k}(\chi(s))^{l}$
(16)
$\displaystyle+\sum_{\Psi}\int_{0}^{\infty}\mathrm{d}s\mu(s)_{ij}\Theta_{ijkl}(\chi(s))_{k}(\chi(s))_{l}$
(17) $\displaystyle+\sum_{\Psi}\mu(s)_{ij}\mu(s)_{kl}\Theta_{ijkl}$ (18)

where the sum over $\Psi$ corresponds to sums over relevant permutations in
the indices $i,j,k,l$ (the method of generation is given in the appendix) and
we have introduced the matrix $\underline{\mu}(s)$ and vector $\bm{\chi}(s)$
that are given by:

$\displaystyle\underline{\mu}(s)$
$\displaystyle=\left(\beta\underline{\kappa})^{-1}.(\underline{I}-\exp(-2\beta
s\underline{\kappa})\right)$ (19) $\displaystyle\bm{\chi}(s)$
$\displaystyle=\exp\left(-\beta\underline{\kappa}s\right).\mathbf{x}$ (20)

$\exp$ here refers to the matrix exponential, and $I$ is the identity matrix.
The derivations of these are also given in the appendix.

From the above, it can be seen that given a base harmonic Hamiltonian, the
Green’s function can be applied mechanically onto any kind of source function
$\rho(\mathbf{x})$, reminiscint of Feynman rules. For example, to source
function going as $\sim x^{6}$, one could apply the following rule: contract
the six index tensor with the matrix $\underline{\mu}$, zero, one, two and
three times, and for any remaining indices contract them over the vector
$\bm{\chi}$. The only remaining facet is to be careful over the indices being
contracted.

We can combine the above discussion with the expression for the Green’s
function in a generic system, not just for a harmonic Hamiltonian (see
appendix C). We imagine the equilibrium Hamiltonian is given by:
$H(\mathbf{x})=H_{h}(\mathbf{x})+V(\mathbf{x})$ where $H_{h}$ are the parts of
the Hamiltonian that go at most quadratically in $\mathbf{x}$.

In such a scenario, the Green’s function is given by[36]:

$G(\mathbf{x}^{\prime\prime},\mathbf{x}^{\prime})=e^{\frac{\beta}{2}\left(V(\mathbf{x}^{\prime\prime})-V(\mathbf{x}^{\prime})\right)}\int_{0}^{\infty}\mathrm{d}sG_{h}(\mathbf{x}^{\prime\prime},\mathbf{x}^{\prime},s)\exp\left[\int_{0}^{s}\mathrm{d}s^{\prime}V_{a}\left(\frac{\delta}{\delta\mathbf{J}(s^{\prime})}\right)\right]\exp(-W[\mathbf{J}])\big{|}_{\mathbf{J}=0}$
(21)

Where the function $V_{a}$ is given by:

$V_{a}(\mathbf{x})=\nabla V(\mathbf{x}).\nabla V(\mathbf{x})+2\nabla
V(\mathbf{x}).\nabla H_{h}(\mathbf{x})-\nabla^{2}V(\mathbf{x})$ (22)

, and $W[\mathbf{J}]$ is a functional of an arbitrary function $\mathbf{J}$
and is given by:

$W[J]=\int\mathrm{d}s\mathbf{J}(s).\mathbf{x}_{cl}(s)+\frac{1}{2}\int\mathrm{d}s\int\mathrm{d}s^{\prime}\mathbf{J}(s).\underline{G}_{D}(s,s^{\prime}).\mathbf{\underline{missing}}\mathbf{J}(s^{\prime})$
(23)

where $\underline{G}_{D}(s,s^{\prime})$ is an analytic functions given in the
appendix and $G_{h}(\mathbf{x},\mathbf{x}^{\prime},s)$ is given by:

$G_{h}(\mathbf{x},\mathbf{x}^{\prime},s)=e^{\frac{\beta}{2}\left(H_{h}(\mathbf{x}^{\prime\prime})-H_{h}(\mathbf{x}^{\prime})\right)}\left(\frac{1}{2\pi}\right)^{D/2}\sqrt{\det\left(\frac{1}{2}\underline{\kappa}\text{csch}\left(\underline{\kappa}s\right)\right)}\exp\left(-S_{cl}\right)$
(24)

with $S_{cl}$ as the classical action integrated over the classical path,
which can be found in the appendix.

The above introduced machinery appears formidable, but can be utilized in
practice to derive corrections for anharmonic systems subject to non-
conservative forces. As the applications of the Green’s functions integrals
are mechanical, all the above taken together can be used to derive leading
order terms for anharmonic systems subjected to non-conservative forces. We
shall see that interesting effects don’t require the employment of the
entirety of this machinery. Equation (21) is exact, nevertheless it is
difficult to use in practice, as is the usual story with anharmonic systems.
To that end we investigate progressively more accurate representations of
equation (21), firstly by ignoring the terms going as $V(\mathbf{x})$ in the
exponential, the harmonic approximation:

$G(\mathbf{x}^{\prime\prime},\mathbf{x}^{\prime})\approx\int_{0}^{\infty}\mathrm{d}sG_{h}(\mathbf{x}^{\prime\prime},\mathbf{x}^{\prime},s)$
(25)

then doing integrations over this kernel with the endpoint values (anharmonic
integral approximation):

$G(\mathbf{x}^{\prime\prime},\mathbf{x}^{\prime})\approx
e^{\frac{\beta}{2}\left(V(\mathbf{x}^{\prime\prime})-V(\mathbf{x}^{\prime})\right)}\int_{0}^{\infty}\mathrm{d}sG_{h}(\mathbf{x}^{\prime\prime},\mathbf{x}^{\prime},s)$
(26)

and finally, by including the first order effect of the dissipation along the
path at equilibrium (anharmonic integral with dissipative correction
approximation)::

$G(\mathbf{x}^{\prime\prime},\mathbf{x}^{\prime})\approx
e^{\frac{\beta}{2}\left(V(\mathbf{x}^{\prime\prime})-V(\mathbf{x}^{\prime})\right)}\int_{0}^{\infty}\mathrm{d}sG_{h}(\mathbf{x}^{\prime\prime},\mathbf{x}^{\prime},s)\left(1+\int_{0}^{s}\mathrm{d}s^{\prime}V_{a}\left(\frac{\delta}{\delta\mathbf{J}(s^{\prime})}\right)\right)\exp(-W[\mathbf{J}])\bigg{|}_{\mathbf{J}=0}$
(27)

###  III.1 Single Brownian particle in harmonic or anharmonic trap subject to
active noise

A useful application of the previous equations is in the calculation of the
probability distribution of a particle bound in some anharmonic or harmonic
potential subject to some active noise $A$. The harmonic form of this problem
is a canonical model in prior literature[37]. Such a system can be described
by the following equations:

$\displaystyle\gamma\frac{\mathrm{d}x}{\mathrm{d}t}$
$\displaystyle=-q^{(\alpha)}x^{\alpha}+aA(t)+\xi(t)$ (28)
$\displaystyle\gamma\frac{\mathrm{d}A(t)}{\mathrm{d}t}$
$\displaystyle=-\frac{\gamma}{\tau_{a}}A(t)+\xi(t)$ (29)

where $\gamma$ is the viscous damping, the white noise $\xi(t)$ obeys the
following relationship
$\langle\xi(t)\xi(t^{\prime})\rangle=2k_{b}T\gamma\delta(t-t^{\prime})$ with a
temperature $k_{b}T$, $\tau_{a}$ is some timescale of relaxation of the active
noise $A$, and $a$ characterizes the magnitude of the impact of that noise on
the spatial variable $x$. The integer $\alpha$ is $\alpha=1$ for a harmonic
system and $\alpha=3$ for an anharmonic one. The parameter $q^{(\alpha)}$
gives the magnitude of the restoring force. We include the superscript
$(\alpha)$ to make clear that the units of $q^{(\alpha)}$ will vary depending
on $\alpha$, $[q^{(\alpha)}]=[E][L]^{-\alpha+1}$. The conservative and non-
conservative parts of the system 28 can be found via symmetrizing the forces,
yielding
$\mathbf{F}_{c}=\left(-q^{(\alpha)}x^{\alpha},-\frac{\gamma}{\tau_{a}}A\right)$
and $\mathbf{F}_{nc}=\left(aA,0\right)$. Using the source function defined in
equation (2), we have for system (28)

$\rho(\mathbf{x})=q^{(\alpha)}ax^{\alpha}A$ (30)

We can go through the process of applying the machinery of the previous
section to such a source term. We can also ask what level of approximations of
the Green’s function is necessary for the calculation of accurate results. For
a harmonic system ($\alpha=1$), the matrix $\kappa$ quoted in the previous
section is given by $\kappa=\begin{bmatrix}&q^{(1)},&0\\\
&0,&\gamma/\tau_{a}\end{bmatrix}$, after performing the integral given in (14)
we obtain for the first order perturbation:

$P_{\text{harmonic}}(x,A)\sim\exp\left(-\frac{\beta}{2}q^{(1)}x^{2}-\frac{\beta}{2}\frac{\gamma}{\tau_{a}}A^{2}+\beta\frac{aq^{(1)}\tau_{a}}{q^{(1)}\tau_{a}+\gamma}xA\right)$
(31)

where $\beta=1/k_{b}T$. As expected, the introduction of a driving noise $A$
onto the spatial variable $x$ leads to a term in the probability that
correlates $x$ and $A$. When viewed in the space of only $x$ (integrating out
over $A$), this would correspond to increased variance of $x$, as would be
expected and is known from an active noise acting on a particle.

A similar calculation (see appendix) can be performed for an active particle
in an anharmonic trap ($\alpha=3$), yielding:

$P_{\text{anharmonic}}(x,A)\sim\exp\left(-\frac{\beta}{4}q^{(3)}x^{4}-\frac{\beta}{2}\frac{\gamma}{\tau_{a}}A^{2}+cxA+dx^{3}A\right)$
(32)

which is the first term for $q^{(3)}$ is small. Full integration over the
anharmonic terms can be performed, however, a simple approximation results if
we take only the leading order correction over the harmonic form of the
Green’s function (which we call the harmonic approximation introduced in the
previous section), which leads to the form of the variables $c$ and $d$, which
goes as:

$\displaystyle c$
$\displaystyle=\frac{3aq^{(3)}\tau_{a}^{2}}{\gamma^{2}}+\mathcal{O}(q^{2}\tau_{a}^{4}/\gamma^{4})$
(33) $\displaystyle d$ $\displaystyle=\frac{\beta
aq^{(3)}\tau_{a}}{2\gamma}+\mathcal{O}(q^{2}\tau_{a}^{3}/\gamma^{3})$ (34)

Whereupon the rest of the series can be seen as including terms of higher
powers of $\tau_{a}/\gamma$

Comparison of analytic results with simulations of a Brownian particle under
an active force are presented in figure 1. For ease of comparison, we compare
both harmonic and anharmonic trapping. It is plainly observable that a term of
the form $xA$ is present in the probability density. Presented are also the
probability densities integrated out over the active noise $A$ in figures 1B
and 1C. The full microscopic probability densities are shown for both theory
and simulation in D-G, and display good agreement.

Figure 1:  The impact of an active noise $A$ on a Brownian particle in a
harmonic or anharmonic well. A) schematic showing the particle in two
different kinds of wells. B) Simulations of a Brownian particle in a harmonic
well compared to theoretical predictions of the same arising from the first
order perturbation for $q^{(1)}=1,\gamma=10,\tau_{a}=10,a=0.5$ C) The same as
in B but for an anharmonic well with $q^{(3)}=0.1,\gamma=10,\tau_{a}=10,a=0.5$
(taking the anharmonic integral approximation). D) and E) comparison of the
full stationary probability in the space of $x$ and $A$ for the harmonic
system with the parameters just mentioned, measured in simulation (E) and in
theory (F). A coupling going as $xA$ is introduced by the non-conservative
force. F) and G), the same but for the anharmonc system. As in the harmonic
system, a coupling going as $xA$ is visible.

At first glance, it is somewhat surprising that we can obtain reasonable
results in this regime, as we are perturbing around the harmonic system where
there is no harmonic contribution to the variable $x$ at equilibrium, and
$\gamma/\tau_{a}$ is not very large in comparison to $a$. This points to some
robustness in the perturbative scheme developed thus far.

One important qualitative result is the term going as $xA$, which appears in
both harmonic and anharmonic traps, but we have yet to discuss its impact.
Generally, the Green’s function integrated out over a quartic source function
will bring out a harmonic contribution (i.e. $xA$). If we use the following
integral identity:
$\int\mathrm{d}A\exp(-c_{1}x^{4}-c_{2}A^{2}+c_{3}xA)=\text{constant}\times\exp\left(\frac{c_{3}^{2}x^{2}}{4c_{2}}-c_{1}x^{4}\right)$
it can be observed that an anharmonic system subject to non-conservative
forces will begin to accumulate probability away from $x=0$ (which is its
equilibrium minima) as the maxima of the probability will now occur at
$x=c_{3}/\sqrt{8c_{1}c_{2}}$. Such an effect isn’t observable for the
magnitudes of non-conservative forces in figure 1, but if the non-conservative
force is made stronger, it can be seen; as in figure 2B where probability
begins to accumulate away from $x=0$, in contrast to the equilibrium system.
Stronger drives also increase the importance of the term going as $x^{3}A$, as
can be observed in figure 2A.

What matters here is not so much the magnitude of the factor $c$, but the form
of the coupling $xA$. This suggests less exact methods and/or approximate
Green’s funcions may get the precise magnitude of the coupling terms
incorrect, but still obtain qualitatively correct couplings. More generally,
Green’s function integrals over any anharmonic system should bring out terms
of lower order whose magnitude can be important for larger amounts of non-
conservative force. Therefore, accumulation away from the equilibrium mean
should be expected as a generic result of anharmonic systems under non-
conservative forces[30], as contrasted to harmonic systems.

The approximate expressions given in equation (33) can be compared against
simulations with strong drives, by increasing $\tau_{a}$. The resulting
measured distributions can then be fit to the form
$P(x)=\exp(-f_{1}(\beta/4)q^{(\alpha)}x^{4}-f_{2}\beta(\gamma/2\tau_{a})A^{2}+cxA+dx^{3}A)$
and the values of $c$ and $d$ can be measured and compared to against
expectations from the theory, which is seen in figures 2C (we also include
additional fitting parameters $f_{1}$ and $f_{2}$ to account for higher than
first order effects in the simulations). The value of the parameters can be
computed not only for the simple approximate expressions presented above, but
for more detailed approximations to the Green’s function
$G(\mathbf{x},\mathbf{x}^{\prime})$. We focus on the two further
approximations given previously, the one given in equation (26) and the result
for the first order perturbation of the dissipation-along-the-path
contribution to the path integral arising in equation (27). This allows us to
probe the importance of different contributions to the Green’s function in the
magnitude of the coupling constants $c,d$. Figure 2C shows that the
approximation (9) is excellent even up to quite large values of $\tau_{a}$.
This is notable as at large $\tau_{a}$ $(>10)$, the non-conservative force
exceeds the conservative one. When the additional dissipation-along-the-path
contribution is included, the value of $c$ measured in simulation is better
captured, though the prediction for $d$ actually becomes worse. We mentioned
previously that the contribution of the dissipation along the path for the
Green’s function becomes important only when the Hamiltonian varies quickly,
which is why it’s introduction is important for capturing the value of $c$ at
large $\tau_{a}$, for under such noise the system is probing the quartic
region of the Hamiltonian and therefore the terms going as $\nabla.H.\nabla
H-\nabla^{2}H$ become non-negligible.

As a final point, the above illustrates another far simpler method towards
discovery of modifications of equilibrium distributions in the presence of
non-conservative forces. By symmetry, the contributions to the probability
distribution of the non-conservative force have odd parity $f(-x)=-f(x)$.
Therefore, we can reason out the additional terms in a quartic system as by
necessity being either $xA,x^{3}A$ or $A^{3}x$. Dimensional analysis can yield
the approximate scalings of the prefactors. Such a method would not be able to
derive numerical prefactors, however when attempting to construct effective
theories these could simply be parameterized.

Figure 2:  Probability distributions under larger magnitudes of the active
noise. A) the probability distribution in both $x$, $A$, from which a coupling
going as $x^{3}A$ is visible. B) Integration over $A$ yields the probability
distribution of the particle in the trap. The maxima of this distribution is
now away from $x=0$. C) Simulations of an active Brownian particle in an
anharmonic trap where the measured distributions are fit to equation 32 for
the parameters $c$ and $d$. The values of these parameters can also be
calculated in theory under different approximations (dashed lines) where it
can be seen that approximate forms of the Green’s function nevertheless
capture the correct scaling up to quite large $\tau_{a}$.

##  IV Larger systems and mode space analysis: Application to Polymers

Using some of the results from the previous section, we can briefly try to
learn about the properties of larger systems based on the modifications of the
microscopic probability that occurs due to non-conservative forces. Larger
systems pose quite general challenges for microscopic descriptions as when the
number of degrees of freedom is large, it is difficult to even represent
microscopic probabilities. Thus, we want to describe some more macroscopic
observables through consideration of the microscopic probability distributions
derived in the previous section. A natural candidate system for larger systems
are active polymers [38, 39, 40, 41], given that these have both non-
conservative forces and harmonic (or anharmonic) backbone potentials. Such a
system can be thought of as a collection of coupled particles in wells of the
previous section. We imagine a polymer of beads each of whose monomers $i$ and
$i+1$ with position vectors $\mathbf{x}_{i}$ and $\mathbf{x}_{i+1}$ are
connected by a spring potential of the form:

$\phi(\mathbf{x}_{i+1},\mathbf{x}_{i})=\frac{1}{\alpha+1}q^{(\alpha)}(\mathbf{x}_{i+1}-\mathbf{x}_{i})^{\alpha+1}$
(35)

where once more $\alpha=1$ is a harmonic system and $\alpha=3$ is an
anharmonic one. The full Hamiltonian of this polymer can be written as
$H=\sum_{i}^{N}\phi(\mathbf{x}_{i+1},\mathbf{x}_{i})$. We subject each monomer
in this system to an active noise of the form given in the previous section,
such that the time evolution equation of each individual monomer will be given
by:

$\gamma\frac{\mathrm{d}\mathbf{x}_{i}}{\mathrm{d}t}=-\nabla_{\mathbf{x}_{i}}\left(\phi(\mathbf{x}_{i+1},\mathbf{x}_{i})+)\phi(\mathbf{x}_{i},\mathbf{x}_{i-1})\right)+a\mathbf{A}_{i}(t)+\bm{\xi}(t)$
(36)

where each entry of vector $\mathbf{A}_{i}$ obeys equation (29). This is an
extended version of the system in the previous example. We may now ask about
the stationary properties of such a system, but it would be useful to compare
the results to a similar equilibrium system. We already saw that the tendency
of the addition of active noise is to bring forth a harmonic term for an
anharmonic Hamiltonian that leads to accumulation away from the equilibrium
minima. We therefore compare these driven harmonic and anharmonic polymers
with an equilibrium polymer with the following potential:

$\Phi(\mathbf{x}_{i+1},\mathbf{x}i)=-\frac{1}{2}k(\mathbf{x}_{i+1}-\mathbf{x}_{i})^{2}+\frac{1}{4}q(\mathbf{x}_{i+1}-\mathbf{x}_{i})^{4}$
(37)

where $k,q$ are positive constants. We concentrate on two different aspects of
the polymer. One particular macroscopic quantity which is natural to study in
polymers is the average end to end distance $\langle R^{2}\rangle$ which is
the average distance between the two endpoints of the polymer. Another
quantity we wish to analyze is how the modes of the polymer are modified in
the presence of drive. Mode space analysis is standard in polymer systems,
defining a set of collective coordinates of the entire polymer, i.e.
$\tilde{\mathbf{x}}=\underline{Q}.\mathbf{x}$ where $\underline{Q}$ is a
matrix of eigenvectors. Each mode is an independent degree of freedom with
some associated eigenvalue $\lambda$. We know that in equilibrium that each of
these modes should follow the distribution
$\sim\exp(-\beta\lambda_{i}\tilde{x}^{2})$. In other words, in a harmonic
system under equilibrium, each mode will have a variance that goes as the
reciprocal of its eigenvalue. By analyzing the fluctuations of modes in an
active system we can see roughly how “unusual” such a system is compared to an
equilibrium one if we divide through by
$\exp(-\beta\lambda_{i}\tilde{x}^{2})$. By equipartition, for an equilibrium
system, the variance of each mode divided by $1/\lambda_{i}$ should be a
constant, however this no longer holds for non-equilibrium systems. Moreover,
anharmonic polymers can be linearized around their minima, yielding similar
mode space distributions as harmonic polymers. The analysis of all these
quantities is presented in figure LABEL:fig:fig3.

Figure 3:  Simulations of an active polymer. A) The polymer is being driven by
some active noise $A$, with $\tau_{a}=10.$,$\gamma=10$,$k_{b}T=2.$, and we
test our results for both a harmonic and anharmonic backbone potential, as
well as against an equilibrium (no active noise) polymer containing both
harmonic and anharmonic backbone potentials (equation (37)) (see text for
reasoning). B) Measured end to end distance for harmonic ($q^{(1)}=1$) and
anharmonic($q^{(3)}=0.1$) polymers under active noise. Also shown, an
equilibrium polymer with $k=a$ and ($q^{(3)}=0.1$). The size of the active
harmonic polymer scales with noise as $a^{2}$, and the equilibrium polymer
with $k^{5/4}$. The anharmonic active polymer will scale with a crossover
between these two regimes due to the interplay of backbone potential and
active noise. C) Variances of modes of the polymer measured against their
eigenvalues. In the absence of active noise, these lines would be flat,
however, when active noise is introduced the modes with smaller eigenvalues
are excessively affected. D) Same as C) but at higher values of the drive.

In figure 3B we can observe several features that arise directly from the
consideration of our microscopic model in the previous section. Firstly, the
interaction between two neighboring monomers being driven by an active force
$A$ of magnitude $a$ should have a microscopic distribution that has the form
$P(x,A)\sim\exp(-c_{1}x^{4}-c_{2}A^{2}+ac_{3}xA)$ which leads to an effective
microscopic probability arising between the monomers of form
$P(x)\sim\exp({\frac{a^{2}c_{3}^{2}x^{2}}{4c_{2}}-c_{1}x^{4}})$. This means
that for low $a$, the effective interaction between monomers is unmodified,
but as $a$ is increased to $a>1$, the monomers should start to accumulate at a
preferred length given by the interplay of the potential and the active noise.
This can be compared to the equilibrium model defined by equation (37). In
such a model, the total end to end distance scales against the parameter $k$
as $k^{5/4}$, as can be seen in figure 3B. This is the same scaling observed
in the active model against $a$ when $a>1$, indicating that the accumulative
effect of the active noise has a similar effect as the interaction in the
equilibrium model. We may contrast this with the active harmonic model, where
the end to end distance of the polymer chain always swells as $a^{2}$, as
there is no accumulative effect away from the equilibrium minimum in harmonic
models.

Differences to the equilibrium model can be observed when looking at the
variances of the modes of the polymer. It can be shown (see appendix D) how
each mode is affected by the presence of the non-conservative force. As the
modes are orthogonal degrees of freedom, in an equilibrium system we should
expect by equipartition that the total energy per mode is the same, as is
observed in simulations of the equilibrium model (and can be seen in figure
3D). In contrast, for active systems equipartition no longer holds. What can
be derived in the harmonic system and also holds in the anharmonic system is
that modes with smaller eigenvalues (i.e. modes that have greater fluctuations
under equilibrium) are more affected by the drive relative to their
equilibrium values than the modes with large eigenvalues. For an active
polymer, this means that the fluctuations of each mode should follow a form
given by approximately $1/(1+\lambda^{2})$ (see appendix D) where $\lambda$ is
the eigenvalue of the mode. As anharmonic systems also adopt harmonic
distributions around the minima in the free energy, it is observed empirically
(in both figure 3C and 3D). This effect becomes more pronounced the stronger
the active noise $a$ is.

##  V Discussion and Conclusion

In the preceding, we have demonstrated that the microscopic probability for a
system perturbed by a deterministic non-conservative force or due to coupling
with different heat baths can be shown to be equal to the equilibrium system
with an additional term appearing as an integral over a microscopic stationary
response (or Green’s) function. This response function is the same for both of
the above non-equilibrium perturbations. The microscopic response function is
a path integral in the unperturbed system. In equilibrium, any two states
$\mathbf{x}$ and $\mathbf{x}^{\prime}$ are connected by the sum of all paths
between them weighted by the path length, work and dissipated power associated
with each of the paths. The modification of the microscopic probability is
analytically solvable for harmonic systems. This formulation can be applied,
though not without some difficulty, to anharmonic systems as well. It can be
shown that purely anharmonic systems should see harmonic terms arising in
their microscopic probability, from which it is apparent that accumulation of
probability away from their equilibrium minima is a generic feature of
anharmonic but not harmonic systems. Finally, we used the results we obtained
for a Brownian particle in an anharmonic trap to understand the properties of
anharmonic polymers under active noise.

The microscopic stationary response function was constructed through
consideration of the dynamical processes involved, and yet its final form has
a more physical interpretation. It would be interesting to start from the
effective Lagrangian formulation as it appears in (8) and apply it to a
completely different system without considering any dynamics and observe
whether it can recapitulate the stationary distribution successfully. This
would be suggestive as to the importance of the response function we have
shown in this paper.

The path integral that represents the response function is summed over all
path lengths, where the path length takes a form that looks suspiciously like
time. In this paper, we have ignored dynamic effects, this can prompt us to
ask the extent to which the effect of time appears as the path length
representation that is already present in the over-damped system. Following on
from this, we have ignored inertial effects here, but the theory could be
modified to account for them. One could then construct dynamical response
functions with explicit time, and test for modifications to the fluctuation-
dissipation relation.[18, 42]

The fact that the soft modes (modes with the smallest eigenvalues) of the
equilibrium system are most affected by the presence of non-equilibrium
driving is intriguing. Interesting phenomena due to dissipative forces would
arise in the case where the eigenspectrum of the original system was quite
broad. In this case, some modes are basically unaffected by the drive, however
the modes which were most strongly fluctuating in the equilibrium system are
those which will be most strongly affected by the dissipative forces. This
could lead to potentially quite exotic states due to the modification of the
microscopic probabilities.

In terms of connection with experimental systems, active matter systems can be
mapped onto this framework when the proper degrees of freedom the system are
considered[25], to this we would add experiments that involve heating
different parts of a system to different temperatures [43]. There has been
recent interest in active matter with what has been deemed odd viscosity[44,
45], as these systems also have non-conservative forces they would be a
natural candidate for investigating how the microscopic phenomena we have
discussed here would bridge length scales to more macroscopic quantities.

## References

  * Ramaswamy [2010] S. Ramaswamy, The mechanics and statistics of active matter, Annual Review of Condensed Matter Physics 1, 323 (2010).
  * Needleman and Dogic [2017] D. Needleman and Z. Dogic, Active matter at the interface between materials science and cell biology, Nature Reviews Materials 2, 10.1038/NATREVMATS.2017.48 (2017).
  * Stillinger and Debenedetti [2013] F. H. Stillinger and P. G. Debenedetti, Glass Transition Thermodynamics and Kinetics, http://dx.doi.org/10.1146/annurev-conmatphys-030212-184329 4, 263 (2013).
  * Netz [2020] R. R. Netz, Approach to equilibrium and nonequilibrium stationary distributions of interacting many-particle systems that are coupled to different heat baths, Physical Review E 101, 22120 (2020).
  * Cates and Tailleur [2015] M. E. Cates and J. Tailleur, Motility-induced phase separation, Annual Review of Condensed Matter Physics 6, 219 (2015), arXiv:1406.3533 .
  * Weber  _et al._ [2019] C. A. Weber, D. Zwicker, F. Jülicher, and C. F. Lee, Physics of active emulsions, Reports on Progress in Physics 82, 10.1088/1361-6633/ab052b (2019), arXiv:1806.09552 .
  * Cates [2019] M. E. Cates, Active Field Theories, Particles 3, 245 (2019), arXiv:1904.01330 .
  * Toner  _et al._ [2005] J. Toner, Y. Tu, and S. Ramaswamy, Hydrodynamics and phases of flocks, Annals of Physics 318, 170 (2005).
  * Callegari and Volpe [2019] A. Callegari and G. Volpe, Numerical Simulations of Active Brownian Particles,  , 211 (2019).
  * Osmanović and Rabin [2019] D. Osmanović and Y. Rabin, Chemically active nanodroplets in a multi-component fluid, Soft Matter 15, 5965 (2019).
  * Seifert [2012] U. Seifert, Stochastic thermodynamics, fluctuation theorems and molecular machines, Reports on Progress in Physics 75, 126001 (2012), arXiv:1205.4176 .
  * Speck [2016] T. Speck, Stochastic thermodynamics for active matter, EPL 114, 10.1209/0295-5075/114/30006 (2016), arXiv:1603.03195 .
  * Netz [2018] R. R. Netz, Fluctuation-dissipation relation and stationary distribution of an exactly solvable many-particle model for active biomatter far from equilibrium, The Journal of Chemical Physics 148, 185101 (2018).
  * Loi  _et al._ [2008] D. Loi, S. Mossa, and L. F. Cugliandolo, Effective temperature of active matter, Physical Review E - Statistical, Nonlinear, and Soft Matter Physics 77, 051111 (2008), arXiv:0711.4946 .
  * Loi  _et al._ [2011] D. Loi, S. Mossa, and L. F. Cugliandolo, Non-conservative forces and effective temperatures in active polymers, Soft Matter 7, 10193 (2011), arXiv:1105.0806 .
  * Bechinger  _et al._ [2016] C. Bechinger, R. D. Leonardo, C. Reichhardt, G. G. Volpe, G. G. Volpe, R. Di Leonardo, H. Löwen, C. Reichhardt, G. G. Volpe, and G. G. Volpe, Active particles in complex and crowded environments, Reviews of Modern Physics 88, 045006 (2016), arXiv:1602.00081 .
  * De Karmakar and Ganesh [2020] S. De Karmakar and R. Ganesh, Phase transition and emergence of active temperature in an active Brownian system in underdamped background, PHYSICAL REVIEW E 101, 32121 (2020).
  * Caprini and Marini Bettolo Marconi [2021] L. Caprini and U. Marini Bettolo Marconi, Inertial self-propelled particles, The Journal of Chemical Physics 154, 024902 (2021), arXiv:2009.14032 .
  * Moyses  _et al._ [2015] H. W. Moyses, R. O. Bauer, A. Y. Grosberg, and D. G. Grier, Perturbative theory for Brownian vortexes, Physical Review E - Statistical, Nonlinear, and Soft Matter Physics 91, 10.1103/PHYSREVE.91.062144 (2015), arXiv:1504.04607 .
  * Amarouchene  _et al._ [2019] Y. Amarouchene, M. Mangeat, B. V. Montes, L. Ondic, T. Guérin, D. S. Dean, and Y. Louyer, Nonequilibrium Dynamics Induced by Scattering Forces for Optically Trapped Nanoparticles in Strongly Inertial Regimes, Physical Review Letters 122, 10.1103/PHYSREVLETT.122.183901 (2019), arXiv:1812.06804 .
  * Mangeat  _et al._ [2019] M. Mangeat, Y. Amarouchene, Y. Louyer, T. Guérin, and D. S. Dean, Role of nonconservative scattering forces and damping on Brownian particles in optical traps, Physical Review E 99, 10.1103/PHYSREVE.99.052107 (2019), arXiv:1812.09188 .
  * Noh and Lee [2015] J. D. Noh and J. Lee, On the steady-state probability distribution of nonequilibrium stochastic systems, Journal of the Korean Physical Society 66, 544 (2015), arXiv:1411.3211 .
  * Wedemann  _et al._ [2016] R. S. Wedemann, A. R. Plastino, and C. Tsallis, Curl forces and the nonlinear Fokker-Planck equation, Physical Review E 94, 10.1103/PHYSREVE.94.062105 (2016).
  * Bertin [2017] E. Bertin, Theoretical approaches to the steady-state statistical physics of interacting dissipative units, Journal of Physics A: Mathematical and Theoretical 50, 10.1088/1751-8121/AA546B (2017), arXiv:1608.08507 .
  * Osmanović [2021] D. Osmanović, Spatial distributions of nonconservatively interacting particles, Physical Review E 103, 022610 (2021).
  * Cowan [2005] B. Cowan,  _Topics in Statistical Mechanics_ (PUBLISHED BY IMPERIAL COLLEGE PRESS AND DISTRIBUTED BY WORLD SCIENTIFIC PUBLISHING CO., 2005) https://www.worldscientific.com/doi/pdf/10.1142/p392 .
  * Cengio  _et al._ [2019] S. D. Cengio, D. Levis, and I. Pagonabarraga, Linear response theory and green-kubo relations for active matter 10.1103/PhysRevLett.123.238003 (2019).
  * [28] (Annals of Mathematics Studies) Robion C.
  * Kwon  _et al._ [2011] C. Kwon, J. D. Noh, and H. Park, Nonequilibrium fluctuations for linear diffusion dynamics, PHYSICAL REVIEW E 83, 61145 (2011).
  * Das  _et al._ [2018] S. Das, G. Gompper, and R. G. Winkler, Confined active brownian particles: theoretical description of propulsion-induced accumulation, New J. Phys 20, 15001 (2018).
  * Kedia  _et al._ [2019] H. Kedia, D. Pan, J.-J. Slotine, and J. L. England, Drive-specific adaptation in disordered mechanical networks of bistable springs, (2019), arXiv:1908.09332 .
  * Su  _et al._ [2022] J. Su, Z. Cao, J. Wang, H. Jiang, and Z. Hou, Dynamic and thermodynamic origins of motility-induced phase separation, (2022), arXiv:2211.16084 .
  * Gersten and Nitzan [1987] J. I. Gersten and A. Nitzan, Path integral approach to electrostatic problems, The Journal of Chemical Physics 86, 3557 (1987).
  * Tomé and de Oliveira [1997] T. Tomé and M. J. de Oliveira, Stochastic mechanics of nonequilibrium systems, Brazilian Journal of Physics 27, 525 (1997).
  * Tome [2006] T. Tome, Entropy production in nonequilibrium systems described by a Fokker-Planck Equation, Brazilian journal of physics 36, 1285 (2006).
  * Wipf [2013] A. Wipf,  _Path Integrals in Quantum and Statistical Mechanics_ (Springer, Berlin, Heidelberg, 2013) pp. 5–23.
  * Dabelow  _et al._ [2019] L. Dabelow, S. Bo, and R. Eichhorn, Irreversibility in Active Matter Systems: Fluctuation Theorem and Mutual Information, Physical Review X 9, 10.1103/PHYSREVX.9.021009 (2019), arXiv:1806.04956 .
  * Eisenstecken  _et al._ [2017] T. Eisenstecken, G. Gompper, and R. G. Winkler, Internal dynamics of semiflexible polymers with active noise., The Journal of Chemical Physics 146, 154903 (2017).
  * Vandebroek and Vanderzande [2015] H. Vandebroek and C. Vanderzande, Dynamics of a polymer in an active and viscoelastic bath., Physical Review. E, Statistical, Nonlinear, and Soft Matter Physics 92, 60601 (2015).
  * Winkler  _et al._ [2017] R. G. Winkler, J. Elgeti, and G. Gompper, Active polymers —emergent conformational and dynamical properties: A brief review, Journal of the Physical Society of Japan 86, 101014 (2017).
  * Osmanović and Rabin [2017] D. Osmanović and Y. Rabin, Dynamics of active rouse chains., Soft matter 13, 963 (2017).
  * Caprini [2021] L. Caprini, Generalized fluctuation–dissipation relations holding in non-equilibrium dynamics, Journal of Statistical Mechanics: Theory and Experiment 2021, 063202 (2021).
  * Grosberg and Joanny [2015] A. Y. Grosberg and J. F. Joanny, Nonequilibrium statistical mechanics of mixtures of particles in contact with different thermostats, Physical Review E - Statistical, Nonlinear, and Soft Matter Physics 92, 032118 (2015).
  * Banerjee  _et al._ [2017] D. Banerjee, A. Souslov, A. G. Abanov, and V. Vitelli, Odd viscosity in chiral active fluids, Nature Communications 2017 8:1 8, 1 (2017), arXiv:1702.02393 .
  * Markovich and Lubensky [2021] T. Markovich and T. C. Lubensky, Odd viscosity in active matter: Microscopic origin and 3d effects, Physical Review Letters 127, 048001 (2021).
  * Risken [1989] H. Risken, Fokker-planck equation for several variables; methods of solution, in  _The Fokker-Planck Equation_, Springer series in synergetics, Vol. 18, edited by H. Haken (Springer Berlin Heidelberg, Berlin, Heidelberg, 1989) pp. 133–162.
  * Grosche and Steiner [1995] C. Grosche and F. Steiner, How to solve path integrals in quantum mechanics, Journal of Mathematical Physics 36, 2354 (1995).
  * Zee [2010] A. Zee,  _Quantum field theory in a nutshell_ , Vol. 7 (Princeton university press, 2010).

##  Appendix A Microscopic probability of overdamped system with multiple
temperature baths

The stationary state of the Smoluchowski equation for the log probability
$\phi$ of a system with temperature $k_{b}T$ is given as the solution to[46]:

$\nabla.\left(\exp(-\phi)\left(-\mathbf{F}-k_{b}T\nabla\phi\right)\right)=0$
(38)

in order to account for different temperatures, we can modify the equation:

$\nabla.\left(\exp(-\phi)\left(-\mathbf{F}-\underline{\beta}^{-1}.\nabla\phi\right)\right)=0$
(39)

where $\underline{\beta}^{-1}$ is a diagonal matrix where each entry on the
diagonal is given by a value $k_{b}T_{i}$ corresponding to the temperature of
that degree of freedom. This equation can be manipulated to give:

$\nabla.\left(\exp(-\phi)\underline{\beta}^{-1}\left(-\underline{\beta}.\mathbf{F}-\nabla\phi\right)\right)=0$
(40)

Splitting the vector field
$\underline{\beta}.\mathbf{F}=\mathbf{f}^{c}+\mathbf{f}^{a}$ where
$\mathbf{f}_{c}=-\nabla\phi$ is a solution if
$\nabla.\left(e^{-\phi}\underline{\beta}\mathbf{f}^{a}\right)=0$. Assuming
that $\mathbf{F}=-\nabla H_{0}$ this leads to the following equation:

$-\nabla
H.\nabla\phi+\nabla\phi.\underline{\beta}^{-1}.\nabla\phi+\nabla^{2}H-\nabla.(\underline{\beta}^{-1}.\nabla\phi)=0$
(41)

where if every entry into the matrix $\beta$ is equivalent, this reproduces
the Boltzmann distribution.

We can write the temperature matrix as a modification around a base
temperature $\beta_{1}^{-1}$:

$\underline{\beta}^{-1}=\beta_{1}^{-1}\underline{I}+\beta_{2}^{-1}\underline{D}$
(42)

we redefine $H$ as including the base temperature $\beta_{1}^{-1}$

$H\coloneqq\beta_{1}H_{0}$ (43)

We therefore have the following equation in $\phi$:

$-\nabla
H.\nabla\phi+|\nabla\phi|^{2}+\alpha\nabla\phi.\underline{D}.\nabla\phi+\nabla^{2}H-\nabla^{2}\phi-\alpha\nabla.(\underline{D}.\nabla\phi)$
(44)

where $\alpha=\beta_{1}\beta_{2}^{-1}$.

A perturbative expansion in $\alpha$,
$\phi=\sum_{n=0}^{\infty}\alpha^{n}\phi_{(n)}$ gives for the term
$\phi_{(k)}$:

$-\nabla
H.\nabla\phi_{(k)}+\sum_{|n+m|=k}\nabla\phi_{(n)}.\nabla\phi_{(m)}+\sum_{|n+m+1|=k}\nabla\phi_{(n)}.\underline{D}.\nabla\phi_{(m)}-\nabla^{2}\phi_{(k)}-\nabla.(\underline{D}.\nabla\phi_{((k-1)})=0$
(45)

for $k>0$ and where $\phi_{(0)}=H$ which leads to the first order
perturbation:

$(\nabla^{2}-\nabla H.\nabla)\phi_{(1)}=\nabla.(\underline{D}.\nabla H)+\nabla
H.\underline{D}.\nabla H$ (46)

It can be seen through repeated application of 45 that the Green’s function is
of the form specified in the main text for every perturbation.

##  Appendix B Evaluation of harmonic path integral for non-conservative
system

We use the generic solution to a path integral with a harmonic Lagrangian for
a system with a Hamiltonain given by
$H_{0}(\mathbf{x})=\frac{1}{2}\mathbf{x}.\underline{\kappa}.\mathbf{x}$, which
is given through the following formula[47]:

$\int_{\mathbf{r}(0)=\mathbf{x}^{\prime}}^{\mathbf{r}(t)=\mathbf{x}^{\prime\prime}}\\!\\!\\!\\!\\!\\!\mathcal{D}\left[\mathbf{r}(t)\right]e^{\frac{i}{\hbar}\int_{0}^{t}\mathrm{d}t^{\prime}L(\dot{x},x)}=\left(\frac{1}{2\pi
i\hbar}\right)^{D/2}\sqrt{\det\left(-\frac{\partial^{2}S_{cl}[\mathbf{x}^{\prime\prime},\mathbf{x}^{\prime}]}{\partial
x_{a}^{\prime\prime}\partial
x_{b}^{\prime}}\right)}\exp\left(\frac{i}{\hbar}S_{cl}[\mathbf{x}^{\prime\prime},\mathbf{x}^{\prime}]\right)$
(47)

where we perform the wick rotation $\hbar\to-i$ (note the lack of Boltzmann
factor here, it is in fact not introduced into this analysis due to the
fictitious nature of the time variable we have). Our Lagrangian is given by:

$L=\frac{1}{4}\dot{\mathbf{x}}^{2}+\frac{1}{4}\mathbf{x}.\underline{\kappa}.\underline{\kappa}.\mathbf{x}-\frac{1}{2}\text{Tr}(\underline{\kappa})$
(48)

The quantity $S_{cl}$ is the classical action evaluation along the classical
solution $\mathbf{x}_{cl}$ satisfying the boundary conditions
$\mathbf{x}_{cl}(0)=\mathbf{x}^{\prime}$ and
$\mathbf{x}_{cl}(s)=\mathbf{x}^{\prime\prime}$

The Lagrangian leads to the following Euler Lagrange equation:

$\frac{1}{2}\ddot{\mathbf{x}}=\frac{1}{2}\underline{\kappa}.\underline{\kappa}.\mathbf{x}$
(49)

which has the generic solution:

$\mathbf{x}_{cl}(\tau)=e^{\underline{\kappa}\tau}\mathbf{c}_{1}+e^{-\underline{\kappa}\tau}\mathbf{c}_{2}$
(50)

where we are free to choose the value of $v_{0}$ to satisfy the boundary
conditions. This can be done so, yielding the following classical path:

$\mathbf{x}_{cl}(\tau)=\text{csch}(\underline{\kappa}s)(\mathbf{x}^{\prime\prime}\sinh(\underline{\kappa}\tau)+\mathbf{x}^{\prime}\sinh(\underline{\kappa}(s-\tau)))$
(51)

The evaluation of the classical action evaluated along this path can be
performed (using the fact that all the hyperbolic functions are power series
in $\underline{\kappa}$, and thus commute):

$S_{cl}=\frac{1}{4}\mathbf{x}^{\prime}.\underline{\kappa}\text{coth}\left(\underline{\kappa}s\right).\mathbf{x}^{\prime}+\frac{1}{4}\mathbf{x}^{\prime\prime}.\underline{\kappa}\text{coth}\left(\underline{\kappa}s\right).\mathbf{x}^{\prime\prime}-\frac{1}{2}\mathbf{x}^{\prime\prime}.\underline{\kappa}\text{csch}\left(\underline{\kappa}s\right).\mathbf{x}^{\prime}-\frac{1}{2}\text{Tr}(\underline{\kappa})s$
(52)

this allows us also the evaluate the Von Vleck determinant, which is given by:

$\sqrt{\det\left(-\frac{\partial^{2}S_{cl}[\mathbf{x}^{\prime\prime},\mathbf{x}^{\prime}]}{\partial
x_{a}^{\prime\prime}\partial
x_{b}^{\prime}}\right)}=\sqrt{\det\left(\frac{1}{2}\underline{\kappa}\text{csch}\left(\underline{\kappa}s\right)\right)}$
(53)

Finally, we can write the full expression for
$K(\mathbf{x}^{\prime\prime},\mathbf{x}^{\prime})$:

$\displaystyle K(\mathbf{x}^{\prime\prime},\mathbf{x}^{\prime})$
$\displaystyle=\int_{0}^{\infty}\mathrm{d}s\left(\frac{1}{2\pi}\right)^{D/2}\sqrt{\det\left(\frac{1}{2}\underline{\kappa}\text{csch}\left(\underline{\kappa}s\right)\right)}\exp\left(-S_{cl}\right)$
(54) $\displaystyle K(\mathbf{x}^{\prime\prime},\mathbf{x}^{\prime})$
$\displaystyle=\int_{0}^{\infty}\mathrm{d}s\left(\text{PFS}(s)\exp\left(-S_{cl}\right)\right)$
(55)

where
$\text{PFS}(s)=\left(\frac{1}{2\pi}\right)^{D/2}\sqrt{\det\left(\frac{1}{2}\underline{\kappa}\text{csch}\left(\underline{\kappa}s\right)\right)}$
Furthermore, the prefactors can be added leading to the full Green’s function:

$G(\mathbf{x}^{\prime\prime},\mathbf{x}^{\prime})=\int_{0}^{\infty}\mathrm{d}s\text{PFS}(s)\times\exp\left(-S_{cl}+\frac{1}{4}\mathbf{x}^{\prime\prime}.\underline{\kappa}.\mathbf{x}^{\prime\prime}-\frac{1}{4}\mathbf{x}^{\prime}.\underline{\kappa}.\mathbf{x}^{\prime}\right)$
(56)

As an example, we can use standard Gaussian integrals to find the perturbation
associated with a Harmonic source term (such as would arise from the action of
different temperature baths)
$\left(\mathbf{x}^{\prime}.\underline{B}.\mathbf{x}^{\prime}\right)$. We
denote the difference of the log probability from equilibrium as $\delta(\log
P)$. We have to calculate the following integral:

$\delta(\log
P)\sim\int_{0}^{\infty}\mathrm{d}s\text{PFS}(s)\int\mathrm{d}\mathbf{x}^{\prime}\exp\left(-\frac{1}{2}\mathbf{x}^{\prime}.(\underline{M}+\frac{1}{2}\underline{\kappa}).\mathbf{x}^{\prime}-\frac{1}{2}\mathbf{x}^{\prime\prime}.(\underline{M}-\frac{1}{2}\underline{\kappa}).\mathbf{x}^{\prime\prime}+\mathbf{b}.\mathbf{x}^{\prime}\right)\left(\mathbf{x}^{\prime}.\underline{B}.\mathbf{x}^{\prime}\right)$
(57)

where

$\displaystyle\underline{M}(s)$
$\displaystyle=\frac{1}{2}\underline{\kappa}\text{coth}\left(\underline{\kappa}s\right)$
(58) $\displaystyle\mathbf{b}(s)$
$\displaystyle=\mathbf{x^{\prime\prime}}\frac{1}{2}\underline{\kappa}\text{csch}\left(\underline{\kappa}s\right)$
(59)

By completing the square, we can obtain the following equivalent expression
for the exponential:

$\displaystyle\delta(\log
P)\sim\int_{0}^{\infty}\mathrm{d}s\text{PFS}(s)\int\mathrm{d}\mathbf{x}^{\prime}\exp\bigg{(}-$
$\displaystyle\frac{1}{2}\left(\mathbf{x}^{\prime}-(\underline{M}(s)+\frac{1}{2}\underline{\kappa})^{-1}\mathbf{b}(s))\right).(\underline{M}(s)+\frac{1}{2}\underline{\kappa}).\left(\mathbf{x}^{\prime}-(\underline{M}(s)+\frac{1}{2}\underline{\kappa})^{-1}\mathbf{b}(s))\right)$
(60)
$\displaystyle+\frac{1}{2}\mathbf{b}(s).(M(s)+\frac{1}{2}\kappa)^{-1}.\mathbf{b}(s)-\frac{1}{2}\mathbf{x}^{\prime\prime}.(\underline{M}(s)-\frac{1}{2}\underline{\kappa}).\mathbf{x}^{\prime\prime}\bigg{)}$

The substitution
$\mathbf{y}(s)=\mathbf{x}^{\prime}-(\underline{M}(s)+\frac{1}{2}\underline{\kappa})^{-1}\mathbf{b}(s)$
allows us to compute the integral in equation 57. This leads to the following
expression:

$\displaystyle\delta(\log
P)\sim\int_{0}^{\infty}\mathrm{d}s\text{PFS}(s)\int\mathrm{d}\mathbf{y}(s)$
$\displaystyle\exp\left(-\frac{1}{2}\mathbf{y}(s).(\underline{M}(s)+\frac{1}{2}\underline{\kappa}).\mathbf{y}(s)+\frac{1}{2}\mathbf{b}(s).(\underline{M}(s)+\frac{1}{2}\underline{\kappa})^{-1}.\mathbf{b}-\frac{1}{2}\mathbf{x}^{\prime\prime}.(\underline{M}(s)-\frac{1}{2}\underline{\kappa}).\mathbf{x}^{\prime\prime}\right)$
(61)
$\displaystyle\times\left(\mathbf{y}(s).\underline{B}.\mathbf{y}(s)+\mathbf{b}(s).(\underline{M}(s)+\frac{1}{2}\underline{\kappa})^{-1}.\underline{B}.(\underline{M}(s)+\frac{1}{2}\underline{\kappa})^{-1}\mathbf{b})\right)$

Where we have suppressed the terms linearly dependent on $\mathbf{y}$ outside
the exponential, as these will integrate out to zero. This integral can be
performed over $\mathbf{y}$ by using the standard properties of Gaussian
integrals, leading to the following expression:

$\displaystyle\sim\int_{0}^{\infty}\mathrm{d}s\text{PFS}(s)\left(\frac{(2\pi)^{N}}{\det(\underline{M}(s)+\frac{1}{2}\underline{\kappa})}\right)^{1/2}\left(\text{Tr}\left((\underline{M}(s)+\frac{1}{2}\underline{\kappa})^{-1}.\underline{B}\right)+\mathbf{b}(s).(\underline{M}(s)+\frac{1}{2}\underline{\kappa})^{-1}.\underline{B}.(\underline{M}(s)+\frac{1}{2}\underline{\kappa})^{-1}\mathbf{b}(s))\right)$
(62)
$\displaystyle\times\exp\left(\frac{1}{2}\mathbf{b}(s).(\underline{M}(s)+\frac{1}{2}\underline{\kappa})^{-1}.\mathbf{b}(s)-\frac{1}{2}\mathbf{x}^{\prime\prime}.(\underline{M}(s)-\frac{1}{2}\underline{\kappa}).\mathbf{x}^{\prime\prime}\right)$

This expression is then multiplied by the von Vleck determinant and other
prefactors whose dependence we suppresssd earlier to lead to the full
resultant expression for a Harmonic perturbation, once the integral over all
times is performed. Though this expression appears difficult, remarkable
simplifications occur. After all the experssions are reinserted, the term in
the exponent is actually equal to zero.

The von vleck determinant $\text{PFS}(s)$ multiplied by the determinant
arising from the integral is given by:

$\left(\frac{1}{2\pi}\right)^{N/2}\sqrt{\det\left(\frac{1}{2}\underline{\kappa}\text{csch}\left(\underline{\kappa}s\right)\right)}*\left(\frac{(2\pi)^{N}}{\det(\underline{M}+\frac{1}{2}\underline{\kappa})}\right)^{1/2}$
(63)

The hyperbolic expressions can again be simplified, leading to:

$\left(\frac{1}{2\pi}\right)^{N/2}\sqrt{\det\left(\frac{1}{2}\underline{\kappa}\text{csch}\left(\underline{\kappa}s\right)\right)}*\left(\frac{(2\pi)^{N}}{\det(\underline{M}+\frac{1}{2}\underline{\kappa})}\right)^{1/2}=\sqrt{\det\left(\exp(-\underline{\kappa}s)\right)}=\exp(-\text{Tr}(\underline{\kappa})s/2)$
(64)

Where we use the property that the determinant of the matrix exponential being
equal to the exponential of the trace. We can observe that this term cancels
out with the term in the Lagrangian proportional to
$\text{Tr}(\underline{\kappa})$

Further simplifications exist:

$\mathbf{b}.(\underline{M}+\frac{1}{2}\underline{\kappa})^{-1}.\underline{B}.(\underline{M}+\frac{1}{2}\underline{\kappa})^{-1}\mathbf{b})=\mathbf{x}^{\prime\prime}\exp(-\underline{\kappa}s)\underline{B}\exp(-\underline{\kappa}s)\mathbf{x}^{\prime\prime}$
(65)

The term
$\text{Tr}\left((\underline{M}(s)+\frac{1}{2}\underline{\kappa})^{-1}.\underline{B}\right)$
appearing in equation 62 doesn’t have any dependence on
$\mathbf{x}^{\prime\prime}$. In other words, it would change the probability
uniformly up or down, however, this dependence disappears when we divide
through by the partition function. It is important to bear this term in mind
for higher order corrections, however.

We therefore have the following expression:

$\int_{0}^{\infty}\mathrm{d}s\mathbf{x}^{\prime\prime}\exp(-\underline{\kappa}s)\underline{B}\exp(-\underline{\kappa}s)\mathbf{x}^{\prime\prime}$
(66)

for the modification in the probability. As seen in the main text.

###  B.1 Integration over non-harmonic sources

The preceding discussion can be applied to source functions of anharmonic form
fairly straightforwardly. If we use the facts that we saw in the preceding
discussion, that

$(\underline{M}(s)+\frac{1}{2}\underline{\kappa})^{-1}.\mathbf{b}(s)=\exp(-\underline{A}s)\mathbf{x}^{\prime\prime}\equiv\bm{\chi}(s)$
(67)
$\left(\frac{1}{2}\mathbf{b}(s).(\underline{M}(s)+\frac{1}{2}\underline{\kappa})^{-1}.\mathbf{b}(s)-\frac{1}{2}\mathbf{x}^{\prime\prime}.(\underline{M}(s)-\frac{1}{2}\underline{\kappa}).\mathbf{x}^{\prime\prime}\right)=0$
(68)
$(\underline{M}(s)+\frac{1}{2}\underline{\kappa})^{-1}=\left(\beta\underline{\kappa})^{-1}.(\underline{I}-\exp(-2\beta
s\underline{\kappa})\right)\equiv\underline{\mu}(s)$ (69)

We can take an arbitrary function of $\mathbf{x}^{\prime}$ at equation 61
$f(\mathbf{x}^{\prime})$, giving:

$\delta(\log
P)\sim\int_{0}^{\infty}\mathrm{d}s\text{PFS}(s)\int\mathrm{d}\mathbf{y}(s)\exp\left(-\frac{1}{2}\mathbf{y}(s).\underline{\mu}(s)^{-1}.\mathbf{y}(s)\right)f\left(\mathbf{y}(s)+\exp(-\underline{A}s)\mathbf{x}^{\prime\prime}\right)$
(70)

where such formulas can be evaluated using the identity[48]:

$\int\mathrm{d}\mathbf{y}(s)\exp\left(-\frac{1}{2}\mathbf{y}(s).\underline{\mu}(s)^{-1}.\mathbf{y}(s)\right)f\left(\mathbf{y}(s)+\bm{\chi}(s)\right)=\left(\sqrt{\frac{(2\pi)}{\det(\underline{\mu}(s)^{-1})}}\right)^{n}\times\exp\left(\frac{1}{2}\mu(s)_{i,j}\frac{\partial}{\partial
y_{i}}\frac{\partial}{\partial
y_{j}}\right)f\left(\mathbf{y}(s)+\bm{\chi}(s)\right)\bigg{|}_{\mathbf{y}=0}$
(71)

for a generic function $f(\mathbf{y})$ equal to some power series in
$\mathbf{y}$, the above expression can be expanded to yield the sets of
indices necessary for the summations in the main text.

##  Appendix C Anharmonic path integral for non-conservative system

In order to calculate the effect of arbitrary perturbations, we must include
the effect of anharmonicity, much of this discussion bears similarity with
similar problems arising in quantum field theory[48], we recall that the
definition of the Green’s function is given by:

$\displaystyle G(\mathbf{x}^{\prime\prime},\mathbf{x}^{\prime})$
$\displaystyle=\int_{0}^{\infty}\mathrm{d}s\int_{\mathbf{r}(0)=\mathbf{x}^{\prime}}^{\mathbf{r}(s)=\mathbf{x}^{\prime\prime}}\\!\\!\\!\\!\\!\\!\mathcal{D}\left[\mathbf{r}(s)\right]e^{-\int_{0}^{s}\mathrm{d}s^{\prime}L(\dot{x},x)}$
(72)
$\displaystyle=e^{\frac{\beta}{2}\left(H(x)-H(x^{\prime})\right)}\int_{0}^{\infty}\mathrm{d}s\int_{\mathbf{r}(0)=\mathbf{x}^{\prime}}^{\mathbf{r}(s)=\mathbf{x}^{\prime\prime}}\\!\\!\\!\\!\\!\\!\mathcal{D}\left[\mathbf{r}(s)\right]e^{-\int_{0}^{s}\mathrm{d}s^{\prime}L(\dot{x},x)}$
(73)

Our difficulty lies in evaluation of the path integral in the presence of
anharmonic conservative forces in $L$. One way beyond this difficult is
splitting the action into a baseline contribution which we can evaluate, and
another contribution which we cannot.
$\int_{0}^{s}\mathrm{d}s^{\prime}L(\dot{x},x)=S_{h}+S_{a}$. We choose the
baseline to be the harmonic path integral which we have just evaluated in the
previous section. Then, conducting a Taylor series in $S_{a}$ we obtain:

$G(\mathbf{x}^{\prime\prime},\mathbf{x}^{\prime})=e^{\frac{\beta}{2}\left(H(x)-H(x^{\prime})\right)}\int_{0}^{\infty}\mathrm{d}s\sum_{n=0}^{\infty}\frac{1}{n!}\int_{\mathbf{r}(0)=\mathbf{x}^{\prime}}^{\mathbf{r}(s)=\mathbf{x}^{\prime\prime}}\\!\\!\\!\\!\\!\\!\mathcal{D}\left[\mathbf{r}(s)\right]e^{-S_{h}}(-S_{a})^{n}$
(74)

In order to evaluate the path integrals in eq. 74 we introduce the following
Lagrangian:

$L_{J}=\int_{0}^{s}\mathrm{d}s^{\prime}\frac{1}{4}\dot{\mathbf{x}(s^{\prime})}^{2}+\frac{1}{4}\mathbf{x}(s^{\prime}).\underline{\kappa}.\underline{\kappa}.\mathbf{x}(s^{\prime})-\frac{1}{2}\text{Tr}(\underline{\kappa})+\mathbf{J}(s^{\prime}).\mathbf{x}(s^{\prime})$
(75)

By taking repeated functional integrals with respect to $J$, we can reproduce
the terms appearing in eq. 74. Crucially, with the Lagrangian given in 75, the
path integral can be evaluated analytically, giving:

$K(\mathbf{x}^{\prime\prime},\mathbf{x}^{\prime},s)=\int_{\mathbf{r}(0)=\mathbf{x}^{\prime}}^{\mathbf{r}(s)=\mathbf{x}^{\prime\prime}}\\!\\!\\!\\!\\!\\!\mathcal{D}\left[\mathbf{r}(s)\right]e^{-\int_{0}^{s}\mathrm{d}s^{\prime}L_{J}(\dot{x},x)}=K_{h}(\mathbf{x}^{\prime\prime},\mathbf{x}^{\prime},s)\exp(-W[\mathbf{J}])$
(76)

Where we have introduced the Green’s function neglecting endpoint Hamiltonian
values:

$K_{h}(\mathbf{x}^{\prime\prime},\mathbf{x}^{\prime},s)=\left(\frac{1}{2\pi}\right)^{D/2}\sqrt{\det\left(\frac{1}{2}\underline{\kappa}\text{csch}\left(\underline{\kappa}s\right)\right)}\exp\left(-S_{cl}\right)$
(77)

where $S_{cl}$ is the same classical action appearing in the previous section,
that of a $J=0$ harmonic system evaluated along the classical path
$\mathbf{x}_{cl}(s^{\prime})$. The Schwinger functional, which accounts for
the terms that depend on $J$ is given by:

$W[J]=\int\mathrm{d}s\mathbf{J}(s).\mathbf{x}_{cl}(s)+\frac{1}{2}\int\mathrm{d}s\int\mathrm{d}s^{\prime}\mathbf{J}(s).\underline{G}_{D}(s,s^{\prime}).\mathbf{\underline{missing}}\mathbf{J}(s^{\prime})$
(78)

Where $G_{D}$ is (unfortunately) usually also called a many-body Green’s
function, which for a harmonic system is given by:

$\displaystyle\underline{G}_{D}(s,s^{\prime})=\frac{1}{2}\underline{\kappa}^{-1}\text{csch}(\underline{\kappa}s)\bigg{(}$
$\displaystyle\theta\left(s-s^{\prime}\right)\left(\cosh\left(\underline{\kappa}\left(s-s^{\prime}-s\right)\right)-\cosh\left(\underline{\kappa}\left(s+s^{\prime}-s\right)\right)\right)$
(79)
$\displaystyle+\theta\left(s^{\prime}-s\right)\left(\cosh\left(\underline{\kappa}\left(s-s^{\prime}-s\right)\right)-\cosh\left(\underline{\kappa}\left(s-s^{\prime}+s\right)\right)\right)\bigg{)}$

We can calculate an arbitrary perturbation accounting for a potential term
going as $V$ according to the formula:

$G(\mathbf{x}^{\prime\prime},\mathbf{x}^{\prime})=e^{\frac{\beta}{2}\left(H(\mathbf{x}^{\prime\prime})-H(\mathbf{x}^{\prime})\right)}\int_{0}^{\infty}\mathrm{d}sK_{h}(\mathbf{x}^{\prime\prime},\mathbf{x}^{\prime},s)\exp\left[\int\mathrm{d}s^{\prime}V\left(\frac{\delta}{\delta\mathbf{J}(s^{\prime})}\right)\right]\exp(-W[\mathbf{J}])\big{|}_{\mathbf{J}=0}$
(80)

splitting the Hamiltonian into harmonic and anharmonic terms gives us the
expression in the main text, which goes as:

$G(\mathbf{x}^{\prime\prime},\mathbf{x}^{\prime})=e^{\frac{\beta}{2}\left(V(\mathbf{x}^{\prime\prime})-V(\mathbf{x}^{\prime})\right)}\int_{0}^{\infty}\mathrm{d}sG_{h}(\mathbf{x}^{\prime\prime},\mathbf{x}^{\prime},s)\exp\left[\int\mathrm{d}s^{\prime}V_{a}\left(\frac{\delta}{\delta\mathbf{J}(s^{\prime})}\right)\right]\exp(-W[\mathbf{J}])\big{|}_{\mathbf{J}=0}$
(81)

Where we already know the mechanical trick of repetively applying
$G_{h}(\mathbf{x}^{\prime\prime},\mathbf{x}^{\prime})$ to any function $\rho$,
as was expressed in the previous section. The term $V_{a}$ can potentially
contain terms going as order $x^{2}$, but it’s best to keep them in the
perturbative expansion so that $G_{h}$ can be applied simply each time.

###  C.1 Application to active brownian particle in anharmonic potential

We can calculate the application of the above for the active anharmonic system
described in the main text. Such a system has
$\rho(\mathbf{x})=\beta^{2}aq^{(3)}x^{3}A$ where the full state
$\mathbf{x}=(x,A)$ For brevity let us define $q=q^{(3)}$. For ease of
nomenclature, we shall define the integral over the harmonic green’s function
as an operator $\hat{O}$:

$\int_{0}^{\infty}\mathrm{d}s\int\mathrm{d}\mathbf{x}^{\prime}G_{h}(\mathbf{x},\mathbf{x}^{\prime},s)\rho(\mathbf{x}^{\prime})=\hat{O}\rho$
(82)

for this system $V(\mathbf{x})=\frac{1}{4}qx^{4}$ and
$V_{a}(\mathbf{x})=3qx^{2}+q^{2}x^{6}$

The first terms in the anharmonic series then go as:

$\hat{O}[\rho]+\hat{O}\left[\beta/2\right(V(\mathbf{x})-V(\mathbf{x}^{\prime}))\rho]+\hat{O}\left[\left(\int\mathrm{d}s^{\prime}V_{a}\left(\frac{\delta}{\delta\mathbf{J}(s^{\prime})}\right)\exp(-W[\mathbf{J}])\big{|}_{\mathbf{J}=0}\right)\rho\right]$
(83)

the action of an operator $B_{ij}\partial_{i}\partial_{j}$ on the term
$\exp(-W[J])$ is
$\int_{0}^{s}\mathrm{d}s^{\prime}\text{Tr}(\underline{G}_{D}(s^{\prime},s^{\prime}).\underline{B})+\frac{1}{2}\mathbf{x}_{cl}(s^{\prime}).\underline{B}.\mathbf{x}_{cl}(s^{\prime})$.
Using the mechanical rules, each of these terms can be evaluated, for now let
us see how terms going as $\mathbf{x}^{2}$ look:

$\displaystyle\hat{O}\rho$
$\displaystyle=\frac{3aq\tau_{a}^{2}}{\gamma^{2}}xA$ (84)
$\displaystyle\hat{O}(\beta/2)V(\mathbf{x})\rho)$
$\displaystyle\sim\mathcal{O}(x^{4})\text{ as $\hat{O}$ doesn't act on $x$ and
$V\sim x^{4}$}$ (85)
$\displaystyle\hat{O}(\beta/2)V(\mathbf{x}^{\prime})\rho)$
$\displaystyle=\frac{630aq^{2}\tau_{a}^{4}}{\beta\gamma^{4}}xA$ (86)

and to evaluate the Green’s terms we use the fact that
$\int_{0}^{s}\mathrm{d}s^{\prime}\text{Tr}(\underline{G}_{D}(s^{\prime},s^{\prime}).\underline{B})=-(1/2)qs^{2}\beta$
where in this example $B=\begin{bmatrix}&3q,0\\\ &0,0\end{bmatrix}$. The
classical path, which depends on both $\mathbf{x}^{\prime}$ and
$\mathbf{x}^{\prime\prime}$ is given in (51). Taken together, these can be
evaluated, yielding:

$\displaystyle\hat{O}\left(\int_{0}^{s}\mathrm{d}s^{\prime}\text{Tr}(\underline{G}_{D}(s^{\prime},s^{\prime}).\underline{B})\rho\right)$
$\displaystyle=\frac{18aq^{2}\tau_{a}^{4}}{\beta\gamma^{4}}xA$ (87)
$\displaystyle\hat{O}\left(\int_{0}^{s}\mathrm{d}s^{\prime}\frac{1}{2}\mathbf{x}_{cl}(s^{\prime}).\underline{B}.\mathbf{x}_{cl}(s^{\prime})\rho\right)$
$\displaystyle=\frac{432aq^{2}\tau_{a}^{4}}{\beta\gamma^{4}}xA$ (88)

From which the subequent structure of the series seen in the main text becomes
apparent. A similar process can be repeated for terms going as
$\sim\mathbf{x}^{4}$ etc.

However, a problem arises when the anharmonic series diverges (as is usually
the case). In such a scenario, it is better not to perform a series expansion
and instead calculate the following integrals directly, which leads to the
anharmonic integral and anharmonic integral with dissipative correction
approximations respectively:

$\int_{0}^{\infty}\mathrm{d}s\int\mathrm{d}\mathbf{x}^{\prime}e^{(\beta/2)(V(\mathbf{x})-V(\mathbf{x}^{\prime}))}G_{h}(\mathbf{x},\mathbf{x}^{\prime},s)\rho(\mathbf{x}^{\prime})$
(90)
$\int_{0}^{\infty}\mathrm{d}s\int\mathrm{d}\mathbf{x}^{\prime}e^{(\beta/2)(V(\mathbf{x})-V(\mathbf{x}^{\prime}))}G_{h}(\mathbf{x},\mathbf{x}^{\prime},s)\rho(\mathbf{x}^{\prime})\left(1+\int_{0}^{s}\mathrm{d}s^{\prime}\left[\text{Tr}(\underline{G}_{D}(s^{\prime},s^{\prime}).\underline{B}))+\frac{1}{2}\mathbf{x}_{cl}(s^{\prime},\mathbf{x},\mathbf{x}^{\prime}).\underline{B}.\mathbf{x}_{cl}(s^{\prime},\mathbf{x},\mathbf{x}^{\prime})\right]\right)$
(91)

both of these integrals can be performed ”analytically” for an active brownian
system with the help of symbolic software such as Mathematica, however the
resulting equations are too long to be reproduced here, but are available upon
request.

##  Appendix D Mode space

Some generic features of the modification of equilibrium systems to non-
equilibrium perturbations can be discerned if we analyze the microscopic
probability in terms of the modes. This is easily done for a harmonic system.
Recalling that the log probability in these systems goes as:

$\log(P(\mathbf{x}))=-\frac{1}{2}\beta\mathbf{x}.\underline{\kappa}.\mathbf{x}-\frac{1}{2}\beta^{2}\mathbf{x}.\int_{0}^{\infty}\mathrm{d}s\exp(-\beta\underline{\kappa}s).\underline{C}.\exp(-\beta\underline{\kappa}s).\mathbf{x}+c_{1}$
(92)

Where $\underline{C}$ is some dissipation matrix and $c_{1}$ is a
normalization constant.

The matrix $\underline{\kappa}$ will have certain properties. For the system
to be conservative in the absence of any non-equilibrium perturbation, the
matrix $\kappa$ must be symmetric $\underline{\kappa}=\underline{\kappa}^{T}$.
A symmetric matrix can always be diagonalized
$\underline{\kappa}=\underline{Q}^{-1}\underline{L}\underline{Q}$ where the
matrix $\underline{L}$ will be a diagonal matrix with each diagonal entry
being equal to one of the eigenvalues of the original matrix
$L_{ii}=\lambda_{i}$ and the matrix $\underline{Q}$ is the matrix of
eigenvectors. The matrix exponential for diagonalizable matrices is also well
known to have the simplification
$e^{\underline{\kappa}}=\underline{Q}^{-1}e^{\underline{L}}\underline{Q}$. The
matrix $Q$ defines a linear transformation, therefore we can define a new
space for our problem given by $\mathbf{\tilde{x}}=\underline{Q}\mathbf{x}$.
Each entry in the vector $\tilde{\mathbf{x}}$ is a mode of the equilibrium
system. This allows re-expression of eq. 92 as:

$\log(P(\mathbf{\tilde{x}}))=-\frac{1}{2}\beta\mathbf{\tilde{x}}.\underline{L}.\mathbf{\tilde{x}}-\frac{1}{2}\beta^{2}\mathbf{\tilde{x}}.\int_{0}^{\infty}\mathrm{d}t\exp(-\beta\underline{L}t).\underline{\tilde{C}}.\exp(-\beta\underline{L}t).\mathbf{\tilde{x}}$
(93)

where $\underline{\tilde{C}}=\underline{Q}.\underline{C}.\underline{Q}^{-1}$.
In mode space, the integral over all path lengths can be performed, leaving
the final probability as:

$\log(P(\mathbf{\tilde{x}}))=-\frac{1}{2}\beta\mathbf{\tilde{x}}.\underline{L}.\mathbf{\tilde{x}}-\frac{1}{2}\beta^{2}\mathbf{\tilde{x}}.\underline{\tilde{D}}.\mathbf{\tilde{x}}$
(94)

where we have a new matrix $\underline{\tilde{D}}$ whose elements are given
by:

$\tilde{D}_{ij}=\frac{\tilde{C}_{ij}}{\lambda_{i}+\lambda_{j}}$ (95)

As the matrix $C$ for driving under a non-conservative force goes as
$\underline{C}=\underline{\kappa}.\underline{M}=\underline{Q}^{-1}.\underline{L}.\underline{Q}.\underline{M}$,
which then means
$\tilde{\underline{C}}=\underline{L}\underline{Q}.\underline{M}.\underline{Q}^{-1}$
and defining
$\tilde{\underline{M}}=\underline{Q}.\underline{M}.\underline{Q}^{-1}$

$\tilde{D}_{ij}=\frac{\lambda_{i}\tilde{M}_{ij}}{\lambda_{i}+\lambda_{j}}$
(96)

we get the property that the non-conservative force affects the variance of
each mode by an amount that is roughly independent of eigenvalue of the mode.

For an active polymer, half of each ”modes” refers to the degrees of freedom
of the active forces $A$. When integrating out over these degrees of freedom,
we can obtain the effective $1/(1+\lambda^{2})$ expression that is seen in the
main text. The true expressions can be found analytically, but are rather
long.

##  Appendix E Simuation details

The simulated harmonic system of $N$ particles where the form of potential
between particles $i$ and $j$ is given by:

$\phi(\mathbf{x}_{i},\mathbf{x}_{j})=\frac{1}{2}k(\mathbf{x}_{i}-\mathbf{x}_{j})^{2}$
(97)

where $k$ is the spring constant. For simplicity, we consider a system where
particles are only bound to their 1D neighbors, leading to a polymer. The
system is thermalized through the Langevin thermostat, where the noise is
uncorrelated in time:

$\langle\bm{\xi}_{i}(t)\bm{\xi}_{i}(t^{\prime})\rangle=2k_{b}T_{i}\gamma\delta(t-t^{\prime})$
(98)

where $k_{b}T$ is the temperature. For the cases where an anharmonic potential
is introduced, we modify the relationship in (99) to

$\Phi(\mathbf{x}_{i},\mathbf{x}_{j})=\frac{1}{2}k(\mathbf{x}_{i}-\mathbf{x}_{j})^{2}+\frac{1}{4}q(\mathbf{x}_{i}-\mathbf{x}_{j})^{4}$
(99)

and we integrate numerically the following equations

$\displaystyle\gamma\frac{\mathrm{d}\mathbf{x}_{i}}{\mathrm{d}t}$
$\displaystyle=\mathbf{F}_{i}(\\{\mathbf{x}\\})+a\mathbf{A}_{i}+\bm{\xi}_{i}(t)$
(100) $\displaystyle\gamma\frac{\mathrm{d}\mathbf{A}_{i}}{\mathrm{d}t}$
$\displaystyle=-(\gamma/\tau_{a})\mathbf{A}_{i}+\bm{\xi}_{i}(t)$ (101)

where
$\mathbf{F}_{i}(\\{\mathbf{x}\\})=-\nabla_{i}\sum_{j}^{N-1}\Phi(\mathbf{x}_{j},\mathbf{x}_{j+1})$