# Locality of critical percolation on expanding graph sequences

Michael Ren⋆∗ and  Nike Sun⋆

###### Abstract.

We study the locality of critical percolation on finite graphs: let $G_{n}$ be
a sequence of finite graphs, converging locally weakly to a (random, rooted)
infinite graph $G$. Consider Bernoulli edge percolation: does the critical
probability for the emergence of an infinite component on $G$ coincide with
the critical probability for the emergence of a linear-sized component on
$G_{n}$? In this short article we give a positive answer provided the graphs
$G_{n}$ satisfy an expansion condition, and the limiting graph $G$ has finite
expected root degree. The main result of Benjamini, Nachmias, and Peres
(2011), where this question was first formulated, showed the result assuming
the $G_{n}$ satisfy a uniform degree bound and uniform expansion condition,
and converge to a deterministic limit $G$. Later work of Sarkar (2021)
extended the result to allow for a random limit $G$, but still required a
uniform degree bound and uniform expansion for $G_{n}$. Our result replaces
the degree bound on $G_{n}$ with the (milder) requirement that $G$ must have
finite expected root degree. Our proof is a modification of the previous
results, using a pruning procedure and the second moment method to control
unbounded degrees.

⋆Massachusetts Institute of Technology. ∗University of Cambridge.

##  1\. Introduction

In this article we study the question of locality of critical percolation, as
formulated by [BNP11]. Informally, if $G_{n}$ is a graph sequence converging
locally weakly to the random graph $(G,\rho)$, then does Bernoulli bond
percolation have the same critical probability on $G$ as on $G_{n}$? The main
result of [BNP11] gives a positive answer in the case that the $G_{n}$ have a
uniform Cheeger constant $h>0$ and a uniform degree bound $d<\infty$, and
converges locally weakly to a deterministic limiting graph $(G,\rho)$. A later
work [Sar21] extends this result to the more general case of a random limiting
graph $(G,\rho)$, but still requires that the graphs $G_{n}$ satisfy uniform
expansion and a uniform degree bound. In this article we extend the result of
[Sar21] by further relaxing the degree condition: more precisely, we show that
the degree restriction on $G_{n}$ can be replaced with the requirement that
the limit $(G,\rho)$ must have finite expected root degree.

###  1.1. Statement of main result

To formally state our main result, we set some notations and review some
standard definitions (see [AL07]). If $v$ is a vertex in graph $H$, we let
$B_{R}(v;H)$ be the $R$-neighborhood of $v$ in $H$; we can regard $B_{R}(v;H)$
as a graph rooted at $v$.

######  Definition 1 (space of rooted graphs).

Let $\mathscr{G}_{\star}$ denote the space of connected locally finite rooted
graphs $(G,\rho)$, modulo rooted isomorphism. A metric on
$\mathscr{G}_{\star}$ is given by

$\bm{d}\Big{(}(G_{1},\rho_{1}),(G_{2},\rho_{2})\Big{)}=\inf\bigg{\\{}\frac{1}{1+R}:R\geq
0,B_{R}(\rho_{1};G_{1})\cong B_{R}(\rho_{2};G_{2})\bigg{\\}}\in[0,1]\,.$ (1.1)

The space $\mathscr{G}_{\star}$ is separable and complete in this metric.

######  Definition 2 (local weak convergence).

Suppose $G_{n}\equiv(V_{n},E_{n})$ is a sequence of finite graphs, and let
$U_{n}$ be a uniformly random vertex in $G_{n}$. We say $G_{n}$ converges
locally weakly to the (random) element $(G,\rho)\in\mathscr{G}_{\star}$ if
$B_{R}(U_{n};G_{n})$ converges in law to $(G,\rho)$ in the topology of
$\mathscr{G}_{\star}$. This will be denoted $G_{n}\to_{\textit{lwc}}(G,\rho)$;
we will often write simply $G_{n}\to_{\textit{lwc}}G$ with the understanding
that $G$ is a (random) rooted graph.

######  Definition 3 (expanding graph sequences).

Given a finite graph $G=(V,E)$ and $0<\delta\leq 1/2$, define

$h_{\delta}(G)\equiv\min\bigg{\\{}\frac{|E_{G}(A,V\setminus
A)|}{|A|}:A\subseteq V,\delta\leq\frac{|A|}{|V|}\leq\frac{1}{2}\bigg{\\}}\,,$

where $E_{G}(A,B)$ denotes the set of all edges in $G$ between the vertex sets
$A$ and $B$. We say that a sequence $G_{n}$ of finite graphs is expanding if
for all $\delta>0$ we have $\liminf_{n}h_{\delta}(G_{n})>0$. That is to say,
for all $\delta>0$ there exist $n_{\delta}<\infty$ and $c_{\delta}>0$ such
that $h_{\delta}(G_{n})\geq c_{\delta}$ for all $n\geq n_{\delta}$.

######  Definition 4 (critical percolation).

Let $G$ be an infinite graph, and let $\mathbb{P}_{p}$ denote the law of a
Bernoulli bond percolation (or edge percolation) configuration
$\omega=\omega(G)$ on $G$. Then

$p_{c}(G)=\sup_{p}\bigg{\\{}p\in[0,1]:\mathbb{P}_{p}\Big{(}\textup{$\omega$
has no infinite component}\Big{)}=1\bigg{\\}}$

is the critical probability for Bernoulli bond percolation on $G$.

######  Theorem 1.

Suppose that $G_{n}=(V_{n},E_{n})$ is an expanding sequence of finite graphs
such that $G_{n}\to_{\textit{lwc}}(G,\rho)$ for some locally finite random
infinite rooted graph $(G,\rho)$ with $\mathbb{E}[\deg\rho]<\infty$. If
$\omega(G_{n})$ denotes edge percolation on $G_{n}$ with probability $p$, then
the following hold:

  * (i)

The critical edge percolation probability $p_{c}(G)$ is almost surely
constant.

  * (ii)

If $p<p_{c}(G)$, then for any $\alpha>0$, the probability that the largest
component in $\omega(G_{n})$ has size more than $\alpha|V_{n}|$ tends to $0$
as $n\rightarrow\infty$.

  * (iii)

If $p>p_{c}(G)$, then there exists $\alpha>0$ such that the probability that
the largest component in $\omega(G_{n})$ has size more than $\alpha|V_{n}|$
tends to $1$ as $n\rightarrow\infty$.

Items (i)–(iii) can be summarized by saying that “$p_{c}((G_{n})_{n\geq
0})=p_{c}(G)$.”

Theorem 1 removes the uniform degree bound assumption from [BNP11] and
[Sar21]. Further, [BNP11, Sar21] assume uniform expansion
($h_{\delta}(G_{n})\geq c_{\delta}$ for all $n\geq 1$), while Theorem 1 uses
only the weaker expansion assumption of Definition 3 ($h_{\delta}(G_{n})\geq
c_{\delta}$ for all $n\geq n_{\delta}$). In fact, our argument for the
unbounded degree setting involves pruning high-degree vertices, which can
break a uniform expansion assumption.

Additionally, we remark that the locality result can fail without some kind of
expansion condition. As an example,111We learned of this example from Elchanan
Mossel. let $h_{n}\to 0$, and let $H$ be a graph on $d/h_{n}$ vertices,
uniformly random subject to the requirement that all vertices have degree $2d$
except for one vertex $v_{*}(H)$ of degree $d$. Let $H_{1},\ldots,H_{n}$ be
disjoint copies of $H$, and let $J$ be a random $d$-regular graph on the
vertices $\\{v(H_{1}),\ldots,v(H_{n})\\}$. Let $G_{n}$ be the union of
$H_{1},\ldots,H_{n},J$, so it is a $(2d)$-regular graph. It has poor
expansion: if we take $A=H_{1}\cup\ldots\cup H_{k}$, then there are at most
$kd$ edges between $A$ and its complement, so

$\frac{|E_{G}(A,V\setminus A)|}{|A|}\leq\frac{kd}{kd/h_{n}}=h_{n}\,.$

Then $G_{n}$ converges locally weakly to the $(2d)$-regular tree $G$, which
has critical probability $p_{c}(G)=1/(2d-1)$. On the other hand, the sequence
of finite graphs has critical probability $p_{c}((G_{n})_{n\geq 0})\asymp
1/(d-1)$, as the emergence of a giant component on $G_{n}$ is determined by
percolation on $J$.

###  1.2. Related work

The locality question, as formulated by [BNP11], was motivated by a related
conjecture of Oded Schramm (see [BNP11, Conjecture 1.2]) for the setting of
infinite, vertex transitive graphs $G_{n}$. This conjecture has since been
proven in a few special cases, e.g., for uniformly nonemenable transitive
graphs [SXZ14], for abelian Cayley graphs [MT17], for transitive graphs with
exponential growth [Hut20], and for transitive graphs with exponential growth
[CMT22].

Percolation on general finite graphs, initiated by [ABS04], has also been
well-studied. Beyond the existence of a percolation threshold for sequences of
graphs [Fou07, Sha13], work has also been done to address the uniqueness of
the giant component [ABS04, EH21] and scaling windows around the critical
probability [NP10, BCvdH+05]. We remark that many of these results which give
detailed descriptions of behavior above, below, and near the threshold are
about specific graph sequences such as random regular graphs or rely on strong
assumptions about the graphs such as having uniformly bounded degrees. On the
other hand, the results from [BNP11] and our results apply more generally but
only show the existence of a threshold for locally weakly convergent graph
sequences.

In the case that $G_{n}$ converges locally weakly to a (random, rooted) tree
$G$, more can be said about the percolation threshold: it is well known that
$p_{c}(G)$ is the reciprocal of the tree branching number
$\operatorname{br}(G)$ [Lyo90]. More recently, physicists have proposed
[KNZ14] that $p_{c}((G_{n})_{n\geq 0})$ can be estimated by the reciprocal of
$\lambda^{\textup{NB}}(G_{n})$, the top eigenvalue of the nonbacktracking
matrix of $G_{n}$. Our main result Theorem 1, for the case the local limit $G$
is a tree, implies $p_{c}((G_{n})_{n\geq 0})=1/\operatorname{br}(G)$. This
naturally suggests that one approach for proving the conjecture of [KNZ14]
would be to relate the asympotics of $\lambda(B_{n})$ to the branching number
of $G_{n}$.

###  1.3. Organization

In Section 2 we prove some preliminary consequences of the local weak
convergence assumption, and give our pruning argument. In Section 3 we
complete the proof of Theorem 1.

### Acknowledgements

