# ALIFE2022 template

Eleni Nisioti1 and Clément Moulin-Frier1
1 Flowers Team, Inria and Ensta ParisTech, Bordeaux, France
<EMAIL_ADDRESS>
# Dynamics of niche construction
in adaptable populations evolving in diverse environments

Eleni Nisioti1 and Clément Moulin-Frier1
1 Flowers Team, Inria and Ensta ParisTech, Bordeaux, France
<EMAIL_ADDRESS>
###### Abstract

In both natural and artificial studies, evolution is often seen as synonymous
to natural selection. Individuals evolve under pressures set by environments
that are either reset or do not carry over significant changes from previous
generations. Thus, niche construction (NC), the reciprocal process to natural
selection where individuals incur inheritable changes to their environment, is
ignored. Arguably due to this lack of study, the dynamics of NC are today
little understood, especially in real-world settings. In this work, we study
NC in simulation environments that consist of multiple, diverse niches and
populations that evolve their plasticity, evolvability and niche-constructing
behaviors. Our empirical analysis reveals many interesting dynamics, with
populations experiencing mass extinctions, arms races and oscillations 111We
provide an online repo for reproducing our simulations at
https://github.com/eleninisioti/NicheConstructionModel/tree/main. To
understand these behaviors, we analyze the interaction between NC and
adaptability and the effect of NC on the population’s genomic diversity and
dispersal, observing that NC diversifies niches. Our study suggests that
complexifying the simulation environments studying NC, by considering multiple
and diverse niches, is necessary for understanding its dynamics and can lend
testable hypotheses to future studies of both natural and artificial systems.

## Introduction

Biological organisms and their environments share a reciprocal relationship:
organisms survive and reproduce under selection pressures present in their
habitats and environments are modified by their inhabitants, with changes
being inherited to the next generation and accumulating with evolutionary time
(Odling-Smee et al.,, 2003). The first process, natural selection, became the
cornerstone of evolutionary theory in the early 20th century (Provine,, 2001)
and, in both natural and artificial life studies, is often seen as synonymous
to evolution. The second, niche construction (NC), was characterized as the
neglected process in evolution in the early 2000s (Odling-Smee et al.,, 2003),
as evolutionary theory assumed that NC’s effect on selection pressures is
negligible. Since then, studies of natural populations have shown that NC
often affects selection pressures by helping species protect themselves from
environmental uncertainty and accelerates evolution by complexifying
environments (Odling-Smee et al.,, 2003; Schwilk,, 2003; Boivin et al.,,
2016). By now, NC is, at best, the elephant in the room of evolutionary
synthesis: despite evidence for its existence, our understanding of it is too
limited to enable its study in settings capturing the complexity of the real
world.

Recent hypotheses studying major events in evolution point to environmental
complexity (Maslin et al.,, 2015; Potts,, 2013; d’Errico et al.,, 2017). For
example, the birth of our own lineage in East Africa occurred amidst large
climatic instability that fragmented the landscape into patches of land
differing in resource availability and separated by large lakes (Maslin et
al.,, 2014). Other hypotheses further emphasize the ability of natural
landmarks such as deserts and oceans to form barriers that isolate populations
(Larrasoaña et al.,, 2013). A similar story is unfolding in artificial life:
under Quality-Diversity optimization (Pugh et al.,, 2016), the search space is
divided into behavioral niches to evolve diverse solutions, a paradigm that
challenged the dominant approach of training a single agent, as it proved more
robust in real-world settings (Cully et al.,, 2015).

When in a heterogeneous environment, organisms can survive by: a) specializing
in a niche to out-compete others, paying the cost that they are out-competed
in other niches Nisioti and Moulin-Frier, (2022); Grove, (2014) b) becoming
phenotypically plastic: a plastic individual adapts its phenotype to its
environment without genetic change and survive in diverse niches within its
lifetime. On the downside, plasticity comes with fitness costs, so that
plastic individuals are out-competed by specialists in their preferred niche
(Grove,, 2014; Nisioti and Moulin-Frier,, 2022) c) becoming highly evolvable:
a high mutation rate enables quick adaptation at an evolutionary scale but
comes at the cost of increased deleterious mutations (Lynch et al.,, 2016;
Nisioti and Moulin-Frier,, 2022). As both plasticity and evolvability enable
adaptation we jointly refer to them as adaptability.

Although limited, our understanding of NC shows that it is influenced both by
the adaptability of populations and the diversity of environments. For
example, processes such as the invention of agriculture (Zohari,, 1986) and
cultural innovation (Mesoudi and Thornton,, 2022) are both influenced by
geography and our impressive social learning abilities (Diamond,, 1998;
Migliano et al.,, 2020; Migliano and Vinicius,, 2022). Explanations for this
rely on the inter-play between a population’s connectivity and its ability to
accumulate solutions: spatial heterogeneity enforces isolation and, therefore,
diversity within a population, while adaptability ensures that environmental
barriers can be crossed and solutions spread, causing an intensification of NC
(Cantor and Whitehead,, 2013; Derex and Boyd,, 2016; Nisioti et al.,, 2022).

In this work, we study the evolution of NC in environments divided in diverse
niches where populations evolve their adaptability and niche-constructing
behavior. Each niche is characterized by its environmental state, which
determines how many agents it can fit, and consists of two components: the
intrinsic state (which can for example model climatic variations
uncontrollable by the agents) and the niche-constructed state, which is the
product of niche-constructing individuals inhabiting the niche. The niche-
constructed state is also inherited from the previous generation with some
decaying applied to capture the fact that ecologically-inherited artifacts
cannot persist indefinitely.  Modeling NC as a process that changes the
capacity of the environment is common (Laland et al.,, 1999; Krakauer et al.,,
2009) and finds inspiration in natural behaviors, such as nest-constructing
species. Genomes consist of four genes: preferred environmental state (the
state at which their fitness is highest), plasticity (the variance in
environmental state they can tolerate without a large impact on fitness),
evolvability (the mutation rate) and niche construction (the amount by which
they change the environmental state of a niche they inhabit). We consider two
different mechanisms for selecting which agents will reproduce: a) under
global competition we select them with a probability proportional to their
average fitness across the niches they can survive in and reproduction stops
when the environment’s capacity is reached b) under local competition we
consider reproduction in each niche independently: an agent is selected based
on its fitness in each niche and reproduction stops when the capacity of that
niche is filled. Thus, agents benefit from surviving in multiple niches.  Our
models of local and global competition have their equivalents both in previous
studies in evolutionary computation (with Quality-Diversity algorithms
representing local competition and classical algorithms global) and in the
study of natural populations (where local adaptation considers that
populations of the same species inhabiting different niches experience
different natural selection pressures (Leimu and Fischer,, 2008)).

An agent may benefit from NC by: a) increasing the capacity of a niche, for
example through nest-building, thus, reducing competition for resources b)
reducing the capacity of a niche thus making it less desirable to others.
This behavior, which we refer to as negative niche construction, is
exemplified in nature by pines, which spread their needles to increase the
probability of fire and out-compete not fire-resistant species (Schwilk,,
2003) c) bringinf the environmental state closer to its preferred niche and
further away from others  For example, aquatic earthworms can inhabit earth
soil only by changing its consistency (Turner,, 2000). But NC also comes at a
cost: a) increasing the capacity of a niche may invite other species in and
increase competition b) decreasing its capacity may make it uninhabitable c)
changing the environmental state will require adaptability, which comes at its
own cost d) multiple agents co-inhabiting a niche can create high
environmental uncertainty that can lead to mass extinctions if adaptability
cannot increase enough to respond to it, a case that feels most relevant to
our own evolution (Boivin et al.,, 2016; West,, 2017).

Taking into account these complex interactions between environmental
heterogeneity, population adaptability and niche construction, we believe that
studies of NC should depart from simplified settings with a single niche
(Laland et al.,, 1999; Suzuki and Arita,, 2009; Chiba et al.,, 2020) and non-
plastic populations (Laland et al.,, 1999; Krakauer et al.,, 2009). As we show
in this work, rich models can offer interesting hypotheses about the dynamics
of NC, such as that:

  * •

having multiple niches is necessary for avoiding mass extinctions. We observed
that, in a single niche, environmental uncertainty induced by NC becomes too
high for adaptability to cope with;

  * •

NC promotes adaptability and populations adapt differently depending on the
selection mechanism: under global competition NC leads to higher plasticity
while under local to higher evolvability;

  * •

the population may niche-construct negatively for prolonged periods of time,
making niches uninhabitable. Agents deal with this by evacuating the niche
until it becomes inhabitable due to the decay of NC;

  * •

under local competition NC increases genomic and environmental diversity

## Related works

Previous studies have primarily modeled NC in two different ways: as an
increase in the amount of resources that increases the environment’s capacity
(Laland et al.,, 1999; Krakauer et al.,, 2009) and as a direct increase or
decrease in fitness that does not affect capacity (Suzuki and Arita,, 2009).
We follow the first approach but also allow reductions in capacity. NC studies
were at first primarily theoretical, employing differential equations to
predict how NC affects evolution (Laland et al.,, 1999; Krakauer et al.,,
2009). Such models lent intuitive insights, such that NC can emerge without
direct selection (Laland et al.,, 1999) and that the ability to monopolize
niches enables the emergence of NC (Krakauer et al.,, 2009) but come with
limitations, such as the fact that niches are identical and agents are not
plastic.

Suzuki and Arita, (2009) study the co-evolution of plasticity and NC in an
agent-based model and show that populations alternate between phases of high
plasticity and NC, provided that agents niche-construct sequentially. This
model differs from ours as the environment consists of a single niche, NC does
not modify its capacity and evolvability is constant. Interestingly, when
Suzuki and Arita, (2009) applied NC in parallel these patterns disappeared due
to agents canceling out each other’s behavior. Instead, the patterns in our
analysis occur with NC applied in parallel. As we show, this becomes possible
due to the presence of multiple niches that stabilize NC. Nisioti and Moulin-
Frier, (2022) studie the co-evolution of plasticity and evolvability in an
environment with multiple niches, but do not consider NC.

Recent advents in deep reinforcement learning (DRL) and neuroevolution have
enabled the study of complex behaviors in simulations, such as foraging and
tool-use (Pérolat et al.,, 2017; Hamon et al.,, 2023; Baker et al.,, 2020).
Chiba et al., (2020) leverage such techniques to study NC in a 2D environment
where agents can construct artifacts useful for avoiding predation. Our study
can lend insights to this direction, as we can view DRL as the mechanism that
underlies the behavioral plasticity considered in our model.

## Modeling and methodology

We now separately discuss our model of the environment and genomes, the
evolutionary algorithm and the set of metrics we use to monitor evolution.

