††thanks<EMAIL_ADDRESS>
# Phase diagram of quantum generalized Potts-Hopfield neural networks

Eliana Fiorelli  Institute for Theoretical Nanoelectronics (PGI-2),
Forschungszentrum Jülich, 52428 Jülich, Germany  JARA-Institute for Quantum
Information, RWTH Aachen University, 52056 Aachen, Germany  Igor Lesanovsky
School of Physics and Astronomy, University of Nottingham, Nottingham, NG7
2RD, UK  Centre for the Mathematics and Theoretical Physics of Quantum Non-
equilibrium Systems, University of Nottingham, Nottingham NG7 2RD, UK
Institut für Theoretische Physik, Universität Tübingen, Auf der Morgenstelle
14, 72076 Tübingen, Germany  Markus Müller  Institute for Theoretical
Nanoelectronics (PGI-2), Forschungszentrum Jülich, 52428 Jülich, Germany
JARA-Institute for Quantum Information, RWTH Aachen University, 52056 Aachen,
Germany

###### Abstract

We introduce and analyze an open quantum generalization of the q-state Potts-
Hopfield neural network, which is an associative memory model based on multi-
level classical spins. The dynamics of this many-body system is formulated in
terms of a Markovian master equation of Lindblad type, which allows to
incorporate both probabilistic classical and coherent quantum processes on an
equal footing. By employing a mean field description we investigate how
classical fluctuations due to temperature and quantum fluctuations effectuated
by coherent spin rotations affect the ability of the network to retrieve
stored memory patterns. We construct the corresponding phase diagram, which in
the low temperature regime displays pattern retrieval in analogy to the
classical Potts-Hopfield neural network. When increasing quantum fluctuations,
however, a limit cycle phase emerges, which has no classical counterpart. This
shows that quantum effects can qualitatively alter the structure of the
stationary state manifold with respect to the classical model, and potentially
allow one to encode and retrieve novel types of patterns.

##  I Introduction

Machine Learning (ML) is considered the core of artificial intelligence and
data science, and it is nowadays an expanding research area. The interest in
giving systems the ability to learn how to accomplish a task, without being
explicitly programmed, ranges from computer science to neurobiology [1, 2],
with applications spreading from science to commerce [3, 4]. One of the most
successful architectures for ML is represented by neural networks (NNs) [5],
which are artificial systems that are able to mimic what is known about the
inner workings of the brain in relation to assimilation and comprehension
processes. While in general these models may be too simple for the pursuit of
understanding the computational properties of the brain itself [6], they have
nevertheless powerful applications in artificial intelligence. NNs are systems
described in terms of interconnected artificial neurons with learning
capabilities. Among the several instances of NNs, many efforts are focused on
attractor NNs, which are well suited to be described via classical spin
systems subject to thermal fluctuations [7]. This formulation enables one to
adopt concepts commonly used in the context of statistical mechanics.

In the last decades, it has been shown that genuine quantum features such as
superposition of states and entanglement provide remarkable advantages for
solving many computational problems [8, 9, 10, 11], which has led to
significant developments in the field of quantum computation. Within this
research area, a number of current developments focus on exploring possible
advantages in the design and realization of near-term quantum devices by
applying concepts from classical NN computing to quantum systems [12]. One
example is given by suitably constructed classical NNs, which are used to
approximate many-body quantum states, such as Restricted Boltzmann Machines
[13, 14, 15]. Other recent works deal with more sophisticated architectures,
such as deep NNs, for efficiently modelling quantum many-body states [16, 17,
18], allowing to implement representations of mixed states, rather than pure
states only [19]. Also noteworthy in this regard are recent advances in
relating the NNs description of many-body quantum states to tensor network
states, with the intention to exploit the latter as a tool for ML tasks [20,
21, 22, 23].

In addition to the aforementioned paradigms, various efforts aim at harnessing
the potential computational power of quantum generalizations of classical NNs
[24, 25]. Here the goal is to first understand how to model such quantum
versions, and subsequently to analyze the emergence of novel features embodied
in their quantum nature [26]. Current instances of quantum formulations of NNs
range from quantum algorithms and quantum circuit settings [27, 28, 29, 30,
31, 32, 33, 34, 35], which include, e.g., the so-called feed-forward quantum
NNs, to condensed matter systems [36, 37, 38] and biological settings [39,
40]. One approach for generalizing NNs into the quantum domain — which we will
also pursue here — views a task-performing NN as a dynamical process that
takes place within an ensemble of neurons that are represented by quantum
systems [41, 42, 43]. This perspective is inspired by Hopfield-type NNs [44],
which are, amongst attractor NNs, the simplest instances of associative
memories, permitting the recognition, or retrieval, of patterns. Here, the
retrieval of information is described by means of a classical non-equilibrium
dynamics, where memories are stored as long-time solutions of a stochastic
dynamics [7, 45, 46]. Such a classical out-of-equilibrium dynamics can be
generalized into the quantum domain via the formalism of open quantum systems
[47], which allows to introduce quantum generalizations of Hopfield-type NNs.
This idea can be further generalized, as recently done in Ref. [48], where
quantum NNs are represented as dynamical maps. They evolve initial quantum
states, which play the role of inputs of the quantum NN, and take them towards
stored memories, which correspond to stationary states of the maps.

In this work we introduce a quantum generalized NN, which is derived from the
classical q-state Potts-Hopfield neural network (qPHN). Its basic constituents
are multi-level spins which are subject to an all-to-all interaction, as shown
in Fig. 1(a). We introduce an open quantum dynamics that evolves this system
by means of a Markovian dynamical map, using the Gorini-Kossakowski-Sudarshan-
Lindblad (GKSL) quantum master equation [49]. This allows us to include the
following competing effects: (i) dissipative, non-unitary processes represent
the analogy with the classical out-of-equilibrium pattern-retrieval dynamics;
(ii) coherent, unitary processes, which generate quantum fluctuations due to
the creation of superpositions between spin states. In the context of quantum
generalization of classical dynamical maps, the competition between coherent
and dissipative dynamics has been explored in different settings, often
providing the access to new classes of non-equilibrium many-body quantum
systems and/or to novel phases of matter [50, 51]. Here, in order to
investigate the impact of (ii) on the retrieval mechanism of the NN, we focus
on analyzing the long-time state of the dynamics by means of a mean field
approach. The derived phase diagram shows that the purely dissipative case,
where coherent effects are absent, is indeed consistent with the classical
pattern retrieval phenomenology. In the presence of quantum fluctuations the
NN retains the capability to retrieve patterns. However, we find that due to
the competition between dissipation and coherent processes a new limit cycle
phase emerges. A similar phenomenon was previously observed within a quantum
generalized Hopfield Neural Network (HNN) [41], and we discuss similarities
and differences with respect to the qPHN. Our study corroborates that quantum
generalized Hopfield-type NNs feature novel phases which could potentially be
exploited to define new types of retrieval. Moreover, they provide a setting
for the realization of quantum maps with multiple stationary states, and
therefore yield a physical platform for addressing questions concerning a
possibly enhanced storage capacity of quantum NNs, as discussed in Ref. [48].

Figure 1: Quantum effects in the q-Potts-Hopfield neural network. Scheme of
the quantum generalization of the qPHN model. (a) All-to-all connected network
of $N=4$ Potts spins. The $i$-th and $j$-th spins (neurons) are in the state
$\alpha$, $\beta$, respectively, and interact with coefficient
$J_{ij}^{\alpha\beta}$. Here we set $q=3$, so that each neuron can assume the
values $\\{1,e^{i2\pi/3},e^{-i2\pi/3}\\}$. A classical stochastic dynamics
leading to an equilibrium state at temperature $T=1/\beta$ can be implemented
via a purely dissipative Lindblad master equation. The thermal transitions
between the internal states of a Potts spin occur with rate
$\Gamma_{\alpha,\pm}^{i}$, which depends on the configuration of all the other
neurons (see main text for details). Quantum effects are included by allowing
coherent transitions (parametrized by the rate $\lambda$) between states of
the Potts spins. (b) Sketch of the energy landscape of the fully qPHN. Memory
states, i.e. patterns labeled as $\vec{\xi}^{\mu}$, are global minima of the
energy function. In case of a purely dissipative dynamics, i.e. $\lambda=0$,
the stationary state is a thermal one, $\rho_{ss}$, with respect to the energy
function $E$ [see Eq. (2)]. When $\lambda\neq 0$, non-classical stationary
states and also non-stationary long time behavior, such as limit cycle
solutions, may emerge. (c) As a function of the temperature, the qPHN passes
from a disordered, paramagnetic phase to an ordered one where patterns are
retrieved. The transition occurs discontinuously at $T_{c}$, and it is
signalled by the behavior of the overlap between Potts spin configurations and
patterns [see Eq. (8)]. Within the temperature range $T_{m}\leq T\leq T_{M}$ a
phase coexistence of paramagnetic and retrieval solutions takes place.

##  II q-state Potts-Hopfield neural network

In this section we review the main features of the classical model that
describes a multi-level neural network referred to as q-state Potts-Hopfield
neural network (qPHN) [52, 53]. The basic constituents of this NN are $N$
multi-level classical spins, $\sigma_{i}$, $i=1,...,N$, referred to as Potts
spins. A single Potts spin can assume $q$ possible states which may be labeled
by the roots of unity:

$\sigma_{i}=\omega^{k_{i}},\;\omega=e^{i\frac{2\pi}{q}},\;k_{i}=0,...,q-1,$
(1)

$\forall i=1,...,N$. Information is stored in the form of special Potts spin
configurations, or patterns, denoted by the vectors
$\vec{\xi}^{\mu}=(\xi^{\mu}_{1},...,\xi^{\mu}_{N})$, with $\mu=1,...,p$. Each
component $\xi^{\mu}_{i}$ of a pattern can assume the values
$\\{\omega^{k_{i}^{\mu}}\\}_{k_{i}^{\mu}=0,...,q-1}$. The interaction between
Potts spins is governed by the energy functional

$E=-\frac{1}{2}\sum_{i\neq
j}^{N}\sum_{\alpha,\beta=1}^{q}J_{ij}^{\alpha\beta}u_{\sigma_{i},\alpha}u_{\sigma_{j},\beta},$
(2)

where $u_{\sigma_{i},\sigma_{j}}\equiv q\delta_{\sigma_{i},\sigma_{j}}-1$. The
interaction strengths are parametrized by the symmetric connectivity matrix

$J_{ij}^{\alpha\beta}=\frac{1}{q^{2}N}\sum_{\mu=1}^{p}u_{\xi_{i}^{\mu},\alpha}u_{\xi_{j}^{\mu},\beta}.$
(3)

A scheme of such an all-to-all network is given in Fig. 1(a). With these
definitions, and in the limit $p/N\ll 1$, i.e. a small number of patterns with
respect to the network size, the patterns $\vec{\xi}^{\mu}$ are minima of the
energy functional (2), as illustrated in Fig. 1(b). The ratio of the number of
patterns and Potts spins — $\alpha=p/N$ — is referred to as load parameter.
When it increases above a critical threshold $\alpha_{c}$ spurious minima
emerge, and the hence minima of the energy function can no longer be
associated with stored patterns. The critical load paramater $\alpha_{c}$ thus
yields the storage capacity, i.e. the maximum number of patterns that can be
stored in a network of $N$ Potts spins.

For $\alpha<\alpha_{c}$, a pattern retrieval mechanism can be constructed via
a classical discrete-time gradient descent dynamics with respect to the energy
(2). Under this dynamics the $i$-th Potts spin, in a generic configuration
$\sigma_{i}(t)$ at time $t$, is updated to a lower energy configuration,
$\sigma^{\prime}_{i}(t+1)$, at time $t+1$:

$\sigma_{i}(t)\rightarrow\sigma^{\prime}_{i}(t+1)\Leftrightarrow\Delta
E_{i}=h_{\sigma^{\prime}_{i}}-h_{\sigma_{i}}<0.$ (4)

Here $\Delta E_{i}$ is the energy difference associated with the transition
$\sigma_{i}\rightarrow\sigma^{\prime}_{i}$, which depends on the potential
energy of the $i$-th Potts spin in the state $\sigma_{i}$:

