# Temporal Graph Rewiring with Expander Graphs

Katarina Petrović  Shenyang Huang  Farimah Poursafaei  Petar Veličković

###### Abstract

Evolving relations in real-world networks are often modelled by temporal
graphs. Graph rewiring techniques have been utilised on Graph Neural Networks
(GNNs) to improve expressiveness and increase model performance. In this work,
we propose Temporal Graph Rewiring (TGR), the first approach for graph
rewiring on temporal graphs. TGR enables communication between temporally
distant nodes in a continuous time dynamic graph by utilising expander graph
propagation to construct a message passing highway for message passing between
distant nodes. Expander graphs are suitable candidates for rewiring as they
help overcome the oversquashing problem often observed in GNNs. On the public
tgbl-wiki benchmark, we show that TGR improves the performance of a widely
used TGN model by a significant margin. Our code repository is accessible at:
https://github.com/kpetrovicc/TGR.git.

Temporal Graphs, Dynamic Graphs, Graph Neural Network, Graph Representation
Learning

##  1 Introduction

Graph representation learning (Hamilton, 2020) aims to learn node
representations on graph structured data to solve tasks such as node property
prediction (Hamilton et al., 2017), link prediction (Ying et al., 2018) and
graph property prediction (Gilmer et al., 2017). Graph neural networks (GNN)
(Kipf & Welling, 2016; Veličković et al., 2017; Xu et al., 2018) capture
relationships between nodes following message passing paradigm (Veličković,
2022). GNNs have been successfully applied to model many real-world networks
such as biological networks (Johnson et al., 2023; Zitnik et al., 2018) and
social networks (Ying et al., 2018).

Oversquashing in GNNs. GNNs operate through message passing mechanism which
aggregates information over a node’s direct neighbourhood at each GNN layer.
Hence, propagating information between nodes that are at $k$-hop distance
requires $k$ GNN layers. This leads to information bottleneck, since the
node’s receptive field in most underlying graph topologies grows exponentially
with an increase of number of GNN layers. This information is further stored
in fixed length embedding vectors, leading to a phenomena called
oversquashing, causing a severe loss of information. In static GNNs,
bottlenecks and oversquashing are addressed using graph rewiring, where
underlying graph topology is altered such that it connects distant nodes.
There are several methods applied to rewire static GNNs such as diffusion-
based graph rewiring (Gasteiger et al., 2019) or reducing negative Ricci
curvature (Topping et al., 2021), leading to an impressive boost in
performance.

Temporal graph learning (TGL) (Kazemi et al., 2020) emerged from graph
representation learning to study evolving relationships in temporal graph
data. Given the dynamic nature of temporal graphs, temporal graph neural
networks (TGNNs) (Longa et al., 2023) are developed to capture the evolution
of graph topology introducing novel components such as temporal memory (Kazemi
et al., 2020) and time-encoding (Xu et al., 2020). TGNNs are underpinned by a
static GNN message passing mechanism, meaning they are vulnerable to
bottlenecks and oversquashing. Given additional temporal dimension, receptive
field of a node in a temporal graph grows faster than in static graphs,
leading to a higher presence of oversquashing and bottlenecks in the
structure. A natural question arises: Can we improve TGNN performance by
introducing temporal graph rewiring aimed to relieve oversquashing and ensure
global propagation between temporally distant nodes?

##### Why rewiring on temporal graphs?

Rewiring fundamentally calls for a dynamic nature of how the data exists in
the real world. There might be a process different from the input graph
structure that governs the message passing and information exchange in a
temporal graph. Additionally, temporal graph rewiring can be seen as a way to
resolve so-called memory staleness problem (Kazemi et al., 2020). By memory
staleness we refer to a process occurring in TGNN temporal memory. Temporal
memory is only updated if a node of interest interacts with another node in a
graph, causing inactive nodes to appear as stale. This set-up poses
significant issues in e.g social networks where if a user or a node is
inactive for a period of time, it loses connections to its active neighbours.
We propose temporal graph rewiring to add connections between temporally
distant nodes, allowing information flow between nodes that might be inactive
for a longer period of time.

Figure 1: Node Bank stores all nodes observed up until batch $t$. In figure,
timestamp $t^{-}$ represents timestamp before $t$. Temporal memory stores
observed node states mixed through expander message passing layer to generate
expander embeddings $\mathbf{H}(t^{-})$. TGNN input features are constructed
by combining expander embeddings $\mathbf{H}(t^{-})$ for previously observed
nodes and temporal memory node states $\mathbf{S}(t^{-})$ for unobserved
nodes. The green node is a node of interest shown alongside its 1-hop
neighbourhood.

In this work, we propose a novel TGR framework for Temporal Graph Rewiring
with expander graphs. We leverage recent work by Deac et al. (2022) on
expander graph propagation in combination with a TGNN base model to ensure
global message passing between temporally distant nodes. It is shown that
expander graphs satisfy four desirable criteria for graph rewiring: 1). the
ability to propagate information globally within the input graph, 2).
relieving bottlenecks and oversquashing, 3). subquadratic space and time
complexity and 4). no additional pre-processing of the input graph (Deac et
al., 2022).

##### Main Paper Contributions.

Our main paper contributions include:

  * •

First Rewiring Method on Temporal Graphs. In this work, we propose _Temporal
Graph Rewiring_ or TGR. To the best of our knowledge, TGR is the first
approach which applies graph rewiring techniques on temporal graphs. TGR
combines expander graph propagation (Deac et al., 2022) to rewire a base TGNN
such as TGN (Rossi et al., 2020), while operating with minimal computational
overhead.

  * •

Model-agnostic. TGR is agnostic to the chosen expander message passing layer
and helps alleviate the _over-squashing_ and _memory staleness_ problems in
TGNNs. TGR is also agnostic to the base TGNN model as TGR provides dynamic
node features thus has the potential to be applied to a wide range of TGL
models.

  * •

Improved Performance. We test TGR on a temporal link prediction task using
tgbl-wiki data taken from a publicly available Temporal Graph Benchmark (TGB)
(Huang et al., 2023). In section 5 we show that using TGR across four
different expander message passing layers outperforms the base TGNN model.

##  2 Related Work

In this section we discuss some related work on TGNNs and static graph
rewiring, highlighting their key advantages and disadvantages in relation to
temporal graph rewiring.

###  2.1 Temporal Graph Neural Networks

TGL has been a popular topic in graph machine learning community focusing on
modelling evolutionary behaviour of complex dynamic interaction networks and
leading to development of temporal graph neural network (TGNN) architectures
(Longa et al., 2023). Following terminology introduced in (Longa et al., 2023)
TGNNs can be classified into snapshot-based architectures operating over
discrete time dynamic graphs (DTDGs) and event-based architectures operating
over continuous time dynamic graphs (CTDGs). CTDGs capture events that can
happen at any time, while DTDGs capture events that happen in discrete time
intervals. In this work we focus on CTDGs, as most real-world networks operate
over a continuous time domain.

Event-based TGNNs such as Neighborhood-aware Scalable Temporal Network (NAT)
(Johnson et al., 2023) or Temporal Graph Network (TGN) (Rossi et al., 2020)
are characterised by use of graph-based operators in combination with novel
model components such as temporal memory (Kazemi et al., 2020) and time-
encoding (Xu et al., 2020). TGNNs use temporal memory to store long-standing
relationships between observed node states used to generate dynamic node
embeddings. Similarly to static GNNs, TGNNs operate through aggregating
information from a local temporal neighbourhood to generate temporal
embeddings following a message passing paradigm. This leads to their
vulnerability to bottlenecks and oversquashing.

###  2.2 Static Graph Rewiring

In static GNNs, bottlenecks and oversquashing have been addressed with graph
rewiring leading to an impressive boost in performance. There has been an
extensive work done on applying graph rewiring to static graphs. For example
diffusion-based graph rewiring (Gasteiger et al., 2019) diffuses additional
edges in the graph with use of kernels such as PageRank (Brin, 1998). However
as stated in (Topping et al., 2021), these models generally fail to reduce
bottlenecks in the input graph. Topping et al. (2021) introduces a method to
modify a portion of edges with negative Ricci curvature to reduce
oversquashing. However, this is shown to come with higher pre-processing cost
and it is sub-optimal for analysing large graphs such as continuous time
dynamic graphs which are the focus of our paper.

##### Rewiring with Expander Graphs.

In this work we opt to work with expander graph propagation (Deac et al.,
2022) as a starting point. Unlike other rewiring methods, expander graph
propagation involves propagating messages over a graph that is independent of
the input graph topology, by replacing input graph with an expander graph. In
section 3.4, we describe in detail set-up by Deac et al. (2022), where
expander and input graph GNN layers are alternated to compute node embeddings.
In recent work by Giovanni et al. (2024), it is formally validated that such
approach is effective in relieving oversquashing by decreasing overall
effective resistance or commute time (Chandra et al., 1989). In fact, commute
time is shown to be closely aligned with oversquashing, and it is stated that
commute time in expander graphs grows linearly with the number of edges. In
section 3.3, we provide further information on expander graph properties,
showing they have high spectral gap while maintaining high sparsity making
them an ideal candidate for graph rewiring. We are aware that expander graph
propagation may not be the only way to do temporal graph rewiring. In this
work, we aim to bring perspective of using expander graph propagation in
combination with a temporal graph learning model and create pathway for
further avenues to be explored which we leave for future work.

##  3 Preliminary Background

This section provides notation and overview of interesting expander graph
topological properties. It discusses expander graph propagation (Deac et al.,
2022), a framework for static graph rewiring using expander graphs, and
temporal graph network (Rossi et al., 2020), a widely used TGNN model. We
combine these models later in section 4 where we describe TGR method.

###  3.1 Temporal Graphs

Temporal graphs are characterised by an evolving underlying graph topology,
thus making them suitable for modelling complex dynamic networks and evolving
interactions between nodes. In this work we focus on continuous time dynamic
graphs as they model events which can occur at any point in time, thus making
them suitable for modelling a majority of real-world networks. Continuous time
dynamic graphs (CTDGs) are defined as a sequence of chronologically sorted
events that can be split into node-wise or edge-wise events. Node-wise events
$\mathbf{v}_{i}(t)$ include addition and deletion of node $i$ as well as
change in node features at time $t$. Edge-wise events $\mathbf{e}_{i,j}(t)$
relate to an edge $ij$ with endpoint nodes $i$ and $j$ and include formation
and deletion of edges in a graph at time $t$.

###  3.2 Temporal Graph Network

Temporal Graph Network (TGN) (Rossi et al., 2020) serves as a base model to
TGR framework to generate temporal embeddings $\mathbf{Z}(t)$. TGN proposes a
combination of temporal memory and graph-based operators for learning on
continuous time dynamic graphs represented as chronologically sorted sequences
of time-based events.