Figure 1: Our model of the environment: each niche $n$ is characterized by its
environmental state $e_{n}^{g}$ which is the sum of the intrinsic state
$i_{n}$ (in blue) and the niche-constructed state $b_{n}^{g}$ (in red).
$i_{n}$ depends on the state of the reference niche $(e_{0})$, the offset
$\epsilon$ and $n$. $b_{n}^{g}$ is the sum of the niche-constructed state at
the previous generation discounted by $\gamma$ plus the sum of the niche-
constructing behavior of agents that reproduced in this niche in the current
generation.

### Modeling the environment

The environment, illustrated in Figure 1, is divided into $N$ niches arranged
in a simple latitudinal model: we consider a reference niche at $n=0$, $N/2$
“northern” niches indicated with positive indexes $n\in(0,N/2]$ and $N/2-1$
“southern” niches with negative indexes $n\in(-N/2,0)$ . Each niche is
characterized by its environmental state $e_{n}^{g}$ which is the sum of two
elements:

  * •

an intrinsic state $i_{n}$ that remains constant with evolutionary time and
depends on the location of the niche. Specifically, the state of niche $n$ at
generation $g$ is $i_{n}=i_{0}+\epsilon\cdot n$, where $i_{0}$ is the state of
the reference niche and $\epsilon$ is a constant capturing the difference
between adjacent niches.

  * •

the niche-constructed state $b_{n}^{g}$, capturing the modifications that
agents inhabiting the niche cause. These modifications are carried over
generations but are discounted by a factor $\gamma$. Formally, the niche-
constructed state $b_{n}^{g}$ of niche $n$ is equal to
$b_{n}^{g-1}\cdot\gamma$, with $\gamma<1$, plus the total amount of NC applied
by all agents that reproduced in this niche in generation $g$. We denote this
latter amount as $\sum_{k\in\mathcal{K}_{n}}a_{k}^{g}$, where
$\mathcal{K}_{n}$ is the subset of all agents that reproduced in niche $n$ and
$a_{k}^{g}$ is the amount of niche construction by a single agent (we later
explain how agents niche-construct).

Thus, the general equation describing the evolution of niche $n$ with
generations $g$ is:

$\displaystyle b_{n}^{g}$
$\displaystyle=b_{n}^{g-1}\cdot\gamma+\sum_{k\in\mathcal{K}_{n}}a_{k}^{g}$
$\displaystyle e_{n}^{g}$ $\displaystyle=i_{0}+\epsilon\cdot n+b_{n}^{g}$ (1)

As we explain later, the environmental state determines the fitness of an
agent based on its genome and the capacity of the niche $c_{n}^{g}$ as:
$c_{n}^{g}=e_{n}^{g}C_{N}$, where $C_{N}$ is the reference niche capacity.
Thus, higher environmental states can support larger populations and are
termed “high-quality”. To ensure that the maximum population size is
independent of the number of niches we define $C_{N}=C_{\text{ref}}/N$, where
$C_{\text{ref}}$ is equal to the desirable maximum population size. An
assumption of this model is that there is spatial smoothness, i.e, nearby
niches are similar. Note that $C_{\text{ref}}$ is independent of niche-
constructing behavior: by niche-constructing the population can exceed this
capacity. We, therefore, further bound the population size by randomly
discarding agents if the population exceeds a value $K_{\text{max}}$.

### Modeling the genome

We model plasticity through tolerance curves, a tool developed in ecology
(Lynch and Gabriel,, 1987) and previously employed in simulation environments
(Grove,, 2014; Nisioti and Moulin-Frier,, 2022). A tolerance curve is a normal
distribution with mean $\mu_{k}^{g}$, indicating the environmental state of
highest fitness for an individual $k$ at generation $g$, called the preferred
state, and standard deviation $\sigma_{k}^{g}$ that captures how quickly the
fitness of the genome drops as the environmental state varies from its
preferred state. Genomes with large $\sigma_{k}^{g}$ are indicative of plastic
individuals (we illustrate the tolerance curves of a plastic and a non-plastic
individual in Figure 2).

Figure 2: Modeling plasticity as a normal distribution
$\mathcal{N}(\mu_{k},\sigma_{k})$. A non-plastic individual ($k$) has small
$\sigma_{k}$ and a high peak at their preferred niche, while a plastic
individual ($k^{{}^{\prime}}$) has large $\sigma_{k}$ and a lower peak at
their preferred niche. Fitness in a given niche $n$ is computed as the
probability density function of the distribution at the environmental state
$e_{n}$. This figure also illustrates the cost and benefit of plasticity,
assuming that $\mu_{k}=\mu_{k}^{{}^{\prime}}$. If $e_{n}=\mu_{k}$ (the actual
environmental state is identical to the preferred niche of both individuals)
the plastic individual has lower fitness (cost of plasticity). If $e_{n}$
differs significantly from $\mu_{k}$ (the actual environmental state differs
from the preferred one) the plastic individual has higher fitness (benefit of
plasticity).

To model niche construction we introduce a niche-constructing gene
$a_{k}\in[-N/K_{\text{max}},N/K_{\text{max}}]$ where we bound the amount an
agent can niche-construct to ensure that agents inhabiting large environments
have limited capabilities. The niche-constructing gene is expressed in the
niche an agent reproduces in by modifying its environmental state. Thus, at
each generation, niche $n$ is constructed by an amount
$\sum_{k\in\mathcal{K}_{n}}a_{k}$, where $\mathcal{K}_{n}$ denotes the subset
of agents that reproduced in it. The genome $o_{k}^{g}$ also includes the
mutation rate $r_{k}^{g}$. Thus the complete form of a genome is
$o_{k}^{g}=[\mu_{k}^{g},\sigma_{k}^{g},r_{k}^{g},a_{k}^{g}]$ and, upon
reproduction, it mutates as:

$\displaystyle\mu_{k}^{g+1}=\mu_{k}^{g}+\mathcal{N}(0,r_{k}^{g})$
$\displaystyle\sigma_{k}^{g+1}=\sigma_{k}^{g}+\mathcal{N}(0,r_{k}^{g})$
$\displaystyle
a_{k}^{g+1}=[c_{k}^{g}+\mathcal{N}(0,r_{a})]_{-N/K_{\text{max}}}^{N/K_{\text{max}}}$
$\displaystyle r_{k}^{g+1}=r_{k}^{g}+\mathcal{N}(0,r_{k}^{g})$ (2)

where $\mathcal{N}(x,y)$ denotes a normal distribution with mean $x$ and
variance $y$ and we have highlighted that the niche-constructing gene is
bounded. Note also that, for stability reasons, we employ a different mutation
rate for the niche-constructing gene, $r_{a}$, that remains constant.

### Global and local selection

At the end of a generation agents are selected for crossover based on their
fitness. Each chosen individual has two offspring and the next generation
consists only of offspring. To compute the fitness of an individual $k$ in
generation $g$ we first detect the niches in which it can survive as:

$\displaystyle{n\in\\{1,\cdots,N\\}\quad|\quad
e_{n}^{g}\in[\mu_{k}^{g}-2\sigma_{k}^{g},\mu_{k}^{g}+2\sigma_{k}^{g}]}$ (3)

and compute its fitness in each one of them as
$f_{k,n}^{g}=pdf(\mu_{k}^{g},\sigma_{k}^{g},e_{n}^{g})$, where $pdf$ denotes
the value of the normal probability density function with mean $\mu_{k}^{g}$,
and variance $\sigma_{k}^{g}$ at location $e_{n}^{g}$. We study two selection
mechanisms:

  * •

under global selection all agents in the population are ranked based on their
average fitness across the niches they can survive in and reproduce with a
probability proportional to it.  Agents reproduce until the capacity of the
environment is filled.

  * •

under local selection we apply the same criterion but independently for each
niche: we detect which agents survive in a niche and reproduce them with a
probability proportional to their fitness. Thus, agents that can survive in
multiple niches have higher chances of reproduction under this mechanism.
Again, agents reproduce in a niche until its capacity is filled.

Thus, at the end of a generation a new population is formed that consists of
the offspring of agents from the previous population that were chosen for
reproduction, based on their fitness and the environment’s capacity. Note that
the niche is not inherited from parent to offspring: the offspring will
inhabit the niches its genome describes (as in Eq. 3). We present the
pseudocode of our algorithm in Algorithm 1 and 2, provided in our online repo
due to limited space.

### Metrics

In addition to population-wide averages of the genome values and environment-
wide averages of the intrinsic and niche-constructed states we monitor the
following metrics:

  1. 1.

$X^{g}=\sum_{k}X_{k}^{g}$, the number of extinctions. We denote the survival
of individual $k$ in niche $n$ at generation $g$ as a binary variable:

$\displaystyle
s_{k,n}^{g}=(e_{n,g}\in[\mu_{k}^{g}-2\sigma_{k}^{g},\mu_{k}^{g}+2\sigma_{k}^{g}])$
(4)

Thus, an individual goes extinct ($X_{k}^{g}=1$) if $\sum_{n}^{N}s_{k,n}^{g}$
is zero and survives ($X_{k}^{g}=0$) if $\sum_{n}^{N}s_{k,n}^{g}$ is positive.

  2. 2.

$V^{g}_{\mu}$, the diversity of the population defined as the standard
deviation of the population’s preferred state, namely

$\displaystyle V^{g}_{\mu}=\sigma_{\mu^{g}}$ (5)

This metric captures the genetic diversity of the population computed for the
gene of preferred state.

  3. 3.

$D^{g}$, the dispersal of the population, computed as the number of niches
over which at least one individual survives for a temporal window of at least
$w$ generations. Formally, $D^{g}=\sum_{n=1}^{N}d_{n,w}^{g}$, where
$d_{n,w}^{g}$ denotes the persistence of the population in a given niche for
the required time window and is computed as

$\displaystyle
d_{n,w}^{g}=\begin{cases}1\quad\text{if}\sum_{g^{\prime}=g-w}^{g}s_{n}^{g^{\prime}}=w\\\
\text{0}\quad\text{otherwise}\end{cases}$ (6)

where $s_{n}^{g}$ is indicates the survival of at least one individual in a
given niche and is computed as

$\displaystyle s_{n}^{g}=\begin{cases}1\quad{if}\sum_{k}^{K}s_{k,n}^{g}>1\\\
0\quad\text{otherwise}\end{cases}$ (7)

with $s_{k,n}^{g}$ defined in Eq. (4).

  4. 4.

$H^{g}$, the competition for reproduction within the population. We count, for
each niche, the number of agents that survived in it and were fit enough to be
chosen for reproduction but did not reproduce because its capacity was reached
and, then, sum over all niches.

## Results