$h_{\sigma_{i}}=-\sum_{j\neq
i,j=1}^{N}\sum_{\alpha,\beta}J_{ij}^{\alpha\beta}u_{\sigma_{i},\alpha}u_{\sigma_{j},\beta}.$
(5)

It is worth noting that $h_{\sigma_{i}}$ depends on both the states of all the
other neurons, $j\neq i$, and the state of the $i$-th neuron itself.

The pattern retrieval via gradient descent can be interpreted as a dynamics
that takes the system of Potts spins to the equilibrium state with respect to
the energy (2) at zero temperature. Noise can be included in the NN by
introducing a finite inverse temperature $\beta=1/T$ and a stochastic dynamics
that leads the system to an equilibrium distribution $P_{eq}=e^{-\beta
E}/Z(\beta)$, with $Z(\beta)$ denoting the partition function. Let
$p(\vec{\sigma},t)$ be the probability to find the system in the configuration
$\vec{\sigma}=(\sigma_{1},...,\sigma_{N})^{T}$ at time $t$. The stochastic
dynamics with the stationary equilibrium state
$P_{eq}=\lim_{t\rightarrow\infty}p(\vec{\sigma},t)$ is then given by

$\dot{p}(\vec{\sigma})=\sum_{i=1}^{N}\sum_{\sigma_{i}^{\prime}\neq\sigma_{i}}w_{\sigma_{i}^{\prime}\rightarrow\sigma_{i}}p(\vec{\sigma}_{i}^{\prime})-w_{\sigma_{i}\rightarrow\sigma_{i}^{\prime}}p(\vec{\sigma}_{i}).$
(6)

The transition rate for changing the state of the $i$-th Potts spin is

$w_{\sigma_{i}\rightarrow\sigma_{i}^{\prime}}=\frac{1}{Z_{i}}e^{-\beta\Delta
E_{i}},$ (7)

with
$Z_{i}=\sum_{\sigma_{i}=1}^{q}\sum_{\sigma^{\prime}\neq\sigma}e^{-\beta(h_{\sigma^{\prime}_{i}}-h_{\sigma_{i}})/2}$
[54].

To analyze the pattern retrieval dynamics of the qPHN it is useful to
introduce the overlaps between the Potts spins and the patterns. These play
the role of order parameters, and their introduction drastically reduces the
complexity of the problem (instead of the $q_{N}$-dimensional probability
vector one considers an order parameter with $p$ components). The overlap of
the spin configuration with the $\mu$-th pattern is given by

$m^{\mu}=\frac{1}{N(q-1)}\sum_{i=1}^{N}u_{\xi_{i}^{\mu},\sigma_{i}}.$ (8)

Positive values of the overlap correspond to the retrieval of patterns, which
is maximum when the overlap assumes its largest value. For instance,
$m^{\mu}=1$ signals the complete retrieval of the $\mu$-th pattern. On the
contrary, negative values of the overlap do not correspond to the storage of
any of the patterns [55, 54]. For sufficiently low temperatures, i.e. below a
temperature $T_{c}$, thermal fluctuations do not affect the ability of the
network to retrieve patterns. Here, positive values of the overlap are stable
stationary solutions, allowing the retrieval of one pattern, when the
corresponding overlap, say $m^{\mu}$, assumes a positive value in the
stationary state. This characterizes the retrieval phase of the network.
Negative solutions of the overlap can also appear as stable stationary
solutions. However, for the cases that have been analysed, they disappear
first than the retrieval solutions as the temperature is increased [54]. For
$T>T_{c}$ thermal fluctuations prevent patterns to be stored as long time
solutions, and the stationary state is characterized by vanishing overlaps,
$m^{\mu}=0$, identifying a paramagnetic phase. At the critical temperature
$T_{c}$ a first order phase transition between the retrieval and the
paramagnetic phase takes place [56], where the overlap changes discontinuously
from $m^{\mu}>0$ to $m^{\mu}=0$. Moreover, for temperatures within the
interval $T\in[T_{m},T_{M}]$, with $T_{m}<T_{c}<T_{M}$, both retrieval
solutions and the paramagnetic one are present, this giving rise to a
coexistence region. Here, retrieval (paramagnetic) solutions are stable
(unstable) for $T\in[T_{m},T_{c}]$, and unstable (stable) for
$T\in[T_{c},T_{M}]$, as highlighted in Fig. 1(c). In the case of a qPHN with
$q=3$, it is $T_{m}=2$, $T_{c}=2.14$, $T_{M}=2.18$ [55], with negative overlap
solutions disappearing at $T=T_{m}$ [54], and becoming unstable at $T=0$ for
$p>1$ [55].

To conclude this section, let us focus on the particularly simple case of a
qPHN storing one pattern, i.e. $p=1$. Here, one can apply the gauge
transformation $\sigma_{i}\rightarrow\xi_{i}^{1}\sigma_{i}$, which aligns all
the components of the pattern. This situation corresponds to the standard
Potts model [57] described by the Hamiltonian
$H=-\sum_{i,j=1}^{N}J_{i,j}u_{\sigma_{i},\sigma_{j}}$. The latter further
reduces to the all-to-all Ising model for the choice $q=2$. For a generic
number $p$ of memories, the two-state Potts-Hopfield NN, $q=2$, corresponds to
the Hopfield Neural Network (HNN) [44],
$E_{H}=-\sum_{i,j=1}^{N}J_{i,j}\sigma_{i},\sigma_{j}$, with $\sigma_{i}=\\{\pm
1\\}$ and $J_{ij}$ derived from Eq. (3) by setting $q=2$. It is worth noticing
that, despite the fact that the HNN can be derived from the qPHN, and that it
features the associative memory behavior, some differences can be highlighted.
Firstly, the energy function of the HNN is characterized by a
$\mathbb{Z}_{2}$-symmetry. As a consequence, given a pattern
$\\{\xi^{\mu}_{1},...,\xi_{N}^{\mu}\\}$, its opposite one,
$\\{-\xi^{\mu}_{1},...,-\xi_{N}^{\mu}\\}$, is also a stored memory,
corresponding to negative value of the overlap. Contrarily, this symmetry is
absent in the qPHN case for $q\geq 3$, and indeed negative values of the
overlap do not correspond to storage of patterns. Secondly, when the thermal
fluctuations are considered in the HNN, it undergoes a second order phase
transition form a paramagnetic to a retrieval phase when the temperature is
lowered, instead of a first order transition that is characteristic for the
qPHN.

In the next sections, we will introduce a quantum generalization of the qPHN.
The type of patterns we will focus on in this work are $(i)$ i.i.d. random
variables, with probability distribution
$\mathbf{P}(\xi_{i}^{\mu}=\omega^{k_{i}^{\mu}})=1/q$; $(ii)$ almost orthogonal
patterns, satisfying
$\lim_{N\rightarrow+\infty}\vec{\xi}^{\mu}\cdot\vec{\xi}^{\nu}=\delta_{\mu,\nu}$,
so that they can be considered distinguishable; furthermore we will focus on
the limit of zero load, $p/N\ll 1$.

##  III Quantum model

In order to include quantum effects in the qPHN, hereafter referred to as
”quantum qPHN”, we first promote the classical variables $\sigma_{i}$ to
quantum operators $\hat{\Omega}_{i}$. Their eigenstates are $\ket{k_{i}}$,
with $k_{i}=0,...,q-1$, and the eigenvalue equation is
$\hat{\Omega}_{i}\ket{k_{i}}=\omega^{k_{i}}\ket{k_{i}}$. Transitions between
states of a q-Potts spin are effectuated by the operators $\hat{T}^{\pm}_{i}$,
such that $\hat{T}^{\pm}_{i}\ket{k_{i}}=\ket{k_{i}\pm 1}$, with
$\hat{T}_{i}^{+}\ket{(q-1)_{i}}=\ket{0_{i}}$ and
$\hat{T}_{i}^{-}\ket{0_{i}}=\ket{(q-1)_{i}}$. It is worth noting that the
Potts spin operators obey the $\mathbb{Z}_{q}$-algebra

$\displaystyle\hat{T}_{i}^{+}\hat{\Omega}_{i}=\omega\hat{\Omega}_{i}\hat{T}_{i}^{+},$
(9)
$\displaystyle\hat{T}_{i}^{-}\hat{\Omega}_{i}=\omega^{*}\hat{\Omega}_{i}\hat{T}_{i}^{-},$
(10)
$\displaystyle(\hat{T}_{i}^{+})^{q}=(\hat{T}_{i}^{-})^{q}=\hat{\Omega}_{i}^{q}=\mathbf{1}_{q}.$
(11)

The first step towards including quantum effects into the qPHN is to formulate
the classical stochastic dynamics of Eq. (6) in terms of a purely dissipative
Markovian evolution, as given by the GKSL equation

$\dot{\rho}=\mathcal{L}(\rho)=\sum_{l}\hat{L}_{l}\rho\hat{L}^{\dagger}_{l}-\frac{1}{2}\\{\hat{L}^{\dagger}_{l}\hat{L}_{l},\rho\\}.$
(12)

Here, $\rho$ is the state (density operator) of $N$ quantum q-Potts spins, and
$\hat{L}_{l}$ are the jump operators. The latter are chosen such that they
give rise to the same local processes described by the aforementioned
classical dynamics. This means that the diagonal of the density operator
$\rho$ evolves under exactly the same dynamics as the probability distribution
of the classical qPHN. The jump operators that achieve this are

$\begin{split}&\hat{L}_{\alpha,s}^{i}=\sqrt{\gamma}~{}\hat{\Gamma}_{\alpha,s}^{i}\hat{T}_{\alpha,s}^{i},\\\
&\hat{\Gamma}_{\alpha,s}^{i}=\frac{1}{\sqrt{\hat{Z}_{i}}}e^{\frac{\beta}{2}\Delta\hat{E}_{\alpha,s}^{i}},\end{split}$
(13)

for $i=1,...,N$, $\alpha=1,...,q$, $s=\pm 1$. Here the operator
$\hat{T}_{\alpha,s}^{i}$ implements the transition of the $i$-th Potts spin
from the state $\ket{\alpha}$ to the state $\ket{\alpha+s}$, and the operator
$\hat{\Gamma}_{\alpha,s}^{i}$ represents the operatorial form of the classical
transition rate defined by Eq. (7). The expression in the exponent reads

$\Delta\hat{E}_{\alpha,s}^{i}=-\frac{1}{2}\sum_{j\neq
i}^{N}\sum_{\eta,\eta^{\prime}=1}^{q-1}\left[\mathcal{J}^{\eta,\eta^{\prime}}_{ij}\hat{\Omega}^{\eta^{\prime}}_{j}\omega^{(\alpha-1)\eta}(\omega^{s\eta}-1)+\mathrm{h.c.}\right],$
(14)

with
$\mathcal{J}_{ij}^{\eta\eta^{\prime}}=\frac{1}{N}\sum_{\mu=1}^{p}(\xi_{i}^{*,\mu})^{\eta}(\xi_{j}^{*,\mu})^{\eta^{\prime}}$,
and can be regarded as the corresponding change of energy under such
transition. Analogously, the operatorial formulation of the classical
partition function reads
$\hat{Z}_{i}=\sum_{\alpha}\sum_{s}e^{-\beta\Delta\hat{E}_{\alpha,s}^{i}}$.
With these definitions, we can see that the operators
$\hat{\Gamma}_{\alpha,s}^{i}$ are non-local ones, as they depend on the Potts-
operators $\hat{\Omega}_{j}$ $\forall j\neq i$, i.e. the transition of the
$i$-th Potts-spin from the state $\ket{\alpha}$ to the state $\ket{\alpha+s}$
is ruled by a rate which depends on the state of all the rest of the network.
The expression defined by Eq. (14) can be understood as the energy difference
related to the operator

$\hat{E}=-\frac{1}{2N}\sum_{\mu=1}^{p}\left\\{\left[\sum_{i=1}^{N}\sum_{\alpha=1}^{q-1}(\xi_{i}^{*,\mu}\hat{\Omega}_{i})^{\alpha}\right]^{2}+\mathrm{h.c.}\right\\},$
(15)