A unique advantage of TGN is use of temporal memory to store and update node
states. TGN memory module $\mathbf{S}(t)$ stores memory state
$\mathbf{s}_{i}(t)$ for each node $i$ in the graph. When a node interacts with
another node through a node-wise or edge-wise event, node states of the nodes
involved in the event are updated in the memory through a recurrent neural
network such as LSTM (Hochreiter & Schmidhuber, 1997) or GRU (Cho et al.,
2014). The purpose of a memory module is to store historical information and
long-term dependencies of a node in a compressed format.

###  3.3 Expander Graphs

Following terminology laid out in (Deac et al., 2022), we provide a non-
exhaustive list of most important expander graph topological properties in
relation to temporal graph rewiring. Expander graphs are a fundamentally
sparse family of graphs with number of edges scaling linearly with number of
nodes ($|E|=O(n)$). They are characterised by two main properties: (P1) they
have no bottlenecks and maintain high connectivity or spectral gap and (P2)
they are efficiently precomputable through use of group operators and act as
an independent graph topology over which input graph is rewired. Furthermore,
their commute time or effective resistance scales linearly with the number of
edges as stated in (Giovanni et al., 2024) showing that expander graphs are
effective in reducing oversquashing.

##### Expander graphs have no bottlenecks.

Spectral gap of the graph $G$ with $n$ nodes is related to the connectivity
and existence of bottlenecks within its structure. We observe the eigenvalues
of a graph-Laplacian $L$ given as:

$\lambda_{0}<\lambda_{1}\leq\lambda_{2}\leq...\leq\lambda_{n},$ (1)

where the first eigenvalue $\lambda_{1}$ is referred to as spectral gap or
measure of connectivity. From Cheeger constant, it follows that a larger
$\lambda_{1}$ leads to a better connectivity within graph $G$. Expander graph
eigenvalues are extensively discussed in (Alon & Milman, 1985; Alon, 1986;
Dodziuk, 1984; Tanner, 1984). From (Deac et al., 2022) it is known that
Cheeger constant of expander graph families $\\{G_{i}\\}$ is lower-bounded by
a positive constant $\epsilon>0$, leading to high connectivity.

##### Expander graphs are efficiently precomputable.

We follow construction given in (Deac et al., 2022) and we consider Cayley
graph expander family. Cayley graphs are constructed through a use of the
special linear group $\text{SL}(2,\mathbb{Z}_{n})$ as a generating set. Full
details on pre-computing Cayley graphs can be found in the work done by (Deac
et al., 2022; Kowalski, 2019; Selberg, 1965; Davidoff et al., 2003). It is
known that a lower bound on the first eigenvalue of Cayley graph Laplacian is
at least $3/16$ with full proofs given in (Kowalski, 2019; Davidoff et al.,
2003). This means that Cayley graphs do not exhibit bottlenecks in their
structure (satisfying P1).

Computing Cayley graphs comes with a low computational cost and no additional
pre-processing of the input. It is shown in (Deac et al., 2022) that for an
expander family $\\{G_{i}\\}$ of finite graphs with uniform upper bound on
their vertex degree, the following inequality holds:

$\text{diam}(G_{i})\leq k\log V(G_{i}),$ (2)

where constant $k>0$ and $\text{diam}(G_{i})$ is the diameter of graph
$G_{i}$, yielding subquadratic time complexity. Due to their low-diameter, two
expander nodes can reach each other within a small number of hops, relieving
bottlenecks and offering efficient structure for global information
propagation.

It is also important to note that the number of nodes in Cayley graphs grows
cubically. As shown in (Deac et al., 2022), the number of nodes is computed
as:

$|V(\text{Cay}(SL(2,\mathbb{Z}_{n});S_{n}))|=n^{3}\prod_{\begin{subarray}{c}\text{prime
}p|n\end{subarray}}\left(1-\frac{1}{p^{2}}\right).$ (3)

###  3.4 Expander Graph Propagation

Expander Graph Propagation (EGP) framework is simplistic in nature and shows
incredibly favourable properties for rewiring on static graphs: relieving
information bottleneck and oversquashing, while maintaining low computational
cost and subquadratic time. For a set of input node features
$\mathbf{X}^{n\times d}$, where $n$ is the number of nodes and $d$ is the
feature vector size, EGP operates by alternating GNN layers and propagating
them over the adjacency matrix of the input graph and Cayley graph in a
following manner:

$\mathbf{H}=\text{GNN}(\text{GNN}(\mathbf{X},\mathbf{A}),\mathbf{A}^{Cay}).$
(4)

Here, $\mathbf{H}$ represents node embedding vector generated after EGP
forward pass, and $\mathbf{A}$ and $\mathbf{A}^{Cay}$ are the adjacency
matrices for the input and Cayley graph respectively. A GNN can be any
classical GNN layer, such as graph attentional network (GAT) (Veličković et
al., 2017):

$\mathbf{h}_{i}=\underset{k=1}{\overset{K}{{\bigg{\|}}}}\sigma\left(\sum_{j\in\mathcal{N}_{i}}\alpha_{ij}^{k}\mathbf{W}^{k}\mathbf{x}_{j}\right).$
(5)

In the equation above for multi-head-attention layer, $i$ is a node of
interest and $\mathcal{N}_{i}$ is its direct neighbourhood, $\alpha_{ij}^{k}$
are the attention weights of the $k$-th attention head and $\mathbf{W}^{k}$ is
corresponding input linear transformation’s weight matrix.

##  4 Temporal Graph Rewiring

TGR combines expander graph propagation with use of memory module to reduce
memory staleness and alleviate oversquashing in a base TGNN model. A detailed
mechanism behind our framework is given in Figure 2 and algorithm in appendix
A.1. During training, TGR stores observed nodes in Node Bank. After TGNN
forward pass, TGR extracts observed node states from TGNN temporal memory and
passes them to the expander message passing layer, which can be any classic
static GNN layer. This process is called memory mixing and it is used to
reduce memory staleness. Expander graphs introduce edges between temporally
distant nodes in TGNN temporal memory, regardless whether they interacted with
other nodes in the batch. Embeddings generated through memory mixing are then
passed as input features to the base TGNN model.

As TGR essentially provides dynamic node features, it is agnostic to the
underlying TGNN architecture thus it can be easily extended to other TGL
models. In this work, we demonstrate how TGR integrates with a widely used TGN
architecture. It is important to note that in our method, TGNN construction
remains intact. We add memory mixing as a separate process to TGNN after TGNN
temporal memory has been updated. This a difference to most static graph
rewiring set-ups, where rewiring alters individual GNN layers in the
architecture.

##### Node Bank.

Node Bank stores all nodes observed up to a temporal batch $t$. We define two
classes of nodes in the batch: previously observed and unobserved nodes. By
previously observed nodes we refer to nodes that have been observed by the
TGNN model in previous batches. By an unobserved node, we refer to nodes that
have not yet been seen by the TGNN model. Node Bank module is dynamically
updated with new nodes at every temporal batch.

Figure 2: TGR batch processing. Green nodes are nodes that have been
previously observed and blue nodes are the new nodes. Input node features
$\mathbf{X}(t)$ are constructed by concatenating expander embeddings
$\mathbf{H}(t^{-})$ (for previously observed nodes), and node states
$\mathbf{S}(t^{-})$ (for new nodes). After computing temporal embeddings
$\mathbf{Z}(t)$, temporal memory is updated and mixed with other previously
observed node states to compute expander embeddings $\mathbf{H^{\prime}}(t)$
and update expander memory.

##### Memory Mixing.

TGR introduces expander graph propagation through memory mixing to reduce
memory staleness. After TGNN forward pass, TGNN temporal memory module is
updated with new node states. Consider $\mathbf{S^{\prime}}(t)$ which consists
of updated node states for all nodes in the Node Bank. Expander embeddings
$\mathbf{H^{\prime}}(t)$ are generated through an expander message passing
layer:

$\mathbf{H^{\prime}}(t)=\text{GNN}(\mathbf{S^{\prime}}(t),\mathbf{A}^{Cay}).$
(6)

Here, $\mathbf{A}^{Cay}$ is the expander graph adjacency matrix and GNN can be
any classical GNN layer, such as graph attentional network (GAT) (Veličković
et al., 2017):

$\mathbf{h}_{i}=\underset{k=1}{\overset{K}{{\bigg{\|}}}}\sigma\left(\sum_{j\in\mathcal{N}_{i}}\alpha_{ij}^{k}\mathbf{W}^{k}\mathbf{x}_{j}\right).$
(7)

In the equation above for multi-head-attention layer, $i$ is a node of
interest and $\mathcal{N}_{i}$ is its direct neighbourhood, $\alpha_{ij}^{k}$
are the attention weights of the $k$-th attention head and $\mathbf{W}^{k}$ is
the corresponding input linear transformation’s weight matrix. In our set-up,
we pre-compute a large expander graph from Cayley graph family with $N$ nodes,
where $N$ is the number of nodes in the training dataset. We use entire
expander graph for passing messages in the expander message passing layer.
When expander graph is smaller than number of nodes in the node bank, we
compute a new expander graph of size $(1+m)N$, where $m<1$.

##### Expander Memory.

Expander memory $\mathbf{H^{\prime}}(t)$ stores computed expander embeddings
and provides input features for previously observed nodes for TGNN forward
pass.

##### Input Features for TGNN.

The input feature vector $\mathbf{X}(t)$ for TGNN forward pass is constructed
through a concatenation ($\|$) of expander embeddings $\mathbf{H}(t^{-})$ (for
previously observed nodes) from expander memory and TGNN node states
$\mathbf{S}(t^{-})$ (for new nodes) from temporal memory:

$\mathbf{X}(t)=\mathbf{H}(t^{-})\|\mathbf{S}(t^{-}).$ (8)

##### TGNN Embedding.

TGNN embedding module generates temporal embeddings $\mathbf{Z}(t)$ by
aggregating information from node’s direct temporal neighbourhood:

$\mathbf{Z}(t)=\text{TGNN}(\mathbf{X}(t),\mathbf{A}(t)).$ (9)

We define temporal neighbourhood $N_{i}^{t}$, of node $i$ at time $t$ as a set
of temporal events where every edge $(i,j)\in N_{i}^{t}$. In this work we opt
for a widely used TGN (Rossi et al., 2020) for a choice of a TGNN embedding.
For node $i$ its TGN embedding $\mathbf{z}_{i}(t)$ is computed as:

$\displaystyle\mathbf{z}_{i}(t)=\text{MLP}(\mathbf{x}_{i}(t)\|\tilde{\mathbf{z}}_{i}(t)),$
(10)
$\displaystyle\tilde{\mathbf{z}}_{i}(t)=\text{MultiHeadAttention}(\mathbf{q}(t),\mathbf{K}(t),\mathbf{V}(t)),$
(11) $\displaystyle\mathbf{q}(t)=\mathbf{x}_{i}(t)\|\phi(0),$ (12)
$\displaystyle\mathbf{K}(t)=\mathbf{V}(t)=\mathbf{C}(t),$ (13)
$\displaystyle\mathbf{C}(t)=[\mathbf{x}_{1}\|\mathbf{e}_{i1}\|\phi(t-t_{1}),...,\mathbf{x}_{N}\|\mathbf{e}_{iN}\|\phi(t-t_{N})]$
(14)

MultiHeadAttention represents multi-head-attention layer (Vaswani et al.,
2023), and $\phi$ is the generic time encoding (Xu et al., 2020). Query
$\mathbf{q}(t)$ relates to node of interest and the keys $\mathbf{K}(t)$ and
values $\mathbf{V}(t)$ relate to its direct neighbours.

Table 1: Results for _dynamic link property prediction_ on tgbl-wiki.

Method | MRR (%) | Improvement (%) |  Compute Time (Sec)
---|---|---|---
Val. | Test | Val. | Test | Val. | Test
TGN | 51.5$\pm$2.9 | 46.8$\pm$2.2 | - | - | 87 | 86
TGR-GCN | 56.5$\pm$1.2 | 50.3$\pm$1.2 | 5.0 | 3.6 | 92 | 94
TGR-GAT | 57.7$\pm$1.2 | 53.1$\pm$1.1 | 6.3 | 6.3 | 95 | 99
TGR-GIN | 55.9$\pm$1.5 | 51.3$\pm$1.9 | 4.5 | 4.5 | 93 | 94
TGR-GATv2 | 54.6$\pm$0.9 | 46.7$\pm$3.4 | 3.2 | -0.1 | 96 | 97

##  5 Experiments

##### Evaluation Setting.

In this section, we evaluate TGR framework on the dynamic link property
prediction task from TGB (Huang et al., 2023) on the tgbl-wiki dataset.
Predicting links in continuous time dynamic graphs consists of computing
temporal node embeddings and predicting a probability of an edge formation and
associated edge features. We compare TGR performance to that of the TGN base
model using four different baselines: _TGR-GCN_ (Kipf & Welling, 2016), _TGR-
GAT_ (Veličković et al., 2017), _TGR-GIN_ (Xu et al., 2018) and _TGR-GATv2_
(Brody et al., 2022).

##### TGR Implementation.

We implement TGR on top of the TGN base model used in TGB (Huang et al., 2023)
to include memory mixing, expander memory and node bank module. Expander
memory is initialised as an empty tensor which matches temporal memory size in
(Rossi et al., 2020) set-up. We also initialise node bank and expander graph
of size equal to the size of training dataset. In all experiments we set
learning rate $lr=5\cdot 10^{-4}$ and run models with a tolerance of 50 epochs
for a maximum of 100 epochs during training, test and validation. We maintain
same batch size as in (Rossi et al., 2020) which contains 200 temporal events.

##### TGB Benchmark.

We leverage availability of temporal graph data in the TGB Benchmark (Huang et
al., 2023) to validate TGR performance. TGB Benchmark collects a variety of
real-world dynamic datasets which contain temporal interactions in many real-
world networks such as flights, transactions and beyond. In this work we
explore tgbl-wiki dataset which stores dynamic information about a co-editing
network on Wikipedia pages over a span of one month. The data is stored in a
bi-partite temporal graph where nodes represent either editors or wiki-pages
they interact with. An edge represents an action user takes when editing a
Wikipedia page and edge features contain textual information about a page of
interest. The goal is to predict existence and nature of links between editors
and Wikipedia pages at a future timestamp.

##### Results.

We show that TGR significantly outperforms TGN on a temporal link prediction
task. The best performing result is achieved using TGR-GAT showing $6\%$
improvement across Validation and Test MRR comparison to TGN trained with our
hyperparameters. Full table of results is given in table 1. Best results using
GAT reiterate importance of prioritising most important neighbours in random
selection of expander edges between observed nodes.

##  6 Conclusion

In this work, we proposed TGR to address gaps in TGNNs: bottlenecks,
oversquashing and memory staleness. We show that using expander graphs to
rewire temporal graphs is an optimal method to relieve oversquahsing while
adding minimal cost overhead.

##### Limitations and Future Work.

We are aware that expander graph propagation may not be the best way of doing
temporal graph rewiring. The purpose of this work is to conceptually propose
temporal graph rewiring as a solution to oversquashing using method we know is
efficient in relieving oversquashing in static GNNs. In future work, we hope
to explore other graph rewiring methods and compare their performance to TGR.
We have empirically tested TGR on tgbl-wiki dataset and shown that TGR
significantly outperforms base TGN model. In future work, we intend to test
TGR on tasks such as node and and graph property prediction and larger
datasets available on TGB Benchmark.

## Acknowledgements

We thank João Araújo and Simon Osindero for their comments and for reviewing
the paper prior to submission. We thank the TGB team for creating the platform
which allowed us to rapidly test and benchmark our method. We also thank
Emanuele Rossi, author of TGN, for his useful comments and feedback while we
were developing TGR.

## References

  * Alon (1986) Alon, N.  Eigenvalues and expanders.  _Combinatorica_ , 6(2):83–96, 1986.
  * Alon & Milman (1985) Alon, N. and Milman, V. D.  $\lambda$1, isoperimetric inequalities for graphs, and superconcentrators.  _Journal of Combinatorial Theory, Series B_ , 38(1):73–88, 1985.
  * Brin (1998) Brin, S.  The pagerank citation ranking: bringing order to the web.  _Proceedings of ASIS, 1998_ , 98:161–172, 1998.
  * Brody et al. (2022) Brody, S., Alon, U., and Yahav, E.  How attentive are graph attention networks?, 2022.
  * Chandra et al. (1989) Chandra, A. K., Raghavan, P., Ruzzo, W. L., and Smolensky, R.  The electrical resistance of a graph captures its commute and cover times.  In _Proceedings of the twenty-first annual ACM symposium on Theory of computing_ , pp. 574–586, 1989.
  * Cho et al. (2014) Cho, K., Van Merriënboer, B., Gulcehre, C., Bahdanau, D., Bougares, F., Schwenk, H., and Bengio, Y.  Learning phrase representations using rnn encoder-decoder for statistical machine translation.  _arXiv preprint arXiv:1406.1078_ , 2014.
  * Davidoff et al. (2003) Davidoff, G. P., Sarnak, P., and Valette, A.  _Elementary number theory, group theory, and Ramanujan graphs_ , volume 55.  Cambridge university press Cambridge, 2003.
  * Deac et al. (2022) Deac, A., Lackenby, M., and Veličković, P.  Expander graph propagation.  In Rieck, B. and Pascanu, R. (eds.), _Proceedings of the First Learning on Graphs Conference_ , volume 198 of _Proceedings of Machine Learning Research_ , pp. 38:1–38:18. PMLR, 09–12 Dec 2022.  URL https://proceedings.mlr.press/v198/deac22a.html.
  * Dodziuk (1984) Dodziuk, J.  Difference equations, isoperimetric inequality and transience of certain random walks.  _Transactions of the American Mathematical Society_ , 284(2):787–794, 1984.
  * Gasteiger et al. (2019) Gasteiger, J., Weißenberger, S., and Günnemann, S.  Diffusion improves graph learning.  _Advances in neural information processing systems_ , 32, 2019.
  * Gilmer et al. (2017) Gilmer, J., Schoenholz, S. S., Riley, P. F., Vinyals, O., and Dahl, G. E.  Neural message passing for quantum chemistry.  In _International conference on machine learning_ , pp. 1263–1272. PMLR, 2017.
  * Giovanni et al. (2024) Giovanni, F. D., Rusch, T. K., Bronstein, M. M., Deac, A., Lackenby, M., Mishra, S., and Veličković, P.  How does over-squashing affect the power of gnns?, 2024.
  * Hamilton et al. (2017) Hamilton, W., Ying, Z., and Leskovec, J.  Inductive representation learning on large graphs.  _Advances in neural information processing systems_ , 30, 2017.
  * Hamilton (2020) Hamilton, W. L.  _Graph representation learning_.  Morgan & Claypool Publishers, 2020.
  * Hochreiter & Schmidhuber (1997) Hochreiter, S. and Schmidhuber, J.  Long Short-Term Memory.  _Neural Computation_ , 9(8):1735–1780, 11 1997.  ISSN 0899-7667.  doi: 10.1162/neco.1997.9.8.1735.  URL https://doi.org/10.1162/neco.1997.9.8.1735.
  * Huang et al. (2023) Huang, S., Poursafaei, F., Danovitch, J., Fey, M., Hu, W., Rossi, E., Leskovec, J., Bronstein, M., Rabusseau, G., and Rabbany, R.  Temporal graph benchmark for machine learning on temporal graphs, 2023.
  * Johnson et al. (2023) Johnson, R., Li, M. M., Noori, A., Queen, O., and Zitnik, M.  Graph ai in medicine.  _arXiv preprint arXiv:2310.13767_ , 2023.
  * Kazemi et al. (2020) Kazemi, S. M., Goel, R., Jain, K., Kobyzev, I., Sethi, A., Forsyth, P., and Poupart, P.  Representation learning for dynamic graphs: A survey.  _Journal of Machine Learning Research_ , 21(70):1–73, 2020.
  * Kipf & Welling (2016) Kipf, T. N. and Welling, M.  Semi-supervised classification with graph convolutional networks.  _arXiv preprint arXiv:1609.02907_ , 2016.
  * Kowalski (2019) Kowalski, E.  _An introduction to expander graphs_.  Société mathématique de France Paris, 2019.
  * Longa et al. (2023) Longa, A., Lachi, V., Santin, G., Bianchini, M., Lepri, B., Lio, P., Scarselli, F., and Passerini, A.  Graph neural networks for temporal graphs: State of the art, open challenges, and opportunities, 2023.
  * Rossi et al. (2020) Rossi, E., Chamberlain, B., Frasca, F., Eynard, D., Monti, F., and Bronstein, M.  Temporal graph networks for deep learning on dynamic graphs, 2020.
  * Selberg (1965) Selberg, A.  On the estimation of fourier coefficients of modular forms.  In _Proceedings of Symposia in Pure Mathematics_ , pp. 1–15. American Mathematical Society, 1965.
  * Tanner (1984) Tanner, R. M.  Explicit concentrators from generalized n-gons.  _SIAM Journal on Algebraic Discrete Methods_ , 5(3):287–293, 1984.
  * Topping et al. (2021) Topping, J., Di Giovanni, F., Chamberlain, B. P., Dong, X., and Bronstein, M. M.  Understanding over-squashing and bottlenecks on graphs via curvature.  _arXiv preprint arXiv:2111.14522_ , 2021.
  * Vaswani et al. (2023) Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A. N., Kaiser, L., and Polosukhin, I.  Attention is all you need, 2023.
  * Veličković et al. (2017) Veličković, P., Cucurull, G., Casanova, A., Romero, A., Lio, P., and Bengio, Y.  Graph attention networks.  _arXiv preprint arXiv:1710.10903_ , 2017.
  * Veličković (2022) Veličković, P.  Message passing all the way up, 2022.
  * Xu et al. (2020) Xu, D., Ruan, C., Korpeoglu, E., Kumar, S., and Achan, K.  Inductive representation learning on temporal graphs, 2020.
  * Xu et al. (2018) Xu, K., Hu, W., Leskovec, J., and Jegelka, S.  How powerful are graph neural networks?  _arXiv preprint arXiv:1810.00826_ , 2018.
  * Ying et al. (2018) Ying, R., He, R., Chen, K., Eksombatchai, P., Hamilton, W. L., and Leskovec, J.  Graph convolutional neural networks for web-scale recommender systems.  In _Proceedings of the 24th ACM SIGKDD international conference on knowledge discovery & data mining_, pp. 974–983, 2018.
  * Zitnik et al. (2018) Zitnik, M., Agrawal, M., and Leskovec, J.  Modeling polypharmacy side effects with graph convolutional networks.  _Bioinformatics_ , 34(13):i457–i466, 2018.