We now examine the behavior of agents under different settings. First, we
compare behaviors under global competition between environments with different
number of niches. After, we consider only heterogeneous environments ($N=100$)
niches and we contrast the behavior across two dimensions: a) whether NC takes
place or not (where we force all niche-constructing genes to be zero) and b)
whether selection is local or global. In all simulations we set the discount
factor $\gamma$ to 0.5, the intrinsic state of the reference niche $i_{0}$ to
0.6, $\epsilon$ (the difference between adjacent niches) to 0.01, the
reference capacity to $C_{\text{ref}}=1000$, the maximum population size to
$K_{\text{max}}=5000$ and the mutation rate of the niche-constructing gene,
$r_{c}$, to 0.0003.  As we discuss later it would be interesting to study the
effect of some of these hyper-parameters. For this study we note that the
values of $C_{\text{ref}},K_{\text{max}}$ where chosen to limit computational
complexity, increasing them further should not qualitatively change our
conclusions. Also regarding $r_{c}$, setting it to too low a value will
disable NC and setting it too high may destabilize NC. We performed ten
independent trials and present in plots median values and $95\%$ confidence
intervals. When studying niche-constructing behaviors we do not average across
trials but present them individually, as differences would average out and
conceal information.

### Spatial heterogeneity stabilizes niche construction

