# Propensity matrix method for age dependent stochastic infectious disease
models111This work was done in the framework of the Hungarian National
Development, Research, and Innovation (NKFIH) Fund 2020-2.1.1-ED-2020-00003.

P. BOLDOG  N. BOGYA222Work partially supported by grant NKFIH-1279-2/2020 of
the Ministry for Innovation and Technology, Hungary and ZS. VIZI333Work
partially supported by grant Functional Differential Equations in Mathematical
Epidemiology National Research, Development and Innovation Office NKFI FK
124016 University of Szeged, Bolyai Institute,
Aradi vértanúk tere 1,
6720, Szeged, Hungary
E-mail<EMAIL_ADDRESS>
###### Abstract

Mathematical modeling is one of the key factors of the effective control of
newly found infectious diseases, such as COVID-19. Our knowledge about the
parameters and the course of the infection is highly limited in the beginning
of the epidemic, hence computer implementation of the models have to be quick
and flexible. The propensity matrix - update graph method we discuss in this
paper serves as a convenient approach to efficiently implement age structured
stochastic epidemic models. The code base we implemented for our forecasting
work is also included in the attached GitHub repository[11].

##  1 Introduction

Quick and effective response to an emerging infection requires a large amount
of information about the course of the epidemic (such as infection and
recovery rate, latency period, etc.) and state of the art mathematical models.
As it has been observed during the COVID-19 pandemic, the information at our
disposal is highly limited, especially at the very beginning of the disease
outbreak – when we have to act quickly. Thus, our models, and therefore their
computer implementations have to be modified drastically on a daily basis.

Experts from the field of mathematical epidemiology, such as Rost et al.[1]
and Barbarossa[2], usually apply deterministic models to predict the spread of
the disease as a first approximation. Nevertheless, realistic models (with
several classes of infection, age structure or spatial patches) tend to get
analytically intractable. Furthermore, deterministic models do not produce
information about several important aspects of the epidemic, such as the
variance of the state variables or the probability of certain events (like
extinction) that are particularly important at the beginning of the outbreak.
Stochastic modeling approach offers feasible alternatives for tackling the
above mentioned problems. Although the governing equation (the so called
stochastic master equation) is often mathematically intractable, numerical
simulations of the corresponding Markov process and the average of the
generated time series provides valuable information of several aspects of the
epidemic.

Gillespie’s stochastic simulation algorithm[4] (SSA) was originally designed
to produce exact realisations of the stochastic master equation in case of
coupled chemical reaction systems. Since then, the scope of the SSA has been
extended by extensive research on this field[7] in the last 40 years to
stochastically simulate other chemical or biological phenomena such as
diffusion processes, migration of cells and animals, or disease spread.
Gillespie’s algorithm yields a convenient method of simulation to gather
information of an epidemic in the stochastic approach.

During the COVID-19 pandemic, one important lesson was that, quick
implementation of complicated stochastic epidemic models with several age
groups is just as crucial as running time or resource efficacy. The aim of the
present paper is to provide a framework that enables researchers to quickly
build and modify stochastic epidemic models with age structure.

The structure of the paper is as follows: in the remaining subsections of the
Introduction, we introduce Gillespie’s first reaction method and provide a
Python code for the case of the SIR model. In Section 2, we introduce our
propensity matrix approach and extend our model implementation with age
structure. Then, we demonstrate the flexibility of our approach by
generalizing the model with other state variables to investigate an
epidemiologically more feasible age-stratified model containing the E
(exposed) and the D (dead) classes additionally. In Section 3, we introduce
the update graph to make the algorithm much more efficient in terms of
computation. Then in Section 4, we demonstrate the strength of stochastic
epidemic models via some experiments with a hypothetical population. Finally,
we summarize and discuss our results in Section 5.

###  1.1 Gillespie’s stochastic simulation algorithms

Considering a group of chemical species ($S_{1},S_{2},\ldots,S_{n}$) reacting
in a coupled system with reaction channels ($R_{1},R_{2},\ldots,R_{m}$) in a
well mixed environment, Gillespie’s stochastic simulation algorithms are based
on answering two questions[4]:

What kind of reaction ($R_{1},R_{2},\ldots,R_{m}$) will the next reaction be?

In what time ($\tau$) will it occur? Gillespie formulated two distinct, but
mathematically identical versions of the SSA: the so called ’direct method’
(DM) and the ’first-reaction method’ (FRM). The main difference is in the way
of deciding which reaction channel to fire. The first-reaction method
generates $m$ random time values $(\tau_{1},\ldots,\tau_{m})$ for the $m$
possible reactions from the corresponding exponential distributions and
selects the reaction channel with the least time to fire. Alternatively, the
direct method only requires two random numbers (one from a uniform
distribution and one from an exponential distribution). Considering the fact
that if there are more than 2 possible reactions then the FRM requires more
computation and memory than the DM, moreover drawing $m$ random numbers from
exponential distributions requires calculating $m$ logarithms, thus the
computational cost of the first-reaction method is even greater[3]. As we are
about to work with robust systems of reactions, we only review and make use of
the direct method in the following.

#### The direct method algorithm

The direct method SSA creates stochastic realizations of the corresponding
Markov chain that is continuous in time and discrete in state variables.
Technically, starting from the vector of initial values
$(X_{1}(0),\dots,X_{n}(0))\in\mathbb{N}_{0}^{n}$, the time series of the state
variables $(X_{1}(t),\dots,X_{n}(t))\in\mathbb{N}_{0}^{n}$ is generated by
constantly updating the state variables in properly generated subsequent times
($\tau$) according to the answers to question (1.1) and (1.1).

To this end, the reaction probability density function, $P(\tau,\mu)d\tau$, is
defined[3]. That is, at time $t$, the next reaction in the reaction chamber
will occur in the differential time interval $(t+\tau,t+\tau+d\tau)$ with the
probability $P(\tau,\mu)d\tau$, and will be an $R_{\mu}$ reaction
($\mu\in{1,\ldots,m}$).

It can be shown that[3], with the procedure called conditioning, the two
variable density function $P(\tau,\mu)$ can be written as the product of two
one-variable probability density functions:

$P(\tau,\mu)=P_{1}(\tau)\cdot P_{2}(\mu|\tau).$ (1)

In particular, from the paper of Gillespie[3], it turns out to be

$P(\tau,\mu)=a_{\mu}\cdot e^{-\tau\cdot a},$

where the symbol $a_{\mu}$ stands for the so called propensity function that
characterizes reaction $R_{\mu}$, and may depend on the quality, and actual
number of the reactants or the environment, etc. For convenience, we use the
notation $a=\sum_{\mu=1}^{m}a_{\mu}$ for the sum of the propensity functions.

Let $P_{1}(\tau)$ be the probability that the next reaction will occur between
times $t+\tau$ and $t+\tau+d\tau$, independent of which reaction it might be.
Similarly, $P_{2}(\mu|\tau)$ is the probability that the next reaction will be
the $R_{\mu}$ reaction, given that the next reaction occurs at $t+\tau$. The
probability $P_{1}(\tau)$ is obtained by summing the reaction probability
density function over all possible $\mu$ values:

$P_{1}(\tau)=\sum_{\mu=1}^{m}P(\tau,\mu)=a\cdot e^{-\tau\cdot a}.$ (2)

Substituting $P_{1}(\tau)$ into (1) and solving for $P_{2}(\mu|\tau)$ we
obtain

$P_{2}(\mu|\tau)=\frac{a_{\mu}}{a}.$ (3)

It is clear that $\int_{0}^{\infty}P_{1}(\tau)d\tau=\int_{0}^{\infty}a\cdot
e^{-\tau\cdot a}=1$. Moreover,
$\sum_{\mu=1}^{m}P_{2}(\mu|\tau)=\sum_{\mu=1}^{m}\frac{a_{\mu}}{a}=1.$

With these notations and probability distributions, the direct method SSA can
be described as follows.

  1. 1.

Initialization:

  2. 2.

set $t\leftarrow 0$,

  3. 3.

set initial values $(X_{1}(0),\dots,X_{n}(0))\in\mathbb{N}_{0}^{n}$,

  4. 4.

prescribe halting conditions $C_{H}$.

  5. 5.

Calculate propensity functions $a_{\mu}$ for all $\mu\in{1,2,\ldots,m}$.

  6. 6.

Decide when the next reaction will occur: choose $\tau$ according to eq. (2),

  7. 7.

choose $r_{1}$ from $(0,1)$ with a uniform distribution,

  8. 8.

obtain $\tau=(1/a)\ln(1/r_{1})$.

  9. 9.

Decide which reaction occurs: choose $\mu$ according to eq. (3),

  10. 10.

choose $r_{2}$ from $(0,1)$ with a uniform distribution,

  11. 11.

take $\mu$ to be the integer for which

$\sum_{j=1}^{\mu-1}a_{j}<r_{2}a\leq\sum_{j=1}^{\mu}a_{j}$

.

  12. 12.

Update time and state variables:

  13. 13.

calculate change $(\Delta X_{1},\dots,\Delta X_{n})\in\mathbb{Z}^{n}$ in
number of reactants according to $\mu$,

  14. 14.

change the number of molecules according to

$(\\!X_{1}\\!(t+\tau),\\!\dots\\!,X_{n}\\!(t+\tau)\\!)\\!=\\!\\!(\\!X_{1}(t)\\!,\\!\dots\\!,\\!X_{n}(t)\\!)\\!+\\!(\\!\Delta\\!X_{1}\\!,\\!\dots\\!,\\!\Delta
X_{n}\\!),$

  15. 15.

set $t\leftarrow t+\tau$.

  16. 16.

Halt if $C_{H}=True$ else continue the process with Step (2).

