1,21,2affiliationtext: Department of Aerospace Engineering, University of
Maryland College Park, MD, 2074211affiliationtext: Department of Computer
Science, University of Maryland College Park, MD, 20742

# High-Entanglement Capabilities for Variational Quantum Algorithms: The
Poisson Equation Case

Fouad Ayoub James D. Baeder

The discretized Poisson equation matrix (DPEM) is vital to the field of
computational fluid dynamics (CFD), and an algorithm that solves it on a
quantum computer could potentially grant exponential space and time complexity
speedups. However, the DPEM in 1D has been shown to have trouble being
represented as a decomposition in the Pauli basis. Additionally, traditional
ansatz for Variational Quantum Algorithms (VQAs) that are used to
heuristically solve linear systems (such as the DPEM) have large numbers of
parameters, making them harder to train. This research attempts to resolve
these problems by utilizing the IonQ Aria quantum computer capabilities that
boast all-to-all connectivity of qubits. We propose a decomposition of the
DPEM that is based on 2- or 3-qubit entanglement gates and is shown to have
$O(1)$ terms with respect to system size, with one term having an $O(n^{2})$
circuit depth and the rest having only an $O(1)$ circuit depth (where $n$ is
the number of qubits defining the system size). To test these new
improvements, we ran numerical simulations to examine how well the VQAs
performed with varying system sizes, showing that the new setup offers an
$O(n)$ scaling of the number of iterations required for convergence, providing
an exponential speedup over their classical computing counterparts. This
project shows that the future of computational fluid dynamics may involve
quantum computers to provide significant time and space complexity speedups.

##  1 Introduction

The Poisson equation is a partial differential equation that appears in many
areas of physics and engineering, but one important application of the Poisson
equation is its usage in fluid dynamics. This equation appears as a relation
between the pressure and velocity field of an incompressible flow. As such,
this equation holds an essential truth to any sort of analysis involving
incompressible and viscous flows, for example in supernovae explosions [1] or
blood flow in arteries [2]. One of the more common classical methods of
solving the Poisson is by discretizing the surface into a mesh grid and
solving for the conditions at these points, reducing the problem to solving a
system of linear equations [3, 4, 5]. There are also some methods of solving
the Poisson in recent years that include the use of modern machine-learning
techniques [6, 7].

The focus of this paper is to develop methods of solving the one-dimensional
discretized Poisson equation matrix (DPEM) on a quantum computer. This
equation is one of many central pieces to the solutions of the Naiver-Stokes
equation, and is key to solving a multitude of fluid-flow problems. As such,
the push for greater and greater computational capabilities for simulating
such fluid flows using the Poisson equation requires new techniques to be
developed and improved upon. Our overall task is to solve some system of
equations $A\mathbf{x}=\mathbf{b}$, where $A$ is the DPEM. There exists an
algorithm known as the HHL algorithm [8] that can explicitly solve for
$\mathbf{x}$, even providing an exponential speedup; however, that method
requires many noiseless qubits to run on reasonably sized systems. In recent
literature, it has been shown that the HHL algorithm can solve the Poisson,
but only for small systems [9, 10]. In the current era of quantum computing,
we only have access to a limited number of qubits, each with some noise
guarantees. This era is referred to as the Noisy Intermediate Quantum (NISQ)
era of computing [11]. Recently, Variational Quantum Algorithms (VQAs) have
had favor in the literature for near-term use, as they provide some resilience
to quantum error as well as some possibilities of speedups [12]. These VQAs
use the quantum computer as the main method of computing a cost function,
while a classical computer performs optimization across the parameter space of
the quantum circuit. The Variational Quantum Linear Solver (VQLS), as proposed
in [13], provides a framework for which any linear system of equations can be
solved.

Typically, the VQLS is solved by using some version of a Local-Entangling
Ansatz (LEA), also known as a Hardware-Efficient Ansatz (HEA), due to the
nature of architecture of most currently available quantum computers, such as
the IBM quantum computers [14]. These ansatz take advantage of the fact that
most quantum computers can much more easily entangle qubits with their
immediate neighbors, without requiring any decomposition into native gates.
However, the usage of ansatz involving high-entanglement schemes, what this
paper introduces as Globally-Entangling Ansatz (GEA) has yet to be put into
practice, and with the advent of IonQ quantum computers that boast “all-to-all
interconnectivity” of qubits [15], it may be wise to explore these avenues. We
will show that these ansatz require much less rotation and entanglement layers
than HEA, allowing for optimization algorithms to converge faster and ensuring
circuit depth is reduced. The resultant scheme offers an approximately linear
scaling of the number of iterations required for convergence. Since $n$ is the
number of qubits defining the system size, which scales exponentially with
respect to $n$, this is an exponential speedup compared to the classical
algorithm for solving linear systems.

Additionally, in typical analysis of solving the DPEM with VQLS, the matrix
itself is decomposed into a linear combination of sub-matrices, which has a
number of terms that scales either exponentially [16] or polynomially [17]
with respect to system size. This means that as the size of the matrix gets
larger, the number of quantum circuits required to evaluate the cost function
grows exponentially or polynomially. In this paper, we propose the usage of a
high-entanglement decomposition (HED), which provides a constant scaling of
sub-matrices with respect to system size, containing only four terms in the
final decomposition. The decomposition is inspired by the HED proposed in
[16]. However, their decomposition is slightly incorrect and not optimal.
Additionally, three of the terms have an $O(1)$ circuit depth, and one term
has an $O(n^{2})$ circuit depth, meaning that the HED does not require the
evaluation of quantum circuits with circuit depth longer than the maximum on
IonQ quantum computers [18].

There are some papers that claim that VQA schemes for solving system of
equations cannot be done effectively in the NISQ era of quantum computers
[19]. While their papers may hold weight for algorithms that use traditional
ansatz and local entangling circuits, this paper shows that taking advantage
of new types of quantum computer architectures may be the key to allowing VQAs
to realize their full potential. As shown in our numerical simulations,
assuming one has access to high-entanglement quantum computers allows for
faster convergence to the solution and less jobs per cost evaluation. VQLS is
only one of the many potential algorithms that may be benefited from the fully
connected qubit assumption. This paper shows that imagining algorithms that
use high-entanglement schemes may be the path forward for improving the
current state of quantum computing, bringing us closer to an age of practical,
commercial use of quantum computers.

Encode matrix $A=\sum_{l}c_{l}A_{l}$ and $b=\sum_{l}U_{l}$Prepare
$\ket{\psi(\boldsymbol{\theta})}$ in quantum variational ansatzCalculate cost
function value $C$ with Hadamard Test$C<\gamma$?Find new $\boldsymbol{\theta}$
with optimizer$\ket{\psi(\boldsymbol{\theta})}\approx\ket{x}$noyes Figure 1:
Depiction of the general workflow of the algorithm. This algorithm is based on
the one proposed in [13].

##  2 Variational Quantum Linear Solver

The following section will be rephrasing and summarizing the theory behind the
VQLS [13], and will specifically describe how the VQLS structure is being
implemented in this paper.

The VQLS’ goal is to solve some arbitrary linear system of equations, namely,
finding $\mathbf{x}$ given some unitary matrix $A$ and some vector
$\mathbf{b}$ that solves the equation $A\mathbf{x}=\mathbf{b}$. The goal of
any algorithm that solves a Quantum Linear Systems Problem (QLSP) is to find
some quantum circuit that prepares a state $\ket{x}$ such that
$\ket{x}\propto\mathbf{x}$.

###  2.1 Algorithm Overview

The VQLS uses a heuristic approach to solving the linear system
$A\mathbf{x}=\mathbf{b}$. The schematic of the algorithm can be seen in Figure
1. It begins by assuming that the matrix $A$ can be decomposed in some finite
number of sub-matrices that can be represented as a combination and tensor
product of quantum logic gates. Additionally, we assume that the vector
$\mathbf{b}$ can be decomposed into a short-depth quantum circuit denoted as
$U$, comprised of some sub-unitaries $U_{l}$, such that a quantum state
$\ket{b}\propto\mathbf{b}$ can be prepared in constant circuit depth.

$\displaystyle A=\sum_{l}c_{l}A_{l}$ (1) $\displaystyle b=\sum_{l}U_{l}$ (2)

We initially prepare some guess of the solution as a quantum state, which we
call $\ket{\psi(\boldsymbol{\theta})}$, parameterized by a vector of rotations
$\boldsymbol{\theta}$. This guess can be prepared using a parameterized
quantum circuit known as a variational ansatz, in which the parameters
$\boldsymbol{\theta}$ will determine certain rotation values of the circuit.
The ansatz we use in this paper is a combination of parameterized rotation
layers (using $R_{y}$ gates) and entanglement layers (using $CZ$ gates).

Once this guess is prepared, we evaluate our cost function (described below),
and use the result to perform some type of minimization across the parameters
of $\boldsymbol{\theta}$. Once the cost function has reached some convergence
threshold, the optimization loop will exit, and the probability vector of the
resulting guess $\ket{\psi(\boldsymbol{\theta})}$ may be read out as a
proportional solution to the system.

###  2.2 The Poisson Equation

The DPEM in 1D is as follows. It will be the matrix of interest in all of our
numerical simulations. Its derivation is not the importance of this paper, but
may be found in [16].

$\displaystyle A=\begin{bmatrix}2&-1&0&\cdots&0\\\ -1&2&-1&\cdots&\vdots\\\
0&\ddots&\ddots&\ddots&0\\\ \vdots&\cdots&-1&2&-1\\\
0&\cdots&0&-1&2\end{bmatrix}$ (3)

The most difficult part about the VQLS algorithm is the decomposition of the
desired matrix into sub-matrices that can be represented on a quantum
computer. For the DPEM, there are 2 decompositions we will discuss in this
paper. The first utilizes the Pauli basis of quantum logic gates, which is
easy to understand and has short circuit lengths when applied to the actual
algorithm. The second decomposition is more complicated and requires quantum
circuits that entangle almost all qubits with each other. However, this
decomposition will be favorable in later use because its decomposition has an
$O(1)$ scaling of sub-matrices and a generally well-behaved scaling of circuit
depth, which is $O(n^{2})$. With currently available hardware (based on
current circuit depth limitations of the IonQ Aria machine [20]), one could
already simulate a one-million particle mesh with only twenty qubits using the
HED method, as shown in Figure LABEL:fig:HED.

####  2.2.1 Decomposition in Pauli Basis

The first, and simplest decomposition of the Poisson matrix is a linear
combination of multiple Pauli gates. The derivation of a recursive formula
based on system size can be found in [16], but here we list the decompositions
for 4x4 and 8x8 DPEM.

N=4:

$\displaystyle\begin{split}A=2I\otimes I-I\otimes X\\\ -0.5X\otimes
X-0.5Y\otimes Y\end{split}$ (4)

N=8:

$\displaystyle\begin{split}A=2I\otimes I\otimes I-I\otimes I\otimes X\\\
-0.5I\otimes X\otimes X-0.25X\otimes X\otimes X\\\ -0.25Y\otimes Y\otimes
X-0.25Y\otimes X\otimes Y\\\ -0.5I\otimes Y\otimes Y+0.25X\otimes Y\otimes
Y\end{split}$ (5)

