# Statistical mechanics of passive Brownian particles in a fluctuating
harmonic trap

Derek Frydel  Department of Chemistry, Universidad Técnica Federico Santa
María, Campus San Joaquin, Santiago, Chile

###### Abstract

We consider passive Brownian particles trapped in an "imperfect" harmonic
trap. The trap is imperfect because it is randomly turned off and on, and as a
result particles fail to equilibrate. In other words, the harmonic trap is
time-dependent on account of its strength evolving stochastically in time.
Particles in such a system are passive and activity arises through external
control of a trapping potential, therefore, no internal energy is used to
power particle motion. A stationary Fokker-Planck equation of this system can
be represented as a third-order differential equation, and its solution, a
stationary distribution, can be found by representing it as a superposition of
Gaussian distributions for different strengths of a harmonic trap. The
resulting stationary distributions tend to be spread out with an exponential
asymptotic tail. A pure Gaussian distribution is recovered only in the limit
of a very fast alternating on and off cycle. The average potential energy of
the system is found to be the same as that for a system in equilibrium, that
is, for a system with a time independent harmonic trap. The entropy production
rate, which quantifies the distance from equilibrium, is found to be
independent of the rate with which the trap fluctuates.

###### pacs:

##  I Introduction

A typical active particle system consists of particles that derive motion from
an internal energy source. Standard models that represent such particles are
the run-and-tumble (RTP) Tailleur08 ; Tailleur09 ; Evans18 ; Dhar18 ; Dhar19 ;
Basu20 ; Soto20 ; Razin20 ; Frydel21b ; Frydel21c ; Farago22 ; Frydel22b ;
Connor23 ; Farago24 ; Ybert24 ; Loewe24  and active-Brownian particle (ABP)
models Dhar20 ; Cargalio22 , one more suitable to biological and the other to
chemical systems. Another active scenario arises when a passive Brownian
particle is immersed in a bath of active particles. The model that addresses
this situation is the active Ornstein-Uhlenbeck particle model (AOUP) Szamel14
; Shankar18 ; Sevilla19 ; Cates21 ; Marconi18 ; Lowen22 .

A different and less explored design of active dynamics is realized by placing
passive Brownian particles in a fluctuating external potential, thereby
preventing particles from equilibrating. In this setup, there is no need for a
special type of particles and the only requirement is that an external
potential varies stochastically with time. A possible experimental realization
of such a system could be attained using tweezer instruments and techniques
Yael11 ; Brady16 ; Brady21 ; Lowen22b ; Genet22 ; Yael23 .

This work considers passive Brownian particles trapped in a harmonic potential
$u=Kr^{2}/2$ with time-dependent strength $K\equiv K(t)$. We consider a
specific evolution of $K(t)$ in which $K(t)$ changes discontinuously between
two discrete values, zero and another finite value of the strength, therefore,
the trap is either off or on. We refer to this external potential as
"imperfect" harmonic trap. The trap fluctuations are implemented by random
sampling of the times during which a potential is in a given state, $t_{p}$,
from an exponential probability distribution $\sim e^{-t_{p}/\tau}$. This
means that on average the the times during which a trap is on and off are the
same and what changes is the rate, given by $\tau^{-1}$, with which the trap
fluctuates between the two states.

A stationary situation of this system gives rise to a third-order differential
equation. A third-order differential equation has previously been found to
arise for run-and-tumble particles in a harmonic trap in one- and two-
dimensions Frydel22c ; Frydel23b . The solution was found to be a convolution
of two distributions, each one separately obeying a first-order differential
equation. The convolution form of the solution reflects the fact that the two
random processes involved, Brownian motion and active dynamics, remain
independent in a harmonic potentia Frydel24a .

In the present system, a solution of the corresponding third-order equation
can be represented as a superposition of Gaussian distributions for different
strengths of a harmonic trap. This decomposes the solution into two
contributions, a Gaussian distribution and the probability distribution of
effective trap strengths, each distribution satisfying a first-order
differential equation. The superposition of Gaussian distributions has been
known in quantum theory and financial markets Kleinert04 ; Kleinert08 . In
non-equilibrium statistical mechanics, those distributions have been
considered in Beck01  in the context of non-extensive statistical mechanics.

The current work could be viewed as related to the systems that explore
control of external forces and information feedback to perform work Sagawal15
; Sagawa08 ; Cao09 ; Horowitz10 ; Sagawa10 ; Pon10 ; Suzuki10 ; Kundu12 .
Although the model that is analyzed does not involve information feedback and
the only control comes from regulating the time duration when particles are
either trapped or released, it could be considered a small first step in
theoretical understanding of this class of problems. In a more general
context, by proposing a non-equilibrium model and solving it, the current work
contributes an item to the collection of exactly solvable non-equilibrium
problems.

This work is organized as follows. In Sec. (II) we introduce the model and
obtain a third-order differential equation for a stationary state. In Sec.
(III) we represent the solution as a superposition of Gaussian distributions
for different trap strengths K. In Sec. (IV) we consider particles trapped in
a harmonic potential with a general time-dependent trap strength. It is shown
that the time-dependent distribution at any given moment is a Gaussian
distribution. This result justifies the use of the superposition formula. In
Sec. (V) we extend all the exact results to a system for an arbitrary
dimension. In Sec. (VI) we consider quantities of physical interest. The work
is concluded in Sec. (VII).

##  II The model

The focus of this work is a conceptually simple model: an ideal gas trapped in
an "imperfect" harmonic potential. The potential is imperfect because it is
turned off and on at random time intervals. The alternate cycle of trapping
and releasing prevents particles from attaining equilibrium and as a result
gives rise to a non-equilibrium situation. The times $t_{p}$ during which a
particle is either trapped or released are sampled from the exponential
distribution $\sim e^{-t_{p}/\tau}$, where $\tau$ is the average time during
which a particle persists in a given state. At the end of each time $t_{p}$, a
particle switches to another state with probability one.

Since the times when a trap is in the "on" and "off" state are sampled from
the same distribution, the system spends on average the same amount of time in
both states. What changes is the rate, given by $\tau^{-1}$, with which the
trap fluctuates between the two states. The model represents a specific case
of a harmonic potential $u=Kr^{2}/2$ with the time-dependent strength $K\equiv
K(t)$.

For a system in one-dimension, the Fokker-Planck formulation that describes
such a model might be written as

$\displaystyle\dot{\rho}_{+}=\mu\left(\frac{K}{2}+\Delta\right)[x\rho_{+}]^{\prime}+D\rho_{+}^{\prime\prime}+\frac{1}{\tau}\left(\rho_{-}-\rho_{+}\right)$
$\displaystyle\dot{\rho}_{-}=\mu\left(\frac{K}{2}-\Delta\right)[x\rho_{-}]^{\prime}+D\rho_{-}^{\prime\prime}+\frac{1}{\tau}\left(\rho_{+}-\rho_{-}\right),$
(1)

where $\mu$ is the mobility, $D=\mu k_{B}T$ is the diffusion constant, $T$ is
the temperature, and $k_{B}$ is the Boltzmann constant. To simplify the
expressions, we use dot notation to represent the time derivative,
$\dot{\rho}_{+}\equiv\frac{\partial\rho_{+}}{\partial t}$, and the prime
notation to represent derivatives with respect to position,
$[x\rho_{+}]^{\prime}\equiv\frac{\partial}{\partial x}[x\rho_{+}]$ and
$\rho_{+}^{\prime\prime}\equiv\frac{\partial^{2}}{\partial x^{2}}$. $\rho_{+}$
and $\rho_{-}$ are the distributions of particles in a harmonic potential with
the respective strength $K/2+\Delta$ and $K/2-\Delta$. The last term on the
right-hand side of each equation represents the conversion of one type of
particle into another and provides coupling between the two equations and also
prevents particles to equilibrate. The remaining terms represent the usual
flux given by

$j_{\pm}=-\mu(K/2\pm\Delta)x\rho_{\pm}-D\rho_{\pm}^{\prime}.$

In this work we are interested in the situation $\Delta=K/2$, in which case
the Fokker-Planck equation becomes

$\displaystyle\dot{\rho}_{+}=\mu
K[x\rho_{+}]^{\prime}+D\rho_{+}^{\prime\prime}+\frac{1}{\tau}\left(\rho_{-}-\rho_{+}\right)$
$\displaystyle\dot{\rho}_{-}=D\rho_{-}^{\prime\prime}+\frac{1}{\tau}\left(\rho_{+}-\rho_{-}\right).$
(2)

The distribution $\rho_{-}$ represents unconstrained particles corresponding
to the trap being turned off, and $\rho_{+}$ represents particles subject to a
confining potential.

To simplify the notation, we introduce dimensionless parameters. The position
of particles on the $x$-axis and the rate at which a harmonic trap fluctuates
in dimensionless units are then given by

$z=x\sqrt{\frac{\mu K}{D}},~{}~{}~{}~{}~{}~{}~{}~{}\alpha=\frac{1}{\tau\mu
K}.$

At a stationary state and dimensionless parameters Eq. (2) becomes

$\displaystyle
0=[z\rho_{+}]^{\prime}+\rho_{+}^{\prime\prime}+\alpha\left(\rho_{-}-\rho_{+}\right)$
$\displaystyle
0=\rho_{-}^{\prime\prime}+\alpha\left(\rho_{+}-\rho_{-}\right).$ (3)

The two equations can be merged into a single differential equation for a
distribution $\rho=(\rho_{+}+\rho_{-})/2$. This results in a third-order
differential equation given below

$0=-\alpha z^{3}\rho+2(1-\alpha
z^{2})\rho^{\prime}-(2-z^{2})z\rho^{\prime\prime}+z^{2}\rho^{\prime\prime\prime}.$
(4)

By eliminating coupling, we increase order of the equation. See Appendix (A)
for details of the derivation.

In the limit of a fast appearing-disappearing trap, $\alpha\to\infty$, Eq. (4)
reduces to

$0=-\frac{z\rho}{2}-\rho^{\prime}.$ (5)

$\rho$ in this limit converges to a Gaussian function and it can be identified
with a system in equilibrium in a harmonic trap half the strength of the
original trap.

Finding the complete solution to Eq. (4) is more challenging on account of it
being third-order. While first- and second-order differential equations are
relatively common in physics, third-order differential equations are
encountered less frequently Pati13 . They are also more challenging to solve
and analyze and require more advanced and creative approach. However, before
solving Eq. (4), it is possible to analyze it with other techniques.

Moments of a stationary distribution $\rho$ can be obtained directly from Eq.
(4), without solving it Frydel22c ; Frydel23b ; Frydel24a . To do this, we
transform Eq. (4) into a recurrence relation by operating on it with
$\int_{-\infty}^{\infty}dz\,z^{2n-3}$. Followed by the integration by parts
this yields

$m_{n}=\frac{4n-2}{\alpha}\bigg{[}(n+\alpha)m_{n-1}-n(2n-3)m_{n-2}\bigg{]},$
(6)

where $m_{n}=\langle z^{2n}\rangle$ are even moments of $\rho$, all odd
moments being zero.

Given the initial condition $m_{0}=1$, all subsequent moments can be
calculated, and the first two terms of the sequence are

$\displaystyle\langle z^{2}\rangle=\left(2+\frac{1}{\alpha}\right)$
$\displaystyle\langle
z^{4}\rangle=6\left(2+\frac{3}{\alpha}+\frac{2}{\alpha^{2}}\right).$ (7)

In the limit $\alpha=\infty$, all the moments reduce to the moments of a
Gaussian distribution $\sim e^{-z^{2}/4}$. The moments increase as $\alpha$
decreases, indicating the spreading out of $\rho$.