which can be obtained from the classical Potts-Hopfield energy, defined by Eq.
(2), by substituting classical Potts-spin variables with their quantum
operators, furthermore requiring to satisfy hermiticity.

Using the master equation formulation of the qPHN dynamics immediately allows
the inclusion of quantum effects through a Hamiltonian. More specifically, we
consider the latter being a non-commuting term with respect to the set of jump
operators (13). This leads to the population of off-diagonal terms in the
density operator $\rho$, when written with respect to the classical basis
$\\{\otimes_{i=1}^{N}\ket{k_{i}}\\}$. For the prupose of the present work, we
choose a simple Hamiltonian [58], which reads

$\hat{H}=\lambda\sum_{i=1}^{N}\hat{T}_{i},$ (16)

with $\hat{T}^{i}=\hat{T}_{+}^{i}+\hat{T}_{-}^{i}$. Given the non-
commutativity of such a term and the rate components
$\hat{\Gamma}_{\alpha,s}^{i}$ of the jump operators (16), it acts on the Potts
spins as a ”transverse field”, i.e. it will create superpositions between
classical basis states. This leads to the formation of coherences and
concomitant quantum fluctuations. Combining the adjoint Eq. (12) with the
adjoint von Neumann equation of the Hamiltonian (16) yields the equation of
motion (EoM) of a generic operator $\hat{O}$:

$\dot{\hat{O}}=i[\hat{H},\hat{O}]+\sum_{l}\hat{L}_{l}^{\dagger}\hat{O}\hat{L}_{l}-\frac{1}{2}\\{\hat{L}_{l}^{\dagger}\hat{L}_{l},\hat{O}\\}.$
(17)

Equipped with the operatorial evolution provided by Eq. (17), we can analyze
how the competition between the classical qPHN dynamics and the Hamiltonian
dynamics modifies the retrieval mechanism of the network. As such, we consider
the functional form of the classical order parameter (8) and promote it to be
a quantum operator, getting

$\hat{m}^{\mu}=\frac{1}{2N(q-1)}\sum_{i=1}^{N}\sum_{\alpha=1}^{q-1}(\xi_{i}^{*,\mu}\hat{\Omega}_{i})^{\alpha}+\mathrm{h.c.}$
(18)

We will thus investigate the expectation value of such a quantity,
$\braket{\hat{m}^{\mu}}$, that can thus be interpreted as the overlap between
the pattern configurations and the Potts-spin one. To do so, we will consider
the stationary solutions of EoMs of the type of Eq. (17), at varying of the
coherent control parameter $\lambda$ and the temperature $T$, and setting
$\gamma=1$ in the definition Eq. (13) of the jump operators.

##  IV q=3 quantum Potts Neural Network

We now specialize the generic quantum qPHN described in Sec. III by fixing
$q$, the number of levels for each Potts spin state. In this respect, we
consider the simplest model of the quantum qPHN that goes beyond the $q=2$
case, which is equivalent to a HNN and has been explored in [58]. Thus, we
will set $q=3$, focusing on this case throughout the remainder of this work.
In this section, we explain the main steps for deriving a set of EoMs that
allow us to study the retrieval properties of the $q=3$ quantum qPHN, leaving
a more detailed description of the derivation in Appendix B.

###  IV.1 Collective operators

We aim at obtaining a closed set of EoMs for the $N$ Potts spins, by employing
the smallest possible number of operators. Such a set is given by the Potts
spin operators $\hat{\Omega}_{i}$ and $\hat{\Omega}_{i}^{\dagger}$, previously
defined in Sec. III, and the following additional ones

$\begin{split}&\hat{Y}_{\alpha}^{i}=\frac{1}{i}[\hat{T}_{\alpha,+}^{i}-(\hat{T}_{\alpha,+}^{i})^{\dagger}],\\\
&\hat{X}_{\alpha}^{i}=\hat{T}_{\alpha,-}^{i}+(\hat{T}_{\alpha,-}^{i})^{\dagger}.\end{split}$
(19)

Referring to the EoM (12) for the density operator, written in the classical
basis $\\{\otimes_{i=1}^{N}\ket{k_{i}}\\}$, the dynamics of the operators
$\hat{\Omega}_{i}$, $\hat{\Omega}_{i}^{\dagger}$ corresponds to the evolution
of the populations, whereas the operators $\hat{X}_{\alpha}^{i}$ and
$\hat{Y}_{\alpha}^{i}$ capture the off-diagonal elements, i.e. the coherences.
The EoM for the operator $\hat{\Omega}_{i}$ reads

$\dot{\hat{\Omega}}_{i}=-\frac{\gamma}{2}\hat{\Omega}_{i}-\gamma[\hat{f}_{1}^{i}\hat{\Omega}_{i}+\hat{f}_{2}^{i}\hat{\Omega}_{i}^{\dagger}+\hat{f}_{3}^{i}]-\lambda\sum_{\alpha=1}^{3}\omega^{\alpha-1}(\omega-1)\hat{Y}_{\alpha}^{i},$
(20)

where $\hat{f}^{i}_{\alpha}$, $i=1,...,N$ and $\alpha=1,2,3$ are given as
combinations of the operators $\hat{\Gamma}^{i}_{\alpha,s}$ and defined by
Eqs. (38)-(40). The EoMs for the operators $\hat{X}_{\alpha}^{i}$ and
$\hat{Y}_{\alpha}^{i}$ are more involved, given that such operators do not
commute with the rates $\hat{\Gamma}^{j}_{\alpha^{\prime},s}$, and crossing
terms with $j\neq i$ are in principle present in the corresponding Lindblad
equation. However, such additional terms are expected to scale as $O(p/N)$
[41], and can thus be neglected in the thermodynamic limit (see also Appendix
B) that we enforce in the following. This yields

$\begin{split}&\dot{\hat{X}}_{\alpha}^{i}=-\frac{\gamma}{2}h_{i}^{-}(\alpha)\hat{X}_{\alpha}^{i}-\lambda(\hat{Y}_{\alpha-2}^{i}-\hat{Y}_{\alpha}^{i}),\\\
&\dot{\hat{Y}}_{\alpha}^{i}=-\frac{\gamma}{2}h_{i}^{+}(\alpha)\hat{Y}_{\alpha}^{i}+\lambda[\hat{X}_{\alpha+2}^{i}-\hat{X}_{\alpha}^{i}+2(\hat{P}_{\alpha+1}-\hat{P}_{\alpha})],\end{split}$
(21)

where
$h^{s}_{i}(\alpha)=\sum_{s^{\prime}=\pm}(\Gamma_{\alpha-s^{\prime},s^{\prime}}^{i,2}+\Gamma_{\alpha-s^{\prime}+s,s^{\prime}}^{i,2})$.
The operators $\hat{P}_{\eta}^{i}$ are projectors on the classical basis
$\\{\ket{(\eta-1)_{i}}\\}$ of the $i$-th Potts spin, and are given by Eq. (36)
in terms of the operators $\hat{\Omega}_{i}$, $\hat{\Omega}_{i}^{\dagger}$.

###  IV.2 Equations of motion

We are interested in investigating the behavior of the network by analyzing
the expectation values of the overlap operator (18). The latter will allow us
to identify the presence or the absence of pattern retrieval. In the
considered case, $q=3$, this is given by

$\displaystyle\hat{m}^{\mu}=\frac{1}{2N}\sum_{i=1}^{N}(\xi^{*,\mu}_{i}\hat{\Omega}_{i}+\xi^{\mu}_{i}\hat{\Omega}_{i}^{\dagger}).$
(22)

In the presence of quantum coherences we have to consider four further
collective observables, which dynamically couple to $\hat{m}^{\mu}$. These are
constructed from the single spin operators (LABEL:eq_X_Y):

$\begin{split}&\hat{x}^{\mu}=\frac{1}{6N}\sum_{\alpha=1}^{3}\sum_{i=1}^{N}\omega^{\alpha}\xi_{i}^{\mu}\hat{X}_{\alpha}^{i}+\mathrm{h.c.},\\\
&\hat{\bar{x}}^{\mu}=i\left[\frac{1}{6N}\sum_{\alpha=1}^{3}\sum_{i=1}^{N}\omega^{-\alpha}\xi_{i}^{*,\mu}\hat{X}_{\alpha}^{i}-\mathrm{h.c.}\right],\\\
&\hat{y}^{\mu}=\frac{1}{6N}\sum_{\alpha=1}^{3}\sum_{i=1}^{N}\omega^{\alpha}\xi_{i}^{\mu}\hat{Y}_{\alpha}^{i}+\mathrm{h.c.},\\\
&\hat{\bar{y}}^{\mu}=i\left[\frac{1}{6N}\sum_{\alpha=1}^{3}\sum_{i=1}^{N}\omega^{-\alpha}\xi_{i}^{*,\mu}\hat{Y}_{\alpha}^{i}-\mathrm{h.c.}\right].\end{split}$
(23)

These collective operators are constructed such that their EoMs, together with
that of $\hat{m}^{\mu}$, form a closed set, provided that the following
approximations are performed (see Appendix B for details): $(i)$ We neglect
correlations amongst operators, which amounts to a mean field approximation.
This means that, given a collective operator $\hat{O}^{\mu}$, we derive the
EoMs for the quantity $O^{\mu}\equiv\braket{\hat{O}^{\mu}}$, using
$\braket{\hat{O}^{\mu}\hat{O}^{\nu}}\approx\braket{\hat{O}^{\mu}}\braket{\hat{O}^{\nu}}$.
Given the all-to-all coupling of the Potts spin model this approximation is
expected to yield exact results in the thermodynamic limit [59]. $(ii)$ Since
the stationary state is built upon the classical stochastic dynamics, we
expect the dynamics toward stationarity to be mostly determined by the
properties of the overlaps $m^{\mu}$, and therefore we replace
$\braket{\hat{\Omega}_{i}}\approx\xi_{\mu}^{i}m^{\mu}.$ $(iii)$ We assume that
the system exhibits the property of self averaging, i.e. that, given a generic
function $g(\xi_{i}^{\mu})$ of the patterns, for large $N$ we can perform the
substitution

$\frac{1}{N}\sum_{i}g(\xi_{i}^{\mu})\rightarrow\sum_{\\{\xi^{\mu}\\}}\mathbf{P}(\xi)g(\xi^{\mu})\equiv\braket{\braket{g(\xi^{\mu})}},$

where $\mathbf{P}(\xi)$ identifies the pattern probability distribution,
chosen such that $\mathbf{P}(\xi_{i}^{\mu}=\omega^{k_{i}^{\mu}})=1/q$.

Under these approximations, the complete set of EoMs reads

$\begin{split}\dot{m}^{\mu}=&-\frac{\gamma}{2}m^{\mu}\left\\{1+2\braket{\braket{\mathrm{Re}[\xi^{\mu}f_{2}(\\{\xi\\})]}}\right\\}\\\
&-\gamma\braket{\braket{\mathrm{Re}[\xi^{*,\mu}f_{3}(\\{\xi\\})]}}\\\
&-3\lambda\left\\{[1-\mathrm{Re}(\omega)]y^{\mu}+\mathrm{Im}(\omega)\bar{y}^{\mu}\right\\},\\\
\dot{x}^{\mu}=&-\frac{\gamma}{3}x^{\mu}-\lambda[(\mathrm{Re}(\omega^{2})-1)y^{\mu}+\mathrm{Im}(\omega^{2})\bar{y}^{\mu}],\\\
\dot{\bar{x}}^{\mu}=&-\frac{\gamma}{3}\bar{x}^{\mu}-\lambda[(\mathrm{Re}(\omega^{2})-1)\bar{y}^{\mu}-\mathrm{Im}(\omega^{2})y^{\mu}]\\\
\dot{y}^{\mu}=&-\frac{\gamma}{3}y^{\mu}+\lambda[(\mathrm{Re}(\omega^{2})-1)x^{\mu}-\mathrm{Im}(\omega^{2})\bar{x}^{\mu}+m^{\mu}],\\\
\dot{\bar{y}}^{\mu}=&-\frac{\gamma}{3}\bar{y}^{\mu}+\lambda[(\mathrm{Re}(\omega^{2})-1)\bar{x}^{\mu}+\mathrm{Im}(\omega^{2})x^{\mu}+\frac{m^{\mu}}{\sqrt{3}}].\end{split}$
(24)