###  1.2 Gillespie’s algorithm for the SIR model

The well known deterministic SIR model is a system of ordinary differential
equations (4), and it is one of the simplest models to capture the spread of
an infection in a population without demography. It separates the population
into three classes: $S(t)$ \- the number of susceptible, $I(t)$ \- the number
of infected, and $R(t)$ \- the number of recovered individuals in the
population at time $t$, the symbol $()^{\prime}$ stands for the usual time
derivative. The model defines two reactions that we will call events from now
on. Susceptible individuals get infected with rate $\beta$ by making effective
contact with infected individuals in the population. Infected individuals
recover at a rate $\gamma$, and the per capita contact between $S$ and $I$ is
$\frac{SI}{N}$. We consider a constant population size $N=S(t)+I(t)+R(t)$ and
the following rate of change in the classes:

$\begin{cases}S^{\prime}&=-\beta\frac{SI}{N}\\\
I^{\prime}&=\beta\frac{SI}{N}-\gamma I\\\ R^{\prime}&=\gamma I.\end{cases}$
(4)

In this section, by following the work of Gillespie[4] and Allen[5], we apply
the direct method SSA to the SIR model (4). The code we implemented can be
found in our GitHub repository[11]. It is important to note that, as we
prepare to handle robust models in a flexible way, we will not follow the
common practice of reducing the state variables with decoupling (4) by using
the fact that $N$ is constant in time.

In the stochastic SIR model, state variables
$S(t),I(t),R(t)\in\\{0,1,\ldots,N\\}$, and $t\in[0,\infty)$. In Table (1.2),
we summarise the possible events with the corresponding transitions, change in
state variables, and we also define the belonging propensity functions and
probabilities. The intervent time $\tau$ is generated from the exponential
distribution (2) with parameter $1/a$, where $a=\beta\frac{SI}{N}+\gamma I$
and the probability distribution of the possible events is $p_{1}=a_{1}/a$ and
$p_{2}=a_{2}/a$ (with propensities $a_{1}=\beta\frac{SI}{N},a_{2}=\gamma I$).

The table shows the possible events in the SSA in case of the SIR epidemic
model.  Event Transition Change ($\Delta S,\Delta I,\Delta R$) Propensity
Probability Infection $S\rightarrow I$ $(-1,+1,0)$ $a_{1}=\beta\frac{SI}{N}$
$p_{1}=\frac{a_{1}}{a}$ Recovery $I\rightarrow R$ $(0,-1,+1)$ $a_{2}=\gamma I$
$p_{2}=\frac{a_{2}}{a}$

##  2 The propensity matrix method for age dependent stochastic epidemic
models

In order to extend the SSA algorithm to handle the required age structure, we
make the following definitions and assumptions:

  * •

assume a population of $N$ individuals,

  * •

the population is divided into $l$ classes ($l=3$ in the SIR model) according
to the characteristics of the epidemic: $X_{1},\dots,X_{l}$,

  * •

the population is stratified into $n$ age groups,

  * •

this partition of the population leads us to the $n\cdot l$ compartments (that
we will call state variables as well) at time $t$:
$X_{1}^{1}(t),X_{2}^{1}(t),\dots,X_{l-1}^{n}(t),X_{l}^{n}(t)$,

  * •

an event (or reaction) during the simulation corresponds to a transition of an
individual from one compartment to another in the same age group,

  * •

we consider the epidemic in such a short time scale when aging is not
significant, thus we assume there is no transition between age groups,

  * •

we assume that there is no demography (birth or natural death) during this
short time – if death is incorporated in the model, it only occurs due to the
infection,

  * •

the same transitions occur in every age group – possibly with zero rate.

Now, we define the propensity matrices as follows: for age group $i$, we
define propensity matrix $M^{i}$ with index $i$ in the following way: element
of the propensity matrix $M^{i}(j,k)=a_{j,k}^{i}$ represents the propensity
function corresponding to the transition of an individual from compartment
$X_{j}\rightarrow X_{k}$ at age group $i$.

During the simulation, we calculate the propensity matrix for every age group.
Since, we observe the same transitions in every age group, we can outsource
the repeating calculations into $n$ cycles, obtaining a cleaner and more
compact code.

Then we choose the intervent time $\tau$ from distribution (2), where $a$ is
the sum of all propensities over all propensity matrices:

$a=\sum_{g=1}^{n}\sum_{f=1}^{l}\sum_{h=1}^{l}a_{f,h}^{g},$ (5)

and choose the next event from distribution

$P_{2}\left((\text{transition }j\rightarrow k\text{ at age group
}i)|\tau\right)=p_{j,k}^{i}=\frac{a_{j,k}^{i}}{a}.$ (6)

Basically, one can obtain a transition from distribution (6) by choosing a
random real $r\in(0,1)$ with uniform distribution and summing the propensities
of all transitions over all age groups, then choosing transition
$X_{j}\rightarrow X_{k}$ at age group $i$ whenever

$\sum_{g=1}^{n}\sum_{f=1}^{l}\sum_{h=1}^{l}a_{f,h}^{g}<r\cdot
a\leq\sum_{g=1}^{i}\sum_{f=1}^{j}\sum_{h=1}^{k}a_{f,h}^{g}.$

Eventually, we have to update time: $t\leftarrow t+\tau$ and the state
variables with the change vector $(\Delta X_{1}^{1},\ldots,\Delta X_{l}^{n})$,
where only $\Delta X_{j}^{i}=-1$ and $\Delta X_{k}^{i}=+1$ are nonzero
elements:

$(X_{1}^{1}(t+\tau),\ldots,X_{l}^{n}(t+\tau))=(X_{1}^{1}(t),\ldots,X_{l}^{n}(t))+(\Delta
X_{1}^{1},\ldots,\Delta X_{l}^{n}).$

With these assumptions and notations we are prepared to present the propensity
matrix algorithm that derives the time evolution of the state variables.

### Algorithm propensity matrix method

  1. 0.

Fix the order of the state variables: $(X_{1}^{1},\ldots,X_{l}^{n})$.

  2. 1.

Initialize: Set $t\leftarrow 0$, initial values
$(X_{1}^{1}(0),\ldots,X_{l}^{n}(0))$, and halting condition $C_{H}$.

  3. 2.

Calculate propensity matrices $M_{i}$ for all age groups
$i\in\\{1,\ldots,n\\}$.

  4. 3.

Choose intervent time $\tau$ from $P(\tau)=a\cdot e^{-\tau\cdot a}$ where $a$
comes from (5),

  5. 4.

Choose the next reaction from the distribution
$p_{j,k}^{i}=\frac{a_{j,k}^{i}}{a}$ (6),

  6. 5.

Update state variables and time.

  7. 6.

Halt if $C_{H}=True$ else continue with step 2.

Propensity matrices of epidemic models are sparse, since, from one
compartment, there are usually only one or very few transitions to other
compartments. We also emphasise that this data structure provides a convenient
way to easily modify our model by including new transitions (like waning
immunity), or deleting transitions between age groups. Also, by extending the
number of state variables, we can further introduce new classes (like exposed,
latent, hospitalized, dead or recovered). In the following, we present the
propensity matrices of two epidemic models.

###  2.1 The age structured SIR model with waning immunity

When we extend model (4) with age structure we have to consider the effect of
every infected age group on the susceptible population of age group $i$,
furthermore to include waning immunity to (4) we consider the
$R_{i}\rightarrow S_{i}$ transitions with age dependent rate $\omega_{i}$, in
the general case. Naturally, with $\omega_{i}=0$ for $i\in\\{1,\ldots,n\\}$ we
get back to the case where waning immunity is not applied in the process.
Hence, with age dependent transmission rates $\beta_{i}$, recovery rates
$\gamma_{i}$ and waning rates $\omega_{i}$ we obtain the system of ODEs for
age group $i$:

$\begin{cases}\displaystyle
S_{i}{{}^{\prime}}=-\frac{1}{N}S_{i}\cdot\sum_{f=1}^{n}\beta_{f}I_{f}+\omega_{i}S_{i}\\\
\displaystyle
I_{i}{{}^{\prime}}=\frac{1}{N}S_{i}\cdot\sum_{f=1}^{n}\beta_{f}I_{f}-\gamma
I_{i}\\\ \displaystyle R_{i}{{}^{\prime}}=\gamma
I_{i}-\omega_{i}S_{i}\end{cases}$ (7)

In the stochastic simulation formulation, let the order of state variables be
$(X_{1}^{i}=S_{i},X_{2}^{i}=I_{i},X_{3}^{i}=R_{i})$, thus let and with this
notation, $a_{jk}^{i}$ is the propensity of transition from compartment
$X_{j}$ to $X_{k}$ at age group $i$. Thus, the possible transitions in age
group $i$ (and the belonging propensities) in this model are infection
($a_{SI}^{i}$), waning immunity ($a_{SR}^{i}$) and recovery ($a_{IR}^{i}$).
The propensity matrix and the corresponding propensities can be seen in Figure
1, where ’$\bullet$’ stand for a zero element of the matrix. Python code for
this model is included in the repository[11].

$\begin{cases}\displaystyle
a_{SI}^{i}=\frac{1}{N}S_{i}\cdot\sum_{f=1}^{n}\beta_{f}I_{f}\\\ \displaystyle
a_{SR}^{i}=\alpha_{i}S_{i}\\\ \displaystyle
a_{IR}^{i}=\gamma_{i}I_{i}\end{cases}$

Figure 1: Propensity matrix and the corresponding propensities of the SIR
model

###  2.2 SEIRD model with age structure and waning immunity