##  Appendix A Appendix

###  A.1 TGR: Algorithm

This section provides full algorithm for TGR batch processing.

##### Get new and observed node IDs.

At the beginning of batch processing, we extract 200 node IDs and store them
in $n_{id}$. We then compare $n_{id}$ to node bank node IDs $n^{obs}_{id}$ and
extract $n^{seen}_{id}$ where we store previously observed nodes and
$n^{new}_{id}$ where we store unobserved nodes. After extracting
$n^{seen}_{id}$ and $n^{new}_{id}$, we update node bank $n^{obs}_{id}$ to use
in the next batch.

##### Compute input node features $\mathbf{X}(t)$.

We compute input node features $\mathbf{X}(t)$ as previously discussed in
section 3 through concatenating node $n^{seen}_{id}$ features from expander
memory and $n^{new}_{id}$ from TGNN temporal memory.

##### TGNN Forward Pass.

In TGN implementation $\mathbf{X}(t)$ includes node states for all node IDs in
$n_{id}$. The main difference between TGN implementation and ours is a
different set of input node features $\mathbf{X}(t)$. After forward pass, TGN
temporal memory is updated with new node states.

##### Extract node states from Temporal Memory.

We extract updated node states $\mathbf{S^{\prime}}(t)$ for node IDs
$n^{obs}_{id}$ from node bank. We pad $\mathbf{S^{\prime}}(t)$ to match
expander graph size. Note that $\mathbf{A}^{Cay}$ remains the same throughout
batch processing unless the node bank size grows larger than the size of
expander graph, as previously described in section 3.

##### Memory Mixing

Memory mixing module computes expander embeddings $\mathbf{H^{\prime}}(t)$ by
passing $\mathbf{S^{\prime}}(t)$ through an expander message passing layer.
After memory mixing expander memory is updated to include new expander
embeddings.

Algorithm 1 TGR Batch Processing

0: ExpanderMemory, TemporalMemory, Node Bank

1\. Get new and observed node IDs

$n_{id}^{new}\leftarrow n_{id}\setminus n_{id}^{obs}$ {Identify new nodes.}

$n_{id}^{seen}\leftarrow n_{id}\cap n_{id}^{obs}$ {Identify observed nodes.}

$n_{id}^{obs}\leftarrow n_{id}^{obs}\cup n_{id}^{new}$ {Update the Node Bank.}

2\. Compute input node features $\mathbf{X}(t)$

$\mathbf{X}(t)[n_{id}^{new}]$:
$\mathbf{S}(t^{-})\leftarrow\text{TemporalMemory}(n_{id}^{new})$

$\mathbf{X}(t)[n_{id}^{seen}]$:
$\mathbf{H}(t^{-})\leftarrow\text{ExpanderMemory}(n_{id}^{seen})$

3\. TGNN Forward Pass

$\mathbf{Z}(t)\leftarrow\text{TGNN}(\mathbf{X}(t),\mathbf{A}(t))$

4\. Extract node states from Temporal Memory

$\mathbf{S^{\prime}}(t)\leftarrow\text{TemporalMemory}(n_{id}^{obs})$

Pad $\mathbf{S^{\prime}}(t)$ to match $\mathbf{A}^{Cay}$ size:
$\mathbf{S^{\prime}}(t)\leftarrow\text{Pad}(\mathbf{S^{\prime}}(t))$

5\. Memory Mixing

$\mathbf{H^{\prime}}(t)\leftarrow\text{GNN}(\mathbf{S^{\prime}}(t),\mathbf{A}^{Cay})$

Update ExpanderMemory $\leftarrow\mathbf{H^{\prime}}(t)$

###  A.2 Experiments

####  A.2.1 Model Parameters

We do not perform extensive hyperparameter tuning in comparison to TGN
implementation (Huang et al., 2023) (see reported table LABEL:tab:parameters).
We match expander memory and embedding dimension to temporal memory.

Table 2: Model Hyperparameters. | Value
---|---
Temporal Memory Dimension | 100
Node Embedding Dimension | 100
Time Embedding Dimension | 100
Expander Memory Dimension | 100
Expander Embedding Dimension | 100
# Attention Heads | 2
Dropout | 0.1


# Mean-field approach to Rydberg facilitation in a gas of atoms at high and
low temperatures.

Daniel Brady  Department of Physics and Research Center OPTIMAS, University of
Kaiserslautern-Landau, D-67663 Kaiserslautern, Germany  Michael Fleischhauer
Department of Physics and Research Center OPTIMAS, University of
Kaiserslautern-Landau, D-67663 Kaiserslautern, Germany

###### Abstract

The excitation spread caused by Rydberg facilitation in a gas of laser driven
atoms is an interesting model system for studying epidemic dynamics. We derive
a mean-field approach to describe this facilitation process in the limits of
high and low temperatures, which takes into account Rydberg blockade and the
network character of excitation spreading in a low-temperature gas. As opposed
to previous mean-field models, our approach accurately predicts all stages of
the facilitation dynamics from the initial fast epidemic growth, an extended
saturation period, to the final relaxation phase.

##  I Introduction

Rydberg atoms have gained a lot of interest in the last few decades due to
their strongly exaggerated properties. In particular, they have very long life
times and strong interactions over distances covering several µm [1]. These
features allow Rydberg systems to be especially useful in a multitude of
applications such as quantum information processing [2, 3, 4, 5, 6] or the
study of many-body spin physics [7, 8, 9, 10, 11, 12, 13, 14].

One interesting process in many-body Rydberg systems is Rydberg facilitation,
which has been used to study dissipative spin systems [15], transport and
localization phenomena in disordered systems [16], or self-organized
criticality [17].

In this type of many-body system, atoms are coupled off-resonantly to a
Rydberg state. As a result of the Rydberg dipole interaction, however, atoms
near an already excited Rydberg atom can be excited resonantly. Thus an
initial seed excitation can lead to a cascade of excitations. It has been
shown experimentally that this type of system bears close similarities to
epidemic dynamics [18].

The most simple description of epidemic type systems is given by susceptible-
infected-susceptible (SIS) models. Here, each individual has two internal
states, susceptible (S) or infected (I). Susceptible individuals are infected
with rate $\kappa$ in the presence of other, already infected ones, and an
infected individual can return to the susceptible state with rate $\gamma$
[19, 20, 21].

Under the assumption of homogeneous mixing, where all individuals interact
with each other completely at random, all information about the epidemic
dynamics is contained in the total fractions $\rho^{\nu}$ in the susceptible
($\nu$=S) and infected state ($\nu$=I), which obey simple homogeneous mean-
field equations given by [22]

$\displaystyle\frac{\text{d}}{\text{d}t}\rho^{\textrm{I}}$
$\displaystyle=\kappa\rho^{\textrm{I}}\rho^{\textrm{S}}-\gamma\rho^{\textrm{I}},$
(1a) $\displaystyle\frac{\text{d}}{\text{d}t}\rho^{\textrm{S}}$
$\displaystyle=-\kappa\rho^{\textrm{I}}\rho^{\textrm{S}}+\gamma\rho^{\textrm{I}}.$
(1b)

These systems feature an absorbing-state phase transition between two
dynamical phases, namely an absorbing phase in which all infections die out,
and an active phase where - in the thermodynamic limit - infections last
forever. A suitable order parameter to distinguish these phases is the steady-
state active (infected) density $\rho_{\mathrm{ss}}^{\textrm{I}}$. From eq.
(1a) one recognizes that this phase transition occurs when

$\displaystyle\rho^{\textrm{S}}\kappa=\gamma.$ (2)

Figure 1: (a): Level scheme of internal atomic states with ground $\ket{g}$
(susceptible) state, Rydberg $\ket{r}$ (infected) state, and inert $\ket{0}$
(recovered) state. An external laser drives the transition between $\ket{g}$
and $\ket{r}$ and spontaneous decay occurs from $\ket{r}$ to either $\ket{g}$
or $\ket{0}$, modulated by branching parameter ${b\in[0,1]}$. (b): Monte-Carlo
data of the Rydberg density $\rho$ (red) and total active density $n$ (in
states $|g\rangle$ and $|r\rangle$) (blue) over time, with ${b=0.3}$ and
initial condition ${\rho(t=0)=0}$, showing the typical epidemic stages. The
blue dashed line corresponds to ${n/n_{\mathrm{crit}}=1}$. (c): Schematic of
Rydberg atoms (red dots) spanning facilitation shells (yellow region) and
blockade spheres (red region). Ground state atoms (grey dots) in the blue
region are subject to off-resonant laser coupling. (d): Schematic phase
diagram for the high (top) and low (bottom) temperature regimes depending on
the total gas density $n$.

This model can be expanded to include a recovered state, from which an
individual cannot be infected again, which corresponds to a susceptible-
infected-recovered (SIR) model [23, 24, 25]. For Rydberg systems this can be
approximately achieved by allowing an additional decay channel from the
Rydberg state to a dark state $\ket{0}$ (see Fig. 1a). A concrete mapping of
the Rydberg facilitation system to SIS and SIR models will be discussed later.