The benefit to using the Pauli basis for decomposition is that we are left
with relatively simple sub-matrices to work with when calculating expectation
values of the quantum state. Also, the circuit depth is guaranteed to be
$O(1)$ for each quantum circuit required to calculate the cost function.
However, the Pauli basis decomposition scales exponentially with system size
[16], meaning that solving the Poisson matrix in this manner would not offer
an asymptotic speedup when compared to classical methods. For this reason, the
Pauli decomposition was not favored in our numerical simulations and runs on
the IonQ quantum computers.

####  2.2.2 High Entanglement Decomposition

A second decomposition of the Poisson equation exists, although it requires
more complicated procedures for its generation, and it requires that the
quantum computer for which it runs on can complete two or three qubit
entanglement gates between almost all qubits. IonQ’s latest quantum computers
indeed are fully capable of such operations [20], so their quantum computers
are a prime candidate for testing this decomposition of the matrix on a real
quantum computer. The HED shortens the number of expectation values required
to be evaluated per iteration, has a constant number of sub-matrices, and has
reasonable circuit depth for each quantum circuit required to evaluate the
cost function, as will be shown in the following derivation and can be seen in
Figure LABEL:fig:HED. This decomposition takes much less computational
resources overall compared to the Pauli decomposition, so we use it in all our
numerical simulations.

The HED proposed in this paper is composed of only four sub-matrices. The
linear combination of the matrices are as follows:

$\displaystyle A=2.5I-L_{1}-L_{2}-0.5L_{3}$ (6)

The matrix $I$ is the identity matrix acting on all qubits. The remaining
matrices are defined below.

$L_{1}$ is simply the $X$ gate applied to the least significant qubit.

$\displaystyle L_{1}=I^{\otimes n}\otimes X$ (7)

This gives us the matrix

$\displaystyle L_{1}=\begin{bmatrix}0&1&&&\\\ 1&0&&&\\\ &&\ddots&&\\\
&&&0&1\\\ &&&1&0\end{bmatrix}$ (8)

which clearly has an $O(1)$ circuit depth. A.1(a) shows the quantum circuit
for $L_{1}$ for the $n=6$ case.

$L_{2}$ is more complicated to craft, and it requires the definition of an
extra sub-matrix, $C_{i}$, which is composed of $CX$s and multi-$CX$s. The
definition used the following notation of $CX$ gates: $CX^{b}_{a}$ describes
an $X$ gate acting on $a$ controlled by $b$. $mCX$ gates simply can have more
controls $b_{0},b_{1},\ldots$

$\displaystyle\begin{split}C_{i}&=CX^{i}_{i-1}CX^{i}_{i-2}\cdots CX^{i}_{0}\\\
&mCX^{0,1,\ldots,i-1}_{i}\\\ &CX^{i}_{0}\cdots
CX^{i}_{i-2}CX^{i}_{i-1}\end{split}$ (9)

A.LABEL:fig:Cgates depicts the quantum circuits representing
$C_{1},C_{2},\ldots,C_{5}$ as an example for the $n=6$ case.

$L_{2}$ is thus formed in the following manner

$\displaystyle L_{2}=C_{1}C_{2}\ldots C_{n-1}$ (10)

where $n$ is the number of qubits. Figure LABEL:fig:Cgates(f) shows the final
result for $n=6$. This gives us the matrix

$\displaystyle L_{2}=\begin{bmatrix}1&&&&&\\\ &0&1&&&\\\ &1&0&&&\\\
&&&\ddots&&\\\ &&&&0&1\\\ &&&&1&0\\\ &&&&&&1\end{bmatrix}$ (11)

The $L_{2}$ circuit depth is

$\displaystyle 3+5+7+\ldots+2(n-1)+1$ (12)
$\displaystyle=\sum_{i=1}^{n-1}2i+1$ (13) $\displaystyle=O(n^{2})$ (14)

Finally, we just need the $L_{3}$ matrix. This one is much more easily
implemented. It is simply applying multi-controlled $Z$ gates and an $X$ gate
on the least-significant qubit in alternating order.

$\displaystyle L_{3}=mCZ^{0,1,\ldots n-1}\>X\>mCZ^{0,1,\ldots n-1}\>X$ (15)

The resultant $L_{3}$ matrix is

$\displaystyle L_{3}=\begin{bmatrix}-1&&&&\\\ &1&&&\\\ &&\ddots&&\\\ &&&1&\\\
&&&&-1\end{bmatrix}$ (16)

An example of the implementation of this matrix as a quantum circuit is shown
for the $n=6$ case in Appendix A.1.

###  2.3 The Cost Function

We define a cost function, which will determine how good the current guess
$\ket{\psi}$ is by finding out how close $A\ket{\psi}$ is to $\ket{b}$. This
cost function is inherently hybrid, in that it will be partly evaluated with a
series of quantum circuits that evaluate circuit matrix arithmetic, as well as
some classical matrix arithmetic.

####  2.3.1 Global Entangling Ansatz

In order to generate our $\ket{\psi(\boldsymbol{\theta})}$ state, we use a
Global-Entangling Ansatz (GEA), which utilizes the capabilities of the IonQ
hardware. The quantum computers available through IonQ’s quantum computing
cloud services boast all-to-all connectivity between qubits, allowing one to
easily encode two-qubit gates between any two qubits [15]. We define the GEA
with multiple layers of parameterized rotation and entanglement gates, where
the parameters are defined by $\boldsymbol{\theta}$. Additionally, we pre-
conditioned the ansatz with the unitary that prepares the b state, denoted by
$U_{b}$. We decided to see whether or not this type of entangling scheme would
benefit the performance of the cost functions, and as shown through our
numerical results, we found that it offered great improvements with respect to
number of qubits and number of iterations when compared to the LEA (or HEA)
proposed in [13]. There are certain factors going into how many layers (and
thus the number of parameters) the ansatz should have, but simply put, if the
parameter space is too small, the solution might not exist in it at all, but
if it is too large, it might take far too long to find the solution [12]. As
can be seen in Figure 3, the GEA has far fewer parameters than the HEA, which
will make it easier to train. The high-entanglement nature of the GEA
preserves its expressibility, thus making the solution to each system still
well-accessible to the ansatz. The GEA used in this paper is defined with
three repeated layers of rotation gates ($R_{y}$) and entangling gates ($CZ$),
with a unique rotation parameter attributed to each $R_{y}$ gate. Thus, the
number of parameters per variational ansatz is $3n$ where $n$ is the number of
qubits defining the system size. An example of the GEA for $n=6$ can be found
in Figure 4.

Figure 3: Parameters versus qubits depending on the ansatz chosen.

Figure 4: The global-entangling variational ansatz, defined with $3n$
parameterized gates and $CZ$ entangling gates, which prepares the
$\ket{\psi(\boldsymbol{\theta})}$ state. This figure depicts the $n=6$ case.
Here, we pre-condition the ansatz with the gate sequence that prepares the b
state, $U_{b}$; in this case, that gate sequence is $H^{\otimes n}\ket{0}$

####  2.3.2 The Hadamard Test

The Hadamard Test is a tool used to actually determine the expectation values
of evaluating vector products of unitary matrices, or compositions of unitary
matrices, with prepared states [21]. The general idea is that, given some
initially prepared state, $\ket{\psi}$, and some unitary matrix, $U$, the
value of $\bra{\psi}U\ket{\psi}$ can be determined by applying the Hadamard
gate to an auxiliary qubit, then applying a $CU$ gate (a controlled version of
that unitary gate) with the auxiliary qubit as the control, and finally
applying another Hadamard gate on the auxiliary qubit. The resulting
probabilities of measuring the auxiliary qubit as a 0 or 1 can reconstruct the
expectation value as follows:

$\displaystyle P(0)-P(1)=Re\bra{\psi}U\ket{\psi}$ (17)

This uses a phenomenon known as “phase kickback” and is further explained in
other works [13]. Since we do not use or expect any gates, matrices, or
vectors to contain imaginary numbers, we neglect the calculation of the
imaginary part of the expectation values.

####  2.3.3 Locally and Globally Defined Cost

The cost function must signify how close we are to the real solution, where
the larger it is, the further away we are, and the smaller it is, the closer
we are. Here we inherit the “projection” Hamiltonian from [13], defined as

$\displaystyle H_{P}=\mathbb{I}-\ket{b}\bra{b},$ (18)

so our cost function $C_{G}$ may be defined as the expectation value of the
projection Hamiltonian, given that $\Phi=A\ket{\psi(\boldsymbol{\theta})}$:

$\displaystyle C_{G}$ $\displaystyle=\bra{\Phi}H_{P}\ket{\Phi}$ (19)
$\displaystyle=\bra{\Phi}(\mathbb{I}-\ket{b}\bra{b})\ket{\Phi}$ (20)
$\displaystyle=\braket{\Phi}{\Phi}-\braket{\Phi}{b}\braket{b}$

So when $\ket{\Phi}$ and $\ket{b}$ are orthogonal, the inner product
$\braket{b}{\Phi}$ gets smaller, making the cost greater, and when the vectors
get closer to each other, that inner product gets larger, making the cost
smaller. However, we can normalize the cost to give it a better operational
meaning [13]. Taking the norm with respect to $\Phi$, we get the following new
cost function:

$\displaystyle\hat{C}_{G}$
$\displaystyle=\frac{\braket{\Phi}{\Phi}}{\braket{\Phi}{\Phi}}-\frac{\braket{\Phi}{b}\braket{b}{\Phi}}{\braket{\Phi}{\Phi}}$
(21)
$\displaystyle=1-\frac{\braket{\Phi}{b}\braket{b}{\Phi}}{\braket{\Phi}{\Phi}}$
(22) $\displaystyle=1-\frac{|\braket{b}{\Phi}|^{2}}{\braket{\Phi}{\Phi}}$ (23)

To break this down into quantum gate components (where $U$ is some unitary
preparing the state $\ket{b}$, and $V$ is some variational unitary preparing
the state $\ket{\psi(\boldsymbol{\theta})}$):

$\displaystyle\bra{b}$ $\displaystyle=(U\ket{0})^{\dagger}=\bra{0}U^{\dagger}$
(24) $\displaystyle\ket{\Phi}$
$\displaystyle=A\ket{\psi(\boldsymbol{\theta})}=AV(\boldsymbol{\theta})\ket{0}$
(25)

Since we know that the $A$ matrix is a linear combination of some decomposed
sub-matrices, we generalize $A$ to being

$\displaystyle A=\sum_{l}c_{l}A_{l},$ (26)

where each $A_{l}$ is a sub-matrix that can be represented as a quantum gate.
This gives us

$\displaystyle\ket{\Phi}$
$\displaystyle=\sum_{l}c_{l}A_{l}\ket{\psi(\boldsymbol{\theta})}=\sum_{l}c_{l}A_{l}V(\boldsymbol{\theta})\ket{0}$
(27)

Now we can rewrite $\hat{C}_{G}$ as follows:

$\displaystyle\hat{C}_{G}=$ (28) $\displaystyle
1-\frac{\sum_{l,l^{\prime}}c_{l}c_{l^{\prime}}\bra{0}V(\boldsymbol{\theta})^{\dagger}A_{l}^{\dagger}U\ket{0}\bra{0}U^{\dagger}A_{l^{\prime}}V(\boldsymbol{\theta})\ket{0}}{\sum_{l,l^{\prime}}c_{l}c_{l^{\prime}}\bra{0}V(\boldsymbol{\theta})^{\dagger}A_{l}^{\dagger}A_{l^{\prime}}V(\boldsymbol{\theta})\ket{0}}$

