# Efficient On-Chip Communication for Parallel Graph-Analytics on Spatial
Architectures

Khushal Sethi<EMAIL_ADDRESS>Stanford University

###### Abstract.

Large-scale graph processing has drawn great attention in recent years. Most
of the modern-day datacenter workloads can be represented in the form of Graph
Processing such as MapReduce etc. Consequently, a lot of designs for Domain-
Specific Accelerators have been proposed for Graph Processing. Spatial
Architectures have been promising in the execution of Graph Processing, where
the graph is partitioned on several nodes and each node works in parallel.
We conduct experiments to analyze the on-chip movement of data in graph
processing on a Spatial Architecture. Based on the observations, we identify a
data movement bottleneck, in the execution of such highly-parallel processing
accelerators. To mitigate the bottleneck we propose a novel power-law aware
Graph Partitioning and Data Mapping scheme to reduce the communication latency
by minimizing the hop counts on a scalable network-on-chip. The experimental
results on popular graph algorithms show that our implementation makes the
execution $2-5\times$ faster and $2.7-4\times$ energy-efficient by reducing
the data movement time in comparison to a baseline implementation.

Graph processing, Application-Specific Accelerators, On-chip Communication,
Spatial Architecture, Content-Addressable Memory

††ccs: Computer systems organization Spatial Architecture††ccs: Networks
Network on chip††ccs: Computer systems organization Domain Specific
Accelerators

##  1\. INTRODUCTION

Graphs are widely used to model both data and relationships in real-world
problems, such as knowledge mining, modelling complex systems, network
analytics etc.

Popular graph algorithms such as BFS, SSSP, PageRank are widely used in real
world server workloads. The demand of graph processsing workloads is
increasing rapidly with the increase in internet usage.

Hence, a large number of systems for efficient execution of large graph
workloads on CPU, GPUs, FPGAs, Hybrid Memory Cube (HMC) and also custom
accelerators have been proposed in the past (gui2019survey, ).

In this paper, we analyze a Spatial architecture based domain-specific
accelerator that utilizes Content Addressable Memories to accelerate Graph
Applications efficiently with high storage density and efficiency by modeling
both computation and memory in graph processing. Although, several works for
graph processing in-memory have been investigated, we recognize a data
movement bottleneck, in the highly-parallel execution of such accelerators. To
mitigate the bottleneck, we propose a novel ”power-law” aware Graph
Partitioning and Data Mapping scheme to reduce the communication latency. Our
mapping scheme minimizes hop counts on a scalable Network on Chip. The
experimental results on popular graph algorithms, show that our implementation
achieves $2-5\times$ speedup and $2.7-4\times$ energy reduction in comparison
to the baselines.
The paper is organized as follows : Section 2 provides a background on the
Graph Frameworks (Vertex-Centric Processing) and Spatial Architectures used in
this paper. Section 3 describes our Experimental Setup, Section 4 contains the
discussion for analyzing data movement and our proposed data mapping scheme.
Section 5 lists our key results and performance of our implementation. Section
6 and 7 contain previous related work and conclusions respectively.

##  2\. Background

Table 1. Vertex-centric Programming Model Application | Property | Process | Reduce | Apply
---|---|---|---|---
Breadth First Search | Depth | eProp=u.Prop+1 | u.Temp=min(e.Prop,u.Temp) | u.Prop=min(u.Prop,u.Temp)
SSSP | Distance | eProp=u.Prop+weight | u.temp=min(eProp, u.Temp) | u.Temp=min(u.Prop,u.Temp)
Page Rank | Page rank score | eProp=u.prop+weight | u.Temp=u.Temp+eProp | u.Prop=a*u.Prop+base

###  2.1. Graph Frameworks

Algorithm 1 Vertex-Centric Model

_Process-Phase_ :

for $u:ActiveList$ do

for $v:OutNeighbour(u)$ do

$eProp(u,v)=process(v.Prop,edge(v,u))$

_Reduce-Phase_ :

for $u:AllVertices$ do

for $v:InNeighbour(u)$ do

$eProp(u,v)=reduce(u.Temp,edge(v,u))$

_Apply-Phase_ :

for v: $AllVertices$ do

$v.Prop=Apply(v.Temp,v.Prop))$

The vertex-centric programming model for graph processing consists of three
phases : Process, Reduce and Apply. This functioning takes several iterations
for processing to converge, in each iteration the search frontier of vertices
increases. Algorithm 1 describes the vertex-centric programming model in
graph-processing. The execution is stopped when a converged solution is found.
The Process-Reduce-Apply functions for processing graphs for popular
algorithms in shown in Table 1. The process and apply phases do not contain
dependencies and can be executed in parallel where-as reduce phase should be
executed serially, but can be parallelized by addition of extra fields in the
Graph (zhou2019gram, ).

Figure 1. Spatial Architecture for Graph Processing. Each Processing Engine is
connected with a NOC (a) 2D Mesh (b) Flattened Butterfly. Figure 2. In-Memory
Graph Processing Flow

###  2.2. Spatial Architectures based Domain Specific Accelerators

Many spatial graph processing architectures/accelerators based have been
proposed in the literature (e.g., ReRAM based RPBFS (han2018novel, ), GraphR
(song2018graphr, ), GraphSAR (dai2019graphsar, ) and HyVE (huang2018hyve, )),
achieving great speedup and energy efficiency improvement. The implementations
in GraphR, GraphSAR and HyVE utilize the matrix-vector multiplication
capabilitites of ReRAM in analog domain. In the design of RPBFS, GraphR, and
HyVE, edges in the graph are stored using the edge list format. GraphR
converts stores the graph in edge list format and converts it into adjacency
matrix for MVM operation and stores back into the ReRAM crossbar. GraphSAR
proposed a sparsity-aware data partitioning for improving the performance of
GraphR. HyVE proposed a design for a memory hierarchy of hetergenous storage
over various levels of caches and ReRAM to optimize its execution, which
involves off-chip data movement. HyVE uses conventional CMOS circuits to
process the subgraphs edge by edge. GRAM (zhou2019gram, ) uses 2-D ReCAM based
architecture, which operates search operations in the digital domain and all
the blocks are serially connected to a memory controller which dictates the
execution and data movement on a shared bus. RADAR(huangfu2018radar, ) similar
to us analyzes 3D-ReCAM for genome sequencing.
The high-level architecture of a spatial architecture is shown in Fig. 1. It
typically consists of multiple Graph Engines, the fundamental computational
blocks. Each Graph Engine can execute in parallel which provides a high degree
of parallelism for execution in different graph processing stages. The On-Chip
Controller is responsible for broadcasting operations to each Graph Engine.
Each Graph Engine may consist of input and output buffers, an on-node storage
(for storing the partitioned Graph), Arithmetic Logic Units (ALUs) for post
processing the data.
The Graph Engines are connected in a ”memory-centric” NOC (such as 2D Mesh,
DragonFly, Torus) structure by which they can route packets of data.

###  2.3. Data Flow

The processing flow in Fig 2. follows through an in-memory layout. The graph
structure in memory is represented by a Edge Table ET and a Vertex Table VT.
The edge table and vertex table are arranged column-wise in the dense memory
to provide fast content addressable search during the execution (zhou2019gram,
). For graph G, with M edges and N vertices, properties related to Edges and
Vertices have to be stored in the memory. The Process-Reduce-Apply execution
involves movement of data between these Graph Engines shown in Fig. 1. A
parallel search in the Edge list looks up the data and neighbours. Thus large
data movement happens between Edge List and Vertex Prop (vprop), Vertex Temp
(vtemp) data structures. No movement of data is between the Edge Property
(eprop) and Edge list during the execution. Similarly, data movement happens
between the Edge Property (eprop) and Vertex Prop (vprop), Vertex Temp (vtemp)
data structures. The data-structures will be distributed on the spatial
architecture, on-chip content addressable memory at each Graph Engine
containing only a portion of the data structures.

Table 2. Graph Processing Workloads Graph |  $\\#$Vertices |  $\\#$Edges | Description
---|---|---|---
amazon | 304K | 4.3M | Purchasing Network
soc-pokec | 1.6M | 30.6M | Social Network
wiki-topcats (wiki) | 1.8M | 28.5M | Hyperlinks of Wikepedia
ljournal | 5.4M | 78M | Live Journal

##  3\. Experimental Setup

We study the design in a trace-driven simulator to model the functionality for
the spatial architecture considered. We modified open-source GraphMAT
(sundaram2015graphmat, ) to generate traces for our simulation. We demonstrate
on three commonly used graph-processing algorithms: Breadth First Search
(BFS), Single Shortest Source Path (SSSP), and Page-Rank(PR). The vertex-
centric model for each is shown in Table 1. The experiments are run on real
world graphs obtained from the SNAP Datasets (snap, ).

##  4\. Analysis of Data Movement

The amount of data movement (normalized w/ size of graph) is shown in Fig. 3
for 3 different algorithms considered in our experiments. The amount of data
transferred is nearly equal for process and reduce phase, whereas data
movement in apply phase is negligible compared to the first two phases.
PageRank requires more data-movement because it takes more number of
iterations to converge, hence data is updated several times. The bulk of data
movement in Process phase occurs from reading the Edge Table and
correspondingly transferring data to Vertex Property (vprop) and then to Edge
Property (eprop) for updates (Figure 1). In reduce phase also, the data
movement between arises from updating Temporarily stored vertex data (vTemp)
by Edge Property (eprop) data structure and reading Edge Table (ET) for
neighbours (Figure 1). If additional book-keeping is done for parallel-reduce
operations (as suggested in (zhou2019gram, )) the network traffic in
considerably increased further.
Figure 4 shows the skewed distribution of edges on different vertices, where
sometimes even less than 10% of vertices are connected in 90% of the edges.
This skewness is even greater in larger graph-databases. This skewness is
explained by power-law which is that the degree distribution of vertices
vector n(d) with non-zero entries that follows the relation

(1) $n(d)\propto 1/d^{\alpha}$

where $\alpha>0$ is the slope of the power law (in log scale), d is the number
of edges in the vertex of a graph, and n(d) is the number vertices with a
specific degree d.

Figure 3. Total On-Chip Data Movement (Normalized with the size of Graph) in
Graph Processing of different algorithms for Process and Reduce phase is
shown.

##  5\. Data Mapping Optimization