One strength of the propensity matrix technique is the flexibility in terms of
computer implementation. For instance several disease models require an
exposed class E as right after infection susceptibles usually do not show
symptoms during the latency period of length $1/\epsilon$. Thus with rate
$\epsilon$ people move from class $E_{i}$ to class $I_{i}$. The corresponding
propensity is $a_{EI}^{i}=\epsilon E_{i}$. We can further add class D that
counts the death cases caused by the infection - with age dependent mortality
rate $p_{i}\gamma$ people move from class $I_{i}$ to $D_{i}$ with
corresponding propensity, $a_{ID}^{i}=p_{i}\gamma I_{i}$. We also assume that
recovered people loose immunity against the disease with rate $\omega$, thus,
the belonging propensity is $a^{i}_{RS}=\omega R_{i}$.

Let the order of state variables at age group $i$ be
$X_{1}^{i}=S_{i},X_{2}^{i}=E_{i},X_{3}^{i}=I_{i},X_{4}^{i}=R_{i},X_{5}^{i}=D_{i}$
and $a_{jk}^{i}$ stands for the propensity of transition from compartment $j$
to $k$ at age group $i$. For flowchart of the process and the governing ODE
model with the propensity matrix and the remaining propensities see Figure 2.

$\begin{cases}S_{i}{{}^{\prime}}=-\frac{1}{N_{i}}S_{i}\cdot\sum_{f=1}^{n}\beta_{f}I_{f}+\omega
R_{i}\\\
E_{i}{{}^{\prime}}=\frac{1}{N_{i}}S_{i}\cdot\sum_{f=1}^{n}\beta_{f}I_{f}-\epsilon
E_{i}\\\ I_{i}{{}^{\prime}}=\epsilon E_{i}-\gamma I_{i}\\\
R_{i}{{}^{\prime}}=(1-p_{i})\gamma I_{i}-\omega R_{i}\\\
D_{i}{{}^{\prime}}=p_{i}\gamma I_{i}\end{cases}$

Figure 2: The SEIRD model with age structure and waning immunity.
Possible transitions with the belonging propensities are: infection
($a_{SE}^{i}=\frac{1}{N_{i}}S_{i}\cdot\sum_{f=1}^{n}\beta_{f}I_{f}$), going
infectious from exposed ($a_{EI}^{i}=\epsilon E_{i}$), recovery
($a_{IR}^{i}=(1-p_{i})\gamma I_{i}$), death due to infection
($a_{ID}^{i}=p_{i}\gamma I_{i}$) and loss of immunity ($a_{RS}^{i}=\omega
R_{i}$).

##  3 Updating the propensities

Gillespie[4] suggests that ”…it is necessary to recalculate only those
quantities $a_{\nu}$, corresponding to reactions $R_{\nu}$, whose reactant
population levels were just altered…” in the reaction selection step.

Probably the most well known attempt to get around this concern belongs to
Gibson and Bruck[8]. Their Next Reaction Method (NRM) may be regarded as an
extension of Gillespie’s original first reaction method. They define the so
called reaction dependency graph that contains information about which
propensity function ($a_{\nu}$) needs to be updated according to the chosen
reaction in the selection step. However, in the detailed comparison of the DM,
FRM and NRM by Cao et al.[6] it is found that ”even with the best data
structure, the NRM is less efficient than the DM except for a very specialized
class of problems”. As Scvehm[9] points out, this is mostly due to the fact
that in case of the Gibson-Bruck method ”the simulator engine spends most of
its execution time for maintaining the priority queue of the tentative
reaction times”.

In this section we provide a dependency graph like method to make use of the
fact that upon a transition event mostly only a small number of propensity
functions have to be recalculated. For instance in case of the age structured
SEIRD model in Section 2.2 the waning immunity event at age group $i$ only
changes state variable $R_{i}$ and $S_{i}$ hence it only influences propensity
function $a^{i}_{RS}$ and $a^{i}_{SE}$. All other propensities in age group
$i$ and every other propensities in the other age groups remain unchanged,
that is $3+(n-1)\cdot 5=5n-2$ number of events in case of $n$ age groups.

In our methodology, in every iteration step after choosing the transition from
distribution (6), we selectively update the propensity values based on this
transition. We do this with the help of a suitable data structure that is a
directed bipartite graph, called the Update Graph (UG) and is defined the
following way: let $C$ be the set of nodes associated with the classes (cf.
Sec. 2) and let $T$ be the other set of nodes representing the possible
transitions between the classes. The edge $c\rightarrow t_{X_{j},X_{k}}$
(where $t_{X_{j},X_{k}}$ is the transition from $X_{j}$ to $X_{k}$) exists
only if the class of a modified state variable $c\in C$ updates the propensity
of transition $t_{X_{j},X_{k}}\in T$. Also edge $t_{X_{j},X_{k}}\rightarrow c$
exists only if $c=X_{j}$ or $c=X_{k}$. For convenience we note that from any
node $t_{X_{j},X_{k}}\in T$ there are always two edges pointing out, one to
class $X_{j}$ and the other one to class $X_{k}$.

However from class $c\in C$ there may be several edges pointing to different
transitions, for instance in the $SEIRD$ model infection, recovery and death
are all depending on the number of $I$ individuals in a given age group. We
also point out that (in case of the $SEIRD$ model) $S_{i}$ individuals can be
infected not only by individuals from compartment $I_{i}$ but from all
compartment $I_{g}$ at all age groups, $g\in\\{1,\dots,n\\}$. Therefore,
change in any infectious compartment affects the propensity for transitions in
all age groups. Thus upon an infection event all propensities
$a_{SE}^{g},g\in\\{1,\dots,n\\}$ have to be updated. In more complicated
models there may be several infectious classes (such as latent or in case of
ebola models even dead individuals may cause infection). To this end, in
general, we flag (*) the classes $c^{*}$ that influence propensities across
all age groups and also flag transitions $t^{*}$ that have to be updated for
all age group whenever a compartment from a flagged class changed. We handle
this issue in the computer implementation.

For instance the two part of the UG of the SEIRD model with waning immunity is
shown on Fig.3. The left part of the figure shows edges
$t_{X_{j}X_{k}}\rightarrow c$ and the right part shows edges $c\rightarrow
t_{X_{j}X_{k}}$. Class $I^{*}$ and transition $S\xrightarrow[]{*}E$ is flagged
thus whenever an $I^{i}$ compartment changes all propensities $a_{SE}^{i}$
need to be updated.

Figure 3: Update graph of the SEIRD model The figure shows the bipartite
update graph that contains information about the classes and propensity
functions to be updated upon event selection.

During the simulation, after the event selection step we end up with the
coordinate triple $(i,j,k)$ that selects transition $t_{X_{j},X_{k}}$ at age
group $i$. By using the $t_{X_{j},X_{k}}\rightarrow c$ edges of the graph
(left figure) we update the state variables according to
$X_{j}^{i}(t+\tau)=X_{j}^{i}(t)-1$, and $X_{k}^{i}(t+\tau)=X_{k}^{i}(t)+1$ and
all other $X_{f}^{i}(t+\tau)=X_{f}^{i}(t)$ remain unchanged. Next we update
the propensity $a_{jk}^{i}$ according to the $c\rightarrow t_{X_{j},X_{k}}$
edges. If we encounter and edge that points from a flagged class to a flagged
transition then we updated all $a_{jk}^{i}$ for $i\in\\{1,\dots,n\\}$.

###  3.1 Propensity matrix method extended with the update graph

  1. 0.

Initialisation step:

     * •

Fix the order of the state variables: $(X_{1}^{1},\ldots,X_{l}^{n})$

     * •

Construct the Update Graph according to Sec. 3

     * •

Set $t\leftarrow 0$, initial values $(X_{1}^{1}(0),\ldots,X_{l}^{i}(0))$, and
halting condition $C_{H}$.

     * •

Calculate propensity matrices $M_{i}$ for all age groups
$i\in\\{1,\ldots,n\\}$

  2. 1.

Selection step:

     * •

Choose intervent time from $P(\tau)=a\cdot e^{-\tau\cdot a}$ where $a$ comes
from (5),

     * •

Choose the next reaction $t_{X_{j},X_{k}}$ at age group $i$ from the
distribution $p_{j,k}^{i}=\frac{a_{j,k}^{i}}{a}$ (6),

  3. 2.

Update step:

     * •

Time: $t\leftarrow t+\tau$

     * •

Update state variables: by using the $t_{X_{j},X_{k}}\rightarrow c$ edges of
the graph update $X_{j}^{i}(t+\tau)=X_{j}^{i}(t)-1$, and
$X_{k}^{i}(t+\tau)=X_{k}^{i}(t)+1$ and all other
$X_{f}^{i}(t+\tau)=X_{f}^{i}(t)$ remain unchanged,

     * •

Update propensities: update the propensity $a_{jk}^{i}$ according to the
$c\rightarrow t_{X_{j},X_{k}}$ edges. If we encounter and edge that points
from a flagged class to a flagged transition then we updated all $a_{jk}^{i}$
for $i\in\\{1,\dots,n\\}$.

  4. 3.

Halt if $C_{H}=True$ else continue with step 2.

##  4 Experiments

In this section, we demonstrate the advantages of the stochastic approach and
the flexibility of the propensity matrix - update graph method by
experimenting with some real life problems.

During the following simulations, we consider a population of 200,000
individuals divided into three age groups (0-14, 15-59 and 60+), with the
aggregated contact matrix from Prem et al.[10] and with the population
distribution of Hungary from the Hungarian Central Statistical Office (KSH)
age-stratification[12] (cf. Figure 4). We symmetrised the contact matrix
according to Sec. 2.3.3. in Röst et al.[1].