Moments for stationary distributions of particles in a different state,
$\rho_{\pm}$, can be calculated by operating on the two equations in Eq. (3)
with $\int dz\,z^{2n+2}$. After some manipulation, the details of which are
given in Appendix B, we get two coupled recurrence relations, from which we
get

$\displaystyle\langle z^{2}\rangle_{+}=2$ $\displaystyle\langle
z^{4}\rangle_{+}=6\left(2+\frac{1}{\alpha}\right),$ (8)

for particles in a trapped state, and

$\displaystyle\langle z^{2}\rangle_{-}=2\left(1+\frac{1}{\alpha}\right)$
$\displaystyle\langle
z^{4}\rangle_{-}=6\left(2+\frac{5}{\alpha}+\frac{4}{\alpha^{2}}\right),$ (9)

for particles in a released state. An interesting observation is that the
second moment corresponding to trapped particles does not depend on $\alpha$.
This will have interesting consequences on various physical quantities
discussed later in this work.

##  III Solution as a superposition of Gaussian distributions

In this section we develop a method that would permit us to solve Eq. (4). As
mentioned earlier, third-order differential equations are less common in
physics. In the context of active particles, a third-order differential
equation arises for RTP particles in a harmonic trap in one- and two-
dimensions Frydel22c . For higher dimensions, a differential equation of the
same system has a more complex structure and to date it is not possible to
obtain it.

To solve Eq. (4) we represent $\rho$ as a superposition of Gaussian
distributions for different effective strengths,

$\rho(z)=\int_{0}^{1}d\lambda\,p(\lambda)\rho_{G}(z;\lambda),$ (10)

where

$\rho_{G}(z;\lambda)=\sqrt{\frac{\lambda}{2\pi}}e^{-\lambda z^{2}/2}$

is a Gaussian distribution for a dimensionless strength $\lambda$. The
superposition trick will be justified later in this work where we consider a
time-dependent system. Here we assume that the superposition is correct and
proceed with a technical part.

Since $p(\lambda)$ is normalized, it follows that
$\int_{-\infty}^{\infty}dz\,\rho=1$. Also, note that $p(\lambda)$ is defined
on the interval $[0,1]$, where $\lambda=1$ is the actual strength of a trap in
reduced units, so that contributions for $\lambda>1$ have no physical meaning.

Next, we attempt to determine $p(\lambda)$. Since $p(\lambda)$ is only one
contribution of $\rho$, the other being the Gaussian distribution via the
relation in Eq. (10), it is reasonable to assume that $p(\lambda)$ is
mathematically less complex than $\rho$, and so an exact expression for
$p(\lambda)$ is plausible.

To proceed, we distinguish between the probability distributions
$p_{+}(\lambda)$ and $p_{-}(\lambda)$, such that
$\rho_{\pm}=\int_{0}^{1}d\lambda\,p_{\pm}(\lambda)\rho_{G}(z;\lambda)$. We
infer a differential equation for $p_{\pm}$ from Eq. (3). We use the word
"infer" intentionally as there is no standard technique for obtaining a
differential equation for $p_{\pm}$ from Eq. (3). One rather proceeds using a
trial and error approach. We find that the equations

$\displaystyle
0=\frac{\partial}{\partial\lambda}[\lambda(\lambda-1)p_{+}]+\frac{\alpha}{2}(p_{-}-p_{+})$
$\displaystyle
0=\frac{\partial}{\partial\lambda}[\lambda^{2}p_{-}]+\frac{\alpha}{2}(p_{+}-p_{-})$
(11)

recover Eq. (3) if we operate on each one of them with
$\int_{0}^{1}d\lambda\,\rho_{G}(z;\lambda)$. Consequently, Eq. (11) are
correct differential equations for $p_{\pm}(\lambda)$. An important fact is
that the resulting equations for $p_{\pm}$ are of a lower order than their
counterpart in Eq. (3) — as conjectured earlier.

Merging the two equations in Eq. (11) leads to

$0=\left(3\lambda^{2}-2\lambda-\alpha\lambda+\frac{\alpha}{2}\right)p+\left(\lambda-1\right)\lambda^{2}p^{\prime},$
(12)

where $p=(p_{+}+p_{-})/2$. See Appendix (C) for details of the derivation. Eq.
(12) is easily solved and the solution is

$p(\lambda)=\frac{1}{2}\frac{(\alpha/2)^{\alpha/2}e^{\alpha/2}}{\Gamma(\alpha/2)}\left(\frac{1}{\lambda}-1\right)^{\frac{\alpha}{2}-1}\frac{e^{-\frac{\alpha}{2}\frac{1}{\lambda}}}{\lambda^{3}}.$
(13)

The solutions for separate distributions $p_{\pm}$ are given by

$\displaystyle p_{+}=2\lambda p$ $\displaystyle p_{-}=2(1-\lambda)p,$ (14)

with $p$ given in Eq. (13). Both solutions can be verified by inserting into
Eq. (11).

In Fig. (1) we plot $p(\lambda)$ in Eq. (13) for different parameters
$\alpha$.

|  |
---|---|---
|  |
|  |
Figure 1:  Probability distribution $p(\lambda)$ given in Eq. (13) for
different values of $\alpha$. Distributions on the right-hand side column are
for $p_{+}$ and $p_{-}$. Blue colored distributions represent released
particles $p_{-}$ and red colored distributions represent trapped particles
$p_{+}$.

The probability distribution $p(\lambda)$ is defined on the interval
$\lambda\in[0,1]$, where $\lambda=1$ is the strength of an actual harmonic
trap in dimensionless units and $\lambda=0$ represents the absence of a
harmonic trap. The fact that $p(\lambda=0)=0$ for any $\alpha$ implies that
there is zero probability for particles to spread out into infinity. In other
words, particles always remain confined to a finite region despite there being
periods when an external potential is turned off.

The behavior at $\lambda=1$ is more interesting. The probability at
$\lambda=1$ exhibits a crossover at $\alpha=2$ such that

$p(1)=\begin{cases}0&\text{$\alpha>2$}\\\ \frac{1}{2}&\text{$\alpha=2$}\\\
\infty&\text{$\alpha<2$}.\end{cases}$

At $\lambda=1$, the probability is zero if fluctuations of a trap strength are
sufficiently fast, such that $\alpha>2$. Then for slower fluctuations, such as
$\alpha<2$, the trap is turned on for a sufficiently long time which permits
particles to attain something like a temporal equilibrium, as a result the
probability distribution at $\lambda=1$ diverges.

To visualize the tendency of distributions $p_{\pm}(\lambda)$ to develop
different properties as $\alpha\to 0$, in Fig. (2) we plot the mean $\lambda$
corresponding to each distribution, $\langle\lambda\rangle_{\pm}$. For the
limiting case $\alpha\to 0$, $\langle\lambda\rangle_{+}\to 1$ and
$\langle\lambda\rangle_{-}\to 0$, and in the limit $\alpha\to\infty$,
$\langle\lambda\rangle_{+}=\langle\lambda\rangle_{-}$.

|
---|---
Figure 2:  Average $\lambda$ for distributions $p_{+}(\lambda)$ and
$p_{-}(\lambda)$.

Despite a very different shape of each $p_{\pm}(\lambda)$, the mean value of
$\lambda$ for a total distribution $p(\lambda)$ is fixed at
$\langle\lambda\rangle=\frac{1}{2}$ for any value of $\alpha$. This reflects
the fact that on average the amount of time in which a system is in an "on"
and "off" state is the same, regardless of $\alpha$. The parameter $\alpha$
only controls the rate with which the system changes between the two states
but not the preference for being in a particular state.

Probability distributions for particles in different states $p_{\pm}(\lambda)$
are shown on the right-hand-side in Fig. (1). The two probabilities exhibit
very different properties. These differences become greater for a slowly
fluctuation harmonic trap, $\alpha\to 0$. As $\alpha$ decreases, distribution
for particles in a trapped state $p_{+}$ shifts toward $\lambda=1$, resulting
in a divergence, and distributions for particles in a released state shift
toward $\alpha=0$, without ever reaching it. In the limit $\alpha\to\infty$,
very fast fluctuation trap, the two distributions $p_{\pm}$ shift toward
$\lambda=1/2$ and converge to the same functional form.

Using the superposition formula in Eq. (10) together with the formula for
$p(\lambda)$ in Eq. (13), the solution of the third-order differential
equation in Eq. (4) is given by

$\rho(z)=\frac{(\alpha
e/2)^{\alpha/2}}{2\Gamma(\alpha/2)}\left(\frac{1}{2\pi}\right)^{\frac{1}{2}}\int_{0}^{1}d\lambda\,\left(\frac{1}{\lambda}-1\right)^{\frac{\alpha}{2}-1}\frac{e^{-\frac{\alpha}{2\lambda}}e^{-\frac{\lambda
z^{2}}{2}}}{\lambda^{3-\frac{1}{2}}}.$ (15)

The integral in Eq. (15) cannot be evaluated except for specific discrete
values of $\alpha$ corresponding to $\alpha=2n$ where $n$ is the positive
integer. For $\alpha=2$, corresponding to a crossover point, the term
$\left(\frac{1}{\lambda}-1\right)^{\frac{\alpha}{2}-1}$ in Eq. (15) becomes
unity and $p(\lambda)$ evaluates to the following analytical form

$\displaystyle\rho$ $\displaystyle=$
$\displaystyle\frac{1}{2}\frac{e^{-z^{2}/2}}{\sqrt{2\pi}}$ (16)
$\displaystyle+$
$\displaystyle\frac{1}{8\sqrt{2}}\left[\left(1-\sqrt{2z^{2}}\right)\text{erfc}\left[1+\frac{\sqrt{2z^{2}}}{2}\right]e^{1+\sqrt{2z^{2}}}\right]$
$\displaystyle+$
$\displaystyle\frac{1}{8\sqrt{2}}\left[\left(1+\sqrt{2z^{2}}\right)\text{erfc}\left[1-\frac{\sqrt{2z^{2}}}{2}\right]e^{1-\sqrt{2z^{2}}}\right].$

Note that only the first term has a Gaussian form. The remaining half of
contributions comes from non-Gaussian terms.

In Fig. (3) we plot $\rho(z)$ for different values of $\alpha$. For those
values of $\alpha$ where no analytical expression is available, the integral
in Eq. (15) is evaluated numerically. All probability distributions are in
addition compared with distributions obtained from simulation based on
numerical integration of the Langevin equation to confirm the correctness of
analytical results.

---
Figure 3:  Distributions $\rho(z)$ calculated using Eq. (15) for different
values of $\alpha$. Red data points are obtained from a simulation based on
numerical integration of the Langevin equation.

The Langevin equation is integrated using the Euler method,

$\displaystyle x_{+}(t+\Delta t)=x_{+}(t)-\mu
Kx_{+}(t)+\xi_{+}(t)\sqrt{2D\Delta t}$ $\displaystyle x_{-}(t+\Delta
t)=x_{-}(t)+\xi_{-}(t)\sqrt{2D\Delta t}$ (17)

where $\xi_{\pm}(t)$ is the white noise with zero mean and unity variance, and
$x_{+}$ and $x_{-}$ is the position of a particle in a trapped and released
state, respectively. An individual particle changes from $x_{+}$ to $x_{-}$
and vice versa at the end of the time $t_{p}$ drawn from the exponential
distribution $\sim e^{-t_{p}/\tau}$.

The primary non-Gaussian feature of the distributions is the spreading-out
effect, giving the impression that some particles are capable of leaking from
the trap. This feature becomes more pronounced for small values of $\alpha$
where particles are permitted to remain in a given state for a longer time,
allowing particles to remain in a released state for longer times, which
results in a greater spread.