This is in fact our “global” cost function. A “local” cost function is defined
as such:

$\displaystyle\hat{C}_{L}=$ (29)
$\displaystyle\frac{1}{2}-\frac{1}{2n}\frac{\sum_{l,l^{\prime},j}c_{l}c_{l^{\prime}}\bra{0}V(\boldsymbol{\theta})^{\dagger}A^{\dagger}_{l}UZ_{j}U^{\dagger}A_{l^{\prime}}V(\boldsymbol{\theta})\ket{0}}{\sum_{l,l^{\prime}}c_{l}c_{l^{\prime}}\bra{0}V(\boldsymbol{\theta})^{\dagger}A_{l}^{\dagger}A_{l^{\prime}}V(\boldsymbol{\theta})\ket{0}}$

where the $\ket{0}\bra{0}$ from $\hat{C}_{G}$ has been replaced by the
positive operator

$\displaystyle P=\frac{1}{2}+\frac{1}{2n}\sum_{j=0}^{n-1}Z_{j}.$ (30)

The local cost function is shown to avoid barren plateau problems that appear
in the global cost function [13]. Both of these cost functions were tested in
this paper, but due to the severe barren plateau problems that appear for the
global cost function, all results shown are from utilizing the local cost
function.

The distance between the quantum state $\ket{\psi(\boldsymbol{\theta})}$ and
the real solution (as represented by a proportional quantum state, $\ket{x}$),
is denoted by $\epsilon$, which is defined by

$\displaystyle\epsilon=\frac{1}{2}Tr|\ket{x}\bra{x}-\ket{\psi}\bra{\psi}|.$
(31)

We have inherited this definition as well as the operational relation between
the local and global cost functions that is introduced in [13], which includes
the following bounds:

$\displaystyle\hat{C}_{G}\geq\frac{\epsilon^{2}}{\kappa^{2}},C_{G}\geq\frac{\epsilon^{2}}{\kappa^{2}},\hat{C}_{L}\geq\frac{\epsilon^{2}}{n\kappa^{2}},C_{L}\geq\frac{\epsilon^{2}}{n\kappa^{2}},$
(32)

where $\kappa$ is the condition number of the matrix, and $n$ is the number of
qubits defining the size of the matrix. We use these bounds in our numerical
simulations, where the optimization routine attempts to reduce the cost
function below these thresholds in order to reach a desired $\epsilon$ error
(set by the user), given the $n$ and $\kappa$ attributed to the matrix.

For the sake of simplicity, we took $\ket{b}=H^{\otimes n}\ket{0}$.

####  2.3.4 Quantum Circuit Evaluations per Local Cost Function Evaluation

Due to the nature of the cost function definitions including multiple
summations over expectation values, there is a need to quantify exactly how
many quantum circuit evaluations are required to complete a single cost
function evaluation. In most literature, number of iterations are usually
presented in the main results (for example, in [13, 16], etc). However, the
number of quantum circuits required to reach a solution is a much more
practically important convergence metric. To illustrate this point, one could
imagine a cost function that has a blowup in the number of quantum circuits
required to evaluate the cost function with respect to system size. A
contrived example of this is if one were to use the Pauli decomposition of the
DPEM, which would require the number of circuit evaluations for a cost
function evaluation to scale exponentially with respect to system size. Most
of the time, the dependence on system size is linear at worst, but the example
still demonstrates the importance of communicating exactly how one’s cost
function computational scales.

As such, we now will show the number of quantum circuits required to run the
local cost function, since the global cost function is not the focus of this
paper. In order to reduce the number of quantum circuit required to evaluate
the local cost function, one can use the following equivalences:

$\displaystyle
l=l^{\prime}\Longrightarrow\bra{0}V(\boldsymbol{\theta})^{\dagger}A_{l}^{\dagger}A_{l^{\prime}}V(\boldsymbol{\theta})\ket{0}=1$
(33) $\displaystyle\bra{0}V(\boldsymbol{\theta})^{\dagger}$ $\displaystyle
A_{l}^{\dagger}A_{l^{\prime}}V(\boldsymbol{\theta})\ket{0}=$ (34)
$\displaystyle\bra{0}V(\boldsymbol{\theta})^{\dagger}A_{l^{\prime}}^{\dagger}A_{l}V(\boldsymbol{\theta})\ket{0}$
$\displaystyle\bra{0}V(\boldsymbol{\theta})^{\dagger}$ $\displaystyle
A^{\dagger}_{l}UZ_{j}U^{\dagger}A_{l^{\prime}}V(\boldsymbol{\theta})\ket{0}=$
(35)
$\displaystyle\bra{0}V(\boldsymbol{\theta})^{\dagger}A^{\dagger}_{l^{\prime}}UZ_{j}U^{\dagger}A_{l}V(\boldsymbol{\theta})\ket{0}$

After doing these reductions, one finds that the number of quantum circuits
required to evaluate the local cost function is

$\displaystyle N_{q}$ $\displaystyle=n\sum_{i=1}^{c}i+\sum_{i=1}^{c-1}i$ (36)
$\displaystyle=\frac{1}{2}c[n(c+1)+c-1],$ (37)

where $c$ is the number of sub-matrices in the decomposition of $A$. For the
HED, which all of our numerical results use (to avoid the exponential blowup
of the Pauli decomposition), $c=4$, so we see that $N_{q}$ has a mild linear
dependence on system size.

$\displaystyle N_{q}$ $\displaystyle=10n+6$ (38) $\displaystyle=O(n)$ (39)

One might now understand how $c=O(2^{n})$ would cause an exponential blowup of
the time required to evaluate a single cost function. While we know $n$ is
still exponentially related to system size, this would result in the
cancellation of any exponential speedups we might hope to achieve with quantum
computing. We will later show how the number of quantum circuits required per
cost function evaluation affects the time-to-solution on real quantum
computers.

Brief Note on Parallelization: One might notice that these quantum circuits
may be parallelized to improve predicted time-to-solutions metrics that we
will show below (following a method similar to [22]). However, according to
IonQ Application Scientist Dr. Daiwei Zhu [23], doing so would introduce more
noise into the circuits, making it much more difficult to train the ansatz. It
is best to not operate under the parallel assumption at all in this case, so
in this paper, we make no mention of a speedup introduced with
parallelization.

##  3 Numerical Results

In the following numerical results, we utilized the Pennylane library for
quantum circuit representation. Since we found that the Pauli decomposition
required an exponential scaling of the sub-matrices, we used the HED of the
Poisson in all cases. System sizes were only considered from 3 to 9 qubit
systems. We noticed that, although both the global and local cost functions
could successfully minimize and output the correct solution vector, for higher
qubit systems, the global cost function often encounters barren plateau
problems. For this reason, all our numerical simulations used the local cost
function in order to give both ansatz the chance to converge to a solution.
Also, we tested two different variances in the randomized weight
initialization, $q_{\Delta}=0.01$ and $q_{\Delta}=0.1$. We found the most
success with $q_{\Delta}=0.01$, so the conclusions about the results in this
section are drawn based on those simulations. However, all remaining
simulations done with $q_{\Delta}=0.1$ are shown in Appendix B.

###  3.1 Practical Constraints

The goal of this paper is to provide truly near-term quantum algorithms that
could be used to solve complex systems of equations, such as the DPEM, within
the next 5-10 years. As such, we believe that the practical constraints (that
will be explained below) we placed on our numerical simulations were
justified.

Since we used Pennylane’s fastest numerical quantum simulator, the Lightning
Qubit, our tests were done at a much faster pace than any state-of-the-art
quantum computer (which would be far slower due to the sampling and validation
requirements completed on the back-end). For this reason, we imposed a
practical limitation on our numerical testing, namely that we would not use
super-computers or high-performance GPUs to further speed up the computational
process. All the simulations were run on a Google Colab notebook, which had a
time-out period of 12 hours. If a run of the simulator did not complete
sufficient iterations before the 12 hour mark (i.e., getting past more than
one or two iterations), we removed those results from the data. In reality,
simulations that we see take near 8-10 hours on the Pennylane Lightning Qubit
simulator would take upwards of days on a real quantum computer. At this order
of magnitude, allowing the use of super-computers or high-performance GPUs to
speed-up the simulations would be pointless; even if these results worked, the
practical implementation of such enormous numbers of quantum circuit
evaluations would take unreasonably long to run on quantum computers
scientists have access to within the next 5-10 years.

To better understand this problem, we have shown the time-to-solution metrics
as both number of iterations, cost function evaluations, and estimated days
required to converge to a solution with respect to system size. As per IonQ’s
website [24] and the advice of IonQ Application Scientist Dr. Daiwei Zhu, each
quantum circuit takes $1\pm 0.5$ minutes to run. Using this metric, we
calculated the estimated number of days it would require to complete one run
of the VQLS for solving the DPEM. We also compared these numbers to the less
daunting iteration counts per run with respect to system size in Figure
LABEL:fig:0_01_metrics. It is worth noting that these estimated times assume
that the result of each quantum circuit evaluation offers a sufficient
precision that is required to minimize below the convergence threshold. More
on this is examined in Section 4.

###  3.2 Comparison of the HEA and GEA

We tested the difference in convergence metrics between using the HEA and the
GEA with numerical, state-vector simulations. In Figure
LABEL:fig:0_01_metrics(a), we see the iterations scaling is linear for both
the HEA and GEA, but the GEA has a far better scaling compared to the
predicted scaling of the HEA. If we only examined the iteration scaling,
however, we would miss out on a more interesting metrics. For example, in
Figure LABEL:fig:0_01_metrics(b), we can see that what we initially thought
was a benign linear scaling with respect to qubits turns out to have
significant consequences for the number of quantum circuits required to
converge to a solution. For the HEA, the simulated and predicting number is in
the millions, with a large-coefficient linear scaling, while the GEA remains
below 0.5e6 quantum circuits and has a much better coefficient for its linear
scaling. Finally, when we examine the estimated days of runtime required for
solving each system, we can see again the significant consequences of a bad
linear coefficient. While the GEA is constrained to no more than 5 days for
system sizes $n=3$ to $n=9$, the predicted scaling of the HEA exposes how
unreasonably long the time-to-solutions would be if one were to run them on
real quantum computers. Ideally, our algorithms would be best if the number of
quantum circuits required to run was kept at a minimum, since there are many
time-consuming components of running a quantum circuit on a real quantum
computer that cannot be represented in quantum simulators like Pennylane.

###  3.3 Scaling of $\kappa$

In the special case of the DPEM, the condition number scales exponentially
with the number of qubits, as can be seen in Figure 7 (note that this simply
means that its condition number scales linearly with respect to system size).
As shown in [25], linear scaling of iterations with respect to condition
number is optimal in the general case. Since condition number is exponentially
related to qubits, and the iterations scale linearly with with respect to
qubits for both the HEA and GEA, both methods are sub-linear, indicating at
least optimally in both cases. This sub-linear scaling is found also in [13],
and is thought to be due to the special case of the problem being examined.