In Figure 3 we compare the behavior under global competition for different
number of niches ($N\in\\{1,50,100\\}$ for a single trial. We observe that,
for $N=1$, the agents are niche-constructing positively (third row), which
pushes the environment state to higher values and more variability (first
row). We see that the the population reacts by adapting its preferred niche
(second row), increasing its plasticity (fourth row) and keeping its
evolvability high (fifth row). Despite that, the population goes extinct
around generation 250, after experiencing many booms and busts (sixth row).
This behavior was consistent across all trials for $N=1$, with half of the
trials experiencing a mass extinction due to negative NC and half due to
positive NC (as in Figure 3).

Figure 3: Comparison between a homogeneous ($N=1$) and a heterogeneous
environment ($N=100$) under global competition

In the heterogeneous environments, on the other hand, the population survives
for the whole simulation. We see that the average NC among agents
($\mathbb{E}(a)$) is similar to that of the homogeneous environment but that
this leads to less niche construction in the environment (evident through the
environmental state). This is because the population in the heterogeneous
environment spreads over multiple niches and often migrates collectively out
of a niche (we will take a look at this behavior in the next section), leading
to a decrease to the accumulation of NC due to its discounting. In contrast,
when $N=1$ NC accumulates in a single niche. We also observe that NC in
$N=100$ does not increase continuously but experiences oscillations, which
suggests that some agents are niche-constructing negatively or reducing their
positive NC. As a result, the environment experiences lower oscillations which
enables the population to reduce its evolvability. Plasticity remains
relatively high but much lower than in the homogeneous environment.  We also
measured the number of mass extinctions for environments with different number
of niches and observed that even a small number of niches is useful and
increasing the number of niches continuously decreases the probability of
extinction (10/10 went extinct for $N=1$, 4/10 for 20, 3/10 for 50 and 2/10
for 100). Thus, the presence of multiple niches was necessary for the survival
of the niche-constructing population. Next, we will look into how adaptability
interacts with NC to give rise to this behavior. As we hypothesize that the
presence of niches matters, we will focus on how agents disperse and niche-
construct differently in different niches.

### Niche construction promotes adaptability

Figure 4: Comparison between populations with ($R_{NC}$) and without NC ($R$)
under global competition.

We now compare a niche-constructing population (denoted as $R_{NC}$) and a
population that cannot niche -construct ($R$). We examine populations under
global selection in Figure 4 and under local selection in Figure 5. We observe
that both populations increase their adaptability but do so differently.

Figure 5: Comparison between populations with ($R_{NC}$) and without NC ($R$)
under local competition.

Under global selection, the $R_{NC}$ population increases its plasticity but
keeps its evolvability low. Dispersal is a bit higher and genomic diversity is
initially high but, by the middle of the simulation, reaches the low level of
the $R$ population. We also see that $R_{NC}$ has a smaller population size
early in the simulation. As extinctions are low, this is due to NC reducing
the capacity of the environment. Near the end of the simulation we see that
$R_{NC}$ reaches and, in some trials, surpasses the size of the $R$
population. This indicates that positive NC has increased the environment’s
capacity.

Under local selection $R_{NC}$ has very high evolvability but plasticity is at
similar level to $R$. We observe that some extinctions persist throughout
evolution but the population size is relatively stable.  Why did the niche-
constructing population prefer to adapt through evolvability rather than
plasticity? As we will see more closely in our analysis of specific trials,
agents in a niche-constructing population manage to coordinate their NC
behavior within a niche. This would not have been possible under high
plasticity, as agents would stochastically niche-construct in different niches
and increase environmental uncertainty.

### Niche construction can cause an arms race

Figure 6: Analyzing one trial for global competition with niche-constructing
population where an arms race emerged (Top) evolution of metrics (Bottom)
Heatmap with rows corresponding to niches, columns to generations and
value/color of cell to the sum of niche construction of all agents in a given
niche and generation.

In Figure 6 we analyze one of the trials for the niche-constructing population
under global competition through heatmaps that show the amount of NC at each
generation and niche. Populations are initialized with diverse genomes, so
that at the beginning they are dispersed in all niches. Then, global
competition leads to a decrease in plasticity and quickly wipes out the
diversity of the population, gathering all agents in the same narrow stripe of
niches. Around generation 50, the population starts positively niche-
constructing in the middle niche, then starts moving to the north and, once it
reaches the end of the environment (around generation 150), stays there and
oscillates between phases of positive and negative NC. On the top plot, we see
that this population has adapted its preferred state to the last niche and
that as soon it reaches it, plasticity stops decreasing and the population
size stops increasing. The trial depicted in Figure 6 is representative of 4
out of the 10 trials. In another five trials the same pattern occured, but in
the south instead of the north. In contrast, we observed that populations
without NC do not gather at the edge of the environment but stay in a random
niche near the center. Thus, NC under global competition leads to a
”geographic” arms race due to the following dynamics: once the agents have
been gathered to a narrow stripe of niches, they begin niche-constructing,
initially randomly. Then, NC will, by random chance become slightly positive
or negative. As all agents are competing in the same niche, this will force
them to adapt their preferred state towards the direction that the sum of
agents is constructing to. This enables inhabiting the adjacent niche, where
the agents will niche-construct towards the same direction. Eventually, the
population will reach the end of the environment. Why doesn’t the population
continue increasing its niche construction to the same direction instead of
experiencing cycles? We believe that this is not possible as the mutation rate
is low and, therefore, the preferred niche cannot be adapted anymore.

### Negative niche construction reduces competition

Figure 7: Analyzing one trial for global competition with niche-constructing
population where negative niche construction emerged.

In Figure 7 we analyze another trial for the niche-constructing population
under global competition, where a different behavior emerged: the population
niche-constructs negatively, which pushes it to the low-index niches. There we
see an interesting pattern: the population is gathered in a narrow stripe of
niches where it niche-constructs negatively and with increasing intensity.
Then the population leaves this stripe and moves to the adjacent northern and
southern niches. This is possible due to the maintained plasticity and does
not require an adaptation of the preferred niche. Then, the same behavior
happens until the population moves back to the previous stripe and we see many
cycles of moving back and forth. This switching behavior is caused by the fact
that, once the population niche constructs negatively in a niche for some
time, then this niche becomes uninhabitable, so that the population needs to
move out until it becomes inhabitable again. Eventually, some positive NC
happens, which leads to an extinction, as evolvability is too low to enable
adaptation. This indicates that negative niche construction is not stable in
the long-term.

### Niche construction diversifies the environment

We now move to local selection, where we saw in Figure 5 that NC leads to a
more evolvable population. In Figure 8, we analyze a typical trial where we
see that NC is positive and high in northern niches and negative in more
southern niches. Thus, the majority lives in the north (this can be inferred
from the preferred state $\mathbb{E}(\mu)=0.8$ surpassing the state of the
reference niche $i_{0}=0.6$). The population experiences many changes, with
niches often switching between being positively and negatively niche-
constructed. We also measured the variance of NC within a niche and found that
it is two orders of magnitude smaller than the one under global selection.
This suggests that local selection enables agents to coordinate their niche
construction within a niche, which diversifies the environment and leads to
the higher genomic diversity observed in Figure 5.

Figure 8: Analyzing a typical trial for local competition with a niche-
constructing population.

An intriguing question is why populations under global selection adapt by
increasing their evolvability but keep their plasticity at even lower levels
than without NC. We believe that this is due to the interplay between NC and
plasticity: very plastic agents could niche-construct in many niches, which
would make it difficult to maintain low variance in NC within a niche. By
monopolizing a limited number of niches and quickly changing niche through
mutations these populations can better coordinate their NC.

## Discussion

We have shown that the evolution of niche construction is contingent on the
number of niches in the environment and the selection mechanism. We explained
this by looking into the interplay between adaptability and NC, in particular
their effect on the population’s and environment’s diversity and dispersal
patterns. Populations are impressive in their ability to self-regulate
themselves: NC remains bounded, even though we do not introduce an explicit
cost for it. Populations fail to coordinate their NC and go extinct only when
the environment has a single niche.

Our empirical study could be extended in multiple ways. First, we could study
the effect of additional hyper-parameters. For example, we could allow the
intrinsic state to vary with time, following a noisy or periodic signal
(Nisioti and Moulin-Frier,, 2022). We could also consider that NC is plastic:
currently an agent’s gene determines its NC independently of the environmental
state while we could imagine scenarios where agents niche-construct
differently for different states. Finally, we believe that our empirical
observations should be tested in grounded environments with RL agents, such as
grid-worlds where a population forages or avoids predation (Hamon et al.,,
2023; Chiba et al.,, 2020).

We hope that our work will contribute towards incorporating NC in future
studies of collective adaptation. We think this is particularly relevant for
the Artificial Intelligence community, which recently focused on agents that
can generalize to diverse environments, the importance of environmental
complexity and multi-agent dynamics (Reed et al.,, 2022; Jaderberg et al.,,
2022; Nisioti et al.,, 2021; Moulin-Frier,, 2022). As we show here, niche
construction, which can be seen as a meta-learning mechanism (Constant et
al.,, 2018), can prove promising in complexifying environments and population
dynamics and bring us closer to artificial agents with behaviors reminiscent
of natural ones.

#### Acknowledgements

This research was partially funded by the French National Research Agency
(https://anr.fr/, project ECOCURL, Grant ANR-20-CE23-0006). This work also
benefited from access to the HPC resources of IDRIS under the allocation
2020-[A0091011996] made by GENCI.

## References

  * Baker et al., (2020) Baker, B., Kanitscheider, I., Markov, T., Wu, Y., Powell, G., McGrew, B., and Mordatch, I. (2020).  Emergent Tool Use From Multi-Agent Autocurricula.  In International Conference on Learning Representations.
  * Boivin et al., (2016) Boivin, N. L., Zeder, M. A., Fuller, D. Q., Crowther, A., Larson, G., Erlandson, J. M., Denham, T., and Petraglia, M. D. (2016).  Ecological consequences of human niche construction: Examining long-term anthropogenic shaping of global species distributions.  Proceedings of the National Academy of Sciences, 113(23):6388–6396.
  * Cantor and Whitehead, (2013) Cantor, M. and Whitehead, H. (2013).  The interplay between social networks and culture: Theoretically and among whales and dolphins.  Philosophical Transactions of the Royal Society B: Biological Sciences, 368(1618):20120340.
  * Chiba et al., (2020) Chiba, N., Suzuki, R., and Arita, T. (2020).  Evolution of Complex Niche-Constructing Behaviors and Ecological Inheritance of Adaptive Structures in a Physically Grounded Environment.  Frontiers in Robotics and AI, 7:45.
  * Constant et al., (2018) Constant, A., Ramstead, M. J. D., Veissière, S. P. L., Campbell, J. O., and Friston, K. J. (2018).  A variational approach to niche construction.  Journal of The Royal Society Interface, 15(141):20170685.
  * Cully et al., (2015) Cully, A., Clune, J., Tarapore, D., and Mouret, J.-B. (2015).  Robots that can adapt like animals.  Nature, 521(7553):503–507.
  * Derex and Boyd, (2016) Derex, M. and Boyd, R. (2016).  Partial connectivity increases cultural accumulation within groups.  Proceedings of the National Academy of Sciences, 113(11):2982–2987.
  * d’Errico et al., (2017) d’Errico, F., Banks, W. E., Warren, D. L., Sgubin, G., van Niekerk, K., Henshilwood, C., Daniau, A.-L., and Sánchez Goñi, M. F. (2017).  Identifying early modern human ecological niche expansions and associated cultural dynamics in the South African Middle Stone Age.  Proceedings of the National Academy of Sciences, 114(30):7869–7876.
  * Diamond, (1998) Diamond, J. M. (1998).  Guns, Germs, and Steel: the Fates of Human Societies.  W. W. Norton & Co., New York.
  * Grove, (2014) Grove, M. (2014).  Evolution and dispersal under climatic instability: A simple evolutionary algorithm.  Adaptive Behavior, 22(4):235–254.
  * Hamon et al., (2023) Hamon, G., Nisioti, E., and Moulin-Frier, C. (2023).  Eco-evolutionary Dynamics of Non-episodic Neuroevolution in Large Multi-agent Environments.
  * Jaderberg et al., (2022) Jaderberg, M., Mathieu, M., McAleese, N., Bradley-Schmieg, N., Wong, N., Porcel, N., Hughes-Fitt, S., Dalibard, V., and Czarnecki, W. M. (2022).  Open-Ended learning leads to Generally Capable Agents.  page 54.
  * Krakauer et al., (2009) Krakauer, D. C., Page, K. M., and Erwin, D. H. (2009).  Diversity, Dilemmas, and Monopolies of Niche Construction.  The American Naturalist, 173(1):26–40.
  * Laland et al., (1999) Laland, K. N., Odling-Smee, F. J., and Feldman, M. W. (1999).  Evolutionary consequences of niche construction and their implications for ecology.  Proceedings of the National Academy of Sciences, 96(18):10242–10247.
  * Larrasoaña et al., (2013) Larrasoaña, J. C., Roberts, A. P., and Rohling, E. J. (2013).  Dynamics of Green Sahara Periods and Their Role in Hominin Evolution.  PLOS ONE, 8(10):e76514.
  * Leimu and Fischer, (2008) Leimu, R. and Fischer, M. (2008).  A Meta-Analysis of Local Adaptation in Plants.  PLoS ONE, 3(12):e4010.
  * Lynch et al., (2016) Lynch, M., Ackerman, M. S., Gout, J.-F., Long, H., Sung, W., Thomas, W. K., and Foster, P. L. (2016).  Genetic drift, selection and the evolution of the mutation rate.  Nature Reviews Genetics, 17(11):704–714.
  * Lynch and Gabriel, (1987) Lynch, M. and Gabriel, W. (1987).  Environmental tolerance.  The American Naturalist, 129(2):283–303.
  * Maslin et al., (2014) Maslin, M. A., Brierley, C. M., Milner, A. M., Shultz, S., Trauth, M. H., and Wilson, K. E. (2014).  East African climate pulses and early human evolution.  Quaternary Science Reviews, 101:1–17.
  * Maslin et al., (2015) Maslin, M. A., Shultz, S., and Trauth, M. H. (2015).  A synthesis of the theories and concepts of early human evolution.  Philosophical Transactions of the Royal Society B: Biological Sciences, 370(1663):20140064.
  * Mesoudi and Thornton, (2022) Mesoudi, A. and Thornton, A. (2022).  What is cumulative cultural evolution?  Proceedings of the Royal Society B: Biological Sciences, 285(1880):20180712.
  * Migliano et al., (2020) Migliano, A. B., Battiston, F., Viguier, S., Page, A. E., Dyble, M., Schlaepfer, R., Smith, D., Astete, L., Ngales, M., Gomez-Gardenes, J., Latora, V., and Vinicius, L. (2020).  Hunter-gatherer multilevel sociality accelerates cumulative cultural evolution.  Science Advances, 6(9):eaax5913.
  * Migliano and Vinicius, (2022) Migliano, A. B. and Vinicius, L. (2022).  The origins of human cumulative culture: From the foraging niche to collective intelligence.  Philosophical Transactions of the Royal Society B: Biological Sciences, 377(1843):20200317.
  * Moulin-Frier, (2022) Moulin-Frier, C. (2022).  The Ecology of Open-Ended Skill Acquisition.  Habilitation à diriger des recherches, Université de Bordeaux (UB).
  * Nisioti et al., (2021) Nisioti, E., Jodogne-del Litto, K., and Moulin-Frier, C. (2021).  Grounding an Ecological Theory of Artificial Intelligence in Human Evolution.  In NeurIPS 2021 - Conference on Neural Information Processing Systems / Workshop: Ecological Theory of Reinforcement Learning, virtual event, France.
  * Nisioti et al., (2022) Nisioti, E., Mahaut, M., Oudeyer, P.-Y., Momennejad, I., and Moulin-Frier, C. (2022).  Social Network Structure Shapes Innovation: Experience-sharing in RL with SAPIENS.
  * Nisioti and Moulin-Frier, (2022) Nisioti, E. and Moulin-Frier, C. (2022).  Plasticity and evolvability under environmental variability: The joint role of fitness-based selection and niche-limited competition.  In GECCO 2022 - The Genetic and Evolutionary Computation Conference, Boston / Hybrid, United States.
  * Odling-Smee et al., (2003) Odling-Smee, F. J., Laland, K. N., and Feldman, M. W. (2003).  Niche Construction: The Neglected Process in Evolution (MPB-37).  Princeton University Press.
  * Pérolat et al., (2017) Pérolat, J., Leibo, J. Z., Zambaldi, V., Beattie, C., Tuyls, K., and Graepel, T. (2017).  A multi-agent reinforcement learning model of common-pool resource appropriation.  In Guyon, I., Luxburg, U. V., Bengio, S., Wallach, H., Fergus, R., Vishwanathan, S., and Garnett, R., editors, Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc.
  * Potts, (2013) Potts, R. (2013).  Hominin evolution in settings of strong environmental variability.  Quaternary Science Reviews, 73:1–13.
  * Provine, (2001) Provine, W. B. (2001).  The Origins of Theoretical Population Genetics: With a New Afterword.  University of Chicago Press, Chicago, IL.
  * Pugh et al., (2016) Pugh, J. K., Soros, L. B., and Stanley, K. O. (2016).  Quality Diversity: A New Frontier for Evolutionary Computation.  Frontiers in Robotics and AI, 3.
  * Reed et al., (2022) Reed, S., Zolna, K., Parisotto, E., Colmenarejo, S. G., Novikov, A., Barth-Maron, G., Gimenez, M., Sulsky, Y., Kay, J., Springenberg, J. T., Eccles, T., Bruce, J., Razavi, A., Edwards, A., Heess, N., Chen, Y., Hadsell, R., Vinyals, O., Bordbar, M., and de Freitas, N. (2022).  A generalist agent.
  * Schwilk, (2003) Schwilk, D. W. (2003).  Flammability Is a Niche Construction Trait: Canopy Architecture Affects Fire Intensity.  The American Naturalist, 162(6):725–733.
  * Suzuki and Arita, (2009) Suzuki, R. and Arita, T. (2009).  Effects of Temporal Locality of Ecological Processes on Coevolution of Learning and Niche Construction.
  * Turner, (2000) Turner, J. S. (2000).  The extended organism.
  * West, (2017) West, G. (2017).  Scale.
  * Zohari, (1986) Zohari, D. (1986).  The Origin and Early Spread of Agriculture in the Old World.  In Barigozzi, C., editor, Developments in Agricultural and Managed Forest Ecology, volume 16 of The Origin and Domestication of Cultivated Plants, pages 3–20. Elsevier.


# Robust quantum many-body scars in the one-dimensional spin-$1$ Kitaev model

Sashikanta Mohapatra<EMAIL_ADDRESS>Ajit C. Balram<EMAIL_ADDRESS>Institute of Mathematical Sciences, CIT Campus, Chennai 600113, India  Homi
Bhabha National Institute, Training School Complex, Anushaktinagar, Mumbai
400094, India

###### Abstract

Experimental observation of coherent oscillations in a Rydberg atom chain
[Bernien _et al._ , Nature 551, 579 (2017)] has led to the discovery of
quantum many-body scars (QMBS) which is a new paradigm for ergodicity-
breaking. The experimental findings in the Rydberg chain can be well captured
by a kinetically constrained model called the “PXP” model, which has been
shown to host the Eigenstate Thermalization Hypothesis (ETH)-violating scar
states in the middle of the spectrum. Much effort has been put into
identifying similar kinetically restricted systems that show a violation of
ETH. In this work, we study the QMBS that can arise in one such model, namely
the spin-$1$ Kitaev chain, where owing to some conserved quantities, the
Hilbert space gets fragmented into unequal disconnected subspaces. Recently,
You _et. al_ [Phys. Rev. Research 4, 013103 (2022)] showed that the ground
state sector of this chain can be mapped exactly onto the prototypical PXP
model and thus hosts QMBSs. Here, we demonstrate that the phenomenon of
scarring is also present in other sectors, and in particular, we identify a
sector that exhibits substantially more scarring than the ground state one. We
propose an initial state and numerically demonstrate that its fidelity
revivals are robust and longer-lived than those in the PXP model.

##  I Introduction

Rapid improvements in the platforms for realizing and controlling non-
equilibrium dynamics of closed quantum systems, such as ultracold atoms [1],
trapped ions [2], nitrogen-vacancy centers [3], etc., has enabled a study of
the thermalization of quantum systems isolated from external baths. A generic
isolated quantum system is expected to be ergodic, i.e., under the unitary
dynamics of its Hamiltonian, any initial state would eventually evolve into a
featureless thermal state. This loss of information on the initial state’s
configuration presents a barrier to protecting quantum information. As a
result, it is crucial to search for non-ergodic systems that resist
thermalization. The Eigenstate Thermalization Hypothesis (ETH) [4, 5]
regulates the characteristics of ergodic quantum systems and describes how
far-from-equilibrium initial states evolve in time to reach a final state that
is described by a thermal ensemble. ETH suggests _all_ the eigenstates of
ergodic systems are thermal and thus any initial state evolves into a thermal
state at long times. Two well-known exceptions to the ETH paradigm are
integrable and many-body localized systems [6, 7, 8]. In integrable systems,
the presence of an extensive number of conserved quantities prevents an
initial state from fully exploring all the allowed configurations in the
Hilbert space. In MBL systems, the presence of interactions [9] and strong
disorder [10] leads to an emergent integrability that prevents thermalization.
These two ergodicity-breaking mechanisms are of the strong form in that
_every_ eigenstate exhibits features of an athermal state.

Recently, experimental findings in an ultracold Rydberg atom chain [11]
revealed a new mechanism for weak ergodicity-breaking. When the Rydberg atoms
were initialized in a particular state, the so-called Néel state, they do not
thermalize and instead display long-lived coherent oscillations. On the other
hand, certain other initial states do exhibit thermal behavior. The
theoretical description of the Rydberg chain is captured by the kinetically
constrained “PXP” model [12]. Since the Rydberg atoms are quite large, it is
energetically prohibitive to simultaneously excite two nearest neighboring
atoms [13, 14]. The ‘P’ in the PXP is a projector that exactly projects out
these configurations in which the nearest neighboring sites are both in
excited states. This Rydberg blockade constraint imposes a restriction on the
allowed configurations for the system which results in a constrained Hilbert
space that the system can access. Numerical studies of the PXP model [15, 16,
17, 18] have revealed the presence of anomalous states at equidistant energies
that have sub-extensive entanglement entropy (EE) in the otherwise thermal
bulk spectrum. These special eigenstates obey the area-law of EE rather than
the volume-law of EE as anticipated by ETH and have substantial overlap with
the Néel state which results in the observed coherent revivals. This
phenomenon is dubbed quantum many-body scars (QMBS) [15]. These scar states
are vanishingly rare and typically their number grows only algebraically with
system size while the Hilbert space dimension grows exponentially with system
size. As a result, these scars only lead to a weak or incomplete breach of
ETH.

In recent years, substantial theoretical effort [19, 20, 21, 22, 23, 24] has
been put in, in tandem with experiments [22, 25], to look for systems that can
support QMBS. QMBS have been identified in many-body systems, such as the
Affleck-Kennedy-Lieb-Tasaki model [26, 27, 28], integer spin XY model [29,
30], $\eta$-pairing Hubbard model [31, 32, 33], thin torus limit of quantum
hall phases [34], tilted 1D Fermi-Hubbard model [22], etc. In this work, we
study QMBS in the spin-$1$ Kitaev chain [35, 36], where previous studies [37]
have demonstrated that the PXP model is embedded in one of its subspaces,
thereby making it an ideal candidate system to support QMBS. We study other
subspaces (besides the one that has the PXP in it) of this model and see if
they too can support QMBS. We identify a sector where the scarring is
considerably stronger than that observed in the PXP model. Analogous to the
Néel state of the PXP model, we propose an initial state in this sector that
shows remarkably persistent fidelity oscillations.

The remainder of this paper is organized as follows. We give a brief overview
of the one-dimensional (1D) spin-$1$ Kitaev model in Sec. II. In Sec. III.1 we
study a particular sector of this model and its associated constrained
dynamics and find that this subspace hosts anomalous scarred states. We
identify an initial state in this subspace and numerically demonstrate that it
has robust and long-lived coherent oscillations. We show that the forward
scattering approximation nicely captures these scarred states. In Sec. III.2
we consider some other subspaces of the Kitaev chain and show that the
fidelity oscillations of analogous initial states in these subspaces decay
rapidly. Finally, we summarize our results in Sec. IV and present an outlook
for the future.

##  II The one-dimensional Kitaev model

The spin-$1$ Kitaev chain can be obtained as a single row of the two-
dimensional Kitaev model [35]. We start with the general spin-$S$ Kitaev model
on the honeycomb lattice that is described by the Hamiltonian

$H_{K}^{\rm 2D}=J_{x}\sum_{{\langle
i,j\rangle}_{x}}S_{i}^{x}S_{j}^{x}+J_{y}\sum_{{\langle
i,j\rangle}_{y}}S_{i}^{y}S_{j}^{y}+J_{z}\sum_{{\langle
i,j\rangle}_{z}}S_{i}^{z}S_{j}^{z},$ (1)

where operators $S_{j}^{a}$ (with $a{=}x,y,z$) are the spin-$S$ operators at
site $j$ and $\langle i,j\rangle_{a}$ denotes nearest neighbors in the
$a$-direction. The spin operators satisfy the usual $SU(2)$ algebra i.e.,
$[S_{i}^{a},S_{j}^{b}]{=}i\delta_{ij}\epsilon_{abc}S_{j}^{c}$, where
$\epsilon_{abc}$ is the totally anti-symmetric Levi-Civita tensor. Setting
$J_{z}{=}0$ in Eq. (1), we get a set of decoupled 1D chains any one of which
of $N$ sites is described by the Hamiltonian [36]

$H^{\rm
1D}_{K}(\\{J\\})=\sum_{j=1}^{N/2}(J_{2j-1}S_{2j-1}^{x}S_{2j}^{x}+J_{2j}S_{2j}^{y}S_{2j+1}^{y}).$
(2)

In general, the coupling constants $J$’s could be different from each other.
However, throughout this work, we will consider the simple case where all
$J$’s are equal and set to unit strength i.e., $J_{l}=1~{}\forall l$. Thus, we
end up with the following Hamiltonian for the spin-$S$ Kitaev chain

$H^{\rm
1D}_{K}=\sum_{j=1}^{N/2}(S_{2j-1}^{x}S_{2j}^{x}+S_{2j}^{y}S_{2j+1}^{y}),$ (3)

which is the model that we will work with throughout this paper. Next, we
would like to find the symmetries of the Hamiltonian of Eq. (3). To do so, we
define site parity operators $\mathcal{P}_{j}^{a}$ on every site as

$\mathcal{P}^{a}=e^{i\pi S_{j}^{a}}.$ (4)

The Ising-like terms in Eq. (3) change the value of total $S_{z}$ at the sites
adjoining a link i.e., $S^{z}_{2j{-}1}{+}S^{z}_{2j}$ at the $x$-link
($2j{-}1,2j$) and the value of $S^{z}_{2j}{+}S^{z}_{2j{+}1}$ at the $y$-link
($2j,2j{+}1$), by either $0$ or ${\pm}2$. Therefore the bond parity operators
$\mathcal{B}_{j}$ on odd and even bonds defined by

$\mathcal{B}_{2j-1}=\mathcal{P}_{2j-1}^{y}\mathcal{P}_{2j}^{y},\hskip
5.69046pt\text{and}\hskip
5.69046pt\mathcal{B}_{2j}=\mathcal{P}_{2j}^{x}\mathcal{P}_{2j+1}^{x}$ (5)

remain invariant under the action of Hamiltonian. Thus we have

$[\mathcal{B}_{j},H]=0,~{}\forall j,$ (6)

and these constitute symmetries of the spin-$S$ Kitaev chain. By performing
the following unitary transformation on the even sites [36]

$S_{2j}^{x}\rightarrow S_{2j}^{y},\hskip 14.22636ptS_{2j}^{y}\rightarrow
S_{2j}^{x}\hskip 14.22636ptS_{2j}^{z}\rightarrow-S_{2j}^{z},$ (7)

the Hamiltonian can be cast into the following convenient translationally
invariant form

$H=\sum_{j=1}^{N}S_{j}^{x}S_{j+1}^{y}.$ (8)

Upon the unitary transformation of Eq. (7), the bond parity operators take the
universal form (independent of whether bond $j$ is even or odd)

$\mathcal{B}_{j}=\mathcal{P}_{j}^{y}\mathcal{P}_{j+1}^{x}.$ (9)

From here on, we shall restrict ourselves to the spin-$1$ case of our interest
and work with its natural representation given by the orthonormal basis states
$\\{|x\rangle,|y\rangle,|z\rangle\\}$ defined as

$\begin{split}\ket{x}\equiv\frac{1}{\sqrt{2}}(\ket{-1}-\ket{1}),\ket{y}\equiv\frac{i}{\sqrt{2}}(\ket{-1}+\ket{1}),\ket{z}\equiv\ket{0},\end{split}$
(10)

where $|m\rangle$ is the eigenstate of the spin-$1$ operator $S_{i}^{z}$ with
eigenvalue $m{=}{-}1,0,1$. In this representation the spin-$1$ operators can
be written as $S_{bc}^{a}{=}i\epsilon_{abc}$ and their matrix representation
is

$\begin{split}S^{x}=\begin{pmatrix}0&0&0\\\ 0&0&-i\\\
0&i&0\end{pmatrix},S^{y}=\begin{pmatrix}0&0&i\\\ 0&0&0\\\
-i&0&0\end{pmatrix},S^{z}=\begin{pmatrix}0&-i&0\\\ i&0&0\\\
0&0&0\end{pmatrix}.\end{split}$ (11)

Furthermore, the $3\times 3$ matrices corresponding to the site parity
operators $\mathcal{P}^{a}$ of Eq. (4) are diagonal and given by

$\begin{split}\mathcal{P}^{x}{=}\begin{pmatrix}1&0&0\\\ 0&{-}1&0\\\
0&0&{-}1\end{pmatrix},\mathcal{P}^{y}{=}\begin{pmatrix}{-}1&0&0\\\ 0&1&0\\\
0&0&{-}1\end{pmatrix},\mathcal{P}^{z}{=}\begin{pmatrix}{-}1&0&0\\\ 0&{-}1&0\\\
0&0&1\end{pmatrix}.\end{split}$ (12)

From this matrix representation, we can readily read off that the eigenvalues
of the operators $\mathcal{P}_{j}^{a}$ are ${\pm}1$ with the eigenvalue ${-}1$
being doubly degenerate. Therefore, the eigenvalues of bond parity operators
$\mathcal{B}_{j}$ defined in Eq. (9) are also $b_{j}{=}{\pm}1$ since they are
just products of the site parity operators. Moreover, as the site-parity
operators $\mathcal{P}_{j}^{a}$ are diagonal, they commute with each other.
The bond operators $\mathcal{B}_{j}$ being products of diagonal site-parity
operators are also diagonal and commute with each other [along with the fact
that they commute with the Hamiltonian as shown in Eq. (6)]. This implies the
Hilbert space can be decomposed into $2^{N}$ sectors (of unequal sizes since
the eigenvalue ${-}1$ is doubly degenerate) and each sector can be represented
by a set of bond invariants $\vec{b}{\equiv}\\{b_{1},b_{2},{\cdots},b_{N}\\}$.

Projection into these sectors imposes restrictions on the allowed
configurations of two neighboring sites. For the nearest neighbor sites
$\langle j,j{+}1\rangle$ there are a total of $3{\times}3{=}9$ allowed states
which, based on the eigenvalue of the bond operator $\mathcal{B}_{j}$, get
fragmented into the following two sets

$\ket{xy},\hskip 5.69046pt\ket{xz},\hskip 5.69046pt\ket{yx},\hskip
5.69046pt\ket{zy}\hskip 5.69046pt\text{and}\hskip 5.69046pt\ket{zz}\hskip
5.69046pt\text{have}\hskip 5.69046ptb_{j}{=}1,$ (13)

and

$|{xx}\rangle,\hskip 5.69046pt|{yy}\rangle,\hskip 5.69046pt|{yz}\rangle\hskip
5.69046pt\text{and}\hskip 5.69046pt|{zx}\rangle\hskip
5.69046pt\text{have}\hskip 5.69046ptb_{j}{=}{-}1.$ (14)

The existence of these constrained subspaces makes the spin-$1$ Kitaev chain a
viable candidate to host QMBS.

##  III QMBS in the spin-$1$ Kitaev chain

The authors of Ref. [36] showed that the ground state of the Hamiltonian of
Eq. (8) lies in the subspace with all $b_{j}{=}1$. The restriction on the
neighboring sites in this sector exactly mimics the Rydberg blockade
constraint [36, 37]. Thus, the $\vec{b}{=}\\{1,1,\cdots,1\\}$ subspace can
exactly be mapped into the PXP model (see App. A) and therefore hosts QMBS
[37]. The corresponding Néel state for the spin-$1$ chain is given by
$|Z_{2}\rangle_{\rm Kitaev}{=}|yx\rangle{\equiv}|yxyx\cdots yx\rangle$ and the
fidelity for this state
$F(t){=}|\langle{Z_{2}}|\exp(-iHt)|{Z_{2}}\rangle|^{2}$ gives rise to coherent
oscillation as shown in Fig. 1.

Figure 1: Fidelity of the $|Z_{2}\rangle_{\rm Kitaev}$ state showing periodic
revivals. Data are for $N{=}22$ sites with periodic boundary conditions. The
dimension of the sector $D(1,1,\cdots,1){=}39,603$.

We will show in the subsequent sections that some other subspaces of the
spin-$1$ Kitaev chain also harbor scarred eigenstates, though in general, it
is difficult to find the corresponding spin-$1/2$ Hamiltonian like the PXP one
as it involves complicated forms with long-range interactions. In particular,
we find that the $\vec{b}{=}\\{1,1,{-}1,1,1,{-}1,\cdots,1,1,{-}1\\}$ sector
exhibits a more pronounced scarring effect than the ground state one and we
will discuss the fate of QMBS in this sector next.

###  III.1 The $\vec{b}{=}\\{1,1,{-}1,1,1-1,{\cdots},1,1,{-}1\\}$ sector

We first unravel the structure of the constrained Hilbert space of this
sector. There are two types of constraints on the states of nearest
neighboring sites: i) since $b_{3j}{=}{-}1$, there are four possible states
given in Eq. (14) that the neighboring sites $\langle 3j,3j{+}1\rangle$ can be
in, whereas, ii) since $b_{k}{=}1~{}\forall k{\neq 3j}$, there are five
possible states given in Eq. (13) that the neighboring sites $\langle
k,k{+}1\rangle$, $k{\neq 3j}$ can be in. The dimension of Hilbert space
$\mathcal{H}$ of this sector is known to be
$D(1,1,{-}1,1,1,{-}1,{\cdots},1,1,{-}1){\approx}1.55113^{N}$ for a system of
size $N$ with periodic boundary condition (PBC) [36]. In this sector, the
Hamiltonian transforms the state $|{yz}\rangle$ to $|{zx}\rangle$ and vice-
versa over the bond with $b_{j}{=}{-}1$ and $|{zz}\rangle$ to $|{yx}\rangle$
and vice-versa over the bond with $b_{j}{=}1$. Therefore the Hilbert space of
this sector can be constructed by taking any initial state of the sector as
root state (call it $|{R}\rangle)$ and successively applying the Hamiltonian
on it, i.e.,

