††thanks: Corresponding Auhtor: Dibakar Ghosh

# Swarmalators on a ring with uncorrelated pinning

Gourab Kumar Sar  Physics and Applied Mathematics Unit, Indian Statistical
Institute, 203 B. T. Road, Kolkata 700108, India  Kevin O’Keeffe  Senseable
City Lab, Massachusetts Institute of Technology, Cambridge, MA 02139  Dibakar
Ghosh  Physics and Applied Mathematics Unit, Indian Statistical Institute, 203
B. T. Road, Kolkata 700108, India<EMAIL_ADDRESS>(D Ghosh);
<EMAIL_ADDRESS>(G.K.Sar)

###### Abstract

We present a case study of swarmalators (mobile oscillators) which move on a
1D ring and are subject to pinning. Previous work considered the special case
where the pinning in space and the pinning in the phase dimension were
correlated. Here we study the general case where the space and phase pinning
are uncorrelated, both being chosen uniformly at random. This induces several
new effects, such as pinned async, mixed states, and a first order phase
transition. These phenomena may be found in real world swarmalators such as
systems of vinegar eels, Janus matchsticks, electrorotated Quincke rollers or
Japanese tree frogs.

> Pinning, the sticking of flow to defects in the underlying medium, is a
> central topic in statistical physics and nonlinear dynamics. It is key to
> the dynamics of charge density waves, for instance, which give rise to
> useful forms of nonlinear conduction Balandin, Zaitsev-Zotov, and Grüner
> (2021). It is also the mechanism behind the flux pinning of type-II
> superconductors Blatter  _et al._ (1994). Here we study pinning in a new
> context: phase oscillators which swarm around in space and they synchronize
> in time. The competition between the pinning, synchronization, and spatial
> self-assembly of such “swarmalators” gives rise to rich collective behavior
> which we investigate with a mix of numerics and analysis.

##  I Introduction

Synchronization and swarming have been at the forefront of the research
avenues while studying the collective behaviors of coupled systems. After the
pioneering work of Kuramoto in 1975, studies on synchronization took a steep
acceleration and have skyrocketed ever since Kuramoto (1975). Such have been
the impact of these studies that its application is ubiquitously found in
technology from cardiac pacemaker cells Michaels, Matyas, and Jalife (1987);
Verheijck  _et al._ (1998) to Josephson junction arrays Wiesenfeld, Colet, and
Strogatz (1996) and power grid networks Motter  _et al._ (2013); Dorfler and
Bullo (2012) to name a few, apart from its overflowing examples in nature like
flashing of fireflies Buck and Buck (1976); Buck (1988), chorusing frogs
Aihara  _et al._ (2008), movement of crowds Ma  _et al._ (2021), clapping in
unison Néda  _et al._ (2000) etc. On a parallel front, studies on swarming
have also gained attention of the researchers and have not lagged far behind,
especially after the seminal work of Vicsek et al. Vicsek  _et al._ (1995).
Particularly, they have been given core attention when investigating the
movements of animal groups, bird flocks, fish schools, insect swarms etc. Buhl
_et al._ (2006); Couzin (2007); Ballerini  _et al._ (2008); Sumpter (2010);
Herbert-Read (2016). The understanding from these studies has been applied to
technological fields like swarm robotics Ducatelle  _et al._ (2011, 2014),
aerial drones Tahir  _et al._ (2019); Abdelkader  _et al._ (2021) and even for
traffic controlling Ma and He (2019).

Although there is quite a vast understanding of synchronization and swarming
as independent fields, one can not say so when these two fields are brought
into contact. The interplay between swarming and synchronization which is
recently termed as swarmalation lacks theoretical understanding. The co-
occurrence of these two phenomena can be found in biological microswimmers
Yang, Elgeti, and Gompper (2008); Riedel, Kruse, and Howard (2005); Quillen
_et al._ (2021, 2022); Peshkov, McGaffigan, and Quillen (2022); Tamm,
Sonneborn, and Dippell (1975); Verberck (2022); Belovs, Livanovičs, and Cēbers
(2017), driven colloids Yan  _et al._ (2012); Hwang  _et al._ (2020); Zhang,
Sokolov, and Snezhko (2020); Bricard  _et al._ (2015); Zhang  _et al._ (2021);
Manna, Shklyaev, and Balazs (2021); Li  _et al._ (2018); Chaudhary  _et al._
(2014), embryonic cells Tan  _et al._ (2022); Petrungaro, Morelli, and Uriu
(2019) etc. This necessitates careful and rigorous research on entities which
swarm in space and sync over time, commonly known as swarmalators.

Striding on that path, Tanaka et al. came up with a model of chemotactic
oscillators where reaction among the particles were mediated by the
surrounding chemicalsTanaka (2007). However, swarmalator as an entity
representing the bidirectional interaction between spatial and internal
dynamics was first introduced by O’Keeffe et al. O’Keeffe, Hong, and Strogatz
(2017) which paved the way for future works. This two dimensional (2D) model
has been explored with various mechanisms thus far O’Keeffe, Evers, and
Kolokolnikov (2018); Ha  _et al._ (2021); Hong, Yeo, and Lee (2021); Jiménez-
Morales (2020); Ha  _et al._ (2019); Sar  _et al._ (2022); Lizarraga and de
Aguiar (2020); Lee, Yeo, and Hong (2021); Degond, Diez, and Walczak (2022);
Schilcher  _et al._ (2021); Yadav  _et al._ (2023); Ceron, O’Keeffe, and
Petersen (2023); O’Keeffe and Bettstetter (2019); Sar and Ghosh (2022); Kongni
_et al._ (2023); Yadav  _et al._ (2023); Ghosh  _et al._ (2023). Very recently
a one dimensional (1D) counterpart of the 2D swarmalator model has been
proposed by O’Keeffe et al. O’Keeffe, Ceron, and Petersen (2022). This model
is essentially a pair of Kuramoto-like models which is analytically tractable
and thus has captured the well deserved attention of the research community
Yoon  _et al._ (2022); O’Keeffe and Hong (2022); Sar, Ghosh, and O’Keeffe
(2023a); Hong  _et al._ (2023); Lizárraga and de Aguiar (2023); Sar, Ghosh,
and O’Keeffe (2023b); Hao, Zhong, and O’Keeffe (2023); Hong  _et al._ (2023).

In this paper, we aim to study the dynamics of swarmalators on a 1D ring when
random pinning is introduced in the system. Pinning generally refers to a
substance’s propensity to adhere to the irregularities of the underlying
medium, necessitating external force to produce flow. Examples can be laid out
in form of charge density waves Gor’kov and Grüner (2012); Grüner and Zettl
(1985) or a pair of magnetic domain walls Hrabec  _et al._ (2018); Haltz  _et
al._ (2021). In the first example, the phases of the density waves which stand
for $\theta_{i}$, located at positions $x_{i}$ in the underlying lattice,
depin from their designated (pinned) values at a critical forcing amplitude.
For the second case, the center of mass $x_{i}$ and the magnetic dipole vector
$\theta_{i}$ of the magnetic domain walls overcome the pinned configuration
and begin to interact when they are sufficiently forced. These implications
motivate us to employ random pinning in systems consisting of position and
phase interaction.

Swarmalators with spatial and phase pinning have been recently studied by Sar
et al. Sar, Ghosh, and O’Keeffe (2023a) and diverse collective states
including chaos and quasi-periodicity were reported. However, the spatial and
phase pinning sites denoted by $\alpha_{i}$ and $\beta_{i}$, respectively were
chosen linearly as $\alpha_{i}=\beta_{i}=2\pi i/N$, $N$ denoting the number of
swarmalators. This form of pinning sites although helped in the analytical
treatment, is unable to capture the broader scenario where the impurities of
the underlying systems are more often than not randomly distributed Reichhardt
and Reichhardt (2016). To fill this gap, here we have chosen both the spatial
and phase pinning sites randomly and shown the robustness of the previously
reported states as well as the emergence of new ones. We have also bailed
ourselves out of the symmetric spatial and phase coupling strengths and
explored the model when these strengths are not necessarily equal. Different
phase transition phenomena have also been shown varying the driving/forcing
amplitude and coupling strengths. The highlight is a first order phase
transition from a pinned async state to a sync state. We believe our work adds
to the existing results of swarmalators on disordered environment and hope
that it injects impetus to further such explorations.

##  II Model Description

The system we set out to study is a pair of Kuramoto-like models

$\displaystyle\dot{x_{i}}$
$\displaystyle=E+b\sin(\alpha_{i}-x_{i})+\frac{J}{N}\sum_{j}^{N}\sin(x_{j}-x_{i})\cos(\theta_{j}-\theta_{i}),$
(1) $\displaystyle\dot{\theta_{i}}$
$\displaystyle=E+b\sin(\beta_{i}-\theta_{i})+\frac{K}{N}\sum_{j}^{N}\sin(\theta_{j}-\theta_{i})\cos(x_{j}-x_{i}),$
(2)

where $(x_{i},\theta_{i})\in(\mathbb{S}^{1},\mathbb{S}^{1})$ are the position
and phase of the $i$-th swarmalator for $i=1,2,\ldots,N$ and $N$ is the number
of swarmalators. $E$ represents external forcing or strength of driving. $J$,
$K$ are the spatial and phase coupling strengths which control the phase
dependent aggregation and position dependent synchronization, $b$ is the
strength of the pinning. $\alpha_{i}$ and $\beta_{i}$ are the position and
phase pinning sites, respectively which are chosen randomly from $[-\pi,\pi]$.
$b$ $(>0)$ keeps the swarmalators pinned to their pinning sites, whereas, the
coupling strengths $J$ $(>0)$ and $K$ $(>0)$ minimize the differences in both
positions and phase. When $b\gg J,K$, swarmalators remain pinned and when
$J,K\gg b$, they move close to each other and get synchronized in phase. The
driving strength $E$ determines whether the solution will be static or moving
(non-stationary). Without loss of generality, we set the pinning strength
$b=1$ and study our model with three parameters $E$, $J$, and $K$. The
$\alpha_{i}=\beta_{i}$ model with $J=K$ case has been previously studied in
Ref. Sar, Ghosh, and O’Keeffe (2023a). Here we focus on studying the model
where the pinning locations are asymmetric, i.e., $\alpha_{i}\neq\beta_{i}$.
We explore both the cases when the coupling strengths are symmetric ($J=K$)
and asymmetric ($J\neq K$).

We define the following order parameters which are instrumental for studying
the collective states

$\displaystyle W_{\pm}$
$\displaystyle=S_{\pm}e^{{i}\Phi_{\pm}}=\frac{1}{N}\sum_{j}e^{{i}(x_{j}\pm\theta_{j})}.\hskip
10.0pt({i}=\sqrt{-1})$ (3)

$S_{\pm}$ measure the phase-space correlation of the swarmalators. If the
phases and positions are perfectly correlated as $x_{i}=\pm\theta_{i}+C$ for
some constant $C$, then we have $S_{\mp}=1$. Otherwise $S_{\pm}<1$. In the
complete asynchronous state where both the positions and phases are uniformly
distributed between $-\pi$ to $\pi$, we get $S_{\pm}=0$. We define two more
order parameters to measure if the solutions are pinned or not,

$\displaystyle Z_{x}$
$\displaystyle=R_{x}e^{{i}\Psi_{x}}=\frac{1}{N}\sum_{j}e^{{i}(x_{j}-\alpha_{j})},$
(4) $\displaystyle Z_{\theta}$
$\displaystyle=R_{\theta}e^{{i}\Psi_{\theta}}=\frac{1}{N}\sum_{j}e^{{i}(\theta_{j}-\beta_{j})}.$
(5)

By definition, $0<R_{x},R_{\theta}<1$. When the positions of the swarmalators
$x_{i}$ are pinned to the assigned locations $\alpha_{i}$, we get
$R_{x}\approx 1$. Similarly, $R_{\theta}\approx 1$ when
$\theta_{i}\rightarrow\beta_{i}$. We rewrite the model in terms of the order
parameters $S_{\pm}$ as

$\displaystyle\dot{x}_{i}$
$\displaystyle=E+\sin(\alpha_{i}-x_{i})+\frac{J}{2}S_{+}\sin(\Phi_{+}-(x_{i}+\theta_{i}))$
$\displaystyle+\frac{J}{2}S_{-}\sin(\Phi_{-}-(x_{i}-\theta_{i})),$ (6)
$\displaystyle\dot{\theta}_{i}$
$\displaystyle=E+\sin(\beta_{i}-\theta_{i})+\frac{K}{2}S_{+}\sin(\Phi_{+}-(x_{i}+\theta_{i}))$
$\displaystyle-\frac{K}{2}S_{-}\sin(\Phi_{-}-(x_{i}-\theta_{i})).$ (7)