Figure 7: Condition number versus qubits defining solution vector. Note that
the y-axis is logarithmic scaling, meaning that the DPEM scales exponentially
in condition number with respect to qubits. However, this does not mean it
scales exponentially with respect to system size; since qubits are related to
system size exponentially, condition number scales linearly with respect to
system size.

##  4 Discussion

The results of our work show that crafting VQAs under the assumption of all-
to-all connectivity of qubits, utilizing their high-entanglement capabilities
for both ansatz and encoding on quantum computers, may be a promising path
forward towards the goal of near-term, practical quantum algorithms. Our paper
proposes a novel ansatz, the GEA, which offers many computational improvements
compared to most ansatz derived in the literature. With the entanglement of
all qubits with each other at each layer, the expressibility of our ansatz is
greatly increased without the need for more rotation gates. This makes it
easier for an optimizer to converge to a solution, since the parameter space
has far less dimensions compared to typical hardware-efficient ansatz. We saw
that, while both the GEA and HEA had linear scaling of parameter space with
respect to qubits, the GEA had a far better linear coefficient, as seen in
Figure 3. The numerical simulations we conducted confirms our conjecture that
the reduction in parameter space enables the optimizer to converge far quicker
with the GEA, as shown in Figure LABEL:fig:0_01_qubits. While both the GEA and
HEA had linear scaling of iterations and quantum circuit evaluations, the GEA
has a far superior linear coefficient, meaning that the time-to-solution would
be considerably quicker when running the algorithm with a GEA on a real
quantum computer (see Figure LABEL:fig:0_01_metrics).

Additionally, our example of the HED for the DPEM shows that using high-
entanglement methods to represent matrices or other, more complicated
structures on a quantum computer may be the key to unlocking better quantum
algorithms, with faster speedups than before. This paper showed that the HED
offers a constant-scaling decomposition of the DPEM compared to the Pauli
decomposition, which is $O(2^{n})$, while remaining under a reasonable circuit
depth scaling of $O(n^{2})$. Note that some other papers have been written on
reducing the decomposition size of the DPEM, but our paper offers unique
insights and better decompositions than these other papers. [26] proposes a
decomposition that has $O(logn)$ sub-matrices, while our decomposition is
$O(1)$. [27] proposes another $O(1)$ decomposition, but the authors do not
discuss the circuit depth requirements of their decomposition. For example,
the shift operator is mentioned, but the consequences on the required circuit
depth for cost function evaluations and time-to-solution metrics are not
examined closely. Our paper indeed examines the bounds on circuit depth that
come as a trade-off when implementing our $O(1)$ decompositions and the
resultant effect on the practicality of this algorithm.

One thing we did not consider in this paper (nor was considered in [13],
though it was mentioned) is the difficulties that arise due to precision
limits on cost function evaluations due to the nature of quantum computing
being rooted in probability. A direct result of the Chebyshev inequality is
that shot count scales as $O(\delta^{-2})$, which is discussed in [28], where
$\delta$ denotes the precision of the cost function evaluation. This means
that in order to allow the cost function to converge to desired $\hat{C}_{L}$
values on a real quantum computer, one would need a significant number of
shots to reach a sufficient precision required for such convergence. The
scaling of cost function precision with respect to digits of precision in
base-10 is $\delta=O(10^{-d})$ where $d$ is the number of digits of precision,
so the scaling of shot count with respect to digits of precision is
$O(10^{2d})$. Since the scaling of digits of precision required for
$\hat{C}_{L}$ scales roughly linearly with respect to number of qubits
defining a system size, the scaling of shot count with respect to required
digits of precision is $O(10^{2O(n)})$. While the exponential system size
dependence on $n$ reduces this scaling back to a polynomial, this still means
that to get a trace precision on the order of $\epsilon=0.01$, as we have done
in this paper, our convergence threshold value, $\hat{C}_{L}$, would be on the
order of $10^{-8}$. Thus, the necessary cost function precision would be
$\delta<10^{-8}$ for any system with $n>3$, meaning at the least $10^{16}$
shots per cost function evaluation would be required to complete a
theoretically perfect run of the algorithm.

This can be remedied using a procedure known as Quantum Mean Estimation
proposed in [29], which would turn the scaling of shot count from exponential
to constant, at the cost of an exponential (but quadratically better) scaling
of circuit depth. Given these two directions, the authors of this paper
suggest that the applicability of this algorithm (and other VQAs that follow a
similar structure) will depend on the whether at least one of the following
capabilities of quantum computers are improved: coherence times in quantum
circuits or shot evaluation time. In the event that one or both of these
capabilities are realized, many other quantum algorithms would also see their
practicality greatly increased [8, 30, 31, 32].

These are only a few of the many potential ideas that may unlock the
capabilities of high-entanglement quantum computers such as the IonQ Aria.
Other applications of high-entanglement capabilities include the usage of GEA
in other quantum machine learning tasks as outlined in [33, 34] or reducing
the number of ancilla qubits required for decomposing Quantum Phase Estimation
(QPE) routines into native gates [32] (such reductions are examined in papers
such as [35] without the use of high-entanglement). A similar procedure may be
done with routines such as Quantum Amplitude Estimation [36] or Quantum Mean
Estimation [29]. Future work is still necessary to be done in order to examine
more closely other ways in which high-entanglement capabilities may be used to
speed-up quantum computing algorithms, whether through a reduction of circuit
depth, simplification of encoding quantum states, or shortening convergence
times of other heuristic algorithms.

###  4.1 Acknowledgements

The authors of this paper would like to thank the team at IonQ and QLab that
gave us the resources to support our work, both for computing capabilities and
knowledge from their staff. I would specifically like to thank Dr. Franz J.
Klein and Dr. Daiwei Zhu for the invaluable guidance given to the authors
about this project. We would also like to thank the Johns Hopkins
Undergraduate Society for Applied Math for inviting us to present our work at
the 2024 Mid-Atlantic Research Exchange (MATRX) in both a lecture and poster-
presentation.

###  4.2 Statement of Funding

The research conducted that culminated into the contents of this paper was
funded by the IonQ Seed Grant for $50k, which was granted to Dr. James Baeder
(PI) and Rahul Babu Koneru (Co-PI) for use in the project titled “Variational
Quantum Linear Solver for Computational Fluid Dynamics”.

## References

  * [1] RP Drake, CC Kuranz, AR Miles, HJ Muthsam, and T Plewa.  Stellar explosions, instabilities, and turbulence.  Physics of Plasmas, 16(4):041004, 2009.
  * [2] Siamak N Doost, Dhanjoo Ghista, Boyang Su, Liang Zhong, and Yosry S Morsi.  Heart blood flow simulation: a perspective review.  Biomedical engineering online, 15(1):1–28, 2016.
  * [3] Joe D. Hoffman and Steven Frankel.  Numerical Methods for Engineers and Scientists.  CRC Press, 2nd edition, May 2001.
  * [4] Mohammad Asif Zaman.  Numerical solution of the poisson equation using finite difference matrix operators.  Electronics, 11(15), 2022.
  * [5] Kun Zhang, Liangbi Wang, and Yuwen Zhang.  Improved finite difference method with a compact correction term for solving poisson’s equations.  Numerical Heat Transfer, Part B: Fundamentals, 70(5):393–405, September 2016.
  * [6] Zhongyang Zhang, Ling Zhang, Ze Sun, Nicholas Erickson, Ryan From, and Jun Fan.  Solving poisson’s equation using deep learning in particle simulation of pn junction.  In 2019 Joint International Symposium on Electromagnetic Compatibility, Sapporo and Asia-Pacific International Symposium on Electromagnetic Compatibility (EMC Sapporo/APEMC), pages 305–308, 2019.
  * [7] Lionel Cheng, Ekhi Ajuria Illarramendi, Guillaume Bogopolsky, Michael Bauerheim, and Benedicte Cuenot.  Using neural networks to solve the 2d poisson equation for electric field computation in plasma fluid simulations, 2021.
  * [8] Aram W. Harrow, Avinatan Hassidim, and Seth Lloyd.  Quantum algorithm for solving linear systems of equations.  Phys. Rev. Lett., 15(103), 2008.
  * [9] Kamal K. Saha, Walter Robson, Connor Howington, In-Saeng Suh, Zhimin Wang, and Jaroslaw Nabrzyski.  Advancing algorithm to scale and accurately solve quantum poisson equation on near-term quantum hardware, 2022.
  * [10] Walter Robson, Kamal K. Saha, Connor Howington, In-Saeng Suh, and Jaroslaw Nabrzyski.  Advanced quantum poisson solver in the nisq era, 2022.
  * [11] Antonio D. Corcoles, Abhinav Kandala, Ali Javadi-Abhari, Douglas T. McClure, Andrew W. Cross, Kristan Temme, Paul D. Nation, Matthias Steffen, and Jay M. Gambetta.  Challenges and opportunities of near-term quantum computing systems.  Proceedings of the IEEE, 108(8):1338–1352, August 2020.
  * [12] M. Cerezo, Andrew Arrasmith, Ryan Babbush, Simon C. Benjamin, Suguru Endo, Keisuke Fujii, Jarrod R. McClean, Kosuke Mitarai, Xiao Yuan, Lukasz Cincio, and Patrick J. Coles.  Variational quantum algorithms.  Nature Reviews Physics, 3(9):625–644, August 2021.
  * [13] Carlos Bravo-Prieto, M. Cerezo Ryan LaRose, Yiğit Subaşı, Lukasz Cincio, and Patrick J. Coles.  Variational quantum linear solver.  Quantum, 7:1188, 2023.
  * [14] Lorenzo Leone, Salvatore F. E. Oliviero, Lukasz Cincio, and M. Cerezo.  On the practical usefulness of the hardware efficient ansatz, 2022.
  * [15] IonQ Staff.  Ionq aria: Past and future, 2024.
  * [16] Enrico Cappanera.  Variational Quantum Linear Solver for Finite Element Problems: a Poisson equation test case.  PhD thesis, TU Delft, 2021.
  * [17] Hai-Ling Liu, Yu-Sen Wu, Lin-Chun Wan, Shi-Jie Pan, Su-Juan Qin, Fei Gao, and Qiao-Yan Wen.  Variational quantum algorithm for the poisson equation.  Physical Review A, 104(2), August 2021.
  * [18] IonQ Staff.  Ionq aria: Practical performance, 2024.
  * [19] Mazen Ali and Matthias Kabel.  Performance study of variational quantum algorithms for solving the poisson equation on a quantum computer.  Phys. Rev. Appl., 20:014054, Jul 2023.
  * [20] IonQ Staff.  Best practices for using ionq hardware, 2024.
  * [21] Dorit Aharonov, Vaughan Jones, and Zeph Landau.  A polynomial quantum algorithm for approximating the jones polynomial.  STOC 2006, 2005.
  * [22] Marek Kowalik.  Parallelization for quantum computers — multi-programming:.  Medium, 2022.
  * [23] Daiwei Zhu.  Personal Communication.  Various conversations about VQLS, VQAs, and IonQ hardware.
  * [24] IonQ Staff.  How to run your first program on ionq cloud, 2024.
  * [25] Aram W. Harrow, Avinatan Hassidim, and Seth Lloyd.  Quantum algorithm for linear systems of equations.  Physical Review Letters, 103(15), October 2009.
  * [26] Hai-Ling Liu, Yu-Sen Wu, Lin-Chun Wan, Shi-Jie Pan, Su-Juan Qin, Fei Gao, and Qiao-Yan Wen.  Variational quantum algorithm for the poisson equation.  Physical Review A, 104(2), August 2021.
  * [27] Yuki Sato, Ruho Kondo, Satoshi Koide, Hideki Takamatsu, and Nobuyuki Imoto.  Variational quantum algorithm based on the minimum potential energy for solving the poisson equation.  Physical Review A, 104(5), November 2021.
  * [28] Alberto Peruzzo, Jarrod McClean, Peter Shadbolt, Man-Hong Yung, Xiao-Qi Zhou, Peter J. Love, Alán Aspuru-Guzik, and Jeremy L. O’Brien.  A variational eigenvalue solver on a photonic quantum processor.  Nature Communications, 5(1), July 2014.
  * [29] Prasanth Shyamsundar.  Non-boolean quantum amplitude amplification and quantum mean estimation, 2021.
  * [30] Peter W. Shor.  Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer.  SIAM Journal on Computing, 26(5):1484–1509, October 1997.
  * [31] Lov K. Grover.  A fast quantum mechanical algorithm for database search.  In Proceedings of the twenty-eighth annual ACM symposium on Theory of computing - STOC ’96, STOC ’96. ACM Press, 1996.
  * [32] A. Yu. Kitaev.  Quantum measurements and the abelian stabilizer problem, 1995.
  * [33] Jacob Biamonte, Peter Wittek, Nicola Pancotti, Patrick Rebentrost, Nathan Wiebe, and Seth Lloyd.  Quantum machine learning.  Nature, 549(7671):195–202, September 2017.
  * [34] Carlo Ciliberto, Mark Herbster, Alessandro Davide Ialongo, Massimiliano Pontil, Andrea Rocchetto, Simone Severini, and Leonard Wossnig.  Quantum machine learning: a classical perspective.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 474(2209):20170551, January 2018.
  * [35] Hongkang Ni, Haoya Li, and Lexing Ying.  On low-depth algorithms for quantum phase estimation.  Quantum, 7:1165, November 2023.
  * [36] Gilles Brassard, Peter Høyer, Michele Mosca, and Alain Tapp.  Quantum amplitude amplification and estimation, 2002.