We thank Elchanan Mossel, Lenka Zdeborová, Jiaoyang Huang, and Fan Wei for
helpful conversations. Research supported in part by NSF CAREER grant
DMS-1940092.

##  2\. Preliminaries

We fix some basic notation to be used throughout. In a graph $H$, let
$d_{H}(u,v)$ be the distance between vertices $u$ and $v$; we will sometimes
write simply $d(u,v)$ if $H$ is clear from context. Recall that $B_{R}(v;H)$
denotes the $R$-neighborhood of $v$ in $H$, i.e., the subgraph of $H$ induced
by all vertices within distance $R$ of $v$. More generally, if $S$ is any
subset of vertices in $H$, we write $B_{R}(S;H)$ for the subgraph of $H$
induced by all vertices within distance $R$ of $S$. We write $\partial
B_{R}(v;H)$ for the set of vertices in $H$ at exactly distance $R$ from $v$.

We will also be working with a few different probability measures throughout
the proof. Let $\mathbb{P}_{p}$ and $\mathbb{P}_{n,p}$ denote edge percolation
on $G$ and $G_{n}$, respectively, with probability $p$. We use $\omega(G)$ and
$\omega(G_{n})$ to denote the percolation configurations. Let $\mathbf{L}_{n}$
be the law of a uniform random vertex $U_{n}$ of $G_{n}$, and let $\mu$ be the
law of $(G,\rho)$.

###  2.1. Continuity lemma

It follows directly from the definition of local weak convergence (Definition
2) that if $f$ is any bounded continuous function on the space
$\mathscr{G}_{\star}$ (Definition 1), and $G_{n}\to_{\textit{lwc}}(G,\rho)$,
then

$\mathbb{E}_{\mathbf{L}_{n}}f(G_{n},U_{n})\stackrel{{\scriptstyle
n\to\infty}}{{\longrightarrow}}\int f(G,\rho)\,d\mu\,.$ (2.1)

An immediate consequence is the following:

######  Lemma 2.

Suppose $f(G,\rho)$ is a bounded function on $\mathscr{G}_{\star}$ that
depends only on $B_{R}(\rho;G)$ for some finite $R$. Then $f$ is continuous on
$\mathscr{G}_{\star}$, so (2.1) holds whenever
$G_{n}\to_{\textit{lwc}}(G,\rho)$.

###### Proof.

Recall from Definition 1 that a metric on $\mathscr{G}_{\star}$ is given by
$\bm{d}$ from (1.1). If

$\bm{d}\Big{(}(G_{1},\rho_{1}),(G_{2},\rho_{2})\Big{)}\leq\frac{1}{1+R}$

then $B_{R}(\rho_{1};G_{1})$ and $B_{R}(\rho_{2};G_{2})$ agree, in which case
the assumption on $f$ implies $f(G_{1},\rho_{1})=f(G_{2},\rho_{2})$. This
shows that $f$ is (uniformly) continuous on $\mathscr{G}_{\star}$. ∎

######  Corollary 3.

If $G_{n}\to_{\textit{lwc}}(G,\rho)$, then for any measurable event $E$ we
have

$\mathbf{L}_{n}\otimes\mathbb{P}_{n,p}\bigg{(}\omega(B_{R}(U_{n};G_{n}))\in
E\bigg{)}\stackrel{{\scriptstyle
n\to\infty}}{{\longrightarrow}}\mu\otimes\mathbb{P}_{p}\bigg{(}\omega(B_{R}(\rho;G))\in
E\bigg{)}\,.$

###### Proof.

Define a function $f$ on $\mathscr{G}_{\star}$ by setting

$f(G,\rho)=\mathbb{P}_{p}\bigg{(}\omega(B_{R}(\rho;G))\in E\bigg{)}\,.$

Then $f$ clearly depends only on $B_{R}(\rho;G)$, so the claim follows from
Lemma 2. ∎

Lemma 2 and Corollary 3 will be used repeatedly throughout the proofs below.

###  2.2. Pruning

We will assume that the average degree of $G_{n}$ is uniformly bounded, i.e.,
that there exists a constant $d$ such that the average degree of $G_{n}$ is at
most $d$ for all $n$. In the case that this does not hold, the result can be
recovered from the former case via a pruning argument, which we now give. The
main result of this subsection is the following:

######  Proposition 4.

Suppose $(G,\rho)$ is a (random) infinite rooted graph, locally finite, with
$\mathbb{E}(\deg\rho)<\infty$. Suppose $G_{n}$ is an expanding sequence of
finite graphs with $G_{n}\to_{\textit{lwc}}(G,\rho)$. Then, for any
$\varepsilon>0$, we can choose a sequence of subgraphs $\bar{G}_{n}$ of
$G_{n}$ such that (i) each $\bar{G}_{n}$ has average degree at most
$\mathbb{E}(\deg\rho)+\varepsilon$; and (ii) the $\bar{G}_{n}$ also form an
expanding sequence, with $\bar{G}_{n}\to_{\textit{lwc}}(G,\rho)$.

######  Lemma 5.

Suppose $(G,\rho)$ is a (random) infinite rooted graph, locally finite, with
$\mathbb{E}(\deg\rho)<\infty$. Suppose $G_{n}$ is an expanding sequence of
finite graphs with $G_{n}\to_{\textit{lwc}}(G,\rho)$. Let $k_{n}$ be any
sequence of integers with $k_{n}\to\infty$, and let $\bar{G}_{n}$ be the
subgraph of $G_{n}$ formed by removing all edges incident to vertices of
degree at least $k_{n}$. Then $\bar{G}_{n}$ is also an expanding sequence with
$\bar{G}_{n}\to_{\textit{lwc}}(G,\rho)$.

###### Proof.

Let $S_{n}$ denote the subset of all vertices in $G_{n}$ with degree at least
$k_{n}$. We first claim that for any fixed $R<\infty$,

$\lim_{n\to\infty}\frac{|B_{R}(S_{n};G_{n})|}{|V_{n}|}=0\,.$ (2.2)

Indeed, if a vertex is within distance $R$ of $S_{n}$, then it must have at
least $k_{n}$ vertices in its $(R+1)$-neighborhood, so

$\frac{|B_{R}(S_{n};G_{n})|}{|V_{n}|}\leq\mathbf{L}_{n}\bigg{(}\Big{|}B_{R+1}(U_{n};G_{n})\Big{|}\geq
k_{n}\bigg{)}\,.$

Since $k_{n}\to\infty$ by assumption, it follows that for any finite $k$ we
have

$\limsup_{n\to\infty}\mathbf{L}_{n}\bigg{(}\Big{|}B_{R+1}(U_{n};G_{n})\Big{|}\geq
k_{n}\bigg{)}\leq\limsup_{n\to\infty}\mathbf{L}_{n}\bigg{(}\Big{|}B_{R+1}(U_{n};G_{n})\Big{|}\geq
k\bigg{)}$

Finally, by the assumption $G_{n}\to_{\textit{lwc}}(G,\rho)\sim\mu$ together
with Lemma 2, we have

$\lim_{n\to\infty}\mathbf{L}_{n}\bigg{(}\Big{|}B_{R+1}(U_{n};G_{n})\Big{|}\geq
k\bigg{)}=\mu\bigg{(}\Big{|}B_{R+1}(\rho,G)\Big{|}\geq k\bigg{)}\,.$

Combining the above calculations gives

$\limsup_{n\to\infty}\frac{|B_{R}(S_{n};G_{n})|}{|V_{n}|}\leq\mu\bigg{(}\Big{|}B_{R+1}(\rho,G)\Big{|}\geq
k\bigg{)}$

for any finite $k$. In the above, the right-hand side can be made arbitrarily
small by taking $k\to\infty$, while the left-hand side does not depend on $k$,
so this proves the claim (2.2).

We next argue that $\bar{G}_{n}\to_{\textit{lwc}}(G,\rho)$. The graphs $G_{n}$
and $\bar{G}_{n}$ have the same vertex set, so if $U_{n}$ is a uniformly
random vertex in $G_{n}$, then it is also a uniformly random vertex in
$\bar{G}_{n}$. Moreover, it is clear that the subgraphs $B_{R}(U_{n};G_{n})$
and $B_{R}(U_{n};\bar{G}_{n})$ are the same as long as $U_{n}$ does not belong
to $B_{R}(S_{n};G_{n})$. We see from (2.2) that $U_{n}$ lies in
$B_{R}(S_{n};G_{n})$ with probability $o_{n}(1)$, so we conclude that
$B_{R}(U_{n};\bar{G}_{n})$ also converges in law to $B_{R}(\rho;G)$. This
proves that $\bar{G}_{n}\to_{\textit{lwc}}(G,\rho)$, as claimed.

We now argue that $\bar{G}_{n}$ is an expanding sequence. Given $\delta>0$,
let $A\subseteq V_{n}$ with $\delta\leq|A|/|V_{n}|\leq 1/2$. Since the
original sequence $G_{n}$ is assumed to be expanding, there exists
$c_{\delta}>0$ such that

$\frac{|E_{G_{n}}(A,V_{n}\setminus A)|}{|A|}\geq c_{\delta}$

as long as $n\geq n_{\delta}$. The graph $\bar{G}_{n}$ differs from $G_{n}$
only in the deletion of edges incident to $S_{n}$, and the total number of
deleted edges can be upper bounded by $|B_{1}(S_{n};G_{n})|$. It follows using
(2.2) that

$\frac{|E_{\bar{G}_{n}}(A,V_{n}\setminus
A)|}{|A|}\geq\frac{|E_{{G}_{n}}(A,V_{n}\setminus
A)|}{|A|}-\frac{|B_{1}(S_{n};G_{n})|}{|V_{n}|}\frac{|V_{n}|}{|A|}\geq\frac{|E_{{G}_{n}}(A,V_{n}\setminus
A)|}{|A|}-\frac{o_{n}(1)}{\delta}\,,$

which is at least $c_{\delta}/2$ for $n$ large enough. This shows that
$\bar{G}_{n}$ is also an expanding sequence, and this concludes the proof of
the lemma. ∎

######  Lemma 6.

Let $D_{n}$ be a nonnegative integer-valued random variable converging in law
to an integrable random variable $D$. Then, for any $\varepsilon>0$, it is
possible to choose a sequence $k_{n}\to\infty$ such that
$\mathbb{E}(D_{n};D_{n}<k_{n})\leq\mathbb{E}D+\varepsilon<\infty$.

###### Proof.

First we note that for any finite integer $k$,