Next, we investigate the emerging collective states of the model and their
properties by changing three parameters $E$, $J$ and $K$.

Figure 1: Collective states in absence of driving ($E=0$). (a)-(d) Pinned
async for $J=K=-1$. (e)-(h) Phase wave of type $x_{i}=-\theta_{i}+C$ for
$J=-1,K=5$. (i)-(l) Phase wave of type $x_{i}=\theta_{i}+C$ for $J=5,K=-1$.
(m)-(p) Sync state for $J=K=5$. (q)-(t) Mixed state for $J=5,K=0.5$. Here,
$(dt,T,N)=(0.1,200,2000)$ and initially $x_{i},\theta_{i}$ were drawn randomly
from $[-\pi,\pi]$. The cyan and magenta circles in the first column refer to
the complex order parameter $W_{+}$ and $W_{-}$, respectively. The pinning
sites are also chosen randomly from $[-\pi,\pi]$. In the first column, the
positions of the swarmalators on the ring are shown where they are colored
according to their phases. The $x$-$\theta$, $x$-$\alpha$ and $\theta$-$\beta$
correlations are depicted in the second, third, and fourth columns,
respectively for all the five states.

##  III Case I: $E=0$

First, we set $E=0$, i.e., explore the collective states when there is no
driving force present in the system. Numerics suggest four different
collective states when the driving strength is zero and we vary the coupling
strengths $J$ and $K$.

1\. Pinned async state: Swarmalators’ positions and phases are pinned to their
pinning locations and we have $x_{i}\rightarrow\alpha_{i}$,
$\theta_{i}\rightarrow\beta_{i}$ for all $i=1,2,\ldots N$. As a result
$(R_{x},R_{\theta})=(1,1)$. Since the pinning locations are randomly
distributed inside $[-\pi,\pi]$, the positions and phases of the swarmalators
are uncorrelated in this state. Which gives $(S_{+},S_{-})=(0,0)$ here. See
Fig. 1 first row for visual demonstration of this state. In Fig. 1(a), the
positions of the swarmalators on the 1D ring is delineated where they are
colored according to their phases. Figure 1(b) shows the position-phase
correlation in the $x$-$\theta$ plane. The relations between the positions and
phases to the pinning sites are portrayed in Figs. 1(c)-(d).

2\. Phase wave state: In this state, the positions and phases of the
swarmalators are found to be correlated. They can either be like
$x_{i}=-\theta_{i}+C$ for some constant $C$ which yields
$(S_{+},S_{-})=(1,0)$, or like $x_{i}=\theta_{i}+C$ where the order parameters
$(S_{+},S_{-})=(0,1)$. These two instances of the phase wave state are
depicted in the second and third columns of Fig. 1, respectively. More
generally, the phase wave state is characterized as the $(S,0)$ or $(0,S)$
state where $S\neq 0$. One of $R_{x}$ and $R_{\theta}$ takes value close to
$1$ and the other is almost zero. The emergence of these two configurations of
the phase wave state depends on the choice of initial conditions.

3\. Sync state: As the name suggests, in here, we observe that the coupling
strengths $J$ and $K$ overcome the pinned arrangement and swarmalators start
to synchronize both their positions and phases. In Fig. 1(m), they can be
found forming two spatial clusters on the ring. The width of the clusters
decreases with increment in the coupling strengths. The $x$-$\theta$
correlation in Fig. 1(n) highlights the fact that both the positions and
phases are synchronized here. Figures. 1(o)-(p) depict the $x$-$\alpha$ and
$\theta$-$\beta$ correlations, respectively. In terms of the order parameters,
the sync state can be represented as $(S_{+},S_{-})=(S,S)$ where $S\neq 0$ and
$R_{x},R_{\theta}\ll 1$.

4\. Mixed state: Alike in the sync state, here also both the order parameters
$S_{\pm}$ are nonzero. But, they are not equal and that is why we call this as
the mixed state. Order parameters characterize this state as
$(S_{+},S_{-})=(S_{1},S_{2})$ where $S_{1},S_{2}\neq 0$ and $S_{1}\neq S_{2}$.
The fifth column of Fig. 1 illustrates this state. Either $x$ or $\theta$ are
completely shifted from their pinning sites ($x$ in Fig. 1(s)) and the other
component is on the verge of losing the pinning configuration ($\theta$ in
Fig. 1(t)). As a result, one of $R_{x}$ and $R_{\theta}$ is much less than $1$
and the other is almost $1$.

The system with no driving ($E=0$) is further studied under two scenarios.
One, when the coupling strengths $J$ and $K$ are equal and two, when they are
unequal.

###  III.1 Symmetric coupling strengths ($J=K$)

With symmetric coupling strengths, we find the occurrence of only those states
where the values of the order parameters $S_{+}$ and $S_{-}$ are equal. These
are the pinned async state and the sync state. Here, we are eventually left
with only one parameter $K(=J)$ and we vary it to study the transition of
states. In Fig. 2, $S_{+}$ and $S_{-}$ are plotted as functions of $K$ by
filled blue circles and void red square markers, respectively. They are joined
by blue and red solid lines and we find that these lines overlap indicating
the emergence of those states with $S_{+}=S_{-}$. From the figure we can see
that there is a phase transition from the pinned async state where
$(S_{+},S_{-})=(0,0)$ to the sync state $(S_{+},S_{-})=(S,S)$ ($S>0$) around
$K=2$. The nature of the transition reveals that there is an abrupt jump of
$S_{\pm}$ near this transition point. Which essentially instructs us to check
for hysteresis which is associated with first order phase transition. In the
inset of Fig. 2, we plot the forward and backward transition curves for
$S_{+}$ ($=S_{-}$) by filled blue circles and magenta diamonds, respectively
over a zoomed in portion near the transition point $K_{c}=2$. In the forward
transition, $S_{+}$ takes a sudden jump from the $(0,0)$ state to the $(S,S)$
state around $K\approx 1.97$. For the backward transition, the $(S,S)$ state
persists till the value of $K$ drops down to $1.86$. The interval between
these $K$ values observes the co-existence of the pinned async and the sync
state resulting in a bistable hysteretic region.

Figure 2: Phase transition with symmetric coupling $J=K$ in absence of
driving. $S_{+}$ and $S_{-}$ are calculated as $K$ is being varied over an
interval. $S_{+}$ (filled blue circles) and $S_{-}$ (void red squares) are
equal in magnitude for all $K$. There is a transition from the
$(S_{+},S_{-})=(0,0)$ (pinned async) state to the $(S,S)$ (sync) state near
$K_{c}=2$. Simulation parameters are $(dt,T,N)=(0.5,200,10^{4})$. In the
inset, hysteretic phase transition is shown. Filled blue circles and magenta
diamonds indicate the values of $S_{+}$ for forward and backward transitions,
respectively. The shaded light blue region inside the forward and backward
transition curves signify the hysteresis region. Simulation parameters:
$(dt,T,N)=(0.5,200,10^{5})$. Time average is taken over last $10\%$ data which
is further averaged over $10$ realizations for all the data points.

Next, we proceed with the analytical treatment of the pinned async state. In
the $N\rightarrow\infty$ limit, $x_{\alpha}$ is assigned to the pinning site
$\alpha$ where we get rid of the suffix in $\alpha_{i}$. Similarly,
$\theta_{\alpha}$ is assigned to the phase pinning site $\Gamma(\alpha)$ where
$\Gamma$ is any map on $\mathbb{S}^{1}$. If $\Gamma(\alpha)=\alpha$ then we
get the symmetric pinning. Under this set up, we reform the governing
equations as

$\displaystyle\dot{x}_{\alpha}$ $\displaystyle=E+\sin(\alpha-
x_{\alpha})+\frac{K}{2}S_{+}\sin(\Phi_{+}-(x_{\alpha}+\theta_{\alpha}))$
$\displaystyle+\frac{K}{2}S_{-}\sin(\Phi_{-}-(x_{\alpha}-\theta_{\alpha})),$
(8) $\displaystyle\dot{\theta}_{\alpha}$
$\displaystyle=E+\sin(\Gamma(\alpha)-\theta_{\alpha})+\frac{K}{2}S_{+}\sin(\Phi_{+}-(x_{\alpha}+\theta_{\alpha}))$
$\displaystyle-\frac{K}{2}S_{-}\sin(\Phi_{-}-(x_{\alpha}-\theta_{\alpha})),$
(9)

where $S_{\pm}$ are rewritten as

$S_{\pm}e^{{i}\Phi_{\pm}}=\frac{1}{2\pi}\int_{0}^{2\pi}e^{{i}(x_{\alpha}\pm\theta_{\alpha})}\,d\alpha.$
(10)

$x:\alpha\mapsto x_{\alpha}$ and $\theta:\alpha\mapsto\theta_{\alpha}$ are
self-maps of the unit circle. In the async state, $S_{\pm}=0$ and the fixed
points are of the form

$\displaystyle x_{\alpha}$ $\displaystyle=\alpha+\sin^{-1}(E),$ (11)
$\displaystyle x_{\alpha}$ $\displaystyle=\Gamma(\alpha)+\sin^{-1}(E).$ (12)

Here, we will study the local stability of the pinned async state. It can be
achieved by diagonalizing the second variation of the potential function
Strogatz  _et al._ (1989) which is of the form

$\displaystyle
P(x_{\alpha},\theta_{\alpha})=-E\int_{0}^{2\pi}x_{\alpha}\,d\alpha\;-E\int_{0}^{2\pi}\theta_{\alpha}\,d\alpha$
$\displaystyle-\int_{0}^{2\pi}\cos(\alpha-
x_{\alpha})\,d\alpha-\int_{0}^{2\pi}\cos(\Gamma(\alpha)-\theta_{\alpha})\,d\alpha$
$\displaystyle-\frac{K}{2\pi}\int_{0}^{2\pi}\int_{0}^{2\pi}\cos(x_{\beta}-x_{\alpha})\cos(\theta_{\beta}-\theta_{\alpha})\,d\alpha
d\beta.$ (13)

We consider a perturbation $\eta:\alpha\mapsto\eta_{\alpha}$ about the pinned
async state

$\displaystyle x_{\alpha}(\epsilon)$
$\displaystyle=\alpha+\sin^{-1}(E)+\epsilon\eta_{\alpha},$ (14)
$\displaystyle\theta_{\alpha}(\epsilon)$
$\displaystyle=\Gamma(\alpha)+\sin^{-1}(E)+\epsilon\eta_{\alpha}.$ (15)

Here $\epsilon$ is small in magnitude. We plug Eqs. (14)-(15) into Eq. (13)
and get

$\displaystyle P(x_{\alpha}$
$\displaystyle(\epsilon),\theta_{\alpha}(\epsilon))=-E\int_{0}^{2\pi}(\alpha+\sin^{-1}(E)+\epsilon\eta_{\alpha})\,d\alpha\;$
$\displaystyle-E\int_{0}^{2\pi}(\Gamma(\alpha)+\sin^{-1}(E)+\epsilon\eta_{\alpha})\,d\alpha\;$
$\displaystyle-2\int_{0}^{2\pi}\cos(\sin^{-1}(E)+\epsilon\eta_{\alpha})\,d\alpha$
$\displaystyle-\frac{K}{4\pi}\int_{0}^{2\pi}\int_{0}^{2\pi}\cos(\beta-\alpha+\epsilon\eta_{\beta}-\epsilon\eta_{\alpha})$
$\displaystyle\qquad\times\cos(\Gamma(\beta)-\Gamma(\alpha)+\epsilon\eta_{\beta}-\epsilon\eta_{\alpha})\,d\alpha
d\beta.$ (16)

We calculate the second variation of $P$ which is a quadratic form

$\Gamma(\eta)=\frac{d^{2}}{d\epsilon^{2}}P(x_{\alpha}(\epsilon))\Big{|}_{\epsilon=0}.$
(17)

After performing the differentiation, we get

$\displaystyle\Omega($
$\displaystyle\eta)=2\sqrt{1-E^{2}}\int_{0}^{2\pi}\eta^{2}_{\alpha}\,d\alpha$
$\displaystyle+\frac{K}{2\pi}\int_{0}^{2\pi}\int_{0}^{2\pi}(\eta_{\beta}-\eta_{\alpha})^{2}\cos(\gamma(\beta)-\gamma(\alpha))\,d\alpha
d\beta,$ (18)