Note that, by getting the EoMs for the collective operators (23) by exploiting
Eqs. (21), we first replace the damping term $h_{i}^{s}(\alpha)$ by its
average over the pattern disorder distribution,
$h_{i}^{s}(\alpha)\rightarrow\braket{\braket{h^{s}(\alpha)}}=\frac{2}{3}$ (see
Appendix B for more details).

Before going ahead with the analysis of the quantum qPHN as described by Eqs.
(24), it is worth commenting upon their structure. First of all, all the $5p$
differential equations (24) display a damping term, which will enable the
system to reach a stationary state that we will characterize in the next
section. Moreover, we can see that the EoMs of the overlaps and the ones of
the coherences are coupled via terms which depend on the coherent control
parameter $\lambda$. Such a feature already suggests that the addition of the
quantum Hamiltonian (16) to the classical qPHN may play a non-trivial role
with respect to the retrieval properties of the network, as we are going to
further highlight in the remaining part of this work.

##  V Results

In this section we investigate the long-time behavior of a quantum qPHN which
evolves under Eqs. (24), hereafter setting $\gamma=1$. We first focus on the
case of one stored pattern, $p=1$. As pointed out in Sec. II, a gauge
transformation allows one to describe such single memory qPHN by means of the
standard Potts model. For this reason, we eventually focus on the case $p=2$
(two memories), which displays a more interesting phenomenology.

###  V.1 Purely dissipative case

Let us first consider the case $\lambda=0$, i.e. the absence of coherent
effects, in order to establish a baseline and to investigate the role of
temperature. The results, displayed in Fig. 2, show the value of the overlap
at sufficiently long times. By the latter expression we mean that our results
are obtained from numerically integrating the EoMs (24), up to times such that
the overlap converges to the same value, named $m_{\mathrm{stat}}$. We refer
to this type of solution as long-time one, and it reasonably captures the
behaviour of the system sufficiently close to stationarity. In Fig. 2 the
stationary solution is reached starting form an initial configuration that can
partially overlap with the pattern, i.e. we select $m(0)\gtrsim 0$. Here we
see that at low temperatures $m_{\mathrm{stat}}\approx 1$, signaling the
retrieval of the pattern, whereas at high temperatures the overlap vanishes,
$m_{\mathrm{stat}}\approx 0$, characterizing a paramagnetic phase. The
transition between the two phases occurs discontinuously, consistently with a
first order phase transition, as already highlighted in Fig. 1(c).
Furthermore, due to the existence of a coexistence region — typical for first
order transitions — in the interval $T\in[T_{m},T_{M}]$, the crossing is
characterized by hysteresis, which manifests when starting from different
initial conditions. Indeed, the red, thicker line in Fig. 2 has been obtained
as long-time solutions of the EoMs (24) by choosing the initial condition near
the stored pattern, $m(0)\approx 1$. In contrast, for the blue, thinner line
we set $m(0)\approx 0$.

The dependence on the choice of the initial conditions within the coexistence
region is further highlighted in the insets. At relatively low and high
temperatures, with respect to the critical one that emerges in the classical
qPHN, ($T_{c}$ in Sec. II), the results are consistent in both cases. To show
an instance, we set $T=1$ and $T=3$ in the insets, highlighting retrieval of
the pattern and non-retrieval, respectively. When setting the temperature to a
value inside the coexistence region, e.g. $T=2.14$, the system shows retrieval
of the pattern if the initial condition is close to the solution corresponding
to the stored pattern (top right inset). Conversely, it displays non-retrieval
if the initial condition is chosen so that the overlap between the network
configuration and the pattern is almost vanishing (bottom left inset).

Figure 2: Phase diagram of the qPHN in the absence of quantum effects. We set
$q=3$, $\lambda=0$, $p=1$, $\gamma=1$. Stationary value of the overlap
$m_{\mathrm{stat}}$ as a function of the temperature $T$. The values of the
overlap are derived as long-time solutions of the EoMs, as explained in the
text. Initial conditions have been chosen $0<m(0)<0.2$ for the red, thicker
line, and $m(0)>0.5$ for the blue, thinner line. In both cases a paramagnetic
phase, with a vanishing overlap, occurs at high temperature, and a retrieval
phase, characterized by $m_{\mathrm{stat}}\approx 1$, takes place at low
temperatures. The transition between the two phases occurs discontinuously,
consistently with a first order phase transition. Furthermore, it is
characterized by hysteresis: the transition from retrieval to paramagnetic
occurs at lower temperatures if the initial condition is near the paramagnetic
solution (blue, thinner line), and at higher temperatures if the initial
condition is near the retrieval solution (red, thicker line). This feature is
additionally stressed in the insets, which show the evolution of the overlap
for different initial conditions. We set $m(0)>0.5$ for the top-right inset,
and $0<m(0)<0.2$ for the bottom-left one. In both the insets, at the
temperature $T=3$ (black dashed line) the overlap approaches the paramagnetic
solution, and at $T=1$ (black solid line) the overlap approaches the retrieval
solution. At the temperature $T=2.14$, the stationary solution depends on the
initial condition, due to the coexistence of the retrieval and paramagnetic
solutions, namely they are both stable fixed points of the dynamics.

###  V.2 Quantum effects

In the following we consider the long-time behaviour in the presence of
quantum effects, i.e. $\lambda\neq 0$. The EoMs (24) are numerically
integrated at sufficiently long times, as explained at the beginning of Sec.
V.1. They thus provide us what we refer to as stationary solutions. In this
case quantum fluctuations compete with thermal fluctuations, which gives rise
to a rich phase diagram as shown in Fig. 3 for the case $p=1$. The coherent
control parameter $\lambda$ is here varied in the region $[0,1]$, which we
refer to as ”low” $\lambda$ regime. We will eventually consider the ”high”
$\lambda$ case, ($\lambda>1$), that gives rise to some additional physical
effects. Referring to Fig. 3(a), at high temperature (light-blue region), we
find one stationary solution of the EoMs. This is a paramagnetic solution,
characterized by $m_{p}=y_{p}=\bar{y}_{p}=x_{p}=\bar{x}_{p}=0$. At low
temperature (gray region), stationary solutions with finite overlap, $m\neq
0$, appear. We refer to this regime as retrieval phase. The two distinct
phases are separated by a coexistence region (dark blue strip) of paramagnetic
solution and positive retrieval solutions, $m>0$.

Figure 3: Phase diagram of the qPHN in the presence of thermal and quantum
fluctuations. We set $\gamma=1$, $q=3$, and $p=1$. (a) Retrieval and
paramagnetic phase. At high temperature (light blue region) there is one
stationary solution of Eqs. (24) with vanishing order parameters,
corresponding to a paramagnetic phase. At low temperature (grey region) the
stationarity solutions are characterized by finite overlap values, both
positive, $m>0$, and negative, $m<0$. This phase is referred to as retrieval
phase. The dark blue region within dashed black lines divides paramagnetic
phase and retrieval one, and it corresponds to a coexistence region. Here,
both paramagnetic solution, $m=0$, and positive retrieval one, $m>0$, are
present. (b, c) Stationary value of the overlap when the dynamics is
initialized in the BA of the positive solution, and of the negative one,
respectively. In the former case, the transition from the paramagnetic phase
to the retrieval one is discontinuous, whereas in the latter case, such a
transition is a continuous one.

Also in the presence of quantum effects, the retrieval phase admits both
positive and negative solutions, $m_{\mathrm{stat}}>0$, and
$m_{\mathrm{stat}}<0$, respectively. This is highlighted in Fig. 3(b,c), where
the transition from the paramagnetic to the retrieval phase is shown for
varying $T$ and $\lambda$. As commented in Sec. II, negative solutions
correspond to network configurations that are uncorrelated with respect to the
pattern, at variance with the positive ones. These distinct fixed points are
characterized by different basins of attraction (BAs). The latter determines
which stationary solution will be reached, once the system is initialized
within its BA. Thus, the value taken by the overlap towards stationarity,
either signaling pattern retrieval ($m_{\mathrm{stat}}>0$) or not
($m_{\mathrm{stat}}<0$), depends on the chosen initial conditions. It is worth
noticing that the transition between the paramagnetic solution and the
retrieval one occurs discontinuously [see panel (b)], whereas it occurs
continuously when negative solutions are reached [see panel (c)]. Such a
phenomenology holds also in the purely dissipative case $\lambda=0$, where
negative solutions vanish at $T=2$. Thus, consistently with our construction,
the quantum qPHN reproduces the the classical phenomenology in absence of the
Hamiltonian term (16). Finally, referring to the retrieval phase in panel (b),
as the control parameter $\lambda$ is increased, the actual capacity of the
network to store the pattern is diminished, as the overlap reduces from
$m\approx 0.8-1.$ to $0.4$.

From the phase diagram shown in Fig. 3, we see that the addition of the
quantum coherent term to the classical dynamics still enables the network to
retrieve patterns, but it diminishes its capability to fully recover patterns.
Beyond that, quantum effects can, however, also lead to a qualitative change
of the network’s dynamics: for a certain range of parameters the overlap and
the coherence variables display self-sustained oscillations at long times, as
shown in Figs. 4 and 5. Let us now analyze in more detail these periodic long-
time solutions, which occur in the ”high”-$\lambda$ parameter regime and are
referred to as limit cycle (LC) solutions. First we analyze the latter in the
case $p=1$, and eventually we set $p=2$. The former case is highlighted in
Fig. 4. It shows the standard deviation $\sigma_{m}$ of the overlap $m(t_{i})$
with respect to the time average of $m(t)$, say $\bar{m}$, at sufficiently
long times. The values $m(t_{i})$ are obtained at times $t_{i}\in I$, with $I$
an arbitrary time interval chosen at long times. Thus the quantity displayed
reads

$\sigma_{m}=\sqrt{\frac{1}{N_{I}}\sum_{i\in I}[m(t_{i})-\bar{m}]^{2}},$ (25)

being $N_{I}$ the number of points $t_{i}\in I$. We set $N_{I}=2000$, and
$I=[9\times 10^{3},10^{4}]$. The light-blue region where $\sigma_{m}$ is
finite corresponds to the LC phase. In the inset, we display the corresponding
oscillating behavior of the overlap. As a matter of fact, the amplitude of the
oscillations varies on a considerably small scale, $\sigma_{m}\lesssim
10^{-4}$.

As further detailed in App. B, we also explore the LC phase via a linear
stability analysis. To this end, we linearize Eqs. (24) with respect to their
stationary solutions, and study the eigenvalues of the corresponding Jacobian
matrix. In this analysis we fix $\lambda$ to values well inside the region
displaying a LC phase, and decrease the temperature $T$ starting from
$T$-values corresponding to the paramagnetic phase. With changing $T$ and
while approaching the LC phase, a complex conjugate pair of eigenvalues
crosses the complex plane imaginary axes. In non-linear systems such a
phenomenology is referred to as Hopf bifurcation [60]. It identifies a
switching in the stability of the system, often giving rise to periodic
solutions. For instance, Hopf bifurcations and LCs occur in the known Lotka-
Volterra dynamical system [61], which is frequently used to model the dynamics
of predator-prey interactions in biological systems [62]. In two-dimensional
dynamical systems, periodic long-time solutions characterizing LCs are
displayed in phase space as isolated orbits [60]. In higher dimensional
dynamical systems, as it is our case, one can still highlight isolated orbits
by considering two out of the total number of dynamical variables involved, as
we are going to show for the case of two stored patterns $p=2$.

Figure 4: One memory limit cycle. We set $\gamma=1$, $q=3$, $p=1$. Standard
deviation $\sigma_{m}$ of the overlap $m(t)$ at long times, with respect to
the fixed point solutions, as defined by Eq. (25). Initial conditions are
chosen such that $m(0)\approx 0$. The dashed lines identify the separation
among paramagnetic phase, coexistence region and retrieval phase. In the blue
region $\sigma_{m}$ is finite, and such a portion of the phase diagram
corresponds to a limit cycle phase. The inset shows the evolution of the
overlap $m(t)$ at long times, highlighting an oscillatory behavior. We set
$\lambda=4.0$ , $T=0.86$.