We can get a better sense of what is happening by plotting distributions for
separate states $\rho_{\pm}$. In Fig. (4) we plot the distributions
$\rho_{\pm}$, calculated using the superposition formula in Eq. (10) and the
distributions $p_{\pm}$ in Eq. (14).

|
---|---
Figure 4:  Distributions of trapped and released particles, $\rho_{+}$ and
$\rho_{-}$, obtained from the superposition formula in Eq. (10) and the
distributions $p_{\pm}$ in Eq. (14) for different parameters $\alpha$.

Distributions for particles in a trapped state $\rho_{+}$ show weaker
dependence on $\alpha$. Distributions $\rho_{+}$ and $\rho_{-}$ also respond
differently to decreasing $\alpha$. While $\rho_{+}$ converges to a Gaussian
distribution for the trap strength $\lambda=1$, $\rho_{-}$ becomes more spread
out and deviates increasingly from a Gaussian form. The spreading out of
$\rho$ seen in Fig. (3) comes from particles in a released state $\rho_{-}$.

###  III.1 characteristic function

Even though the integral in Eq. (15) cannot be evaluated for an arbitrary
$\alpha$, the characteristic function
$\hat{\rho}(q)=\int_{-\infty}^{\infty}dz\,e^{-iqz}\rho(z)$ has a relatively
simple algebraic form given by

$\hat{\rho}(q)=e^{-{q^{2}}/{2}}\left[\frac{2\alpha+q^{2}}{2\alpha}\left(\frac{\alpha}{\alpha+q^{2}}\right)^{\frac{\alpha}{2}+1}\right].$
(18)

Since the factor $e^{-{q^{2}}/{2}}$ is the characteristic function of the
Boltzmann distribution, we can represent $\hat{\rho}(q)$ as

$\hat{\rho}(q)=\hat{\rho}_{eq}(q)\hat{\rho}_{neq}(q),$

where $\hat{\rho}_{eq}(q)=e^{-{q^{2}}/{2}}$, and $\hat{\rho}_{neq}(q)$
captures non-equilibrium contributions. As the product of two Fourier
transformed function corresponds to the convolution in the real space, we
could represent $\rho$ as

$\rho(z)=\int dz^{\prime}\,\rho_{eq}(z-z^{\prime})\rho_{neq}(z^{\prime}),$
(19)

and because the convolution construction implies the presence of two
independent random processes, we could interpret $\rho_{neq}$ as some type of
random process. At the crossover $\alpha=2$, $\rho_{neq}$ is found to have a
simple form given by

$\rho_{neq}=\frac{e^{-|z|\sqrt{2}}}{8}\left(3\sqrt{2}+2z\right),$ (20)

and the convolution formula in Eq. (19) smears out this result into the
formula in Eq. (16). Despite Gaussian smearing, the asymptotic behavior of
$\rho$ should be dominated by the exponential form of $\rho_{neq}$, which
indicates that $\rho$ has different asymptotic behavior than that of a
Boltzmann distribution.

##  IV Time-dependent harmonic trap

In this section we provide justification for the superposition formula in Eq.
(10). To provide such a justification, we consider a harmonic potential with a
general time-dependent strength,

$u(x,t)=\frac{K(t)x^{2}}{2}.$ (21)

It might seem that the time-dependent distribution of such a system is tricky
to calculate, but we claim that it is a Gaussian function at all times,

$\rho(x,t)=\sqrt{\frac{\mu K_{eff}}{2\pi D}}\exp\left[-\frac{\mu
K_{eff}}{D}\frac{x^{2}}{2}\right].$ (22)

If this claim is correct, then all we need to calculate to derive the
expression for $\rho(x,t)$ is the time-dependent effective strength
$K_{eff}(t)$ such that $K_{eff}(t)\neq K(t)$ (unless $K(t)$ varies very
slowly). The only constraint we introduce is that at $t=0$, $\rho(x,t)$ is a
Gaussian function corresponding to some initial effective strength
$K_{0}=K_{eff}(0)$.

To calculate $K_{eff}$, we insert the Gaussian distribution in Eq. (22) into
the corresponding time-dependent Fokker-Planck equation,

$\dot{\rho}=\mu K(t)[x\rho]^{\prime}+D\nabla^{2}\rho,$ (23)

which yields the following equation

$\mu K-\mu K_{eff}=\frac{1}{2}\frac{d\ln K_{eff}}{dt},$ (24)

for which the solution is

$K_{eff}(t)=\frac{K_{0}}{e^{-2\int_{0}^{t}dt^{\prime}\,\mu K(t^{\prime})}+2\mu
K_{0}\int_{0}^{t}dt^{\prime}\,e^{-2\int_{t^{\prime}}^{t}dt^{\prime\prime}\,\mu
K(t^{\prime\prime})}}.$ (25)

The fact that Eq. (24) can be solved implies that the Gaussian distribution in
Eq. (22) is a correct time-dependent distribution and a solution of Eq. (23).

If $K(t)$ changes in some periodic, quasi-periodic, or any other repetitive
fashion, then it should be possible to obtain an average time-independent
distribution defined as

$\rho(x)=\lim_{t\to\infty}\frac{1}{t}\int_{0}^{t}ds\,\rho(x,s).$ (26)

And since it was determined that the distribution at all times has a Gaussian
form, we could represent $\rho(x)$ defined above as a superposition of
Gaussian distributions,

$\rho(x)=\lim_{t\to\infty}\frac{1}{t}\int_{0}^{t}ds\,\rho_{G}(z,s)=\int_{0}^{1}d\lambda\,p(\lambda)\rho_{G}(z;\lambda),$
(27)

where the probability distribution $p(\lambda)$ depends on a specific
evolution of $K_{eff}(t)$. The superposition formula in Eq. (10) is simply a
reflection of this underlying behavior.

In the case that $K(t)$ changes discontinuously as $K_{0}\to K_{1}$, Eq. (25)
reduces to

$K_{eff}(t)=\frac{K_{0}K_{1}}{K_{0}-(K_{0}-K_{1})e^{-2\mu K_{1}t}}.$ (28)

In the model considered in this work, the strength of the harmonic potential
changes between two discrete values, $0$ and $K$, corresponding to the
strength of the actual trap.

Using Eq. (28), we could calculate the evolution of $K_{eff}(t)$ as $K(t)$
fluctuates between zero and the actual strength of a trap, and where the trap
persists in a given state for the length of time drawn from the exponential
distribution. The resulting evolution of $K_{eff}$ is shown in Fig. (5) for
different values of $\alpha$.

---
Figure 5:  $K_{eff}$ of a fluctuating trap for different $\alpha=1/\mu K\tau$.
The initial value is $K_{eff}(0)=K/2$ for all figures.

For $\alpha=4$, which is above the point of crossover, $K_{eff}$ fails to get
close to $K$ due to rapid fluctuations of a trap. This explains why the
probability distribution $p(\lambda)$ in Fig. (1) for the same value of
$\alpha$ is zero at $K_{eff}/K=\lambda=1$. For $\alpha=1$, which is below the
crossover, $K_{eff}$ can stay very close to $K$ for extended periods. This
gives rise to a divergence in $p(\lambda)$ at $\lambda=1$ seen in Fig. (1) for
the same value of $\alpha$.

It is possible to obtain the distribution $p(\lambda)$ from the evolution of
$K_{eff}(t)$ using
$p(\lambda)\sim\lim_{t\to\infty}\frac{1}{t}\frac{dt}{d\lambda}$, where
$\lambda(t)=K_{eff}(t)/K$.

##  V Extension to higher dimensions

It is straightforward to extend the results in Sec. (IV) to an arbitrary
dimension $d$. Given a time-dependent harmonic potential,

$u(x,t)=\frac{K(t)r^{2}}{2},$ (29)

where $r^{2}=\sum_{i=1}^{d}x_{i}^{2}$ is the radial distance from a trap
center for a system in dimension $d$. We next assume that the time-dependent
distribution has a Gaussian form at all times,

$\rho(r,t)=\left(\frac{\mu K_{eff}}{2\pi D}\right)^{d/2}\exp\left[-\frac{\mu
K_{eff}}{D}\frac{r^{2}}{2}\right],$ (30)

and to obtain an expression for $K_{eff}(t)$, we insert the Gaussian
distribution above into the corresponding Fokker-Planck equation,

$\dot{\rho}=\mu K(t)\mbox{\boldmath$\nabla$}\cdot({\bf
r}\rho)+D\nabla^{2}\rho.$ (31)

Such a procedure recovers the relation in Eq. (24) for which the solution is
the expression in Eq. (25). This proves that $K_{eff}(t)$ in Eq. (25) is valid
for any dimension.

Other results in this work can also be extended to an arbitrary dimension. The
Fokker-Planck equation for the fluctuating potential model, analogous to Eq.
(2) but for an arbitrary dimension, is given by

$\displaystyle\dot{\rho}_{+}=\mu K\mbox{\boldmath$\nabla$}\cdot[{\bf
r}\rho_{+}]+D\nabla^{2}\rho_{+}+\frac{1}{\tau}\left(\rho_{-}-\rho_{+}\right)$
$\displaystyle\dot{\rho}_{-}=D\nabla^{2}\rho_{-}-\frac{1}{\tau}\left(\rho_{-}-\rho_{+}\right).$
(32)

For a stationary state and in reduced units the two equations above become

$\displaystyle 0=\mbox{\boldmath$\nabla$}\cdot[{\bf
s}\rho_{+}]+\nabla^{2}\rho_{+}+\alpha\left(\rho_{-}-\rho_{+}\right)$
$\displaystyle 0=\nabla^{2}\rho_{-}-\alpha\left(\rho_{-}-\rho_{+}\right),$
(33)

where ${\bf s}={\bf r}\sqrt{\frac{\mu K}{D}}$.

The moments of stationary distributions can be obtained by operating on both
equations in Eq. (33) with $\int d{\bf r}\,r^{d-1}\,r^{2n}$. The resulting
expressions for the second moment are

$\displaystyle\langle s^{2}\rangle_{+}=2d$ $\displaystyle\langle
s^{2}\rangle_{-}=2d\left(1+\frac{1}{\alpha}\right)$ $\displaystyle\langle
s^{2}\rangle=d\left(2+\frac{1}{\alpha}\right).$ (34)

By combining the two equations in Eq. (33) we obtain a third-order
differential equation,

$0=-\alpha
s^{3}\rho+[3-d-(2\alpha+1-d)s^{2}]\rho^{\prime}-[3-d-s^{2}]s\rho^{\prime\prime}+s^{2}\rho^{\prime\prime\prime}.$
(35)

Note that by setting $d=1$ we recover Eq. (4). To find the solution of the
above equation, we represent a stationary distribution as a superposition of
Gaussian distributions,

$\rho(s)=\int_{0}^{1}d\lambda\,p(\lambda)\rho_{G}(s;\lambda),$ (36)

where

$\rho_{G}(s;\lambda)=\left(\frac{\lambda}{2\pi}\right)^{d/2}e^{-\lambda
s^{2}/2}.$ (37)

Since it was determined that $K_{eff}(t)$ is independent of $d$, and since it
is possible to obtain the probability distribution $p(\lambda)$ from the
evolution of $K_{eff}(t)$, we conclude that the expression for $p(\lambda)$ in
Eq. (13) is valid for all $d$. This means that any dependence on $d$ comes
from the Gaussian distribution in Eq. (37), and the solution of the third-
order differential equation in Eq. (35) is given by