where $\gamma(\alpha)=\alpha+\Gamma(\alpha)$. We work with
$\Gamma(\alpha)=(r-1)\alpha$, $r\in\mathbb{Z}$ so that
$\gamma(\alpha)=r\alpha$. The $(\eta_{\beta}-\eta_{\alpha})^{2}$ term on the
RHS of Eq. (18) is expanded the and further simplified as

$\displaystyle\int_{0}^{2\pi}$
$\displaystyle\int_{0}^{2\pi}(\eta_{\beta}-\eta_{\alpha})^{2}\cos(r\beta-r\alpha)\,d\alpha
d\beta$
$\displaystyle=-2\int_{0}^{2\pi}\int_{0}^{2\pi}\eta_{\beta}\eta_{\alpha}\cos(r\beta-r\alpha)\,d\alpha
d\beta.$ (19)

Now, the $\alpha$ and $\beta$ integrals separate on expanding the term
$\cos(r\beta-r\alpha)$:

$\displaystyle\int_{0}^{2\pi}\int_{0}^{2\pi}\eta_{\beta}\eta_{\alpha}\cos(r\beta-r\alpha)\,d\alpha
d\beta$
$\displaystyle=\int_{0}^{2\pi}\int_{0}^{2\pi}\eta_{\beta}\eta_{\alpha}(\cos{r\alpha}\cos{r\beta}+\sin
r\alpha\sin r\beta)\,d\alpha d\beta$
$\displaystyle=\Bigg{[}\int_{0}^{2\pi}\eta_{\alpha}\cos
r\alpha\,d\alpha\Bigg{]}^{2}+\Bigg{[}\int_{0}^{2\pi}\eta_{\alpha}\sin
r\alpha\,d\alpha\Bigg{]}^{2}$
$\displaystyle=\Bigg{|}\int_{0}^{2\pi}\eta_{\alpha}e^{{i}r\alpha}\,d\alpha\Bigg{|}^{2}=4\pi^{2}\big{|}\hat{\eta}_{\alpha}(-r)\big{|}^{2},$
(20)

where $\hat{\eta}_{\alpha}$ denotes the Fourier transform of $\eta_{\alpha}$
which is defined by

$\hat{\eta}_{\alpha}(m)=\frac{1}{2\pi}\int_{0}^{2\pi}\eta_{\alpha}e^{-{i}m\alpha}\,d\alpha.$
(21)

We move to the Hilbert space $L^{2}(\mathbb{S}^{1})$ with the inner product
defined on it as,

$\mu_{\alpha}\cdot\nu_{\alpha}=\frac{1}{2\pi}\int_{0}^{2\pi}\mu_{\alpha}\nu_{\alpha}\,d\alpha.$
(22)

Let us choose $\mu_{\alpha}=\cos r\alpha$, $\nu_{\alpha}=\sin r\alpha$ and we
have $||\mu_{\alpha}||^{2}=||\nu_{\alpha}||^{2}=1/2$ and
$\mu_{\alpha}\cdot\nu_{\alpha}=0$. We consider $\eta_{\alpha}^{\perp}$, a
function which is orthogonal to both $\mu_{\alpha}$ and $\nu_{\alpha}$. This
means
$\mu_{\alpha}\cdot\eta_{\alpha}^{\perp}=\nu_{\alpha}\cdot\eta_{\alpha}^{\perp}=0$.
Now, $\eta_{\alpha}$ can be expressed as a linear combination of
$\mu_{\alpha}$, $\nu_{\alpha}$, $\eta_{\alpha}^{\perp}$ as:

$\eta_{\alpha}=p\frac{\mu_{\alpha}}{||\mu_{\alpha}||}+q\frac{\nu_{\alpha}}{||\nu_{\alpha}||}+\eta_{\alpha}^{\perp},$
(23)

Also, we get

$||\eta_{\alpha}||^{2}=p^{2}+q^{2}+||\eta_{\alpha}^{\perp}||^{2}.$ (24)

With these tools, it is elementary to see that

$\displaystyle\big{|}\hat{\eta}_{\alpha}(-r)\big{|}^{2}$
$\displaystyle=\Bigg{|}\frac{1}{2\pi}\int_{0}^{2\pi}\eta_{\alpha}e^{{i}r\alpha}\,d\alpha\Bigg{|}^{2}$
$\displaystyle=\Bigg{[}\int_{0}^{2\pi}\eta_{\alpha}\cos
r\alpha\,d\alpha\Bigg{]}^{2}+\Bigg{[}\int_{0}^{2\pi}\eta_{\alpha}\sin
r\alpha\,d\alpha\Bigg{]}^{2}$
$\displaystyle=(\eta_{\alpha}\cdot\mu_{\alpha})^{2}+(\eta_{\alpha}\cdot\nu_{\alpha})^{2}$
$\displaystyle=(p||\mu_{\alpha}||)^{2}+(q||\nu_{\alpha}||)^{2}=\frac{p^{2}+q^{2}}{2}.$
(25)

After substitutions and simplification, from Eq. (18) we finally get

$\displaystyle\Omega(\eta)$
$\displaystyle=4\pi\Bigg{[}\sqrt{1-E^{2}}(p^{2}+q^{2}+||\eta_{\alpha}^{\perp}||^{2})-K\Big{(}\frac{p^{2}+q^{2}}{2}\Big{)}\Bigg{]}$
$\displaystyle=4\pi\Bigg{[}(p^{2}+q^{2})\Big{(}\sqrt{1-E^{2}}-\frac{K}{2}\Big{)}+\sqrt{1-E^{2}}||\eta_{\alpha}^{\perp}||^{2}\Bigg{]}.$
(26)

So, $\Omega$ is positive definite when $\sqrt{1-E^{2}}-K/2>0$. From this we
get the stability boundary

$\displaystyle E_{c}(K)=\sqrt{1-\frac{K^{2}}{4}}.$ (27)

When $E=0$, we get $K_{c}=2$ which matches the transition point found in Fig.
2. The small deviation from the analytically predicted value is due to the
finite system size and reduces with increasing $N$.

###  III.2 Asymmetric coupling strengths ($J\neq K$)

We found only two collective states out of the four mentioned at the beginning
of the section with symmetric coupling strengths $J=K$. The rest of the two
states which are the phase wave and mixed state can only be seen when $J$ and
$K$ are unequal. For a particular case, we fix $J=4$ and inspect the different
phase transition behaviors while varying $K$. The result is displayed in Fig.
3. Starting from the left with negative $K$, first we encounter the pinned
async state where both $S_{+}$ and $S_{-}$ lie near zero. This is highlighted
by the light red region in Fig. 3. Moving to the right with increasing $K$, we
find that pinned async state becomes unstable near $K\approx-0.5$ where
$S_{+}$ starts to increase from zero while $S_{-}$ continue to stay near zero
(look at the light yellow region in Fig. 3). This is the $(S,0)$ state or the
phase wave as we have discussed earlier. It is to be noted that, the phase
wave state in this region can be either of the $(S,0)$ or $(0,S)$ types
depending on the initial condition. We set the maximum of $S_{+}$ and $S_{-}$
to $S_{+}$ so that it is always of type $(S,0)$ during numerical simulation.
With increasing $K$, the phase wave state loses its stability around $K\approx
0.3$. From the phase wave state, $S_{+}$ decreases and $S_{-}$ increases
resulting in the occurrence of the $(S_{1},S_{2})$ case where $S_{1}>S_{2}$.
This is the mixed state which is indicated by the light pink region in Fig. 3.
From here, $S_{+}$ further decreases while $S_{-}$ increases until they
collide around $K\approx 0.5$. Finally, we observe the sync state
characterized by $(S,S)$ where $S\neq 0$. The interval of occurrence of the
sync state is delineated by the light blue region in Fig. 3.

Figure 3: Phase transition scenario with asymmetric coupling strengths $J\neq
K$ for $E=0$. Filled blue circles and void red markers represent the values of
$S_{+}$ and $S_{-}$, respectively. The regions corresponding to the emergence
of different states are colored by light red (pinned async), light yellow
(phase wave), light pink (mixed) and light blue (sync). Simulation parameters
used here are $(dt,T,N)=(0.1,200,10^{4})$. Both $S_{+}$ and $S_{-}$ are time
averaged over last $10\%$ data and further averaged over $10$ realizations. We
have fixed $J=4$ here. The phase transition here follows the route pinned
async $(0,0)$ $\rightarrow$ phase wave $(S,0)$ $\rightarrow$ mixed
$(S_{1},S_{2})$ $\rightarrow$ sync $(S,S)$ with increasing $K$.

We have also looked for the bifurcation in the $J$-$K$ parameter plane. We
demonstrate the diagram in Fig. 4. We numerically calculate this by dividing
the plane into $100$ by $100$ mesh points. At each of these point, the values
of the order parameters $S_{+}$, $S_{-}$ are calculated. We color the mesh
points by red, yellow, pink and blue if the emerging states at those points
are pinned async, phase wave, mixed and sync, respectively. We have
distinguished the states as below

$\displaystyle\text{Pinned async}:\\{S_{+}<0.2\\}\cap\\{S_{-}<0.2\\}.$
$\displaystyle\text{Phase wave}:\\{S_{+}>0.2\\}\cap\\{S_{-}<0.2\\}.$
$\displaystyle\text{Mixed}:\\{S_{+}>0.2\\}\cap\\{S_{-}>0.2\\}\cap\\{S_{+}-S_{-}>0.1\\}.$
$\displaystyle\text{Sync}:\\{S_{+}>0.2\\}\cap\\{S_{-}>0.2\\}\cap\\{S_{+}-S_{-}<0.1\\}.$

Based on the above conditions, we compute the bifurcation diagram in Fig. 4.
We find the existence of pinned async state (red region) when $J+K\lessapprox
3$. The phase wave state (yellow region) is encountered when one of $J$ and
$K$ is positively large and the other is small. The mixed state (pink region)
is seen at the boundary between the phase wave and the sync state and their
region of occurrence is small compared to the other states. When both the
coupling strengths $J$ and $K$ are positive and large, the sync state (blue
region) prevails.

Figure 4: Phase diagram in the $J$-$K$ plane for $E=0$. Red, yellow, pink and
blue indicate the pinned async, phase wave (PW), mixed and sync states,
respectively. We use the following simulation parameters:
$(dt,T,N)=(0.1,500,1000)$. First $90\%$ data were discarded as transients.
Each data point is an average of $5$ realizations. We have characterized the
states in terms of the order parameters $S_{\pm}$. It can be seen that one of
$J$ and $K$ has to be positively large for the phase wave state to be
achieved.

##  IV Case II: $E\neq 0$

With a sound understanding of the collective states in the $E=0$ limit, we now
move to explore the dynamical features of the system in presence of driving,
i.e., in the $E\neq 0$ limit. In the undriven limit, all the states were
stationary. The mean or average velocity $\langle
V\rangle=N^{-1}\sum_{i}\sqrt{v_{x,i}^{2}+v_{\theta,i}^{2}}$ was zero there.
But, in presence of driving in the system, we expect moving solutions where
$\langle V\rangle\neq 0$. Looking at the fixed points given by Eqs. (11)-
(12), we can say that they exist when $E<1$. From numerics, we observe that
all the solutions for $E>1$ are nonstationary in nature, i.e., $\langle
V\rangle>0$.

Another stark contrast with the $E=0$ case is that even with the symmetric
coupling strengths $J=K$ we find the existence of phase wave and mixed states.
Moreover, a new state is born where the positions and phases of the
swarmalators are distributed inside $[-\pi,\pi]$ but unlike in the pinned
async state they lose their connection with the pinning locations. In terms of
the order parameter, in this state we have $S_{\pm}\approx 0$ and
$R_{x},R_{\theta}\ll 1$. This state is called as the async state. To capture
the entire scenario, we move to the $K$-$E$ parameter space. Figure 5 captures
the overall behavior of the system in the undriven limit with symmetric
coupling strengths. Along with the aforementioned conditions to distinguish
the pinned async, phase wave, mixed and sync states, here we write down the
condition for realization of the async state as

$\displaystyle\text{Async}:\\{S_{+}<0.2\\}\cap\\{S_{-}<0.2\\}\cap\\{R_{x}<0.2\\}\cap\\{R_{\theta}<0.2\\}.$