The dynamics of both SIS and SIR models are well understood in the homogeneous
mean-field regime. While the homogeneous mixing assumption is well justified
in systems where the infection spreading occurs on a regular lattice in high
spatial dimensions, it fails in many relevant cases, for example if the the
SIS/SIR dynamics take place on real-life networks. For such network systems,
including e.g. random Erdős–Rényi (ER) [26] and scale-free (SF) networks [27],
a large body of more sophisticated approximation methods have been developed
[28, 29, 30], but many questions remain unsolved. Here, Rydberg atoms provide
a platform to experimentally investigate the epidemic dynamics on a variety of
complex networks, which can e.g. be engineered by the use of tweezer arrays
[31]. Moreover, in a gas of atoms contained in some macroscopic trapping
potential one can investigate the transition between a random ER network at
very low temperatures, where the motion of atoms on the relevant time-scales
of the facilitation process can be ignored, to the homogeneous mixing limit
for a gas of high temperatures [32]. The latter corresponding to an annealed
random network.

With the addition of a slow decay channel into the inert (recovered) state,
the loss of individuals in the population (corresponding to the sum of Rydberg
and ground state densities) gradually moves the system into an absorbing-
state. This gives rise to three typical epidemic stages, which have been
experimentally observed in Rydberg facilitation systems [18]. Following an
initial infection, there is a rapid epidemic growth in infected individuals,
or Rydberg atoms, in the system. This is followed by a saturation and an
eventual relaxation at long times as a result of the system reaching an
absorbing-state on these time scales (see Fig. 1b).

In order to describe the macroscopic dynamics of the Rydberg facilitation
process in a gas, a simple mean-field model has been put forward in [17],
which however fails to provide a quantitative prediction of the microscopic
dynamics, accurately calculated by Monte-Carlo simulations [32] (see Fig. 2).
This discrepancy results from the mean-field model not regarding (i) Rydberg
blockade, which prevents the excitation of any atom closer than some radius
$r_{\mathrm{blockade}}$ to a Rydberg atom, and (ii) the emergent ER network at
low temperatures.

In the following, we will develop a mean-field description of the dynamics of
Rydberg excitations in a many-body facilitated gas that accounts for both of
these effects and provides accurate predictions of the full facilitation
dynamics, which we demonstrate by comparing our predictions with Monte-Carlo
simulations. The network structure of the cold gas leads to a higher total gas
density at very long times and subsequently a higher Rydberg density in the
saturation phase in comparison to mean-field predictions. Rydberg Blockade
causes a significant modification of the facilitation (infection) rate if the
density of Rydberg excited atoms (infected individuals) reaches some threshold
value. Similarly to the effect of regulatory measures on the dynamics of
epidemics (”lockdown”), it limits the maximum density of Rydberg-excited atoms
(infected individuals), but at the same time leads to a substantial
prolongation of the slow transition into the absorbing (recovery) phase.

##  II Microscopic Model of Rydberg facilitation

A microscopic description of Rydberg facilitation in a gas can be achieved
from a Lindblad master equation of the density matrix $\hat{\rho}$ which takes
the form

$\displaystyle\frac{\text{d}}{\text{d}t}\hat{\rho}=i[\hat{\rho},\hat{\mathcal{H}}]+\sum_{l}\hat{L}_{l}\hat{\rho}\hat{L}_{l}^{\dagger}-\frac{1}{2}\\{\hat{L}_{l}^{\dagger}\hat{L}_{l},\hat{\rho}\\},$
(3)

with the atom-light interaction Hamiltonian $\hat{\mathcal{H}}$ given by

$\displaystyle\hat{\mathcal{H}}$
$\displaystyle=\sum_{i}\Big{[}\Omega(\hat{\sigma}^{gr}_{i}+\hat{\sigma}^{rg}_{i})+\Big{(}\sum_{j\neq
i}\frac{c_{6}}{r_{ij}^{6}}\hat{\sigma}^{rr}_{j}-\Delta\Big{)}\hat{\sigma}^{rr}_{i}\Big{]}.$
(4)

Here $\hat{\sigma}_{i}^{\mu\nu}=|\mu\rangle_{ii}\langle\nu|$ is the projection
operator of the $i$-th atom from the internal state $\nu$ to $\mu$. The
external driving field is described by the Rabi-frequency $\Omega$ and the
detuning $\Delta$, and the van-der-Waals interaction energy between Rydberg
atoms $i$ and $j$ is given by ${c_{6}/r_{ij}^{6}}$, with $r_{ij}$ being the
distance between the atoms. Finally, in eq. (3) dissipation is described by
the jump operators $\hat{L}_{l}$. These take the form
${\hat{L}_{1}^{(i)}=\sqrt{(1-b)\gamma}\hat{\sigma}_{i}^{gr}}$,
${\hat{L}_{2}^{(i)}=\sqrt{b\gamma}\hat{\sigma}_{i}^{0r}}$ for spontaneous
decay from $\ket{r}$ to $\ket{g}$ and $\ket{0}$, as well as
${\hat{L}_{3}^{(i)}=\sqrt{\gamma_{\perp}}\hat{\sigma}_{i}^{rr}}$ for the
dephasing of the Rydberg state. The parameter ${b\in[0,1]}$ corresponds to the
percentage of Rydberg atoms which spontaneously decay to $\ket{0}$ and are
thereby removed from the system. Therefore, if ${b=0}$ the system corresponds
to a two-level system and resembles an SIS epidemic (see Fig. 1a).

Dephasing results from e.g. Doppler broadening or the spread of the atomic
wave packet over the van-der-Waals potential [17, 33]. In the large dephasing
limit, the dynamics of a many-body Rydberg gas are effectively governed by
classical rate equations [34]. As a result, this system can be simulated to
great accuracy using Monte-Carlo simulations. Starting from eq. (3), after
adiabatic elimination of coherences, one can formulate a set of rate equations
for the probabilities of atom $i$ being in the Rydberg state with
$P_{\mathrm{r}}^{(i)}$ or ground state with $P_{\mathrm{g}}^{(i)}$ as

$\displaystyle\frac{\text{d}}{\text{d}t}P_{\mathrm{r}}^{(i)}$
$\displaystyle=\Gamma_{\text{f}}(\Sigma)P_{\mathrm{g}}^{(i)}-\bigl{(}\Gamma_{\text{f}}(\Sigma)+\gamma\bigr{)}P_{\mathrm{r}}^{(i)}$
(5a) $\displaystyle\frac{\text{d}}{\text{d}t}P_{\mathrm{g}}^{(i)}$
$\displaystyle=\bigl{(}\Gamma_{\text{f}}(\Sigma)+(1-b)\gamma\bigr{)}P_{\mathrm{r}}^{(i)}-\Gamma_{\text{f}}(\Sigma)P_{\mathrm{g}}^{(i)},$
(5b)

with the stimulated excitation rate given by

$\displaystyle\Gamma_{\text{f}}(\Sigma)=\frac{2\Omega^{2}\gamma_{\perp}}{\gamma_{\perp}^{2}+\Delta^{2}\big{(}\sum_{\begin{subarray}{c}j\neq
i\\\
j\in\Sigma\end{subarray}}\frac{r_{\text{f}}^{6}}{r_{ij}^{6}}-1\big{)}^{2}}.$
(6)

Here $\Sigma$ is the set of indices of Rydberg-excited atoms. If no other
Rydberg atom exists in the gas or their distance is much larger than
$r_{\text{f}}$, $\Gamma_{\text{f}}(\Sigma)$ reduces to the off-resonant
excitation rate of an isolated atom

$\displaystyle\tau=\frac{2\Omega^{2}\gamma_{\perp}}{\gamma_{\perp}^{2}+\Delta^{2}}.$
(7)

If a Rydberg atom is present in the system, atoms located around a certain
distance to it, called facilitation distance $r_{\text{f}}$, are shifted into
resonance and can be excited on a much faster time-scale, given by the
facilitation rate ${\Gamma_{\text{f}}=2\Omega^{2}/\gamma_{\perp}}$. The
facilitation distance is given by

$\displaystyle r_{\text{f}}=\sqrt[6]{\frac{c_{6}}{\Delta}}.$ (8)

Rydberg facilitation can be observed when ${\Delta\gg\Omega}$, as this
naturally gives rise to a hierarchy in time-scales such that

$\displaystyle\Gamma_{\text{f}}\gg\gamma\gg\tau.$ (9)

In this case, off-resonant excitations and the decay of Rydberg atoms are
effectively static on the time-scale of facilitated excitations. Each Rydberg
atom spans a spherical shell, with distance $r_{\text{f}}$ and approximate
width $\delta r_{\text{f}}$, in which atoms are resonantly coupled to the
driving laser field. The width of the facilitation shell is thereby determined
by the effective linewidth of the excitation transition and reads

$\delta r_{\text{f}}=\frac{\gamma_{\perp}}{2\Delta}r_{\text{f}}.$ (10)

Atoms closer than $r_{\text{f}}-\delta r_{\text{f}}/2$ to a Rydberg atom are
subject to Rydberg blockade [35] and cannot be excited, since they are shifted
out of resonance again. These three regions around a Rydberg atom (off-
resonant coupling, facilitation, and blockade) can be schematically seen in
Fig. 1c.

In [17], a mean-field equation for a macroscopic description of the many-body
Rydberg facilitation dynamics has been derived. For a homogeneous gas this
reads

$\displaystyle\frac{\text{d}}{\text{d}t}\rho$
$\displaystyle=-\kappa(2\rho^{2}-\rho n)-\gamma\rho-\tau(2\rho-n),$ (11a)
$\displaystyle\frac{\text{d}}{\text{d}t}n$ $\displaystyle=-b\gamma\rho.$ (11b)

Here $\rho$ corresponds to the coarse grained Rydberg density (in a small
volume $\Delta V$)

$\displaystyle\rho(\vec{r},t)=\frac{1}{\Delta V}\sum_{i:\vec{r}_{i}\in\Delta
V}\langle\hat{\sigma}^{rr}_{i}\rangle,$ (12)

and $n$ is the density of ground and Rydberg state atoms

$\displaystyle n(\vec{r},t)=\frac{1}{\Delta V}\sum_{i:\vec{r}_{i}\in\Delta
V}\Bigl{(}\langle\hat{\sigma}^{rr}_{i}\rangle+\langle\hat{\sigma}^{gg}_{i}\rangle\Bigr{)}.$
(13)

Note that $n$ does not count $\ket{0}$ state atoms and therefore decreases
over time if ${b>0}$. The spreading rate of Rydberg excitations in the many-
body gas is given by the two-body facilitation rate integrated over the
facilitation shell

$\displaystyle\kappa=\Gamma_{\text{f}}V_{s},$ (14)