$\rho(s)=\frac{(\alpha
e/2)^{\alpha/2}}{2\Gamma(\alpha/2)}\left(\frac{1}{2\pi}\right)^{\frac{d}{2}}\int_{0}^{1}d\lambda\,\left(\frac{1}{\lambda}-1\right)^{\frac{\alpha}{2}-1}\frac{e^{-\frac{\alpha}{2\lambda}}e^{-\frac{\lambda
s^{2}}{2}}}{\lambda^{3-\frac{d}{2}}}.$ (38)

The solution is not limited to integer values of $d$ and applies to any value
of $d$.

##  VI Quantities of physical interest

###  VI.1 potential energy

In this section we look into physical quantities of the model, starting with
the average potential energy $\langle u\rangle_{\pm}=K\langle
r^{2}\rangle_{\pm}/2$. Using the moments in Eq. (34), the potential energy for
particles in each state is found to be

$\displaystyle\langle u\rangle_{+}=dk_{B}T,$ $\displaystyle\langle
u\rangle_{-}=0.$ (39)

The potential energy of particles in a released state is zero since the trap
is switched off. Another interesting observation is that the potential energy
of particles in a trapped state does not depend on the rate parameter
$\alpha$. Adding the two contributions, the total potential energy is found to
be the same as that for a system in equilibrium,

$\langle u\rangle=\frac{dk_{B}T}{2}.$ (40)

The imperfect harmonic potential might considerably alter stationary
distributions but the average potential energy is unaffected.

###  VI.2 flux

We next consider flux, which for particles in a each state and for $d=1$ is
given by

$j_{-}=-D\rho^{\prime}_{-},~{}~{}~{}~{}~{}j_{+}=-\mu
Kx\rho_{+}-D\rho^{\prime}_{+},$ (41)

and is such that the total flux is zero everywhere, $j=j_{+}+j_{-}=0$.

In Fig. (6) we plot flux for particles in a released state, $j_{-}$, in
dimensionless units given by $j_{-}(z)=-\rho^{\prime}_{-}(z)$. (No insight is
gained by plotting $j_{+}$ since $j_{+}=-j_{-}$). The first observation is
that released particles move away from the trap center, which is expected. The
second observation is that the magnitude of the flux increases with increased
rate $\alpha$. This makes sense as released particles tend to be more
compressed when $\alpha$ is large and so diffusion away from the trap center
when particles are released should be greater.

|
---|---
Figure 6:  Flux for particles released from the trap. Flux for particles in a
trapped state, not shown in the figure, is $j_{+}=-j_{-}$.

As both fluxes have the same functional form, a better insight could be gained
by looking into the local velocity, in dimensionless units given by

$v_{+}=\frac{j_{+}}{\rho_{+}}=-[\ln\rho_{+}]^{\prime}-z,~{}~{}~{}~{}~{}~{}v_{-}=\frac{j_{-}}{\rho_{-}}=-[\ln\rho_{-}]^{\prime}.$
(42)

Note that the velocity in a trapped state has an additional linear term
resulting from an external force. The other term in each expression comes from
diffusion.

Local velocities are shown in Fig. (7). For particles in a trapped state we
plot $v_{+}+z$ to subtract contributions of an external force. Consequently,
the two plots amount to $-[\ln\rho_{\pm}]^{\prime}$ and represent diffusional
component of a local velocity.

|  |
---|---|---
Figure 7: Local velocity for different values of $\alpha$. For trapped
particles we plot $v_{+}+z$ to subtract contributions of an external force, as
a result, both quantities represent diffusional component of a local velocity.
In both cases this corresponds to plotting $-[\ln\rho_{\pm}]^{\prime}$.

The most striking feature is that the local velocity does not vanish at
infinity but saturate to a non-zero constant. Since the quantity that is
plotted is $-[\ln\rho_{\pm}]^{\prime}$, this implies that distributions
$\rho_{\pm}$ decay exponentially. Such an asymptotic behavior has been hinted
at in Sec. (III.1) and was proven for $\alpha=2$ in Eq. (20).

###  VI.3 entropy production rate

One feature of an active system is that it generates heat that is dissipated
into reservoir. By measuring the rate of heat that is dissipated, we could
quantify a distance of how far from equilibrium a system is. The rate of heat
dissipation, $\langle\dot{q}\rangle$, is, in addition, related to the entropy
production rate $\Pi$ via the relation $T\Pi\equiv\langle\dot{q}\rangle$
Frydel23a .

A formula for the dissipation of heat can be obtained from the drag force a
particle experiences when it moves through a dissipating medium, ${\bf
F}_{d}=-\mu^{-1}{\bf v}$, where ${\bf v}$ is the velocity vector. An
instantaneous rate of heat dissipation is then given by
$\dot{q}=\mu^{-1}v^{2}$. However, when calculating the average rate of heat
dissipation, we need to subtract equilibrium contributions,
$\langle\dot{q}\rangle=\mu^{-1}\langle v^{2}\rangle-\mu^{-1}\langle
v^{2}\rangle_{eq}$, since these contributions are not truly dissipated but are
recovered at some point in the form of thermal fluctuations as a result of the
fluctuation-dissipation relation. Since $\langle v^{2}\rangle_{eq}=dk_{B}T/m$,
we can write

$\langle\dot{q}\rangle=\frac{1}{\mu}\left[\langle
v^{2}\rangle-\frac{dD}{\tau_{r}}\right],$ (43)

where $\tau_{r}=m\mu$ is the inertial relaxation and $m$ is the mass of a
particle. For overdamped dynamics there is no inertia and $\tau_{r}=0$.

To calculate $\langle v^{2}\rangle$ we need to formulate our model in the
underdamped regime. The overdamped regime can then be recovered by setting
$\tau_{r}=0$. We find the following result for the overdamped regime:

$\langle v^{2}\rangle-\frac{dD}{\tau_{r}}=\frac{d\mu KD}{2}.$ (44)

The derivation of the above formula requires that we formulate the problem
using the Kramer’s equation, instead of the Fokker-Planck equation Frydel23a .
The details of the derivation can be found in Appendix D.

Since $\langle v^{2}\rangle$ is proportional to the kinetic energy, the term
on the right-hand-side in Eq. (44) can be seen as the excess kinetic energy
that results due to trap fluctuations. The formula, however, does not depend
on the rate of the trap fluctuations $\alpha$. This means that the kinetic
energy also does not depend on $\alpha$, similar to the potential formula
energy in Eq. (40).

Inserting Eq. (44) in Eq. (43) yields

$T\Pi\equiv\langle\dot{q}\rangle=\frac{dDK}{2}.$ (45)

The result in Eq. (45) shows that the dissipation of heat is a function of the
temperature and the trap strength and is completely independent of the rate
parameter $\alpha$ that controls fluctuations of the harmonic trap. This
independence implies that even in the limit $\alpha\to\infty$, where a
stationary distribution converges to the Gaussian distribution for the trap
strength $K/2$, the heat continues to be dissipated in the same way as for any
other value of $\alpha$ and the system is always at the same distance from
equilibrium.

##  VII Conclusion

This work considers an ideal gas confined to a harmonic trap $u=Kr^{2}/2$ with
time dependent strength $K\equiv K(t)$. It is demonstrated that the time
dependent distribution of such a system $\rho(x,t)$ has a Gaussian form at all
times, corresponding to some effective strength of the harmonic potential,
$K_{eff}(t)\neq K(t)$.

We are interested in a specific time evolution of $K(t)$ such that it changes
discontinuously between two discrete values, $0$ and $K$, and where the time
during which a trap remains in a given state is drawn from an exponential
distribution $\sim e^{-t_{p}/\tau}$.

A stationary distribution of this system can be obtained from a third-order
differential equation, the solution of which is represented as a superposition
of Gaussian distributions for different strengths of a harmonic trap,
$\rho(x)=\int_{0}^{K}dK^{\prime}\,p(K^{\prime})\rho_{G}(x;K^{\prime})$. The
reason for the superposition can be traced back to the fact that a time-
dependent distribution at all times has a Gaussian form.

The probability distribution $p(K)$ can be obtained and analyzed. The
resulting algebraic expression for $p(K)$ exhibits a crossover at $\tau=2/\mu
K$. For $\tau>2\mu K$, the distribution $p(K^{\prime})$ diverges at
$K^{\prime}=K$, indicating that a fraction of particles comes close to an
equilibrium-like behavior. The divergence disappears for $\tau<2/\mu K$, and
for very small $\tau$, a stationary distribution converges to a Gaussian form
for a trap with the strength which is half the strength of the physical trap,
$K/2$.

The primary feature of the resulting stationary distribution
$\rho(x)=\int_{0}^{K}dK^{\prime}\,p(K^{\prime})\rho_{G}(x;K^{\prime})$ is the
spreading-out effect, which is understood to be the result of a trap being
turned off. Since the period during which a trap is turned off increases with
$\tau$, the distributions tend to be more spread out for larger $\tau$. As a
result of spreading-out, the asymptotic behavior of $\rho(x)$ is dominated by
an exponential decay.

The spreading-out feature of the current model is different from typical
active particles that tend, rather than to spread-out and penetrate the trap
boundaries, to accumulate at the trap boundaries. This behavior is seen for
RTP and ABP particles. The accumulation of particles at the boundaries leads
to bimodal stationary distributions, with excess of particles close to
boundaries and depletion of particles near a trap center. These models exhibit
crossover when a stationary distribution changes from a bimodal to unimodal
distribution.

###### Acknowledgements.

D.F. acknowledges financial support from FONDECYT through grant number
1241694.

##  VIII DATA AVAILABILITY

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

##  Appendix A Derivation of Eq. (4)

In the section we provide details how the two coupled equations in Eq. (3),
which we reproduce below,

$\displaystyle
0=[z\rho_{+}]^{\prime}+\rho_{+}^{\prime\prime}+\alpha\left(\rho_{-}-\rho_{+}\right),$
$\displaystyle
0=\rho_{-}^{\prime\prime}+\alpha\left(\rho_{+}-\rho_{-}\right),$ (46)

are combined to yield the third-order differential equation in Eq. (4). Adding
and subtracting the two equations yield

$\displaystyle
0=\frac{1}{2}[z\rho]^{\prime}+\frac{1}{2}[z\sigma]^{\prime}+\rho^{\prime\prime},$
$\displaystyle
0=\frac{1}{2}[z\rho]^{\prime}+\frac{1}{2}[z\sigma]^{\prime}+\sigma^{\prime\prime}-2\alpha\sigma,$
(47)

where $\rho=(\rho_{+}+\rho_{-})/2$ and $\sigma=(\rho_{+}-\rho_{-})/2$. From
the first equation we get $z\sigma=-z\rho-2\rho^{\prime}+\text{const},$ where
$\text{const}=0$ as a result of the even symmetry of $\rho(z)$. From this we
can generate a sequence of expressions

$\displaystyle\sigma=-\rho-\frac{2}{z}\rho^{\prime},$
$\displaystyle\sigma^{\prime}=-\rho^{\prime}+\frac{2}{z^{2}}\rho^{\prime}-\frac{2}{z}\rho^{\prime\prime},$
$\displaystyle\sigma^{\prime\prime}=-\rho^{\prime\prime}-\frac{4}{z^{3}}\rho^{\prime}+\frac{4}{z^{2}}\rho^{\prime\prime}-\frac{2}{z}\rho^{\prime\prime\prime}.$
(48)

Inserting these expressions into the second equation in Eq. (47) recovers Eq.
(4).

##  Appendix B Recurrence relation

This section shows how even moments, given by

$\langle z^{2n}\rangle_{\pm}=\int_{-\infty}^{\infty}dz\,z^{2n}\rho_{\pm}(z)$

can be obtained from the two equations in Eq. (3), which we repeat below for
clarity:

$\displaystyle
0=[z\rho_{+}]^{\prime}+\rho_{+}^{\prime\prime}+\alpha\left(\rho_{-}-\rho_{+}\right)$
$\displaystyle
0=\rho_{-}^{\prime\prime}+\alpha\left(\rho_{+}-\rho_{-}\right).$ (49)

To convert the two equations into the recurrence relation, we operate on both
equation with $\int_{-\infty}^{\infty}dz\,z^{2n+2}$. This results in the
following two coupled recurrence relations:

$\displaystyle
0=-(2n+2)m_{n+1}^{+}+(2n+2)(2n+1)m_{n}^{+}+\alpha\left(m_{n+1}^{-}-m_{n+1}^{+}\right)$
$\displaystyle
0=(2n+2)(2n+1)m_{n}^{-}+\alpha\left(m_{n+1}^{+}-m_{n+1}^{-}\right),$ (50)

where $m_{n}^{+}=\langle z^{2n}\rangle_{+}$ and $m_{n}^{-}=\langle
z^{2n}\rangle_{-}$. After rearrangement, the two equations become

$\displaystyle m_{n+1}^{+}=(2n+1)(m_{n}^{+}+m_{n}^{-})$ $\displaystyle
m_{n+1}^{+}=(2n+1)\left(m_{n}^{+}-m_{n}^{-}\right)-\frac{2\alpha}{(2n+2)}\left(m_{n+1}^{+}-m_{n+1}^{-}\right).$

The initial terms of both sequences are $m_{0}^{+}=m_{0}^{-}=1$ and the
remaining terms are obtained from the recurrence relations above.

##  Appendix C Derivation of Eq. (11)

In this section we demonstrate how combining the two equations in Eq. (11),
which we reproduce below,

$\displaystyle
0=[\lambda(\lambda-1)p_{+}]^{\prime}+\frac{\alpha}{2}(p_{-}-p_{+})$
$\displaystyle 0=[\lambda^{2}p_{-}]^{\prime}+\frac{\alpha}{2}(p_{+}-p_{-}),$
(51)

leads to the first order differential equation in Eq. (14). Adding and
subtracting the two equations yield

$\displaystyle 0=[\lambda(2\lambda-1)p]^{\prime}-\lambda\sigma^{\prime}$
$\displaystyle 0=2[\lambda^{2}\sigma]^{\prime}-\lambda
p^{\prime}-\lambda\sigma^{\prime}-2\alpha\sigma,$ (52)

where $p=(p_{+}+p_{-})/2$ and $\sigma=(p_{+}-p_{-})/2$. From the first
equation we get $\sigma=(2\lambda-1)p$ and from which we generate the
consecutive expressions

$\displaystyle\sigma=(2\lambda-1)p,$
$\displaystyle\sigma^{\prime}=2p+(2\lambda-1)p^{\prime},$ (53)

Inserting these expressions into the second equation in Eq. (52) recovers Eq.
(14).

##  Appendix D  Derivation of $\langle v^{2}\rangle$

To obtain the expression for $\langle v^{2}\rangle$, we need to formulate our
model within underdamped dynamics. This means that we need to formulate the
system within Kramer’s equation, which is the type of Fokker-Planck equation
for the underdamped regime Frydel23a . The evolution of the distributions
$\rho_{\pm}\equiv\rho_{\pm}({\bf r},{\bf v},t)$ is governed by the following
Kramer’s equations:

$\displaystyle\frac{\partial\rho_{+}}{\partial t}$ $\displaystyle=$
$\displaystyle-{\bf
v}\cdot\mbox{\boldmath$\nabla$}_{r}\rho_{+}+\frac{1}{\tau_{r}}\mbox{\boldmath$\nabla$}_{v}\cdot\left[\left(\mu
K{\bf r}+{\bf
v}\right)\rho_{+}\right]+\frac{D}{\tau_{r}^{2}}\nabla_{v}^{2}\rho_{+}$
$\displaystyle+$ $\displaystyle\frac{1}{\tau}(\rho_{-}-\rho_{+})$
$\displaystyle\frac{\partial\rho_{-}}{\partial t}$ $\displaystyle=$
$\displaystyle-{\bf
v}\cdot\mbox{\boldmath$\nabla$}_{r}\rho_{-}+\frac{1}{\tau_{r}}\mbox{\boldmath$\nabla$}_{v}\cdot\left[{\bf
v}\rho_{-}\right]+\frac{D}{\tau_{r}^{2}}\nabla_{v}^{2}\rho_{-}+\frac{1}{\tau}(\rho_{+}-\rho_{-}),$

and we recall that $\tau_{r}=m\mu$ is the inertial relaxation, $m$ is a
particle mass, and the overdamped regime is recovered for $\tau_{r}=0$.

To calculate $\langle v^{2}\rangle$, or other average quantities of interest,
we operate on the stationary Kramer’s equation, $\dot{\rho}=0$, with an
integral operator $\hat{O}_{g}=\int d{\bf v}\int d{\bf r}\,\,g({\bf r},{\bf
v}),$ where the function $g$ is going to be defined later. Using integration
by parts, the terms of a transformed Kramer’s equation can be represented as
average quantities, designated by the angular brackets
$\langle\dots\rangle=\int d{\bf v}\int d{\bf r}\,\,(\dots)$. The resulting two
equations are

$\displaystyle 0$ $\displaystyle=$ $\displaystyle\langle{\bf
v}\cdot\mbox{\boldmath$\nabla$}_{r}g\rangle_{+}-\frac{1}{\tau_{r}}\langle\left(\mu
K{\bf r}+{\bf
v}\right)\cdot\mbox{\boldmath$\nabla$}_{v}g\rangle_{+}+\frac{D}{\tau_{r}^{2}}\langle\nabla_{v}^{2}g\rangle_{+}$
$\displaystyle+$ $\displaystyle\frac{\langle g\rangle_{-}}{\tau}-\frac{\langle
g\rangle_{+}}{\tau}$ $\displaystyle 0=\langle{\bf
v}\cdot\mbox{\boldmath$\nabla$}_{r}g\rangle_{-}-\frac{1}{\tau_{r}}\langle{\bf
v}\cdot\mbox{\boldmath$\nabla$}_{v}g\rangle_{-}+\frac{D}{\tau_{r}^{2}}\langle\nabla_{v}^{2}g\rangle_{-}+\frac{\langle
g\rangle_{+}}{\tau}-\frac{\langle g\rangle_{-}}{\tau}.$

Using $g=v^{2}$, $g=r^{2}$, and $g={\bf r}\cdot{\bf v}$, we can generate from
the relation above six equations involving six unknown average quantities,

$\displaystyle 0=-\frac{2}{\tau_{r}}\langle
v^{2}\rangle_{+}+\frac{2dD}{\tau_{r}^{2}}-\frac{1}{\tau}\langle
v^{2}\rangle_{+}+\frac{1}{\tau}\langle v^{2}\rangle_{-}-\frac{2\mu
K}{\tau_{r}}\langle{\bf r}\cdot{\bf v}\rangle_{+}$ $\displaystyle
0=-\frac{2}{\tau_{r}}\langle
v^{2}\rangle_{-}+\frac{2dD}{\tau_{r}^{2}}-\frac{1}{\tau}\langle
v^{2}\rangle_{-}+\frac{1}{\tau}\langle v^{2}\rangle_{+}$ $\displaystyle
0=\langle v^{2}\rangle_{+}+\frac{1}{\tau_{r}}\langle{\bf v}\cdot{\bf
r}\rangle_{+}-\frac{1}{\tau}\langle{\bf v}\cdot{\bf
r}\rangle_{+}+\frac{1}{\tau}\langle{\bf v}\cdot{\bf r}\rangle_{-}-\frac{\mu
K}{\tau_{r}}\langle r^{2}\rangle_{+}$ $\displaystyle 0=\langle
v^{2}\rangle_{-}-\frac{1}{\tau_{r}}\langle{\bf v}\cdot{}\bf
r\rangle_{-}-\frac{1}{\tau}\langle{\bf v}\cdot{\bf
r}\rangle_{-}+\frac{1}{\tau}\langle{\bf v}\cdot{\bf r}\rangle_{+}$
$\displaystyle 0=2\langle{\bf v}\cdot{\bf r}\rangle_{+}-\frac{1}{\tau}\langle
r^{2}\rangle_{+}+\frac{1}{\tau}\langle r^{2}\rangle_{-}$ $\displaystyle
0=2\langle{\bf v}\cdot{\bf r}\rangle_{-}-\frac{1}{\tau}\langle
r^{2}\rangle_{-}+\frac{1}{\tau}\langle r^{2}\rangle_{+}.$ (54)

After solving the system of coupled equations above, we get

$\langle v^{2}\rangle-\frac{dD}{\tau_{r}}=-\frac{d\tau D\mu
K(\tau_{r}+\tau)}{\mu
K\tau_{r}^{2}\tau-4\tau_{r}^{2}-6\tau_{r}\tau-2\tau^{2}},$ (55)

which for the overdamped regime $\tau_{r}=0$ reduces to

$\langle v^{2}\rangle-\frac{dD}{\tau_{r}}=\frac{dD\mu K}{2}.$ (56)

Since $\langle v^{2}\rangle_{eq}=\frac{dD}{\tau_{r}}$, and since $\langle
v^{2}\rangle$ is proportional to the kinetic energy, we can think of the above
quantity as representing the excess kinetic energy that arises due to
fluctuations of a harmonic trap.