$\mathbb{E}(D;D<k)=\sum_{\ell\geq
1}\mathbb{P}\bigg{(}D\mathbf{1}\\{D<k\\}\geq\ell\bigg{)}=\sum_{\ell=1}^{k-1}\mathbb{P}(\ell\leq
D<k)\leq\sum_{\ell=1}^{k-1}\mathbb{P}(D\geq\ell)\,.$

It follows from the assumption that for any finite $\ell$ we have
$\mathbb{P}(D_{n}\geq\ell)\to\mathbb{P}(D\geq\ell)$ as $n\to\infty$, and
consequently
$\mathbb{P}(D_{n}\geq\ell)\leq\mathbb{P}(D\geq\ell)+\varepsilon/2^{\ell}$ for
all $n\geq m_{\ell}$, where $m_{\ell}$ is a finite integer depending on $\ell$
and $\varepsilon$. We then set $k_{n}\equiv n$ if
$\mathbb{P}(D_{n}\geq\ell)\leq\mathbb{P}(D\geq\ell)+\varepsilon/2^{\ell}$ for
all $\ell\geq 1$ and otherwise

$k_{n}\equiv\min\bigg{\\{}\ell\geq
1:\mathbb{P}(D_{n}\geq\ell)>\mathbb{P}(D\geq\ell)+\frac{\varepsilon}{2^{\ell}}\bigg{\\}}\geq\min\bigg{\\{}\ell\geq
1:m_{\ell}>n\bigg{\\}}.$

Our chosen $k_{n}$ must tend to infinity as $n\to\infty$ since the $m_{\ell}$
are finite. The definition of $k_{n}$ implies

$\mathbb{E}(D_{n};D_{n}<k_{n})\leq\sum_{\ell=1}^{k_{n}-1}\mathbb{P}(D_{n}\geq\ell)\leq\sum_{\ell=1}^{k_{n}-1}\bigg{\\{}\mathbb{P}(D\geq\ell)+\frac{\varepsilon}{2^{\ell}}\bigg{\\}}\leq\mathbb{E}D+\varepsilon\,,$

as required. This finishes the proof. ∎

###### Proof of Proposition 4.

This follows immediately by combining Lemmas 5 and 6. ∎

##  3\. Proof of Theorem 1

In this section we give the proof of Theorem 1. We prove parts (i), (ii), and
(iii) respectively in §3.1, §3.2, and §3.3. As noted above, our proofs are
based on arguments of [BNP11, Sar21], with modifications to handle unbounded
degrees. With the aim of being self-contained, we present a full proof below,
but will repeatedly point out similarities and differences with the prior
works.

###  3.1. Rerooting invariant functions

The proof of Theorem 1 part (i), after applying Proposition 4, is very similar
to the proof of the corresponding result in [Sar21], but we have made some
simplifications to the main lemma. As in [Sar21], the proof shows that any
rerooting invariant function of $G$ is almost surely constant, so the local
weak limit of an expanding sequence of finite graphs with finite expected root
degree is in fact an extremal unimodular graph (see Corollary 11 below).

######  Proposition 7.

Suppose that $\mathcal{G}_{1}$ and $\mathcal{G}_{2}$ are disjoint subsets of
$\mathscr{G}_{\star}$ such that each $\mathcal{G}_{i}$ is closed under
rerooting. Then, for any compact subsets $H_{i}\subseteq\mathcal{G}_{i}$ and
any finite integer $K$, there exists a finite integer $R$ such that for any
$(G_{1},\rho_{1})\in H_{1}$ and $(G_{2},\rho_{2})\in H_{2}$, we have that
$B_{R}(u_{1},G_{1})$ is not isomorphic to $B_{R}(u_{2},G_{2})$ for all
$u_{1}\in B_{K}(\rho_{1},G_{1}),u_{2}\in B_{K}(\rho_{2},G_{2})$.

In other words, given a rooted graph from one of $\mathcal{G}_{1}$ and
$\mathcal{G}_{2}$, we can determine which set it came from by looking at the
$R$-neighborhood of an arbitrary vertex in the $K$-neighborhood of the root.
This was proven in [Sar21] by considering finite open covers of the $H_{i}$.
We now give a conceptually cleaner proof of this fact.

######  Lemma 8.

Recall the metric $\bm{d}((G_{1}.\rho_{1}),(G_{2},\rho_{2}))$ on
$\mathscr{G}_{\star}$ defined by (1.1). For any finite $K$, let
$D_{K}:\mathscr{G}_{\star}\times\mathscr{G}_{\star}\to[0,1]$ be defined by

$D_{K}\Big{(}(G_{1},\rho_{1}),(G_{2},\rho_{2})\Big{)}\equiv\min\bigg{\\{}\bm{d}\Big{(}(G_{1},u_{1}),(G_{2},u_{2})\Big{)}:u_{i}\in
B_{K}(\rho_{i};G_{i})\bigg{\\}}\,,$

where $(G_{i},u_{i})$ denotes the graph $(G_{i},\rho_{i})$ rerooted at
$u_{i}$. The function $D_{K}$ is continuous.

###### Proof.

Since the graphs $G_{i}$ are locally finite, the definition of $D_{K}$ is a
minimum over finitely many quantities. Therefore we can choose $u_{i\star}\in
B_{K}(\rho_{i};G_{i})$ such that

$D_{K}\Big{(}(G_{1},\rho_{1}),(G_{2},\rho_{2})\Big{)}=\bm{d}\Big{(}(G_{1},u_{1\star}),(G_{2},u_{2\star})\Big{)}\,.$
(3.1)

Now suppose that $(H,\sigma)$ is any element of $\mathscr{G}_{\star}$ that
lies within distance $\delta=1/(1+R)$ of $(G_{2},\rho_{2})$, where we assume
$R\geq K$. From the definition (1.1), it implies that there is a rooted
isomorphism $\varphi$ that maps $B_{R}(\rho_{2};G_{2})$ to $B_{R}(\sigma;H)$.
Let $v_{2\star}\equiv\varphi(u_{2\star})$. We then have

$\displaystyle
D_{K}\Big{(}(G_{1},\rho_{1}),(H,\sigma)\Big{)}\leq\bm{d}\Big{(}(G_{1},u_{1\star}),(H,v_{2\star})\Big{)}$
$\displaystyle\qquad\leq\bm{d}\Big{(}(G_{1},u_{1\star}),(G_{2},u_{2\star})\Big{)}+\bm{d}\Big{(}(G_{2},u_{2\star}),(H,v_{2\star})\Big{)}$
$\displaystyle\qquad\leq
D_{K}\Big{(}(G_{1},\rho_{1}),(G_{2},\rho_{2})\Big{)}+\frac{1}{1+(R-K)}$

where the first step is by the definition of $D_{K}$, the second step is by
the triangle inequality, and the last step is by (3.1) together with the
observation that $B_{R-K}(u_{2\star};G_{2})$ must be isomorphic to
$B_{R-K}(u_{2\star};H)$. By exchanging the roles of $(G_{2},\rho_{2})$ and
$(H,\sigma)$ we must also have

$D_{K}\Big{(}(G_{1},\rho_{1}),(G_{2},\rho_{2})\Big{)}\leq
D_{K}\Big{(}(G_{1},\rho_{1}),(H,\sigma)\Big{)}+\frac{1}{1+(R-K)}\,,$

and combining the last two inequalities gives

$\bigg{|}D_{K}\Big{(}(G_{1},\rho_{1}),(G_{2},\rho_{2})\Big{)}-D_{K}\Big{(}(G_{1},\rho_{1}),(H,\sigma)\Big{)}\bigg{|}\leq\frac{1}{1+(R-K)}\,.$

This shows that $D_{K}:\mathscr{G}_{\star}\times\mathscr{G}_{\star}\to[0,1]$
is uniformly continuous in the second coordinate. By symmetry it is also
uniformly continuous in the first coordinate, which proves the claim. ∎

###### Proof of Proposition 7.

By the assumption that the $\mathcal{G}_{i}$ are disjoint and closed under
rerooting, the function $D_{K}$ must be strictly positive on
$\mathscr{G}_{1}\times\mathscr{G}_{2}$, and hence also on $H_{1}\times H_{2}$.
Since $D_{K}$ is continuous by Lemma 8, we conclude that the minimum of
$D_{K}$ on the compact set $H_{1}\times H_{2}$ must be lower bounded by some
$\delta>0$. Then, as long as we have $\delta>1/(1+R)$, it follows from the
definition of $D_{K}$ that for any $(G_{i},\rho_{i})\in H_{i}$, the graphs
$B_{R}(u_{1};G_{1})$ and $B_{R}(u_{2};G_{2})$ must be non-isomorphic for all
$u_{i}\in B_{K}(\rho_{i};G_{i})$. ∎

As in [Sar21], the rest of the proof proceeds by showing that if $p_{c}(G)$
can take on different values, then the root vertices in $G_{n}$ that produce
these different values in the limit can be connected by a short path,
contradicting Proposition 7. This relies crucially on the expansion condition,
which will be used in the application of the following classical result:

######  Theorem 9 (Menger’s theorem).

In a finite graph $G=(V,E)$, for any two disjoint subsets $A,B\subseteq V$,
the minimum size of an $A$-$B$ edge cut equals the maximum number of pairwise
edge-disjoint paths from $A$ to $B$.

######  Corollary 10.

Suppose $G_{n}$ is an expanding sequence of finite graphs with
$G_{n}\to_{\textit{lwc}}(G,\rho)$, and such that each $G_{n}$ has average
degree at most $d<\infty$. Then $p_{c}(G)$ is almost surely constant.

###### Proof.

This proof is essentially identical to that of [Sar21, Theorem 1.1]. The
uniform degree bound is only used to bound the number of edges in the graph,
so the argument can be directly adapted to when there is a uniform average
degree bound.

Suppose to the contrary that $p_{c}(G)$ is not almost surely constant. Then
there exist $0\leq a<b\leq 1$ such that, if we define

$\mathcal{G}^{1}=\Big{\\{}G:p_{c}(G)\leq
a\Big{\\}}\,,\quad\mathcal{G}^{2}=\Big{\\{}G:p_{c}(G)\geq b\Big{\\}}\,,$

then we will have a positive constant $p_{0}$ such that

$\min\Big{\\{}\mu(\mathcal{G}^{1}),\mu(\mathcal{G}^{2})\Big{\\}}\geq p_{0}\,.$

It is clear that $\mathcal{G}^{1}$ and $\mathcal{G}^{2}$ are disjoint, and
invariant under rerooting. Recall that $\mathscr{G}_{\star}$ is a Polish space
(see Definition 1), and that all Borel measures on Polish spaces are inner
regular. It follows that there exist compact subsets
$H^{i}\subseteq\mathcal{G}^{i}$ such that