Figure 5: Phase diagram in the $K$-$E$ plane for $J=K$. Red, magenta, yellow,
pink and blue colored regions correspond to the occurrence of pinned async,
async, phase wave, mixed and sync states, respectively. Black dashed line is
the $E=1$ line. Solid black line is the analytical prediction Eq. (27). Here,
we have used $(dt,T,N)=(0.1,500,1000)$. First $90\%$ data were discarded as
transients. Each data point is an average of $5$ realizations. We see that
with negative $K$ the pinned async state becomes unstable as $E$ crosses 1
from below and the nonstationary async state takes over. The phase wave and
mixed states are observed over relatively small regions near the boundary
between the pinned async and sync states.

In Fig. 5, the black dashed line indicates $E=1$ beyond which the stationary
pinned async state (red region) loses its stability and the nonstationary
async state (magenta region) takes place. The solid black line signifies the
analytical boundary Eq. (27) between the pinned async state and the sync state
(blue region). We also observe the scattered presence of the phase wave
(yellow) and the mixed (pink) states.

Figure 6: Phase transition with asymmetric coupling $J\neq K$ with $E\neq 0$.
Filled blue circles and void red squares correspond to the values of $S_{+}$
and $S_{-}$, respectively. The transition phenomena is shown in (a) for
$E=0.5$, in (b) for $E=1.5$ and in (c) for $E=2.5$. Simulation parameters:
$(dt,T,N)=(0.5,200,10^{4})$. Both the order parameters are time averaged over
last $10\%$ data and averaged over $5$ realizations. Here, we have used
$J=4.0$. In (a), we find pinned async $\rightarrow$ phase wave $\rightarrow$
mixed $\rightarrow$ sync, whereas in (b) and (c) we observe async
$\rightarrow$ phase wave $\rightarrow$ sync. We also observe first-order
transition taking place from phase wave state to sync state when $E$ is large.

We continue our study with asymmetric coupling strengths $J\neq K$. Again for
instance we fix $J=4$. We accumulate the transition phenomena for different
$E$ values in Fig. 6 where we vary $K$ over an interval $[-6,2]$. For $E=0.5$
in Fig. 6(a), the transition behavior is similar to the one we found in Fig. 3
with $E=0$. The route with increasing $K$ is noted down as: pinned async
$\rightarrow$ phase wave $\rightarrow$ mixed $\rightarrow$ sync. An increase
to the value of $E$ to $1.5$ shows the disappearance of the stationary pinned
async state and arrival of the nonstationary async state in place of it.
Moreover, in Fig. 6(b) we observe that the interval for phase wave state
expands compared to the $E=0.5$ case in Fig. 6(a). The mixed state which
exists over a tiny interval for $E=0.5$ vanishes completely for $E=1.5$. As a
result, the transition route is: async $\rightarrow$ phase wave $\rightarrow$
sync. With $E=2.5$ in Fig. 6(c), the dynamical features do not change much
from the $E=1.5$ scenario. However, the phase wave state is more prominent
here. We also see that with higher driving strength as the mixed state
vanishes, the transition from the phase wave state to the sync state becomes
abrupt.

##  V Discussion

In this work, we have broken the extreme symmetries of the pinning model by
taking both the pinning locations randomly. This resulted in the emergence of
a pinned async state where $x_{i}\rightarrow\alpha_{i}$ and
$\theta_{i}\rightarrow\beta_{i}$ which trivially transforms into the phase
wave state if we assume $\alpha_{i}=\beta_{i}$. On top of it, with the
asymmetric random pinning sites, we encounter phase wave states that arise
when one of the position and phase variables remains pinned and the other
depins. Under sufficiently large coupling strengths, sync state is observed
where positions and phases of the swarmalators are synchronized. At an
intermediate stage between the phase wave and sync state, the mixed state
develops for small $E$ and disappears when $E$ is large.

Randomness is always difficult to deal with while treating it analytically
compared to systematic patterns. While the random pinning locations make our
model realistic and produce unknown collective states, they take away from us
the luxury of deterministic mathematical framework. That is why we had to keep
ourselves limited to numerical treatment of majority of the collective states.
This being said, the collective states along with their counterparts in real
world have enough to offer to this line of study.

Our model captures the behavior of real world swarmalators like circularly
confined spermatazoa Creppy  _et al._ (2016) and bordertaxic vinegar eels
Quillen  _et al._ (2021, 2022); Peshkov, McGaffigan, and Quillen (2022). When
combined with the findings of symmetric pinning Sar, Ghosh, and O’Keeffe
(2023a), the model’s dynamical states can analogously be found in models of
Japanese tree frogs Aihara  _et al._ (2014) and Janus matchsticks Chaudhary
_et al._ (2014). The 1D ring model also has implications in groups of sperm
Creppy  _et al._ (2016) and forced colloids Yan  _et al._ (2012); Yan, Bae,
and Granick (2015); Zhang, Sokolov, and Snezhko (2020).

In a nutshell, our model establishes dynamics of swarmalators with random
pinning locations with both symmetric and asymmetric coupling strengths.
Future research works can look for an analytical treatment of these states in
terms of the variation of order parameters, $\dot{S}_{\pm}$. The driving
strength $E$ can also be chosen asymmetrically in the position and phase
components, so is the pinning strength $b$. Local coupling can be introduced
in the model. This work can also be extended to swarmalators where the spatial
movements are in the two dimensional plane. The coupling strengths, driving
force and the pinning strength which are identical (i.e., chosen from a delta
distribution) here for all the swarmalators, can be chosen from some other
kind of distributions as well. Periodic forcing/driving strengths can also be
considered in place of the constant driving.

## Data Availability

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