Figure 4: Population related parameters of the demonstrated experiments
Agregated contact matrix (left) and age stratification (right) of the
investigated model population. Both data collection corresponds to the
Hungarian society.

Using the stochastic SEIRD model with immunity waning from Sec.2.2 we
investigate the course of a COVID-19 epidemic with aggregated parameters from
Rost et al.[1] (cf. Table 4). Every simulation starts with $I_{1}(0)=10$
infected individuals in age group 1 at $t=0$ – all other individuals
considered to be susceptible according to the age partition
($S_{1}(0)=29119,S_{2}(0)=117637,S_{3}(0)=53234$). The governing propensity
functions can be found in Sec. 2.2.

Parameters of the demonstrated epidemic  Parameter Value/age-dependent vector
Incubation period ($\epsilon^{-1}$) $5.2$ (days) Infectious period
($\gamma^{-1}$) $5.0$ (days) Time spent immunized ($\omega^{-1}$) $180$ (days)
Infection rate ($\beta$) $0.05$ Probabilty of death ($p_{i}$) ($0.0000451$,
$0.00117$, $0.0281$)

The left part of figure 5 shows the time series of the age aggregated state
variables $\sum_{i=1}^{3}S_{i},\sum_{i=1}^{3}E_{i}$, etc. from a single
simulation outcome. We can observe that during the ”first wave” the exposed
class peaks before the infected class, and class $S$ and $R$ shows oscillation
due to loss of immunity. The right part of Figure 5 shows the time series of
the three age groups of class $S$ from the same simulation.

Figure 5: Stochastic realisation of the SEIRD model with immunity waning
Time series of the age agregated state variables (left), and the three
susceptible age groups (right) from the same realisation of the stochastic
SEIRD model. For the parameters cf. Fig.4 and Tab.4.

On Figure 6 we focus on the peak size of the epidemic and show the effect of
contact reduction in the early stage of the disease spread (on day 45 in this
case). By $10\%,20\%$ and $30\%$ of contact reduction we mean that we multiply
every element of the contact matrix with $0.9,0.8$ and $0.7$ respectively. The
simulation was halted 10 days after the peak. We can notice that a uniform
contact reduction of $30\%$ in every age group may decrease the peak size by
$1/3$, and in the meanwhile it also delays the time of the peak.

Figure 6: Effect of contact reduction in the early stage of the epidemic
In this simulation we applied a $10\%,20\%$ and $30\%$ uniform contact
reduction in the total contact structure of the population. The simulation
stops 10 days after the peak.

Montecarlo simulations, such as any version of Gillespie’s SSA can only serve
information about the process if we run several stochastic realisations and
investigate the statistics of them. Thus it is always important to know how
many simulations we need for plausible conclusions.

In the following experiments we run a large number of simulations the
following way: we obtain 100 simulations (for having enough data for
statistical analysis) and for the daily sampled time series of every single
compartment $X^{i}_{j}(t)$ we evaluate the mean $\mu_{X_{j}^{i}}(t)$ and
variance $\sigma_{X_{j}^{t}}(t)$. We also calculate the mean value and
variance for the time of the peak $(\mu_{t},\sigma_{t})$ and for the size of
the peak $(\mu_{s},\sigma_{s})$. Then after each simulation we check whether
95% of the collected first peak sizes are in the confidence interval
$[\mu_{s}-2\sigma_{s},\mu_{s}+2\sigma_{s}]$ as well as in
$[\mu_{t}-2\sigma_{t},\mu_{t}+2\sigma_{t}]$ and stop the routine, if this
condition is fulfilled.

On both part of Fig. 7, with solid red curve, we plotted the mean of the daily
sampled time aggregated infected compartments $\sum_{i=1}^{3}\mu_{I_{i}}(t)$.
The blue band on the left plot shows the minimal and maximal values we
calculated for each days from the ensemble of the independent realisations of
the process as well. Let us remark that the band at the top (around the peak)
is almost flat, meaning that the different simulations produced similar peak
sizes but different peak times. On the right part of Fig. 7 the highlighted
rectangle shows the $2\sigma_{t}$ wide and $2\sigma_{s}$ high confidence
interval around the maximum of the mean. Thus we may conclude that the peak of
the epidemic occurs between day 63 and 79 and it is expected to be between
approximately 27,900 and 28,600 with 95% probability.

Figure 7: Time series of the daily sampled aggregated infected classes
The red curve shows the daily sampled $\sum_{i=1}^{3}\mu_{I_{i}}(t)$ in case
of 110 simulations around the peak of the first wave. The blue band on the
left shows the daily sampled minimal and maximal values of the stochastic
ensemble, the blue rectangle on the right represents the $2\sigma_{t}$ wide
and $2\sigma_{s}$ high confidence interval around the maximum of the mean.

##  5 Summary and discussion

In this paper we introduced our approach that we use in our forecasting work
to efficiently implement age structured stochastic epidemic models. The core
idea, the so called propensity matrix method, serves a data structure to
handle the propensities and select the upcoming reaction/event to execute in a
convenient way, in case of robust epidemic models. In Section 2 we introduced
this method and the algorithm that obtains the time evolution of the state
variables. In Section 3 we further improved the algorithm by introducing the
so called update graph that helps us to speed up the algorithm by updating
only the minimal number of propensities that are required to be updated due to
the change in the state variables. Finally in Section 4 we showed some real
life experiments to demonstrate the strength of the stochastic approach and
the flexibility of our method.

The speed of the simulation is usually a limiting factor of stochastic
simulations. Generations of single trajectorie requires several independent
runs and therefore are typically expensive in terms of computational time.
From coding perspective interpreted languages like Python may have low
performance. The simulation time may be reduced by several magnitudes using
low-level languages like C or C++. The propensity matrix algorithm can be
easily applied to such languages. We provide our sample code in Python to
emphasize the algorithm and help readability as much as possible.

However, from the modeling perspective simulation of whole trajectories during
years, from the outbreak untill the very end of the epidemic, is usually
unnecessary in practice. Information provided by stochastic simulations is
particularly valuable i) in the beginning of the epidemic - when only a small
portion of the population is infected, ii) in case of parameter estimation,
iii) when we want to make short-term forecasts, or iv) we want to gain
information of the variance of state variables in a short time scale - such
like in Section 4 where we estimated the number of infected individuals near
the peak of the epidemic. Thus, reducing the scope of the simulation to
shorter time intervals is usually a convenient way of dealing with speed.

Incorporating demography and aging into our models may change the process
drastically in a long term and may lead to much more realistic models.
However, in case of short term simulations it usually leads to unnecessary
complications and a great number of (aging) events that occupies the
simulation engine and it leads to slower simulations. Moreover in case of a
small number of age groups the change between two age groups during a short
term simulation remains negligible. However, if unavoidable, we suggest that
instead of Gillespie’s[3] exact stochastic algorithm consider an approximating
algorithm when aging events occur in discrete time steps, say on every 30th
day of the simulation, and the appropriate portion of the population in
compartment $X_{j}^{i}$ (for $j\in\\{0,\dots,l\\}$) moves to compartment
$X_{j}^{i+1}$ if $i\neq n$ or dies whenever $i=n$. In the meanwhile
appropriate portion of newborns need to be added into the susceptible,
immunised or the infected compartments of age group 1.

## References

  * [1] G. Röst, F. Bartha, N. Bogya, P. Boldog, A. Dénes, T. Ferenci, K. Horváth, A. Juhász, C. Nagy, T. Tekeli, Z. Vizi. and B. Oroszi, Viruses. 12 (2020)
  * [2] M. V. Barbarossa, PLOS ONE. 15, 1-22 (2020,9)
  * [3] D. Gillespie, Journal Of Computational Physics. 22, 403-434 (1976)
  * [4] D. Gillespie, The Journal Of Physical Chemistry. 81 (1977)
  * [5] L. Allen, Infectious Disease Modelling. 2, 128-142 (2017)
  * [6] L. Yang Cao, The Journal Of Chemical Physics. 121, 4059-4067 (2004)
  * [7] D. Gillespie, A. Hellander and L. R. Petzold, The Journal of Chemical Physics 138, 170901 (2013)
  * [8] M Gibson and J. Bruck, The Journal Of Physical Chemistry A. 104 (2000,3)
  * [9] M. Schwehm, Poster in German Conference on Bioinformatics. (2001)
  * [10] K. Prem, PLOS Computational Biology. 13, 1-21 (2017,9)
  * [11] Z. Vizi, GitHub Repository. (2021), https://github.com/zsvizi/propensity-matrix-epidemics
  * [12] Hungarian Central Statistical Office, website. Accessed: 2021-09-20, www.ksh.hu/stadat_files/nep/en/nep0003.html


# Exploring thermal equilibria of the Fermi–Hubbard model
with variational quantum algorithms

Jack Y. Araz  0000-0001-8721-8042<EMAIL_ADDRESS>Jefferson Lab, Newport
News, VA 23606, United States  Michael Spannowsky  0000-0002-8362-0576
<EMAIL_ADDRESS>Institute for Particle Physics Phenomenology,
Durham University, Durham, United Kingdom  Matthew Wingate  0000-0001-6568-988
<EMAIL_ADDRESS>DAMTP, University of Cambridge, Cambridge, CB3 0WA,
United Kingdom

###### Abstract

This study investigates the thermal properties of the repulsive Fermi-Hubbard
model with chemical potential using variational quantum algorithms, crucial in
comprehending particle behaviour within lattices at high temperatures in
condensed matter systems. Conventional computational methods encounter
challenges, especially in managing chemical potential, prompting exploration
into Hamiltonian approaches. Despite the promise of quantum algorithms, their
efficacy is hampered by coherence limitations when simulating extended
imaginary time evolution sequences. To overcome such constraints, this
research focuses on optimizing variational quantum algorithms to probe the
thermal properties of the Fermi-Hubbard model. Physics-inspired circuit
designs are tailored to alleviate coherence constraints, facilitating a more
comprehensive exploration of materials at elevated temperatures. Our study
demonstrates the potential of variational algorithms in simulating the thermal
properties of the Fermi-Hubbard model while acknowledging limitations stemming
from error sources in quantum devices and encountering barren plateaus.