$\min\Big{\\{}\mu(H^{1}),\mu(H^{2})\Big{\\}}\geq\frac{p_{0}}{2}\,.$ (3.2)

Recall that $G_{n}$ is assumed to be an expanding sequence in the sense of
Definition 3, and let $h$ be a positive constant such that

$\frac{1}{2}\liminf_{n\to\infty}h_{p_{0}/4}(G_{n})\geq h\,.$ (3.3)

Recall that $G_{n}$ is also assumed to have average degree at most $d$, and
define

$K=\frac{4d}{hp_{0}}\,.$

Now apply Proposition 7 with $H^{1}$, $H^{2}$, and $K$ as defined above, and
let $R$ be the integer that results from the conclusion of the proposition. If
we define the events

$E^{i}\equiv\bigg{\\{}(G,\rho):B_{R+K}(\rho;G)\cong
B_{R+K}(\rho^{i};G^{i})\textup{ for some }(G^{i},\rho^{i})\in
H^{i}\bigg{\\}}\,,$

then $E^{1}$ and $E^{2}$ are disjoint, because Proposition 7 guarantees in
particular that $B_{R}(\rho^{1};G^{1})\not\cong B_{R}(\rho^{2};G^{2})$ for all
$(G^{i},\rho^{i})\in H^{i}$. Let $A^{i}$ be the set of vertices in $G_{n}$
such that $(G_{n},v)$ belongs to event $E^{i}$. Since the events $E^{1}$ and
$E^{2}$ are disjoint, the vertex sets $A^{1}$ and $A^{1}$ are disjoint. Since
membership in $E^{i}$ depends only on the $(R+K)$-neighborhood of the root
vertex, Lemma 2 gives

$\lim_{n\to\infty}\frac{|A^{i}|}{|V_{n}|}=\lim_{n\to\infty}\mathbf{L}_{n}\Big{(}(G_{n},U_{n})\in
E^{i}\Big{)}=\mu\Big{(}(G,\rho)\in
E^{i}\Big{)}\geq\mu(H^{i})\stackrel{{\scriptstyle\eqref{e:choice.of.H}}}{{\geq}}\frac{p_{0}}{2}>0\,.$

Thus, for large enough $n$ we have $|A^{i}|/|V_{n}|\geq p_{0}/4$. It follows
from the expansion assumption (3.3) that, again for $n$ large enough, the
minimum edge cut separating $A^{1}$ from $A^{2}$ has at least
$\kappa_{n}=hp_{0}|V_{n}|/4$ edges. It follows by Menger’s theorem (Theorem 9)
that there are at least $\kappa_{n}$ edge-disjoint paths joining $A^{1}$ to
$A^{2}$. Since $G_{n}$ has at most $|V_{n}|d/2$ edges in total, at least half
of the paths must have length at most

$\frac{|V_{n}|d/2}{\kappa_{n}/2}=\frac{4d}{hp_{0}}=K\,.$

However, even the existence of a single such path joining $v^{1}\in A^{1}$ to
$v^{2}\in A^{2}$ results in a contradiction: it implies that the
$R$-neighborhood of $v^{2}$ lies inside the $(R+K)$-neighborhood of $v^{1}$ in
$G_{n}$. The definition of $A^{i}$ implies $B_{R+K}(v^{i};G_{n})\cong
B_{R+K}(\rho^{i};G^{i})$ for some $(G^{i},\rho^{i})\in H_{i}$. This
contradicts the definition of $R$ and concludes the proof. ∎

We remark that the above argument implies the result of [Sar21, Theorem 1.1]
under the milder assumptions of Theorem 1. Any local weak limit of a sequence
of finite graphs is a random rooted graph whose measure is unimodular (see
[AL07, Definition 2.1]). The space of unimodular probability measures on
$\mathscr{G}_{\star}$ is convex, and its extreme points are said to be
extremal or ergodic. A unimodular measure $\mu$ on $\mathscr{G}_{\star}$ is
extremal if and only if $\mu(E)\in\\{0,1\\}$ for any event $E$ that is
invariant under non-rooted isomorphisms ([AL07, Theorem 4.7]). It then follows
from the above that the local weak limit of an expanding sequence of finite
graphs is extremal if its root has finite expected degree:

######  Corollary 11.

Under the conditions of Theorem 1, $G$ is an extremal unimodular random graph.

###### Proof of Theorem 1 part (i).

By Proposition 4, there exists a expanding sequence of subgraphs $\bar{G}_{n}$
of $G_{n}$ such that each $\bar{G}_{n}$ has average degree at most
$\mathbb{E}(\deg\rho)+1$ and $\bar{G}_{n}\to_{\textit{lwc}}(G,\rho)$. Applying
Corollary 10 to $\bar{G}_{n}$ gives the result. ∎

###  3.2. Subcritical percolation regime

In this subsection we prove Theorem 1 part (ii), which says in short that
$p_{c}(G_{n})\geq p_{c}(G)$. The result is a straightforward consequence of
Lemmas 2 and Corollary 3. Let us point out that for this part of the proof we
work with the original sequence $G_{n}$, without pruning, so we do not assume
that $G_{n}$ has bounded average degree. The reason for this is that (a
priori) we only know that $p_{c}(\bar{G}_{n})\geq p_{c}(G_{n})$, so it is not
sufficient to only prove $p_{c}(\bar{G}_{n})\geq p_{c}(G)$.

Recall that we let $\omega=\omega(G)$ denote an edge percolation configuration
on graph $G$; we also view $\omega$ as a (random) subgraph of $G$. Let
$C(v)\equiv C(v;\omega)$ denote the connected component of $\omega$ that
contains $v$. For any vertex subsets $A$ and $B$, let
$A\leftrightarrow_{\omega}B$ indicate that $A$ and $B$ are connected by a path
using only edges in $\omega$. We write simply $A\leftrightarrow B$ when
$\omega$ is clear from context.

###### Proof of Theorem 1 part (ii).

The proof of [Sar21, Corollary 1.2] (which is based on the proof of [BNP11,
Theorem 1.3]) applies here verbatim. To be self-contained, we also give the
argument here. Under the assumptions of Theorem 1, we will show that for any
$p<p_{c}(G)$ and any $\alpha>0$, the probability that the largest component in
$\omega(G_{n})$ contains more than $\alpha|V_{n}|$ vertices tends to zero as
$n\to\infty$. Let $C_{n}$ denote the vertices in the largest component of
$\omega(G_{n})$. For any $R>0$, we have the chain of inequalities

$\displaystyle\alpha\mathbb{P}_{n,p}\Big{(}|C_{n}|\geq\alpha|V_{n}|\Big{)}$
$\displaystyle\qquad\leq\mathbf{L}_{n}\otimes\mathbb{P}_{n,p}\Big{(}|C_{n}|\geq\alpha|V_{n}|,U_{n}\in
C_{n}\Big{)}$
$\displaystyle\qquad\leq\mathbf{L}_{n}\otimes\mathbb{P}_{n,p}\Big{(}|C(U_{n})|\geq\alpha|V_{n}|\Big{)}$
$\displaystyle\qquad\leq\mathbf{L}_{n}\otimes\mathbb{P}_{n,p}\Big{(}|C(U_{n})|\geq|B_{R}(U_{n},G_{n})|\Big{)}+\mathbf{L}_{n}\Big{(}|B_{R}(U_{n},G_{n})|\geq\alpha|V_{n}|\Big{)}$
$\displaystyle\qquad\leq\mathbf{L}_{n}\otimes\mathbb{P}_{n,p}\Big{(}U_{n}\leftrightarrow\partial
B_{R}(U_{n},G_{n})\Big{)}+\mathbf{L}_{n}\Big{(}|B_{R}(U_{n},G_{n})|\geq\alpha|V_{n}|\Big{)}\,.$
(3.4)

We will consider separately the two terms on the right-hand side of (3.4). For
the first term, Corollary 3 directly implies

$\lim_{n\to\infty}\mathbf{L}_{n}\otimes\mathbb{P}_{n,p}\Big{(}U_{n}\leftrightarrow\partial
B_{R}(U_{n},G_{n})\Big{)}=\mu\otimes\mathbb{P}_{p}\Big{(}\rho\leftrightarrow\partial
B_{R}(\rho,G)\Big{)}\,.$

For the second term on the right-hand side of (3.4), for any finite $K$ we
have

$\limsup_{n\to\infty}\mathbf{L}_{n}\Big{(}|B_{R}(U_{n};G_{n})|\geq\alpha|V_{n}|\Big{)}\leq\lim_{n\to\infty}\mathbf{L}_{n}\Big{(}|B_{R}(U_{n};G_{n})|\geq
K\Big{)}=\mu\Big{(}|B_{R}(\rho;G)|\geq K\Big{)}\,,$

where the last step is by Lemma 2. Since $|B_{R}(\rho;G)|$ is almost surely
finite, the right-hand side tends to zero as $K\to\infty$. The left-hand side
does not depend on $K$, so we conclude that in fact

$\lim_{n\to\infty}\mathbf{L}_{n}\Big{(}|B_{R}(U_{n};G_{n})|\geq\alpha|V_{n}|\Big{)}=0\,.$

Substituting these limits back into (3.4) gives

$\limsup_{n\to\infty}\alpha\mathbb{P}_{n,p}\Big{(}|C_{n}|\geq\alpha|V_{n}|\Big{)}\leq\mu\otimes\mathbb{P}_{p}\Big{(}\rho\leftrightarrow\partial
B_{R}(\rho,G)\Big{)}\,.$

Since $p<p_{c}(G)$, the right-hand side tends to zero as $R\to\infty$. The
left-hand side does not depend on $R$, so we finally conclude that

$\lim_{n\to\infty}\mathbb{P}_{n,p}\Big{(}|C_{n}|\geq\alpha|V_{n}|\Big{)}=0\,,$

as claimed. ∎

###  3.3. Supercritical percolation regime

In this subsection we prove Theorem 1 part (iii), which says in short that
$p_{c}(G_{n})\leq p_{c}(G)$. We make use of Proposition 4 to reduce to the
case where the $G_{n}$ have uniformly bounded average degree. The rest of the
proof is adapted from that of [BNP11, Theorem 1.3], where the idea is to first
percolate with a smaller probability $p-\varepsilon>p_{c}(G)$, then connect
the resulting clusters into a linear-sized component with the remaining
probability. While [BNP11] used the bounded degree assumption to obtain a
strong concentration bound for the number of percolation clusters, we will use
Chebychev’s inequality to obtain a weaker (but sufficient) bound for our more
general setting. This is based on the following lemma:

######  Lemma 12.

Under the assumptions of Theorem 1, we have for any fixed $R$ that

$\lim_{n\to\infty}\frac{1}{|V_{n}|^{2}}\bigg{|}\Big{\\{}(u,v)\in(V_{n})^{2}:d(u,v)\leq
R\Big{\\}}\bigg{|}=0\,.$

###### Proof.

We can bound the quantity of interest by

$\displaystyle r(n)$
$\displaystyle\equiv\frac{1}{|V_{n}|^{2}}\bigg{|}\Big{\\{}(u,v)\in(V_{n})^{2}:d(u,v)\leq
R\Big{\\}}\bigg{|}=\sum_{v\in V_{n}}\frac{|B_{R}(v;G_{n})|}{|V_{n}|^{2}}$
$\displaystyle=\frac{\mathbb{E}_{\mathbf{L}_{n}}|B_{R}(U_{n};G_{n})|}{|V_{n}|}\leq\frac{K}{|V_{n}|}+\mathbf{L}_{n}\Big{(}|B_{R}(U_{n},G_{n})|\geq
K\Big{)}\,,$

where the last inequality holds for any finite $K$. It follows using Lemma 2
that

$\limsup_{n\to\infty}r(n)\leq\mu\Big{(}|B_{R}(\rho,G)|\geq K\Big{)}\,.$

Since $|B_{R}(\rho,G)|$ is finite almost surely, the right-hand side tends to
zero as $K\to\infty$. The left-hand side does not depend on $K$, so we
conclude $r(n)\to 0$ as claimed.∎

The rest of the proof closely follows that of [BNP11, Theorem 1.3], which uses
a sprinkling argument (see also [AKS82, ABS04]). Given $p>p_{c}(G)$, let
$p_{1}>p_{c}(G)$ and $\varepsilon>0$ such that

$1-p=(1-p_{1})(1-\varepsilon)\,.$

Let $\omega_{1}$ denote an edge percolation configuration with edge
probability $p_{1}$, and let $\omega_{\varepsilon}$ denote an edge percolation
configuration with edge probability $\varepsilon$. Then the union of
$\omega_{1}$ and $\omega_{\varepsilon}$ is equidistributed as a
$p$-percolation. The first step is to show that a positive proportion of the
vertices in $\omega_{1}$ percolate far and hence are contained in sizable
components.

######  Lemma 13.

If $p>p_{c}(G)$, then there exists $\delta>0$ such that for any $R>0$,

$\lim_{n\to\infty}\mathbb{P}_{n,p}\bigg{(}\frac{1}{|V_{n}|}\bigg{|}\Big{\\{}v\in
G_{n}:v\leftrightarrow\partial
B_{R}(v;G_{n})\Big{\\}}\bigg{|}\leq\delta\bigg{)}=0\,.$

###### Proof.

The first part of the proof establishes a positive lower bound on the
expectation of the quantity of interest and is essentially identical to the
argument in the proof of [Sar21, Corollary 1.2]. Given $(G,\rho)$, let
$f(G,\rho)$ be the probability under a $p$-percolation that the root
percolates:

$f(G,\rho)\equiv\mathbb{P}_{p}\bigg{(}\Big{|}C(\rho;\omega(G))\Big{|}=\infty\bigg{)}\,,$

where $C(\rho;\omega(G))$ is the connected component of $\omega(G)$ that
contains $\rho$. This is a measurable function on $\mathscr{G}_{\star}$, since
it is the decreasing limit of the local (hence measurable) functions

$f_{R}(G,\rho)\equiv\mathbb{P}_{p}\bigg{(}\rho\leftrightarrow\partial
B_{R}(\rho;G)\bigg{)}\,.$

Since $p>p_{c}(G)$, the quantity $f(G,\rho)$ must be positive, $\mu$-almost
surely. Averaging over $\mu$ gives

$\mu\otimes\mathbb{P}_{p}\Big{(}\rho\leftrightarrow\partial
B(\rho,R)\Big{)}\geq\mu\otimes\mathbb{P}_{p}\bigg{(}\Big{|}C(\rho;\omega(G))\Big{|}=\infty\bigg{)}\geq\delta>0$

uniformly over all $R$. It follows by Corollary 3 that

$\lim_{n\to\infty}\mathbf{L}_{n}\otimes\mathbb{P}_{n,p}\Big{(}U_{n}\leftrightarrow\partial
B_{R}(U_{n};G_{n})\Big{)}=\mu\otimes\mathbb{P}_{p}\Big{(}\rho\leftrightarrow\partial
B(\rho,R)\Big{)}\geq\delta\,.$

Let $Y_{v}=Y_{v}(R)$ be the indicator that $v$ percolates at least distance
$R$ under the $p$-percolation,

$Y_{v}=\mathbf{1}\Big{\\{}v\leftrightarrow\partial B_{R}(v;G_{n})\Big{\\}}\,.$

Let $X_{n}=X_{n}(R)$ be the sum of $Y_{v}$ over all $v\in V_{n}$, so $X_{n}$
is the quantity of interest. If $\mathbb{E}_{n,p}$ denotes expectation with
respect to $\mathbb{P}_{n,p}$, then for all $n$ large enough we have

$\frac{\mathbb{E}_{n,p}X_{n}}{|V_{n}|}=\mathbf{L}_{n}\otimes\mathbb{P}_{n,p}\Big{(}U_{n}\leftrightarrow\partial
B_{R}(U_{n};G_{n})\Big{)}\geq\frac{\delta}{2}\,.$

We next apply the second moment method to show that $X_{n}$ concentrates, as
the method used in [BNP11] no longer applies due to potential unboundedness of
the degrees. Note that $\operatorname{Cov}(Y_{u},Y_{v})\leq 1$ always (since
$Y_{u},Y_{v}$ are indicators). If $u$ and $v$ lie at distance more than $2R$
in $G_{n}$, then $Y_{u}$ and $Y_{v}$ are clearly independent, so in this case
$\operatorname{Cov}(Y_{u},Y_{v})=0$. Hence, we have that

$\frac{\operatorname{Var}_{n,p}X_{n}}{|V_{n}|^{2}}\leq\frac{1}{|V_{n}|^{2}}\bigg{|}\Big{\\{}(u,v)\in(V_{n})^{2}:d(u,v)\leq
2R\Big{\\}}\bigg{|}=o_{n}(1)\,,$

by Lemma 12. It follows by Chebyshev’s inequality that

$\mathbb{P}_{n,p_{1}}\bigg{(}X_{n}\leq\frac{\delta|V_{n}|}{4}\bigg{)}\to 0\,.$

The claim follows by redefining $\delta$. ∎

Once we prune the $G_{n}$, the remainder of the proof is essentially the same
as in [BNP11], where the uniform degree bound is only used to bound the number
of edges in the graph. For the sake of being self-contained, we reproduce the
following lemma from the end of the proof of [BNP11, Theorem 1.3]:

######  Lemma 14.

Suppose $G_{n}$ is an expanding sequence of finite graphs, each with average
degree at most $d$. Let $\omega_{1},\omega_{\varepsilon}$ be two independent
edge percolation configurations on $G_{n}$, with edge probabilities $p_{1}$
and $\varepsilon$ respectively. Let $X_{n}(R)$ be the number of vertices in
$G_{n}$ that percolate to distance $R$ under $\omega_{1}$. Let $C_{n}$ be the
largest connected component in $\omega_{1}\cup\omega$ and $\delta\in(0,1)$.
Then there is a large enough constant $R$ such that

$\lim_{n\to\infty}\mathbb{P}_{n,p}\bigg{(}|C_{n}|\geq\frac{\delta|V_{n}|}{3}\,\bigg{|}\,X_{n}(R)>\delta|V_{n}|\bigg{)}=1\,.$

The value of $R$ depends on $d$, $\delta$, $\varepsilon$, and the expansion
rate, but not on $n$.

###### Proof.

We condition on $\omega_{1}$ and use the randomness of $\omega_{\varepsilon}$
only. Consider the connected components of $\omega_{1}$ of size at least $R$,
of which there are at most $m=|V_{n}|/R$. By the bound on $X_{n}(R)$, these
components cover at least a $\delta$ fraction of the vertices in $G_{n}$. Let
$A,B$ be any partition of these components into two sets, such that

$\min\bigg{\\{}|A|,|B|\bigg{\\}}\geq\frac{\delta|V_{n}|}{3}\,.$ (3.5)

Let $h$ be a positive constant such that

$\frac{1}{2}\liminf_{n\rightarrow\infty}h_{\delta/3}(G_{n})\geq h.$

By the same reasoning as in the proof of Corollary 10, $A$ and $B$ must be
joined by at least

$\lambda_{n}=\frac{h\delta|V_{n}|}{3}$

edge-disjoint paths for sufficiently large $n$ , at least half of which must
be of length at most

$\frac{|V_{n}|d/2}{\lambda_{n}/2}=\frac{3d}{h\delta}\equiv L\,.$

Thus, the probability that $\omega_{\varepsilon}$ fails to have a path from
$A$ to $B$ is at most

$(1-\varepsilon^{L})^{\lambda_{n}/2}\leq\exp\bigg{\\{}-\frac{\lambda_{n}}{2}\varepsilon^{L}\bigg{\\}}$

Let $\bm{E}$ be the event that $\omega_{\varepsilon}$ has a path from $A$ to
$B$ for every partition $(A,B)$ satisfying (3.5). Then

$1-\mathbb{P}(\bm{E})\leq
2^{m}\exp\bigg{\\{}-\frac{\lambda_{n}}{2}\varepsilon^{L}\bigg{\\}}\leq\exp\bigg{\\{}|V_{n}|\bigg{(}\frac{1}{R}-\frac{h\delta}{6}\varepsilon^{L}\bigg{)}\bigg{\\}}\,,$

which tends to zero by taking $R$ a large enough constant. On the event
$\bm{E}$, the configuration $\omega=\omega_{1}\cup\omega_{\varepsilon}$
contains a component that occupies at least a $\delta/3$ fraction of the
vertices in $G_{n}$, as desired. ∎

###### Proof of Theorem 1 part (iii).