$\mathcal{H}_{\\{1,1,-1,1,1-1,\cdots,1,1-1\\}}\equiv\text{Span}\\{|{R}\rangle,H|{R}\rangle,H^{2}|{R}\rangle,\cdots\\}.$
(15)

Figure 2: The action of the Kitaev Hamiltonian of Eq. (8) on the Hilbert space
of the sector $\vec{b}{=}\\{1,1,-1,1,1,-1\\}$ for $N{=}6$ sites with periodic
boundary conditions.

Fig. 2 shows the constrained Hilbert space and the action of the Hamiltonian
in this sector for $N{=}6$ sites with PBC. In the graph, each node corresponds
to a product state of the subspace and the edges connect the configurations
that result from a given product state due to the action of the Hamiltonian.
This graph representation will be helpful in the forward scattering
approximation (FSA) defined later in this section. We now study the dynamics
of the basis states of this subspace using the exact diagonalization of the
Hamiltonian. We find that initial states of the kind
$|yxy\rangle{\equiv}|yxyyxy\cdots yxy\rangle$ and
$|xyx\rangle{\equiv}|xyxxyx\cdots xyx\rangle$ show long-lived revivals. Fig.
3(a) depicts the evolution of the initial state $|yxy\rangle$ and a randomly
chosen product state under the Kitaev Hamiltonian with $N{=}24$ sites. The
$|yxy\rangle$ state shows the hallmark of QMBS wherein fidelity oscillations
are robust and long-lived. In particular, the fidelity oscillations for this
state are more robust (peaks heights are higher as evidenced by the fact that
the first revival peak displays ${>}80\%$ return probability to the initial
state) and longer-lived (persist for a longer time) as compared to that of the
$|Z_{2}\rangle_{\text{Kitaev}}$ state shown in Fig. 1. In sharp contrast, a
random state thermalizes rapidly [see the green curve shown in Fig. 3(a)].