††preprint: IPPP/23/72, JLAB-THY-23-3968

## Introduction

Quantum physics presents profound challenges in unravelling the intricate
dynamics governing particle interactions within condensed matter systems.
Among the pivotal models illuminating such phenomena, the Fermi-Hubbard model
[1, 2] encapsulates correlated particles’ behaviour within a lattice
framework. A crucial facet in comprehending these systems lies in probing
their thermal properties.

However, the computational scrutiny of models encounter substantial hurdles,
especially when traditional methodologies confront facets like the
introduction of a chemical potential [3] or topological terms [4]. Notably,
Monte Carlo algorithms falter in handling such characteristics, compelling the
exploration of Hamiltonian methodologies for their robustness in navigating
these complexities.

Despite the promise inherent in Hamiltonian methods, challenges endure. Due to
their Hamiltonian-based approach, Tensor Networks do not suffer from sign
problems. Employing imaginary time evolution through techniques like Time
Evolving Block Decimation (TEBD) allows for accurate simulation of thermal
states for lattice field theories [5, 6, 7, 8]. However, while potent, TEBD
confronts computational bottlenecks; the escalating computational demand due
to large bond dimensions associated with extensive time steps restricts the
efficacy of such methods in capturing the thermal equilibria of these systems.

To address these challenges, attention has pivoted towards leveraging the
capabilities of quantum computing. Quantum algorithms, immune to specific
computational roadblocks faced by classical algorithms, offer a promising
avenue for exploring complex condensed matter systems. Although it is possible
to implement imaginary time evolution for quantum algorithms to study the
thermal properties of lattice systems without suffering from large bond
dimensions [9], the finite coherence time in quantum devices imposes
constraints on simulating prolonged time evolution sequences, necessitating
alternative strategies.

This study focuses on employing variational quantum algorithms to investigate
the thermal characteristics of the spin-$\tfrac{1}{2}$ Fermi-Hubbard model.
These techniques find their roots in the Variational Quantum Eigensolver (VQE)
[10, 11, 12], which utilizes a trainable unitary gate sequence to approximate
the diagonalization of a given Hamiltonian. Previous studies have applied this
method to analyze the ground state properties of the Fermi-Hubbard model [13,
14]. Furthermore, references [15, 16, 17] have extended VQE techniques to
learn thermal state density matrices using hybrid classical-quantum network
systems. Despite their success, integrating classical and quantum computing
methods creates unavoidable bottlenecks, as evidenced in data analysis
applications [18]. These challenges have prompted exploring and developing
alternative quantum-based algorithms [19, 20, 21, 22].

We employ the quantum Variational Quantum Thermaliser (qVQT) algorithm,
proposed in ref. [20]. Through efficient circuit designs, we aim to surmount
the limitations of quantum coherence constraints. By doing so, we aim to glean
more profound insights into the behaviour of these materials at elevated
temperatures, thereby contributing to the foundational understanding of
quantum condensed matter systems.

Subsequent sections delve into the methodological intricacies deployed in our
study; the model will be introduced in section I, and in section II, we will
present our methodology, delineating the tailored circuit designs. Finally, in
sec. III, we will discuss the insights gleaned from our investigations. This
research endeavours to provide a nuanced understanding of materials at a
quantum scale, potentially paving the way for advancements in material science
and technological innovation.

##  I The repulsive Fermi-Hubbard model with chemical potential

The Fermi-Hubbard model belongs to the category of “tight-binding models”,
representing effective electron-electron interactions by a contact term. In
its second-quantized form, the Fermi-Hubbard Hamiltonian is expressed as:

$\displaystyle\hat{H}$ $\displaystyle=$ $\displaystyle-t\sum_{\langle
i,j\rangle,s}c^{\dagger}_{i,s}c_{j,s}+U\sum_{i}\hat{n}_{i,\uparrow}\hat{n}_{i,\downarrow}-\mu\sum_{i}\hat{n}_{i}\
.$ (1)

Here, $\hat{n}_{i,s}=c^{\dagger}_{i,s}c_{i,s}$ and
$\hat{n}_{i}=\sum_{s}\hat{n}_{i,s}$ represent the number operator. The
parameter $U$ denotes the on-site repulsive coupling constant, while $\mu$
stands for the chemical potential [1, 2]. The notation $\langle i,j\rangle$
denotes nearest-neighbor interaction, and the $s$ index accounts for spin
orientations, where $s\in\\{\uparrow,\downarrow\\}$. The creation
(annihilation) operators, $c_{is}^{\dagger}$ ($c_{is}$), follow
anticommutation relations:

$\displaystyle\left\\{c_{is},c_{js^{\prime}}^{\dagger}\right\\}=\delta_{ij}\delta_{ss^{\prime}}\
,\
\left\\{c_{is},c_{js^{\prime}}\right\\}=\left\\{c_{is}^{\dagger},c_{js^{\prime}}^{\dagger}\right\\}=0\
.$

For $\mu\neq 0$, eq. (1) suffers from the so-called sign problem due to the
complex terms in the partition function, leading MC-based algorithms to become
ill-posed [23, 24, 25].

In order to embed such a Hamiltonian on a quantum device, one needs to go
through a qubitization process via the Jordan-Wigner transformation [26, 27]

$\displaystyle
c^{\dagger}_{j}=\left[\prod_{k=1}^{j-1}\left(-Z_{k}\right)\right]S^{+}_{j}\quad,\quad
S^{\pm}_{j}=\frac{1}{2}\left(X_{j}\pm iY_{j}\right)\ ,$ (2)

where $X,\ Y$ and $Z$ represent Pauli operators. Notice that for nearest-
neighbour interactions, the chain of $Z$ operators vanishes and only reappears
at the periodic boundaries. By applying eq. (2) to eq. (1) one can rewrite the
terms of the Hamiltonian $\hat{H}=H_{0}+H_{1}+H_{2}$ as

$\displaystyle-\frac{H_{0}}{t}$ $\displaystyle=$ $\displaystyle\sum_{\langle
i,j\rangle,s}^{N-1}S^{+}_{i,s}S^{-}_{j,s}$ (3) $\displaystyle+$
$\displaystyle\left[\prod_{i=1}^{N-1}\left(-Z_{i}\right)\right](S^{+}_{N,s}S^{-}_{0,s}+{\rm
h.c.})$ $\displaystyle H_{1}$ $\displaystyle=$
$\displaystyle\frac{U}{4}\sum_{i}^{N}\left(Z_{i,\uparrow}Z_{i,\downarrow}+Z_{i,\uparrow}+Z_{i,\downarrow}+\mathds{1}_{i}\right)\
,$ (4) $\displaystyle H_{2}$ $\displaystyle=$
$\displaystyle-\frac{\mu}{2}\sum_{i,s}\left(Z_{i,s}+\mathds{1}_{i,s}\right)\
,$ (5)

where $N$ stands for the number of sites. Notice that in order to simulate
both spin flavours, one needs to have two sets of $N$ qubits, one for spin-up
and one for spin-down, which increases the number of qubits to $N_{q}=2N$.

##  II Quantum Simulation

In this section, we outline the methodology employed to investigate the
thermal equilibria of the Fermi-Hubbard model at a fixed temperature using
variational quantum algorithms. Our primary focus will be on the method
illustrated in ref. [20]. This exercise represents the thermal mixed state
induced by the Hamiltonian within quantum circuits.

The density matrix of a given Hamiltonian $\hat{H}$ at a fixed temperature can
be written as:

$\displaystyle\rho_{\beta}=\frac{e^{-\beta\hat{H}}}{Z}\quad,\quad
Z=\mathrm{Tr}\left[e^{-\beta\hat{H}}\right]\ ,$ (6)

where $\beta$ is the inverse temperature ($1/T$), and $Z$ is the partition
function. For a variational ansatz, the thermal state at a fixed temperature
can be found by optimising the system with respect to the Helmholtz free
energy $F=E-TS$, where $E$ corresponds to the expectation value of the
Hamiltonian at a fixed temperature
$\langle\hat{H}\rangle_{\beta}=\mathrm{Tr}\left[\rho_{\beta}\hat{H}\right]$,
and entropy, $S$, can be written as
$S=-\mathrm{Tr}\left[\rho_{\beta}\log\rho_{\beta}\right]$. The challenging
part of the quantum simulation arises from the fact that a quantum computer is
a pure-state simulator, and one does not have direct access to the density
matrix within the circuit.

Our approach’s key components involve utilising two distinct quantum circuits,
each designed to address specific aspects of the problem. Our procedure is
summarised in the following three subsections.

###  II.1 Entropy estimation

The first variational quantum circuit (VQC1) determines the probabilities
associated with each pure state within the ensemble under investigation. This
circuit calculates the probability of finding the system in each pure state.
Since it is not possible to extract the density matrix directly to compute the
entropy, one can compute the Shannon entropy after measurement:

$\displaystyle S=-\sum_{i}p_{i}(\theta_{1})\log p_{i}(\theta_{1})\ ,$ (7)

where $p_{i}(\theta_{1})=\left\lVert\langle
0|\hat{U}_{1}^{\dagger}(\theta_{1})|\phi_{i}\rangle\right\rVert^{2}$
corresponds to the probability of the $i$-th pure state, $|\phi_{i}\rangle$,
within the ensemble. Here, $\theta_{1}$ represents the variational parameters
of the variational quantum circuit represented by $\hat{U}_{1}(\theta)$.