## References

  * Balandin, Zaitsev-Zotov, and Grüner (2021) A. A. Balandin, S. V. Zaitsev-Zotov, and G. Grüner, “Charge-density-wave quantum materials and devices—new developments and future prospects,” Applied Physics Letters 119 (2021).
  * Blatter  _et al._ (1994) G. Blatter, M. V. Feigel’man, V. B. Geshkenbein, A. I. Larkin, and V. M. Vinokur, “Vortices in high-temperature superconductors,” Reviews of Modern Physics 66, 1125 (1994).
  * Kuramoto (1975) Y. Kuramoto, “International symposium on mathematical problems in theoretical physics,” Lecture notes in Physics 30, 420 (1975).
  * Michaels, Matyas, and Jalife (1987) D. C. Michaels, E. P. Matyas, and J. Jalife, “Mechanisms of sinoatrial pacemaker synchronization: a new hypothesis.” Circulation Research 61, 704–714 (1987).
  * Verheijck  _et al._ (1998) E. E. Verheijck, R. Wilders, R. W. Joyner, D. A. Golod, R. Kumar, H. J. Jongsma, L. N. Bouman, and A. C. v. Ginneken, “Pacemaker synchronization of electrically coupled rabbit sinoatrial node cells,” The Journal of General Physiology 111, 95–112 (1998).
  * Wiesenfeld, Colet, and Strogatz (1996) K. Wiesenfeld, P. Colet, and S. H. Strogatz, “Synchronization transitions in a disordered josephson series array,” Physical Review Letters 76, 404 (1996).
  * Motter  _et al._ (2013) A. E. Motter, S. A. Myers, M. Anghel, and T. Nishikawa, “Spontaneous synchrony in power-grid networks,” Nature Physics 9, 191–197 (2013).
  * Dorfler and Bullo (2012) F. Dorfler and F. Bullo, “Synchronization and transient stability in power networks and nonuniform kuramoto oscillators,” SIAM Journal on Control and Optimization 50, 1616–1642 (2012).
  * Buck and Buck (1976) J. Buck and E. Buck, “Synchronous fireflies,” Scientific American 234, 74–85 (1976).
  * Buck (1988) J. Buck, “Synchronous rhythmic flashing of fireflies. ii.” The Quarterly Review of Biology 63, 265–289 (1988).
  * Aihara  _et al._ (2008) I. Aihara, H. Kitahata, K. Yoshikawa, and K. Aihara, “Mathematical modeling of frogs’ calling behavior and its possible application to artificial life and robotics,” Artificial Life and Robotics 12, 29–32 (2008).
  * Ma  _et al._ (2021) Y. Ma, E. W. M. Lee, M. Shi, and R. K. K. Yuen, “Spontaneous synchronization of motion in pedestrian crowds of different densities,” Nature Human Behaviour 5, 447–457 (2021).
  * Néda  _et al._ (2000) Z. Néda, E. Ravasz, T. Vicsek, Y. Brechet, and A.-L. Barabási, “Physics of the rhythmic applause,” Physical Review E 61, 6987 (2000).
  * Vicsek  _et al._ (1995) T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet, “Novel type of phase transition in a system of self-driven particles,” Physical Review Letters 75, 1226 (1995).
  * Buhl  _et al._ (2006) J. Buhl, D. J. Sumpter, I. D. Couzin, J. J. Hale, E. Despland, E. R. Miller, and S. J. Simpson, “From disorder to order in marching locusts,” Science 312, 1402–1406 (2006).
  * Couzin (2007) I. Couzin, “Collective minds,” Nature 445, 715–715 (2007).
  * Ballerini  _et al._ (2008) M. Ballerini, N. Cabibbo, R. Candelier, A. Cavagna, E. Cisbani, I. Giardina, V. Lecomte, A. Orlandi, G. Parisi, A. Procaccini, _et al._ , “Interaction ruling animal collective behavior depends on topological rather than metric distance: Evidence from a field study,” Proceedings of the National Academy of Sciences 105, 1232–1237 (2008).
  * Sumpter (2010) D. J. Sumpter,  _Collective animal behavior_ (Princeton University Press, 2010).
  * Herbert-Read (2016) J. E. Herbert-Read, “Understanding how animal groups achieve coordinated movement,” Journal of Experimental Biology 219, 2971–2983 (2016).
  * Ducatelle  _et al._ (2011) F. Ducatelle, G. A. Di Caro, C. Pinciroli, and L. M. Gambardella, “Self-organized cooperation between robotic swarms,” Swarm Intelligence 5, 73–96 (2011).
  * Ducatelle  _et al._ (2014) F. Ducatelle, G. A. Di Caro, A. Förster, M. Bonani, M. Dorigo, S. Magnenat, F. Mondada, R. O’Grady, C. Pinciroli, P. Rétornaz, _et al._ , “Cooperative navigation in robotic swarms,” Swarm Intelligence 8, 1–33 (2014).
  * Tahir  _et al._ (2019) A. Tahir, J. Böling, M.-H. Haghbayan, H. T. Toivonen, and J. Plosila, “Swarms of unmanned aerial vehicles—a survey,” Journal of Industrial Information Integration 16, 100106 (2019).
  * Abdelkader  _et al._ (2021) M. Abdelkader, S. Güler, H. Jaleel, and J. S. Shamma, “Aerial swarms: Recent applications and challenges,” Current Robotics Reports 2, 309–320 (2021).
  * Ma and He (2019) C. Ma and R. He, “Green wave traffic control system optimization based on adaptive genetic-artificial fish swarm algorithm,” Neural Computing and Applications 31, 2073–2083 (2019).
  * Yang, Elgeti, and Gompper (2008) Y. Yang, J. Elgeti, and G. Gompper, “Cooperation of sperm in two dimensions: synchronization, attraction, and aggregation through hydrodynamic interactions,” Physical Review E 78, 061903 (2008).
  * Riedel, Kruse, and Howard (2005) I. H. Riedel, K. Kruse, and J. Howard, “A self-organized vortex array of hydrodynamically entrained sperm cells,” Science 309, 300–303 (2005).
  * Quillen  _et al._ (2021) A. Quillen, A. Peshkov, E. Wright, and S. McGaffigan, “Metachronal waves in concentrations of swimming turbatrix aceti nematodes and an oscillator chain model for their coordinated motions,” Physical Review E 104, 014412 (2021).
  * Quillen  _et al._ (2022) A. Quillen, A. Peshkov, B. Chakrabarti, N. Skerrett, S. McGaffigan, and R. Zapiach, “Fluid circulation driven by collectively organized metachronal waves in swimming t. aceti nematodes,” Physical Review E 106, 064401 (2022).
  * Peshkov, McGaffigan, and Quillen (2022) A. Peshkov, S. McGaffigan, and A. C. Quillen, “Synchronized oscillations in swarms of nematode turbatrix aceti,” Soft Matter 18, 1174–1182 (2022).
  * Tamm, Sonneborn, and Dippell (1975) S. L. Tamm, T. Sonneborn, and R. V. Dippell, “The role of cortical orientation in the control of the direction of ciliary beat in paramecium.” The Journal of Cell Biology 64, 98–112 (1975).
  * Verberck (2022) B. Verberck, “Wavy worms,” Nature Physics 18, 131–131 (2022).
  * Belovs, Livanovičs, and Cēbers (2017) M. Belovs, R. Livanovičs, and A. Cēbers, “Synchronized rotation in swarms of magnetotactic bacteria,” Physical Review E 96, 042408 (2017).
  * Yan  _et al._ (2012) J. Yan, M. Bloom, S. C. Bae, E. Luijten, and S. Granick, “Linking synchronization to self-assembly using magnetic janus colloids,” Nature 491, 578–581 (2012).
  * Hwang  _et al._ (2020) S. Hwang, T. D. Nguyen, S. Bhaskar, J. Yoon, M. Klaiber, K. J. Lee, S. C. Glotzer, and J. Lahann, “Cooperative switching in large-area assemblies of magnetic janus particles,” Advanced Functional Materials 30, 1907865 (2020).
  * Zhang, Sokolov, and Snezhko (2020) B. Zhang, A. Sokolov, and A. Snezhko, “Reconfigurable emergent patterns in active chiral fluids,” Nature Communications 11, 1–9 (2020).
  * Bricard  _et al._ (2015) A. Bricard, J.-B. Caussin, D. Das, C. Savoie, V. Chikkadi, K. Shitara, O. Chepizhko, F. Peruani, D. Saintillan, and D. Bartolo, “Emergent vortices in populations of colloidal rollers,” Nature Communications 6, 1–8 (2015).
  * Zhang  _et al._ (2021) B. Zhang, H. Karani, P. M. Vlahovska, and A. Snezhko, “Persistence length regulates emergent dynamics in active roller ensembles,” Soft Matter (2021).
  * Manna, Shklyaev, and Balazs (2021) R. K. Manna, O. E. Shklyaev, and A. C. Balazs, “Chemical pumps and flexible sheets spontaneously form self-regulating oscillators in solution,” Proceedings of the National Academy of Sciences 118, e2022987118 (2021).
  * Li  _et al._ (2018) M. Li, M. Brinkmann, I. Pagonabarraga, R. Seemann, and J.-B. Fleury, “Spatiotemporal control of cargo delivery performed by programmable self-propelled janus droplets,” Communications Physics 1, 1–8 (2018).
  * Chaudhary  _et al._ (2014) K. Chaudhary, J. J. Juárez, Q. Chen, S. Granick, and J. A. Lewis, “Reconfigurable assemblies of janus rods in ac electric fields,” Soft Matter 10, 1320–1324 (2014).
  * Tan  _et al._ (2022) T. H. Tan, A. Mietke, J. Li, Y. Chen, H. Higinbotham, P. J. Foster, S. Gokhale, J. Dunkel, and N. Fakhri, “Odd dynamics of living chiral crystals,” Nature 607, 287–293 (2022).
  * Petrungaro, Morelli, and Uriu (2019) G. Petrungaro, L. G. Morelli, and K. Uriu, “Information flow in the presence of cell mixing and signaling delays during embryonic development,” in  _Seminars in cell & developmental biology_, Vol. 93 (Elsevier, 2019) pp. 26–35.
  * Tanaka (2007) D. Tanaka, “General chemotactic model of oscillators,” Physical Review Letters 99, 134103 (2007).
  * O’Keeffe, Hong, and Strogatz (2017) K. P. O’Keeffe, H. Hong, and S. H. Strogatz, “Oscillators that sync and swarm,” Nature Communications 8, 1–13 (2017).
  * O’Keeffe, Evers, and Kolokolnikov (2018) K. P. O’Keeffe, J. H. M. Evers, and T. Kolokolnikov, “Ring states in swarmalator systems,” Phys. Rev. E 98, 022203 (2018).
  * Ha  _et al._ (2021) S.-Y. Ha, J. Jung, J. Kim, J. Park, and X. Zhang, “A mean-field limit of the particle swarmalator model,” Kinetic and Related Models 14, 429–468 (2021).
  * Hong, Yeo, and Lee (2021) H. Hong, K. Yeo, and H. K. Lee, “Coupling disorder in a population of swarmalators,” Physical Review E 104, 044214 (2021).
  * Jiménez-Morales (2020) F. Jiménez-Morales, “Oscillatory behavior in a system of swarmalators with a short-range repulsive interaction,” Physical Review E 101, 062202 (2020).
  * Ha  _et al._ (2019) S.-Y. Ha, J. Jung, J. Kim, J. Park, and X. Zhang, “Emergent behaviors of the swarmalator model for position-phase aggregation,” Mathematical Models and Methods in Applied Sciences 29, 2225–2269 (2019).
  * Sar  _et al._ (2022) G. K. Sar, S. N. Chowdhury, M. Perc, and D. Ghosh, “Swarmalators under competitive time-varying phase interactions,” New Journal of Physics 24, 043004 (2022).
  * Lizarraga and de Aguiar (2020) J. U. Lizarraga and M. A. de Aguiar, “Synchronization and spatial patterns in forced swarmalators,” Chaos: An Interdisciplinary Journal of Nonlinear Science 30, 053112 (2020).
  * Lee, Yeo, and Hong (2021) H. K. Lee, K. Yeo, and H. Hong, “Collective steady-state patterns of swarmalators with finite-cutoff interaction distance,” Chaos: An Interdisciplinary Journal of Nonlinear Science 31, 033134 (2021).
  * Degond, Diez, and Walczak (2022) P. Degond, A. Diez, and A. Walczak, “Topological states and continuum model for swarmalators without force reciprocity,” Analysis and applications 20, 1215–1270 (2022).
  * Schilcher  _et al._ (2021) U. Schilcher, J. F. Schmidt, A. Vogell, and C. Bettstetter, “Swarmalators with stochastic coupling and memory,” in  _2021 IEEE International Conference on Autonomic Computing and Self-Organizing Systems (ACSOS)_ (2021) pp. 90–99.
  * Yadav  _et al._ (2023) A. Yadav, V. Chandrasekar, W. Zou, J. Kurths, D. Senthilkumar, _et al._ , “Exotic swarming dynamics of high-dimensional swarmalators,” arXiv preprint arXiv:2308.03803 (2023).
  * Ceron, O’Keeffe, and Petersen (2023) S. Ceron, K. O’Keeffe, and K. Petersen, “Diverse behaviors in non-uniform chiral and non-chiral swarmalators,” Nature Communications 14, 940 (2023).
  * O’Keeffe and Bettstetter (2019) K. O’Keeffe and C. Bettstetter, “A review of swarmalators and their potential in bio-inspired computing,” in  _Micro-and Nanotechnology Sensors, Systems, and Applications XI_ , Vol. 10982 (International Society for Optics and Photonics, 2019) p. 109822E.
  * Sar and Ghosh (2022) G. K. Sar and D. Ghosh, “Dynamics of swarmalators: A pedagogical review,” Europhysics Letters 139, 53001 (2022).
  * Kongni  _et al._ (2023) S. J. Kongni, V. Nguefoue, T. Njougouo, P. Louodop, F. F. Ferreira, R. Tchitnga, and H. A. Cerdeira, “Phase transitions on a multiplex of swarmalators,” Physical Review E 108, 034303 (2023).
  * Ghosh  _et al._ (2023) S. Ghosh, G. K. Sar, S. Majhi, and D. Ghosh, “Antiphase synchronization in a population of swarmalators,” Physical Review E 108, 034217 (2023).
  * O’Keeffe, Ceron, and Petersen (2022) K. O’Keeffe, S. Ceron, and K. Petersen, “Collective behavior of swarmalators on a ring,” Phys. Rev. E 105, 014211 (2022).
  * Yoon  _et al._ (2022) S. Yoon, K. P. O’Keeffe, J. F. F. Mendes, and A. V. Goltsev, “Sync and swarm: Solvable model of nonidentical swarmalators,” Phys. Rev. Lett. 129, 208002 (2022).
  * O’Keeffe and Hong (2022) K. O’Keeffe and H. Hong, “Swarmalators on a ring with distributed couplings,” Phys. Rev. E 105, 064208 (2022).
  * Sar, Ghosh, and O’Keeffe (2023a) G. K. Sar, D. Ghosh, and K. O’Keeffe, “Pinning in a system of swarmalators,” Physical Review E 107, 024215 (2023a).
  * Hong  _et al._ (2023) H. Hong, K. P. O’Keeffe, J. S. Lee, and H. Park, “Swarmalators with thermal noise,” Physical Review Research 5, 023105 (2023).
  * Lizárraga and de Aguiar (2023) J. U. Lizárraga and M. A. de Aguiar, “Synchronization of sakaguchi swarmalators,” Physical Review E 108, 024212 (2023).
  * Sar, Ghosh, and O’Keeffe (2023b) G. K. Sar, D. Ghosh, and K. O’Keeffe, “Solvable model of driven matter with pinning,” arXiv preprint arXiv:2306.09589 (2023b).
  * Hao, Zhong, and O’Keeffe (2023) B. Hao, M. Zhong, and K. O’Keeffe, “Mixed sign interactions in the 1d swarmalator model,” arXiv preprint arXiv:2309.02342 (2023).
  * Gor’kov and Grüner (2012) L. P. Gor’kov and G. Grüner,  _Charge density waves in solids_ (Elsevier, 2012).
  * Grüner and Zettl (1985) G. Grüner and A. Zettl, “Charge density wave conduction: A novel collective transport phenomenon in solids,” Physics Reports 119, 117–232 (1985).
  * Hrabec  _et al._ (2018) A. Hrabec, V. Křižáková, S. Pizzini, J. Sampaio, A. Thiaville, S. Rohart, and J. Vogel, “Velocity enhancement by synchronization of magnetic domain walls,” Physical Review Letters 120, 227204 (2018).
  * Haltz  _et al._ (2021) E. Haltz, S. Krishnia, L. Berges, A. Mougin, and J. Sampaio, “Domain wall dynamics in antiferromagnetically coupled double-lattice systems,” Physical Review B 103, 014444 (2021).
  * Reichhardt and Reichhardt (2016) C. Reichhardt and C. O. Reichhardt, “Depinning and nonequilibrium dynamic phases of particle assemblies driven over random and ordered substrates: a review,” Reports on Progress in Physics 80, 026501 (2016).
  * Strogatz  _et al._ (1989) S. H. Strogatz, C. M. Marcus, R. M. Westervelt, and R. E. Mirollo, “Collective dynamics of coupled oscillators with random pinning,” Physica D: Nonlinear Phenomena 36, 23–50 (1989).
  * Creppy  _et al._ (2016) A. Creppy, F. Plouraboué, O. Praud, X. Druart, S. Cazin, H. Yu, and P. Degond, “Symmetry-breaking phase transitions in highly concentrated semen,” Journal of The Royal Society Interface 13, 20160575 (2016).
  * Aihara  _et al._ (2014) I. Aihara, T. Mizumoto, T. Otsuka, H. Awano, K. Nagira, H. G. Okuno, and K. Aihara, “Spatio-temporal dynamics in collective frog choruses examined by mathematical modeling and field observations,” Scientific Reports 4, 1–8 (2014).
  * Yan, Bae, and Granick (2015) J. Yan, S. C. Bae, and S. Granick, “Rotating crystals of magnetic janus colloids,” Soft Matter 11, 147–153 (2015).