## References

  * (1) J. Tailleur and M. E. Cates, Statistical Mechanics of Interacting Run-and-Tumble Bacteria, Phys. Rev. Lett. 100, 218103 (2008).
  * (2) J. Tailleur and M. E. Cates, Sedimentation, trapping, and rectification of dilute bacteria, Europhys. Lett. 86, 60002 (2009).
  * (3) M. R. Evans and S. N. Majumdar, Run and tumble particle under resetting: a renewal approach, J. Phys. A: Math. Theor. 51 475003 (2018).
  * (4) K. Malakar, V. Jemseena, A. Kundu, K. V. Kumar, S. Sabhapandit, S. N. Majumdar, S. Redner, and A. Dhar, Steady state, relaxation and first-passage properties of a run-and-tumble particle in one-dimension, J. Stat. Mech.: Theory Exp. 043215 (2018).
  * (5) A. Dhar, A. Kundu, S. N. Majumdar, S. Sabhapandit, and G. Schehr, Run-and-tumble particle in one-dimensional confining potentials: Steady-state, relaxation, and first-passage properties, Phys. Rev. E 99, 032132 (2019).
  * (6) U. Basu, S. N. Majumdar, A. Rosso, S. Sabhapandit, and G. Schehr, Exact stationary state of a run-and-tumble particle with three internal states in a harmonic trap, J. Phys. A 53, 09LT01 (2020).
  * (7) A. Villa-Torrealba, C.l Chávez-Raby, P. de Castro, and R. Soto, Run-and-tumble bacteria slowly approaching the diffusive regime, Phys. Rev. E 101, 062607 (2020),
  * (8) N. Razin, Entropy production of an active particle in a box, Phys. Rev. E (R) 102, 030103(R) (2020).
  * (9) D. Frydel, Generalized run-and-tumble model for an arbitrary distribution of velocities in 1D geometry, J. Stat. Mech. 103, 052603 (2021).
  * (10) D. Frydel, Kuramoto model with run-and-tumble dynamics, Phys. Rev. E 104, 024203 (2021).
  * (11) N. R. Smith and O. Farago, Nonequilibrium steady state for harmonically confined active particles, Phys. Rev. E 106, 054118 (2022).
  * (12) D. Frydel, The four-state RTP model: exact solution at zero temperature, Phys. Fluids 34, 027111 (2022).
  * (13) C. Roberts and Z. Zhen, Run-and-tumble motion in a linear ratchet potential: Analytic solution, power extraction, and first-passage properties, Phys. Rev. E 108, 014139 (2023).
  * (14) O. Farago and N. R. Smith, Confined run-and-tumble particles with non-Markovian tumbling statistics, Phys. Rev. E 109, 044121 (2024),
  * (15) T. Pietrangeli, C. Ybert, C. Cottin-Bizonne, and F. Detcheverry, Optimal run-and-tumble in slit-like confinement, Phys. Rev. Research 6, 023028 (2024).
  * (16) B. Loewe, T. Kozhukhov and, T. N. Shendruk Anisotropic run-and-tumble-turn dynamics, Soft Matter 20, 1133 (2024).
  * (17) K. Malakar, A. Das, A. Kundu, K. V. Kumar, and A. Dhar Steady state of an active Brownian particle in a two-dimensional harmonic trap, Phys. Rev. E 101, 022610 (2020).
  * (18) M. Caraglio, T. Franosch, Analytic solution of an active brownian particle in a harmonic well, Phys. Rev. Lett. 129, 158001, (2022).
  * (19) G. Szamel, Self-propelled particle in an external potential: Existence of an effective temperature, Phys, Rev. E 90, 012111 (2014).
  * (20) S. Shankar and M. C. Marchetti, Hidden entropy production and work fluctuations in an ideal active gas, Phys. Rev. E 98, 020604(R) (2018).
  * (21) F. J. Sevilla, R. F. Rodríguez, and J. R. Gomez-Solano, Generalized Ornstein-Uhlenbeck model for active motion, Phys. Rev. E 100, 032123 (2019).
  * (22) D. Martin, J. O’Byrne, M. E. Cates, É. Fodor, C. Nardini, J. Tailleur, and F. van Wijland Statistical mechanics of active Ornstein-Uhlenbeck particles, Phys. Rev. E 103, 032607 (2021).
  * (23) L. Caprini, U. M. B. Marconi, and A. Vulpiani, Linear response and correlation of a self-propelled particle in the presence of external fields, J. Stat. Mech.: Theory Exp. 033203 (2018).
  * (24) G H Philipp Nguyen, Active Ornstein–Uhlenbeck model for self-propelled particles with inertia, J. Phys.: Condens. Matter 34, 035101 (2022).
  * (25) Y. Sokolov, D. Frydel, D. G. Grier, H. Diamant, Y. Roichman, Hydrodynamic pair attractions between driven colloidal particles, Phys. Rev. Lett. 107, 158302 (2011).
  * (26) S. C. Takatori, R. De Dier, J. Vermant, and J. F. Brady, Acoustic trapping of active matter, Nat. Commun. 7, 10694 (2016).
  * (27) R. Goerlich, M. Li, S. Albert, G. Manfredi, P. Hervieux, and C. Genet, Noise and ergodic properties of Brownian motion in an optical tweezer: Looking at regime crossovers in an Ornstein-Uhlenbeck process, Phys. Rev. E 103, 032132 (2021).
  * (28) I. Buttinoni, L. Caprini, L. Alvarez , F. J. Schwarzendahl, and H. Löwen, Active colloids in harmonic optical potentials, EPL 140, 27001 (2022).
  * (29) R. Goerlich, L.B. Pires, G. Manfredi, P.A. Hervieux, C. Genet, Harvesting information to control nonequilibrium states of active matter, Phys. Rev. E 106, 054617, (2022).
  * (30) O. Chor, A. Sohachi, R. Goerlich, E. Rosen, S. Rahav, and Y. Roichman Many-body Szilárd engine with giant number fluctuations, Phys. Rev. Research 5, 043193 (2023).
  * (31) D. Frydel, Positing the problem of stationary distributions of active particles as third-order differential equation, Phys. Rev. E 106, 024121 (2022).
  * (32) D. Frydel, Run-and-tumble oscillator: Moment analysis of stationary distributions, Phys. Fluids 35, 101905 (2023).
  * (33) D. Frydel, Active oscillator: Recurrence relation approach , Phys. Fluids 36, 011910 (2024).
  * (34) H.Kleinert, Path Integrals in Quantum Mechanics, Statistics, Polymer Physics and Financial Markets (World Scientific, Singapore, 2004).
  * (35) Petr Jizba and Hagen Kleinert, Superpositions of probability distributions, Phys. Rev. E 78, 031122 (2008).
  * (36) C. Beck, Dynamical Foundations of Nonextensive Statistical Mechanics, Phys. Rev. Lett. 87, 180601 (2001).
  * (37) J. M. Parrondo, J. M. Horowitz, and T. Sagawa, Thermodynamics of information, Nat. Phys. 11, 131 (2015).
  * (38) T. Sagawa and M. Ueda, Second law of thermodynamics with discrete quantum feedback control, Phys. Rev. Lett. 100, 080403 (2008).
  * (39) F. J. Cao and M. Feito, Thermodynamics of feedback controlled systems, Phys. Rev. E 79, 041118 (2009).
  * (40) J. M. Horowitz and S. Vaikuntanathan, Nonequilibrium detailed fluctuation theorem for repeated discrete feedback, Phys. Rev. E 82, 061120 (2010).
  * (41) T. Sagawa and M. Ueda, Generalized Jarzynski equality under nonequilibrium feedback control, Phys. Rev. Lett. 104, 090602 (2010).
  * (42) M. Ponmurugan, Generalized detailed fluctuation theorem under nonequilibrium feedback control, Phys. Rev. E 82, 031129 (2010).
  * (43) Y. Fujitani and H. Suzuki, Jarzynski equality modified in the linear feedback system, J. Phys. Soc. Jpn. 79, 104003 (2010).
  * (44) A. Kundu, Nonequilibrium fluctuation theorem for systems under discrete and continuous feedback control, Phys. Rev. E 86, 021107 (2012).
  * (45) Seshadev Padhi and Smita Pati, Theory of Third-Order Differential Equations, (Springer Science and Business Media, 2013).
  * (46) D. Frydel, Entropy production of active particles formulated for underdamped dynamics, Phys. Rev. E 107, 014604 (2023).
  * (47) D. Frydel, Intuitive view of entropy production of ideal run-and-tumble particles, Phys. Rev. E 105, 034113 (2022).


# Interpreting the Residual Stream of ResNet18

André Longon
Division of Computer Science and Engineering
Louisiana State University
Baton Rouge, LA 70803
<EMAIL_ADDRESS>

###### Abstract

A mechanistic understanding of the computations learned by deep neural
networks (DNNs) is far from complete. In the domain of visual object
recognition, prior research has illuminated inner workings of InceptionV1, but
DNNs with different architectures have remained largely unexplored. This work
investigates ResNet18 with a particular focus on its residual stream, an
architectural mechanism which InceptionV1 lacks. We observe that for a given
block, channel features of the stream are updated along a spectrum: either the
input feature skips to the output, the block feature overwrites the output, or
the output is some mixture between the input and block features. Furthermore,
we show that many residual stream channels compute scale invariant
representations through a mixture of the input’s smaller-scale feature with
the block’s larger-scale feature. This not only mounts evidence for the
universality of scale equivariance, but also presents how the residual stream
further implements scale invariance. Collectively, our results begin an
interpretation of the residual stream in visual object recognition, finding it
to be a flexible feature manager and a medium to build scale invariant
representations.

##  1 Introduction

Deep neural networks (DNNs) lead humanity into a quandary: they learn to
perform increasingly sophisticated tasks but leave us ignorant of how they
ultimately solve them. As DNNs are deployed in safety-critical applications
such as autonomous driving and medical diagnostics, understanding their inner
workings is imperative to better predict their strengths and limitations.
Additionally, the leading models of neural response predictions in visual [13]
and language [14] processing are DNNs, so reverse-engineering these models may
uncover the neural mechanisms from which these cognitive abilities emerge.

Toward this end, the field of mechanistic interpretability has progressed in
illuminating DNN circuitry with a concentration on vision [10] and language
[3] models. In the visual domain, InceptionV1 [16] is the most well-studied
DNN where curve detectors [1], equivariances for hue, rotation, and scale
[11], and high-low frequency detectors [15] have been identified. However, a
key issue of mechanistic interpretability is universality: to what extent do
these mechanisms exist in other DNNs, such as those with different
architectures or trained on different datasets?

To expand the mechanistic interpretation of visual models, we zoom into the
residual stream of ResNet18 [5], offering an empirical account of its
behavior. The residual stream, which InceptionV1 lacks, is an architectural
mechanism that permits features to bypass layers of processing via their
summation to the output of a downstream layer. With ResNet18 being the
smallest member of the ResNet family, it is manageable to perform a network-
scale study of its stream, striking a balance between a simplistic “toy model”
while having $\sim$5x less parameters than ResNet152, the largest of the
family.

We present the following results from our investigation:

  * •

Stream channels at a given block update their features either by skipping the
input to the output of the block, erasing the input and overwriting with the
block feature if present, or mixing the input and block features (see Figure
1A). We highlight network trends of how channels are distributed across this
skip-overwrite spectrum and demonstrate some mechanisms of how the spectrum is
implemented.

  * •

With simple criteria, residual stream channels which exhibit scale invariance
are extracted (see Figure 1B). The input channel possesses the smaller scale
version, the block’s pre-sum output possesses the larger scale, and these are
summed together to produce the final scale invariant output. This pattern is
especially pronounced in intermediate blocks. These findings not only mount
evidence for the universality of scale equivariance, but further show how the
residual stream uses these features to construct scale invariance.

Figure 1: Each grid contains feature visualizations of the channel’s center
neuron, the entire channel, and the top 9 center neuron activating natural
images from the ImageNet validation set [2]. All channel indices are zero-
indexed. A: Exemplary overwrite (38), mix (19), and skip (25) channels from
output 1.1 (1.0 input + 1.1.bn2 block). Note the presence or lack of overlap
in top vals between Output and Input and/or Block. B: Exemplary scale
invariant channels (178, 215, 25) from output 3.1.

Together, our results serve as a starting point for exploring the function of
the residual stream in visual processing. We hope to inspire further work on
elucidating the circuitry of these mechanisms, and even to search for similar
computations in biological visual systems. It is known that the macaque [6]
and mouse [4] visual systems contain so-called “bypass connections”, where
neurons project to farther downstream processing areas, skipping over
intermediate areas. This anatomical feature is highly analogous to the DNN
residual stream, so it is conjectured that a function of bypass connections
is, for instance, to compute scale invariance. This possibility is important
for visual neuroscience and to further establish universality in computations
between artificial and biological neural networks.

##  2 Methods

We use the PyTorch [12] torchvision implementation of ResNet18 with the
default ImageNet-trained [2] weights, `IMAGENET1K_V1`. Neuron activations are
obtained via the THINGSvision library [8].

###  2.1 Residual Stream Architecture

After an initial convolutional layer and max pool, the ResNet18 architecture
is an alternating block structure of two types: one with a downsample
submodule and one without, which we refer to as a downsample block and simple
block respectively. In a simple block, its input tensor is element-wise added
with the tensor obtained from its last processing step, the second batch
normalization (BN). Note there is no ReLU nonlinearity before this addition.
In a downsample block, this element-wise addition cannot occur as the block’s
first convolution doubles the number of channels and halves both spatial
dimensions via a stride of 2. Therefore, the input is passed through a
separate convolutional layer with a stride of 2 and doubled channels, followed
by a BN. The result of this downsample submodule is then added with the second
BN tensor. In both downsample and simple blocks, the summed result undergoes a
ReLU nonlinearity to produce the final output activations of the block.

#### Notation