The NoC topology can be represented as a directed graph P(U, F) with each
vertex $u_{i}\in U$ representing a node in the topology and the directed edge
$f_{i,j}\in F$ representing a physical link between the vertices $u_{i}$ and
$u_{j}$ . If the weights of all connections are considered equal
($f_{i,j}\in\\{0,1\\}$) then, latency (T) of a packet through a network can be
summarized as

(2) $T=H*(T_{r}+T_{w})$

where H is the hop count and $T_{r}$ is the per-hop latency through the switch
and $T_{w}$ is the channel link latency. For minimizing the average hop count
an efficient data mapping is required.
The in-memory graph structure denoted by M (V,E) consists of Edge Table (ET)
where edges are denoted by E, Vertex Temp (vtemp), Vertex Property (vprop) and
Edge Property (eprop) data structures, with data movements between them as
described in the previous sections.
The ”power-law” aware graph partitioning and mapping in each of the data
structures and their spatial arrangement is described in Algorithm 2.

###  5.1. Graph Partitioning to Nodes :

The graph is partitioned in a source-cut (edge) fashion. As a first step, the
vertices are sorted by descending order of degree, for ease. The higher degree
vertices will have more number of edges (combined) compared to lower degree
vertices.
For edge-list partitioning, each edge stored in a Node will have one vertex
from a small set of higher degree vertices, from the top of sorted vertex
list. In other words, the edges from these higher degree vertices are
distributed on to the nodes.

Now, the edges/vertices are allocated to nodes in a load-balancing aware
scheme.

Load balancing is carried by using modulo scheduling, which means among the
$v$ vertices in sorted vertex list, are distributed in a cyclic fashion to the
nodes.

Algorithm 2 Graph Partitioning in Nodes

1:Inputs : NoC Topology Graph P (U,F) and M (V,E)

2:Output : M (V,E) $\to$ P (U,F) : $V\to U$, such that, $v_{i}\in V$,
$u_{j}\in U$ and map ($v_{i}$) = $u_{j}$.

3:V’ $\leftarrow$ sort all v $\in V$ in order of out-degree, V $\leftarrow$ V’

4:_Edge-list Partitioning $(\forall u\in U,u.index\in\\{1,4\\})$_:

5:for $v\in V\%{num\\_of\\_nodes}\ and\ e\ni e.src=v$ do

6: while $u.size<u.maxsize$ do

7: $u\leftarrow e$

8: $u.rank=min(v)\forall v\in u$

9:_Vertex Partitioning $(\forall u\in U,u.index\in{2,3})$_

10:for $v\in V\%{num\\_of\\_nodes}$ do

11: while $u.size<u.maxsize$ do

12: $u\leftarrow v$

13: $u.rank=min(v)\forall v\in u$

###  5.2. Scalable Placement of Nodes :

Since each graph engine is an NOC node, each NOC node $u_{i}\in U$ is assigned
a coordinate $(x_{i},y_{i})$ and an additional fields $rank_{i}$ and
$index_{i}$, where the index is integer-allocated (1-4) based on the data
structure (Edge Table, Vertex Prop, Vertex Temp, Edge Prop) in the following
order, and a rank is allocated to link corresponding vertex-cut Nodes with
different index.

Suppose a node $u_{i}$ contains data for Vertex Property, then the $rank_{i}$
indicates the minimum vertex id contained and $index_{i}=0$.

Then with the corresponding node $u_{j}$ containing Edge Table with same
$rank$, we have $f_{i,j}=1$.

The initial arrangement of In-Memory data structures is shown in Fig 4, where
the mapping starts from the center of the cluster.

We define regularity constraints on each on the data-structure, for
predictable and homogeneous data-transfers.

For this, different data-structures are placed in columnar and row fashion in
the spatial architecture, such that positions of the data-structure between
which data-transfer is higher should be closer. The formal constraints for
such are defined in Algorithm 3.

Algorithm 3 Scalable Placement of Nodes

1:Inputs : NoC Topology Graph P (U,F) and M (V,E)

2:_Populating Topology Graph Edges_

3:$\forall u_{i},u_{j}\in U\ni u_{i}.rank=u_{j}.rank$ and
$u_{i}.index\in\\{1,4\\},u_{j}.index\in\\{2,3\\},\ do\ f_{i,j}=1$

4:_Constraints for Regular-Scalable Structure_

5:$\forall u\in U,u.index=1,u.y>0$

6:$\forall u\in U,u.index=4,u.y<k$

7:$\forall u\in U,u.index=2,0<u.y<k,u.x>0$

8:$\forall u\in U,u.index=3,0<u.y<k,u.x>0$

###  5.3. ILP Formulation for Optimal Placement :

For the determmining the 2-D coordinates of each node, on the network on chip
mesh, we formulate it as an Integer Linear Programming (ILP) problem that
finds the coordinates such that it minimizes the total cost (hop count) of the
data transfer between nodes.

As defined above, the data-transfer happens between the nodes $u_{i}$ and
$u_{j}$ for which the edge $f_{i,j}=1$. The cost of data-transfer between
these nodes is equivalent to the topological distance between them. In case of
a 2-D mesh this topological distance would just be the $L_{1}$ norm distance.

This optimization (in Algorithm 4) is solved computationally with the
constraints defined in Algorithm 3 as an ILP.

Algorithm 4 ILP Formulation of the Optimization Objective

1:Inputs : NoC Topology Graph P (U,F)

2:Output : $(x_{i},y_{i})\leftrightarrow u_{i}$

3:H = min
$\sum_{i}\sum_{j}f_{ij}cost((x_{i},y_{i})\leftrightarrow(x_{j},y_{j}))$

4:Where $cost((x_{i},y_{i})\leftrightarrow(x_{j},y_{j}))$ depends on the
Topological Distance,

5:For 2-D Mesh : $cost((x_{i},y_{i})\leftrightarrow(x_{j},y_{j}))$ = $(\mid
x_{i}-x_{j}\mid+\mid y_{i}-y_{j}\mid)$

6:For Flattened Butterfly : $cost((x_{i},y_{i})\leftrightarrow(x_{j},y_{j}))$
= $(\mid x_{i}-x_{j}\mid+\mid y_{i}-y_{j}\mid)$

Fig. 5 shows the decrease in hop-count from the proposed strategy after
finding the node coordinates by solving the ILP, and using those coordinates
in our simulation in a comparison with the randomized mapping strategy for the
nodes.

Figure 4. Schematic of the Placement Strategy of In-Memory Data Structures for
Graph Processing. Figure 5. Decrease in the Average Hop Count with our
proposed placement strategy (Assuming 2D Mesh Topology).

##  6\. Results

###  6.1. Overall Performance

Figure 6. (a) Node Configuration as presented in GRAM (zhou2019gram, ) showing
I/O Buffers, ALUs and Content Addressable Memories (b) 3DReCAM implementation
Figure 7. Speedup from our proposed optimization due to lower data movement :
for a 3D-ReCAM (Resistive Content Addressable Memory) Architecture with
2D-Mesh and F-Butterfly NoC Configuration Figure 8. Energy Consumption
Reduction : for a 3D-ReCAM (Resistive Content Addressable Memory) Architecture
with 2D-Mesh and F-Butterfly NoC Configuration

To characterize the benefit of our optimization compared to previous works, we
consider the spatial architecture (100 Mhz) configuration as described in
(zhou2019gram, ) with a 3D ReCAM in Graph Nodes for doing Parallel Search
(shown in Fig ). We analyze the effect of utilizing 2D-Mesh and F-Butterfly
Network on Chip configuration in the Spatial Accelerator.

The parameters for the simulation are taken from NVSim-CAM (li2016nvsim, ) and
Destiny (poremba2015destiny, ) modelling tools. The parameters of the 2-D Mesh
NOC used are configured from ORION and input-output buffers are modelled using
CACTI 6.5 (muralimanohar2009cacti, ). These are listed in Table 3.