Fig. 5 (b), left-hand side panel, illustrates an instance of the isolated
orbit that characterizes the LCs in our two-memories model. Here we show the
parametric plot of the overlap $m^{1}(t)$ and its corresponding coherence
variable $y^{1}(t)$. Trajectories surrounding the isolated orbit constitute
the flux diagram of the vector field $(\dot{y}^{1},\dot{m}^{1})$, with the
value of its norm color-coded and increasing from purple to yellow. We also
set the parameter value inside the LC phase, e.g. $T=0.8$ and $\lambda=4.5$.
Different choices of parameters lead to a similar phenomenology.

The parameter region where the LC phase takes place is displayed in Fig. 5(a).
It shows the standard deviation $\sigma_{m}$ defined by Eq. (25), averaged
with respect to the two overlaps. At this point, it is worth considering the
comparison with the LC phase of the single memory case shown in Fig. 4. We can
see that the LC phase for the case $p=2$ is $(i)$ extended into the low
temperature region, and $(ii)$ characterized by a larger amplitude of the
oscillations occurring at long times. The latter feature can be further
observed, e.g., from comparing the right-hand side of Fig. 5(b), which
displays the time evolution of $m^{1}(t)$ and $y^{1}(t)$, and the inset of
Fig. 4. Let us then try to provide a qualitative explanation of the features
$(i)$ and $(ii)$. To this end, we contrast our results with those obtained for
the quantum HNN analyzed in Ref. [41]. There, similarly to our model, as a
result of the competition between dissipation and coherent dynamics, a LC
phase emerges. However, it turns out to be independent of the number of stored
patterns. In other words, in this case there is no distinction between the
case $p=1$ and $p=2$, at variance with our quantum qPHN model. A potential
relevant difference between the HNN and the qPHN is the symmetry of the model.
As already mentioned in Sec. II, the former is characterized by a
$\mathbb{Z}_{2}$ symmetry. Direct manifestation of the latter is that once the
pattern $\\{\xi_{1}^{\mu},...,\xi_{N}^{\mu}\\}$ is stored as a memory, so is
the opposite one, $\\{-\xi_{1}^{\mu},...,-\xi_{N}^{\mu}\\}$. Contrarily, the
$q\geq 2$ qPHN does not display the same symmetry, and the notion of
”opposite” pattern is meaningless. Such a property might be the origin of the
dissimilarity between the two models: when fixing $p=1$ and selecting the
proper parameter regime, in the HNN the pattern and its opposite one can give
rise to the LC phase. Instead, the $p=1$ qPHN cannot give rise to the same
type of LC solutions, and it is only when considering at least $p=2$, that our
model features two overlaps variables that can give rise to a LC phase.

The behavior of the overlaps for a $p=2$ qPHN, and within the LC phase, is
shown in Fig. 5(c). As visible in the right-hand side panel, close to the
paramagnetic phase, there exists a narrow region where only one of the two
patterns features a LC solution. Notably, the amplitude of the corresponding
oscillations decreases with respect to the ones displayed in the left-hand
side panel figure. Such a behavior takes place within the coexistence region
of paramagnetic and retrieval fixed point solutions (the region within the
dashed black lines). It is worth stressing that this region coincides also
with the LC phase of the $p=1$ case (see Fig. 4). Moreover, also in this case
the periodic solution is characterized by a relatively small amplitude of
oscillations. These features hint that the coexistence of fixed point
solutions may play a role in determining the properties of LC solutions
occurring in such a parameter regime. In the bulk of the LC phase, as shown in
the left-hand side panel of Fig. 5(c), the two overlaps present out-of-phase,
self-sustained oscillations in the long-time limit. Here the amount of
retrieval is relatively small, i.e. it is at most $m\approx 0.4$. However, as
a consequence of the out-of-phase oscillations between the two overlaps, the
storage of either one of the two patterns is in principle possible. Indeed,
the periodic switching between positive values of the two overlaps, with
$m\lesssim 0.4$, implies that it is in principle possible to alternatively
retrieve either one the two patterns. This phenomenology allows us to identify
the bulk of the LC phase as a new type of ”quantum” retrieval phase, emerging
as a feature of the quantum formulation of the model.

Figure 5: Two memory limit cycle phase. We set $\gamma=1$, $q=3$, $p=2$. (a)
Standard deviation $\sigma_{m}$ of the overlap $m^{\mu}(t)$ at long times,
with respect to the stationary solution. The dashed lines identify the
separation among parmagnetic phase, coexistence region and retrieval phase. In
the blue region $\sigma_{m}$ is finite, corresponding to a limit cycle phase.
(b) The right-hand panel shows the time evolution of the overlap $m^{1}(t)$
and the coherence variable $y^{1}(t)$ at long times. In the left-hand panel
the corresponding parametric plot and flux diagram are displayed, highlighting
the closed orbit that characterizes the limit cycle solution. Colors for the
flux diagram identify the norm of the corresponding vector field, which
increases from purple to yellow. We set $T=0.8$ and $\lambda=4.5$. (c) Time
evolution of the overlaps $m^{1}(t)$ and $m^{2}(t)$ as specified in the
legend. We set $\lambda=3.0$ and $T=0.99$ for the right-hand panel, where the
self-sustained oscillations at long times appear for only one of the two
overlap. For the left-hand panel it is $\lambda=4.5$ and $T=0.6$, and both the
overlaps display a limit cycle solution.

##  VI Conclusions and Outlook

In this work, we have established a method for generalizing associative memory
NNs via open quantum systems dynamics, which goes beyond the previously
studied case of the HNN [41]. Starting point of our approach is the classical
qPHN, whose non-equilibrium dynamics has been formulated in terms of a
Markovian master equation of Lindblad form. This formulation allows to
incorporate quantum fluctuations through a Hamiltonian term whose strength is
given by a control parameter. By exploiting MF techniques, we construct the
phase diagram for a many-body systems made of $N$ q-Potts spins. This shows
that memory retrieval can be achieved for certain temperatures and values of
this control parameter. Furthermore, the competition between quantum coherent
dynamics and thermal fluctuations has shown to give rise to a new non-
equilibrium quantum phase. The latter features LC solutions, which display
persistent oscillations at long times. We analyze the LC phase for a single
memory model and a two memories one, finding a different phenomenology in the
two cases. The parameter regime where the LC phase occurs shrinks in the
single memory case with respect to the two memories one. In the latter, both
the overlaps sustain persistent, out-of-phase oscillations at long times which
display a larger amplitude with respect to the one of the single memory case.

The analysis and characterization of the phase diagram done in this work
represent initial steps in order to eventually explore more sophisticated
questions. These concern, for example, the issue of identifying the storage
capacity $\alpha_{c}$, i.e. how many patterns can be stored in a network with
$N$ constituents. This aspect has been widely investigated in the classical
realm. For instance, it is known that it is $\alpha_{c}=0.14$ for a classical
HNN [46]. For the case of a classical qPHN, the storage capacity depends on
the number $q$ of levels of Potts spins, being $\alpha_{c}=0.41$ if $q=3$ [52,
56, 53]. The main question here is understanding whether quantum formulations
allow for a larger storage capacity than the classical models. One potential
way to tackle such a problem has been pointed out in the recent theoretical
work of Ref. [48]. Here, the authors consider a generic set of attractive
quantum maps, evolving input quantum states towards stationary ones, labeled
as retrieval states. The size of such a stationary manifold gives access to
the storage capacity, which is found to exponentially outperform the classical
counterpart. Our work lays the ground to tackle the question of storage
capacity for the concrete and novel class of qPHNs introduced in our work, and
to also investigate how physically realizable instances of such quantum neural
networks can be engineered in state-of-the-art quantum hardware.

The presence of the novel types of retrieval phases identified in our work,
which stem from the quantum description of the NNs, opens interesting new
directions in relation to more quantum-information oriented questions: A first
interesting question is how to formulate non-classical patterns, i.e. specific
input quantum states, in order to store and access them. Here, the
investigation of the speed of retrieval [42] with respect to the classical
counterpart represents an important aspect. Furthermore, it would be worth
exploring the link of quantum formulations of associative memories as realized
by our quantum neural networks, with the active research field of quantum
error correction. Indeed, one could investigate the capabilities of quantum
generalization of Hopfield-type NNs in protecting and correcting faulty
quantum states.

## Acknowledgments

We gratefully acknowledge useful discussions with L. Bödeker. EF and MM
acknowledge support by the ERC Starting Grant QNets Grant Number 804247. IL
acknowledges support from the “Wissenschaftler Rückkehrprogramm GSO/CZS” of
the Carl-Zeiss-Stiftung and the German Scholars Organization e.V. EF, MM and
IL are grateful for funding from the Deutsche Forschungsgemeinschaft through
Grant No. 449905436.

##  Appendix A Equation of motion for the purely dissipative case

In this section we derive the EoMs for the purely dissipative case
($\lambda=0$) of the qPHN, showing that one can recover the HNN model for
$q=2$.

As a first step, we evaluate the EoMs of the Potts spin operator
$\hat{\Omega}_{i}^{\alpha}$, $i=0,...,N$. To this end, we exploit Eqs.
(9)-(11) for getting the following relations for the operators
$\hat{T}^{i}_{\eta,\pm}$,

$\displaystyle\hat{T}^{i}_{\eta,\pm}\hat{\Omega}_{i}^{\alpha}=\omega^{\pm\alpha}\hat{\Omega}_{i}^{\alpha}\hat{T}^{i}_{\eta,\pm},$
(26)
$\displaystyle\hat{T}^{i}_{\eta,\pm}\hat{T}^{i}_{\eta,\mp}=\hat{P}^{i}_{\eta},$
(27)
$\displaystyle\sum_{\eta=1}^{q}\hat{T}^{i}_{\eta,\pm}\hat{T}^{i}_{\eta,\mp}=\mathbf{1}_{q},$
(28)

where $\hat{P}_{\eta}=\ket{\eta-1}\bra{\eta-1}$. The equations of motion read

$\dot{\hat{\Omega}}_{i}^{\alpha}=\hat{\Omega}_{i}^{\alpha}\left\\{-\gamma[(1-\omega^{-\alpha})\hat{W}_{i,+}+(1-\omega^{\alpha})\hat{W}_{i,-}]\right\\},$
(29)

where we have defined

$\displaystyle\gamma\hat{W}_{i,s}$ $\displaystyle\equiv$
$\displaystyle\sum_{\alpha=1}^{q}\hat{L}_{\alpha,s}^{i,\dagger}\hat{L}_{\alpha,s}^{i}=\gamma\sum_{\alpha=1}^{q}\Gamma_{\alpha,s}^{i,2}(\hat{T}_{\alpha,s}^{i})^{\dagger}\hat{T}_{\alpha,s}^{i}$
$\displaystyle=$
$\displaystyle\gamma\sum_{\alpha=1}^{q}\Gamma_{\alpha,s}^{i,2}\hat{P}^{i}_{\alpha+s}=\gamma\sum_{\alpha=1}^{q}\Gamma_{\alpha+s,s}^{i,2}\hat{P}_{\alpha}^{i}.$

It is straightforward to express the Potts spin operators
$\hat{\Omega}_{i}^{\alpha}$ in terms of the projection operators as
$\hat{\Omega}_{i}^{\alpha}=\sum_{\eta=1}^{q}\omega^{\alpha(\eta-1)}\hat{P}_{\eta}^{\alpha}$.
This relation can in general be inverted, yielding

$\hat{P}_{\eta}^{i}=\sum_{a=1}^{q}g_{\eta a}\hat{\Omega}_{i}^{a}.$ (31)

with some coefficients $g_{\eta\alpha}$ that can be derived for each choice of
$q$. The Potts spin operators evolve according to the equation