This work is only concerned with the activations from three layers for a given
block. We name the layers $B$, $I$, and $O$ for block, input, and output
respectively (the block index is not notated as we analyze them in isolation).
$B$ is the block’s second BN. $I$ is the layer whose output bypasses the main
layers of the block and added to the output of $B$. For the simple block, $I$
is the previous block’s output. In the downsample block, $I$ is the result
from the BN in the downsample submodule. Finally, $O$ is the layer which
outputs the sum of outputs from $I$ and $B$ followed by a ReLU.

As an example, if the given block is simple such as 3.1, then $B$ is the
second BN in 3.1, $I$ is the output of block 3.0, and $O$ is the output of 3.1
itself. For a downsample block such as 4.0, $I$ is the BN of the downsample
submodule in 4.0, $B$ is still the second BN from 4.0, and $O$ is still the
final output of 4.0.

###  2.2 Feature Visualization

Our results rely on feature visualization, an optimization algorithm used to
generate images which maximally activate a target channel (or an individual
neuron) in a given layer. Such images are henceforth referred to as feature
visualizations (FZs). We use the version of feature visualization described by
[9] and implemented in the Lucent
library111https://github.com/greentfrapp/lucent for PyTorch. We use near-
identical regularizations as prescribed by [9] for improved interpretability,
with altered jitter values discussed in the appendix (see Section A.2).

We use a channel’s center-neuron optimized FZs for all subsequent analyses.
While we also display channel-optimized FZs, they are for visualization
purposes only. We also do not study the first residual stream block (1.0) in
this study. This is due to its input resulting from a maxpool layer, making it
is difficult to obtain FZs due to frequent zero-valued gradients222This
problem also arises when optimizing on a post-ReLU neuron. To circumvent this,
Identity layers are added after stream summation to obtain pre-ReLU
activations.. Therefore, we start at the 1.1 block and work our way to the end
of the stream through 4.1.

#### Notation

In a given block of the residual stream, for a given channel $c$, we obtain
the center-neuron FZs for $I_{c},B_{c},\textrm{and }O_{c}$, which we denote as
$\hat{X}_{I_{c}},\hat{X}_{B_{c}},\textrm{and }\hat{X}_{O_{c}}$ respectively.
We denote the output of layer $I$’s center-neuron activation in channel $c$ to
image $X$ as $I_{c}(X)$ (equivalent notation for $B$ and $O$). Thus,
$I_{c}(\hat{X}_{I_{c}})$ is the center-neuron activation in channel $c$
resulting from its optimized FZ.

###  2.3 Mix Ratio

In order to measure where a given block’s channel falls along the Skip-
Overwrite spectrum, we define the mix ratio as:

$M_{c}=\frac{O_{c}(\hat{X}_{I_{c}})}{O_{c}(\hat{X}_{B_{c}})}$ (1)

$M_{c}$ measures the relative activations of the two components that when
summed, make up $O_{c}$: $I_{c}$ and $B_{c}$. The mix ratio informs us of the
importance of the input feature versus the block feature. A high $M$ thus
indicates a skip-like behavior, a low $M$ is overwrite-like, and an $M$ close
to $1$ is a mixture of input and block features (resembling an OR gate). We
max-clamp $M$ to $5$ for our analyses, and as a channel’s $B_{c}$ can be
negative, we also set those $M_{c}$ values to $5$ to indicate a skip.

###  2.4 Scale Invariance Criteria

We assume that if a particular $O_{c}$ at a given block is to be scale
invariant, then $I_{c}$ will represent the smaller-scale copy, and $B_{c}$
will represent the larger-scale copy. This is a reasonable assumption as the
receptive field of a neuron in $I_{c}$ is always smaller than one in $B_{c}$
(or equal in deep enough blocks). We then determine if $O_{c}$ exhibits scale
invariance by checking if all the following conditions are true:

$\frac{2}{3}<M_{c}<\frac{3}{2}$ (2)
$\textrm{ReLU}(B_{c}(\hat{X}_{I_{c}}))<B_{c}(S(\hat{X}_{I_{c}}))$ (3)
$\textrm{ReLU}(I_{c}(\hat{X}_{B_{c}}))<I_{c}(S^{-1}(\hat{X}_{B_{c}}))$ (4)

where

$S(X):=\textrm{Resize}(\textrm{CenterCrop(X)})\quad\textrm{and}\quad
S^{-1}(X):=\textrm{Pad}(\textrm{Resize(X)}).$

$S(X)$ denotes a scale transformation of an image where it is first center-
cropped to 112 pixels (half the spatial resolution), then resized back to 224
pixels using bilinear interpolation. $S^{-1}(X)$ denotes an “inverse”333We
acknowledge that $S^{-1}(X)$ is not a true inverse so there is an asymmetry,
but we find it sufficient to provide an entry point into this study. scale
transformation where the image is resized to 112, then reflection padded to
224.

The rationale for these criteria is as follows:

  * •

Equation 2 ensures the block output $O_{c}$ is invariant (with some tolerance)
to the scale of the input and block features, where both the smaller-scale
$\hat{X}_{I_{c}}$ and larger-scale $\hat{X}_{B_{c}}$ activate $O_{c}$
similarly.

  * •

Equation 3 increases the scale of the smaller-scale neuron’s FZ. It then
checks if the upscaled FZ activates the larger-scale neuron more so than the
unscaled FZ. This ensures that scaling up $\hat{X}_{I_{c}}$ makes it more
closely resemble the larger-scale neuron’s FZ.

  * •

Equation 4 similarly decreases the scale of the larger-scale neuron’s FZ to
see if it increases the activation of the smaller-scale neuron relative to the
unscaled FZ.

We apply a ReLU in the unscaled case to guarantee we are not bringing a
negative activation from the unscaled FZ up to 0 activation (or less negative)
in the scaled case, indicative of destroying inhibitory selectivity without
necessarily increasing the selectivity for the excitatory feature.

#### Scale Metric

When channels that meet our criteria are identified, we next quantify the
degree of their scale invariance by:

$SM_{c}:=\frac{B_{c}(S(\hat{X}_{I_{c}}))-\textrm{ReLU}(B_{c}(\hat{X}_{I_{c}}))}{B_{c}(\hat{X}_{B_{c}})}+\frac{I_{c}(S^{-1}(\hat{X}_{B_{c}}))-\textrm{ReLU}(I_{c}(\hat{X}_{B_{c}}))}{I_{c}(\hat{X}_{I_{c}})}$
(5)

which is the sum of activation deltas $S(X)$ and $S^{-1}(X)$ produce from
equations 3 and 4, scaled by input and block’s respective FZ-evoked
activations to account for different activation ranges across the channels. We
perform this scaling as this metric will be later used to rank the channels.

##  3 Results

###  3.1 Skip-Overwrite Spectrum

To begin, we obtain mix ratios $M_{c}$ for all channels in every block
studied: 1.1 through 4.1. As seen in Figure 2, blocks exhibit a variety of
distributions across the spectrum. All downsample blocks are dominated by
overwrite channels, whereas simple blocks exhibit less skewed distributions.
We note that the final block, 4.1, is an exception to this pattern: it is a
simple block with all of its channels in the lowest bin (we confirm that
nearly all of them are non-zero).

Figure 2: Histograms of mix ratios for all inspected blocks of ResNet18,
starting at 1.1 through 4.1. Bins are from 0 to 5 at 0.25 intervals. Close to
1 is a mixture of input and block features, below 1 indicates an overwrite
behavior, above 1 indicates a skip behavior.

#### Skip Weight Analysis

To explore the mechanisms that implement this spectrum, we first suppose that
the weights of skip channels will have lower magnitude weights relative to the
layer. In Figure 3, we plot each channel’s second BN weight magnitude and the
mean weight magnitude from their second convolutional layer. We do not display
the downsample blocks and 4.1 due to the sharp overwrite skew. In the
remaining simple blocks, we confirm that higher mix ratio channels possess
significantly smaller weight magnitudes in the inspected layers. It seems
ResNet18 simply silences the block output with lower weights, allowing the
input to skip to the output relatively undisturbed.

Figure 3: Weight magnitudes with respect to mix ratio of second batch
normalization (top row) and the second convolutional layer (bottom row, mean
taken). Spearman rank correlations and p-values displayed for each block,
which show that channels with more skip-like behavior have smaller weight
magnitudes.

#### Overwrite Inhibition

We observe that for many channels, $B_{c}(\hat{X}_{I_{c}})<0$, which appears
instrumental to facilitate overwrite. That is, when presented with the input’s
preferred feature, the block’s second BN inhibits the output of the channel.
In Figure 4 (Top Row), we see an abundance of such negative activations
amongst overwrite channels ($M_{c}<1$), significantly lower than the
activations of skip channels ($M_{c}>1$) in blocks 2.1 and 3.1. However, we
find that many skip channels also possess such negative activations, so it is
difficult to make conclusions from this result alone.

Next, we obtained the activations from the input and block to all ImageNet
validation images where the input responded with a positive activation. We
then took a Spearman rank correlation of the block activations with respect to
the input activations for each channel and report the results in 4 (Bottom
Row), plotted against the mix ratio. If input features inhibit the block for
overwrite channels, we expect them to have significantly more negative
correlations (as input activation increases, the block tends to decrease)
compared to skip channels. This is the case for block 3.1, but not 1.1 nor
2.1.

Figure 4: Top Row: Activations of second batch norm (pre-sum block output)
with respect to mix ratio, with activation means for skip (mix > 1, blue) and
overwrite (mix < 1, red) channels. Bottom Row: Spearman rank correlations of
Block activations with respect to Input activations for ImageNet validation
images which positively activate Input, plotted against mix ratio. ANOVA Fs
and p-values between skip and overwrite activations displayed for each block.

From these conflicting results, it is difficult to conclude that inhibition is
a mechanism deployed to erase and/or overwrite features of the stream. It may
be that ResNet18 relies more on higher weight magnitudes to accomplish
overwrite as seen in Figure 3. That said, the negative activations and
correlations in block 3.1 suggest an intuitively appealing strategy by which
DNNs can erase channel features and overwrite them if the block feature is
present.

Figure 5: Percent of channels that meet all three scale invariance criteria
(purple), and those with the mix ratio criterion omitted (yellow).

We believe more fine-grained circuit analyses of the input and block features
could help clarify this issue.

###  3.2 Scale Invariance in the Residual Stream

In this section, the aforementioned three criteria (see Section 2.4) are
applied across all studied blocks to see if any scale invariant channels of
the residual stream rise to the surface. Indeed, in blocks 1.1, 2.0, 2.1, and
3.1, many channels fulfill the criteria. Displayed in Figure 5 are the totals
as a percentage of total block channels. Zero scale invariant channels were
found in blocks 3.0, 4.0, and 4.1. At first it seemed this was due to the
overwrite skew, but only a few were found in 4.0 and 4.1 even with the mix
ratio criterion (Equation 2) removed (see yellow bars in Figure 5). This may
be partly explained by the receptive fields of neurons in 4.0 and 4.1 being
close to the size of the image, which constrains their features to have a
fixed scale. However, 3.0 counters this with an appreciable percentage of
channels, albeit substantially less than its neighboring blocks. It appears
that overall, the reason for the dearth of scale invariance in these blocks
lies more so in the absence of scaled features per se.

We now obtain the top 3 channels of each block with the highest scale metric
$SM$ defined in Equation 5. We visualize these channels in Figure 6 via grids
of their FZs and center-neuron top-activating ImageNet validation images,
together with their mix ratio $M$ and scale metric $SM$. In these examples, we
uncover a variety of features such as small and large black dots (block 2.0
channel 15), rectilinear patterns (block 2.1 channel 108), and oriented
cylindrical curves (block 3.1 channels 113 and 178). Some of these examples
such as the cylinders look like intuitively plausible cases of scale
invariance, but we note others do not appear as obvious, such as those in
block 1.1 (although there are relatively lower scale metrics for channels 19
and 46).