11institutetext: Faculty of Mathematics and Computer Science, University of
Muenster, Germany 11email: {karim.huesmann, luis.garcia, linsen, b.risse}@uni-
muenster.de
∗ Corresponding Author (https://cvmls.uni-muenster.de)

# The Impact of Activation Sparsity on Overfitting in Convolutional Neural
Networks

Karim Huesmann  Luis Garcia Rodriguez  Lars Linsen  Benjamin Risse∗

###### Abstract

Overfitting is one of the fundamental challenges when training convolutional
neural networks and is usually identified by a diverging training and test
loss. The underlying dynamics of how the flow of activations induce
overfitting is however poorly understood. In this study we introduce a
perplexity-based sparsity definition to derive and visualise layer-wise
activation measures. These novel explainable AI strategies reveal a surprising
relationship between activation sparsity and overfitting, namely an increase
in sparsity in the feature extraction layers shortly before the test loss
starts rising. This tendency is preserved across network architectures and
reguralisation strategies so that our measures can be used as a reliable
indicator for overfitting while decoupling the network’s generalisation
capabilities from its loss-based definition. Moreover, our differentiable
sparsity formulation can be used to explicitly penalise the emergence of
sparsity during training so that the impact of reduced sparsity on overfitting
can be studied in real-time. Applying this penalty and analysing activation
sparsity for well known regularisers and in common network architectures
supports the hypothesis that reduced activation sparsity can effectively
improve the generalisation and classification performance. In line with other
recent work on this topic, our methods reveal novel insights into the
contradicting concepts of activation sparsity and network capacity by
demonstrating that dense activations can enable discriminative feature
learning while efficiently exploiting the capacity of deep models without
suffering from overfitting, even when trained excessively.

###### Keywords:

Explainable AI Sparstiy Overfitting Visualisation Technique CNNs

##  1 Introduction

In recent years deep convolutional neural networks (CNNs) achieved state-of-
the-art performances in most computer vision applications [21, 35]. The
ultimate goal of training neural networks networks is to achieve high
performance measures while avoiding the generalisation error, which is defined
by the difference between the training and test set loss [40]. This error
estimate requires that the trained model is independent of the test set
(independence hypothesis) [37]. If the independence hypothesis holds [37] and
generalisation error is high, then the model is suspected to have an
inappropriately high variance and is therefore overfitted to the training data
[41]. In fact, overfitting appears to be one of the fundamental challenges in
training deep CNNs since models often tend to learn too specific features of
the training set [14, 3].

###  1.1 Related Work

In the past, a variety of strategies have been proposed to prevent
overfitting, which can be roughly categorised into: (i) increasing the amount
of training data; (ii) reducing the models’ capacity; (iii) regularising the
model parameters; and (iv) early stopping procedures. Increasing the training
set can either be done by collecting additional data or by augmenting the
existing set [11, 17]. The reduction of the capacity is usually done by
explicitly pruning learnable parameters [16, 24, 22, 28]. In a similar
fashion, regularisation can also be used to decrease the capacity of the model
[29, 20]. For example [29] and [5] both aim to extenuate the model complexity
by using weight decay and it has also been shown that classical L1
regularisation tends to less complex models over time [38, 26]. In contrast,
decorrelation-based regularisers aim to employ the given capacity by reducing
redundancies based on hidden features [7, 4] or activations [6, 10], whereas
entropy-based regularisation uses the output distribution of the network to
reduce overfitting [27, 30]. Alternative regularisation strategies incorporate
additional layers. The two most common examples are dropout [34] and batch
normalisation [20]. Finally, early stopping has been studied for more than a
decade to end the training before the generalisation error gets too high [31]
and has recently been studied in the presence of label noise [23].

An important concept often related to overfitting is neural network sparsity,
which is interpreted as an indicator of success in learning discriminative
features [2] while filtering out irrelevant information, in analogy to
sparsity observed in the neocortex [36]. Consequently, many of the previously
mentioned pruning and regularisation techniques assume that sparse models
suffer less from generalisation error [16, 24, 22, 28, 29, 5, 38, 26]. In
addition, several empirical investigations have been done to demonstrate
beneficial effects of sparsity on the performance of deep neural networks [8,
13, 12]. However, recent experiments suggest that our comprehension of
sparsity is insufficient to fully understand its underlying relationship to
generalisation errors [25, 42]. Interestingly, even though sparsity inevitably
induces an underutilisation of the network’s capacity [4] it has never been
considered to be used to explain overfitting.

###  1.2 Contribution

Usually sparsity is referred to the property of zero-valued weights [13] and
is distinguished from activation sparsity which counts the number of zeros
after applying a non-linearity (in general ReLU) [32, 39]. Most existing work
on (activation) sparsity focuses on its beneficiary effects such as improved
inference performances or the robustness to adversarial attacks and noise [2,
18, 12]. Here we propose a novel activation sparsity definition which reveals
a surprising relationship between sparsity and overfitting. Our analysis
suggests that high activation sparsity is a reliable indicator for an
underutilisation of the potential network capacity and therefore limits the
generalisation capabilities of the model which induces overfitting. To study
this hypothesis we propose:

  1. 1.

a perplexity-based activation sparsity definition which uses network features
across neural receptive fields, yielding a novel layer-wise (i.e. targeted)
sparsity measures, called neural activation sparsity (NAS);

  2. 2.

targeted explainable AI (XAI) visualisation strategies to identify neural
activation sparsity in both, convolutional and fully connected layers in real-
time;

  3. 3.

a regularisation strategy derived from our neural activation sparsity
definition which can be used to penalise or reward neural activation sparsity
in a targeted manner.

Using these techniques, we demonstrate that an increase in neural activation
sparsity in all layers except the prediction layer indeed coincides with test
loss-based overfitting measures. Therefore, our XAI measures can be used as an
indicator for overfitting, which is independent of loss values, thus enables
the detection of generalisation errors even if the independence hypothesis of
the training and test data is violated. To further demonstrate this
relationship, we counteract neural activation sparsity in a layer-specific
fashion using a novel regularisation strategy: by penalising neural activation
sparsity in the feature extraction layers we can prevent the network from
overfitting while increasing the overall accuracy. We elaborate our findings
by performing several experiments including common network architectures
(including VGG16, ResNet50 and Xception) and regularisation strategies (L1,
L2, dropout and batch normalisation). Interestingly, targeted neural
activation sparsity regularisation prevents overfitting even in low training
data / high model size scenarios and outperforms conventional regularisation
techniques. In line with other recent work on activation sparsity [25, 42] we
demonstrate that neural network sparsity must not be a desirable property by
definition. Instead salient and discriminative feature learning while
exploiting the full capacity of the model appears to be another important
dimension to effectively train deep convolutional neural networks.

##  2 Methods

###  2.1 Activation Sparsity Definition

Neural network sparsity is usually defined by the amount of weights which are
exactly or close to zero [13]. Even though this quantification is inspired by
observations made in real brains [36], neuroscientific sparsity definitions
are based on neuronal activation patterns [33]. In a similar fashion sparse
information propagation in the forward path of artificial neural networks can
be the result of zero-valued weights or zero-valued input coefficients.
Motivated by this observation we introduce a novel sparsity indicator that is
purely based on the feature map activations which we refer to as neural
activation sparsity (NAS). We note that our sparsity definition differs from
existing activation sparsity definitions in two ways. Firstly, conventional
activation sparsity measures count the number of zeros after the use of the
non-linearity [32, 39], whereas our metric is calculated before the
nonlinearity is applied. As a consequence NAS is independent of the activation
function (i.e. less biased) and captures both, weight- and data-induced
sparsifications (i.e. more sensitive). Secondly, we measure relative neural
activity values in comparison to all other feature map activations at a given
neural receptive field. The neural receptive field is defined as the region of
a layer’s direct input that a filter is being affected by (see Figure 1).

Figure 1: NAS and NAS Heatmap generation. Left: sample input image. Middle:
schematic network architecture and sketch of an exemplary proximate receptive
field from the conv2 layer and the corresponding receptive field activation
vector in conv3 layer. Right: NAS Heatmap extracted by calculating the
activation sparsity on the depicted receptive field activation vector.

More formally, a layer with $D$ filters creates a $D\times M\times N$ shaped
feature tensor. Thus, the input of this layer consists of $M\cdot N=R$
receptive fields. Let $x_{d^{\prime},i,j,m,n}$ be a pixel at position $(i,j)$
and channel $d^{\prime}$ of receptive field $r_{m,n}$ . The corresponding
weight of filter $f_{d}$ which affects this pixel is given as
$w_{d^{\prime},i,j,d}$. The number of pixels in $r_{m,n}$ is equal to the
number of weights in $f_{d}$. Hence, the linear activation $a_{m,n,d}$ created
by $f_{d}$ when applied to $r_{m,n}$ is defined as

$a_{m,n,d}=\sum_{d^{\prime},i,j}w_{d^{\prime},i,j,d}\cdot
x_{d^{\prime},i,j,m,n}\ .$ (1)

To improve readability, we denote the results of this linear filtering as
receptive field activation vectors
$\textbf{a}_{k}\in\mathbb{R}^{D},k\in[0,R-1]$, where $k$ is a linear index
over all receptive fields computed by $k=m\cdot N+n$.

In order to define activation sparsity based on $\mathbf{a}_{k}$ we first
transform its components into a probability-like distribution by using the
Softmax function:

$p_{l}(\textbf{a}_{k})=\frac{e^{a_{k}^{l}}}{\sum_{m}e^{a_{k}^{m}}}\ ,\
l\in[1,...,D]$ (2)

where $a_{k}^{l}$ is the $l$-th component of $\textbf{a}_{k}$. Subsequently,
we calculate the perplexity for each normalised neural receptive field by

$\rho_{k}=e^{H_{k}}$ (3)

with

$H_{k}=H(\textbf{a}_{k})=-\sum_{l=1}^{D}p_{l}(\textbf{a}_{k})\
ln(p_{l}(\textbf{a}_{k}))\ .$ (4)

By definition, perplexity results in values in a range between $[1,D]$. Next,
we define a perplexity score $\tau$ by rescaling $\rho$ into values $\leq 1$:

$\tau_{k}=\frac{\rho(\textbf{p})_{k}}{D},$ (5)

where $\tau_{k}\in[\frac{1}{D},1]$. Based on this definition, the
smaller/larger $\tau_{k}$ is, the less/more filters are activated at the
neural receptive field $k$. Finally, we can define the neural activation
sparsity by

$s_{k}=1-\tau_{k}\ .$ (6)

This activation sparsity measure reaches its maximum value
$s_{k}=1-\frac{1}{D}$, if only one of the given $D$ feature maps returns an
activation $a_{l}$, $a_{l}\in\textbf{a}_{k}$, $p_{l}(\textbf{a}_{k})\neq 0$,
indicating lowest possible activation density for a given neural receptive
field. The minimum value $s_{k}=0$ can be reached, if all $D$ feature maps
return the same activation $a_{l}=a_{m},l\neq m$, which is equivalent to the
highest possible activation density and thus minimum activation sparsity.

###  2.2 Activation Sparsity Visualisation

Figure 2: Visualisations of layer NAS over time at three levels of aggregation
with corresponding test loss (red): (a) _NAS Heatmaps_ showing spatial
sparsity distributions for 3 layers. (b) _NAS Line Plot_ showing blue lines
(and bands) of layer NAS values. (c) _NAS Stripe Plot_ for 50 layers encoding
layer NAS in 1D colour coded stripes.

In order to enable a temporal visualisation of layer NAS for a variety of CNN
architectures we propose three different approaches that visually represent
the time series at different levels of aggregation using geometric (lines) and
colour (heatmap) encodings. These strategies are a trade-off between the level
of visualised details and the underlying network depth.

The most detailed visualisation can be achieved by spatially localised colour-
encodings of each NAS value for all layers. By applying the newly defined
activation sparsity, we can calculate $s_{k}$ for each position $k$ in all
layers, indicating how many of the layer’s filters produce a comparatively
high activation for the given filter input. Encoding and re-ordering all
$s_{k}$ into a 2D heatmap we can visualise the activation sparsity in a very
detailed manner for a given input as illustrated in Figure 2 (a): High
sparsity measures are given in yellow indicating that only a few filters
produced high activations at the location of the respective neural receptive
field. Likewise, the dark pixels indicate many highly activated filters (i.e.
low NAS measures) which naturally occurs at contrast-rich locations such as
edges and textured regions. Note that the same visualisation strategy can be
used for fully connected layers, with the only difference that the resulting
heatmap consists of a single value.

Even though NAS heatmaps enable spatial reasoning within the activation maps,
this visualisation can become overloaded given many layers or training epochs.
In order to provide a more concise overview of NAS values for all layers we
therefore introduce the NAS Line Plot which visualise the median NAS values
over all neural receptive fields for a given layer into the same illustration
as shown in Figure 2 (b). Each line therefore represents the change in
activation sparsity (y-axis) of an individual layer over time (x-axis), where
brightness of the colour encodes the depth of the layers (early layers in dark
blue and deep layers in bright blue). Furthermore, we are able to visualise
uncertainties between multiple trainings by including a band, which in fact
makes the plot a functional box plot over time (minimum, median, maximum). NAS
Line Plots allow for effectively reading off the NAS values for each layer at
each time step including their uncertainty, but they are limited to a certain
amount of layers since very deep networks will lead to occluded curves.

Therefore we also provide a visualisation technique which further abstracts
from the measured NAS values by means of a colour-coded 1D NAS heatmap over
time (x-axis), which we refer to as NAS Stripe Plot. The NAS Stripe Plots are
sorted by increasing layer depth (y-axis). As can be seen in Figure 2(c) this
visualisation encodes the least details but can be used to illustrate NAS
dynamics over dozens of layers.

In Figure 2 all NAS visualisations are superimposed with the respective test
loss over time (red line) to demonstrate the impact of activation sparsity on
overfitting. For NAS Heatmaps and Stripe Plots, we always applied the same
colour map as shown in Figure 2, encoding NAS values from low (dark blue,
$s_{k}=0$) to high (yellow, $s_{k}\approx 1$). Moreover these visualisations
can be calculated dynamically during training so that these XAI strategies can
offer deep insights into activation sparsity of a given network in real-time.

###  2.3 Activation Sparsity Regularisation

To test our hypothesis, that NAS is related to overfitting, we introduce a
penalty term that prevents the network layers from producing sparse
activations. This term is directly derived from our aforementioned activation
sparsity definition and utilises the fact that perplexity is differentiable in
our given ranges so that it can be used as an activity regulariser. More
formally, our regularisation penalty term $\mathcal{L}_{s}$ is defined by

$\mathcal{L}_{s}=-\sum_{i}\;\lambda_{i}\sum_{k=0}^{r_{i}}\rho_{k}^{i},$ (7)

where $i$ loops through all layers of the network, $r_{i}$ is the number of
receptive fields in the respective layer and $\rho_{k}^{i}$ refers to
$\rho_{k}$ as defined in Equation 3. In the following, this regularisation is
referred to as NASReg.

As perplexity reaches its maximum when all filters are activated in the same
way, this regulariser can have the tendency to produce highly correlated
filters. The trivial solution for $\mathcal{L}_{s}$ to be minimised is
therefore to generate identical filter responses. As identical filters reduce
the predictive power of neural networks, this effect has to be counterbalanced
by preventing high filter correlations. Therefore we have to set
$\lambda_{i}\geq 0$ in a way, that the perplexity is not reaching it’s
theoretical maximum.

###  2.4 Experimental Design

In our first set of experiments we used an intentionally simple base line
architecture called VanillaNet. This basic network comprises two conv-conv-
pool blocks followed by two fully connected layers and the layer sizes of
conv1, conv2, conv3, conv4, and fc1 are 256, 256, 512, 512, 1024 respectively.
ReLU is used as activation function across all layers except the last fully
connected layer.

Next, we analysed our XAI strategies in the context of different
regularisation techniques. In particular, we extended the VanillaNet by
including L1, L2, dropout and batch normalisation (called L1Net, L2Net,
DropNet and NormNet respectively) since these techniques are known to have an
impact on sparsity. For DropNet we set the dropout rate to 0.3 for layers
conv1 to fc1 and 0.5 between fc1 and prediction layer.

To study the generalisability of our metrics and visualisation strategies for
state-of-the-art architectures we evaluated three very deep CNNs, namely VGG16
[1], ResNet50 [19] and Xception [9]. Moreover, the behaviour of ResNet with
and without batch normalisation is analysed to further evaluate the impact of
regularisation on activation sparsity of deep nets.

For a better interpretability and to reduce side effects of the experiments we
have chosen a non-adaptive optimiser. Accordingly, we used vanilla Stochastic
Gradient Decent without momentum and a fixed learning rate of $0.01$ without
learning rate scheduling. This allows us to carry out a uniform training over
the entire course of the experiments. For the classification we are optimising
the categorical cross-entropy loss. The experiments are trained up to $300$
epochs with a batch size of $32$.

During the training, we calculate NAS for each convolution and fully connected
layer individually. For most of our experiments we chose the cifar-100 dataset
($50,000$ train and $10,000$ test images), because it is very susceptible to
overfitting. For one particular experiment we intentionally violated the
independence hypothesis of our test dataset by creating a new testing set
consisting of $5,000$ train and $5,000$ test images, which we mixed to
generate a distorted test loss. We conducted a hyperparameter search to find
appropriate values for $\lambda_{i}$ which resulted in a strong correlation
between $\lambda_{i}$ and the size of the layers. Therefore we reduced the
mean NAS via $\lambda_{i}=\frac{1}{r_{i}}$ for each layer $i$ in all
experiments. In order to add uncertainty to the results, we repeat each
training 3 times.

##  3 Results

###  3.1 Relationship Between Overfitting and Activation Sparsity

To test if classical test loss based generalisation errors coincide with our
NAS definition, we analysed the NAS Line Plots for multiple experiments. In
Figure 3(a) the resultant NAS Line Plots of VanillaNet on cifar-100 can be
seen, which are not affected by any regularisation, thus allowing a straight-
forward interpretation. Since cifar-100 is prone to overfitting if no
augmentation or advanced training protocol is used, it can be seen that the
test loss starts to rise quickly after epoch $4$. As it is apparent in the
plots, NAS also increases almost simultaneously with the test loss: Whilst the
NAS in conv1 rises relatively slowly, a rapid jump can already be observed in
conv2. The jump in NAS is even recognisable shortly before the test loss
diverges. The observation that conv1 NAS is usually low and that higher NAS
can be observed in deeper layers is in line with other work on activation
sparsity (i.e. based on counting zeros after ReLU activation) [32].

Since our NAS overfitting estimation is only based on activation patterns, it
does not depend on the satisfaction of the independence hypothesis [37]. To
illustrate this we created a test loss which violates the independence
hypothesis by including $50\%$ dependent images to this set (Figure 3(a);
orange line). As can be clearly seen, this perturbed test loss does not allow
generalisation error estimations, whereas our NAS indicators still diverge
accordingly.

To further investigate our observation that NAS coincides with overfitting, we
used dropout since this technique is known to reduce generalisation errors. As
shown in Figure 3(b) dropout indeed delays the divergence of the test loss and
also reduces its slope. However, overfitting cannot be avoided which is also
reflected in the median NAS measures throughout the layers. Moreover, all but
the conv1 NAS measures saturate close to $s_{k}=1$ after epoch $18$ indicating
a drastically reduced number of highly activated features propagating through
the network. This observation, combined with the constantly increasing test
loss value suggests, that the remaining features are over-adapted to the
training data.

As theory suggests, applying L1 regularisation encourages sparse weight
tensors and should therefore affect our NAS measures. As shown in Figure 3(c),
we can see a relatively steep increase of NAS in all layers compared to all
other experiments. Furthermore, all layers (even conv1) quickly converge
towards NAS values close to 1. Again, the moment of diverging loss and
increasing NAS coincide. A similar observation can be made for L2 regularised
networks (Figure 3(d)). Even though L2 regularisation is able to shift the
moment of overfitting to a later epoch, the corresponding NAS values start to
rise slowly before strong generalisation errors are observable in the test
loss estimates.

Adding batch normalisation layers to a given network architecture is a widely
used approach to prevent overfitting [20]. In line with other work on this
regularisation technique the loss does not show a noticeable divergence
throughout our training as can be seen in Figure 3(e). Similarly our NAS
estimates also tend towards constant values $<1$ which indicates that batch
normalisation has an attenuating effect on activation sparsity over time. Note
that the first layer (conv1) has no tendency towards sparsified activations
throughout the entire training. We conclude that despite the different effect
of batch normalisation on overfitting, the overall tendency of coinciding NAS
propagation can also be observed for almost all layers.

Figure 3: NAS Line Plots for Test loss (left axis) and median NAS (right axis)
over time. The dashed red lines indicate the moment of overfitting derived
from the test loss. All trainings are carried out on cifar-100 dataset. (a) -
(e) networks are trained for 30 epochs whereas the NASReg network in (f) is
trained for 300 epochs.

###  3.2 Spatial Analysis of Activation Sparstiy

Figure 4: NAS Heatmaps and corresponding test loss values for (a) VanillaNet,
(b) DropNet, (c) NormNet and (d) NASReg trainings. The respective heatmaps
were created after each corresponding epoch. The original input image can be
seen in Figure 1 (colour map corresponds to Figure 2).

Figure 4(a) illustrates the NAS Heatmaps generated for the individual
VanillaNet layers over $30$ epochs of training in comparison to the
corresponding test loss. Similar to the observations made in the previous
section the test loss indicates overfitting after epoch $4$ and the conv1
layer is not very susceptible to changes in NAS, whereas the other layers show
increased sparsity proportional to their depth. The same phenomenon can be
seen in Figure 4(b) in which we used dropout to delay overfitting. In analogy
to the previous results, the NAS measures shift accordingly to later epochs
but the overall tendency of coinciding generalisation error and increasing NAS
is preserved. Moreover, spatial NAS representations further reveal the
location of sparsified features and therefore less activated filters.

In order to further investigate our spatial NAS visualisation strategy, we
visualised these dynamics for a batch normalisation regularised network
(Figure 4(c)). Batch normalisation can effectively counteract overfitting as
can be seen in both, the loss plot and the corresponding NAS heatmaps. The
heatmaps reveal that the overall NAS is low and does not show a significant
change over time. However, subtle features of the original input appear
relatively early in deeper conv layers highlighting the strengths of this more
detailed visualisation technique.

In summary, these observations suggest that a progressive escalation of NAS
tends to coincide with overfitting. To further investigate both, this
hypothesis and the validity of our sparsity definition, we are penalising high
NAS values in a layer-wise fashion in the next experiment.

###  3.3 Activation Sparsity Regularisation Results

To study the impact of activation sparsity on overfitting in more detail, we
applied NASReg to enforce low NAS values in individual layers except the
prediction layer. Figure 3(f) shows the course of a NAS regularised training.
Note that in comparison to Figure 3(a) - (e) we substantially extended
training to $300$ epochs to show potential long-term effects. Since the
regularisation adds a penalty term to the loss, the resultant cross-entropy
loss is shown here. Juxtaposing NASReg with the other trainings, one can
clearly see the effect of the regularisation on the recorded training history.
VanillaNet, DropNet, L1Net and L2Net clearly start overfitting before the
$20$th epoch. In comparison, the NASReg test loss measures does not indicate
any generalisation errors even when trained for several hundreds of epochs.
Moreover this training reaches a significantly lower loss value compared to
the other trainings. One can observe that the NAS regularised layers reveal a
different behaviour in comparison to the other trainings: Despite the initial
increase in NAS estimates the regulariser pushes all NAS measurements towards
values close to zero shortly before epoch $15$.

The differences between NASReg and the other trainings is also shown in Figure
4(d). Here one can see a strong similarity to the NormNet training, which also
does not suffer from overfitting in the given epochs. Since we explicitly
encourage low NAS values using NASReg the heatmaps given in Figure 4(d)
comprise even lower activation sparsity values than NormNet. Moreover, we
regularise fc1 so that only the final prediction layer learns to produce
sparse activations, as desired for conventional classification tasks. Even
though we used intentionally straight forward models the positive influence of
regularised NAS can also be observed in the test accuracy: As can be seen in
Figure 3 the peak accuracy of NASReg ($0.5119$) outperforms VanillaNet
($0.3851$), DropNet ($0.4034$), L1Net ($0.3267$) and L2Net ($0.2866$). The
smaller accuracy difference to NormNet ($0.4947$) follows the intuition that a
certain degree of activation sparsity is acceptable in deeper layers.

###  3.4 Activation Sparsity in Common Deep Architectures

Figure 5: NAS Stripe Plots and corresponding test loss (a) VGG16, (b)
ResNet50-noBN, (c) ResNet50 and (d) Xception. Colour map equivalent to Figure
2.

After having examined the VanillaNet baseline architecture and its regularised
variations, we next investigate further NAS measurements in common deep
architectures, namely VGG16, ResNet50 and Xception. The training was performed
on cifar-100 using the same conditions as described in Section 2.4.

The VGG16 network mainly differs from VanillaNet by having a higher number of
convolutional layers and a more advanced pooling strategy. In line with our
hypothesis we can observe the same NAS behaviour when the network starts
overfitting (Figure 5(a)). Here, overfitting takes place at around epoch 30
and we can observe increased NAS measures especially in the middle layers
shortly before epoch 30.

The ResNet architecture is another common deep CNN, consisting of multiple
layers with several shortcut connections. Since batch normalisation has a
strong influence on activation sparsity (see Section 3.1) we investigated the
original regularised ResNet50 and an unregularised version (ResNet50-noBN).
The main difference between ResNet50-noBN and VGG16 is the higher number of
layers and the presence of shortcuts. Again, as soon as the training drifts
into overfitting, the NAS starts rising significantly in most layers (see
Figure 5(b)).

When applying batch normalisation to ResNet50, we observe an impact on both
the NAS and the overfitting measures (see Figure 5(c)). Similar to the results
presented for the VanillaNet (cf. Figure 3), continuous test loss increases
can be delayed while achieving better overall loss values when compared to the
non-regularised version. Moreover, batch norm again stabilises many layers to
constant NAS values throughout the training (cf. Section 3.1). However, the
loss rises sharply at epoch 5, 19 and 26, which is also reflected by a sudden
NAS increase for most layers. Overall a more chaotic NAS and test loss
behaviour can be observed until the network starts to overfit in epoch 30.

The overall relationship between overfitting and NAS can also be confirmed in
the Xception architecture (see Figure 5(d)): When compared to VGG16 and
ResNet50 the Xception network achieves best overall generalisation
performance, which is also reflected in the low activation sparsity measures.
Since these networks also utilise batch normalisation in most layers we can
also observe stabilised NAS values over time. Interestingly, even in this
complex architecture subtle overfitting tendency starting at about epoch 5 can
be observed in a deeper layer as slight increases in the NAS estimates are
visible at the moment of overfitting which again demonstrates the sensitivity
of our XAI measures.

##  4 Conclusion

In this paper we introduced a more sensitive and less biased neural activation
sparsity (NAS) definition which we used to derive reliable quantitative
measures and layer-wise visualisation techniques to study the dynamics in CNNs
during training. Interestingly, the analysis of NAS for various network types
and regularisation strategies revealed a correlation between activation
sparsity and overfitting. Throughout our experiments we observed a decreasing
amount of discriminative features in almost all layers over time so that the
influence of most neurons for the overall classification task became
negligible. This observation is in line with work on the beneficiary effects
of network pruning [15, 25]. Potential generalisation error tendencies of
heavily sparsified models were however not taken into consideration
previously. Our results suggest that high activation sparsity indeed coincides
with overfitting which we were able to detect on a per-layer basis. Moreover,
the rise in NAS values even occurs shortly before the test loss based
overfitting measures and our generalisation error estimates do not require the
independence between training and test data.

Using our visualisation strategy we were able to show the network performance
with respect to the utilised capacity during training. These visual
investigation gave additional evidence that low NAS can be a desirable
characteristic for the feature extraction layers. In contrast, the final
classification layer benefits from high activation sparsity to ensure distinct
classification decisions. Applying the NAS measures as a penalty term during
training we derived a novel regulariser which prevents CNNs from overfitting
over the entire course of training. In fact, our NAS-based regularisation
technique not just supports hypothetical relationship between overfitting and
activation sparsity, but also outperforms dropout and batch normalisation in
our experiments.

We studied the effects of NAS in a variety of different well-known
architectures and regularisation strategies and demonstrated that comparable
effects can be observed across all experiments. These experiments gave
additional evidence for the relation between sparsity and overfitting and
revealed interesting sparsity-related activation dynamics of batch
normalisation. In the future we will investigate these dynamics in more
detail, experiment with more complex network architectures, analyse NAS values
for adaptive optimisers such as ADAM and study the potential of our novel
regulariser for overall performance improvements.

#### Acknowledgements.

BR would like to thank the Ministeriums für Kultur und Wissenschaft des Landes
Nordrhein-Westfalen for the AI Starter support (ID 005-2010-005). Moreover,
the authors would like to sincerely thank Sören Klemm for his valuable ideas
and input throughout this project and the WWU IT for the usage of the PALMA2
supercomputer. This work was partially supported by the Deutsche
Forschungsgemeinschaft (DFG) under contract LI 1530/21-2.

## References

  * [1] (2020), C.C.C.S.C.: Very deep convolutional networks for large-scale image recognition. Computational and Biological Learning Society (2015)
  * [2] Ahmad, S., Scheinkman, L.: How can we be so dense? the benefits of using highly sparse representations. arXiv preprint arXiv:1903.11257 (2019)
  * [3] Ayinde, B.O., Inanc, T., Zurada, J.M.: On correlation of features extracted by deep neural networks. In: International Joint Conference on Neural Networks (IJCNN), Proc. pp. 1–8 (2019)
  * [4] Ayinde, B.O., Inanc, T., Zurada, J.M.: Regularizing deep neural networks by enhancing diversity in feature extraction. IEEE transactions on neural networks and learning systems (2019)
  * [5] Ayinde, B.O., Zurada, J.M.: Deep learning of constrained autoencoders for enhanced understanding of data. IEEE transactions on neural networks and learning systems 29(9), 3969–3979 (2017)
  * [6] Bao, Y., Jiang, H., Dai, L., Liu, C.: Incoherent training of deep neural networks to de-correlate bottleneck features for speech recognition. In: 2013 IEEE International Conference on Acoustics, Speech and Signal Processing. pp. 6980–6984. IEEE (2013)
  * [7] Bengio, Y., Bergstra, J.S.: Slow, decorrelated features for pretraining complex cell-like networks. In: Advances in neural information processing systems. pp. 99–107 (2009)
  * [8] Changpinyo, S., Sandler, M., Zhmoginov, A.: The power of sparsity in convolutional neural networks. arXiv preprint arXiv:1702.06257 (2017)
  * [9] Chollet, F.: Xception: Deep learning with depthwise separable convolutions. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 1251–1258 (2017)
  * [10] Cogswell, M., Ahmed, F., Girshick, R., Zitnick, L., Batra, D.: Reducing overfitting in deep networks by decorrelating representations. arXiv preprint arXiv:1511.06068 (2015)
  * [11] Cubuk, E.D., Zoph, B., Mane, D., Vasudevan, V., Le, Q.V.: Autoaugment: Learning augmentation strategies from data. In: IEEE Conference on Computer Vision and Pattern Recognition (CVPR), Proc. pp. 113–123 (2019)
  * [12] Frankle, J., Carbin, M.: The lottery ticket hypothesis: Finding sparse, trainable neural networks. In: International Conference on Learning Representations (2018)
  * [13] Gale, T., Elsen, E., Hooker, S.: The state of sparsity in deep neural networks. arXiv preprint arXiv:1902.09574 (2019)
  * [14] Gavrilov, A.D., Jordache, A., Vasdani, M., Deng, J.: Preventing model overfitting and underfitting in convolutional neural networks. International Journal of Software Science and Computational Intelligence (IJSSCI) 10(4), 19–28 (2018)
  * [15] Gomez, A.N., Zhang, I., Kamalakara, S.R., Madaan, D., Swersky, K., Gal, Y., Hinton, G.E.: Learning Sparse Networks Using Targeted Dropout. arXiv (2019)
  * [16] Goodfellow, I.J., Bengio, Y., Courville, A.C.: Deep Learning. Adaptive computation and machine learning, MIT Press (2016)
  * [17] Guo, H., Mao, Y., Zhang, R.: Mixup as locally linear out-of-manifold regularization. In: AAAI Conference on Artificial Intelligence (AAAI), Proc. pp. 3714–3722 (2019)
  * [18] Guo, Y., Zhang, C., Zhang, C., Chen, Y.: Sparse dnns with improved adversarial robustness. In: Advances in neural information processing systems. pp. 242–251 (2018)
  * [19] He, K., Zhang, X., Ren, S., Sun, J.: Deep residual learning for image recognition. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 770–778 (2016)
  * [20] Ioffe, S., Szegedy, C.: Batch normalization: Accelerating deep network training by reducing internal covariate shift. In: International Conference on Machine Learning (ICML), Proc. pp. 448–456 (2015)
  * [21] Khan, A., Sohail, A., Zahoora, U., Qureshi, A.S.: A survey of the recent architectures of deep convolutional neural networks. Artificial Intelligence Review 53(8), 5455–5516 (2020)
  * [22] Klemm, S., Ortkemper, R.D., Jiang, X.: Deploying deep learning into practice: A case study on fundus segmentation. In: Conference in Medical Imaging, Understanding and Analysis (MIUA), Proc. (2019), [in print]
  * [23] Li, M., Soltanolkotabi, M., Oymak, S.: Gradient descent with early stopping is provably robust to label noise for overparameterized neural networks. In: International Conference on Artificial Intelligence and Statistics. pp. 4313–4324. PMLR (2020)
  * [24] Liu, Z., Li, J., Shen, Z., Huang, G., Yan, S., Zhang, C.: Learning efficient convolutional networks through network slimming. In: Proceedings of the IEEE International Conference on Computer Vision. pp. 2736–2744 (2017)
  * [25] Liu, Z., Sun, M., Zhou, T., Huang, G., Darrell, T.: Rethinking the value of network pruning. In: International Conference on Learning Representations (2018)
  * [26] Mehta, D., Kim, K.I., Theobalt, C.: On implicit filter level sparsity in convolutional neural networks. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition. pp. 520–528 (2019)
  * [27] Miller, D.J., Rao, A.V., Rose, K., Gersho, A.: A global optimization technique for statistical classifier design. IEEE Transactions on Signal Processing 44(12), 3108–3122 (1996)
  * [28] Molchanov, P., Tyree, S., Karras, T., Aila, T., Kautz, J.: Pruning convolutional neural networks for resource efficient inference. In: 5th International Conference on Learning Representations, ICLR 2017-Conference Track Proceedings (2019)
  * [29] Nowlan, S.J., Hinton, G.E.: Simplifying neural networks by soft weight-sharing. Neural computation 4(4), 473–493 (1992)
  * [30] Pereyra, G., Tucker, G., Chorowski, J., Kaiser, L., Hinton, G.E.: Regularizing neural networks by penalizing confident output distributions. International Conference on Learning Representations (ICLR), Proc. (2017)
  * [31] Prechelt, L.: Early stopping-but when? In: Neural Networks: Tricks of the trade, pp. 55–69. Springer (1998)
  * [32] Rhu, M., O’Connor, M., Chatterjee, N., Pool, J., Kwon, Y., Keckler, S.W.: Compressing DMA Engine: Leveraging Activation Sparsity for Training Deep Neural Networks. 2018 IEEE International Symposium on High Performance Computer Architecture (HPCA) pp. 78–91 (2018)
  * [33] Seelig, J.D., Chiappe, M.E., Lott, G.K., Dutta, A., Osborne, J.E., Reiser, M.B., Jayaraman, V.: Two-photon calcium imaging from head-fixed drosophila during optomotor walking behavior. Nature methods 7(7), 535 (2010)
  * [34] Srivastava, N., Hinton, G.E., Krizhevsky, A., Sutskever, I., Salakhutdinov, R.: Dropout: a simple way to prevent neural networks from overfitting. Jornal of Machine Learning Research (JMLR) 15(1), 1929–1958 (2014)
  * [35] Tu, Z., Xie, W., Zhang, D., Poppe, R., Veltkamp, R.C., Li, B., Yuan, J.: A survey of variational and cnn-based optical flow techniques. Signal Processing: Image Communication 72, 9–24 (2019)
  * [36] Vinje, W.E., Gallant, J.L.: Sparse coding and decorrelation in primary visual cortex during natural vision. Science 287(5456), 1273–1276 (2000)
  * [37] Werpachowski, R., György, A., Szepesvári, C.: Detecting overfitting via adversarial examples. In: Advances in Neural Information Processing Systems. pp. 7856–7866 (2019)
  * [38] Yaguchi, A., Suzuki, T., Asano, W., Nitta, S., Sakata, Y., Tanizawa, A.: Adam induces implicit weight sparsity in rectifier neural networks. In: IEEE International Conference on Machine Learning and Applications (ICMLA), Proc. pp. 318–325 (2018)
  * [39] Yang, Q., Mao, J., Wang, Z., Li, H.: Dasnet: Dynamic activation sparsity for neural network efficiency improvement. In: 2019 IEEE 31st International Conference on Tools with Artificial Intelligence (ICTAI). pp. 1401–1405. IEEE (2019)
  * [40] Zhang, C., Bengio, S., Hardt, M., Recht, B., Vinyals, O.: Understanding deep learning requires rethinking generalization. In: International Conference on Learning Representations (ICLR), Proc. (2017)
  * [41] Zhang, C., Vinyals, O., Munos, R., Bengio, S.: A study on overfitting in deep reinforcement learning. arXiv preprint arXiv:1804.06893 (2018)
  * [42] Zhou, H., Lan, J., Liu, R., Yosinski, J.: Deconstructing lottery tickets: Zeros, signs, and the supermask. In: Advances in Neural Information Processing Systems. pp. 3592–3602 (2019)