$\begin{split}&\dot{\hat{\Omega}}_{i}^{\alpha}=-\gamma\sum_{a=1}^{q}f_{\alpha,a}^{i}\hat{\Omega}_{i}^{a},\\\
&f_{\alpha,a}^{i}=\sum_{\eta=1}^{q}g_{\eta
a}\omega^{(\eta-1)\alpha}[(1-\omega^{-\alpha})\Gamma_{\eta-1,+}^{i,2}+(1-\omega^{\alpha})\Gamma_{\eta+1,-}^{i,2})].\end{split}$
(32)

####  A.0.1 The case q=2 - Hopfield neural network

Taking $q=2$, the operator $\hat{\Omega}_{i}^{\alpha}$ becomes the spin $1/2$
Pauli operator $\hat{\sigma}^{z}_{i}$. Specializing Eq. (29) to this case
yields

$\dot{\hat{\sigma}}^{z}_{i}=-2\gamma\hat{\sigma}^{z}_{i}(\hat{W}_{i,+}+\hat{W}_{i,-}),$
(33)

which can be further simplified as

$\displaystyle\hat{W}_{i,+}+\hat{W}_{i,-}$ $\displaystyle=$
$\displaystyle\sum_{\alpha=1}^{2}(\Gamma_{\alpha-1,+}^{i,2}+\Gamma_{\alpha+1,-}^{i,2})\hat{P}_{\alpha}$
$\displaystyle=$
$\displaystyle\frac{1}{Z_{i}}\sum_{\alpha=1}^{2}(e^{\beta\Delta\hat{E}_{\alpha-1,+}^{i}}+e^{\beta\Delta\hat{E}_{\alpha+1,-}^{i}})\hat{P}_{\alpha}$
$\displaystyle=$
$\displaystyle\frac{1}{Z_{i}}(\Gamma_{2,+}^{i,2}+\Gamma_{2,-}^{i,2})\hat{P}_{1}+(\Gamma_{1,+}^{i,2}+\Gamma_{1,-}^{i,2})\hat{P}_{2},$

where $\hat{P}_{1,2}=(1\pm\hat{\sigma}_{i}^{z})/2$, and
$\Gamma_{2,+}^{i,2}=\Gamma_{2,-}^{i,2}=e^{\beta\Delta\hat{E}_{H}}/Z$,
$\Gamma_{1,+}^{i,2}=\Gamma_{1,-}^{i,2}=e^{-\beta\Delta\hat{E}_{H}}/Z$, and
$Z=2(e^{-\beta\Delta\hat{E}_{H}}+e^{\beta\Delta\hat{E}_{H}})$. Thus, we get
the equation for the Hopfield NN [41],

$\dot{\hat{\sigma}}_{i}^{z}=-\gamma\hat{\sigma}_{i}^{z}+\gamma\tanh{(\beta\Delta\hat{E}_{H})}$
(35)

with $\Delta\hat{E}_{H}=\frac{1}{N}\sum_{j\neq i}J_{ij}\hat{\sigma}_{j}^{z}$.

##  Appendix B Equation of motion for the q=3 Quantum Potts Neural Network

This appendix contains details on how to derive Eqs. (24) of the main text,
the quantum generalization of the $q=3$ Potts-Hopfield neural network.

We proceed by specifying Eq. (17) for the operators $\hat{\Omega}_{i}$ and
$\hat{\Omega}_{i}^{2}=\hat{\Omega}_{i}^{\dagger}=\hat{\Omega}_{i}^{-1}$. To
this end, it is useful to write the projectors
$\hat{P}_{\eta}=\ket{\eta-1}\bra{\eta-1}$, $\eta=1,2,3$ in terms of
$\hat{\Omega}_{i}$, $\hat{\Omega}_{i}^{\dagger}$, and the identity operator
$\mathbf{1}_{3}$ as follows

$\begin{split}&\hat{P}^{i}_{1}=\frac{1}{(\omega-1)^{2}}[-\omega\hat{\Omega}_{i}+(1+\omega^{2})\hat{\Omega}^{\dagger}_{i}-\omega\mathbf{1}_{3}],\\\
&\hat{P}^{i}_{2}=\frac{1}{(\omega-1)^{2}}[-\hat{\Omega}_{i}+(1+\omega)\hat{\Omega}^{\dagger}_{i}-\omega\mathbf{1}_{3}],\\\
&\hat{P}^{i}_{3}=\frac{1}{(\omega-1)^{2}}[(1+\omega)\hat{\Omega}_{i}-\hat{\Omega}^{\dagger}_{i}-\omega\mathbf{1}_{3}].\\\
\end{split}$ (36)

Let us first derive the EoMs for Potts spin operators $\hat{\Omega}_{i}$,
$\hat{\Omega}_{i}^{\dagger}$ without the Hamiltonian term, i.e. taking
$\lambda=0$. The evolution is governed by

$\begin{split}&\dot{\hat{\Omega}}_{i}=-\gamma[(\frac{1}{2}+\hat{f}_{1}^{i})\hat{\Omega}_{i}+\hat{f}_{2}^{i}\hat{\Omega}_{i}^{\dagger}+\hat{f}_{3}^{i}],\\\
\end{split}$ (37)

where we have defined

$\hat{f}_{1}^{i}=i\frac{1}{2\sqrt{3}}[\hat{\Gamma}_{1,+}^{i,2}+\hat{\Gamma}_{2,+}^{i,2}+\hat{\Gamma}_{3,+}^{i,2}-(\hat{\Gamma}_{1,-}^{i,2}+\hat{\Gamma}_{2,-}^{i,2}+\hat{\Gamma}_{3,-}^{i,2})],$
(38)
$\begin{split}\hat{f}_{2}^{i}&=\frac{1}{2}[\hat{\Gamma}_{3,+}^{i,2}-\hat{\Gamma}_{3,-}^{i,2}-(\hat{\Gamma}_{2,+}^{i,2}-\hat{\Gamma}_{2,-}^{i,2})]\\\
&+\frac{i}{2\sqrt{3}}[\hat{\Gamma}_{3,+}^{i,2}-\hat{\Gamma}_{3,-}^{i,2}+\hat{\Gamma}_{2,+}^{i,2}-\hat{\Gamma}_{2,-}^{i,2}-2(\hat{\Gamma}_{1,+}^{i,2}-\hat{\Gamma}_{1,-}^{i,2})],\end{split}$
(39)
$\begin{split}\hat{f}_{3}^{i}&=\frac{1}{2}(\hat{\Gamma}_{3,+}^{i,2}+\hat{\Gamma}_{2,-}^{i,2}-\hat{\Gamma}_{1,+}^{i,2}-\hat{\Gamma}_{1,-}^{i,2})\\\
&+\frac{i}{2\sqrt{3}}[\hat{\Gamma}_{3,+}^{i,2}-\hat{\Gamma}_{2,-}^{i,2}+\hat{\Gamma}_{1,+}^{i,2}-\hat{\Gamma}_{1,-}^{i,2}-2(\hat{\Gamma}_{2,+}^{i,2}-\hat{\Gamma}_{3,-}^{i,2})].\\\
\end{split}$ (40)

The rates $\hat{\Gamma}_{\alpha,s}^{i}$ are given by Eq. (13) and depend on
the energy difference

$\begin{split}&\Delta\hat{E}_{\alpha,s}^{i}\\\
&=\sum_{\mu=1}^{p}\hat{m}^{\mu}\left\\{3\cos[\frac{2\pi}{3}(k_{i}^{\mu}-\alpha+1)]\right.\\\
&\left.\pm\sqrt{3}\sin[\frac{2\pi}{3}(k_{i}^{\mu}-\alpha+1)]\right\\}\\\
&=3\sum_{\mu=1}^{p}\hat{m}^{\mu}(\delta_{k_{i}^{\mu},\alpha-1}-\delta_{k_{i}^{\mu}\pm
1,\alpha-1}),\end{split}$ (41)

for $\alpha=1,2,3$.

Let us now include the Hamiltonian term (16), setting $\lambda\neq 0$. This
generates a new term in the EoMs of the Potts operators which reads

$\begin{split}i[\hat{H},\hat{\Omega}_{i}]&=-\lambda\sum_{\alpha=1}^{3}\omega^{\alpha-1}(\omega-1)\hat{Y}_{\alpha}^{i}.\end{split}$
(42)

Here we have used that
$\hat{\Omega}_{i}\hat{T}_{\alpha,+}^{i}=\sum_{\alpha^{\prime}=1}^{3}\omega^{\alpha^{\prime}-1}\hat{T}_{\alpha^{\prime},+}^{i}$
and
$\hat{\Omega}_{i}\hat{T}_{\alpha,-}^{i}=\sum_{\alpha^{\prime}=1}^{3}\omega^{\alpha^{\prime}-1}\omega(\hat{T}_{\alpha^{\prime},+}^{i})^{\dagger}$,
and we have defined the operators

$\hat{Y}_{\alpha}^{i}=\frac{1}{i}[\hat{T}_{\alpha,+}^{i}-(\hat{T}_{\alpha,+}^{i})^{\dagger}].$
(43)

Thus the full EoM of the Potts operators $\hat{\Omega}_{i}$ reads

$\dot{\hat{\Omega}}_{i}=-\frac{\gamma}{2}\hat{\Omega}_{i}-\gamma[\hat{f}_{1}^{i}\hat{\Omega}_{i}+\hat{f}_{2}^{i}\hat{\Omega}^{\dagger}+\hat{f}_{3}^{i}]-\lambda\sum_{\alpha=1}^{3}\omega^{\alpha-1}(\omega-1)\hat{Y}_{\alpha}^{i}.$
(44)

To construct a closed set of EoM, we also need to derive the ones of the
operators $\hat{Y}_{\alpha}^{i}$, and the ones of the operators

$\displaystyle\hat{X}_{\alpha}^{i}=\hat{T}_{\alpha,-}^{i}+(\hat{T}_{\alpha,-}^{i})^{\dagger}.$
(45)

We first derive the EoMs of the operator $\hat{T}_{\alpha,s}^{i}$, which are
more involved in comparison to the EoMs of the Potts spin operator
$\hat{\Omega}_{i}$. This is due to the fact that the operators
$\hat{T}_{\alpha,s}^{i}$ do not commute with the rates
$\hat{\Gamma}_{\alpha^{\prime},s^{\prime}}^{j}$ that appear in the jumps
operators (13). This gives rise to a cross-term, with $j\neq i$, in the
Lindblad equation (17), which reads

$\sum_{\alpha^{\prime},s^{\prime}}\sum_{j\neq
i}L_{\alpha^{\prime},s^{\prime}}^{j,\dagger}\hat{T}_{\alpha,s}^{i}L_{\alpha^{\prime},s^{\prime}}^{j}-\frac{1}{2}\\{L_{\alpha^{\prime},s^{\prime}}^{j,\dagger}L_{\alpha^{\prime},s^{\prime}}^{j},\hat{T}_{\alpha,s}^{i}\\}.$
(46)

For simplifying such an expression, it is sufficient to consider the non-
commuting terms
$\hat{\Gamma}_{\alpha^{\prime},s^{\prime}}^{j}\hat{T}_{\alpha,s}^{i}$ and
exploit the relation (26) for moving all the operators
$\hat{T}_{\alpha,s}^{i}$ to the left, obtaining

$\begin{split}&\hat{\Gamma}_{\alpha^{\prime},s^{\prime}}^{j}(\hat{\Omega}_{1},...,\hat{\Omega}_{i},...,\hat{\Omega}_{N})\hat{T}_{\alpha,s}^{i}\\\
&=\hat{T}_{\alpha,s}^{i}\hat{\Gamma}_{\alpha^{\prime},s^{\prime}}^{j}(\hat{\Omega}_{1},...,\omega^{-s}\hat{\Omega}_{i},...,\hat{\Omega}_{N})=\hat{T}_{\alpha,s}^{i}\hat{\Gamma}_{\alpha^{\prime},s^{\prime}}^{j(i)},\end{split}$
(47)