By Proposition 4, there exists a expanding sequence of subgraphs $\bar{G}_{n}$
of $G_{n}$ such that each $\bar{G}_{n}$ has average degree at most
$\mathbb{E}(\deg\rho)+1$ and $\bar{G}_{n}\to_{\textit{lwc}}(G,\rho)$. Since
$\bar{G}_{n}$ is a subgraph of $G_{n}$, the probability that $\omega(G_{n})$
contains a component of linear size is at least the probability that
$\omega(\bar{G}_{n})$ contains a component of linear size, i.e.
$p_{c}(\bar{G}_{n})\geq p_{c}(G_{n})$. Thus, it suffices to prove the
statement for $\bar{G}_{n}$.

Let $X_{n}(R)$ be the number of vertices in $\bar{G}_{n}$ that percolate to
distance $R$ under $\omega_{1}$ and $C_{n}$ be the largest connected component
in $\omega(\bar{G}_{n})$. We have that

$\mathbb{P}_{n,p}\bigg{(}|C_{n}|\geq\frac{\delta|V_{n}|}{3}\bigg{)}\geq\mathbb{P}_{n,p_{1}}(X_{n}(R)>\delta|V_{n}|)\times\mathbb{P}_{n,p}\bigg{(}|C_{n}|\geq\frac{\delta|V_{n}|}{3}\,\bigg{|}\,X_{n}(R)>\delta|V_{n}|\bigg{)}$

for any $R,\delta>0$. Taking $\delta$ to be the constant that results from
applying Lemma 13 to $\omega_{1}(\bar{G}_{n})$ and $R$ to be the constant that
results from applying Lemma 14 to $\bar{G}_{n}$ and $\delta$, we have that

$\lim_{n\rightarrow\infty}\mathbb{P}_{n,p_{1}}(X_{n}(R)>\delta|V_{n}|)=\lim_{n\rightarrow\infty}\mathbb{P}_{n,p}\bigg{(}|C_{n}|\geq\frac{\delta|V_{n}|}{3}\,\bigg{|}\,X_{n}(R)>\delta|V_{n}|\bigg{)}=1.$

It follows that

$\lim_{n\rightarrow\infty}\mathbb{P}_{n,p}\bigg{(}|C_{n}|\geq\frac{\delta|V_{n}|}{3}\bigg{)}=1,$

so we may take $\alpha=\delta/3$. ∎

## References

  * [ABS04] N. Alon, I. Benjamini, and A. Stacey.  Percolation on finite graphs and isoperimetric inequalities.  Ann. Probab., 32(3A):1727–1745, 2004.
  * [AKS82] M. Ajtai, J. Komlós, and E. Szemerédi.  Largest random component of a $k$-cube.  Combinatorica, 2(1):1–7, 1982.
  * [AL07] D. Aldous and R. Lyons.  Processes on unimodular random networks.  Electron. J. Probab., 12:no. 54, 1454–1508, 2007.
  * [BCvdH+05] C. Borgs, J. T. Chayes, R. van der Hofstad, G. Slade, and J. Spencer.  Random subgraphs of finite graphs. I. The scaling window under the triangle condition.  Random Structures Algorithms, 27(2):137–184, 2005.
  * [BNP11] I. Benjamini, A. Nachmias, and Y. Peres.  Is the critical percolation probability local?  Probab. Theory Related Fields, 149(1-2):261–269, 2011.
  * [CMT22] D. Contreras, S. Martineau, and V. Tassion.  Locality of percolation for graphs with polynomial growth.  arXiv:2205.10253, 2022.
  * [EH21] P. Easo and T. Hutchcroft.  Supercritical percolation on finite transitive graphs I: Uniqueness of the giant component.  arXiv:2112.12778, 2021.
  * [Fou07] N. Fountoulakis.  Percolation on sparse random graphs with given degree sequence.  Internet Math., 4(4):329–356, 2007.
  * [Hut20] T. Hutchcroft.  Locality of the critical probability for transitive graphs of exponential growth.  Ann. Probab., 48(3):1352–1371, 2020.
  * [KNZ14] B. Karrer, M. E. J. Newman, and L. Zdeborová.  Percolation on sparse networks.  Phys. Rev. Lett., 113:208702, Nov 2014.
  * [Lyo90] R. Lyons.  Random walks and percolation on trees.  Ann. Probab., 18(3):931–958, 1990.
  * [MT17] S. Martineau and V. Tassion.  Locality of percolation for Abelian Cayley graphs.  Ann. Probab., 45(2):1247–1277, 2017.
  * [NP10] A. Nachmias and Y. Peres.  Critical percolation on random regular graphs.  Random Structures Algorithms, 36(2):111–148, 2010.
  * [Sar21] S. Sarkar.  A note on the local weak limit of a sequence of expander graphs.  Electron. Commun. Probab., 26:Paper No. 32, 6, 2021.
  * [Sha13] Y. Shang.  The sharp threshold for percolation on expander graphs.  Math. Slovaca, 63(5):1141–1152, 2013.
  * [SXZ14] H. Song, K.-N. Xiang, and S.-C.-H. Zhu.  Locality of percolation critical probabilities: uniformly nonamenable case.  arXiv:1410.2453, 2014.


# Early Exploration of a Flexible Framework for Efficient Quantum Linear
Solvers in Power Systems
††thanks: The Pacific Northwest National Laboratory (PNNL) is operated by
Battelle for the U.S. Department of Energy (DOE) under Contract DE-
AC05-76RL01830. This work was supported by DOE Office of Electricity through
its Advanced Grid Modeling (AGM) program, DOE Office of Basic Energy Science
(BES), National Science Foundation CAREER, and Defense Advanced Research
Projects Agency. This research used resources of the National Energy Research
Scientific Computing Center (NERSC), a DOE Office of Science User Facility
located at Lawrence Berkeley National Laboratory. Emails:
<EMAIL_ADDRESS>and<EMAIL_ADDRESS>
Muqing Zheng12, Yousu Chen1, Xiu Yang2, Ang Li1  1 Pacific Northwest National
Laboratory, Richland, WA, 99352, USA
2Lehigh University, Bethlehem, PA 18015, USA

###### Abstract

The rapid integration of renewable energy resources presents formidable
challenges in managing power grids. While advanced computing and machine
learning techniques offer some solutions for accelerating grid modeling and
simulation, there remain complex problems that classical computers cannot
effectively address. Quantum computing, a promising technology, has the
potential to fundamentally transform how we manage power systems, especially
in scenarios with a higher proportion of renewable energy sources. One
critical aspect is solving large-scale linear systems of equations, crucial
for power system applications like power flow analysis, for which the Harrow-
Hassidim-Lloyd (HHL) algorithm is a well-known quantum solution. However, HHL
quantum circuits often exhibit excessive depth, making them impractical for
current Noisy-Intermediate-Scale-Quantum (NISQ) devices. In this paper, we
introduce a versatile framework, powered by NWQSim, that bridges the gap
between power system applications and quantum linear solvers available in
Qiskit. This framework empowers researchers to efficiently explore power
system applications using quantum linear solvers. Through innovative gate
fusion strategies, reduced circuit depth, and GPU acceleration, our simulator
significantly enhances resource efficiency. Power flow case studies have
demonstrated up to a eight-fold speedup compared to Qiskit Aer, all while
maintaining comparable levels of accuracy.

###### Index Terms:

quantum computing, Harrow-Hassidim-Lloyd, high-performance computing, grid
analytics

This is a preprint version of the work.
This work has been submitted to the IEEE for possible publication.
Copyright may be transferred without notice, after which this version may no
longer be accessible.

##  I Introduction

As nations around the world unite in the battle against climate change, clean
energy and decarbonization have taken center stage. For perspective, in 2021,
one-fifth of the electricity in the U.S. came from renewable sources. With
eyes set on a greener future, America is aiming for $100\%$ clean electricity
by 2035 and zero net emissions by 2050. California has a pathway to achieve
carbon neutral by 2045[1]. Europe is on a similar trajectory, with numerous
countries launching strategies to cut their greenhouse gas emissions
significantly.

The energy landscape is indeed transforming. The shift is evident with the
growing reliance on renewable sources like wind and solar, advanced energy-
consuming devices, and other cutting-edge technologies. Taking California as
an example: there’s been a rapid increase in the use of wind and solar energy.
By 2021, $90\%$ of its power came from these sources. And in a milestone
moment on May 8, 2022, this number surged past $100\%$. Such changes has
brought significant challenges to grid management as it demands a much faster,
quicker responses to a much complex, dynamic, and with uncertainty not been
seen before.

High-performance-computing, as an emerging technique, can help power system
engineers to gain faster computational speeds. However, there are still
problems in power systems, particularly those of large scale and high
granularity, are intractable to classical cluster machines due to their
computational complexity, data communication overhead, or both. Some examples
are listed below:

  * •

Stochastic Optimization for Operational Planning with Renewables: Given the
variability and uncertainty of renewable energy sources like wind and solar,
planning for their optimal integration and dispatch involves solving high-
dimensional stochastic optimization problems.

  * •

Very Large-Scale Contingency Analysis: Assessing the grid’s resilience to
numerous possible contingencies (like multiple line outages) requires
analyzing a combinatorial number of scenarios. It’s a problem that can quickly
grow beyond the capabilities of classical computing clusters.

  * •

Security Constrained Unit Commitment: Determining which power plants should be
turned on or off to meet the anticipated electricity demand, while also
ensuring the grid remains stable and secure against potential disturbances.

Classical cluster machines, while powerful, have limitations in terms of
processing capability and data communication bandwidth. Emerging computing
paradigms, like quantum computing, might offer pathways to tackle some of
these intractable problems in the future.

The basic information unit in quantum computing is called qubits. An $n$-qubit
quantum state can be described as a $2^{n}$-size complex vector in Hilbert
space. Quantum state has a unique feature: at most $2^{n}$ number of
distinguishable, i.e., orthogonal, quantum states can form a valid $n$-qubit
quantum state. This property is called superposition. Another property of
quantum states is entanglement. It is when several quantum states are
intricately intertwined so that any transformation on one state will affect
all the other states. Both are important sources of

quantum advantage.  This is a preprint version of the work.
This work has been submitted to the IEEE for possible publication.
Copyright may be transferred without notice, after which this version may no
longer be accessible.

Quantum computing holds great promise for power systems, even though some
details remain unclear. Many experts are eager to start exploring its
application before it is perfected. The areas that they are working on include
AC/DC power flow, contingency analysis, state estimation, transient stability
assessment, fault diagnosis, unit commitment, and facility location–allocation
problem, as shown in [2] [3], [4]. However, there are still gaps that need to
be bridged to enable power system engineers, who typically do not have a
quantum computing background, to evaluate how to use quantum algorithms to
accelerate their applications. To address this, we introduce a flexible
framework, empowered by Northwest Quantum Simulator (NWQSim)[5],[6], to bridge
the gaps between power system engineers and quantum techniques.