## Appendix

##  Appendix A 6-qubit Circuit Representations

A. 0: The final circuit that represents the application of the $L_{2}$ matrix
to the ansatz.

A. 1: Quantum circuit representing $L_{1}$ A. 2: Quantum circuit representing
$L_{3}$

##  Appendix B Remaining Numerical Simulations


# Self-organized bistability on scale-free networks

Nikita Frolov1,2<EMAIL_ADDRESS>Alexander Hramov1,2<EMAIL_ADDRESS>1Center for Neurotechnology and Machine Learning, Immanuel Kant Baltic Federal
University, Kaliningrad, Russia
2Institute of Information Technology, Mathematics and Mechanics, Lobachevsky
State University of Nizhny Novgorod, Nizhny Novgorod, Russia

###### Abstract

A dynamical system approaching the first-order transition can exhibit a
specific type of critical behavior known as self-organized bistability (SOB).
It lies in the fact that the system can permanently switch between the
coexisting states under the self-tuning of a control parameter. Many of these
systems have a network organization that should be taken into account to
understand the underlying processes in detail. In the present paper, we
theoretically explore an extension of the SOB concept on the scale-free
network under coupling constraints. As provided by the numerical simulations
and mean-field approximation in the thermodynamic limit, SOB on scale-free
networks originates from facilitated criticality reflected on both macro- and
mesoscopic network scales. We establish that the appearance of switches is
rooted in spatial self-organization and temporal self-similarity of the
network’s critical dynamics and replicates extreme properties of epileptic
seizure recurrences. Our results, thus, indicate that the proposed conceptual
model is suitable to deepen the understanding of emergent collective behavior
behind neurological diseases.

††preprint: APS/123-QED

##  I Introduction

Many natural and man-made systems deliver diverse functionality by switching
between coexisting stable states, or multistability. In living systems, this
phenomenon regulates the processes on different scales: from the interaction
of inner organ systems to cell cycles and neuronal synchronization on
considerably smaller scales [1, 2, 3]. Specifically, the ability to switch
between the multiple patterns of local coherence on the brain’s cortical
network underpins many aspects of consciousness [4]. Past works have
identified that in normal conditions, emergent neuronal avalanches expose
distinct power-law scaling properties, which are the hallmarks of self-
organized criticality (SOC). [5, 6, 7, 8].

On the contrary, excessive synchronization involving large-scale networks into
coherent motion is an example of undesired and harmful behavior, which is the
case of the brain’s epileptic condition [9, 10, 11]. Such extreme bursts of
coherence lie beyond the power-law distribution and are most probably rooted
in the bistable dynamics of neuronal ensembles. In this regard, Di Santo et
al. have extended the theory of SOC and have described these dynamics as a
self-organized bistability (SOB) [12]. They have proposed a conceptual model
of a system in the vicinity of the first-order transition – a singular point
at which both stable states are simultaneously present. The switching between
these states is provided by self-tuning the control parameter under the impact
of dissipation and driving force. Remarkably, the observed avalanches obey the
statistics of dragon-kings (DKs), which are rare and significantly large
events. At the same time, the occurrence of DKs originates from the system’s
nonlinearity. Therefore these events are not random and even predictable to
some extent [13]. Overall, these facts confirm the fundamental significance of
the SOB concept in non-trivial neural dynamics under epileptic conditions.

However, to provide further insight into the origins of bistability on damaged
neuronal populations, one has to account for their networked organization
[14]. It is well-known that the interplay between structure and dynamics
enriches the available modes of collective behaviors [15, 16]; Thus, it is
essential to understand how the structural properties drive SOB on complex
networks and give birth to excessively large avalanches.

In our latest work [17], we have introduced a conceptual network model under
coupling constraints that should be considered a network extension of the SOB
theory. In this model, the dissipation and driving of the interelement
coupling self-tuning, as suggested by SOB, is biologically motivated by the
general aspects of neuro-astrocyte interaction. Its dysregulation provides a
metabolic mechanism for astrocyte-dependent hyperexcitability [18] and is
recognized as one of the driving factors of epilepsy [19, 18]. Another feature
of the model that needs to be emphasized is the exploitation of a self-similar
scale-free (SF) structure that determines the bistability of collective
behavior, i.e., an explosive transition, under certain conditions [20, 21,
22]. Although the SF organization of local neuronal populations is debatable
and has not been proven so far, there are several indirect signs and arguments
based on the underlying criticality of neural dynamics in support of this
hypothesis [8, 23, 24].

Consistent with SOB, our model has replicated the appearance of epileptic
seizure generation and self-termination of a supercritical synchronized state
due to the interplay between the explosive synchronization and coupling
constraints. We referred to such a state as an extreme synchronization event
due to its spontaneous occurrence and short lifetime. Although we have
demonstrated that the designed model produces solitary synchronization events,
the critical transitions leading to such a hazardous behavior have not been
fully addressed.

In this work, we exploit extensive theoretical analysis to derive the
conditions of SOB for the proposed network under coupling constraints. We
primarily focus on the Barabási-Albert (BA) model, a paradigmatic model of SF
network grown through the preferential attachments [25]. This model yields a
degree distribution with constant exponential decay of $-3$. For this model,
we derive the mean-field approximation in the thermodynamic limit that is in
good agreement with the numerical simulations of a finite-size network. Having
analyzed continuous synchronization diagrams and statistical properties of the
synchronization events occurrence, we show an essential role of self-organized
bistability in producing and terminating short-term synchronization. We
establish that these events reproduce the statistical properties of extreme
value theory and epileptic brain activity. Finally, we generalize our analysis
for SF networks with an arbitrary scaling exponent and reveal that such a
structural property has a particular impact on the emergence of SOB and its
temporal self-similarity.

##  II Model

Consider a network of $N$ Kuramoto phase oscillators, whose rotation is
defined by the system of differential equations

$\displaystyle\frac{d\theta_{i}}{dt}$
$\displaystyle=\omega_{i}+\lambda_{i}\sum_{j=1}^{N}\mathcal{A}_{ij}\sin(\theta_{j}-\theta_{i}),$
(1) $\displaystyle\frac{d\lambda_{i}}{dt}$
$\displaystyle=\alpha(\lambda_{0}-\lambda_{i})-\beta r_{i}.$ (2)

In Eq. (1), $\theta_{i}$ and $\omega_{i}$ denote instantaneous phase and
natural frequency of each $i^{th}$ Kuramoto oscillator, $i=1,...,N$. The
adjacency matrix $\mathcal{A}_{ij}$ defines the connectivity of the SF network
yielding a power-law degree distribution $P(k)\sim k^{-\gamma}$. To construct
$\mathcal{A}_{ij}$ we use the BA model that generates a degree-degree
correlated SF network based on preferential attachment [25] yielding a degree
distribution with fixed scaling exponent $\gamma=3$. In the BA model, $m=3$
new edges preferentially attach between a new vertex and existing vertices at
each step of the growing process. To explore the influence of scaling
exponent, we employ the Chung-Lu (CL) that generates an uncorrelated SF
network with an arbitrary pre-defined value of scaling exponent $\gamma$ [26,
27]. The CL model is described in detail in Appendix A. In both models, we
adjust $\omega_{i}=k_{i}$ to fulfill the frequency-degree correlation.
Macroscopic motion of the ensemble is defined by the complex order parameter
$Re^{i\Psi}=1/N\sum_{j=1}^{N}e^{i\theta_{j}}$.

Following Di Santo et al. [12], Eq. (2) describes the self-tuning of
individual coupling strength $\lambda_{i}$. Here, coupling consumption (a
dissipation) is a function of local order parameter
$r_{i}=1/k_{i}\lvert\sum_{j=1}^{N}\mathcal{A}_{ij}e^{i\theta_{j}}\rvert$,
which implies that the maintenance of a certain level of coherence in the
neighborhood of the $i^{th}$ unit consumes its coupling ability at rate
$\beta$. Besides, each unit diffusely couples to an external source that
provides the recovery of individual coupling ability to the level
$\lambda_{0}$ at rate $\alpha$ (a driving force).

To provide deeper insights into the dynamics of Eq. (1)–(2) concerning the BA
model, we used the mean-field formalism proposed by Ichinomiya [28].
Considering the motion of the ensemble at the mean frequency
$\Omega=\langle\omega\rangle=\langle k\rangle$ with the average phase $\Psi$,
we can replace $\phi=\theta-\Psi$. In a continuous limit, the density of
oscillators having phase $\phi$ and degree $k$ at time moment $t$ is defined
by $\rho(\phi,k,t)$, such that $\int_{0}^{2\pi}d\phi\rho(\phi,k,t)=1$, that
obeys the continuity equation:

$\frac{\partial\rho}{\partial
t}+\frac{\partial}{\partial\phi}\left\\{v_{\phi}\rho\right\\}=0,$ (3)