with the volume of the facilitation shell ${V_{s}\approx 4\pi\delta
r_{\text{f}}r_{\text{f}}^{2}}$. The above equations predict an absorbing-state
phase transition between an active and absorbing phase for the critical gas
density

$\displaystyle n_{\mathrm{crit}}=\frac{\gamma}{\Gamma_{\text{f}}V_{s}}.$ (15)

Equation (11a) for the (mean-field) Rydberg density in the many-body gas
strongly resembles the SIS equation of motion of infected individuals given by
(1a). However, in contrast to SIS/SIR epidemics, Rydberg systems additionally
feature (i) off-resonant excitations with rate $\tau$, (ii) resonant
(facilitated) de-excitations of Rydberg atoms (described by the term
${-2\kappa\rho^{2}}$ in eq. (11a)), and (iii) Rydberg blockade which is not
regarded in eq. (11a). Thus eqs. (11) do not directly map to those of the SIS
dynamics eqs.(1) (with ${b=0}$). However, identifying
$\rho^{\textrm{I}}=2\rho$ and $\rho^{\textrm{S}}=n-2\rho\geq 0$ and setting
$\tau=b=0$ reproduces the SIS equations with conserved total density
$\rho^{\text{I}}+\rho^{\text{S}}=n-2\rho+2\rho=n$. Moreover in the high
temperature limit, the de-excitation can effectively be neglected as it is a
second order process in terms of Rydberg density.

For the low temperature gas, the excitations dynamics are constrained to an ER
network in which the individual nodes are comprised of atoms (either in
$\ket{g}$ or $\ket{r}$), and connections between nodes, say $i$ and $j$ exist
if ${r_{ij}\approx r_{\text{f}}}$. The number of connections a node has (i.e.
the number of atoms with distance $r_{\text{f}}$ to an atom) is called the
degree $k$ of the atom. In an ER network the node degrees follow a Poissonian
distribution with average degree

$\displaystyle\langle k\rangle=nV_{s}.$ (16)

This system features a percolation transition at ${\langle k\rangle=1}$
between an (almost fully) connected network and a fragmented network,
comprised of many, small disconnected clusters. Here, clusters refers to a
group of connected nodes. From eq. (16), we can identify the gas density at
which the percolation transition occurs as

$\displaystyle n_{\mathrm{perc}}=\frac{1}{V_{s}}.$ (17)

This density is a factor ${\Gamma_{\text{f}}/\gamma}$ larger than the critical
density of the phase transition to the absorbing phase for a homogeneous gas
[32]. A schematic phase diagram for the high and low temperature gas can be
seen in Fig. 1d.

In this paper, we model the actual many-body dynamics using Monte-Carlo
simulations of the rate equations (5). We assume a cubic box with length
${L=7r_{\text{f}}}$ and periodic boundary conditions. Atom positions are
chosen randomly and velocities are sampled from a Maxwell-Boltzmann
distribution with temperature parameter $\hat{v}$, corresponding to the most
probable atom velocity in the gas. For the time evolution we utilize a fixed
time step Monte-Carlo algorithm [36], with the time step
${dt=1/400\,\gamma^{-1}}$, as long-range interactions paired with the fast
movement of atoms in the high temperature case results in quickly changing
transitional rates in the system. To ensure numeric stability the $c_{6}$
potential in eq. (6) is truncated at a cutoff value around the singularity
${r_{ij}\to 0}$.

##  III Modified Langevin Description of Epidemic Evolution

In the following section, we develop an effective macroscopic theory of the
Rydberg facilitation process, expanding the Langevin equation (11), starting
from the microscopic model. This new equation takes into account Rydberg
blockade, as well as the network structure in the case of the low temperature
gas.

In Fig. 2, the dynamics of Rydberg excitations predicted by the improved
Langevin equation and by Monte-Carlo simulations are compared for the cases:
(a) high temperature gas with ${n_{0}>n_{\mathrm{crit}}}$, (b) low temperature
gas initially above the percolation threshold with
${n_{0}>n_{\mathrm{perc}}}$, and (c) low temperature gas initially below the
percolation threshold ${n_{0}<n_{\mathrm{perc}}}$. Here $n_{0}$ refers to the
gas density at ${t=0}$. Additionally, we use a branching parameter ${b=0.3}$,
allowing some loss into the recovered state $\ket{0}$. Therefore, for all
cases (a-c), the system drives itself to the absorbing-state and follows the
typical epidemic stages as seen in Fig. 1b.

We start from the microscopic Heisenberg-Langevin equations describing the
quantum many-body dynamics of Rydberg excitations for atoms at given spatial
positions given by

$\displaystyle\frac{\text{d}}{\text{d}t}\hat{\sigma}^{rr}_{i}$
$\displaystyle=-i\Omega(\hat{\sigma}^{rg}_{i}-\hat{\sigma}^{gr}_{i})-\gamma\hat{\sigma}^{rr}_{i}+\hat{\xi}_{1},$
(18) $\displaystyle\frac{\text{d}}{\text{d}t}\hat{\sigma}^{rg}_{i}$
$\displaystyle=-i\Big{(}\Omega(\hat{\sigma}^{rr}_{i}-\hat{\sigma}^{gg}_{i})-\hat{V}_{i}\hat{\sigma}^{rg}_{i}\Big{)}-\gamma_{\perp}\hat{\sigma}^{rg}_{i}+\hat{\xi}_{2}.$
(19)

We have accounted for losses and dephasing and added corresponding fluctuation
operators $\hat{\xi}_{1}$ and $\hat{\xi}_{2}$, which disappear in the quantum
mechanical average and whose properties can be obtained from the fluctuation-
dissipation relation [37]. The operator $\hat{V}_{i}$ describes the detuning
of the $i$-th atom and depends on the states of all other atoms. It is given
by

$\displaystyle\hat{V}_{i}=\Delta\Big{(}-1+\sum_{j\neq
i}\frac{r_{\text{f}}^{6}}{r_{ij}^{6}}\hat{\sigma}^{rr}_{j}\Big{)}.$ (20)

We note that the operator valued quantities are objects in Hilbert space
describing the quantum mechanical evolution and are subject to the classical
statistics of the (time dependent) random positions of the atoms. The dynamics
of the atom positions are treated classically, which is well justified in the
high-dephasing limit, assumed throughout the present paper.

Furthermore, the effect of $c_{6}$ forces acting on the center-of-mass motion
of the atoms due to the distance dependence of $\hat{V}_{i}$ are disregarded
in the present paper. They will be discussed elsewhere in more detail [38],
where we will show that under typical experimental conditions they can be
accounted for by a change of the atoms velocity distribution and, in the case
of a trapped gas, by an additional loss channel.

Assuming high dephasing ${\gamma_{\perp}\geq\Omega}$, the coherences
$\hat{\sigma}^{rg}_{i}$ quickly decay to quasi-stationary values relative to
the relevant many-body time scales. Therefore, we adiabatically eliminate
coherences (${\frac{\text{d}}{\text{d}t}\hat{\sigma}^{rg}_{i}=0}$) and arrive
at

$\displaystyle\frac{\text{d}}{\text{d}t}\hat{\sigma}^{rr}_{i}$
$\displaystyle=-\frac{2\Omega^{2}{\gamma_{\perp}}}{{\gamma_{\perp}}^{2}+\hat{V}_{i}^{2}}(\hat{\sigma}^{rr}_{i}-\hat{\sigma}^{gg}_{i})-\gamma\hat{\sigma}^{rr}_{i}+\hat{\xi}.$
(21)

In the following, we expand the leading fraction in eq. (21) with a full basis
of projection operators ${\Pi}_{i}(m)$ projecting onto $m$ Rydberg atoms in
the facilitation sphere of atom $i$. Each of these atoms have a relative
distance to atom $i$ in the interval $0\leq{r_{ij}\leq r_{\text{f}},\;\forall
j\in[1,m]}$. Using this, the equation of motion takes the form

$\displaystyle\frac{\text{d}}{\text{d}t}\hat{\sigma}^{rr}_{i}$
$\displaystyle=-{\Pi}_{i}(0)\frac{2\Omega^{2}{\gamma_{\perp}}}{{\gamma_{\perp}}^{2}+\Delta^{2}}(\hat{\sigma}^{rr}_{i}-\hat{\sigma}^{gg}_{i})$
$\displaystyle-{\Pi}_{i}(1)\underbrace{\frac{2\Omega^{2}{\gamma_{\perp}}}{{\gamma_{\perp}}^{2}+\Delta^{2}\big{(}(\frac{r_{\text{f}}}{r_{1i}})^{6}-1\big{)}^{2}}}_{(*)}(\hat{\sigma}^{rr}_{1}\hat{\sigma}^{rr}_{i}-\hat{\sigma}^{rr}_{1}\hat{\sigma}^{gg}_{i})$
$\displaystyle+...$ $\displaystyle-\gamma\hat{\sigma}^{rr}_{i}+\hat{\xi}.$
(22)

All rates for more than one Rydberg atom in the facilitation sphere ${m>1}$
are strongly suppressed due to blockade. As a result, we truncate the
expansion at ${m=1}$.

The fraction $(*)$ in eq. (22) is the distance dependent facilitation rate for
an atom in the presence of one Rydberg atom. Here, we approximate this rate as
being non-zero only if atom $i$ is in the facilitation _shell_ of the Rydberg
atom ${r_{1i}\in S_{\text{f}}\equiv[r_{\text{f}}-\frac{\delta
r_{\text{f}}}{2},r_{\text{f}}+\frac{\delta r_{\text{f}}}{2}]}$. We express
this as

$\displaystyle\hat{\sigma}^{rr}_{1}\Pi_{i}(1)\to\hat{\sigma}^{rr}_{1}\,\Xi_{i}(1)=\hat{\sigma}^{rr}_{1}\times\begin{cases}1,\quad
r_{1i}\in S_{\text{f}}\\\ 0,\quad\text{else}\end{cases},$ (23)

and obtain

$\displaystyle\frac{\text{d}}{\text{d}t}\hat{\sigma}^{rr}_{i}$
$\displaystyle=-\tau\,{\Pi}_{i}(0)(\hat{\sigma}^{rr}_{i}-\hat{\sigma}^{gg}_{i})$
$\displaystyle-\Gamma_{\text{f}}\,{\Xi}_{i}(1)(\hat{\sigma}^{rr}_{1}\hat{\sigma}^{rr}_{i}-\hat{\sigma}^{rr}_{1}\hat{\sigma}^{gg}_{i})$
$\displaystyle-\gamma\hat{\sigma}^{rr}_{i}+\hat{\xi},$ (24)

with the maximal facilitation rate
$\Gamma_{\text{f}}=\frac{2\Omega^{2}}{{\gamma_{\perp}}}$ and the off-resonant
excitation rate
$\tau=\frac{2\Omega^{2}\gamma_{\perp}}{{\gamma_{\perp}}^{2}+\Delta^{2}}$.