Figure 3: (a) Fidelity dynamics for the initial state $|yxy\rangle$ and a
randomly chosen product state. The former shows coherent oscillations whereas
the latter thermalizes rapidly. The black dotted lines show the probability of
state transfer between $|yxy\rangle$ and $|xyx\rangle$ (b) Entanglement
entropy S for the same initial states. The Entropy of the $|yxy\rangle$ state
grows linearly with time but it also shows oscillations due to scarring. Data
is shown for N=24 sites with PBC.

We can visualize the scarred dynamics in this sector as the state bouncing
between the two corner states $|yxy\rangle$ and $|xyx\rangle$ of the Hilbert
space graph shown in Fig. 2. The dotted line in Fig. 3(a), where we plot
$|\langle xyx|e^{-iHt}|yxy\rangle|^{2}$, which is the probability of finding
the state in $|xyx\rangle$ following time evolution from the initial state
$|yxy\rangle$, illustrates this back-and-forth motion. In Fig. 3(b) we plot
the growth of EE with time for different initial states. The EE of a subregion
$A$ is defined as the von Neumann entropy of the reduced density matrix of the
subsystem as
$S_{A}{=}{-}\text{Tr${}_{A^{c}}$}\\{\rho_{A}\text{ln}\rho_{A}\\}$, where
$\rho_{A}$ is the reduced density matrix of subsystem $A$ and the trace is
taken over its complement $A^{c}$. For the $|yxy\rangle$ state, along with an
increase as a function of time, the EE mirrors the oscillations seen in the
fidelity. Moreover, the rate at which EE grows in the $|yxy\rangle$ state is
much smaller as compared to that of a randomly chosen state, suggesting that
the initial $|yxy\rangle$ state results in non-ergodic behavior.

Thermalization and its breakdown can also be probed by measuring the spread of
EE of eigenstates. ETH predicts a “volume-law” scaling of EE, i.e., for a 1D
system EE scales linearly with the size of the subsystem. Fig. 4 demonstrates
that the bipartite (equipartitioned) $S$ for the majority of the eigenstates
of the $\vec{b}{=}\\{1,1,{-}1,1,1-1,{\cdots},1,1,{-}1\\}$ sector do exhibit
the volume-law behavior that is consistent with the prediction of ETH.
However, in addition to the bulk of highly entangled states, there are
outliers over the entire range of the spectrum that have much lower entropy
that violates the volume law predicted by ETH.

Figure 4: Bipartite (equipartition) entanglement entropy of the eigenstates of
the $\vec{b}{=}\\{1,1,{-}1,1,1-1,{\cdots},1,1,{-}1\\}$ in the spin-$1$ Kitaev
model as a function of energy. The bulk states satisfy the volume law of EE,
however, there are several states that reside over the entire range of the
spectrum and carry low EE and thereby violate ETH. Data is shown for $N{=}24$
sites with PBC. The color scale on the right indicates the density of the data
points.

The fidelity oscillations observed in Fig. 3(a) arise precisely due to the
existence of these relatively small number of athermal eigenstates that are
spread throughout the bulk spectrum but carry low EE. These scarred states
have anomalously high overlap with the initial product state $|{yxy}\rangle$
as shown in Fig. 5. Like in other models hosting QMBS [15, 34], the projection
onto the $|{yxy}\rangle$ state displays towers of special _equispaced_ (with
the spacing in energy determining the inverse time period of oscillations seen
in the fidelity) eigenstates having anomalously high overlap with the initial
product state $|{yxy}\rangle$. The observed coherent oscillations in fidelity,
sub-thermal entanglement entropy of certain eigenstates, and anomalously large
overlap of these eigenstates with a particular initial product state results
in the non-ergodic dynamics and establishes the existence of QMBS in this
subspace of spin-$1$ Kitaev chain.

Figure 5: Density plot showing the overlap of the initial product state
$|{yxy}\rangle$ with the eigenstates of the
$\vec{b}{=}\\{1,1,{-}1,1,1-1,{\cdots},1,1,{-}1\\}$ sector in the spin-$1$
Kitaev model. The cross marks denote the overlap with the eigenstates of the
FSA Hamiltonian [see Eq. (20)], which very well approximate the topmost states
in the towers of scar states.

Furthermore, as has been demonstrated for the PXP mode, the topmost state in
the towers of scarred eigenstates in this sector can be well-approximated
using the so-called Forward Scattering Approximation (FSA) [15]. The FSA
mechanism involves constructing an approximate Hamiltonian whose eigenstates
reproduce the scarred states. We start by splitting the Hamiltonian into
forward and backward propagating parts as $H{=}H^{{+}}{+}H^{{-}}$, where