Eq. (3) therefore describes the motion of the networked system (1)-(2) in the
continuous limit with $v_{\phi}=d\phi/dt$. Let us define the global order
parameter through the density $\rho(\phi,k,t)$ in the integral form

$Re^{i\Psi}=\frac{1}{\langle k\rangle}\int_{m}^{k_{max}}dk\int_{0}^{2\pi}d\phi
kP(k)\rho e^{i\phi},$ (4)

and rewrite Eq. (1) using the obtained definition (4) to find the mean-field
approximation of $v_{\phi}$:

$v_{\phi}=\frac{d\phi}{dt}=(1-\lambda R\sin\phi)k-\langle k\rangle.$ (5)

In Eq. (5), we suggest that $\lambda_{i}=\lambda$ and $r_{i}=R$. Given (4) and
(5), Peron and Rodrigues found the steady-state solution of Eq. (3) for a
frequency-degree correlated SF network [29]:

$\rho(k,\phi)=\begin{cases}\delta\left(\phi-\arcsin\left[\frac{k-\langle
k\rangle}{\lambda Rk}\right]\right),\quad\lvert k-\langle
k\rangle\rvert\leq\lambda Rk\\\ \frac{A(k)}{\lvert(k-\langle k\rangle)-\lambda
Rk\sin\phi\rvert},\quad\text{otherwise},\end{cases}$ (6)

with normalization term $A(k)$.

The condition $\lvert k-\langle k\rangle\rvert\leq\lambda Rk$ is fulfilled for
the phase-locked fraction of networked oscillators while remaining possess
incoherent drift. Since only phase-locked units contribute to the global
coherence of the ensemble, Eq. (4) becomes

$R=\frac{1}{\langle k\rangle}\int_{\lvert k-\langle k\rangle\rvert\leq\lambda
Rk}kP(k)\sqrt{1-\left(\frac{k-\langle k\rangle}{\lambda Rk}\right)^{2}}dk.$
(7)

To complete a mean-field approximation of the system (1)-(2) we should
complement Eq. (7) with account of coupling constraints. In the steady state
$d\lambda/dt=0$, it becomes

$\lambda=\lambda_{0}-\frac{\beta}{\alpha}R.$ (8)

To provide rigorous conclusions, we present Eqs. (7)-(8) in the thermodynamic
limit $N\rightarrow\infty$ by explicitly taking the exact degree distribution
$P^{\infty}(k)=2m(m+1)/(k(k+1)(k+2))$ and mean degree $\langle
k\rangle^{\infty}=2m$ for BA model [30, 31]. Further, we denote a mean-field
approximation of the order parameter as $R^{\infty}$.

Using mean-field analysis one can also estimate conditions for the trivial
solution of Eq. (7), $R^{\infty}=0$, to lose its stability. By letting
$R^{\infty}\rightarrow+0$, Peron and Rodrigues derived this critical point,
which in thermodynamic limit is $\lambda_{c}=1/(\pi mP^{\infty}(2m))\simeq
1.485$ [29]. Accounting for coupling constraints we define the critical point
for our model (1)-(2) as:

$\lambda_{0c}=\frac{1}{\pi mP^{\infty}(2m)}+\frac{\beta}{\alpha}R^{0},$ (9)

where $R^{0}=1/\sqrt{N}$ is a finite estimation of the order parameter in the
incoherent state.

Figure 1: (a) (a) Forward and backward synchronization diagrams without
excitability constraints $\beta=0$. (b) A forward diagram in the presence of a
finite excitability consumption $\beta=0.04$. Synchronization diagrams in both
panels are computed for fixed $\alpha=0.04$. In panel (b), the median value,
the maximum value, and the inter-quartile interval (IQI) of $R(t)$ were
computed over the time series of $2\times 10^{4}$ time units. In both panels,
the black solid line and vertical dashed line show a mean-field approximation
$R^{\infty}$ and the critical point $\lambda_{0c}$, respectively.

##  III Results

We start with considering forward and backward continuations of a described
network by increasing and decreasing $\lambda_{0}$ using increment
$\delta\lambda_{0}=0.02$. Eqs. (1)-(2) for a finite-size network are
integrated for $N=1000$ unit.

Fig. 1 reports synchronization diagrams $R(\lambda_{0})$ in the absence and
presence of coupling constraint for fixed $\alpha=0.04$. As expected, without
consumption $\beta=0$ (Fig. 1a), our system exhibits an explosive transition
reproducing a hysteresis loop demonstrated in the pioneering work on a
frequency-degree correlated SF network by Gómez-Gardeñez et al. [20]. A mean-
field approximation also predicts an explosive transition indicating a
bistability region, where both $R^{\infty}=0$ and $R^{\infty}\neq 0$ are
stable ($\lambda_{0c}\simeq 1.485$).

Numerical results for a system under coupling consumption $\beta=0.04$ show
that the onset of forward transition is largely delayed concerning the
critical point $\lambda_{0c}$. Fig. 1b displays that both coherent and
incoherent solutions are possible in the region between $\lambda_{0c}$ and the
onset of forward transition. It follows from Eq. (2) that the stable coherent
state is only possible if coupling consumption is compensated by its recovery
$\alpha(\lambda_{0}-\lambda_{i})\geq\beta r_{i}$ for each network unit. Since
the latter condition is not fulfilled in the area between $\lambda_{0c}$ and
the onset of forward transition, the turbulent state $R\simeq 0$ here remains
preferable for the network. However, rare transitions to synchronized state
shown in Fig. 1b) as the outliers of $R(t)$ distribution are possible due to
$\lambda_{0}>\lambda_{0c}$.

Figure 2: (a)-(c) Phase portraits of the macroscopic dynamics of the SF
network under coupling constraints on the plane $(R,\langle\lambda\rangle)$
for fixed values of $\beta=0.04$, $\alpha=0.04$ and different $\lambda_{0}$:
(a) $\lambda_{0}=1.5$ – subcritical dynamics; (b) $\lambda_{0}=2.3$ –
supercritical dynamics; (c) $\lambda_{0}=2.0$ – critical dynamics. Dashed
black and gray lines in (a)-(c) display the forward and backward
synchronization diagrams $R(\lambda)$ for $\beta=0$. Corresponding phase
portraits on the complex plane in the polar coordinate system $(R,\Psi)$ are
presented in (d)-(f).

Based on Fig. 1b, three characteristic areas can be distinguished: (i) the
region of subcritical dynamics $\lambda_{0}<\lambda_{0c}$; (ii) the region of
critical bistable dynamics between $\lambda_{0c}$ and the onset of forward
transition; (iii) the region of supercritical dynamics after forward
transition corresponding to a coherent motion.

Now, we illustrate the network’s dynamics in terms of its macroscopic
parameters under variation of $\lambda_{0}$. Fig. 2 reports macroscopic motion
on the parameter plane $(R,\langle\lambda\rangle)$, where
$\langle\lambda\rangle=1/N\sum_{i=1}^{N}\lambda_{i}$ is the ensemble-average
level of coupling ability, accompanied by corresponding trajectories on the
complex polar plane $(R,\Psi)$. As expected, in a subcritical region, the
system converges to a turbulent drift around the fixed point at a lower branch
of the hysteresis curve (Fig. 2a,d). In a supercritical mode, it is attracted
to an upper branch that corresponds to the limit cycle at mean frequency
$\Omega=\langle k\rangle$ (Fig. 2b,e).

Figure 3: Distribution of mesoscopic parameters – local order parameter
$r_{i}$ (a)-(c) and individual coupling level $\lambda_{i}$ (d)-(e) – versus
node’s degree $k_{i}$. The distributions are presented for fixed
$\alpha=0.04$, $\beta=0.04$ and different values of $\lambda_{0}$: (a),(d)
$\lambda_{0}=1.1$ – subcritical behavior; (b),(e) $\lambda_{0}=1.6$ – weak
criticality; (c),(f) $\lambda_{0}=1.9$ – strong criticality. Black lines
indicate approximations $r_{i}=1/\sqrt{k_{i}}$ in left column and
$\lambda_{i}=\lambda_{0}-\beta/(\alpha\sqrt{k_{i}})$ in right column.

Remarkably, in the region of critical bistable dynamics, the system mostly
drifts near the fixed point $R\simeq 0$ in the vicinity of transition point
$\lambda_{c}\simeq 1.485$. However, the turbulent drift may occasionally push
the trajectory beyond the critical point. Hence, the network explosively
synchronizes as reflected by motion along the repelling trajectory toward the
upper branch (Fig. 2c). Having landed the upper branch the network moves to
backward transition forced by resource consumtion. As the network exhausts the
coupling resource, it undergoes the backward transition. Eventually, the
diffusive mechanism of coupling recovery returns the trajectory to the
starting position near the forward transition point on the lower branch of the
hysteresis curve. The reported process repeats over time as the coupling
becomes excessively strong $\langle\lambda\rangle>\lambda_{c}$ due to chaotic
fluctuations in the vicinity of the critical point.

The criticality of the considered Kuramoto model under coupling constraints
that determines the spontaneous appearance of synchronization events is
reflected on the mesoscopic level, i.e., on the level of local populations.
Log-log plots in Fig. 3 show the distributions of local order parameter
$r_{i}(k_{i})$ and individual coupling ability $\lambda_{i}(k_{i})$ during a
long-term period of turbulent motion ($2\times 10^{4}$ time units). In a
subcritical regime $\lambda_{0}=1.1$, local order parameter is well-
approximated by inverse square root law $r_{i}=1/\sqrt{k_{i}}$ (Fig. 3a). It
indicates that below the critical value $\lambda_{0c}=1.517$, phases within
the local groups are homogeneously distributed on the unit circle, implying
turbulent motion across all the mesoscopic scales. Substituting
$r_{i}=1/\sqrt{k_{i}}$ in Eq. (2) and letting $d\lambda_{i}/dt=0$, one finds
the stationary distribution of individual coupling ability in a subcritical
regime $\lambda_{i}=\lambda_{0}-\beta/(\alpha\sqrt{k_{i}})$ that perfectly
fits the empirical distribution at $\lambda_{0}=1.1$ (Fig. 3d).

Figure 4: Synchronization diagrams $R(\lambda_{0})$ for the fixed value
$\beta=0.04$ and different $\alpha$: (a) $\alpha=0.04$; (b) $\alpha=0.08$; (c)
$\alpha=0.12$. (d) Phase diagram in the $(\alpha,\lambda_{0})$ parameter plane
plotted upon the mean-field approximation given by Eqs. (7)–(8).

At $\lambda_{0}=1.6$ slightly exceeding the critical value $\lambda_{0c}$
(weak criticality, Fig. 3b,e), inverse square root approximations
$r_{i},\lambda_{i}\sim 1/\sqrt{k_{i}}$ are applicable only to large-scale
populations formed around the network’s hubs. At the same time, the local
coherence of peripheral groups around low- and medium-degree units starts
rising above the estimate given by this approximation. With increasing
criticality (Fig. 3,c,f), local coherence around the low- and medium-degree
nodes grows while hubs remain desynchronized. Taking into account that
$R\simeq 1/\sqrt{N}$ in all considered cases, facilitation of local coherence
does not contribute to global synchronization. Instead, the critical dynamics
of this model form a complex of populations synchronized on different spatial
and temporal scales and weakly interacting with each other. Noteworthy, due to
the coupling consumption, the coherence of small groups is not long-term
preserved but increases and decreases over time independently of other groups,
thus resisting a transition to global synchronization. However, as soon as a
large number of local populations with similar frequencies accidentally
synchronize, at the same time, critical dynamics create conditions for
synchronizing key components – the network’s hubs. Involvement of hubs in the
coherent motion launches an avalanche-like process of a large-scale phase-
locking resulting in an abrupt first-order transition to metastable global
synchrony.