Finally, we calculate the expectation value of the operator
$\hat{\sigma}^{rr}_{i}$ with a double averaging over the quantum mechanical
state and the ensemble of the many different atom positions in the gas. We
will denote these double averages as
$\langle\langle\hat{\sigma}^{rr}_{i}\rangle\rangle$ and write

$\displaystyle\frac{\text{d}}{\text{d}t}\langle\langle\hat{\sigma}^{rr}_{i}\rangle\rangle$
$\displaystyle=-\tau\langle\langle{\Pi}_{i}(0)(\hat{\sigma}^{rr}_{i}-\hat{\sigma}^{gg}_{i})\rangle\rangle$
$\displaystyle-\Gamma_{\text{f}}\,\langle\langle{\Xi}_{i}(1)\hat{\sigma}^{rr}_{1}(\hat{\sigma}^{rr}_{i}-\hat{\sigma}^{gg}_{i})\rangle\rangle$
$\displaystyle-\gamma\langle\langle\hat{\sigma}^{rr}_{i}\rangle\rangle$ (25)
$\displaystyle\approx-\tau\langle\langle{\Pi}_{i}(0)\rangle\rangle\Bigl{(}\langle\langle\hat{\sigma}^{rr}_{i}\rangle\rangle-\langle\langle\hat{\sigma}^{gg}_{i}\rangle\rangle\Bigr{)}$
$\displaystyle-\Gamma_{\text{f}}\,\langle\langle{\Pi}_{i}(1)\hat{\sigma}^{rr}_{1}\rangle\rangle\Bigl{(}\langle\langle\hat{\sigma}^{rr}_{i}\rangle\rangle-\langle\langle\hat{\sigma}^{gg}_{i}\rangle\rangle\Bigr{)}$
$\displaystyle-\gamma\langle\langle\hat{\sigma}^{rr}_{i}\rangle\rangle.$ (26)

Assuming a randomly distributed gas, we can approximate the probabilities
$\langle\langle{\Pi}_{i}(m)\rangle\rangle$ as Poissonian with the rate $\rho
V_{\textrm{f}}$, i.e. $\langle\langle{\Pi}_{i}(m)\rangle\rangle=(\rho
V_{f})^{m}e^{-\rho V_{f}}/m!$, resulting in

$\displaystyle\langle\langle{\Pi}_{i}(0)\rangle\rangle$
$\displaystyle=\text{e}^{-\rho V_{\textrm{f}}}$ (27)
$\displaystyle\langle\langle{\Xi}_{i}(1)\hat{\sigma}^{rr}_{1}\rangle\rangle$
$\displaystyle\equiv\frac{V_{s}}{V_{f}}\langle\langle{\Pi}_{i}(1)\rangle\rangle=\rho
V_{\textrm{s}}\,\text{e}^{-\rho V_{\textrm{f}}}.$ (28)

The factor $V_{s}/V_{f}$ in the second line takes into account that only atoms
in the facilitation shell contribute. We then perform the coarse-graining
given by eqs. (12) and (13) and arrive at

$\displaystyle\frac{\text{d}}{\text{d}t}\rho=$
$\displaystyle-\kappa\text{e}^{-\rho V_{\text{f}}}\rho(2\rho-n)$
$\displaystyle-\gamma\rho-\tau(2\rho-n).$ (29)

Furthermore, we have assumed here ${\text{e}^{-\rho
V_{\text{f}}}\tau\approx\tau}$ as the off-resonant rate is only relevant when
${\rho V_{\text{f}}\ll 1}$.

The spreading rate of excitations $\kappa$ is now exponentially damped by the
density of Rydberg atoms. This gives a better description of the spreading of
Rydberg excitations in the epidemic growth stage. However, as all atoms with
distances closer than $r_{\mathrm{blockade}}$ to a Rydberg atom cannot be
excited due to Rydberg blockade (red region in Fig. 1c), there exists a
maximum density of Rydberg atoms $\rho_{\mathrm{max}}$, given by the packing
density of non-overlapping spheres, above which no more excitations are
possible.

To quantify the blockade induced saturation density in the gas, we introduce
the parameter $\eta$ corresponding to the packing density of spheres in a
given volume.

For the high temperature gas, this corresponds to the densest packing of
spheres, given by ${\eta=\frac{\pi}{3\sqrt{2}}\approx 74.0\%}$. In this regime
we can assume this packing density to be achieved, as the high motion of the
atoms allow the system to organize itself to this state.

Figure 2: Rydberg density $\rho$ over time for ${\rho(t=0)=0}$ and ${b=0.3}$
modeled with Monte-Carlo simulations (blue dots), eq. (11) (red, faint), and
eq. (III) (red, solid), for a high temperature gas with starting density
${n_{0}=4.0\,r_{\text{f}}^{-3}>n_{\mathrm{crit}}}$, b low temperature
percolating gas with ${n_{0}=20\,r_{\text{f}}^{-3}>n_{\mathrm{perc}}}$, and c
low temperature non-percolating gas with
${n_{0}=4\,r_{\text{f}}^{-3}<n_{\mathrm{perc}}}$.

For the frozen gas, the packing density is slightly lower, and is given by the
closest density of randomly packed spheres, which is given by ${\eta\approx
63.5\%}$ [39].

As $\delta r_{\text{f}}\ll r_{\text{f}}$ we can approximate the blockade
radius as ${r_{\mathrm{blockade}}\approx r_{\text{f}}}$ and write

$\displaystyle\rho_{\mathrm{max}}=2\frac{\eta}{V_{\mathrm{fac}}},$ (30)

with the approximate volume of the blockade sphere
${V_{\mathrm{fac}}=\frac{4}{3}\pi r_{\text{f}}^{3}}$. The factor of 2 emerges
as when a facilitation event occurs, the facilitated atom is centered on the
blockade sphere of the facilitating Rydberg atom. As a result, the blockade
spheres of these atoms overlap. If, however, a third Rydberg atom is
facilitated (by the second Rydberg atom), its blockade sphere borders the
blockade sphere of the first Rydberg atom with, on average, very little
overlap (see Fig. 1c).

As the laser coupling smoothly changes from resonant, for an atom with
distance ${r=r_{\text{f}}}$ to a Rydberg atom, to strongly suppressed for
${r<r_{\text{f}}}$, this can be regarded as a packing of soft spheres with an
uncertainty in volume of ${\delta V_{\mathrm{fac}}=4\pi\delta
r_{\text{f}}r_{\text{f}}^{2}}$. The result is a smearing out of
$\rho_{\mathrm{max}}$ given by ${\delta\rho_{\mathrm{max}}=\delta
V_{\mathrm{fac}}\frac{2\eta}{V_{\mathrm{fac}}^{2}}}$. We can now add a
heuristic function which sets the facilitation rate to 0 if
${\rho>\rho_{\mathrm{max}}}$ as

$\displaystyle
h(\rho)=\frac{1}{2}\Big{(}1+\tanh\Big{(}\frac{\rho_{\mathrm{max}}-\rho}{\delta\rho_{\mathrm{max}}}\Big{)}\Big{)}.$
(31)

The added factors $\text{e}^{-\rho V_{f}}$ and $h(\rho)$ to the facilitation
rate $\kappa$ suffice to accurately describe the dynamics of the Rydberg
density in the high temperature gas (see Fig. 2a).

For the low temperature gas the finite connectivity greatly reduces the
facilitation rate. Taking into account that facilitation can only occur if the
degree of the atom $k$ is not $0$, we alter the facilitation rate to

$\displaystyle\kappa\to\kappa\,P(k>0).$ (32)

For an ER network with average degree $\langle k\rangle\ll 1$, we can
approximate $P(k>0)\approx\langle k\rangle$. In this case the new infection
rate $\kappa$ corresponds to the Kephart-White model [40, 41].

The full Langevin equation for the Rydberg density reads

$\displaystyle\frac{\text{d}}{\text{d}t}\rho=$
$\displaystyle-\kappa\text{e}^{-\rho V_{\text{f}}}h(\rho)P(k>0)\rho(2\rho-n)$
$\displaystyle-\gamma\rho-\tau(2\rho-n),$ (33)

with ${P(k>0)=1-\text{e}^{-nV_{s}}}$ for the low temperature gas and
${P(k>0)=1}$ at high temperatures.

In Fig. 2 we have compared the predictions from the modified Langevin equation
(III) with Monte Carlo simulations in the high temperature gas, the frozen
percolating gas, and the frozen non-percolating gas. For the high temperature
and the frozen percolating case, eq. (III) has a very good agreement for all
epidemic stages with Monte-Carlo data. In particular, it predicts the correct
density in the saturation stage in the high temperature and the low
temperature percolating gas in contrast to eq. (11).

Furthermore, for the case of the low temperature gas, eq. (III) gives a much
better approximation of the relaxation epidemic stage (i.e. for times ${\gamma
t\gtrsim 10^{2}}$). In this stage, the Rydberg density is much higher than the
expected MF density (predicted by the faint red line), which holds for high
temperatures. In contrast to the high temperature regime, the system leaves
the active phase at much higher gas density due to the finite connectivity of
excitation paths in the gas. The factor ${P(k>0)}$ in the facilitation rate
gives a much better approximation of this increased Rydberg density.

##  IV Conclusion

In conclusion, we have developed a modified mean-field approach to model the
Rydberg density over time in a many-body gas under facilitation conditions for
the limits of high and low temperature. In the low temperature regime, we have
additionally differentiated between a system with initial density
${n_{0}>n_{\mathrm{prec}}}$ and ${n_{0}<n_{\mathrm{perc}}}$, where
$n_{\mathrm{perc}}$ is the percolation density below which heterogeneous
effects play a large role.

Our modelling is similar to that developed in [17], but with three key
improvements to the facilitation (or infection) rate $\kappa$. We consider (i)
random atom positions leading to a Poissonian distribution in the number of
Rydberg atoms closer than $r_{\mathrm{blockade}}$ to a given atom. In this
case, the atom cannot be excited or de-excited due to Rydberg blockade. As a
result, with increasing Rydberg density, the global facilitation rate $\kappa$
exponentially decreases.

Additionally, (ii) excited Rydberg atoms can be seen as soft spheres inside of
which no atoms can be excited due to blockade. Therefore, there exists a
tightest packing of excited atoms beyond which the facilitation rate $\kappa$
vanishes. In the high temperature regime, this packing density corresponds to
the tightest packing of spheres in a given volume, as the high thermal
velocities allow the system to continuously organize itself to this state. In
the low temperature regime the packing density is given by that of randomly
placed spheres in a given volume, which, in comparison, is slightly lower.