###  II.2 Estimation of the expectation value of the Hamiltonian

The second quantum circuit (VQC2) is responsible for computing the expectation
value of the Hamiltonian for each of the pure states identified in the
previous step. The Hamiltonian operator encapsulates the system’s energy
characteristics, and its expectation value reveals the energy associated with
each pure state. The expectation value of the density matrix can be computed
as:

$\displaystyle\langle\hat{H}\rangle_{\beta}=\sum_{i}p_{i}(\theta_{1})\langle\phi_{i}|\hat{U}_{2}(\theta_{2})\hat{H}\hat{U}_{2}^{\dagger}(\theta_{2})|\phi_{i}\rangle,$
(8)

where $\theta_{2}$ represents the variational parameters of the variational
quantum circuit represented by $\hat{U}_{2}(\theta)$.

###  II.3 Objective function and optimization

With the probabilities and energies of pure states at hand, we compute the
system’s total free energy. This is achieved through a combination of the
Shannon entropy, derived from the probabilities obtained from the first
quantum circuit, and the energy obtained from the second quantum circuit. The
free energy is a fundamental thermodynamic quantity that encapsulates the
system’s equilibrium properties.

Having established the total free energy of the system, we employ variational
quantum algorithms to optimize our quantum circuits. The objective is to
minimize the free energy:

$\displaystyle\min_{\theta_{1,2}}F(\theta_{1},\theta_{2})\ ,$

which corresponds to finding the variational parameters that best represent
the system’s thermal equilibrium. The variational parameters, $\theta_{1,2}$,
can be updated with respect to the gradient of the free energy:

$\displaystyle\theta^{\prime}_{2}=f\left(\theta_{2},\frac{\partial\langle\hat{H}\rangle_{\beta}}{\partial\theta_{2}}\right)\quad,\quad\theta^{\prime}_{1}=f\left(\theta_{1},\frac{\partial\langle\hat{H}\rangle_{\beta}}{\partial\theta_{1}}-T\frac{dS}{d\theta_{1}}\right),$

where $f$ represents a function to be employed to update variational
parameters, e.g., gradient descent. This optimization process enables us to
refine our quantum states, ultimately leading to a more accurate
representation of the system’s behaviour.

##  III Results

We tested the algorithm described in the previous section on a four-site
system, which required eight qubits to describe the Fermi-Hubbard model, with
four qubits assigned to each spin configuration.

For VQC1, we selected a strongly entangling layer configuration as prescribed
in ref. [28]. We observed that maximising the entanglement structure within
the circuit was crucial for accurately representing the entropy of the model.
Other ansatze did not yield satisfactory results.

For VQC2, we adopted the variational circuit structure proposed in ref. [14].
This configuration implements a Trotterised Hamiltonian time evolution
sequence as a variational ansatz time step for each term, which is implemented
as a variational parameter. One can write eq. (1) as
$\hat{U}_{2}(\theta_{2})=\left(e^{i\theta_{2}\hat{H}}\right)^{L}$:

$\displaystyle\hat{U}_{T}(\theta_{2})$ $\displaystyle=$
$\displaystyle\prod_{\langle
i,j\rangle,s}e^{i\left(\theta^{0}_{ijs}X_{i,s}\otimes
X_{j,s}+\theta^{1}_{ijs}Y_{i,s}\otimes Y_{j,s}\right)}\ ,$ (9)
$\displaystyle\hat{U}_{U}(\theta_{2})$ $\displaystyle=$
$\displaystyle\prod_{i,s}e^{i\left(\theta^{0}_{is}Z_{i,s}\otimes
Z_{i,s+1}+\theta^{1}_{is}Z_{i,s}\right)}\ ,$ (10)

where $\\{\theta^{0,1}_{ijs},\theta^{0,1}_{is}\\}\in\theta_{2}$ each being
independent, and $X,\ Y$ and $Z$ represent the Pauli operators, and subscripts
indicate the qubits that they are acting on. Notice that $\theta_{2}$
indicates the collective set of parameters included within VQC2, each term has
an independent parameter. The complete ansatz per layer consists of

$\displaystyle\hat{U}_{2}(\theta_{2})=\prod^{L}_{k=1}\hat{U}_{T}(\theta^{(k)}_{2})\hat{U}_{U}(\theta^{(k)}_{2})\
,$

where $L$ refers to the number of independent layers. The difference between
eqns. (9)111Also known as Fermionic SWAP gate, see ref. [29]. and (10) are
coming from hopping terms and repulsive term in the Hubbard model
respectively. Here, subindices $i,\ j$ correspond to the site index, $s$ to
the spin index and $k$ is the layer index. This step ensures the preservation
of the physical properties of spin types and the correlation between spin
flavours induced by on-site repulsive coupling. We refer the reader to
Appendix A for details on the circuit representation.

For system optimisation, we used SciPy’s minimiser based on Sequential Least
Squares Programming [30] (version 1.10.0). The quantum circuit and its
gradient were simulated using PennyLane [31] (version 0.32.0) and Qiskit’s AER
simulator [32] (qiskit.aer plugin version 0.28.0, Qiskit version 0.24.1). To
evaluate the quality of the final density matrix, we computed the fidelity of
the reconstructed density matrix with respect to the one computed via exact
diagonalisation:

$\displaystyle\text{Fidelity}(\rho,\rho_{\text{rec}})=\left(\mathrm{Tr}\left[\sqrt{\sqrt{\rho_{\text{rec}}}\rho\sqrt{\rho_{\text{rec}}}}\right]\right)^{2}\
,$

where $\rho_{\text{rec}}$ represents the reconstructed density matrix.

Figure 1: Quality metrics computed for fixed $U$ and $\mu$ values for four-
site Fermi-Hubbard model. The blue line in each panel shows the results
computed via exact diagonalisation, where dots represent the results
reconstructed via quantum simulation. The coloured dots in the top panel show
the Fidelity value achieved for each point. From top to bottom, the panels
compare free energy, Hamiltonian density, entropy and number operator density.
The error bars represent the statistical uncertainty for 3000 shots, and the
black dashed line shows the expectation value of the operators at
$\beta\to\infty$.

To investigate the influence of temperature on the outcomes, we conducted a
scan across the inverse temperature range $\beta\in[0.05,35]$ while
maintaining fixed values for $U=0.8$ and $\mu=0.2$. The VQC prescription
mentioned earlier represents a single layer, and we implemented four such
layers for each VQC as an initial ansatz. The results of this scan for four
different metrics are illustrated in Fig. 1.

In the top panel of the figure, the computed free energy of the system via
exact diagonalisation is depicted by the blue line, while coloured dots
represent the optimised reconstructed free energy. The colour of these dots
reflects the fidelity between the exact and reconstructed density matrix.

Within the range $\beta\in[0.05,2]$, we observed a minimum fidelity of 93%,
indicating considerable success with our initial four-layer proposal. However,
for higher $\beta$ values, despite accurately reconstructing the free energy,
we noticed a significant decrease in fidelity. This decline is attributed to
the reduced significance of the entropy term in the free energy as $\beta$
increases. Consequently, the system encountered greater difficulty converging
to the correct entropy value with increasing $\beta$. Notably, while some
values matched well, they did not reflect achieving physical equilibrium. In
order to assess the importance of such large $\beta$ values, we also include
the expectation value of the operators at $T=0$, shown via a black dashed
line. This indicates the system behaves like a ground state at $\beta\gtrsim
10$.

To independently assess the physical accuracy of the system, we utilised the
number operator density, $\sum_{i}\hat{n}_{i}$, revealing considerable
deviations from its exact value at higher $\beta$ values. This highlights the
necessity of considering physical equilibrium beyond just achieving specific
numerical accuracies. The results within $\beta>10$ region have been
constructed by increasing the number of layers of the ansatz. However, we did
not observe significant improvement in the fidelity for up to five layers.

Figure 2: The top panels show the results computed with exact diagonalisation,
whereas the bottom panels show the same results reconstructed with the
variational quantum algorithm at a fixed temperature, $T=2$ on the left and
$T=0.\bar{3}$ on the right. The colour in each panel shows the number operator
density per $(U,\mu)$ grid.

In order to assess the possible effects on more generic scans, we computed the
density of the number operator at a fixed temperature for $U,\ \mu\in[0.1,1]$.
Figure 2 shows the number density value (represented by the change in colour)
for fixed $\beta$ values, 0.5 and 3, respectively. The top panel in each
figure shows the result from exact diagonalization, whereas the bottom panel
shows the results reconstructed by quantum simulation. We observe reasonable
results for both $\beta$ values, matching the results from exact
diagonalization.

###  III.1 Cost estimation

The variational system under consideration exhibits two sources of
computational complexity. The first arises from computational errors within
the quantum device, while the second stems from the optimization landscape. We
quantified the CNOT gates within each quantum circuit to estimate the error
source. However, it is important to note that the precision required for
implementing rotation gates would amplify the number of T gate implementations
within the circuit, thereby introducing an additional source of error. For
simplicity, we focused solely on the CNOT count.

The generalized CNOT count for both circuits in terms of the number of sites
($N$) and layers ($L$) can be expressed as follows:

$\displaystyle\\#{\rm CNOT}_{\rm VQC_{1}}$ $\displaystyle=$ $\displaystyle
2LN\ ,$ (11) $\displaystyle\\#{\rm CNOT}_{\rm VQC_{2}}$ $\displaystyle=$
$\displaystyle 2L(5N-4)\ .$ (12)