Figure 5: Long-term time series $R(t)$ (a)-(c) and corresponding log-log
histograms (d)-(f) in the proximity of forward transition for fixed
$\alpha=0.04$, $\beta=0.04$ and different $\lambda_{0}$: (a),(d)
$\lambda_{0}=1.93$; (b),(e) $\lambda_{0}=1.95$; (c),(f) $\lambda_{0}=1.99$.
The orange line in (d)-(f) shows a fitted power law. Figure 6: Long-term time
series $R(t)$ (a)-(c) and corresponding distributions of return intervals
(d)-(f) in the proximity of forward transition for fixed $\alpha=0.04$,
$\beta=0.04$ and different $\lambda_{0}$: (a),(d) $\lambda_{0}=1.97$; (b),(e)
$\lambda_{0}=1.99$; (c),(f) $\lambda_{0}=2.03$. Black lines in (d)-(f) show a
fitted power law $\sim\tau^{-3/2}$. The outcomes of the $\chi^{2}$-test are
presented within each panel.

In addition, we demonstrate how the balance between recovery rate $\alpha$ and
consumption rate $\beta$ affects the critical dynamics of the considered
model. Fig. 4 reports the synchronization diagrams $R(\lambda_{0})$ for
different values of $\alpha$ at fixed $\beta=0.04$. It follows from Eqs.
(7)–(8) that increasing $\alpha$, or equivalently decreasing $\beta$, reduces
the impact of coupling constraints on the system’s dynamics. Indeed, the
region of critical behavior bounded by the critical value $\lambda_{0c}$ from
the left and the onset of forward transition from the right shrinks with
increasing recovery rate $\alpha$ (Fig. 4a,b) until the transition to
coherence becomes of a first-order (Fig. 4c). This observation for a finite-
size network is consistent with the mean-field approximation in the
thermodynamic limit $N\rightarrow\infty$ (Fig. 4d).

Up to now, we have considered the emergence of such synchronization events and
associated critical dynamics of the Kuramoto network under coupling
constraints. Now, we explore the statistics of their appearance in the
proximity of forward transition in terms of extreme events theory. Let us fix
$\alpha=0.04$, $\beta=0.04$ and consider the Kuramoto ensemble under coupling
constraints approaching the point of forward transition. Although the area of
criticality is broad enough, as displayed in Fig. 1b and Fig. 4a, the
generation of solitary synchronization events is quite rare and occasional in
a long-term perspective, not exceeding one event per $2\times 10^{4}$ time
units. However, the probability of their emergence drastically increases in
the proximity of forward transition point $\lambda_{0}\simeq 2.05$ (Fig. 5),
allowing us to explore the statistical properties of these events.

Fig. 5a-c show the time series $R(t)$ and corresponding histograms (Fig. 5d-f)
with increasing $\lambda_{0}$. As $\lambda_{0}$ approaches the forward
transition point, the number of produced metastable coherent events increases.
Analysis of the histograms reveals that the variation of $R$ near the
preferable incoherent state obeys power law $P(R)\sim R^{-\zeta}$. Although
the power law is fitted in a relatively narrow region of $R$ values near the
stable fixed point, this observation could be considered a hallmark of self-
organized criticality [7]. It implies that the occurrence of small- and
medium-size events, i.e., the short-term establishment of coherence in small-
and medium-size groups, obeys scaling of the same physical mechanism. Although
the dynamics seem to be highly turbulent at first glance, the power-law fit of
the tail indicates the presence of spatial order. At the same time, a peak
beyond the power-law fit indicates a non-negligible probability of large
deviations of $R$. Such events lie outside the scaling rule and obey different
physical principles. Such distribution is a hallmark of a specific type of
extreme events – dragon kings, DKs [13, 32]. DKs possess a remarkable property
– these events are significant and non-random; therefore, they are predictable
to some degree. Indeed, DKs are generated in a deterministic system and obey
certain mechanisms underlying critical behavior in the vicinity of the typing
point. The possibility of predicting such states, although tangible, is a
challenging task that requires a deep understanding of the structure and
dynamics of the considered networked system.

Figure 7: Forward and backward synchronization diagrams in the absence of
excitability constraints $\beta=0$ for different values of $\gamma$: (a)
$\gamma=2.6$, (b) $\gamma=2.8$, (c) $\gamma=3.2$, and (d) $\gamma=3.4$.
(e)-(h) Respective forward diagrams in the presence of a finite excitability
consumption $\beta=0.04$. Synchronization diagrams in both panels are computed
for fixed $\alpha=0.04$. In panels (e)-(h), the median value, the maximum
value, and the inter-quartile interval (IQI) of $R(t)$ were computed over the
time series of $2\times 10^{4}$ time units.

To finalize our analysis of the BA model, we explored the distribution of
return intervals between neighboring synchronization events $\tau$. Fig. 6
displays long-term time series $R(t)$ and corresponding return intervals
distributions $P(\tau)$ for different $\lambda_{0}$ in the proximity of
forward transition. Despite an increase of criticality and changes in
characteristic scales of return intervals at different $\lambda_{0}$, their
distribution is well-fitted by the power law $P(\tau)\sim\tau^{-3/2}$. The
Pearson’s $\chi^{2}$-test, which has been used to assess the goodness of fit,
yields a correspondence between the observed and expected power-law
distributions (test outcomes are presented in each subplot). Similar to the
previous discussion on the histograms $P(R)$, a distinct power-law fitting of
$P(\tau)$ evidences the presence of temporal order in the sequence of
synchronization events despite being aperiodic and visually irregular. This
observation is most probably rooted in the concept of “on-off” intermittency –
an alteration of turbulent (“on”) and laminar (“off”) phases of the system’s
motion in the vicinity of synchronization onset [33]. In “on-off”
intermittency, an analogous power-law fit with characteristic exponent $-3/2$
defines the distribution of laminar phases duration and reflects the approach
to ordered dynamics from turbulence as the system evolves towards
synchronization. The $-3/2$ power law is also reported in empirical
observations of epileptic activity in the rodents’ brain [34, 35, 36, 37]. For
instance, Refs. [34, 35, 36] propose a theoretical framework bonding “on-off”
intermittency and occurrence of epileptic seizures. Ref. [37] explicitly
demonstrates both extreme properties of seizures’ amplitudes and temporal
scaling fitted by the $-3/2$ power law.

Lastly, we investigate how the structure of the SF network affects SOB.
Specifically, we are interested in the influence of scaling exponent $\gamma$
since this parameter regulates the prevalence of highly-connected hubs in SF
networks. The latter units play a principal role in the first-order transition
and, consequently, the emergence of an SOB in our model. Earlier, Coutinho et
al. [21] have developed a theoretical framework to show that in the frequency-
degree correlated networks, the first-order transition is present at
$2<\gamma<3$. Recall that this particular range of scaling exponent is the
most common in real SF networks [38]. For $\gamma>3$, Coutinho et al. have
predicted the suppression of hysteresis that determines the continuous type of
transition.

Using the CL model, we have generated the degree-degree uncorrelated SF
networks with $\gamma=\\{2.6,2.8,3.2,3.4\\}$. Fig. 7 displays the
synchronization diagrams $R(\lambda_{0})$ for these networks in the absence
and presence of coupling constraint for fixed $\alpha=0.04$. Consistent with
the work by Coutinho et al. [21], our results confirm that in the absence of
coupling constraints, the network undergoes an explosive transition at
$\gamma<3$ and a continuous transition at $\gamma>3$ (Fig. 7a-d). As in [21],
one can also see that the hysteresis area shrinks as $\gamma$ increases to
$3$. Introduction of coupling constraints at rate $\beta=0.04$ results in a
distinct SOB behavior for $\gamma<3$. In Fig. 7e,f, one can see the area of
critical dynamics followed by an abrupt stabilization of a globally coherent
state, i.e., a discontinuous transition of median $R(t)$. With a further
increase of $\gamma$, the area of critical dynamics associated with SOB
vanishes as the transition to coherence becomes continuous (Fig. 7g,h).
Actually, for $\gamma>3$, the very concept of SOB itself is meaningless since
the network is no longer bistable. However, there are still at least two
reasons for the SF network to exhibit SOB-like behavior at $\gamma=3.2$ as
displayed by Fig. 7g. First, it could be due to a consideration of a finite-
sized and relatively small network ($N=1000$). In this case, the scaling of
$\gamma=3.2$ could still provide a sufficient amount of highly-connected units
for the system to synchronize rapidly yet continuously. Second, it could be
due to the chosen values of $\alpha$ and $\beta$, i.e., coupling recovery and
consumption rates. With $\alpha=\beta=0.04$, a continuous transition runs much
faster than the self-tuning of coupling strength $\lambda$ that determines
sufficiently rapid but not bistable switches.

Besides the above-discussed effect, the variation of scaling exponent $\gamma$
within the SOB-associated range $2<\gamma<3$ impacts the slope of distribution
of return intervals. Fig. 8 displays such distributions for $\gamma$ equal to
$2.6$, $2.8$ and $3.0$ along with the respective power-law fits. One can see
that decreasing SF scaling exponent $\gamma$ increases the slope of return
intervals distribution. Although the particular reasons for these observations
require further investigation and are yet to be discovered, one may conclude
that the SF network’s structure is closely related to the properties of its
temporal self-similarity. Indeed, increasing the amount of highly-connected
hubs on the SF network might cause more frequent bistable switches while
maintaining the distribution’s exponential slope. At the same time, the slope
of return intervals distribution converges to the universal constant $-3/2$ as
$\gamma\rightarrow 3$ since $\gamma=3$ is the point of critical singularity,
separating the first- and second-order transitions on the frequency-degree
correlated SF networks.

Figure 8: Distributions of return intervals in the proximity of forward
transition for fixed $\alpha=0.04$, $\beta=0.04$ and different scaling
exponents $\gamma$ in the CL model: $\gamma=2.6$, $\lambda_{0}=2.16$ (blue
dots) and $\gamma=2.8$, $\lambda_{0}=2.06$ (red dots). Green dots repeat the
distribution for the BA model with $\gamma=3.0$ and $\lambda_{0}=2.03$. Solid
lines of respective colors show the power-law fit. The outcomes of the
$\chi^{2}$-test are presented for each line.

##  IV Conclusion