with
$\hat{\Gamma}_{\alpha^{\prime},s^{\prime}}^{j}(\hat{\Omega}_{1},...,\hat{\Omega}_{i},...,\hat{\Omega}_{N})$
highlighting the dependence of $\Gamma_{\alpha^{\prime},s^{\prime}}^{j}$ on
all the $N$ Potts spin operators, and the compact notation
$\hat{\Gamma}_{\alpha^{\prime},s^{\prime}}^{j}(\hat{\Omega}_{1},...,\omega^{-s}\hat{\Omega}_{i},...,\hat{\Omega}_{N})=\hat{\Gamma}_{\alpha^{\prime},s^{\prime}}^{j(i)}$
has been introduced in the second line. Now we can recognize that
$\hat{\Gamma}_{\alpha^{\prime},s^{\prime}}^{j}$ and
$\hat{\Gamma}_{\alpha^{\prime},s^{\prime}}^{j(i)}$ depend on two
configurations which differ by a single Potts spin (the $i$-th), from which we
expect that
$\hat{\Gamma}_{\alpha^{\prime},s^{\prime}}^{j}\approx\hat{\Gamma}_{\alpha^{\prime},s^{\prime}}^{j(i)}$
up to finite-size corrections which scale as $1/N$. As derived in Ref. [41],
the crossing term (47) in the EoMs for operators $\hat{T}_{\alpha,s}^{i}$
should more precisely scale as $O(p/N)$. It can be thus neglected for a
subextensive number of pattern, $p\ll N$, and in the thermodynamic limit,
leading to the EoM

$\dot{\hat{T}}_{\alpha,s}^{i}=-\frac{\gamma}{2}\hat{T}_{\alpha,s}^{i}h^{s}_{i}(\alpha)+i\lambda[\hat{T}^{i}_{\alpha+2s,-s}+\hat{T}^{i}_{\alpha,-s}+\hat{P}^{i}_{\alpha+s}-\hat{P}_{\alpha}^{i}].$
(48)

Here
$h^{s}_{i}(\alpha)=\sum_{s^{\prime}=\pm}(\Gamma_{\alpha-s^{\prime},s^{\prime}}^{i,2}+\Gamma_{\alpha-s^{\prime}+s,s^{\prime}}^{i,2})$,
and the operator $\hat{P}_{\eta}^{i}$ are defined by Eq. (36) in terms of the
Potts operator.

It follows that

$\begin{split}&\dot{\hat{X}}_{\alpha}^{i}=-\frac{\gamma}{2}h_{i}^{-}(\alpha)\hat{X}_{\alpha}^{i}-\lambda(\hat{Y}_{\alpha-2}^{i}-\hat{Y}_{\alpha}^{i}),\\\
&\dot{\hat{Y}}_{\alpha}^{i}=-\frac{\gamma}{2}h_{i}^{+}(\alpha)\hat{Y}_{\alpha}^{i}+\lambda[\hat{X}_{\alpha+2}^{i}-\hat{X}_{\alpha}^{i}+2(\hat{P}_{\alpha+1}-\hat{P}_{\alpha})].\end{split}$
(49)

It is worth noticing that, although Eqs. (44) and (49) form a closed set of
EoMs, they describe microscopic quantities, i.e. individual Potts spins. Our
interest lies, however, in understanding the behavior of macroscopic,
collective observables, such as the overlap operators. When analyzing such
quantities the number of equations is reduced from $O(N)$ to $O(p)$.

To achieve such a reduction, we begin by writing EoMs for the overlaps

$\begin{split}\dot{\hat{m}}^{\mu}&=-\frac{\gamma}{2}\hat{m}^{\mu}-\frac{\gamma}{2N}\sum_{i=1}^{N}[(\xi_{i}^{*,\mu}f_{1}^{i}+\xi_{i}^{\mu}\hat{f}^{i,\dagger}_{2})\hat{\Omega}_{i}+\mathrm{h.c.}]\\\
&-\frac{\gamma}{N}\sum_{i=1}^{N}\mathrm{Re}[\xi_{i}^{*,\mu}\hat{f}_{3}^{i}]\\\
&-\frac{\lambda}{2N}\sum_{i=1}^{N}\sum_{\alpha=1}^{3}\hat{Y}_{\alpha}^{i}[\xi_{i}^{*,\mu}\omega^{\alpha-1}(\omega-1)+\xi_{i}^{\mu}\omega^{-(\alpha-1)}(\omega^{*}-1)].\end{split}$
(50)

The last term suggests to also consider the collective operators

$\begin{split}&\hat{x}^{\mu}=\frac{1}{6N}\sum_{\alpha=1}^{3}\sum_{i=1}^{N}\omega^{\alpha}\xi_{i}^{\mu}\hat{X}_{\alpha}^{i}+\mathrm{h.c.},\\\
&\hat{\bar{x}}^{\mu}=i\left[\frac{1}{6N}\sum_{\alpha=1}^{3}\sum_{i=1}^{N}\omega^{-\alpha}\xi_{i}^{*,\mu}\hat{X}_{\alpha}^{i}-\mathrm{h.c.}\right],\\\
\end{split}$ (51)

and two additional ones, $\hat{y}^{\mu}$ and $\hat{\bar{y}}^{\mu}$, obtained
by replacing $\hat{X}_{\alpha}^{i}$ with $\hat{Y}_{\alpha}^{i}$. This choice
will allow us to find a closed set of EoMs, when providing the system with
some additional approximations.

To achieve this, we firstly perform a mean field (MF) approximation, replacing
collective operators with their averaged quantities, and neglecting
correlations amongst them. For the classical system, MF is exact in the
thermodynamic limit, and we expect such a treatment to be appropriate also in
our quantum generalization [59, 41]. Thus, we derive EoMs for averaged
quantities $O^{\mu}$ such that $\braket{O}^{\mu}\equiv O^{\mu}$.

Secondly, we expect the system near stationarity to be mostly determined by
the behavior of the overlap. Indeed, we build up the model so that it exactly
reproduces the classical case in the limit $\lambda\rightarrow 0$, and quantum
effects are gradually taken into account as $\lambda$ is increased. Moreover,
the classical case is well described by the dynamics of the overlap variables.
We thus expect that also the description of the quantum dynamics near to the
stationary solutions can be mainly embodied in the dynamics of the overlaps.
Therefore, we assume that approximately
$\braket{\hat{\Omega}_{i}}\approx\xi_{i}^{\mu}m^{\mu}$ holds. This allows us
to treat terms of the type
$\sum_{i=1}^{N}\xi_{i}^{\mu}\hat{\Omega}_{i}\hat{f}_{i}^{\dagger}$ which
appear in the Eq. (50) of the overlap. Furthermore, Eqs. (49) give rise to a
set of EoMs which do not immediately form a closed set for the collective
operators $\hat{x}^{\mu},\hat{y}^{\mu}$. This is due to the presence of terms
of the type
$-\frac{1}{6N}\sum_{\alpha=1}^{3}\omega^{\alpha}\sum_{i=1}^{N}h_{i}^{-}(\alpha)\hat{X}_{\alpha}^{i}\xi_{i}^{\mu}$.
Here, we exploit the average of the damping coefficient $h_{i}^{s}(\alpha)$
with respect to the pattern disorder distribution. An analysis of the dynamics
obtained under this approximation and the dynamics without such an average is
displayed in Fig. 6, panel (a) and (b). Here we show the evolution of the
overlap for a single pattern, $p=1$. The dashed black line is obtained by
solving the exact dynamics for $N=2$, and the solid black line represents the
approximated one. We fixed $\lambda=0.1$ and $\gamma=1$. Panel (a) and (b)
show the result for $T=2.5$ and $T=0.5$, respectively. At large temperatures,
the difference between the two results occurs at early stages of the dynamics
only, whereas at small temperatures, the performed approximations yields a
constant, although relatively small shift that holds even at stationarity.
However, the qualitative behavior of the network remains unaltered.

Figure 6: Approximated and exact dynamics. We set $p=1$, $q=3$, $\gamma=1$,
$\lambda=0.1$. Overlap $m(t)$ as a function of the time $t$ for (a) $T=2.5$,
and (b) $T=0.5$. Dashed black line corresponds to the exact solution, the
solid black line represents the mean field, self-averaged case.

The final step towards obtaining the EoMs for the macroscopic observables
consists in exploiting the so-called self-averaging property, typical of
disordered systems [63]. This means, that in the large $N$ limit and for
$p/N\rightarrow 0$, we use that
$\frac{1}{N}\sum_{i}g(\xi_{i}^{\mu})\rightarrow\sum_{\\{\xi^{\mu}\\}}\mathbf{P}(\xi)g(\xi^{\mu})\equiv\braket{\braket{g(\xi^{\mu})}}$.

After the above approximations the EoMs read

$\begin{split}\dot{m}^{\mu}=&-\frac{\gamma}{2}m^{\mu}\left\\{1+2\braket{\braket{\mathrm{Re}[\xi^{\mu}f_{2}(\\{\xi\\})]}}\right\\}\\\
&-\gamma\braket{\braket{\mathrm{Re}[\xi^{*,\mu}f_{3}(\\{\xi\\})]}}\\\
&-3\lambda\left\\{[1-\mathrm{Re}(\omega)]y^{\mu}+\mathrm{Im}(\omega)\bar{y}^{\mu}\right\\},\\\
\dot{x}^{\mu}=&-\frac{\gamma}{3}x^{\mu}-\lambda[(\mathrm{Re}(\omega^{2})-1)y^{\mu}+\mathrm{Im}(\omega^{2})\bar{y}^{\mu}],\\\
\dot{\bar{x}}^{\mu}=&-\frac{\gamma}{3}\bar{x}^{\mu}-\lambda[(\mathrm{Re}(\omega^{2})-1)\bar{y}^{\mu}-\mathrm{Im}(\omega^{2})y^{\mu}]\\\
\dot{y}^{\mu}=&-\frac{\gamma}{3}y^{\mu}+\lambda[(\mathrm{Re}(\omega^{2})-1)x^{\mu}-\mathrm{Im}(\omega^{2})\bar{x}^{\mu}+m^{\mu}],\\\
\dot{\bar{y}}^{\mu}=&-\frac{\gamma}{3}\bar{y}^{\mu}+\lambda[(\mathrm{Re}(\omega^{2})-1)\bar{x}^{\mu}+\mathrm{Im}(\omega^{2})x^{\mu}+m^{\mu}/\sqrt{3}],\end{split}$
(52)

which are $5p$ differential equations for the collective operators. The
presence of a damping term, proportional to $\gamma$, enables the system to
reach a stationary state. The characterization of the latter with varying of
$\lambda$ and $T$ is performed in the main text. However, we can already see
that the EoMs of the overlaps and the ones of the coherences are coupled via
terms proportional to the coherent control parameter $\lambda$. This suggests
that the addition of the quantum Hamiltonian (16) to the classical qPHN may
play a non-trivial role when analyzing the retrieval properties of the
network.

###  B.1 p=1 limit cycle phase

Figure 7: One memory LC phase. We set $\gamma=1$, $q=3$, $p=1$. Standard
deviation $\sigma_{m}$ of the overlap $m(t)$ at long times, with respect to
the fixed point solutions, as defined by Eq. (25). Initial conditions are
chosen (a) $m(0)<0$ and (b) $m(0)>0$. In the blue region $\sigma_{m}$ is
finite, and such a portion of the phase diagram corresponds to a limit cycle
phase.

At stationarity, Eqs. (24) display limit cycle solutions for a certain range
of temperatures $T$ and when the control parameter $\lambda$ is sufficiently
large. As shown in Fig. 4 in the main text, we analyze the standard deviation
of the overlap variables at long times, with respect to the fixed point
solutions of the EoMs. Additional results are reported in Fig. 7, where panels
(a) and (b) show the standard deviation defined by Eq. (25), with initial
conditions $m(0)<0$ and $m(0)>0$, respectively. The shaded region corresponds
to a finite value of the standard deviation $\sigma_{m}$.