These equations demonstrate that both circuits scale the CNOT count linearly
concerning both $N$ and $L$. Here, eq. (11) denotes the CNOT count in VQC1,
while eq. (12) represents VQC2.

Furthermore, an essential aspect of optimizing the system involves estimating
the presence of barren plateaus. This was accomplished by computing the
variance in the free energy based on 500 different $\theta$ values drawn from
a uniform distribution with $\theta_{1,2}\in[-\pi,\pi]$.222It has been shown
that the variance in gradient is equivalent to the variance in the objective
function [33]. Fig. 3 illustrates the variance change concerning the number of
layers and sites, where the number of qubits is twice the number of sites. Our
observations revealed that the variance in the free energy remains around
$\mathcal{O}(10^{-3})$ up to 6 layers and sites without a significant decline.
This suggests that the algorithm can be optimized for a relatively large
number of qubits and layers. It is essential to emphasise that using ansatze,
which is typically used for quantum machine learning applications, has led to
significantly lower values for the variance in free energy. This finding
reveals the importance of using physics-inspired ansatz to prevent or reduce
the risk of barren plateaus in variational quantum simulation applications.

Figure 3: Changes in the variance of the free energy with respect to the
number of sites and layers in the variational circuits where both circuits are
assumed to have the same number of layers. Variance has been estimated by
drawing 500 samples for $\theta_{1}$ and $\theta_{2}$.

##  IV Conclusion

In this study, we investigated the feasibility of variational quantum
algorithms in examining the thermal behaviour of the Fermi-Hubbard model. Our
analysis demonstrates that employing qVQT alongside a physics-inspired ansatz
proves to be a viable approximation technique. Notably, while achieving a
relatively accurate depiction of the free energy value at a given temperature
is possible with any ansatz, our findings underscore that the physics-inspired
ansatz offers superior approximation, effectively capturing the system’s
physical behaviour. To gauge this, we utilized the number density operator as
an independent metric, separate from the optimization process. We showed that,
although higher $\beta$ values pose a significant challenge for the
optimisation algorithm due to the large suppression of the entropy term, it is
possible to increase the fidelity of the system by increasing the number of
layers in the VQC1.

However, the Fermi-Hubbard model’s gapless nature necessitates a substantial
number of sites to capture the lattice’s phase transition behaviour.
Unfortunately, our computational resources limited our ability to explore
larger site numbers. Nonetheless, our estimations indicate that the proposed
ansatz holds promise in providing a relatively barren plateau-free environment
for significantly larger lattices. Implementing such approaches on tensor
networks (TNs) demands substantial resources due to the exponential expansion
in the bond dimension during the state’s imaginary time evolution (e.g., TEBD
algorithm). Despite the significant limitations of current quantum computers,
they are not anticipated to suffer from bond dimension constraints.

## Acknowledgements

JYA acknowledges the hospitality of Cambridge University. JYA is supported by
the U.S. Department of Energy, Office of Science, Contract No. DE-
AC05-06OR23177, under which Jefferson Science Associates, LLC operates
Jefferson Lab and in part by the DOE, Office of Science, Office of Nuclear
Physics, Early Career Program under contract No DE-SC0024358. MW acknowledges
support from STFC consolidated grant ST/X000664/1 and an IPPP Associateship.

##  Appendix A Physics-inspired ansatz

The ansatz that is proposed for VQC2 in eqs. (9) and (10) has been inspired by
the Trotterised time evolution steps of Fermi-Hubbard model in eq. (1). The
circuit version of the three main terms has been represented in Fig. 4. These
terms can be realised using similarity transformations where
$Y=SHZHS^{\dagger}$ and $Z=HXH$. Here $H$ stands for Hadamard gate, and $S$
stands for phase gate with $\lambda=\pi/2$. Notice that $SHS$ can also be
expressed as $R_{X}(-\pi/2)=R_{X}^{\dagger}(\pi/2)$.

Figure 4: Individual terms of the physics-inspired ansatz have been
represented as pieces of the quantum circuit. $\theta\in[-\pi,\pi]$ represents
trainable parameter in each term and $R_{m},\ m\in\\{X,\ Z\\}$ represents
rotation in respective Pauli axis.