We use the metric of performance in terms of speed-up in execution time and
energy consumption, where our proposed optimization ((in 2D Mesh NoC
Configuration) results in a on average (geometric mean) $5\times$, $5.2\times$
and $5.4\times$ faster and $2.7\times$, $2.58\times$ and $2.4\times$ energy-
efficient for BFS, SSSP and PR respectively.

Similarly, for the flattened butterfly configuration, the execution becomes
$1.8\times$, $1.85\times$ and $1.9\times$ faster and $2.9\times$, $4.0\times$
and $4.0\times$ energy-efficient (in 2D Mesh NoC Configuration) for BFS, SSSP
and PR respectively.

Table 3. Parameters Used CAM Parameters ((zhou2019gram, )) |
---|---
MAT Size | 1024*1024*8
Per Engine Capacity | 1MB
Word | 64 bits
NOC Parameters ((zhou2019gram, ; kahng2011orion, )) |
Frequency | 1GHz (1ns)
Packet Size | 8 bytes
Latency of hops | 1ns
Ports | 4
Configuration | 2D Mesh

##  7\. Related Work

GraphP proposed a HMC-based graph processing system that reduces
communication, by two techniques : “Source-cut” partitioning (which utilizes
duplicating data), and Hierarchical ”non-blocking” communication and for
efficiently routing traffic. We have proposed another considerable
optimization for reducing communication based on efficient mapping of data
NOCs by utilizing ”power-law”. Our work is hence, orthogonal to their approach
and can be improved by combining their approaches.

Previous works have explored mapping static and dynamic tasks efficiently to
processor-centric NOCs ((sahu2013survey, ; yang2016application, ;
maqsood2015dynamic, ; sahu2014application, ; wu2015efficient, )) by using
Integer Linear Programming, Heuristic Search and Deterministic Search.
However, this has not been explored for memory-centric PIM systems with high
parallelism which have a high programmable complexity. In such systems a
regular structure is also required in mapping the in-memory data structures to
the nodes which is scalable and easy to implement in the memory controller.
This requires good choice of methodologies for efficient and regular data and
processing mapping to network-on-chip for processing-in-memory technologies.

##  8\. Conclusions

In this work we analyze spatial architectures proposed to accelerate Graph
Applications.

Content Addressable Memories allow faster search, allowing faster execution,
but in the fast execution, the on-chip traffic becomes a bottleneck in the
execution time.

By utilizing the inherent characteristics of such Graph Processing Workloads,
an efficient mapping is designed to reduce the on-chip traffic and minimize
the communication latency between the spatial blocks. The experimental results
show that our method can minimize the hop counts for routing packets of data
and also significantly decrease the execution time and energy consumption on
CAM-based Spatial Architectures.

## References

  * (1) C.-Y. Gui, L. Zheng, B. He, C. Liu, X.-Y. Chen, X.-F. Liao, and H. Jin, “A survey on graph processing accelerators: Challenges and opportunities,” _Journal of Computer Science and Technology_ , vol. 34, no. 2, pp. 339–371, 2019.
  * (2) M. Zhou, M. Imani, S. Gupta, Y. Kim, and T. Rosing, “Gram: graph processing in a reram-based computational memory,” in _Proceedings of the 24th Asia and South Pacific Design Automation Conference_. ACM, 2019, pp. 591–596.
  * (3) L. Han, Z. Shen, D. Liu, Z. Shao, H. H. Huang, and T. Li, “A novel reram-based processing-in-memory architecture for graph traversal,” _ACM Transactions on Storage (TOS)_ , vol. 14, no. 1, p. 9, 2018.
  * (4) L. Song, Y. Zhuo, X. Qian, H. Li, and Y. Chen, “Graphr: Accelerating graph processing using reram,” in _2018 IEEE International Symposium on High Performance Computer Architecture (HPCA)_. IEEE, 2018, pp. 531–543.
  * (5) G. Dai, T. Huang, Y. Wang, H. Yang, and J. Wawrzynek, “Graphsar: a sparsity-aware processing-in-memory architecture for large-scale graph processing on rerams,” in _Proceedings of the 24th Asia and South Pacific Design Automation Conference_. ACM, 2019, pp. 120–126.
  * (6) T. Huang, G. Dai, Y. Wang, and H. Yang, “Hyve: Hybrid vertex-edge memory hierarchy for energy-efficient graph processing,” in _2018 Design, Automation & Test in Europe Conference & Exhibition (DATE)_. IEEE, 2018, pp. 973–978.
  * (7) W. Huangfu, S. Li, X. Hu, and Y. Xie, “Radar: a 3d-reram based dna alignment accelerator architecture,” in _2018 55th ACM/ESDA/IEEE Design Automation Conference (DAC)_. IEEE, 2018, pp. 1–6.
  * (8) N. Sundaram, N. Satish, M. M. A. Patwary, S. R. Dulloor, M. J. Anderson, S. G. Vadlamudi, D. Das, and P. Dubey, “Graphmat: High performance graph analytics made productive,” _Proceedings of the VLDB Endowment_ , vol. 8, no. 11, pp. 1214–1225, 2015.
  * (9) “Stanford snap datasets.” [Online]. Available: https://snap.stanford.edu/data/
  * (10) S. Li, L. Liu, P. Gu, C. Xu, and Y. Xie, “Nvsim-cam: a circuit-level simulator for emerging nonvolatile memory based content-addressable memory,” in _Proceedings of the 35th International Conference on Computer-Aided Design_. ACM, 2016, p. 2.
  * (11) M. Poremba, S. Mittal, D. Li, J. S. Vetter, and Y. Xie, “Destiny: A tool for modeling emerging 3d nvm and edram caches,” in _Proceedings of the 2015 Design, Automation & Test in Europe Conference & Exhibition_. EDA Consortium, 2015, pp. 1543–1546.
  * (12) N. Muralimanohar, R. Balasubramonian, and N. P. Jouppi, “Cacti 6.0: A tool to model large caches.”
  * (13) A. B. Kahng, B. Li, L.-S. Peh, and K. Samadi, “Orion 2.0: A power-area simulator for interconnection networks,” _IEEE Transactions on Very Large Scale Integration (VLSI) Systems_ , vol. 20, no. 1, pp. 191–196, 2011.
  * (14) P. K. Sahu and S. Chattopadhyay, “A survey on application mapping strategies for network-on-chip design,” _Journal of Systems Architecture_ , vol. 59, no. 1, pp. 60–76, 2013.
  * (15) L. Yang, W. Liu, W. Jiang, M. Li, J. Yi, and E. H.-M. Sha, “Application mapping and scheduling for network-on-chip-based multiprocessor system-on-chip with fine-grain communication optimization,” _IEEE Transactions on Very Large Scale Integration (VLSI) Systems_ , vol. 24, no. 10, pp. 3027–3040, 2016.
  * (16) T. Maqsood, S. Ali, S. U. Malik, and S. A. Madani, “Dynamic task mapping for network-on-chip based systems,” _Journal of Systems Architecture_ , vol. 61, no. 7, pp. 293–306, 2015.
  * (17) P. K. Sahu, T. Shah, K. Manna, and S. Chattopadhyay, “Application mapping onto mesh-based network-on-chip using discrete particle swarm optimization,” _IEEE Transactions on Very Large Scale Integration (VLSI) Systems_ , vol. 22, no. 2, pp. 300–312, 2014.
  * (18) C. Wu, C. Deng, L. Liu, J. Han, J. Chen, S. Yin, and S. Wei, “An efficient application mapping approach for the co-optimization of reliability, energy, and performance in reconfigurable noc architectures,” _IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems_ , vol. 34, no. 8, pp. 1264–1277, 2015.


# Structure of the nearly-degenerate manifold of lattice quasiholes on a torus

Z. Zeybek  Universität Hamburg, Institut für Laserphysik, Luruper Chaussee
149, 22761 Hamburg, Germany  R. O. Umucalılar  Department of Physics, Mimar
Sinan Fine Arts University, 34380 Sisli, Istanbul, Turkey

###### Abstract

We study the nearly-degenerate quasihole manifold of the bosonic Hofstadter-
Hubbard model on a torus, known to host the lattice analog of the Laughlin
state at filling fraction $\nu=1/2$. Away from $\nu=1/2$ and in the presence
of both localized and delocalized quasiholes, the ratio between the
numerically calculated many-body Chern number for certain groups of states and
the number of states in the relevant group turns out to be constant for this
manifold, which is also manifested in the density profile as the depleted
charge of localized quasiholes. Inspired by a zero-mode counting formula
derivable from a generalized Pauli principle, we employ a combinatorial scheme
to account for the splittings in the manifold, allowing us to interpret some
groups of states as the quasihole excitations corresponding to filling
fractions lower than $\nu=1/2$. In this scheme, the many-body Chern number of
subgroups appears as a simple combinatorial factor.

##  I Introduction

The interest in topological phases of matter has tremendously increased since
the discovery of the quantum Hall effect [1, 2] and topological insulators
[3]. The rich phenomena revealed in this context, such as conductance
quantization [4], chiral edge states [5], and excitations with fractional
charge [6] and exchange statistics [7, 8, 9], are theorized to be instrumental
for a wide range of applications such as topological quantum computing [10]
and novel material design [11, 12].

Recent developments in quantum simulation with ultracold atomic [13] and
photonic systems [14] have created an alternative environment for
investigating topologically ordered states such as fractional quantum Hall
(FQH) [15, 16, 17, 18, 19] and fractional Chern insulator (FCI) states, the
latter being exclusively studied for lattices [20, 21, 22, 23, 24, 25].
Validating the topological nature of states harbored in a given system is an
essential component of such investigations, in which one typically looks for
specific signatures such as charge fractionalization [26], braiding statistics
[27, 28] and matching of the excitation and/or the entanglement spectra of
ground states [20, 21, 22, 23] with known results. Counting quasihole (zero-
mode) states in an FQH-type system on a torus is considered to be a viable
method highlighting the topological nature of the problem. Counting formulas
for certain Abelian and non-Abelian quasihole states have already been
provided either using a generalized Pauli principle [20, 29] or in the thin-
torus limit [30].

In this paper, focusing on the Hofstadter-Hubbard model of repulsively
interacting bosonic particles in a finite lattice with torus boundary
conditions, which is known to host the discrete version of the $\nu=1/2$
Laughlin state [17, 18], we scrutinize its nearly degenerate quasihole
manifold with the guidance of a zero-mode counting formula valid for a generic
Laughlin filling fraction $\nu=1/r$. Based on numerical calculations, we find
that, different from the continuum limit where the degeneracy is expected to
be exact, there exist certain subgroups of states within a nearly-degenerate
manifold which could be interpreted as the quasihole states existing around
filling fractions lower than $\nu=1/2$. We observe that the ratio of the
numerically calculated many-body Chern number for the states in a subgroup to
the number of states in that group (which we will call Chern number per state
for brevity) is the same for all subgroups in the whole nearly-degenerate
manifold. This ratio also turns out to be equal to the ratio of the number of
particles to the number of unpinned magnetic flux quanta in the system.
Investigation of the density profile reveals that this ratio also represents
the depleted charge at the position of pinned quasiholes. To classify the
subgroups, we propose to use their many-body Chern numbers which can be
expressed in a combinatorial form following a simple distribution rule. We
believe that such a distribution rule and a compact form of the degeneracy
formula we present, with its direct connection to fundamental topological
concepts like the many-body Chern number and the fractional charge of
quasiholes, could be seen as stimulating clues for new investigations towards
a better understanding of the underlying quantum states, which may even
include more exotic non-Abelian ones that can occur in a setting similar to
ours [31].

##  II Generalized counting formula and the many-body Chern number

As quasihole excitations are among the unique defining features of FQH states,
the matching of the number of low-lying states with that of Laughlin-type
quasihole excitations for given system parameters can be considered to be a
strong signal that the observed states resemble FQH states rather than
competing ones like charge density waves [20].

The number of states lying below the many-body gap of an FQH-like system can
be obtained both through a counting procedure in the thin-torus limit [30] and
by using a generalized Pauli principle as detailed in Appendix A. In this
latter approach, the total number of quasihole excitations of a
$\nu_{0}=N/N^{0}_{\phi}=1/r$ state on the torus can be determined by applying
the so-called $(1,r)$ generalized Pauli principle [29], in which one counts
the number of distinct ways of distributing $N$ particles into $N^{0}_{\phi}$
orbitals obeying the condition that no more than a single particle is admitted
in $r$ consecutive orbitals, as well as conforming with the periodic boundary
conditions imposed by the torus. Using this method, we derived a formula to
account for the ground-state degeneracy of a given system at an arbitrary
filling fraction $\nu=N/N_{\phi}$ on the torus, assuming that the ground-state
manifold corresponds to quasihole excitations of a reference Laughlin state
with filling $\nu_{0}=1/r$:

$\displaystyle\mathcal{D}(N,N_{d},N_{\phi})=\frac{(N_{d}+N-1)!}{N_{d}!(N-1)!}\frac{N_{\phi}}{N},$
(1)

where $N$, $N_{\phi}$, and $N_{d}=N_{\phi}-N/\nu_{0}$ are the number of
particles, the number of flux quanta and the number of delocalized (unpinned)
quasiholes [27] in the system respectively. We note that although the above
formula may not be considered new in the sense that it appears in similar
forms in several other works [30, 20, 21], it is still original to the best of
our knowledge, deriving its originality mainly from its simple present form
with a direct connection to the many-body Chern number as we now explain. Our
formula (1) can be compactly written as $\mathcal{D}=\mathcal{C}/\nu$ with
$\mathcal{C}=\binom{N_{d}+N-1}{N_{d}}$ being the binomial coefficient. Most
interestingly, in the cases we investigated for a lattice system via exact
diagonalization, this combinatorial part $\mathcal{C}$ reproduces the
numerically calculated many-body Chern number $C$ of the ground-state manifold
so that the Chern number per $\mathcal{D}$-fold degenerate states is simply
$\nu$. Similar to the case of a system with $\nu_{0}=1/r$ on the torus, for
which the Chern number of the $r$-fold degenerate ground-state manifold equals
one [32, 33], this result suggests that $\nu$ might correspond to the
fractional prefactor of the quantized Hall conductance in this case as well.

We note that for a filling fraction $\nu=N/N_{\phi}=k/r$, $k$ and $r$ being
coprime, general translational symmetry arguments in continuum lead to an
$r$-fold center-of-mass degeneracy on the torus [34] and the remaining
degeneracy can be ascribed to relative translations [22]. For $k=1$, this
latter degeneracy corresponds to the combinatorial part $\mathcal{C}$ of our
formula, hinting at a potential relation between the many-body Chern number
with the total number of states per momentum sector of the relative
translation operators.

##  III Bosonic Hofstadter-Hubbard model

In this work we study bosons on a two-dimensional square lattice under a
uniform perpendicular effective magnetic field, via considering the so-called
Hofstadter-Hubbard (HH) model with the Hamiltonian:

$\displaystyle H=-t\sum_{\langle
ij\rangle}\left(e^{i2\pi\phi_{ij}}c^{\dagger}_{i}c_{j}+\text{h.c.}\right)+\frac{U}{2}\sum_{i}n_{i}(n_{i}-1),$
(2)

where $c^{\dagger}_{i}$ ($c_{j}$) is the bosonic creation (annihilation)
operator at site $i$ ($j$) and $n_{i}=c^{\dagger}_{i}c_{i}$ is the number
operator. The strength of repulsive on-site interactions is given by $U>0$ and
$t>0$ is the hopping amplitude between nearest-neighbor sites. The magnetic
field is characterized by the magnetic flux quantum per unit cell of the
lattice defined as $\phi=\sum_{\square}\phi_{ij}/(2\pi)$, where the sum is
carried out around a unit cell. For the non-interacting case with $\phi=p/q$,
$p$ and $q$ being coprime, the energy band is split into $q$ sub-bands
exhibiting the fractal Hofstadter butterfly [35, 36]. In the literature, there
are quite a number of works (see e.g. [17, 18]) which have shown that the
ground state of the interacting Hamiltonian (2) has a very large overlap with
the Laughlin state for the filling fraction $\nu=N/N_{\phi}=1/2$ under
appropriate boundary conditions and in the so-called continuum limit $\phi\ll
1$. The robustness of this state has also been demonstrated through finite-
size scaling studies which were performed both in the continuum [37] and
thermodynamic continuum [38] limits, yielding an unambiguous many-body gap
that scales as $1/q$.

Although we perform exact diagonalization for small hardcore bosonic systems,
we employ the lowest Landau level (LLL) projection method [33] in real space
to deal with larger systems. Working with larger systems helps us avoid
boundary effects as much as possible, at the same time allowing us to maintain
a certain number of flux quanta through the whole lattice while we reduce the
flux quanta per unit cell to reach the continuum limit. An additional benefit
of a larger system is related to a possible future study of the braiding phase
of quasiholes, which appears more clearly when the quasiholes are sufficiently
far apart. For this purpose, we project the operators into a space spanned by
a group of low-energy single-particle states forming a narrow band for which
the kinetic term is frozen out just analogously as in the LLL approximation in
continuum. In the continuum limit $p/q\ll 1$, this lowest band is composed of
$N_{\phi}$ lowest-energy states. LLL projected site operators are formed via
defining
$\tilde{c}_{i}^{\dagger}\equiv\sum_{k}\psi^{k}_{i}(\vec{\theta})d_{k}^{\dagger}$,
where $\vec{\theta}$ denotes the extra twist angles modifying magneto-periodic
boundary conditions [18] and $d_{k}^{\dagger}$ creates a particle in the $k$th
lattice LLL state, $\psi^{k}_{i}(\vec{\theta})$ being its $i$th component. The
projected interaction Hamiltonian then becomes

$H_{I}=\frac{U}{2}\sum_{i}\tilde{c}_{i}^{\dagger}\tilde{c}_{i}^{\dagger}\tilde{c}_{i}\tilde{c}_{i}=\frac{U}{2}\sum_{i}\sum_{k,l,m,n}d_{k}^{\dagger}d_{l}^{\dagger}d_{m}d_{n}\\\
\times\psi^{k*}_{i}(\vec{\theta})\psi^{l*}_{i}(\vec{\theta})\psi^{m}_{i}(\vec{\theta})\psi^{n}_{i}(\vec{\theta}).$
(3)

We also add to the model Hamiltonian (2) a repulsive impurity potential
$V_{\rm imp}=Vn_{i}$, creating an offset with large enough strength $V>0$, to
pin a quasihole on the $i$th site. Note that the impurity potential must also
be projected into the LLL. Throughout our numerical calculations for an
$L_{x}\times L_{y}$ lattice, we choose the flux density to be $\phi=1/L_{y}$
to fit a full magnetic unit cell along the $y$ direction (for the chosen
Landau gauge vector potential in the $x$ direction), which leads to
$N_{\phi}=L_{x}$.

Figure 1: Energy spectrum for the quasihole states in the presence of hardcore
interactions for $N=2$ particles with $L_{x}=10$, $L_{y}=15$, $\phi=1/15$. The
gap $\sim 0.15t$ between the quasihole manifold and spurious states is much
larger than the energy spread displayed here. In the legend, possible
groupings of quasihole states based on the reference filling fraction
$\nu_{0}$ (without the boson-fermion distinction corresponding to even-odd
denominators of $\nu_{0}=1/r$; see text) are indicated by different
combinations of symbols. $C$ is the Chern number of a subgroup. The grouping
of the first five states for the possible $\nu_{0}=1/5$ reference filling is
not shown as the splitting is very small compared to other groups and the
Chern number cannot be calculated reliably.

##  IV Grouping the quasihole states

For a bosonic system on a lattice with a filling fraction different from the
Laughlin filling $\nu=1/2$, there exists a manifold of nearly-degenerate
quasihole states harboring delocalized quasiholes in the absence of a pinning
potential [39, 27]. These low-lying quasihole excitations are separated by a
relatively large energy gap from spurious states. In addition to this gap, we
numerically find smaller splittings within the manifold itself, which occur
prominently at an integer multiple of $r$ when the filling fraction is
$\nu=N/N_{\phi}=k/r$, $k$ and $r$ being coprime. Our degeneracy formula (1)
for continuum suggests that some of these splittings might mark the grouping
of quasihole states with respect to a reference filling fraction as shown in
Fig. 1. For instance, a system with $N=2$, $N_{\phi}=10$ with the filling
fraction $\nu=1/5$ can be considered to host bosonic Laughlin quasihole states
at different filling fractions such as $\nu_{0}=1/4$ with $N_{d}=2$
($\mathcal{D}=15$) and $\nu_{0}=1/2$ with $N_{d}=6$ ($\mathcal{D}=35$)
quasiholes. Since for a given number of particles and interaction strength the
energy differences become inappreciable as the size of the system grows (see
Appendix B to see an example of this effect), these energy differences may not
be considered unambiguous observables and it is important to develop a method
to distinguish between different types of quasihole excitations.
Investigating, for instance, the particle entanglement spectrum might help one
to make such a distinction [20, 21, 22].

Here, we propose to label certain groups of states by their many-body Chern
number [40, 41] (see Appendix C) in order to resolve the aforementioned
ambiguity. In previous studies (e.g. [18]) many-body Chern number has been
shown to probe topological order by revealing whether a given collection of
states reflects the properties of a Laughlin state of a certain filling
fraction. In Table 1, we show the numerically calculated many-body Chern
numbers $C$ of certain groups of low-lying states for several systems. Each
group for a given $N$ and $N_{\phi}$ can be represented by a filling fraction
of the Laughlin form $\nu_{0}=1/r$, $r$ being an integer (without the even-odd
distinction) starting from $2$, plus a number $N_{d}=N_{\phi}-N/\nu_{0}$ of
delocalized quasiholes so as to yield the physical filling $\nu=N/N_{\phi}$.
As an odd-denominator filling fraction ($\nu_{0}=1/3$) appears in the Table as
well, here is a good point to stress the fact that we do not claim that all
these groups of states correspond to quasihole excitations of bosonic Laughlin
states, for which only even denominators are allowed. Our motivation for
considering all possible $\nu_{0}=1/r$ fractions is to show that odd-
denominator fractions also fit in the combinatorial picture we present in the
following (as in Figs. 1 and 3), without too much speculating on the physical
nature of the states. By localized or delocalized quasihole states, again we
do not mean the usual bosonic Laughlin quasihole states, but just excitations
over a state with a certain (even or odd) number of flux quanta per particle.

Table 1:  Numerical many-body Chern numbers $C$ calculated using $\mathcal{D}$ lowest-energy states for systems of $N=2,3$ hardcore bosons. $N_{\phi}$ is the number of flux quanta, $\nu_{0}$ is the reference filling fraction and $N_{d}$ is the corresponding number of delocalized quasiholes. Flux quantum per unit cell is $\phi=1/12,1/10$ for $N=2,3$, respectively. Degeneracy $\mathcal{D}$ is calculated using Eq. (1). $N$ | $N_{\phi}$ | $\nu_{0}$ | $N_{d}$ | $\mathcal{D}$ | $C,\mathcal{C}$ | $C/\mathcal{D},\nu$
---|---|---|---|---|---|---
$2$ | $5$ | $1/2$ | $1$ | $5$ | $2$ | $2/5$
$2$ | $6$ | $1/3$ | $0$ | $3$ | $1$ | $1/3$
$2$ | $6$ | $1/2$ | $2$ | $9$ | $3$ | $1/3$
$2$ | $7$ | $1/3$ | $1$ | $7$ | $2$ | $2/7$
$2$ | $7$ | $1/2$ | $3$ | $14$ | $4$ | $2/7$
$2$ | $8$ | $1/4$ | $0$ | $4$ | $1$ | $1/4$
$2$ | $8$ | $1/3$ | $2$ | $12$ | $3$ | $1/4$
$2$ | $8$ | $1/2$ | $4$ | $20$ | $5$ | $1/4$
$3$ | $6$ | $1/2$ | $0$ | $2$ | $1$ | $1/2$
$3$ | $7$ | $1/2$ | $1$ | $7$ | $3$ | $3/7$
$3$ | $8$ | $1/2$ | $2$ | $16$ | $6$ | $3/8$

Using the parameters $N$, $N_{\phi}$, and $N_{d}$, we find the corresponding
degeneracy $\mathcal{D}$ through Eq. (1) and take $\mathcal{D}$ nearly-
degenerate lowest-energy states of the system to calculate $C$ for this group
of states. The most striking features observable in Table 1 are: (i) Chern
number $C$ equals the combinatorial factor $\mathcal{C}$ appearing in Eq. (1),
(ii) Chern number per state $C/\mathcal{D}$ equals the physical filling
fraction $\nu$ just as when $\nu$ itself is a Laughlin filling fraction. This
result also suggests that all the states in the lowest-energy manifold share a
common topological feature, namely $C/\mathcal{D}$. The appropriateness of
this quantity in classifying states is apparent for the two systems displayed
in the table with $\nu=2/5$ and $\nu=3/7$, which can both be interpreted as
containing a single quasihole ($N_{d}=1$) with respect to $\nu_{0}=1/2$, but
each has a different $C/\mathcal{D}$ for the degenerate manifold. Next, we
examine the role of pinning a quasihole in counting the degeneracy and suggest
a connection between $C/\mathcal{D}$ and the charge of the pinned quasihole.

##  V The role of pinning

Adding an offset potential on a single site in the HH Hamiltonian (2) has been
shown to be effective in pinning a single quasihole associated with an added
flux quantum, visible as a density depletion around the site [39, 27]. In
general, the energy spectrum turns out to have a smoother variation in the
presence of a pinning potential washing out the subgroup structure. We
calculated the average spatial density $\expectationvalue{n_{i}}$ over the
states that make up the nearly-degenerate manifold and computed the density
depletion around the quasihole center via defining the missing density (or
charge) with respect to its value in the absence of a pinning potential
($V=0$) with one less flux quanta, which also corresponds to the almost
uniform value of the numerical density calculated when a pinning potential is
present, forming an analogy with the continuum case [27, 28]:

$\displaystyle Q_{\rho}\equiv\sum_{i}[\expectationvalue{n_{i}}_{V\neq
0}-\expectationvalue{n_{i}}_{V=0}],$ (4)

where the sum is carried out for sites that are equidistant from the pinning
site in a given radius $\rho$.

Figure 2: Density (charge) depletion profile $Q_{\rho}$ for a system of $N=3$
hardcore bosons with $V=4t$, $L_{x}=8$, $L_{y}=10$ and $\phi=1/10$ (yielding
$N_{\phi}=8$), as a function of the radial distance $\rho$ from the pinning
site. Dashed line indicates the saturation level $N/(N_{\phi}-1)=3/7$ of the
depletion occuring around $\rho\sim 3a$, which determines the quasihole size.
The inset shows the site densities, interpolated for better visualization.

The fractional nature of the quasihole charge can be probed via the depletion
profile. As seen in Fig. 2, the depletion saturates at a radius compatible
with the quasihole size which can be extracted from the two-dimensional
density. The saturation level then reveals the fractional charge. For the
cases we numerically checked, this fractional charge always turned out to be
equal to the filling $\nu=N/(N_{\phi}-N_{\rm loc})$ of a reference system
without a pinning potential, containing a total of $N_{\phi}-N_{\rm loc}$ flux
quanta, where $N_{\rm loc}$ is the number of localized flux quanta (or pinned
quasiholes) in the actual system. We also observed that the number of states
$D$ in the whole nearly-degenerate manifold of the actual system is the same
as that of the reference system and can be given with Eq. (1) by considering
only the unpinned flux quanta in the formula, through the replacement
$N_{\phi}\rightarrow N_{\phi}-N_{\rm loc}$, as exemplified in Table 2.
Moreover, for the systems given in the first two rows of the table, we
calculated the many-body Chern number for the $7$-fold degenerate lowest band
to be $3$, leading to the remarkable result that the Chern number per state
$3/7$ in both cases is equal to the fractional charge of the pinned quasihole
(Fig. 2), which is different from a simple Laughlin fraction $1/r$. Whether
the quasihole (or more properly the localized density depletion) with charge
$\nu=C/\mathcal{D}$ has any anyonic character with regard to braiding
statistics (at least for simple bosonic fractions) is an interesting open
problem to be pursued in a future work. In order to extract information as to
the braiding phase of two quasiholes, one can either perform a real braiding
operation with time-dependent potentials [31, 27, 42] or take an indirect
route which involves the comparison of the density profiles of one- and two-
quasihole states [28].

Table 2:  Numerical ground-state degeneracy $D$ for various systems obtained via exact diagonalization ($N=3$, hardcore) and LLL approximation ($N=4,5$). When there are $N_{\rm loc}$ localized quasiholes, each is centred at a lattice site through a pinning potential with the same strength. The number of delocalized quasiholes $N_{d}=(N_{\phi}-N_{\rm loc})-N/\nu_{0}$ is defined with respect to $\nu_{0}=1/2$. Symbols next to $D$ values are the same for systems with the same $N$ and $N_{\phi}-N_{\rm loc}$. $N$ | $N_{\phi}$ | $N_{\rm loc}$ | $N_{d}$ | $D$ |
---|---|---|---|---|---
$3$ | $7$ | $0$ | $1$ | $7$ | $\ast$
$3$ | $8$ | $1$ | $1$ | $7$ | $\ast$
$4$ | $9$ | $0$ | $1$ | $9$ | $\bullet$
$4$ | $10$ | $0$ | $2$ | $25$ | $\triangle$
$4$ | $10$ | $1$ | $1$ | $9$ | $\bullet$
$4$ | $11$ | $0$ | $3$ | $55$ |
$4$ | $11$ | $1$ | $2$ | $25$ | $\triangle$
$4$ | $11$ | $2$ | $1$ | $9$ | $\bullet$
$5$ | $11$ | $0$ | $1$ | $11$ | $\square$
$5$ | $12$ | $0$ | $2$ | $36$ |
$5$ | $12$ | $1$ | $1$ | $11$ | $\square$

##  VI Exact degeneracy in the Kapit-Mueller model

Here, we briefly discuss whether the splittings we found in our numerical
calculations can be a finite-size artifact, by considering the Kapit-Mueller
model [43] for which the lowest single-particle band is exactly flat. Instead
of the nearest-neighbor hopping of Eq. (2), this model allows for an infinite-
range hopping with properly decaying amplitudes [44], which leads to the
remarkable result that the lowest single-particle band of an infinite lattice
is spanned by the continuum LLL functions sampled at lattice points. We
calculated the energy spectra and many-body Chern numbers for several small
Kapit-Mueller systems with high $\phi$, in the hardcore limit and using
periodic boundary conditions. For all the cases we present in Table 3, the
many-particle groundstate manifold turned out to be exactly degenerate up to
numerical precision even in the presence of a pinning potential. Moreover, the
degeneracy and the many-body Chern number for this manifold were found to be
as expected, conforming to the formulas presented in previous sections.

Table 3:  Numerical ground-state degeneracy $D$ and many-body Chern number $C$ for Kapit-Mueller systems of $N=2,3$ hardcore bosons with respective flux quanta per unit cell $\phi=2/5,1/3$. The lattice size $L_{x}=L_{y}=5,6$ determines the total number of flux quanta as $N_{\phi}=10,12$. The number $N_{\rm loc}=1(0)$ signifies that a quasihole is (not) localized with a pinning potential. The number of delocalized quasiholes $N_{d}=(N_{\phi}-N_{\rm loc})-N/\nu_{0}$ is defined with respect to $\nu_{0}=1/2$. $N$ | $N_{\phi}$ | $N_{\rm loc}$ | $N_{d}$ | $D$ | $C$ | $C/D$
---|---|---|---|---|---|---
$2$ | $10$ | $0$ | $6$ | $35$ | $7$ | $1/5$
$2$ | $10$ | $1$ | $5$ | $27$ | $6$ | $2/9$
$3$ | $12$ | $0$ | $6$ | $112$ | $28$ | $1/4$
$3$ | $12$ | $1$ | $5$ | $77$ | $21$ | $3/11$

The exact degeneracy appearing in these systems, which are even smaller than
the ones investigated in previous sections with the same boundary conditions
and therefore could in principle be prone to the same finite-size issues,
suggests that the small splittings we obtained for the Hofstadter-Hubbard
model in the continuum limit might be ascribed to the nearest-neighbor tight-
binding nature of the model rather than to the small size of the systems we
investigated. From this standpoint, we may argue that the nearest-neighbor
hopping model might help us better understand the competing many-body phases
in the continuum by lifting the exact degeneracy and leading to an intricate
grouping structure of states.

##  VII Possible generalization of the combinatorial scheme

Until this point, we have treated energy splittings that correspond only to
certain groups of quasihole excitations (as in Fig. 1) and provided a simple
combinatorial factor ($\mathcal{C}$) for the Chern number of subgroups. This
factor can be interpreted as the number of ways to distribute $N_{d}$
indistinguishable objects into $N$ distinguishable boxes. Based on numerical
calculations for larger systems (without pinning), we observed that there
exist extra splittings and correspondingly new groupings in the degenerate
quasihole manifold, which can be accounted for by a more general combinatorial
scheme.

Figure 3: For systems with $\nu=N/N_{\phi}=3/12$ and $\nu=N/N_{\phi}=3/13$ at
various $\phi=p/q\ll 1$ and fixed $U/t$, the largest ten of the energy
differences $\Delta E(n)=E_{n+1}-E_{n}$ between consecutive states of
increasing energy are displayed in logarithmic scale.

In Fig. 3, we display the largest ten mini-gaps between consecutive states of
the lowest-energy manifold for two systems in the continuum limit with various
$\phi=p/q\ll 1$. The largest gap (delineated by a circle) separating the whole
nearly-degenerate manifold from spurious states is quite robust. It is also
clearly seen that certain mini-gaps (delineated by ovals and rectangles) close
in a similar fashion as $\phi$ gets smaller. Upon scrutiny, we realized that
all these relatively slowly closing gaps mark groups of states with degeneracy
$\tilde{D}=\tilde{C}/\nu$, where $\nu=N/N_{\phi}$ and $\tilde{C}$ is a new
combinatorial factor encompassing the previously found $\mathcal{C}$. Now,
$\tilde{C}\equiv(N,N_{d},m)$ corresponds to the number of ways to distribute
$N_{d}$ indistinguishable objects into $N$ distinguishable boxes each of which
can accommodate a maximum of $m$ objects. For $m=N_{d}$, $\tilde{C}$ equals
$\mathcal{C}$. In order to exemplify this new counting, we focus on the system
of Fig. 3(a) with $N=3$ and $N_{\phi}=12$. This system can be thought to
possess excitations at different filling fractions: $\nu_{0}=1/2$ with
$N_{d}=6$, $\nu_{0}=1/3$ with $N_{d}=3$, and $\nu_{0}=1/4$ with $N_{d}=0$. The
degeneracy of the groups is then obtained by computing $\tilde{C}$ for certain
$m$ as follows

$\displaystyle\tilde{C}$ $\displaystyle\tilde{D}=\tilde{C}$
$\displaystyle/\nu$ $\displaystyle(N=3,N_{d}=6$ $\displaystyle,m=6)=28$
$\displaystyle 112$ $\displaystyle(N=3,N_{d}=6$ $\displaystyle,m=5)=25$
$\displaystyle 100$ $\displaystyle(N=3,N_{d}=6$ $\displaystyle,m=4)=14$
$\displaystyle 76$ $\displaystyle(N=3,N_{d}=3$ $\displaystyle,m=3)=10$
$\displaystyle 40$ $\displaystyle(N=3,N_{d}=3$ $\displaystyle,m=2)=7$
$\displaystyle 28$ $\displaystyle(N=3,N_{d}=0$ $\displaystyle,m=0)=1$
$\displaystyle 4$

Results for the system of Fig. 3(b) follow analogously. We have also
encountered extra slow-closing gaps in numerical results for various $N=4$
systems, for which we could not find a distinct pattern. Although we were not
able to perform a reliable Chern number calculation using the LLL projection
method, the relation $\nu=\tilde{C}/\tilde{D}$ suggests that this new
combinatorial factor $\tilde{C}$ can be thought to be the many-body Chern
number of the $\tilde{D}$-fold nearly-degenerate states, by analogy with our
previous findings.

##  VIII Conclusion

We presented a compact expression for the degeneracy of Laughlin quasihole
states at different filling fractions on a torus, which can be written as the
inverse physical filling fraction times a simple combinatorial factor. Based
on numerical many-body Chern number calculations for a lattice model, we
observed that the combinatorial part equals the Chern number of certain
subgroups of states suggesting that all states in the lowest-energy manifold
can be characterized by the same Chern number per state, which also equals the
physical filling fraction. In systems containing pinned quasiholes, we showed
that the quasihole charge is controlled by the number of unpinned flux quanta
and is equal to the Chern number per state as well.

In order to study larger systems with slightly larger number of particles and
flux quanta in the continuum limit, which is characterized by small flux
quanta per unit cell, we employed a lowest-band approximation. As one
approaches the continuum limit of an exactly flat manifold of states, the
energy-splittings get smaller. By examining the closing structure of these
mini-gaps, we also proposed a more general combinatorial scheme to account for
additional groupings that do not correspond to quasihole states at simple
$\nu_{0}=1/r$ fractions. We briefly touched upon the possibility of these
energy splittings being a finite-size artifact as well, by studying the Kapit-
Mueller model of long-range hoppings for small systems. Having found the many-
particle ground state manifold to be exactly flat for these small systems, we
speculated that the origin of the small splittings we encounter, with their
specific grouping structure, might be the short-range hopping rather than the
smallness of the system. The grouping structure of the nearly-degenerate
states we identified and the compact form of our degeneracy formula with its
relation to some basic topological invariants can be considered a step forward
towards a better microscopic understanding of the system. As a next step, we
will focus on various FCI systems and the braiding properties of quasiholes in
the direction of expanding our present findings.

###### Acknowledgements.

The authors thank I. Carusotto for inspiring discussions. This work has been
supported through the TÜBİTAK-CNR International Bilateral Cooperation Program
2504 (project no. 119N192). R. O. U. acknowledges support through the TÜBA-
GEBİP Award of the Turkish Academy of Sciences.

##  Appendix A Proof of the generalized counting formula

The counting procedure for the total number of quasihole states of a $\nu=1/r$
Laughlin state of $N$ particles in $N_{\phi}$ orbitals (number of flux quanta)
on the torus can be described in two steps [20]:

  * •

Enumerate the number of ways of occupying $N_{\phi}$ orbitals with $N$
particles with respect to the constraint that the configurations obey the
$(1,r)$ generalized Pauli principle on the sphere [29].

  * •

Discard occupation configurations which violate the torus boundary conditions.

Thus the problem of degeneracy counting can actually be analyzed as two
separate counting problems.

###  A.1 Quasihole Counting on the Sphere

Given $N_{\phi}$ orbitals and $N$ particles we must enumerate all occupation
configurations such that no more than one particle exists in $r$ consecutive
orbitals, according to the $(1,r)$ generalized Pauli principle.

Let $x_{0},x_{1},\dots,x_{N}$ represent the number of unoccupied orbitals
(zeros) between the consecutive occupied orbitals (ones). For instance,
$x_{0}$ is the number of zeros to the left side of the first occupied orbitals
and $x_{1}$ is the number of zeros between the first and second occupied
orbitals, and so forth:

$\underbrace{0}_{x_{0}}1\underbrace{000}_{x_{1}}1\underbrace{000}_{x_{2}}1\underbrace{0000}_{x_{3}}1\dots
1\underbrace{000}_{x_{N}}$

We know there are a total of $N_{\phi}-N$ zeros indicating unoccupied
orbitals. Thus we have

$x_{0}+x_{1}+\dots+x_{N}=N_{\phi}-N.$ (5)

Furthermore, the constraint that no more than one particle exist in $r$
consecutive orbitals implies that for $i\in\\{1,2,\dots,N-1\\}$,
configurations must conform to the condition $x_{i}\geq r-1$. Therefore, the
problem reduces to the number of integer solutions to (5) with the constraint
$x_{i}\geq r-1$ for $i\in\\{1,2,\dots,N-1\\}$:

$\displaystyle\begin{drcases}x_{0}+x_{1}+\dots+x_{N}=N_{\phi}-N\\\ x_{0}\geq
0\\\ x_{1}\geq r-1\\\ x_{2}\geq r-1\\\ \vdots\\\ x_{N-1}\geq r-1\\\ x_{N}\geq
0\end{drcases}$ (6)

In order to have the constraint included inside the integer equation, we make
a change of variables in (6): $a_{0}=x_{0},a_{N}=x_{N},a_{i}=x_{i}+1-r$. Then
the integer equation takes the following form

$\displaystyle\begin{drcases}a_{0}+a_{1}+\dots+a_{N}=N_{\phi}-Nr+r-1\\\
a_{0}\geq 0\\\ a_{1}\geq 0\\\ a_{2}\geq 0\\\ \vdots\\\ a_{N-1}\geq 0\\\
a_{N}\geq 0\end{drcases}$ (7)

with the solution,

$\boxed{\mbox{Number of integer solutions}=\binom{N+n}{n}}$

where $n\equiv N_{\phi}-Nr+r-1$.

###  A.2 Quasihole Counting on the Torus

We now calculate the number of configurations that violate the torus boundary
conditions (TBCs) for a given $r$. TBCs are realized through applying the
1D-chain boundary condition to the $N_{\phi}$-orbital space, i.e., the $i$th
orbital is identified with the $(i+N_{\phi})$th one.

Example: Boundary condition violating patterns for quasihole excitations of
$\nu=\frac{1}{3}$ $1\leavevmode\hbox to6.67pt{\vbox
to6.67pt{\pgfpicture\makeatletter\hbox{\hskip 3.33301pt\lower-3.33301pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\pgfsys@color@gray@fill{0}\pgfsys@invoke{ }\hbox{}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox
to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}00\underline{\hbox
to170.71652pt{\hss(1,3) admissible-patterns\hss}}001\,\leavevmode\hbox
to1.2pt{\vbox to1.2pt{\pgfpicture\makeatletter\hbox{\hskip
0.59999pt\lower-0.59999pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}}{{}}{}
{{}{}}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{1.2pt}\pgfsys@invoke{
}\pgfsys@setdash{3.0pt,3.0pt}{0.0pt}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\,\leavevmode\hbox
to6.67pt{\vbox to6.67pt{\pgfpicture\makeatletter\hbox{\hskip
3.33301pt\lower-3.33301pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\pgfsys@color@gray@fill{0}\pgfsys@invoke{ }\hbox{}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox
to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}1\leavevmode\hbox
to1.5pt{\vbox to12.45pt{\pgfpicture\makeatletter\hbox{\hskip
1.19998pt\lower-10.84619pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{ {}{}{}{}{}}{}{{}}{} {{}{}}{}{
{}{}{}{}{}}{}{}{}{}{{
{\pgfsys@beginscope\pgfsys@setlinewidth{0.32pt}\pgfsys@setdash{}{0.0pt}\pgfsys@roundcap\pgfsys@roundjoin{}
{}{}{} {}{}{}
\pgfsys@moveto{-1.19998pt}{1.59998pt}\pgfsys@curveto{-1.09998pt}{0.99998pt}{0.0pt}{0.09999pt}{0.29999pt}{0.0pt}\pgfsys@curveto{0.0pt}{-0.09999pt}{-1.09998pt}{-0.99998pt}{-1.19998pt}{-1.59998pt}\pgfsys@stroke\pgfsys@endscope}}
}{}{}{{}}\pgfsys@moveto{0.0pt}{-3.533pt}\pgfsys@lineto{0.0pt}{-10.6462pt}\pgfsys@lineto{0.0pt}{-10.6462pt}\pgfsys@lineto{0.0pt}{-3.993pt}\pgfsys@stroke\pgfsys@invoke{
}{{}{{}}{}{}{{}}{{{}}{{{}}{\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{0.0}{1.0}{-1.0}{0.0}{0.0pt}{-3.993pt}\pgfsys@invoke{
}\pgfsys@invoke{ \lxSVG@closescope }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}{{}}}} } \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox
to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}$ $1\leavevmode\hbox
to6.67pt{\vbox to6.67pt{\pgfpicture\makeatletter\hbox{\hskip
3.33301pt\lower-3.33301pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\pgfsys@color@gray@fill{0}\pgfsys@invoke{ }\hbox{}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox
to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}00\underline{\hbox
to170.71652pt{\hss(1,3) admissible-patterns\hss}}0010\,\leavevmode\hbox
to1.2pt{\vbox to1.2pt{\pgfpicture\makeatletter\hbox{\hskip
0.59999pt\lower-0.59999pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}}{{}}{}
{{}{}}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{1.2pt}\pgfsys@invoke{
}\pgfsys@setdash{3.0pt,3.0pt}{0.0pt}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\,1\leavevmode\hbox
to6.67pt{\vbox to6.67pt{\pgfpicture\makeatletter\hbox{\hskip
3.33301pt\lower-3.33301pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\pgfsys@color@gray@fill{0}\pgfsys@invoke{ }\hbox{}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox
to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}0\leavevmode\hbox
to0.4pt{\vbox to7.51pt{\pgfpicture\makeatletter\hbox{\hskip
0.2pt\lower-10.84619pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{ {}{}{}{}{}}{}{{}}{} {{}{}}{}{
{}{}{}{}{}}{}{}{}{}{}{}{}{{}}\pgfsys@moveto{0.0pt}{-3.533pt}\pgfsys@lineto{0.0pt}{-10.6462pt}\pgfsys@lineto{0.0pt}{-10.6462pt}\pgfsys@lineto{0.0pt}{-3.993pt}\pgfsys@stroke\pgfsys@invoke{
}{{}{{}}{}{}{{}}{{{}}{{{}}{\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{0.0}{1.0}{-1.0}{0.0}{0.0pt}{-3.993pt}\pgfsys@invoke{
}\pgfsys@invoke{ \lxSVG@closescope }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}{{}}}} } \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox
to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}$ $0\leavevmode\hbox
to6.67pt{\vbox to6.67pt{\pgfpicture\makeatletter\hbox{\hskip
3.33301pt\lower-3.33301pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\pgfsys@color@gray@fill{0}\pgfsys@invoke{ }\hbox{}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox
to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}100\underline{\hbox
to170.71652pt{\hss(1,3) admissible-patterns\hss}}001\,\leavevmode\hbox
to1.2pt{\vbox to1.2pt{\pgfpicture\makeatletter\hbox{\hskip
0.59999pt\lower-0.59999pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}}{{}}{}
{{}{}}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{1.2pt}\pgfsys@invoke{
}\pgfsys@setdash{3.0pt,3.0pt}{0.0pt}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{0.0pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\,0\leavevmode\hbox
to6.67pt{\vbox to6.67pt{\pgfpicture\makeatletter\hbox{\hskip
3.33301pt\lower-3.33301pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\pgfsys@color@gray@fill{0}\pgfsys@invoke{ }\hbox{}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox
to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}1\leavevmode\hbox
to0.4pt{\vbox to7.51pt{\pgfpicture\makeatletter\hbox{\hskip
0.2pt\lower-10.84619pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{ {}{}{}{}{}}{}{{}}{} {{}{}}{}{
{}{}{}{}{}}{}{}{}{}{}{}{}{{}}\pgfsys@moveto{0.0pt}{-3.533pt}\pgfsys@lineto{0.0pt}{-10.6462pt}\pgfsys@lineto{0.0pt}{-10.6462pt}\pgfsys@lineto{0.0pt}{-3.993pt}\pgfsys@stroke\pgfsys@invoke{
}{{}{{}}{}{}{{}}{{{}}{{{}}{\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{0.0}{1.0}{-1.0}{0.0}{0.0pt}{-3.993pt}\pgfsys@invoke{
}\pgfsys@invoke{ \lxSVG@closescope }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}{{}}}} } \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox
to0.0pt{}{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}$ Above configurations
are all valid on the sphere but not so when TBCs are applied, thus any inner
(1,3) admissible patterns between the edge patterns such as above cannot be
realized on the torus and must be discarded.

As it can be seen from the example, TBCs lead to a reduction in the
degeneracy; therefore, in order to obtain the correct number we must count
invalid configurations and discard them. We will show the counting for the
$\nu=1/3$ case and generalize it for a given $\nu=1/r$.

(a) $\displaystyle 100\text{\text@underline{\hbox to170.71652pt{\hss(1,3)
admissible-patterns\hss}}}001$ $\displaystyle\text{Unoccupied }=N_{\phi}-N-4$
$\displaystyle\text{Occupied }=N-2$

Repeating similar calculations, as we carried out for the sphere, to count the
inner (1,3) admissible-configurations with fixed edge patterns, we obtain

$\displaystyle\begin{rcases}x_{0}+x_{1}+\dots+x_{N-2}=N_{\phi}-N-4\\\
x_{0}\geq 0\\\ x_{1}\geq 2\\\ x_{2}\geq 2\\\ \vdots\\\ x_{N-3}\geq 2\\\
x_{N-2}\geq 0\end{rcases}$ (8)

As we did in the sphere counting, we make a change of variables in (8) to
include the constraint inside the integer equation:
$a_{0}=x_{0},a_{N}=x_{N},a_{i}=x_{i}-2$, which yields

$\displaystyle\begin{rcases}a_{0}+a_{1}+\dots+a_{N-2}=N_{\phi}-3N+2\\\
a_{0}\geq 0\\\ a_{1}\geq 0\\\ a_{2}\geq 0\\\ \vdots\\\ a_{N-3}\geq 0\\\
a_{N-2}\geq 0\end{rcases}$ (9)

with the solution,

$\boxed{\mbox{Number of integer solutions}=\binom{N+n-2}{n}}$

(b) $\displaystyle 100\text{\text@underline{\hbox to170.71652pt{\hss(1,3)
admissible-patterns\hss}}}0010$ $\displaystyle\text{Unoccupied }=N_{\phi}-N-5$
$\displaystyle\text{Occupied }=N-2$

In this case, calculations are same as in the above case (a) except that the
number of unoccupied orbitals changes and the result is

$\boxed{\mbox{Number of integer solutions}=\binom{N+n-3}{n-1}}$

Due to inversion symmetry, 0100(1,3) admissible-patterns001 generates the same
number of configurations as (b). Then the total number of degenerate states
for the quasihole excitations of $\nu=1/3$ with TBCs is

$\mathcal{D_{1/3}}=\binom{N+n}{n}-\binom{N+n-2}{n}-2\binom{N+n-3}{n-1}.$

For a generic fraction $\nu=1/r$, we see by inspection that there exist $r-1$
independent boundary-condition violating terms which need to be multiplied by
appropriate constant factors due to inversion symmetry in order to find the
total number of disallowed configurations. The proposed final formula for the
degeneracy is

$\mathcal{D}_{1/r}=\binom{N+n}{n}-\binom{N+n-2}{n}-2\binom{N+n-3}{n-1}-\dots-(r-1)\binom{N+n-r}{n-r+2}\\\
=(rN+n-r+1)\frac{(N+n-r)!}{N!(n-r+1)!},$ (10)

where $n\equiv N_{\phi}-Nr+r-1$. We will prove the last equality by induction.
We obtained the right-hand side combinatorial form of the degeneracy by
guessing from the numerical degeneracies we encountered in both our
calculations and some other works, e.g. [20, 21, 22]. Later we have become
aware of Ref. [30], which contains essentially the same formula (in its
footnote 49) derived in the thin-torus limit.

Proof of Eq. (10): One can easily show that Eq. (10) is true for $r=2$.
Assuming that the equation holds true for $r=k$:

$\binom{N+n}{n}-\binom{N+n-2}{n}-\dots-(k-1)\binom{N+n-k}{n-k+2}=(kN+n-k+1)\frac{(N+n-k)!}{N!(n-k+1)!},$
(11)

we will show that it holds true for $r=k+1$:

$\binom{N+n}{n}-\binom{N+n-2}{n}-\dots-(k-1)\binom{N+n-k}{n-k+2}-k\binom{N+n-k-1}{n-k+1}\\\
=(kN+N+n-k)\frac{(N+n-k-1)!}{N!(n-k)!}.$ (12)

Inserting Eq. (11) in the left-hand side of Eq. (12) and rearranging terms
yield

$\displaystyle(kN+n-k+1)\frac{(N+n-k)!}{N!(n-k+1)!}-k\frac{(N+n-k-1)!}{(n-k+1)!(N-2)!}$
$\displaystyle=$
$\displaystyle(kN+n-k+1)\frac{(N+n-k)(N+n-k-1)!}{N!(n-k+1)!}-k\frac{N(N-1)(N+n-k-1)!}{N!(n-k+1)!}$
$\displaystyle=$
$\displaystyle\frac{(N+n-k-1)!}{N!(n-k)!}\Bigg{[}\frac{(kN+n-k+1)(N+n-k)-kN(N-1)}{n-k+1}\Bigg{]}$
$\displaystyle=$
$\displaystyle\frac{(N+n-k-1)!}{N!(n-k)!}\Bigg{\\{}\frac{(n-k+1)(N+n-k)+kN\big{[}N+n-k-(N-1)\big{]}}{n-k+1}\Bigg{\\}}$
$\displaystyle=$
$\displaystyle\frac{(N+n-k-1)!}{N!(n-k)!}\Bigg{[}\frac{\cancel{(n-k+1)}(N+n-k+kN)}{\cancel{(n-k+1)}}\Bigg{]}$
$\displaystyle=$ $\displaystyle(kN+N+n-k)\frac{(N+n-k-1)!}{N!(n-k)!}$
$\displaystyle=$ $\displaystyle\blacksquare$

Therefore, for a given $\nu=1/r$ we can obtain the number of degenerate
quasihole excitations, and we can have it in a form that only depends on
system parameters by using the definitions made previously. Since
$n=N_{\phi}-Nr+r-1$,

$\displaystyle\mathcal{D}_{1/r}$
$\displaystyle=(rN+n-r+1)\frac{(N+n-r)!}{N!(n-r+1)!}$
$\displaystyle=N_{\phi}\frac{(N+N_{\phi}-Nr-1)!}{N!(N_{\phi}-Nr)!}$
$\displaystyle\boxed{\mathcal{D}(N,N_{d},N_{\phi})=\frac{(N_{d}+N-1)!}{N_{d}!(N-1)!}\frac{N_{\phi}}{N}},$

where $N_{\phi}=Nr+N_{d}$, $N_{d}$ being the number of delocalized quasiholes
in the system.

##  Appendix B Energy gaps and the system size

In this appendix, we exemplify the dependence of the small energy gaps marking
different groups of states on the system size. This effect is already
implicitly contained in Fig. 3 of the main text, where for fixed interaction
strength and number of particles it is seen that the mini-gaps close as one
approaches the continuum limit of smaller flux quanta per unit cell. As we
mentioned below Eq. (3) of the main text, in our numerical modeling, we choose
to fit full magnetic unit cells along any direction of the lattice in order to
avoid possible boundary effects. As a result, while going to the continuum
limit by decreasing the flux quanta per unit cell $\phi=p/q$ through
increasing $q$ ($=L_{y})$, the lattice size increases automatically. This
relation also ensures that $L_{x}$ flux quanta are contained in the whole
lattice at all times. In Fig. 4, we show this effect for three nearby $q$
values (panel(b) has the same parameters as in Fig. 1) for a system of two
hardcore bosons. It is apparent that the energy gaps decrease with increasing
$q=L_{y}$.

Figure 4: Diminishing of the energy gaps as the system size increases.
Hardcore interactions are assumed for $N=2$ particles with lattice size
$L_{x}=10$, and (a) $L_{y}=14$, (b) $L_{y}=15$, (c) $L_{y}=16$. In each case,
flux is taken to be $\phi=1/q$ with $q=L_{y}$ and $E_{0}$ is the lowest
energy.

##  Appendix C Many-body Chern number calculation

Here, for completeness, we briefly outline the method of Refs. [40, 33, 41]
for the many-body Chern number calculation. In this method, the Berry
curvature is expressed in a discretized parameter space of twisted angles
$\vec{\theta}=(\theta_{x},\theta_{y})$ of the torus boundary conditions [45].
In our numerical calculations twisted angles are discretized as
$\theta_{i}=\frac{2\pi}{N_{\theta}}n_{i}$ ($n_{i}=1,\dots,N_{\theta};i=x,y$),
where $N_{\theta}$ tunes the fineness of the mesh. The discretized Berry
curvature is defined as
$\mathcal{F(\vec{\theta})}\\!=\\!\text{log}[U_{x}(\vec{\theta})U_{y}(\vec{\theta}+\delta_{x})U_{x}(\vec{\theta}+\delta_{y})^{-1}U_{y}(\vec{\theta})^{-1}]$,
where the link variables $U_{i}(\vec{\theta})$ are expressed as

$\displaystyle
U_{i}(\vec{\theta})=\frac{\text{det}[\Phi^{\dagger}(\vec{\theta})\Phi(\vec{\theta}+\delta_{i})]}{|{\text{det}[\Phi^{\dagger}(\vec{\theta})\Phi(\vec{\theta}+\delta_{i})]|}};\>\>i=x,y.$
(13)

Here, $\Phi(\vec{\theta})$ stands for the multiplet matrix
$\Phi(\vec{\theta})=(\ket{G_{1}(\vec{\theta})},\dots,\ket{G_{\mathcal{D}}(\vec{\theta})})$
formed by $\mathcal{D}$-fold degenerate states for which we want to calculate
the Chern number and $\delta_{x}=(\frac{2\pi}{N_{\theta}},0)$,
$\delta_{y}=(0,\frac{2\pi}{N_{\theta}})$ are the incremental units in
parameter space. Calculation of the link variable $U_{i}(\vec{\theta})$
involves finding the overlap between states at different twisted angles
corresponding to different boundary conditions. From this construct, the Chern
number can be calculated by taking the sum over all possible twisted angles in
the mesh:

$\displaystyle\mathcal{C}=\frac{1}{2\pi
i}\sum_{\vec{\theta}}\mathcal{F(\vec{\theta})}.$ (14)

## References

  * [1] K. v. Klitzing, G. Dorda, and M. Pepper, Phys. Rev. Lett. 45, 494 (1980).
  * [2] D. C. Tsui, H. L. Stormer, and A. C. Gossard, Phys. Rev. Lett. 48, 1559 (1982).
  * [3] M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * [4] D. J. Thouless, M. Kohmoto, M. P. Nightingale, and M. den Nijs, Phys. Rev. Lett. 49, 405 (1982).
  * [5] B. I. Halperin, Phys. Rev. B 25, 2185 (1982).
  * [6] R. B. Laughlin, Phys. Rev. Lett. 50, 1395 (1983).
  * [7] J. M. Leinaas and J. Myrheim, Nuovo Cim B 37, 1 (1977).
  * [8] D. Arovas, J. R. Schrieffer, and F. Wilczek, Phys. Rev. Lett. 53, 722 (1984).
  * [9] F. D. M. Haldane, Phys. Rev. Lett. 67, 937 (1991).
  * [10] C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S. Das Sarma, Rev. Mod. Phys. 80, 1083 (2008).
  * [11] Y. Ando, J. Phys. Soc. Jpn. 82, 102001 (2013).
  * [12] I. Carusotto, A. A. Houck, A. J. Kollár, P. Roushan, D. I. Schuster, and J. Simon, Nat. Phys. 16, 268 (2020).
  * [13] I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008)
  * [14] I. Carusotto and C. Ciuti, Rev. Mod. Phys. 85, 299 (2013).
  * [15] N. R. Cooper, N. K. Wilkin, and J. M. F. Gunn, Phys. Rev. Lett. 87, 120405 (2001).
  * [16] B. Paredes, P. Fedichev, J. I. Cirac, and P. Zoller, Phys. Rev. Lett. 87, 010402 (2001).
  * [17] A. S. Sørensen, E. Demler, and M. D. Lukin, Phys. Rev. Lett. 94, 086803 (2005).
  * [18] M. Hafezi, A. S. Sørensen, E. Demler, and M. D. Lukin, Phys. Rev. A 76, 023613 (2007).
  * [19] J. Cho, D. G. Angelakis, and S. Bose, Phys. Rev. Lett. 101, 246809 (2008).
  * [20] N. Regnault and B. A. Bernevig, Phys. Rev. X 1, 021014 (2011).
  * [21] Z. Liu, E. J. Bergholtz, H. Fan, and A. M. Läuchli, Phys. Rev. Lett. 109, 186805 (2012).
  * [22] B. A. Bernevig and N. Regnault, Phys. Rev. B 85, 075128 (2012).
  * [23] Z. Liu, E. J. Bergholtz, and E. Kapit, Phys. Rev. B 88, 205101 (2013).
  * [24] S. A. Parameswaran, R. Roy, and S. L. Sondhi, C R Phys, 14 816 (2013).
  * [25] I. D. Rodríguez and Anne E. B. Nielsen, Phys. Rev. B 92, 125105 (2015).
  * [26] M. Račiūnas, F. N. Ünal, E. Anisimovas, and A. Eckardt, Phys. Rev. A 98, 063621 (2018).
  * [27] Z. Liu, R. N. Bhatt, and N. Regnault, Phys. Rev. B 91, 045126 (2015).
  * [28] E. Macaluso, T. Comparin, R. O. Umucalılar, M. Gerster, S. Montangero, M. Rizzi, and I. Carusotto, Phys. Rev. Research 2, 013145 (2020).
  * [29] B. A. Bernevig and F. D. M. Haldane, Phys. Rev. Lett. 100, 246802 (2008).
  * [30] A. Seidel and K. Yang, Phys. Rev. B 84, 085122 (2011).
  * [31] E. Kapit, P. Ginsparg, and E. Mueller, Phys. Rev. Lett. 108, 066802 (2012).
  * [32] X. G. Wen and Q. Niu, Phys. Rev. B 41, 9377 (1990).
  * [33] K. Kudo, T. Kariyado, and Y. Hatsugai, J. Phys. Soc. Jpn. 86, 103701 (2017).
  * [34] F. D. M. Haldane, Phys. Rev. Lett. 55, 2095 (1985).
  * [35] D. Langbein, Phys. Rev. 180, 633 (1969).
  * [36] D. R. Hofstadter, Phys. Rev. B 14, 2239 (1976).
  * [37] T. S. Jackson, D. Bauer, and R. Roy, Phys. Rev. B 93, 235133 (2016).
  * [38] B. Andrews and G. Möller, Phys. Rev. B 97, 035159 (2018).
  * [39] R. O. Umucalılar, Phys. Rev. A 98, 063629 (2018).
  * [40] T. Fukui, Y. Hatsugai, and H. Suzuki, J. Phys. Soc. Jpn. 74, 1674 (2005).
  * [41] K. Kudo, H. Watanabe, T. Kariyado, and Y. Hatsugai, Phys. Rev. Lett. 122, 146601 (2019).
  * [42] N. S. Srivatsa, X. Li, and A. E. B. Nielsen, Phys. Rev. Research 3, 033044 (2021).
  * [43] E. Kapit and E. Mueller, Phys. Rev. Lett. 105, 215303 (2010).
  * [44] H. Atakişi and M. Ö. Oktel, Phys. Rev. A 88, 033612 (2013).
  * [45] Q. Niu, D. J. Thouless, and Y.-S. Wu, Phys. Rev. B 31, 3372 (1985).