Building upon this potential of quantum computing for power systems, it is
crucial to understand the tools that can help harness this technology
effectively. Quantum solvers, for instance, leverage the principles of quantum
mechanics to address specific types of problems more efficiently than
classical computers. Some of the main quantum solvers include:

  * •

Quantum Linear System Algorithms (QLSAs): QLSAs are used for solving linear
systems of equations efficiently on a quantum computer, one of the
representative algorithms is Harrow-Hassidim-Lloyd (HHL) [7],

  * •

Variational Quantum Eigensolver (VQE): VQE is a quantum solver used for
finding the ground state energy (the smallest eigenvalue) of quantum systems,

  * •

Quantum Approximate Optimization Algorithm (QAOA): QAOA is a quantum solver
designed for solving combinatorial optimization problems, such as MaxCut
problems,

  * •

Quantum Annealing (QA): QA currently solves the problems that can be
reformulated to Ising models and quadratic unconstrained binary optimization
problems.

This paper start with integrating the HHL solver into the framework due to its
ability to handle linear system and its connection to power systems. HHL has
the potential to solve complex linear equations with exponential speedup
compared to classical solvers, which could enable a significant computational
improvement for many power system applications, such as power flow, short-
circuit analysis, contingency analysis, state estimation, and voltage
stability analysis [2, 3].

Nevertheless, applying the HHL algorithm directly to power system problems
presents several challenges. First and foremost, the quantum circuits for the
HHL algorithm tend to be very deep, even for small-scale problems. To
illustrate, consider the HHL circuit that only solves a random $2$-by-$2$
linear system, depicted in Figure 1. It already has 120 one-qubit gates and 90
two-qubit gates [8]. Furthermore, the small experiment illustrated in Figure 2
shows that significant errors persist even in a circuit with one-qubit gate
and a two-qubit gate. In essence, given the current noise level in Noisy
Intermediate-Scale Quantum (NISQ) devices, obtaining meaningful results for
such circuits is difficult. To the best of our knowledge, a $4$-by-$4$ system
is the largest one that can be accurately solved on an actual quantum computer
in existing literature [9], which falls short of the requirements for
practical power system applications.

Figure 1: An example of the transpiled HHL circuit that solves a random
$2$-by-$2$ linear system. It is generated by the package in [8] and transpiled
into one- and two-qubit basis gates. Figure 2: 10,000 measurements of bell
state
$\frac{1}{\sqrt{2}}\left(\left|00\right\rangle+\left|11\right\rangle\right)$
prepared from initial state $\left|00\right\rangle$ using a Hadamard gate and
a Controlled-NOT gate on IBM superconducting quantum computer ibmq_kolkata.
The noise in the quantum computer causes the actual results deviate from the
ideal case where $01$ and $10$ should not appear and the counts for $00$ and
$11$ should be almost equal in the measurements.

Additionally, HHL circuits have certain requirements on qubit connectivity.
This is not an issue for trapped-ion quantum computers as discussed in [9],
but it is challenging for superconducting quantum computers offered by
companies like IBM and Google. Lastly, there is a lack of software tools for
developing and deploying quantum algorithms for power system applications,
especially for non-quantum-background researchers and engineers.

Given all these challenges, the authors propose a flexible simulator-based
framework designed for tackling power system problems. This framework
leverages quantum simulators, which are dedicated to the evaluation of quantum
circuits on classical quantum computers. The use of simulators can facilitate
feasibility studies in this area, targeting more realistic problem settings
and making it more accessible to a broader range of researchers. Once fault-
tolerant quantum computers are available, researchers can then migrate their
existing work

to the real hardware, further advancing the field.  This is a preprint version
of the work.
This work has been submitted to the IEEE for possible publication.
Copyright may be transferred without notice, after which this version may no
longer be accessible.

The remainder of this paper is organized as follows: Section II introduces the
flexible framework along with its key components and setup. Section III
explains the details of applying the HHL algorithm to linear systems in power
flow problems. Section IV illustrates the use of our framework on three power
flow problems. Finally, Section V concludes the paper with a discussion of
future work.

##  II Introduction to the Framework

The motivation is to provide power system engineers with a versatile platform
for testing power system applications by seamlessly combining classical and
quantum algorithms within a hybrid simulation framework. This hybrid
simulation framework allows users to connect with various existing power
system applications with different quantum simulators via the interfaces that
facilitate the smooth exchange of data necessary for the successful execution
of hybrid simulations. The ultimate goal is to equip power system
professionals with the tools and environment they need to explore and advance
quantum computing in power system applications.

###  II-A The proposed hybrid framework

As illustrated in Figure 3, our conceptual framework serves as a bridge,
enhancing the communication between various quantum algorithms, such as HHL
and variational quantum algorithm (VQA), and external simulation tools
compatible with diverse languages and platforms. Our current focus lies in the
integration of the HHL algorithm with the power flow functions in MatPower
[10]. This integration necessitates some data preprocessing to transform the
data into a format suitable for HHL. Subsequently, we are poised to extend our
framework by integrating other quantum algorithms, such as VQA, to further
enrich our hybrid simulation capabilities.

Figure 3: The conceptual view of the proposed framework.

###  II-B Existing software stack choices

On the classical side of the framework, as a starting point, we have
integrated MatPower due to its widespread use in power system problems in
academia. For the quantum computing component of our framework, we have
selected Qiskit[11] developed by IBM, from among other tools such as Cirq from
Google, and Q# from Microsoft. Qiskit offers several advantages that make it a
preferred choice. Foremost is the large and active community supporting Qiskit
as an open-source python library. This helps lower the potential difficulties
on the development of the framework. Additionally, Qiskit provides a wide
range of implementations for popular quantum algorithms, which will assist the
future integration of additional quantum algorithms into the framework. Notice
that Qiskit has provided a complete VQA routine, but its recent update removed
its HHL module. To address this, we have adapted a Qiskit-based HHL package in
[8] and incorporated it into the framework.

Although Qiskit has provided quantum simulators in its Aer submodule, their
efficiency falls short when dealing with complex circuits like those in the
HHL algorithm. To overcome this limitation, our preferred simulator is the SV-
Sim simulator in NWQSim. It is a C++ CPU/GPU statevector simulator includes
various front-end supports such as C++, Qiskit, native Python, and Q#. When
compared to the statevector simulators in Qiskit, Cirq, and Q#, SV-Sim
exhibits exceptional performance in the benchmarking tests in [6]. The speed
comparisons between simulators in Qiskit Aer and NWQSim on power flow use
cases are also shown in Section IV.

###  II-C Embedded techniques in SV-Sim simulator

The majority of time consumption during HHL circuit simulation is attributed
to the process of simulation and measurement. Therefore, it’s crucial to
highlight several techniques employed by SV-Sim (NWQSim) to optimize
computation time for deep circuits. The most significant time sink is the
unitary evolution of the state. For SV-Sim, a noticeable improvement over the
simulators in Aer (from Qiskit) and qsim (from Cirq) is the specialized gate
implementation. While other simulators use general unitary gates for all one-
qubit and two-qubit gates, SV-Sim offers specialized computation for dozens of
supported basis gates and specific architectures of CPUs and GPUs. One
instance of such specialized computation is gate fusion. As illustrated in
Figure 4, gate fusion merges several applicable gates into a single gate based
on four distinct strategies. This substantially reduces both the gate counts
and the depth of the circuits. Taking the 210-gate circuit shown in Figure 1
as an example, the actual circuit executed in SV-Sim only has 67 gates. While
gate fusion is limited to certain basis gates, it does not necessitate the
matrix multiplications. Thus, the overall runtime can be lowered significantly
[6].

Another source of time inefficiency stems from latency arising from
communication between multiple computational units during parallel computing.
On one hand, SV-Sim employs the so-called “PGAS-based SHMEM” communication
model,

which provides enhanced support for software governing direct links between
GPUs compared to the traditional MPI routine. On the other hand, researchers
have also developed a novel method to minimize unnecessary data migration [6].
This is a preprint version of the work.
This work has been submitted to the IEEE for possible publication.
Copyright may be transferred without notice, after which this version may no
longer be accessible.

##  III Power Flow Use Case

Power flow is a fundamental power system application. In this paper, our
primary focus centers on the Newton-Raphson based AC power flow, a widely used
technique for power system analysis.

###  III-A Procedures for applying HHL on practical problems

Let us consider the linear system $\mathbf{A}\mathbf{x}=\mathbf{b}$ in an AC
power flow problem. $\mathbf{A}$ is assumed to be Hermitian and an easy
conversion can be conducted if it is not. HHL algorithm aims to find the
normalized solution vector $|\mathbf{x}\rangle$ by
$|\mathbf{x}\rangle=\mathbf{A}^{-1}|\mathbf{b}\rangle=\sum_{i}\lambda_{i}^{-1}b_{i}|\mathbf{v_{i}}\rangle$,
where $\lambda_{i}$ and $\mathbf{v_{i}}$ are the $i^{th}$ eigenvalue and
eigenvector of $\mathbf{A}$, respectively, and
$|\mathbf{b}\rangle=\sum_{i}b_{i}|\mathbf{v_{i}}\rangle$ is the normalized
vector $\mathbf{b}$ with decomposition. The decomposition of $\mathbf{b}$ in
the eigenbasis of $\mathbf{A}$ is achieved with quantum phase estimation
(QPE), as shown in the circuit in Fig. 5. The complete procedures used to
solve the linear system from the power flow problems are following.