## References

  * Hubbard [1963] J. Hubbard, Math. Phys. Sci. 276, 238 (1963).
  * Anderson [1963] P. W. Anderson, Solid State Physics 14, 99 (1963).
  * Šimkovic  _et al._ [2017] F. Šimkovic, Y. Deng, N. V. Prokof’ev, B. V. Svistunov, I. S. Tupitsyn, and E. Kozik, Phys. Rev. B 96, 081117 (2017).
  * Araz  _et al._ [2023] J. Y. Araz, S. Schenk, and M. Spannowsky, Phys. Rev. A 107, 032619 (2023), arXiv:2210.03679 [quant-ph] .
  * Paeckel  _et al._ [2019] S. Paeckel, T. Köhler, A. Swoboda, S. R. Manmana, U. Schollwöck, and C. Hubig, Annals Phys. 411, 167998 (2019), arXiv:1901.05824 [cond-mat.str-el] .
  * Agasti [2020] S. Agasti, Journal of Physics Communications 4, 015002 (2020).
  * Verstraete  _et al._ [2004] F. Verstraete, J. J. García-Ripoll, and J. I. Cirac, Phys. Rev. Lett. 93, 207204 (2004).
  * Urbanek and Soldán [2016] M. Urbanek and P. Soldán, Computer Physics Communications 199, 170 (2016).
  * Davoudi  _et al._ [2023] Z. Davoudi, N. Mueller, and C. Powers, Phys. Rev. Lett. 131, 081901 (2023), arXiv:2208.13112 [hep-lat] .
  * Peruzzo  _et al._ [2014] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’Brien, Nature Communications 5, 4213 (2014).
  * Tilly  _et al._ [2022] J. Tilly _et al._ , Phys. Rept. 986, 1 (2022), arXiv:2111.05176 [quant-ph] .
  * Cerezo  _et al._ [2021] M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, and P. J. Coles, Nature Reviews Physics 3, 625 (2021).
  * Stanisic  _et al._ [2022] S. Stanisic, J. L. Bosse, F. M. Gambetta, R. A. Santos, W. Mruczkiewicz, T. E. O’Brien, E. Ostby, and A. Montanaro, Nature Commun. 13, 5743 (2022), arXiv:2112.02025 [quant-ph] .
  * Consiglio  _et al._ [2022] M. Consiglio, W. J. Chetcuti, C. Bravo-Prieto, S. Ramos-Calderer, A. Minguzzi, J. I. Latorre, L. Amico, and T. J. G. Apollaro, J. Phys. A 55, 265301 (2022), arXiv:2106.15552 [quant-ph] .
  * Verdon  _et al._ [2019] G. Verdon, J. Marks, S. Nanda, S. Leichenauer, and J. Hidary, Quantum hamiltonian-based models and the variational quantum thermalizer algorithm (2019).
  * Guo  _et al._ [2023] X.-Y. Guo _et al._ , Chin. Phys. B 32, 010307 (2023), arXiv:2107.06234 [quant-ph] .
  * Fromm  _et al._ [2023] M. Fromm, O. Philipsen, M. Spannowsky, and C. Winterowd, (2023), arXiv:2306.06057 [hep-lat] .
  * Araz and Spannowsky [2022] J. Y. Araz and M. Spannowsky, (2022), arXiv:2211.03803 [quant-ph] .
  * Foldager  _et al._ [2022] J. Foldager, A. Pesah, and L. K. Hansen, Scientific Reports 12, 3862 (2022).
  * Selisko  _et al._ [2022] J. Selisko, M. Amsler, T. Hammerschmidt, R. Drautz, and T. Eckl, (2022), arXiv:2208.07621 [quant-ph] .
  * Warren  _et al._ [2022] A. Warren, L. Zhu, N. J. Mayhall, E. Barnes, and S. E. Economou, (2022), arXiv:2203.12757 [quant-ph] .
  * Zhang  _et al._ [2021] D.-B. Zhang, G.-Q. Zhang, Z.-Y. Xue, S.-L. Zhu, and Z. Wang, Physical Review Letters 127, 10.1103/physrevlett.127.020502 (2021).
  * Lee and Schmidt [1992] M. A. Lee and K. E. Schmidt, Computer in Physics 6, 192 (1992), https://pubs.aip.org/aip/cip/article-pdf/6/2/192/12213157/192_1_online.pdf .
  * Hirsch  _et al._ [1982] J. E. Hirsch, R. L. Sugar, D. J. Scalapino, and R. Blankenbecler, Phys. Rev. B 26, 5033 (1982).
  * Ceperley and Alder [1984] D. M. Ceperley and B. J. Alder, The Journal of Chemical Physics 81, 5833 (1984), https://pubs.aip.org/aip/jcp/article-pdf/81/12/5833/14120713/5833_1_online.pdf .
  * Shastry [1986] B. S. Shastry, Phys. Rev. Lett. 56, 1529 (1986).
  * Reiner  _et al._ [2016] J.-M. Reiner, M. Marthaler, J. Braumüller, M. Weides, and G. Schön, Phys. Rev. A 94, 032338 (2016).
  * Schuld  _et al._ [2020] M. Schuld, A. Bocharov, K. M. Svore, and N. Wiebe, Physical Review A 101, 10.1103/physreva.101.032308 (2020).
  * Hashim  _et al._ [2022] A. Hashim, R. Rines, V. Omole, R. K. Naik, J. M. Kreikebaum, D. I. Santiago, F. T. Chong, I. Siddiqi, and P. Gokhale, Phys. Rev. Res. 4, 033028 (2022), arXiv:2111.04572 [quant-ph] .
  * Virtanen  _et al._ [2020] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, Nature Methods 17, 261 (2020).
  * Bergholm  _et al._ [2020] V. Bergholm, J. Izaac, M. Schuld, C. Gogolin, M. S. Alam, S. Ahmed, J. M. Arrazola, C. Blank, A. Delgado, S. Jahangiri, K. McKiernan, J. J. Meyer, Z. Niu, A. Száva, and N. Killoran,  (2020), arXiv:1811.04968 [quant-ph] .
  * ANIS  _et al._ [2021] M. S. ANIS, H. Abraham, AduOffei, R. Agarwal, G. Agliardi, M. Aharoni, I. Y. Akhalwaya, G. Aleksandrowicz, T. Alexander, M. Amy, S. Anagolum, E. Arbel, A. Asfaw, A. Athalye, A. Avkhadiev, C. Azaustre, P. BHOLE, A. Banerjee, S. Banerjee, W. Bang, A. Bansal, P. Barkoutsos, A. Barnawal, G. Barron, G. S. Barron, L. Bello, Y. Ben-Haim, M. C. Bennett, D. Bevenius, D. Bhatnagar, A. Bhobe, P. Bianchini, L. S. Bishop, C. Blank, S. Bolos, S. Bopardikar, S. Bosch, S. Brandhofer, Brandon, S. Bravyi, N. Bronn, Bryce-Fuller, D. Bucher, A. Burov, F. Cabrera, P. Calpin, L. Capelluto, J. Carballo, G. Carrascal, A. Carriker, I. Carvalho, A. Chen, C.-F. Chen, E. Chen, J. C. Chen, R. Chen, F. Chevallier, K. Chinda, R. Cholarajan, J. M. Chow, S. Churchill, CisterMoke, C. Claus, C. Clauss, C. Clothier, R. Cocking, R. Cocuzzo, J. Connor, F. Correa, A. J. Cross, A. W. Cross, S. Cross, J. Cruz-Benito, C. Culver, A. D. Córcoles-Gonzales, N. D, S. Dague, T. E. Dandachi, A. N. Dangwal, J. Daniel, M. Daniels, M. Dartiailh, A. R. Davila, F. Debouni, A. Dekusar, A. Deshmukh, M. Deshpande, D. Ding, J. Doi, E. M. Dow, E. Drechsler, E. Dumitrescu, K. Dumon, I. Duran, K. EL-Safty, E. Eastman, G. Eberle, A. Ebrahimi, P. Eendebak, D. Egger, ElePT, Emilio, A. Espiricueta, M. Everitt, D. Facoetti, Farida, P. M. Fernández, S. Ferracin, D. Ferrari, A. H. Ferrera, R. Fouilland, A. Frisch, A. Fuhrer, B. Fuller, M. GEORGE, J. Gacon, B. G. Gago, C. Gambella, J. M. Gambetta, A. Gammanpila, L. Garcia, T. Garg, S. Garion, J. R. Garrison, T. Gates, L. Gil, A. Gilliam, A. Giridharan, J. Gomez-Mosquera, Gonzalo, S. de la Puente González, J. Gorzinski, I. Gould, D. Greenberg, D. Grinko, W. Guan, D. Guijo, J. A. Gunnels, H. Gupta, N. Gupta, J. M. Günther, M. Haglund, I. Haide, I. Hamamura, O. C. Hamido, F. Harkins, K. Hartman, A. Hasan, V. Havlicek, J. Hellmers, Ł. Herok, S. Hillmich, H. Horii, C. Howington, S. Hu, W. Hu, J. Huang, R. Huisman, H. Imai, T. Imamichi, K. Ishizaki, Ishwor, R. Iten, T. Itoko, A. Ivrii, A. Javadi, A. Javadi-Abhari, W. Javed, Q. Jianhua, M. Jivrajani, K. Johns, S. Johnstun, Jonathan-Shoemaker, JosDenmark, JoshDumo, J. Judge, T. Kachmann, A. Kale, N. Kanazawa, J. Kane, Kang-Bae, A. Kapila, A. Karazeev, P. Kassebaum, J. Kelso, S. Kelso, V. Khanderao, S. King, Y. Kobayashi, Kovi11Day, A. Kovyrshin, R. Krishnakumar, V. Krishnan, K. Krsulich, P. Kumkar, G. Kus, R. LaRose, E. Lacal, R. Lambert, H. Landa, J. Lapeyre, J. Latone, S. Lawrence, C. Lee, G. Li, J. Lishman, D. Liu, P. Liu, L. Madden, Y. Maeng, S. Maheshkar, K. Majmudar, A. Malyshev, M. E. Mandouh, J. Manela, Manjula, J. Marecek, M. Marques, K. Marwaha, D. Maslov, P. Maszota, D. Mathews, A. Matsuo, F. Mazhandu, D. McClure, M. McElaney, C. McGarry, D. McKay, D. McPherson, S. Meesala, D. Meirom, C. Mendell, T. Metcalfe, M. Mevissen, A. Meyer, A. Mezzacapo, R. Midha, D. Miller, Z. Minev, A. Mitchell, N. Moll, A. Montanez, G. Monteiro, M. D. Mooring, R. Morales, N. Moran, D. Morcuende, S. Mostafa, M. Motta, R. Moyard, P. Murali, J. Müggenburg, T. NEMOZ, D. Nadlinger, K. Nakanishi, G. Nannicini, P. Nation, E. Navarro, Y. Naveh, S. W. Neagle, P. Neuweiler, A. Ngoueya, J. Nicander, Nick-Singstock, P. Niroula, H. Norlen, NuoWenLei, L. J. O’Riordan, O. Ogunbayo, P. Ollitrault, T. Onodera, R. Otaolea, S. Oud, D. Padilha, H. Paik, S. Pal, Y. Pang, A. Panigrahi, V. R. Pascuzzi, S. Perriello, E. Peterson, A. Phan, F. Piro, M. Pistoia, C. Piveteau, J. Plewa, P. Pocreau, A. Pozas-Kerstjens, R. Pracht, M. Prokop, V. Prutyanov, S. Puri, D. Puzzuoli, J. Pérez, Quant02, Quintiii, I. R, R. I. Rahman, A. Raja, R. Rajeev, N. Ramagiri, A. Rao, R. Raymond, O. Reardon-Smith, R. M.-C. Redondo, M. Reuter, J. Rice, M. Riedemann, Rietesh, D. Risinger, M. L. Rocca, D. M. Rodríguez, RohithKarur, B. Rosand, M. Rossmannek, M. Ryu, T. SAPV, N. R. C. Sa, A. Saha, A. Ash-Saki, S. Sanand, M. Sandberg, H. Sandesara, R. Sapra, H. Sargsyan, A. Sarkar, N. Sathaye, B. Schmitt, C. Schnabel, Z. Schoenfeld, T. L. Scholten, E. Schoute, M. Schulterbrandt, J. Schwarm, J. Seaward, Sergi, I. F. Sertage, K. Setia, F. Shah, N. Shammah, R. Sharma, Y. Shi, J. Shoemaker, A. Silva, A. Simonetto, D. Singh, D. Singh, P. Singh, P. Singkanipa, Y. Siraichi, Siri, J. Sistos, I. Sitdikov, S. Sivarajah, M. B. Sletfjerding, J. A. Smolin, M. Soeken, I. O. Sokolov, I. Sokolov, V. P. Soloviev, SooluThomas, Starfish, D. Steenken, M. Stypulkoski, A. Suau, S. Sun, K. J. Sung, M. Suwama, O. Słowik, H. Takahashi, T. Takawale, I. Tavernelli, C. Taylor, P. Taylour, S. Thomas, K. Tian, M. Tillet, M. Tod, M. Tomasik, C. Tornow, E. de la Torre, J. L. S. Toural, K. Trabing, M. Treinish, D. Trenev, TrishaPe, F. Truger, G. Tsilimigkounakis, D. Tulsi, W. Turner, Y. Vaknin, C. R. Valcarce, F. Varchon, A. Vartak, A. C. Vazquez, P. Vijaywargiya, V. Villar, B. Vishnu, D. Vogt-Lee, C. Vuillot, J. Weaver, J. Weidenfeller, R. Wieczorek, J. A. Wildstrom, J. Wilson, E. Winston, WinterSoldier, J. J. Woehr, S. Woerner, R. Woo, C. J. Wood, R. Wood, S. Wood, J. Wootton, M. Wright, L. Xing, J. YU, B. Yang, D. Yeralin, R. Yonekura, D. Yonge-Mallo, R. Yoshida, R. Young, J. Yu, L. Yu, C. Zachow, L. Zdanski, H. Zhang, I. Zidaru, C. Zoufal, aeddins ibm, alexzhang13, b63, bartek bartlomiej, bcamorrison, brandhsn, charmerDark, deeplokhande, dekel.meirom, dime10, dlasecki, ehchen, fanizzamarco, fs1132429, gadial, galeinston, georgezhou20, georgios ts, gruu, hhorii, hykavitha, itoko, jessica angel7, jezerjojo14, jliu45, jscott2, klinvill, krutik2966, ma5x, michelle4654, msuwama, ntgiwsvp, ordmoj, sagar pahwa, pritamsinha2304, ryancocuzzo, saswati qiskit, septembrr, sethmerkel, shaashwat, sternparky, strickroman, tigerjack, tsura crisaldo, vadebayo49, welien, willhbang, wmurphy collabstar, yang.luh, and M. Čepulkovskis, Qiskit: An open-source framework for quantum computing (2021).
  * Ragone  _et al._ [2023] M. Ragone, B. N. Bakalov, F. Sauvage, A. F. Kemper, C. O. Marrero, M. Larocca, and M. Cerezo, (2023), arXiv:2309.09342 [quant-ph] .