Figure 6: Columns are blocks with channels that passed the scale invariance
criteria (see Section 2.4). Rows are the top 3 channels in descending order by
the scale metric. Each grid has the same layout as in Figure 1. All channel
indices are zero-indexed. Displayed for each grid are the mix ratio and the
two summed deltas of the scale metric per Equation 5.

These results provide preliminary evidence for the existence of scale
invariant features computed by the residual stream. Our criteria may be a good
starting point, but leaves more to be desired. To add robustness to the
results, we reran all three criteria but replaced FZs with the mean activation
for the top 9 activating ImageNet validation images. For all the channels
which cleared both the FZ and ImageNet versions of the criteria, we visualize
them in the appendix (see Section A.3).

##  4 Discussion

###  4.1 Limitations

As this investigation is in a nascent stage, it leaves further development to
be made. While the basic behavior of the skip-overwrite spectrum is presented,
we have not addressed why particular features are skipped, erased/overwritten,
or mixed with another feature. We also leave the overwrite skew of the
downsample blocks unexplored, although the mixes visualized in appendix
Section A.4 for 4.1 indicate significant feature convergence, making the input
features redundant. The implementation of feature erasure/overwrite also
remains unexplained. While there are hints that inhibition may be a mechanism,
our results show mixed findings. A more circuits-driven analysis may reveal if
negatively-weighted features in the final layer of a block (incorporating the
sign of BN’s weight) bear similarities to the positively-weighted features of
the input.

The criteria used to detect scale invariance might be necessary but is
certainly insufficient, probably leading to false positives. A more definitive
method is needed, such as how the circuit approach reveals scale equivariance
in InceptionV1 [11]. Similarly, if scaled versions of the same feature exist
in the input and block, they must also have correspondingly scaled circuitry.
Additionally, further work could better assess if scale invariance is an
explicit strategy learned in the presence of a residual stream. An approach
would be to test object recognition accuracy robustness to scale
transformations against a DNN with a residual stream and an otherwise
identical DNN without one. One could also detect if scale equivariance occurs
intra-layer in the streamless DNN, which is necessary to form scale invariance
in the subsequent layer in a feedforward architecture.

###  4.2 Future Directions

Aside from rectifying the limitations of this study, new directions of
research can branch from this work. We encourage the continued search for
universality in a variety of ways. First, the skip-overwrite spectrum and
scale invariance could be searched for in other residual stream architectures
such as the other ResNets or vision transformers. We consider recurrent neural
networks such as CORNet-S [7] especially intriguing. As weights are shared
across timesteps, this allows for the detection of the identical feature
across different scales, and its bypass connection can build the invariance.
Second and more speculatively, there may be commonalities in how features are
updated in the residual streams of vision and language models. For instance, a
similar feature erasure/overwrite behavior in transformer-based language
models has been alluded to [3]. This direction contains the possibility for
the universality of residual stream management not only in different
architectures, but across entirely different modalities. Finally, as
conjectured in the introduction, these findings may have relevance for visual
neuroscience. Literature has emphasized the existence of “bypass connections”
in the visual object recognition systems of the macaque [6] and mouse [4], but
little work has been done to identify the functions of these connections. As
bypass connections are analogous to the residual stream, it seems plausible
that one such function is to compute scale invariant representations as
described in this work. Progress in this direction may find parallels between
visual processing in artificial and biological neural networks,
mechanistically explaining why DNNs perform so well at predicting neural
responses to visual stimuli in the ventral stream [13].

###  4.3 Conclusion

We have presented observations regarding the function of the residual stream
in ResNet18. We offer evidence for interpreting the residual stream as a
flexible feature manager through the skip-overwrite spectrum, which also
facilitates the computation of scale invariant representations. These findings
are important to further the mission of mechanistically understanding DNNs,
with exciting possibilities to discover if their computations also take place
in biological neural networks.

## References

  * Cammarata et al. [2020] N. Cammarata, G. Goh, S. Carter, L. Schubert, M. Petrov, and C. Olah.  Curve detectors.  _Distill_ , 2020.  doi: 10.23915/distill.00024.003.  https://distill.pub/2020/circuits/curve-detectors.
  * Deng et al. [2009] J. Deng, W. Dong, R. Socher, L.-J. Li, K. Li, and L. Fei-Fei.  Imagenet: A large-scale hierarchical image database.  In _2009 IEEE Conference on Computer Vision and Pattern Recognition_ , pages 248–255, 2009.  doi: 10.1109/CVPR.2009.5206848.
  * Elhage et al. [2021] N. Elhage, N. Nanda, C. Olsson, T. Henighan, N. Joseph, B. Mann, A. Askell, Y. Bai, A. Chen, T. Conerly, N. DasSarma, D. Drain, D. Ganguli, Z. Hatfield-Dodds, D. Hernandez, A. Jones, J. Kernion, L. Lovitt, K. Ndousse, D. Amodei, T. Brown, J. Clark, J. Kaplan, S. McCandlish, and C. Olah.  A mathematical framework for transformer circuits.  _Transformer Circuits Thread_ , 2021.  https://transformer-circuits.pub/2021/framework/index.html.
  * Gămănuţ and Shimaoka [2021] R. Gămănuţ and D. Shimaoka.  Anatomical and functional connectomes underlying hierarchical visual processing in mouse visual system.  _Brain Struct Funct_ , 227(4):1297–1315, Nov. 2021.
  * He et al. [2016] K. He, X. Zhang, S. Ren, and J. Sun.  Deep residual learning for image recognition.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , June 2016.
  * Kravitz et al. [2012] D. J. Kravitz, K. S. Saleem, C. I. Baker, L. G. Ungerleider, and M. Mishkin.  The ventral visual pathway: an expanded neural framework for the processing of object quality.  _Trends Cogn Sci_ , 17(1):26–49, Dec. 2012.
  * Kubilius et al. [2019] J. Kubilius, M. Schrimpf, K. Kar, R. Rajalingham, H. Hong, N. Majaj, E. Issa, P. Bashivan, J. Prescott-Roy, K. Schmidt, et al.  Brain-like object recognition with high-performing shallow recurrent anns.  _Advances in neural information processing systems_ , 32, 2019.
  * Muttenthaler and Hebart [2021] L. Muttenthaler and M. N. Hebart.  Thingsvision: A python toolbox for streamlining the extraction of activations from deep neural networks.  _Frontiers in Neuroinformatics_ , 15:45, 2021.  ISSN 1662-5196.  doi: 10.3389/fninf.2021.679838.  URL https://www.frontiersin.org/article/10.3389/fninf.2021.679838.
  * Olah et al. [2017] C. Olah, A. Mordvintsev, and L. Schubert.  Feature visualization.  _Distill_ , 2017.  doi: 10.23915/distill.00007.  https://distill.pub/2017/feature-visualization.
  * Olah et al. [2020a] C. Olah, N. Cammarata, L. Schubert, G. Goh, M. Petrov, and S. Carter.  Zoom in: An introduction to circuits.  _Distill_ , 2020a.  doi: 10.23915/distill.00024.001.  https://distill.pub/2020/circuits/zoom-in.
  * Olah et al. [2020b] C. Olah, N. Cammarata, C. Voss, L. Schubert, and G. Goh.  Naturally occurring equivariance in neural networks.  _Distill_ , 2020b.  doi: 10.23915/distill.00024.004.  https://distill.pub/2020/circuits/equivariance.
  * Paszke et al. [2019] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala.  Pytorch: An imperative style, high-performance deep learning library.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, _Advances in Neural Information Processing Systems_ , volume 32. Curran Associates, Inc., 2019.  URL https://proceedings.neurips.cc/paper_files/paper/2019/file/bdbca288fee7f92f2bfa9f7012727740-Paper.pdf.
  * Schrimpf et al. [2020] M. Schrimpf, J. Kubilius, H. Hong, N. J. Majaj, R. Rajalingham, E. B. Issa, K. Kar, P. Bashivan, J. Prescott-Roy, F. Geiger, K. Schmidt, D. L. K. Yamins, and J. J. DiCarlo.  Brain-score: Which artificial neural network for object recognition is most brain-like?  _bioRxiv_ , 2020.  doi: 10.1101/407007.  URL https://www.biorxiv.org/content/early/2020/01/02/407007.
  * Schrimpf et al. [2021] M. Schrimpf, I. A. Blank, G. Tuckute, C. Kauf, E. A. Hosseini, N. Kanwisher, J. B. Tenenbaum, and E. Fedorenko.  The neural architecture of language: Integrative modeling converges on predictive processing.  _Proc Natl Acad Sci U S A_ , 118(45), Nov. 2021.
  * Schubert et al. [2021] L. Schubert, C. Voss, N. Cammarata, G. Goh, and C. Olah.  High-low frequency detectors.  _Distill_ , 2021.  doi: 10.23915/distill.00024.005.  https://distill.pub/2020/circuits/frequency-edges.
  * Szegedy et al. [2015] C. Szegedy, W. Liu, Y. Jia, P. Sermanet, S. Reed, D. Anguelov, D. Erhan, V. Vanhoucke, and A. Rabinovich.  Going deeper with convolutions.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , June 2015.

##  Appendix A Appendix

###  A.1 Code Repository

To promote replication and further exploration, we share our code at:

https://github.com/cest-andre/residual-stream-interp

###  A.2 Feature Visualization Regularization

We use an identical regularization method as described in [9], except we use
different levels of jitter dependent on the depth of the layer (we still use
half the initial jitter value for the second jitter transform). This is done
due to the unreliability high jitter has at producing a high activation on the
center neuron once optimization has completed. We sacrificed some
interpretability that jitter offers for more accurate activation values.

We omitted jitter for the entire 1.1 and 2.0 blocks, we used jitter$=$4 for
2.1 second BN ($B$) and the entire 3.0 block, then used the default
jitter$=$16 value for the remainder of the blocks.

###  A.3 Scale Invariance Robustness

To add robustness to our scale invariance findings, we rerun the criteria in
2.4, but replace FZ stimuli $\hat{X}_{c}$ with $\textbf{V}_{c}$, the batch of
nine ImageNet validation images that maximally activated the center neuron in
the target channel $c$. We then replace activations $B_{c}(\hat{X}_{c})$ with
the mean activation to this batch, $\bar{B}_{c}(\textbf{V}_{c})$ (same for
$I$), for both mix ratio and scale deltas. For all blocks studied, we display
the channels that pass both this mean top validation criteria and the original
FZ criteria in Figure 7 (we order them by the scale metric obtained from FZs).

Figure 7: All channels that passed the FZ and ImagetNet validation scale
invariant criteria for all blocks studied (1.1 through 4.1). The channels are
ordered by the scale metric (with respect to FZ stimuli) left to right and top
to bottom. The channels are as follows (all zero-indexed): 1.1: 19 2.0: 7, 15,
53, 3 2.1: 108, 36, 13, 21, 60, 52, 93, 5, 58 3.1: 113, 178, 22, 215, 204, 25,
29, 154, 200.

###  A.4 Top Mixes

Figure 8: Channels with mix ratio $M_{c}$ closest to 1 in all blocks studied.
The channels are ordered left to right. The channels are as follows (all zero-
indexed): 1.1: 44, 27, 55 2.0: 42, 9, 78 2.1: 22, 116, 72 3.0: 105, 250, 204
3.1: 128, 246, 41 4.0: 163, 3, 51 4.1: 154, 68, 496.