To conclude, this paper reports a theoretical study of the self-organized
bistability of a scale-free Kuramoto network under coupling constraints. Our
analysis shows that the imbalance between the rates of coupling consumption
and recovery determines the delay between the forward synchronization and the
critical point, at which the global incoherence becomes unstable. Facilitated
critical dynamics of the networked ensemble in this area provide self-
consistent switching between coexisting states of incoherence and coherence.
Current observations demonstrate that the appearance of such switches bistable
exhibits the pronounced statistical properties of extreme events and
reproduces the features of epileptic seizures recurrences. It evinces a high
degree of self-similarity and ordered dynamics behind the process looking
uncorrelated at first glance. Finally, we establish that emergence of self-
organized bistability on a scale-free is determined by its structural
properties, e.g., in terms of its scaling exponent $\gamma$. We reveal that
this type of critical behavior is observed at $2<\gamma<3$, which is the most
common in real scale-free networks and suppressed at $\gamma>3$. Our
conceptual and yet biologically motivated network model sheds light on the
macro- and mesoscopic properties of collective behavior behind the disruptive
emergent hypersynchronization of the brain networks paving the way for their
study in dynamical and graph-theoretical contexts.

Figure 9: Constructed uncorrelated scale-free networks according to the Chung-
Lu model. Top row presents expected degree profiles $k^{\prime}_{i}(i)$
(orange curves, determined by Eq. (10)) and actual degree profiles of
generated scale-free networks $k_{i}(i)$ (blue dots) for different values of
scaling exponent $\gamma$. The bottom row presents respective degree
distributions $P(k)$ of generated uncorrelated networks (blue dots) along with
the expected slope $\sim k^{-\gamma}$ (orange curves).

###### Acknowledgements.

The work has been supported by the Russian Science Foundation (Grant No.
21-72-10129).

##  Appendix A Generating a scale-free graph with arbitrary exponent

The Chung-Lu model has been used to produce an uncorrelated scale-free graph
with arbitrary scaling exponent [26, 27]. Specifically, we exploited a recent
algorithm proposed by Fasino et al. [39]. Briefly, it generates a random graph
with an arbitrary power-law degree distribution $P(k)\sim k^{-\gamma}$ by
adjusting a specific degree to each $i^{th}$ unit according to the expected
degree profile:

$k^{\prime}_{i}=c(i_{0}+i)^{-p},\qquad p=\frac{1}{1-\gamma}.$ (10)

In Eq. (10), parameters $c$ and $i_{0}$ are defined as:

$c=(1-p)k^{*}N^{p},\quad i_{0}=\left(\frac{c}{k_{max}}\right)^{1/p}-1.$ (11)

Here, $k^{*}$ and $k_{max}$ are mean and maximal expected degree. For our
calculations, we have set $k^{*}=10$ and $k_{max}=100$. Degree distributions
of constructed scale-free networks for different values of scaling exponent
$\gamma$ are presented in Fig. 9.

## References

  * Angeli  _et al._ [2004] D. Angeli, J. E. Ferrell, and E. D. Sontag, Detection of multistability, bifurcations, and hysteresis in a large class of biological positive-feedback systems, Proceedings of the National Academy of Sciences 101, 1822 (2004).
  * Strogatz [2012] S. H. Strogatz,  _Sync: How order emerges from chaos in the universe, nature, and daily life_ (Hachette UK, 2012).
  * Pisarchik and Feudel [2014] A. N. Pisarchik and U. Feudel, Control of multistability, Physics Reports 540, 167 (2014).
  * Koch  _et al._ [2016] C. Koch, M. Massimini, M. Boly, and G. Tononi, Neural correlates of consciousness: progress and problems, Nature Reviews Neuroscience 17, 307 (2016).
  * Linkenkaer-Hansen  _et al._ [2001] K. Linkenkaer-Hansen, V. V. Nikouline, J. M. Palva, and R. J. Ilmoniemi, Long-range temporal correlations and scaling behavior in human brain oscillations, Journal of Neuroscience 21, 1370 (2001).
  * Beggs and Plenz [2003] J. M. Beggs and D. Plenz, Neuronal avalanches in neocortical circuits, Journal of neuroscience 23, 11167 (2003).
  * Bak [2013] P. Bak,  _How nature works: the science of self-organized criticality_ (Springer Science & Business Media, 2013).
  * Shew and Plenz [2013] W. L. Shew and D. Plenz, The functional benefits of criticality in the cortex, The neuroscientist 19, 88 (2013).
  * Louzada  _et al._ [2012] V. H. Louzada, N. A. Araújo, J. S. Andrade, and H. J. Herrmann, How to suppress undesired synchronization, Scientific reports 2, 1 (2012).
  * Lehnertz  _et al._ [2009] K. Lehnertz, S. Bialonski, M.-T. Horstmann, D. Krug, A. Rothkegel, M. Staniek, and T. Wagner, Synchronization phenomena in human epileptic brain networks, Journal of neuroscience methods 183, 42 (2009).
  * Lüttjohann and van Luijtelaar [2015] A. Lüttjohann and G. van Luijtelaar, Dynamics of networks during absence seizure’s on-and offset in rodents and man, Frontiers in physiology 6, 16 (2015).
  * di Santo  _et al._ [2016] S. di Santo, R. Burioni, A. Vezzani, and M. A. Munoz, Self-organized bistability associated with first-order phase transitions, Physical review letters 116, 240601 (2016).
  * Sornette and Ouillon [2012] D. Sornette and G. Ouillon, Dragon-kings: mechanisms, statistical methods and empirical evidence, The European Physical Journal Special Topics 205, 1 (2012).
  * Kalitzin  _et al._ [2019] S. Kalitzin, G. Petkov, P. Suffczynski, V. Grigorovsky, B. L. Bardakjian, F. L. da Silva, and P. L. Carlen, Epilepsy as a manifestation of a multistate network of oscillatory systems, Neurobiology of disease 130, 104488 (2019).
  * Boccaletti  _et al._ [2006] S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, and D.-U. Hwang, Complex networks: Structure and dynamics, Physics reports 424, 175 (2006).
  * Boccaletti  _et al._ [2018] S. Boccaletti, A. N. Pisarchik, C. I. Del Genio, and A. Amann,  _Synchronization: from coupled systems to complex networks_ (Cambridge University Press, 2018).
  * Frolov and Hramov [2021] N. Frolov and A. Hramov, Extreme synchronization events in a kuramoto model: The interplay between resource constraints and explosive transitions, Chaos: An Interdisciplinary Journal of Nonlinear Science 31, 063103 (2021).
  * Siracusa  _et al._ [2019] R. Siracusa, R. Fusco, and S. Cuzzocrea, Astrocytes: role and functions in brain pathologies, Frontiers in Pharmacology , 1114 (2019).
  * Diaz Verdugo  _et al._ [2019] C. Diaz Verdugo, S. Myren-Svelstad, E. Aydin, E. Van Hoeymissen, C. Deneubourg, S. Vanderhaeghe, J. Vancraeynest, R. Pelgrims, M. I. Cosacak, A. Muto, _et al._ , Glia-neuron interactions underlie state transitions to generalized seizures, Nature communications 10, 1 (2019).
  * Gómez-Gardeñes  _et al._ [2011] J. Gómez-Gardeñes, S. Gómez, A. Arenas, and Y. Moreno, Explosive synchronization transitions in scale-free networks, Physical review letters 106, 128701 (2011).
  * Coutinho  _et al._ [2013] B. Coutinho, A. Goltsev, S. Dorogovtsev, and J. Mendes, Kuramoto model with frequency-degree correlations on complex networks, Physical Review E 87, 032106 (2013).
  * Boccaletti  _et al._ [2016] S. Boccaletti, J. Almendral, S. Guan, I. Leyva, Z. Liu, I. Sendiña-Nadal, Z. Wang, and Y. Zou, Explosive transitions in complex networks’ structure and dynamics: Percolation and synchronization, Physics Reports 660, 1 (2016).
  * Haimovici  _et al._ [2013] A. Haimovici, E. Tagliazucchi, P. Balenzuela, and D. R. Chialvo, Brain organization into resting state networks emerges at criticality on a model of the human connectome, Physical review letters 110, 178101 (2013).
  * Gastner and Ódor [2016] M. T. Gastner and G. Ódor, The topology of large open connectome networks for the human brain, Scientific reports 6, 1 (2016).
  * Barabási and Albert [1999] A.-L. Barabási and R. Albert, Emergence of scaling in random networks, science 286, 509 (1999).
  * Goh  _et al._ [2001] K.-I. Goh, B. Kahng, and D. Kim, Universal behavior of load distribution in scale-free networks, Physical review letters 87, 278701 (2001).
  * Chung and Lu [2002] F. Chung and L. Lu, Connected components in random graphs with given expected degree sequences, Annals of combinatorics 6, 125 (2002).
  * Ichinomiya [2004] T. Ichinomiya, Frequency synchronization in a random oscillator network, Physical Review E 70, 026116 (2004).
  * Peron and Rodrigues [2012] T. K. D. Peron and F. A. Rodrigues, Determination of the critical coupling of explosive synchronization transitions in scale-free networks by mean-field approximations, Physical Review E 86, 056108 (2012).
  * Dorogovtsev  _et al._ [2000] S. N. Dorogovtsev, J. F. F. Mendes, and A. N. Samukhin, Structure of growing networks with preferential linking, Physical review letters 85, 4633 (2000).
  * Krapivsky  _et al._ [2000] P. L. Krapivsky, S. Redner, and F. Leyvraz, Connectivity of growing random networks, Physical review letters 85, 4629 (2000).
  * Mishra  _et al._ [2018] A. Mishra, S. Saha, M. Vigneshwaran, P. Pal, T. Kapitaniak, and S. K. Dana, Dragon-king-like extreme events in coupled bursting neurons, Physical Review E 97, 062311 (2018).
  * Hramov  _et al._ [2006a] A. E. Hramov, A. A. Koronovskii, M. K. Kurovskaya, and S. Boccaletti, Ring intermittency in coupled chaotic oscillators at the boundary of phase synchronization, Physical review letters 97, 114101 (2006a).
  * Hramov  _et al._ [2006b] A. Hramov, A. A. Koronovskii, I. Midzyanovskaya, E. Sitnikova, and C. Van Rijn, On-off intermittency in time series of spontaneous paroxysmal activity in rats with genetic absence epilepsy, CHAOS: An Interdisciplinary Journal of Nonlinear Science 16, 043111 (2006b).
  * Sitnikova  _et al._ [2012] E. Sitnikova, A. E. Hramov, V. V. Grubov, A. A. Ovchinnkov, and A. A. Koronovsky, On–off intermittency of thalamo-cortical oscillations in the electroencephalogram of rats with genetic predisposition to absence epilepsy, Brain research 1436, 147 (2012).
  * Koronovskii  _et al._ [2016] A. A. Koronovskii, A. E. Hramov, V. V. Grubov, O. I. Moskalenko, E. Sitnikova, and A. N. Pavlov, Coexistence of intermittencies in the neuronal network of the epileptic brain, Physical Review E 93, 032220 (2016).
  * Frolov  _et al._ [2019] N. S. Frolov, V. V. Grubov, V. A. Maksimenko, A. Lüttjohann, V. V. Makarov, A. N. Pavlov, E. Sitnikova, A. N. Pisarchik, J. Kurths, and A. E. Hramov, Statistical properties and predictability of extreme epileptic events, Scientific reports 9, 1 (2019).
  * Barabási [2013] A.-L. Barabási, Network science, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 371, 20120375 (2013).
  * Fasino  _et al._ [2021] D. Fasino, A. Tonetto, and F. Tudisco, Generating large scale-free networks with the chung–lu random graph model, Networks 78, 174 (2021).