$\begin{split}H^{+}{=}\sum_{i=1,4,7,\cdots}|{yx}\rangle\langle{zz}|+\sum_{i=2,5,6,\cdots}|{zz}\rangle\langle{yx}|\\\
-\sum_{i=1,4,7,\cdots}|{yz}\rangle\langle{zx}|\end{split}$ (16)
$\begin{split}H^{-}{=}\sum_{i=1,4,7,\cdots}|{zz}\rangle\langle{yx}|+\sum_{i=2,5,6,\cdots}|{yx}\rangle\langle{zz}|\\\
-\sum_{i=1,4,7,\cdots}|{zx}\rangle\langle{yz}|.\end{split}$ (17)

Then we construct the basis vectors
$|{0}\rangle,|{1}\rangle,\cdots,|{N}\rangle$ of the effective Hamiltonian
$H_{\rm FSA}$, where $|{0}\rangle{\equiv}|{yxy}\rangle$ and
$|{n}\rangle{=}(1/\sqrt{c_{n}})(H^{+})^{n}|{yxy}\rangle$ ($c_{n}$ is the
normalization constant). In the Hilbert space graph shown in Fig. 2, the
action of $H^{+}$ corresponds to moving from left to right (right to left for
$H^{-}$) and $H^{+}$ annihilates the $|xyx\rangle$ state ($H^{-}$ annihilates
the $|yxy\rangle$ state). Therefore, starting from the $|yxy\rangle$ state the
FSA recurrence closes after $N{+}1$ steps once the forward propagation reaches
the $|xyx\rangle$ state at the opposite end of the Hilbert space graph shown
in Fig. 2. The approximation in FSA entails that the Hilbert space of basis
states $|0\rangle,|1\rangle,\cdots,|N\rangle$ is closed under the action of
the Kitaev Hamiltonian of Eq. (8). The action of $H$ on these basis states is
given by

$\begin{split}H|n\rangle=&H^{+}|n\rangle+H^{-}|n\rangle\\\
=&\beta_{n+1}|n+1\rangle+H^{-}|n\rangle,\end{split}$ (18)

where $\beta_{n}{=}\sqrt{c_{n}/c_{n-1}}$. Thus to make the Hilbert space
closed under the action of the Hamiltonian we have to approximate

$H^{-}|n\rangle\approx\beta_{n}|n-1\rangle.$ (19)

Using Eqs. (18) and (19) the Hamiltonian takes the form of the following
tridiagonal matrix which is the FSA Hamiltonian

$H_{\rm FSA}=\begin{pmatrix}0&\beta_{1}&&&\\\ \beta_{1}&0&\beta_{2}&&\\\
&\beta_{2}&0&\ddots&\\\ &&\ddots&\ddots&\beta_{N}\\\
&&&\beta_{N}&0\end{pmatrix}.$ (20)

As shown by the cross marks in Fig. 5, the eigenstates of the FSA Hamiltonian
of Eq. (20) provide an excellent approximation to the special scarred
eigenstates of the Kitaev Hamiltonian in the sector
$\vec{b}{=}\\{1,1,{-}1,1,1,{-}1,{\cdots},1,1,{-}1\\}$.

###  III.2 QMBS in other sectors

We have also looked for the possibility of scarring in other sectors of the
Kitaev chain by studying the dynamics from different initial product states.
Amongst all the initial states and sectors we considered, we found that
initial product states $|yyxx\rangle{\equiv}|yyxxyyxx\cdots yyxx\rangle$ in
the sector $\vec{b}{=}\\{{-}1,1,{-}1,1,\cdots,{-}1,1\\}$ and the state
$|yyyx\rangle{\equiv}|yyyxyyyx\cdots yyyx\rangle$ in the sector
$\vec{b}{=}\\{{-}1,{-}1,1,1,{-}1,{-}1,1,1,\cdots,{-}1,{-}1,1,1\\}$ also show
oscillations in the fidelity. In Fig. 6 we plot the return probabilities for
the aforementioned initial states. We note that the oscillations are weaker
(as evidenced by the peak heights) and decay much faster in these sectors. One
way to understand this is that the FSA does not work well in these sectors as
shown in Fig. 7.

Figure 6: Fidelity of the state (a) $|yyxx\rangle$ which lies in the sector
$\vec{b}{=}\\{{-}1,1,{-}1,1,{\cdots},{-}1,1\\}$ (b) $|yyyx\rangle$ which lies
in the sector
$\vec{b}{=}\\{{-}1,{-}1,1,1,{-}1,{-}1,1,1,{\cdots},{-}1,{-}1,1,1\\}$. The
fidelity oscillations are weaker (peak heights are reduced) and decay faster
in these sectors as compared to the ones shown in Fig. 3. Data are shown for
$N{=}24$ sites.

Figure 7: Density plot of the overlap of the initial product states that show
fidelity oscillations in particular sectors (see Fig. 6) of the spin-$1$
Kitaev model with all eigenstates in that sector. The crosses show the overlap
with eigenstates of the forward-scattering approximation Hamiltonian $H_{\rm
FSA}$. The overlap of the scar states with the product states fails to match
the magnitude as predicted by FSA which is consistent with the observation of
a faster decay of fidelity in these states. Data are shown for (a)
$\vec{b}{=}\\{{-}1,1,{-}1,1,\cdots,{-}1,1\\}$ and (b)
$\vec{b}{=}\\{{-}1,{-}1,1,1,{-}1,{-}1,1,1,\cdots,{-}1,{-}1,1,1\\}$ sectors of
the spin-$1$ Kitaev chain of $N{=}24$ sites.

The authors of Ref. [20] showed that Hamiltonians hosting QMBS support an
emergent approximate $SU(2)$ symmetry within a subspace of the Hilbert space.
The revivals from the initial product states can then be thought of as the
coherent rotation of the large $SU(2)$ degree of freedom. In the $SU(2)$
algebra, $H^{+}$ and $H^{-}$ act as an analog of the raising and lowering
operators. Their commutator $H^{z}{=}[H^{+},H^{-}]$ acts as $S^{z}$ operator,
and the FSA states are its eigenstates. However, in Eq. (19) we saw that
$H^{-}$ only approximately inverts the action of $H^{+}$, which is why the
algebra is not exact and perfect revivals are not observed. We expect that if
the FSA gives a good representation of the exact scar states, then one sees
strong revivals. Otherwise, if the FSA does not represent the scar states
well, the fidelity decays quickly. This is consistent with our numerical
observations.

These results also show that it is not necessarily the case that the more
constrained a Hilbert space is the more scarring it shows. In general, the
more the number of $b_{j}{=}{-}1$ larger the number of constraints and the
fewer the number of states in the corresponding Hilbert space. Nevertheless,
as we have shown above, the more constrained sector
$\vec{b}{=}\\{{-}1,{-}1,1,1,{-}1,{-}1,1,1,{\cdots}\\}$ shows less scarring
than the less constrained $\vec{b}{=}\\{1,1,{-}1,1,1,{-}1,{\cdots}\\}$ sector.
The strength of scarring is determined by the structure of the graph of the
corresponding Hilbert space and how well the FSA works there. We note here
that we have checked that all sectors in the spin-$1$ Kitaev chain that we
considered do not show Poison level statistics which rules out an
integrability-based explanation for the athermal behavior we observe.

##  IV Summary and Conclusion

In this paper, we studied the time evolution of initial states in certain
sectors of the spin-$1$ Kitaev chain, where the Hilbert space is fragmented
into $2^{N}$ unequal subspaces. We looked at the dynamics of the initial
states in these constrained subspaces. We found that the $|yxy\rangle$ state
in the $\vec{b}{=}\\{1,1,{-}1,1,1,{-}1,\cdots\\}$ sector showed the most
prominent coherent oscillations in fidelity when evolved under the Kitaev
Hamiltonian, more so than even the $|Z_{2}\rangle_{\rm Kitaev}$ state of the
celebrated PXP model which is embedded in the $\vec{b}{=}\\{1,1,\cdots,1\\}$
sector that hosts the ground state of the spin-$1$ Kitaev chain. The coherent
dynamics in the $\vec{b}{=}\\{1,1,{-}1,1,1,{-}1,\cdots\\}$ sector were
characterized by special eigenstates that have anomalously low entanglement
entropy, and high overlap with the initial $|yxy\rangle$ state. We also showed
that these special scarred states can be well-approximated by the FSA.
Furthermore using the FSA, we showed why certain other sectors do not show
long-lived oscillations in fidelity. It would be interesting to see if the
scarring phenomena we observed can be understood using analytical techniques
such as the recently proposed broken unitary picture of dynamics in QMBS [38],
or interpreting it as a one-dimensional chiral scattering problem [39], or
projector-embedding [40] or commutant algebras [41].

The Kitaev chain provides a model system and framework to study the dynamics
of constrained systems. Here, we only looked at the spin-$1$ chain and it
would be interesting to look at higher spins and see whether they exhibit
QMBS. Another potential avenue that could be worth studying in the future is
to explore the existence of QMBS in higher dimensions and/or in different
geometries.

## Acknowledgments

We acknowledge useful discussions with Diptiman Sen, Kartiek Agarwal, Sanjay
Moudgalya, and Zlatko Papić. Computational portions of this research work were
conducted using the Nandadevi supercomputer, which is maintained and supported
by the Institute of Mathematical Science’s High-Performance Computing Center.

##  Appendix A Mapping the $\vec{b}{=}\\{1,1,{\cdots},1\\}$ sector to the PXP
model

In this appendix, we show that the $\vec{b}{=}\\{1,1,{\cdots},1\\}$ sector of
the spin-$1$ Kitaev chain can be mapped to the PXP model. Since all
$b_{j}{=}1$ in this sector, there are five allowed states for any pair of
nearest neighbor sites $\langle j,j{+}1\rangle$ as shown in Eq. (13). Owing to
these constraints, configurations of neighboring sites can be written in terms
of spin-$1/2$ degrees of freedom on the dual lattice (for each bond between
sites $j$ and $j{+}1$ on the primal lattice, on the dual lattice we define a
site at $\\{j{+}1/2\\}$) of $N$ sites using the following mapping

$\begin{split}&|yx\rangle_{j,j+1}\rightarrow|\uparrow\rangle_{j+\frac{1}{2}}\\\
&|zy\rangle_{j,j+1}\rightarrow|\downarrow\rangle_{j+\frac{1}{2}}\\\
&|xz\rangle_{j,j+1}\rightarrow|\downarrow\rangle_{j+\frac{1}{2}}\\\
&|zz\rangle_{j,j+1}\rightarrow|\downarrow\rangle_{j+\frac{1}{2}}\\\
&|xy\rangle_{j,j+1}\rightarrow|\downarrow\rangle_{j+\frac{1}{2}}.\end{split}$
(21)