Finally, (iii) for the low temperature regime, one has to additionally
consider the finite connectivity of the underlying network along which
facilitated excitations can spread. On a mean-field level, we have described
this by reducing the facilitation rate in correspondence with the portion of
atoms with network degree (i.e. the number of atoms in their facilitation
shell) ${k=0}$. The percentage of these isolated atoms increases as the
network connectivity decreases, and is therefore dependent on the total
density of the gas.

For both the high temperature, as well as the low temperature, high density
case eq. (III) gives excellent correspondence to Monte-Carlo data for all
epidemic stages.

For the low temperature, low density gas the system is characterized by strong
heterogeneity making an accurate mean-field description challenging. However,
for this case we still see a large improvement in the Langevin description of
the dynamics.

### Acknowledgement

The authors thank Simon Ohler and Johannes Otterbach for fruitful discussions.
Financial support from the DFG through SFB TR 185, project number 277625399,
is gratefully acknowledged.

## References

  * Gallagher [2006] T. Gallagher,  _Rydberg atoms_ (2006).
  * Jaksch  _et al._ [2000] D. Jaksch, J. I. Cirac, P. Zoller, S. L. Rolston, R. Côté, and M. D. Lukin, Fast quantum gates for neutral atoms, Phys. Rev. Lett. 85, 2208 (2000).
  * Lukin  _et al._ [2001a] M. D. Lukin, M. Fleischhauer, R. Cote, L. M. Duan, D. Jaksch, J. I. Cirac, and P. Zoller, Dipole blockade and quantum information processing in mesoscopic atomic ensembles, Phys. Rev. Lett. 87, 037901 (2001a).
  * Gaetan  _et al._ [2009] A. Gaetan, Y. Miroshnychenko, T. Wilk, A. Chotia, M. Viteau, D. Comparat, P. Pillet, A. Browaeys, and P. Grangier, Observation of collective excitation of two individual atoms in the Rydberg blockade regime, Nature Physics 5, 115 (2009).
  * Urban  _et al._ [2009] E. Urban, T. A. Johnson, T. Henage, L. Isenhower, D. Yavuz, T. Walker, and M. Saffman, Observation of Rydberg blockade between two atoms, Nature Physics 5, 110 (2009).
  * Saffman  _et al._ [2010] M. Saffman, T. G. Walker, and K. Mølmer, Quantum information with rydberg atoms, Reviews of modern physics 82, 2313 (2010).
  * Weimer  _et al._ [2010] H. Weimer, M. Müller, I. Lesanovsky, P. Zoller, and H. P. Büchler, A Rydberg quantum simulator, Nature Physics 6, 382 (2010).
  * Schauß  _et al._ [2012] P. Schauß, M. Cheneau, M. Endres, T. Fukuhara, S. Hild, A. Omran, T. Pohl, C. Gross, S. Kuhr, and I. Bloch, Observation of spatially ordered structures in a two-dimensional Rydberg gas, Nature 491, 87 (2012).
  * 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, Probing many-body dynamics on a 51-atom quantum simulator, Nature 551, 579 (2017).
  * Browaeys and Lahaye [2020a] A. Browaeys and T. Lahaye, Many-body physics with individually controlled Rydberg atoms, Nature Physics 16, 132 (2020a).
  * Surace  _et al._ [2020] F. M. Surace, P. P. Mazza, G. Giudici, A. Lerose, A. Gambassi, and M. Dalmonte, Lattice gauge theories and string dynamics in Rydberg atom quantum simulators, Phys. Rev. X 10, 021041 (2020).
  * Scholl  _et al._ [2021] P. Scholl, M. Schuler, H. J. Williams, A. A. Eberharter, D. Barredo, K.-N. Schymik, V. Lienhard, L.-P. Henry, T. C. Lang, T. Lahaye, A. M. Läuchli, and A. Browaeys, Quantum simulation of 2D antiferromagnets with hundreds of Rydberg atoms, Nature 595, 233 (2021).
  * de Léséleuc  _et al._ [2019] S. de Léséleuc, V. Lienhard, P. Scholl, D. Barredo, S. Weber, N. Lang, H. P. Büchler, T. Lahaye, and A. Browaeys, Observation of a symmetry-protected topological phase of interacting bosons with Rydberg atoms, Science 365, 775 (2019), https://science.sciencemag.org/content/365/6455/775.full.pdf .
  * Semeghini  _et al._ [2021] G. Semeghini, H. Levine, A. Keesling, S. Ebadi, T. T. Wang, D. Bluvstein, R. Verresen, H. Pichler, M. Kalinowski, R. Samajdar, A. Omran, S. Sachdev, A. Vishwanath, M. Greiner, V. Vuletić, and M. D. Lukin, Probing topological spin liquids on a programmable quantum simulator, Science 374, 1242 (2021).
  * Helmrich  _et al._ [2018] S. Helmrich, A. Arias, and S. Whitlock, Uncovering the nonequilibrium phase structure of an open quantum spin system, Phys. Rev. A 98, 022109 (2018).
  * Marcuzzi  _et al._ [2017] M. Marcuzzi, J. Minář, D. Barredo, S. De Léséleuc, H. Labuhn, T. Lahaye, A. Browaeys, E. Levi, and I. Lesanovsky, Facilitation dynamics and localization phenomena in rydberg lattice gases with position disorder, Physical Review Letters 118, 063606 (2017).
  * Helmrich  _et al._ [2020] S. Helmrich, A. Arias, G. Lochead, T. M. Wintermantel, M. Buchhold, S. Diehl, and S. Whitlock, Signatures of self-organized criticality in an ultracold atomic gas, Nature 577, 481 (2020).
  * Wintermantel  _et al._ [2021] T. M. Wintermantel, M. Buchhold, S. Shevate, M. Morgado, Y. Wang, G. Lochead, S. Diehl, and S. Whitlock, Epidemic growth and griffiths effects on an emergent network of excited atoms, Nature Communications 12, 103 (2021).
  * Anderson and May [1991] R. M. Anderson and R. M. May,  _Infectious diseases of humans: dynamics and control_ (Oxford university press, 1991).
  * Bailey  _et al._ [1975] N. T. Bailey _et al._ ,  _The mathematical theory of infectious diseases and its applications_ (Charles Griffin & Company Ltd, 5a Crendon Street, High Wycombe, Bucks HP13 6LE., 1975).
  * Murray and Murray [1993] J. D. Murray and J. D. Murray, Epidemic models and the dynamics of infectious diseases, Mathematical biology , 610 (1993).
  * Pastor-Satorras  _et al._ [2015] R. Pastor-Satorras, C. Castellano, P. Van Mieghem, and A. Vespignani, Epidemic processes in complex networks, Reviews of modern physics 87, 925 (2015).
  * Harko  _et al._ [2014] T. Harko, F. S. Lobo, and M. Mak, Exact analytical solutions of the susceptible-infected-recovered (sir) epidemic model and of the sir model with equal death and birth rates, Applied Mathematics and Computation 236, 184 (2014).
  * Bjørnstad  _et al._ [2020] O. N. Bjørnstad, K. Shea, M. Krzywinski, and N. Altman, Modeling infectious epidemics, Nature Methods 17, 455 (2020).
  * Beckley  _et al._ [2013] R. Beckley, C. Weatherspoon, M. Alexander, M. Chandler, A. Johnson, and G. S. Bhatt, Modeling epidemics with differential equations, Tennessee State University Internal Report (2013).
  * Erdős  _et al._ [1960] P. Erdős, A. Rényi, _et al._ , On the evolution of random graphs, Publ. Math. Inst. Hung. Acad. Sci 5, 17 (1960).
  * Barabási [2009] A.-L. Barabási, Scale-free networks: a decade and beyond, science 325, 412 (2009).
  * Pastor-Satorras and Vespignani [2001] R. Pastor-Satorras and A. Vespignani, Epidemic spreading in scale-free networks, Physical review letters 86, 3200 (2001).
  * Pastor-Satorras and Vespignani [2002] R. Pastor-Satorras and A. Vespignani, Epidemic dynamics in finite size scale-free networks, Physical Review E 65, 035108 (2002).
  * Durón and Farrell [2022] C. Durón and A. Farrell, A mean-field approximation of sir epidemics on an erdös–rényi network model, Bulletin of Mathematical Biology 84, 70 (2022).
  * Browaeys and Lahaye [2020b] A. Browaeys and T. Lahaye, Many-body physics with individually controlled rydberg atoms, Nature Physics 16, 132 (2020b).
  * Brady  _et al._ [2023] D. Brady, J. Bender, P. Mischke, T. Niederprüm, H. Ott, and M. Fleischhauer, Griffiths phase in a facilitated rydberg gas at low temperature (2023), arXiv:2302.14145 [cond-mat.quant-gas] .
  * Li  _et al._ [2013] W. Li, C. Ates, and I. Lesanovsky, Nonadiabatic motional effects and dissipative blockade for rydberg atoms excited from optical lattices or microtraps, Phys. Rev. Lett. 110, 213005 (2013).
  * Levi  _et al._ [2016] E. Levi, R. Gutiérrez, and I. Lesanovsky, Quantum non-equilibrium dynamics of rydberg gases in the presence of dephasing noise of different strengths, Journal of Physics B: Atomic, Molecular and Optical Physics 49, 184003 (2016).
  * Lukin  _et al._ [2001b] M. D. Lukin, M. Fleischhauer, R. Cote, L. Duan, D. Jaksch, J. I. Cirac, and P. Zoller, Dipole blockade and quantum information processing in mesoscopic atomic ensembles, Physical review letters 87, 037901 (2001b).
  * Ruiz Barlett  _et al._ [2009] V. Ruiz Barlett, J. Bigeón, M. Hoyuelos, and H. Mártin, Differences between fixed time step and kinetic monte carlo methods for biased diffusion, Journal of Computational Physics 228, 5740 (2009).
  * Kubo [1966] R. Kubo, The fluctuation-dissipation theorem, Reports on progress in physics 29, 255 (1966).
  * [38] D. Brady and M. Fleischhauer, in preparation .
  * Wu  _et al._ [2003] Y. Wu, Z. Fan, and Y. Lu, Bulk and interior packing densities of random close packing of hard spheres, Journal of Materials Science 38, 2019 (2003).
  * Kephart and White [1992] J. O. Kephart and S. R. White, Directed-graph epidemiological models of computer viruses, in  _Computation: the micro and the macro view_ (World Scientific, 1992) pp. 71–102.
  * Chakrabarti  _et al._ [2008] D. Chakrabarti, Y. Wang, C. Wang, J. Leskovec, and C. Faloutsos, Epidemic thresholds in real networks, ACM Transactions on Information and System Security (TISSEC) 10, 1 (2008).