Additionally, we perform a stability analysis of the fixed point paramagnetic
solution. Namely, defining $\vec{v}=(m,x,\bar{x},y,\bar{y})^{T}$, we linearize
the Eqs. (24) with respect to the solution
$\vec{v}=\vec{v}_{p}+\delta\vec{v}$, with $\vec{v}_{p}=0$. We obtain the
equations in the form $\delta\dot{\vec{v}}=J\delta\vec{v}$, where $J$ is the
Jacobian matrix. From the eigenvalues of the Jacobian matrix, say
$\vec{\zeta}$, we could identify a Hopf bifurcation [60], signaling a limit
cycle solution. An instance of this feature is shown in Fig. 8, which displays
two complex conjugate eigenvalues of the Jacobian matrix, $\zeta_{1}$ and
$\zeta_{2}$. The coherent control parameter is fixed to the value
$\lambda=2.5$, and the temperature is varied in the interval $T\in[0.5,1.5]$.
We can see that the real part of the eigenvalues changes sign when crossing
the imaginary axis, i.e. with finite imaginary part, highlighting the feature
of a Hopf bifurcation [60].

Figure 8:  Hopf bifurcation.  Complex plane representation of eigenvalues
$\zeta_{1,2}$ of the linearized Eqs. (24) with respect to the paramagnetic
solution. We fix $\lambda=2.5$ and the temperature is varied in the interval
$T\in[0.5,1.5]$. The change in sign of the eigenvalues while keeping a finite
imaginary part signals a Hopf bifurcation. The paramagnetic solution is
stable[unstable] for $\mathrm{Re}(\zeta_{1,2})<[>]0$. However, as the
imaginary axis is crossed by lowering the temperature, a limit cycle solution
appears, as discussed in the main text. Other parameters are $p=1$,
$\gamma=1$.

## References

  * Samuel [1959] A. L. Samuel, IBM Journal of Research and Development 3, 210 (1959).
  * Goodfellow  _et al._ [2016] I. Goodfellow, Y. Bengio, and A. Courville,  _Deep Learning_ (MIT Press, 2016) http://www.deeplearningbook.org.
  * Jordan and Mitchell [2015] M. I. Jordan and T. M. Mitchell, Science 349, 255 (2015).
  * LeCun  _et al._ [2015] Y. LeCun, Y. Bengio, and G. Hinton, Nature 521, 436 (2015).
  * Haykin [1998] S. Haykin,  _Neural Networks: A Comprehensive Foundation_ , 2nd ed. (Prentice Hall PTR, USA, 1998).
  * Crick [1989] F. Crick, Nature 337, 129 (1989).
  * Amit [1992] D. J. Amit,  _Modelling Brain Function: The World of Attractor Neural Networks_ , 1st ed. (Cambridge University Press, USA, 1992).
  * Montanaro [2016] A. Montanaro, npj Quantum Information 2 (2016), 10.1038/npjqi.2015.23.
  * Nielsen and Chuang [2011] M. A. Nielsen and I. L. Chuang,  _Quantum Computation and Quantum Information: 10th Anniversary Edition_ , 10th ed. (Cambridge University Press, USA, 2011).
  * Shor [1999] P. W. Shor, SIAM Review 41, 303 (1999).
  * Grover [1997] L. K. Grover, Phys. Rev. Lett. 79, 325 (1997).
  * Carleo  _et al._ [2019] G. Carleo, I. Cirac, K. Cranmer, L. Daudet, M. Schuld, N. Tishby, L. Vogt-Maranto, and L. Zdeborová, Rev. Mod. Phys. 91, 045002 (2019).
  * Carleo and Troyer [2017] G. Carleo and M. Troyer, Science 355, 602 (2017).
  * Gao and Duan [2017] X. Gao and L.-M. Duan, Nature communications 8, 1 (2017).
  * Huang and Moore [2017] Y. Huang and J. E. Moore, arXiv preprint arXiv:1701.06246 (2017).
  * Choo  _et al._ [2018] K. Choo, G. Carleo, N. Regnault, and T. Neupert, Phys. Rev. Lett. 121, 167204 (2018).
  * Saito [2018] H. Saito, Journal of the Physical Society of Japan 87, 074002 (2018), https://doi.org/10.7566/JPSJ.87.074002 .
  * Sharir  _et al._ [2020] O. Sharir, Y. Levine, N. Wies, G. Carleo, and A. Shashua, Phys. Rev. Lett. 124, 020503 (2020).
  * Torlai and Melko [2018] G. Torlai and R. G. Melko, Phys. Rev. Lett. 120, 240503 (2018).
  * Glasser  _et al._ [2018] I. Glasser, N. Pancotti, M. August, I. D. Rodriguez, and J. I. Cirac, Phys. Rev. X 8, 011006 (2018).
  * Liu  _et al._ [2018] Y. Liu, X. Zhang, M. Lewenstein, and S. J. Ran, arXiv preprint arXiv:1803.09111 (2018).
  * Pastori  _et al._ [2019] L. Pastori, R. Kaubruegger, and J. C. Budich, Phys. Rev. B 99, 165123 (2019).
  * Clark [2018] S. R. Clark, Journal of Physics A: Mathematical and Theoretical 51, 135301 (2018).
  * Schuld  _et al._ [2014] M. Schuld, I. Sinayskiy, and F. Petruccione, Quantum Information Processing 13, 2567 (2014).
  * Deng  _et al._ [2017] D. L. Deng, X. Li, and S. Das Sarma, Phys. Rev. X 7, 021021 (2017).
  * Biamonte  _et al._ [2017] J. Biamonte, P. Wittek, N. Pancotti, P. Rebentrost, N. Wiebe, and S. Lloyd, Nature 549, 195 (2017).
  * Rebentrost  _et al._ [2018] P. Rebentrost, T. R. Bromley, C. Weedbrook, and S. Lloyd, Physical Review A 98, 042308 (2018).
  * Aspuru-Guzik and Cao [2020] A. Aspuru-Guzik and Y. Cao, “Quantum artificial neural networks,” (2020), uS Patent App. 16/647,194.
  * Killoran  _et al._ [2019] N. Killoran, T. R. Bromley, J. M. Arrazola, M. Schuld, N. Quesada, and S. Lloyd, Phys. Rev. Research 1, 033063 (2019).
  * Mangini  _et al._ [2021] S. Mangini, F. Tacchino, D. Gerace, D. Bajoni, and C. Macchiavello, EPL (Europhysics Letters) 134, 10002 (2021).
  * Torrontegui and García-Ripoll [2019] E. Torrontegui and J. J. García-Ripoll, EPL (Europhysics Letters) 125, 30004 (2019).
  * Kristensen  _et al._ [2021] L. Kristensen, M. Degroote, P. Wittek, A. Aspuru-Guzik, and N. Zinner, npj Quantum Information 7, 59 (2021).
  * Cao  _et al._ [2017] Y. Cao, G. G. Guerreschi, and A. Aspuru-Guzik, “Quantum neuron: an elementary building block for machine learning on quantum computers,” (2017), arXiv:1711.11240 [quant-ph] .
  * Cong  _et al._ [2019] I. Cong, S. Choi, and M. D. Lukin, Nature Physics 15, 1273–1278 (2019).
  * Beer  _et al._ [2020] K. Beer, D. Bondarenko, T. Farrelly, T. J. Osborne, R. Salzmann, D. Scheiermann, and R. Wolf, Nature communications 11, 1 (2020).
  * Pons  _et al._ [2007] M. Pons, V. Ahufinger, C. Wunderlich, A. Sanpera, S. Braungardt, A. Sen(De), U. Sen, and M. Lewenstein, Phys. Rev. Lett. 98, 023003 (2007).
  * Gopalakrishnan  _et al._ [2012] S. Gopalakrishnan, B. L. Lev, and P. M. Goldbart, Philosophical Magazine 92, 353 (2012), https://doi.org/10.1080/14786435.2011.637980 .
  * Fiorelli  _et al._ [2020] E. Fiorelli, M. Marcuzzi, P. Rotondo, F. Carollo, and I. Lesanovsky, Phys. Rev. Lett. 125, 070604 (2020).
  * Behrman  _et al._ [2006] E. C. Behrman, K. Gaddam, J. E. Steck, and S. R. Skinner, “Microtubules as a quantum hopfield network,” in  _The Emerging Physics of Consciousness_, edited by J. A. Tuszynski (Springer Berlin Heidelberg, Berlin, Heidelberg, 2006) pp. 351–370.
  * Akazawa  _et al._ [2000] M. Akazawa, E. Tokuda, N. Asahi, and Y. Amemiya, Analog Integrated Circuits and Signal Processing 24, 51 (2000).
  * Rotondo  _et al._ [2018a] P. Rotondo, M. Marcuzzi, J. P. Garrahan, I. Lesanovsky, and M. Müller, Journal of Physics A: Mathematical and Theoretical 51, 115301 (2018a).
  * Fiorelli  _et al._ [2019] E. Fiorelli, P. Rotondo, M. Marcuzzi, J. P. Garrahan, and I. Lesanovsky, Phys. Rev. A 99, 032126 (2019).
  * Diamantini and Trugenberger [2006] M. C. Diamantini and C. A. Trugenberger, Phys. Rev. Lett. 97, 130503 (2006).
  * Hopfield [1982] J. J. Hopfield, PNAS 79, 2554 (1982).
  * Amit  _et al._ [1985] D. J. Amit, H. Gutfreund, and H. Sompolinsky, Phys. Rev. Lett. 55, 1530 (1985).
  * Amit  _et al._ [1987] D. J. Amit, H. Gutfreund, and H. Sompolinsky, Annals of Physics 173, 30 (1987).
  * Garrahan [2018] J. P. Garrahan, Physica A: Statistical Mechanics and its Applications 504, 130 (2018), lecture Notes of the 14th International Summer School on Fundamental Problems in Statistical Physics.
  * Lewenstein  _et al._ [2021] M. Lewenstein, A. Gratsea, A. Riera-Campeny, A. Aloy, V. Kasper, and A. Sanpera, Quantum Science and Technology (2021).
  * Breuer and Petruccione [2002] H. P. Breuer and F. Petruccione,  _The theory of open quantum systems_ (Oxford University Press, Great Clarendon Street, 2002).
  * Müller  _et al._ [2012] M. Müller, S. Diehl, G. Pupillo, and P. Zoller, Advances in Atomic, Molecular, and Optical Physics 61, 1 (2012).
  * Schindler  _et al._ [2013] P. Schindler, M. Müller, D. Nigg, J. T. Barreiro, E. A. Martinez, M. Hennrich, T. Monz, S. Diehl, P. Zoller, and R. Blatt, Nature Physics 9, 361 (2013).
  * Kanter [1988] I. Kanter, Phys. Rev. A 37, 2739 (1988).
  * Bollé  _et al._ [1992a] D. Bollé, P. Dupont, and J. Huyghebaert, Physica A: Statistical Mechanics and its Applications 185, 363 (1992a).
  * Bollé and Mallezie [1989] D. Bollé and F. Mallezie, Journal of Physics A: Mathematical and General 22, 4409 (1989).
  * Bollé  _et al._ [1991] D. Bollé, P. Dupont, and J. van Mourik, Journal of Physics A: Mathematical and General 24, 1065 (1991).
  * Bollé  _et al._ [1992b] D. Bollé, P. Dupont, and J. Huyghebaert, Phys. Rev. A 45, 4194 (1992b).
  * Wu [1982] F. Y. Wu, Rev. Mod. Phys. 54, 235 (1982).
  * Rotondo  _et al._ [2018b] P. Rotondo, M. Marcuzzi, J. P. Garrahan, I. Lesanovsky, and M. Müller, Journal of Physics A: Mathematical and Theoretical 51, 115301 (2018b).
  * Gayrard [1992] V. Gayrard, Journal of Statistical Physics 68, 977 (1992).
  * Strogatz [1994] S. H. Strogatz,  _Nonlinear dynamics and chaos: with applications to Physics, Biology_ (1994).
  * Takeuchi [1996] Y. Takeuchi,  _Global dynamical properties of Lotka-Volterra systems_ (World Scientific, 1996).
  * Goel  _et al._ [1971] N. S. Goel, S. C. Maitra, and E. W. Montroll, Rev. Mod. Phys. 43, 231 (1971).
  * Mézard  _et al._ [1990] M. Mézard, G. Parisi, and M.-A. Virasoro,  _Spin glass theory and beyond._ (World Scientific Publishing Co., Inc., Pergamon Press, 1990).