This mapping does not allow the nearest neighbors on the dual lattice to be in
the configuration $|{\uparrow}{\uparrow}\rangle$ which is precisely the
Rydberg blockade constraint (no two nearest neighbors are simultaneously in
the excited state) that is implemented in the PXP model. Though the mapping in
Eq. (LABEL:Eq:_A1) appears to be many-to-one it is not. The reverse mapping
from the dual lattice to the spin-$1$ primal lattice is given by

$\begin{split}&|\downarrow\uparrow\rangle_{j-\frac{1}{2},j+\frac{1}{2}}\rightarrow|y\rangle_{j}\\\
&|\uparrow\downarrow\rangle_{j-\frac{1}{2},j+\frac{1}{2}}\rightarrow|x\rangle_{j}\\\
&|\downarrow\downarrow\rangle_{j-\frac{1}{2},j+\frac{1}{2}}\rightarrow|z\rangle_{j},\end{split}$
(22)

which ensures that the mapping is one-to-one. Note that a similar mapping was
used in Ref. [34] to map the thin-torus limit of the pair-hoping Hamiltonian
of the $\nu{=}1/3$ fractional quantum hall effect to the PXP model. With this
mapping, the (non-vanishing) action of the spin-$1$ Kitaev model on the primal
lattice leads to the following terms in the Hamiltonian
$\mathcal{H}_{\\{1,1,{\cdots},1\\}}$ in the dual space

$\begin{split}H_{j,j+1}|\ast\overset{j}{y}\overset{j+1}{x}\ast\rangle=|\ast\overset{j}{z}\overset{j+1}{z}\ast\rangle\\\
\implies\mathcal{H}_{j+\frac{1}{2}}|\downarrow\overset{j+\frac{1}{2}}{\uparrow}\downarrow\rangle=|\downarrow\overset{j+\frac{1}{2}}{\downarrow}\downarrow\rangle\\\
H_{j,j+1}|\ast\overset{j}{z}\overset{j+1}{z}\ast\rangle=|\ast\overset{j}{y}\overset{j+1}{x}\ast\rangle\\\
\implies\mathcal{H}_{j+\frac{1}{2}}|\downarrow\overset{j+\frac{1}{2}}{\downarrow}\downarrow\rangle=|\downarrow\overset{j+\frac{1}{2}}{\uparrow}\downarrow\rangle,\end{split}$
(23)

where $|\ast\rangle$ corresponds to any allowed configuration on the sites
that respect the above-mentioned constraints. The terms in the Hamiltonian of
Eq. (23) are exactly the ones that appear in the PXP model, which is given by

$\mathcal{H}_{\\{1,1,\cdots,1\\}}=\sum_{j=1}^{N}P_{j-1}\sigma_{j}^{x}P_{j+1},$
(24)

where the Pauli operators are defined in the usual way with
$\sigma^{x}{=}(|{\downarrow}\rangle\langle{\uparrow}|{+}|{\uparrow}\rangle\langle{\downarrow}|)$
and
$\sigma^{z}{=}(|{\uparrow}\rangle\langle{\uparrow}|{-}|{\downarrow}\rangle\langle{\downarrow}|)$,
and the projectors $P_{j}{=}(1{-}\sigma_{j}^{z})/2$ ensure that no two
spin-$1/2$ nearest neighbors are simultaneously in the excited
$|{\uparrow}\rangle$ state. Thus, the spin-$1$ Kitaev chain Hamiltonian
restricted to the $\vec{b}{=}\\{1,1,\cdots,1\\}$ sector (which hosts its
ground state) exactly maps into the PXP model.

## References

  * Kinoshita  _et al._ [2006] T. Kinoshita, T. Wenger, and D. S. Weiss, Nature 440, 900 (2006).
  * Smith  _et al._ [2016] J. Smith, A. Lee, P. Richerme, B. Neyenhuis, P. W. Hess, P. Hauke, M. Heyl, D. A. Huse, and C. Monroe, Nature Physics 12, 907 (2016).
  * Kucsko  _et al._ [2018] G. Kucsko, S. Choi, J. Choi, P. C. Maurer, H. Zhou, R. Landig, H. Sumiya, S. Onoda, J. Isoya, F. Jelezko, E. Demler, N. Y. Yao, and M. D. Lukin, Phys. Rev. Lett. 121, 023601 (2018).
  * Deutsch [1991] J. M. Deutsch, Phys. Rev. A 43, 2046 (1991).
  * Srednicki [1994] M. Srednicki, Physical review e 50, 888 (1994).
  * Basko  _et al._ [2006] D. Basko, I. Aleiner, and B. Altshuler, Annals of Physics 321, 1126 (2006).
  * Pal and Huse [2010] A. Pal and D. A. Huse, Phys. Rev. B 82, 174411 (2010).
  * Nandkishore and Huse [2015] R. Nandkishore and D. A. Huse, Annual Review of Condensed Matter Physics 6, 15 (2015), https://doi.org/10.1146/annurev-conmatphys-031214-014726 .
  * Schulz  _et al._ [2019] M. Schulz, C. A. Hooley, R. Moessner, and F. Pollmann, Phys. Rev. Lett. 122, 040606 (2019).
  * Kjäll  _et al._ [2014] J. A. Kjäll, J. H. Bardarson, and F. Pollmann, Phys. Rev. Lett. 113, 107204 (2014).
  * Bernien  _et al._ [2017] H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Omran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, Nature 551, 579 (2017).
  * Fendley  _et al._ [2004] P. Fendley, K. Sengupta, and S. Sachdev, Phys. Rev. B 69, 075106 (2004).
  * Urban  _et al._ [2009] E. Urban, T. A. Johnson, T. Henage, L. Isenhower, D. D. Yavuz, T. G. Walker, and M. Saffman, Nature Physics 5, 110 (2009).
  * Lesanovsky [2011] I. Lesanovsky, Phys. Rev. Lett. 106, 025301 (2011).
  * Turner  _et al._ [2018a] C. J. Turner, A. A. Michailidis, D. A. Abanin, M. Serbyn, and Z. Papić, Phys. Rev. B 98, 155134 (2018a).
  * Turner  _et al._ [2018b] C. J. Turner, A. A. Michailidis, D. A. Abanin, M. Serbyn, and Z. Papić, Nature Physics 14, 745 (2018b).
  * Khemani  _et al._ [2019] V. Khemani, C. R. Laumann, and A. Chandran, Phys. Rev. B 99, 161101 (2019).
  * Lin and Motrunich [2019] C.-J. Lin and O. I. Motrunich, Phys. Rev. Lett. 122, 173401 (2019).
  * Moudgalya  _et al._ [2021] S. Moudgalya, B. A. Bernevig, and N. Regnault, arXiv preprint arXiv:2109.00548 https://doi.org/10.1088/1361-6633/ac73a0 (2021).
  * Choi  _et al._ [2019] S. Choi, C. J. Turner, H. Pichler, W. W. Ho, A. A. Michailidis, Z. Papić, M. Serbyn, M. D. Lukin, and D. A. Abanin, Phys. Rev. Lett. 122, 220603 (2019).
  * Chandran  _et al._ [2023] A. Chandran, T. Iadecola, V. Khemani, and R. Moessner, Annual Review of Condensed Matter Physics 14, null (2023), https://doi.org/10.1146/annurev-conmatphys-031620-101617 .
  * Desaules  _et al._ [2021] J.-Y. Desaules, A. Hudomal, C. J. Turner, and Z. Papić, Phys. Rev. Lett. 126, 210601 (2021).
  * Khemani  _et al._ [2020] V. Khemani, M. Hermele, and R. Nandkishore, Phys. Rev. B 101, 174204 (2020).
  * Wildeboer  _et al._ [2022] J. Wildeboer, C. M. Langlett, Z.-C. Yang, A. V. Gorshkov, T. Iadecola, and S. Xu, Phys. Rev. B 106, 205142 (2022).
  * Zhang  _et al._ [2023] P. Zhang, H. Dong, Y. Gao, L. Zhao, J. Hao, J.-Y. Desaules, Q. Guo, J. Chen, J. Deng, B. Liu, W. Ren, Y. Yao, X. Zhang, S. Xu, K. Wang, F. Jin, X. Zhu, B. Zhang, H. Li, C. Song, Z. Wang, F. Liu, Z. Papić, L. Ying, H. Wang, and Y.-C. Lai, Nature Physics 19, 120 (2023).
  * Moudgalya  _et al._ [2018a] S. Moudgalya, S. Rachel, B. A. Bernevig, and N. Regnault, Phys. Rev. B 98, 235155 (2018a).
  * Moudgalya  _et al._ [2018b] S. Moudgalya, N. Regnault, and B. A. Bernevig, Phys. Rev. B 98, 235156 (2018b).
  * Mark  _et al._ [2020] D. K. Mark, C.-J. Lin, and O. I. Motrunich, Phys. Rev. B 101, 195131 (2020).
  * Chattopadhyay  _et al._ [2020] S. Chattopadhyay, H. Pichler, M. D. Lukin, and W. W. Ho, Phys. Rev. B 101, 174308 (2020).
  * Schecter and Iadecola [2019] M. Schecter and T. Iadecola, Phys. Rev. Lett. 123, 147201 (2019).
  * Moudgalya  _et al._ [2020a] S. Moudgalya, N. Regnault, and B. A. Bernevig, Phys. Rev. B 102, 085140 (2020a).
  * Mark and Motrunich [2020] D. K. Mark and O. I. Motrunich, Phys. Rev. B 102, 075132 (2020).
  * Bull  _et al._ [2019] K. Bull, I. Martin, and Z. Papić, Phys. Rev. Lett. 123, 030601 (2019).
  * Moudgalya  _et al._ [2020b] S. Moudgalya, B. A. Bernevig, and N. Regnault, Phys. Rev. B 102, 195150 (2020b).
  * Kitaev [2006] A. Kitaev, Annals of Physics 321, 2 (2006), january Special Issue.
  * Sen  _et al._ [2010] D. Sen, R. Shankar, D. Dhar, and K. Ramola, Phys. Rev. B 82, 195435 (2010).
  * You  _et al._ [2022] W.-L. You, Z. Zhao, J. Ren, G. Sun, L. Li, and A. M. Oleś, Phys. Rev. Res. 4, 013103 (2022).
  * Rozon and Agarwal [2023] P.-G. Rozon and K. Agarwal, Broken unitary picture of dynamics in quantum many-body scars (2023).
  * Windt and Pichler [2022] B. Windt and H. Pichler, Phys. Rev. Lett. 128, 090606 (2022).
  * Shiraishi and Mori [2017] N. Shiraishi and T. Mori, Phys. Rev. Lett. 119, 030601 (2017).
  * Moudgalya and Motrunich [2022] S. Moudgalya and O. I. Motrunich, Phys. Rev. X 12, 011050 (2022).