Pratical HHL procedures Goal: Providing $\mathbf{A}$ and $\mathbf{b}$ in a
linear system $\mathbf{A}\mathbf{x}=\mathbf{b}$, use HHL to solve for solution
vector $\mathbf{x}$. Process: 1. First, normalize $\mathbf{b}$ to get
$|\mathbf{b}\rangle$ and record $\|\mathbf{b}\|$, then check the properties of
matrix $\mathbf{A}$: (a) if the condition number is too high, use a
preconditioner (b) if the dimension is not a power of 2, expand it to the
nearest power of 2 (c) if the matrix is not Hermitian (e.g., due to the
preconditioner), make it Hermitian:
$\displaystyle\left[\begin{array}[]{cc}\mathbf{0}&\mathbf{A}\\\
\mathbf{A}^{\dagger}&\mathbf{0}\end{array}\right]\left[\begin{array}[]{c}\mathbf{0}\\\
|\mathbf{x}\rangle\end{array}\right]=\left[\begin{array}[]{c}|\mathbf{b}\rangle\\\
\mathbf{0}\end{array}\right].$ 2. Generate the HHL circuit as in Figure 5. 3.
Compute $|\mathbf{x}\rangle$ by measuring the final state $|\vec{\psi}\rangle$
from the circuit. • If a simulator is used, the statevector of $|\psi\rangle$
can be read directly. • if a real quantum computer is used, the statevector of
$|\psi\rangle$ has to be reconstructed from state tomography. 4. Retrieve
$\mathbf{x}$ from $|\mathbf{x}\rangle$ and return it. (a) Compute
$\displaystyle\|\mathbf{x}\|=|\lambda_{0}|\cdot P(\text{measure ancilla and
get 1}),$ where $|\lambda_{0}|\leq|\lambda_{i}|$ for all $i\neq 0$, and
$\lambda_{i}$ are eigenvalues of the matrix input in the circuit. (b) Finally,
$\mathbf{x}=\|\mathbf{x}\|\cdot\|\mathbf{b}\|\cdot|\mathbf{x}\rangle$.

Figure 4: Four possible gate fusion strategies for applicable gates. The
individual gates on the left-hand side are combined into one fused gate on the
right-hand side. The priorities among four strategies is from top to the
bottom: first, one-qubit gates are fused; then two different fusion between
one-qubit and two-qubit gates happen; and finally, two-qubit gates are fused
together.  Figure 5: The conceptual circuit of HHL

###  III-B Resources estimation and scalability

Following the settings in [8] and after computing the first step in Section
III-A (expand and Hermitize the matrix $\mathbf{A}$), the estimated total
number of qubits used, $n_{\text{total}}$, is based on the following equation.

$\displaystyle n_{\text{total}}=n_{\text{data}}+n_{\text{QPE}}+n_{\text{neg
val}}$

where

$\displaystyle\centering\begin{cases}n_{\text{data}}&=\log_{2}(\text{number of
rows or columns of matrix})\\\
n_{\text{QPE}}&=\max(n_{\text{data}}+1,\lceil\log_{2}(\kappa+1)\rceil)\\\
n_{\text{neg val}}&=\begin{cases}1&\text{if the input matrix has negative
eigenvalues}\\\ 0&\text{else}\end{cases}\end{cases}\@add@centering$

and $\kappa$ is the condition number of matrix $\mathbf{A}$.

Note that $n_{\text{QPE}}$ governs the accuracy of the eigenvalue computation
and the solution vector $\mathbf{x}$, making it hard to reduce the qubit
requirement in the QPE step, while the computational power needed in quantum
simulators has an exponential dependency on $n_{\text{QPE}}$. To mitigate this
difficulty, a preconditioner is used for systems with relatively high
condition numbers in our routine. The current implementation utilizes the
Gauss-Seidel (GS) preconditioner due to its efficiency. In Section IV, a case
is provided to show the use of a GS preconditioner can improve the accuracy.

With NWQSim, a single NVIDIA A100 GPU can effectively handle around 32 qubits
and this capacity can be expanded to 42 qubits with the utilization of
multiple nodes. Considering memory space constraints, the optimal capability
of solving non-singular non-Hermitian linear systems ranges from
$2^{14}$-by-$2^{14}$ to $2^{20}$-by-$2^{20}$, depending on the number of

available GPU nodes. However, as shown in Table I, the current bottleneck is
in HHL circuit generation.  This is a preprint version of the work.
This work has been submitted to the IEEE for possible publication.
Copyright may be transferred without notice, after which this version may no
longer be accessible.

##  IV Numerical experiments

To demonstrate the capabilities of our framework, we selected the IEEE 14-bus
and IEEE 30-bus systems[12], which are readily available through the MATPOWER
package. These systems serve as our initial testbed, allowing us to showcase
the framework’s functionality.

Experiments were conducted on Perlmutter supercomputer using a single NVIDIA
A100 GPU. We will use multiple GPUs to accommodate the increased computational
demands for future large test cases. In our study, HHL circuits for both the
14-bus and 30-bus cases were executed on both Qiskit Aer and the SV-Sim
simulator. The detailed results are presented in Table I. A highlight is the
impact of the GS preconditioner on the 30-bus case. It can reduce the
condition number of the 30-bus system from 492.5 to 109.3. This reduction had
the effect of lowering the value of $n_{\text{qpe}}$ from 10 to 8, as
elaborated in the entry marked “30-buses∗” in the table.

Although GS preconditioner evidently improved the accuracy of the final
solution, it came at the cost of increased circuit generation time and thus
simulation time. That was because the preconditioner broke the Hermiticity of
the original matrix, resulting in deeper circuits. Additionally, SV-Sim
consistently outperformed Qiskit Aer, achieving simulation speeds up to 8
times faster. Notably, the running time of the C++ segment within SV-Sim was
1.7s, 36.2s, and 34.5s for the three cases, respectively. This observation
indicated that the actual backend computation time for each circuit ranged
from $1/150$ to $1/65$ of the time recorded in Python. This discrepancy can be
attributed to the substantial time overhead involved in casting data in
between C++ and Qiskit. Unfortunately, since it is difficult to obtain the
actual computation time from Qiskit Aer, the authors could not make a direct
comparison on simulation time without type-casting.

##  V Conclusion and future work

This paper presents a flexible hybrid framework, offering expedited problem-
solving through the incorporation of quantum algorithms while integrating with
existing power system tools. The framework is based on an advanced quantum
simulator that significantly minimize resource-intensive communications and
GPU-based parsing through the implementation of specialized computations. Use
case with power flow studies of IEEE test systems showed up to 8 times speedup
compared to Qiskit Aer with similar levels of accuracy. This framework has
potential to empower power system engineers, even those without a quantum
computing background, to explore and advance quantum computing and its
applications in power system domain.

This paper shows an initial exploration with a focus on HHL linear algorithm
in power system applications. As we look ahead, our vision includes several
key improvements to further enhance the framework’s capabilities. First, a C++
implementation of QPE and HHL circuit generation based on NWQSim would greatly
reduce the time overhead for circuit generation and type-casting data between
Qiskit and C++. With this new implementation in place, it will be feasible for
running large-scale linear problems with higher number of qubits and more GPUs
for scalability testing. Furthermore, authors will include more solvers into
the framework, such as various variational quantum solvers, and extend it to a
wide ranges of power system applications. Differential-algebraic equations
related power system dynamics modelings have received considerable attentions
as seen in [13]. Given the underline subroutine still requires linear solvers,
a worthwhile direction is to incorporated those into our framework.

TABLE I: Experiment results on 14-bus case, 30-bus case, and 30-bus case with GS preconditioner (30-bus∗) | 14-bus | 30-bus | 30-bus∗
---|---|---|---
Matrix Size | 13 $\times$ 13 | 29 $\times$ 29 | 29 $\times$ 29
Condition Number | 119.2 | 492.5 | 109.3
$n_{\text{total}}$ and ($n_{\text{data}}$, $n_{\text{QPE}}$) | 13 and (4,8) | 16 and (5,10) | 15 and (6,8)
Circuit Generation Time | 88s | 30.5 min | 37.5 min
Qiskit Aer Simu. Time | 27s | 31.3 min | 37.5 min
\- Error $\|\mathbf{x}_{hhl}-\mathbf{x}_{true}\|_{2}$ | $1.97\cdot 10^{-3}$ | $1.18\cdot 10^{-3}$ | $6.35\cdot 10^{-4}$
SV-Sim Simu. Time | 13s | 3.9min | 4.7 min
\- Error $\|\mathbf{x}_{hhl}-\mathbf{x}_{true}\|_{2}$ | $1.97\cdot 10^{-3}$ | $1.18\cdot 10^{-3}$ | $6.35\cdot 10^{-4}$

## References

  * [1] Edison International. Countdown to 2045: Realizing California’s Pathway to Net Zero September 2023. IEEE Electrification Magazine.
  * [2] Zhou, Y., Tang, Z., Nikmehr, N., et al. Quantum computing in power systems. iEnergy, 2022, 1(2): 170-187.
  * [3] Golestan, S., Habibi, M. R., Mousazadeh Mousavi, S.Y., Guerrero, J. M., and Vasquez, J. C. Quantum computation in power systems: An overview of recent advances. Energy Reports 9 (2023): 584-596.
  * [4] Chen, Y., Stein, S., Li, A., and Huang, Z. H. (2023, July). Is It Coming Soon to Power Systems: Quantum Computing and Its Early Exploration. In 2023 IEEE Power & Energy Society General Meeting (PESGM) (pp. 1-5). IEEE.
  * [5] Li, A., Subasi, O., Yang, X. & Krishnamoorthy, S. Density Matrix Quantum Circuit Simulation via the BSP Machine on Modern GPU Clusters. Proceedings Of The International Conference For High Performance Computing, Networking, Storage And Analysis. (2020). Available at https://github.com/pnnl/NWQ-Sim
  * [6] Li, A. et al. SV-Sim: Scalable PGAS-Based State Vector Simulation of Quantum Circuits. Proceedings Of The International Conference For High Performance Computing, Networking, Storage And Analysis. (2021). Available at https://github.com/pnnl/NWQ-Sim
  * [7] Harrow, A., Hassidim, A. & Lloyd, S. Quantum Algorithm for Linear Systems of Equations. Phys. Rev. Lett.. 103, 150502 (2009,10).
  * [8] Vazquez, A. quantum_linear_solvers. GitHub Repository. (2022), https://github.com/anedumla/quantum_linear_solvers
  * [9] Yalovetzky, R., Minssen, P., Herman, D., and Pistoia, M. Hybrid HHL with Dynamic Quantum Circuits on Real Hardware. ArXiv Preprint ArXiv:2110.15958. (2023).
  * [10] Zimmerman, R., Murillo-Sánchez, C., and Thomas, R. MATPOWER: Steady-State Operations, Planning, and Analysis Tools for Power Systems Research and Education. IEEE Transactions On Power Systems. 26, 12-19 (2011)
  * [11] Qiskit contributors, Qiskit: An Open-source Framework for Quantum Computing. (2023)
  * [12] Power Systems Test Case Archive. Available at: https://labs.ece.uw.edu/pstca/
  * [13] Tran, H. T. T., Nguyen, H. T. , Vu, L. T., and Ojetola, S. T., Applying Quantum Computing to Simulate Power System Dynamics’ Differential-Algebraic Equations, 2023 North American Power Symposium (NAPS), Asheville, NC, USA, 2023, pp. 1-6.