# Dynamical phase transitions in quantum reservoir computing

Rodrigo Martínez-Peña<EMAIL_ADDRESS>Instituto de Física
Interdisciplinar y Sistemas Complejos (IFISC, UIB-CSIC), Campus Universitat de
les Illes Balears E-07122, Palma de Mallorca, Spain  Gian Luca Giorgi
Instituto de Física Interdisciplinar y Sistemas Complejos (IFISC, UIB-CSIC),
Campus Universitat de les Illes Balears E-07122, Palma de Mallorca, Spain
Johannes Nokkala  Instituto de Física Interdisciplinar y Sistemas Complejos
(IFISC, UIB-CSIC), Campus Universitat de les Illes Balears E-07122, Palma de
Mallorca, Spain  Miguel C. Soriano  Instituto de Física Interdisciplinar y
Sistemas Complejos (IFISC, UIB-CSIC), Campus Universitat de les Illes Balears
E-07122, Palma de Mallorca, Spain  Roberta Zambrini<EMAIL_ADDRESS>Instituto de Física Interdisciplinar y Sistemas Complejos (IFISC, UIB-CSIC),
Campus Universitat de les Illes Balears E-07122, Palma de Mallorca, Spain

###### Abstract

Closed quantum systems exhibit different dynamical regimes, like Many-Body
Localization or thermalization, which determine the mechanisms of spread and
processing of information. Here we address the impact of these dynamical
phases in quantum reservoir computing, an unconventional computing paradigm
recently extended into the quantum regime that exploits dynamical systems to
solve nonlinear and temporal tasks. We establish that the thermal phase is
naturally adapted to the requirements of quantum reservoir computing and
report an increased performance at the thermalization transition for the
studied tasks. Uncovering the underlying physical mechanisms behind optimal
information processing capabilities of spin networks is essential for future
experimental implementations and provides a new perspective on dynamical
phases.

– Introduction. Unconventional computing is an interdisciplinary branch of
science that aims to uncover new computing and information processing
mechanisms in physical, chemical and biological systems [1]. In this field,
the challenge is to develop a device theory guaranteeing that a given system,
used as an analog computer, is able to accomplish a computational task. When
it comes to solving temporal tasks, a natural “computer” is represented by a
system exhibiting rich dynamical properties. An example of such approach can
be found in reservoir computing (RC), an unconventional framework belonging to
the broad family of machine learning, derived from recurrent neural networks
[2, 3, 4] but with the major advantage of low training cost and fast learning.
RC is also especially suited for hardware implementations [5, 6, 7, 8].

For big-data processing, an exceptional playground where a rich dynamics can
be exploited is certainly provided by quantum systems, whose exponentially
large number of degrees of freedom pushes them towards computational limits
that are not achievable by classical systems [9]. This is the potential
envisaged in quantum reservoir computing (QRC) [10, 11], as recently explored
in spin-based implementations [10, 12, 13, 14, 15], continuous-variable
bosonic systems [16, 17, 18, 19] and fermionic setups [20]. Efforts to
demonstrate proof-of-principle QRC physical experiments are ongoing [21, 22].
Although all the previous works provide examples of functioning quantum
reservoir computers, the fundamental issue raised at the beginning remains
open: what conditions must a physical system fulfill to be a good quantum
reservoir computer? The aim of this Letter is to establish the relation
between operation regime of complex computing systems and performance of QRC.

Networks of interacting spins enable complex dynamics providing a source of
memory needed for temporal tasks in QRC. The time evolution of these systems
and the conditions for thermalization have been recently debated in the
context of statistical physics. Indeed isolated quantum many-body systems can
display thermalization in local observables, as explained by the Eigenstate
Thermalization Hypothesis, which can be seen as the manifestation of
ergodicity in quantum mechanics [23, 24, 25]. A remarkable case of a dynamical
regime in which this hypothesis is violated is Many-Body Localization (MBL),
where strong disorder causes the emergence of an extensive number of local
integrals of motion that break down the thermalization hypothesis [26].
Indeed, such conserved quantities make local observables retain memory of
their initial states. Transitions between localization and thermalization
manifest in a critical change of a time-averaged order parameter and are
referred to as dynamical phase transitions [27, 26].

The different physical mechanisms underlying the presence or the absence of
thermalization deeply influence the computational capabilities of the
different dynamical phases. For instance, systems presenting MBL can provide
quantum memories at finite temperature [28] and avoid overheating in Floquet
systems [29, 30]. In quantum machine learning, MBL can improve the
trainability of parameterized quantum Ising chains [31]. Contrariwise,
localization can be computationally detrimental in quantum annealing [32, 33]
or quantum random walk algorithms [34, 35]. Our work establishes that optimal
information processing capabilities in QRC not only are favored in the ergodic
phase but also that the onset of this regime can be particularly advantageous.
We uncover the underlying physical mechanisms favoring machine learning and
also provide a new computing perspective on dynamical phases.

– Reservoir layer and dynamical phases. We choose as a reservoir a spin
network described by the transverse-field Ising Hamiltonian plus onsite
disorder:

$H=\sum^{N}_{i>j=1}J_{ij}\sigma_{i}^{x}\sigma_{j}^{x}+\frac{1}{2}\sum^{N}_{i=1}(h+D_{i})\sigma_{i}^{z},$
(1)

where $N$ is the spin number (from now on we will fix $N=10$), $h$ is the
magnetic field, $D_{i}$ is the onsite disorder, $\sigma^{a}_{i}$ ($a=x,y,z$)
are the Pauli matrices and $J_{ij}$ are the spin-spin couplings, randomly
selected from a uniform distribution in the interval $[-J_{s}/2,J_{s}/2]$ as
it is often done in QRC [10, 12, 13]. $D_{i}$ will be also randomly drawn from
the uniform distribution $[-W,W]$, where $W$ is the disorder strength. For
convenience, all the parameters will be expressed in units of $J_{s}$.

To characterize the dynamical phases of Eq. (1), we use as usual the ratio of
adjacent gaps
$r_{n}=\min[\delta_{n+1},\delta_{n}]/\max[\delta_{n+1},\delta_{n}]$, where the
gaps are $\delta_{n}\equiv E_{n}-E_{n-1}$, and $\\{E_{n}\\}$ is the sorted
list in ascending order of the Hamiltonian eigenvalues [36]. We compute the
eigenvalues via exact diagonalization limited to one symmetry sector, as the
model possesses a parity ($\mathcal{Z}_{2}$) symmetry [37]. In the localized
phase, the level spacing is expected to display a Poisson distribution with
$\braket{r}\simeq 0.386$, while in the ergodic phase, according to the random
matrix theory, $\braket{r}\simeq 0.535$ [38].

Figure 1:  Heat map of $\braket{r}$ for different values of the magnetic field
$h$ and the disorder strength $W$ in units of $J_{s}$. Results are averaged
over 1200 realizations.

The full dynamical phase diagram of Eq. (1) depending on magnetic field
strength and disorder is shown in Fig. 1, displaying four different regions:
two localization areas (regions I and III, black) and two ergodicity areas
(regions II and IV, bright yellow). The localized regime in I corresponds to a
MBL paramagnetic phase where eigenvalue statistics is Poissonian, while region
III corresponds to a spin-glass phase, where eigenvalue statistics is also
Poissonian, but it presents a mobility edge [39]. These regimes have been
described in different models as for the transitions III-IV [39, 40, 41], IV-I
[42, 43], or between different localization phases I-II-III [28, 44, 45].
While such transitions are strictly found in the thermodynamic limit [46, 47],
signatures are already evident for finite-size systems.

– Quantum reservoir dynamics. The QRC algorithm can be divided into three
steps associated to the relative system layers: (i) feed an input into the
dynamical system; (ii) let the reservoir, i.e. the spins network (1), evolve;
(iii) extract information from the reservoir, using all or some of its degrees
of freedom via an output layer [4]. Let us assume that our input is given by a
sequence $\\{s_{0},s_{1},\dots,s_{k},\dots\\}$ that is injected into the same
spin (named qubit $1$ for convenience) every time step $k$ [10]. This spin
state is updated every $\Delta t$ as follows:
$\rho_{1}^{(k)}=\ket{\psi_{s_{k}}}\bra{\psi_{s_{k}}}$, where
$\ket{\psi_{s_{k}}}=\sqrt{1-s_{k}}\ket{0}+\sqrt{s_{k}}\ket{1}$, with
$s_{k}\in[0,1]$. The completely positive, trace-preserving map summarizing the
input encoding and information processing is

$\rho(k\Delta t)=e^{-iH\Delta
t}\rho_{1}^{(k)}\otimes\text{Tr}_{1}\left\\{\rho[(k-1)\Delta
t]\right\\}e^{iH\Delta t},$ (2)

where $e^{-iH\Delta t}$ is the operator of the unitary dynamics and
$\text{Tr}_{1}\left\\{\cdot\right\\}$ denotes the partial trace performed over
the first qubit. The output layer will be built using some of the observables
of the system such as the projections $\braket{\sigma^{z}_{i}}$ of each spin
over the $z$-axis or the spin correlations
$\braket{\sigma^{z}_{i}\sigma^{z}_{j}}$ (Sect. V in [37]).

Figure 2:  Dynamics of observables $\braket{\sigma^{z}_{i}}$ with a binary
input ($s_{k}=\\{0,1\\}$). Parameters are (a) $W=0$, and $h/J_{s}=10$; (b)
$W/J_{s}=10$ and $h/J_{s}=1$. We use $J_{s}\Delta t=10$ here and in all the
next figures. Input is fed to the first spin (blue line) and the rest of lines
correspond to the other spins. The initial condition is a random density
matrix in both (a) and (b).

The spin network response to the input injection through the dynamics of the
observables $\braket{\sigma^{z}_{i}}$ provides insightful evidence as shown in
Fig. 2. The evolution of the observables in Fig. 2 (a) corresponds to the
ergodic region (IV in Fig. 1). Color lines represent the different spins,
being the input qubit the blue line. This plot displays that all spin
observables are driven to an input-dependent stationary state within each
$\Delta t$. The explanation behind this behaviour is based on the role of the
conserved quantities of the system. In the ergodic regions, total energy
$\braket{H}$ and parity $\braket{P}=\braket{\prod^{N}_{i}\sigma^{z}_{i}}$ are
the only conserved quantities and both of them are delocalized. As detailed in
Sect. II of [37], their values after the $k$th input injection
($\braket{H(k\Delta t)}$ and $\braket{P(k\Delta t)}$) only depend on the
initial condition of the system $\rho_{0}$ and the input history up to
$s_{k}$. Then, provided that the unitary dynamics is applied during a time
$\Delta t$ allowing for the Eigenstate Thermalization, all $local$
observables, like $\braket{\sigma^{z}_{i}}$ in Fig. 2 (a), only depend on
$\braket{H(k\Delta t)}$ and $\braket{P(k\Delta t)}$ up to finite-size
fluctuations [26]. Therefore, output observables become functions of the input
history through $\braket{H}$ and $\braket{P}$, as for the driven dynamics of
Fig. 2 (a). We anticipate that, resetting the first spin state with the input
protocol defined above implies a partial information erasure. Repeating this
operation several times amounts to losing all trace of the initial conditions
as for the convergence property addressed in the following.

Figure 2 (b) corresponds to the transition from ergodic to MBL (regions IV and
I respectively) and displays a significant change of the system response with
respect to the previous ergodic case. The observables show now little
correlation with respect to the dynamics of the first spin (blue line). This
behavior becomes more evident deep in the localized regimes, where none of the
$\braket{\sigma^{z}_{i}}$ are driven by the input, being instead determined by
the initial condition ([37], e.g. Fig. S2 (d)). The physical reason is that
the presence of an extensive number of local conserved quantities hinders the
information transport across the network. As a prominent effect of MBL, only
those conserved quantities involving the first spin are modified by the input,
while all the others keep memory of their initial conditions.

– Convergence. A fundamental property a system must exhibit to serve as RC is
the convergence or echo state property [2]. This means that, after repeated
input injection, the reservoir forgets its initial condition. This is closely
related to the so-called fading memory, that is the ability of the output
variables to only depend on the recent history of the input sequence [48]. The
convergence property is captured by the distance between two different
reservoir states after several (here $200$) input injections through the
protocol in Eq. (2). The initial conditions considered are two random density
matrices with a typical distance around $||\rho_{A}-\rho_{B}||\sim 0.044$
measured using the Frobenius norm, defined as
$||A||=\sqrt{\text{Tr}(A^{\dagger}A)}$ [37].

Figure 3: Convergence of the system for two different (random) initial
conditions after 200 inputs. Distance values are averaged over $600$
realizations, and those below the threshold of $10^{-8}$ are kept to this
minimum value for clarity. The white-dashed line corresponds to an
intermediate value of $\braket{r}=0.46$ in Fig. 1, used to guide the eye.

Convergence as a function of $h$ and $W$ is shown in Fig. 3 and exhibits an
insightful correspondence with the phase diagram in Fig. 1. Indeed, the
convergence property is enhanced in the ergodic phase enabling then the
realization of QRC with the spins reservoir in this regime. The influence of
different initial conditions persists in the localized phases, which hinders
the QRC performance as it will be shown later. Even if the map (2) is
contractive due to the partial trace, the local conserved quantities make the
contractiveness of this dynamical map much weaker than in the ergodic phase.
Indeed, in the MBL phase the creation of entanglement between the first qubit
and the rest of the network is very weak, as it only grows logarithmically in
time [26]. This provides a different focus of dynamical phases, in terms of
their response to external perturbations.

– RC performance. While convergence can be seen as a necessary condition to
identify suitable RC systems, a characterization of the information processing
capabilities of the spin network is needed in order to determine if the
reservoir computer can accomplish a given task. A convenient known advantage
of RC is that the output layer is the only one that needs to be trained, by
optimizing a linear combination of the responses of the reservoir to the task
at hand [4]. Given the set of observables $\textbf{x}_{k}$ chosen as output,
we write the output layer as ${y_{k}=\textbf{w}^{\top}\textbf{x}_{k}}$, where
w are the weights that are adapted by minimizing the error with respect to a
target function $\bar{y}_{k}$. Training of the output weights w is usually
done using a linear regression [4, 37].

To evaluate the reservoir performance, we will consider two different specific
tasks as well as a more general indicator for the processing capacity [49].
Let us start with the nonlinear auto-regressive moving average (NARMA) model,
which is widely used to characterize recurrent neural networks [50]. The
general NARMA$n$ task is defined as:

$\bar{y}_{k}=0.3\bar{y}_{k-1}+0.05\bar{y}_{k-1}\left(\sum^{n}_{j=1}\bar{y}_{k-j}\right)+1.5s_{k-n}s_{k-1}+0.1,$
(3)

where $n$ is the maximum delay, $\bar{y}_{k}$ is the target and $s_{k}$ is a
random input uniformly distributed, as detailed in Sect. VI of [37]. For this
task, the quantum reservoir computer needs to learn to emulate Eq. (3) from
the random input $s_{k}$, i.e. to reproduce a quadratic nonlinear function of
the input sequence up to a maximum fixed delay. We fix the maximum delay of
Eq. (3) to $n=10$. In order to address performance both in nonlinear and
linear temporal tasks, we also evaluate the linear temporal task
$\bar{y}_{k}=s_{k-\tau}$ fixing $\tau=10$. The value $10$ in both tasks sets
the target memory: indeed the reservoir needs a memory of at least 10 past
inputs at each time step $k$ to solve them. The performance of the spin
network for each task can be measured as
$C={\text{cov}^{2}(\bar{\textbf{y}},\textbf{y})}/[{\sigma^{2}(\textbf{y})\sigma^{2}(\bar{\textbf{y}})]}$
where cov is the covariance, $\sigma^{2}$ is the variance and y and
$\bar{\textbf{y}}$ are the time series of prediction and target respectively.
The performance $C$ is bounded between $0$ (target not approximated, the
system is useless for this task) and $1$ (perfect match between prediction and
target).

Figure 4: Performance covariance $C$ for the NARMA10 and linear memory
$\tau=10$ tasks versus $h$ (a) and $W$ (b). We took $W=0$ for (a) and
$h/J_{s}=10$ for (b). We represent the average value of $C$ over 100
realizations of the random network with a solid line, while the shadows
represent the standard deviation. The grey-dashed line corresponds to the
intermediate value $\braket{r}=0.46$ in Fig. 1, used to guide the eye. The
initial condition for all the realizations is the maximal coherent state
$\rho_{0}=\frac{1}{2^{N}}\sum\ket{i}\bra{j}$.

In this work, we benefit from the high-dimensionality of the Hilbert phase
space of the quantum reservoir by considering an output layer of $O=75$
observables (going beyond classical RC with one-body observables): all the
local spin projections $\braket{\sigma^{x}_{j}}$, $\braket{\sigma^{y}_{j}}$
and $\braket{\sigma^{z}_{j}}$ $(1\leq j\leq N)$, plus the two-spin
correlations along the $z$-axis $\braket{\sigma^{z}_{i}\sigma^{z}_{j}}$
$(1\leq i,j\leq N$, $i\neq j$). The measurement of these observables is
experimentally feasible (Sect. V in [37]). Figure 4 displays the figure of
merit $C$ for the performance of the spin-based reservoir in the NARMA10 and
linear memory tasks. In Fig. 4 (a), we plot $C$ versus $h$ keeping $W=0$,
while in Fig. 4 (b) we fix $h/J_{s}=10$ and vary $W$, in order to detect the
change in performance at the transition between different dynamical phases.
Localized regimes (regions I and III) show the smallest values of $C$: the
ability to reproduce the target is poor due to their slow convergence
properties and actually can be influenced both by the initial elapsed time and
the specific choice of initial conditions for the reservoir. In contrast, the
ergodic regime (region IV) shows a higher performance achieved rapidly and
independently of the reservoir initial state.

Interestingly, in both plots $C$ shows a peak at the transitions between
ergodic and localized phases. To show the generality of our results beyond
specific tasks and to shed light on the possible performance enhancement at
the phase transition, we evaluate [37] the information processing capacity
(IPC) [49]. We find that the linear memory builds up first as we move away
from the localized phase into the ergodic one. Deep into the ergodic region,
nonlinear memory dominates. This trade-off between linear and nonlinear memory
brings the performance enhancement at the transition for the NARMA$10$ task.
Although the IPC results indicate that one could find nonlinear tasks where
the optimal working point is found in the ergodic region and not at the
transition, it is often the case that RC tasks precisely require a combination
of linear memory and nonlinearity, as often reported at the edge of stability
between different dynamical regimes [51].

– Conclusions. High performance in QRC can be achieved thanks to the large
dimensionality of the Hilbert space. Still, the performance of a system to be
used as a quantum reservoir computer crucially depends on its operation
regime. We showed in this work that localization, because of the presence of
local conserved quantities, is detrimental for an optimal information
processing performance due to a slow convergence [52]. We demonstrate this in
specific tasks as well as in the quantification of the reservoir memory
through the IPC. In contrast, the ergodic phase offers a suitable scenario for
the convergence property and facilitates efficient information extraction.
Different tasks can be solved by exploiting the trade-off between linear and
nonlinear memory at the phase transition, and actually the onset of
thermalization can be particularly advantageous for QRC, a feature reminiscent
of the performance enhancement found in classical RC at the edge of stability
[51]. Our QRC study offers an original perspective on thermal and localized
phases in terms of their ability to process information and can be further
explored in the context of quantum correlations, information scrambling, OTOC
(out-of-time-order correlators) [53] and transient real-time evolution of
Loschmidt echoes [54, 27].

Furthermore, our results define the proper conditions for experimental
implementations of QRC. We showed the importance of tuning the reservoir at
the onset of thermalization, which can be easily achieved by controlling the
(average) strength of the magnetic field. Another relevant issue concerns the
network topology. It is already accepted that random connections are necessary
for an optimal performance, avoiding redundancies between different degrees of
freedom. But it is not enough. We show that even those topologies with
disorder leading to an extensive number of conserved quantities are not suited
for RC. Strategies for on-line data processing addressing quantum measurement
[22, 55], need to be further explored [37]. First experiments involve ensemble
computing, obtained by taking many copies of the reservoir [21], or rely on
the use of non-demolition measurements [22]. Several platforms, ranging from
trapped-ion quantum simulators [43, 56, 57], to optical lattices [58, 59] to
superconducting circuits [60, 61], or photonic simulators [62, 63] are mature
to establish the potential of QRC towards applications, both for classical and
quantum time series processing [11].

###### Acknowledgements.

We acknowledge the Spanish State Research Agency, through the Severo Ochoa and
María de Maeztu Program for Centers and Units of Excellence in R&D
(MDM-2017-0711) and through the QUARESC project (PID2019-109094GB-C21 and
-C22/ AEI / 10.13039/501100011033); We also acknowledge funding by CAIB
through the QUAREC project (PRD2018/47). Part of this work has been funded by
MICINN/AEI/FEDER and the University of the Balearic Islands through a “Ramon y
Cajal” Fellowship (RYC-2015-18140) for MCS and a predoctoral fellowship
(MDM-2017-0711-18-1) for RMP. GLG is funded by the Spanish Ministerio de
Educación y Formación Profesional/Ministerio de Universidades and co-funded by
the University of the Balearic Islands through the Beatriz Galindo program
(BG20/00085).

## Supplemental Material

##  I Numerical methods

To characterize the dynamical phases of the system’s Hamiltonian
$H=\sum^{N}_{i>j=1}J_{ij}\sigma_{i}^{x}\sigma_{j}^{x}+\frac{1}{2}\sum^{N}_{i=1}(h+D_{i})\sigma_{i}^{z}$
(Fig. 1 of the main text), we employ the ratio of adjacent gaps
$r_{n}=\min[\delta_{n+1},\delta_{n}]/\max[\delta_{n+1},\delta_{n}]$, where the
gaps are $\delta_{n}\equiv E_{n}-E_{n-1}$, and $\\{E_{n}\\}$ is the sorted
list in ascending order of the Hamiltonian eigenvalues [36]. The transition
from ergodic (thermalizing) to non-ergodic regions in the parameter space can
be tracked down looking at the eigenvalue statistics: in regions where the
system has an extensive number of integrals of motion, the eigenvalues
belonging to different sectors can be seen as independent random variables, so
that they are spaced according to a Poisson distribution ($\braket{r}\sim
0.386$, where $\langle\rangle$ is the average over all the gaps); on the other
hand, in the presence of thermalization, there is level repulsion and the
statistics follows the Wigner-Dyson distribution ($\braket{r}\sim 0.53$). The
eigenvalues of the Hamiltonian are computed via exact diagonalization. For a
system of 10 spins this implies to diagonalize a matrix of $\sim 10^{6}$
elements obtaining a sequence of about 1000 eigenvalues. It is convenient that
as the model possesses a parity ($\mathcal{Z}_{2}$) symmetry, the eigenvalue
calculation is limited to one symmetry sector. We have employed the QuSpin
library for this purpose [64].

When considering this system as the reservoir for learning purposes, the
dynamics is instead governed by the dynamical map we repeat here

$\rho(k\Delta t)=e^{-iH\Delta
t}\rho_{1}^{(k)}\otimes\text{Tr}_{1}\left\\{\rho[(k-1)\Delta
t]\right\\}e^{iH\Delta t},$ (4)

being $\rho_{1}$ the density matrix of the first spin (where we introduce the
input, see main text), $e^{-iH\Delta t}$ the operator of the unitary dynamics
and $\text{Tr}_{1}\left\\{\cdot\right\\}$ denotes the partial trace performed
over the first qubit. This completely positive and trace preserving (CPTP) map
of the QRC algorithm is applied computing again the Hamiltonian eigenvalues
via exact diagonalization, taking the whole spectrum this time. In this way,
we obtained the results of Figs. 2, 3 and 4 of the main text. Some elements of
the code were generated with the help of the QuTiP library [65], such as the
random initial density matrices of Fig. 3 of the main text.

##  II Role of ergodicity

In order to have a theoretical insight about the results presented in the main
text, let us analyze how information enters and leaves the system. First, let
us explain how the input protocol affects the overall picture. After each
input injection, classical information about the input is encoded in the first
spin, which instantaneously modifies the expected value of operators with
support in this spin. For example, observables like $\braket{\sigma^{z}_{1}}$,
$\braket{\sigma^{x}_{1}\sigma^{x}_{2}}$ or the total energy $\braket{H}$
(which has support over all qubits) are instantaneously changed by the input
protocol. In particular, those observables with support on the first spin that
are conserved quantities become functions of both the initial state of the
system and the input sequence, since the unitary dynamics do not affect them.
In the case of the total energy $\braket{H}$, this dependence can be expressed
as:

$\braket{H(k\Delta t)}=\text{Tr}[H\rho(k\Delta
t)]=E(\rho_{0},s_{0}\dots,s_{k-1},s_{k}),$ (5)

where $\rho_{0}$ is the initial state of the system and $s_{k}$ is the input
at time step $k$. In the ergodic regions, total energy $\braket{H}$ and parity
$\braket{P}=\braket{\prod^{N}_{i}\sigma^{z}_{i}}$ are the only conserved
quantities, which are global operators with support over all the qubits.
Provided that the unitary dynamics is applied during a time $\Delta t$ long
enough, the Eigenstate Thermalization Hypothesis guarantees that all local
observables $\mathcal{O}$ only depend on $\braket{H(k\Delta t)}$ and
$\braket{P(k\Delta t)}$ up to finite-size fluctuations [26]. Indeed, after $k$
input injections, the expected value of any of such local observables can be
well approximated as

$\langle\mathcal{O}(k\Delta t)\rangle\simeq\frac{{\rm
Tr}[\mathcal{O}e^{-\beta_{k}H-\gamma_{k}P}]}{\mathcal{Z}_{k}},$ (6)

with the partition function $\mathcal{Z}_{k}={\rm
Tr}[e^{-\beta_{k}H-\gamma_{k}P}]$ and where the Lagrange multipliers
$\beta_{k}$ and $\gamma_{k}$ can in principle be calculated respectively as

$\displaystyle\langle H(k\Delta
t)\rangle=-\frac{\partial\log\mathcal{Z}_{k}}{\partial\beta_{k}},$ (7)
$\displaystyle\langle P(k\Delta
t)\rangle=-\frac{\partial\log\mathcal{Z}_{k}}{\partial\gamma_{k}}.$ (8)

Then, local observables become functions (each of them with its own functional
dependence) of the input history through $\braket{H}$ and $\braket{P}$.

Figure 5: Single realizations of the energy $\braket{H}$ and parity
$\braket{P}$ dynamics under the same input sequence with different initial
conditions. Parameters are $W=0$ and $h/J_{s}=10$ (ergodic region) for (a) and
(c), and $W/J_{s}=100$ and $h/J_{s}=1$ (MBL region) for (b) and (d). The
number of spins is $N=10$ and the time between inputs is $J_{s}\Delta t=10$.
The color lines correspond to different initial conditions: all spins up in
the $z$ basis (blue), all spins down in the $z$ basis (red), half spins up and
half spins down in the $z$ basis (green), half spins up and half spins down in
the $x$ basis (yellow). These particular initial conditions have been chosen
for the sake of visualization.

In Fig. 5, we show an example of the dynamics of $\braket{H(k\Delta t)}$ and
$\braket{P(k\Delta t)}$ in the ergodic phase for different initial conditions
but the same input sequence. As it can be seen in Figs. 5 (a) and (c), all the
trajectories converge to one single trajectory, which is just a consequence of
the convergence property studied in the main text. What we can learn here is
that the partial trace operation (see Eq. 4) that accompanies each input
injection implies a partial information erasure, as the presence of global
conserved quantities implies a significant amount of entanglement between the
first spin and the rest of the network. Repeating this operation enough times
amounts to losing all trace of the initial conditions and inputs in the
distant past, which ensures the convergence observed in Fig. 3 of the main
text. Then, after several input injections, energy and parity can be expressed
as

$\begin{split}&\braket{H(k\Delta t)}\simeq E(s_{k-\alpha},\dots,s_{k}),\\\
&\braket{P(k\Delta t)}\simeq p(s_{k-\alpha},\dots,s_{k}),\end{split}$ (9)

where both $E$ and $p$ are only functions of the input history up to a time
step $\alpha$ in the past. As a consequence of Eq. (9), our reservoir is
equipped with the so-called fading memory property, which means that two input
sequences that are close in the recent past produce outputs that are close in
the present. This fading memory observed in the collective conserved
quantities $\braket{H}$ and $\braket{P}$ is responsible for the good
functioning of the reservoir observed in Fig. 4 of the main text. Indeed, the
collectiveness of these degrees of freedom propagates such fading memory to
all degrees of freedom of the reservoir through Eq. (6).

The underlying physics is completely different for the localized regimes,
where the presence of local integrals of motion, that is, of conserved
quantities that have support only on parts of the reservoir, hinders the
information transport across the network (Fig. 5, right panels (b) and (d)).
Let us consider the Many-Body Localization (MBL) paramagnetic phase (region I
in Fig. 1 of the main text). As in the ergodic phases, local observables in a
MBL phase arrive at a stationary state after a transient time of the unitary
dynamics [26]. The main difference with respect to the ergodic phase is that
now the stationary value of local observables is deeply influenced by the
conserved quantities with support on them, and will be weakly affected by the
partial trace. Then, information erasure will be much less effective, which
makes the reservoir unable to efficiently forget its initial state. Once we
introduce the input protocol, only the observables with support in the first
spin will be significantly affected after each input injection and after each
partial trace application. This slows down the flow of information in the
system: the input enters but hardly arrives at the other spins, and
information contained in local conserved quantities that are not connected to
the first spin does not easily flow back to it and persists in the network.

As a consequence, from the point of view of the learning algorithm, not all
local observable allow for an efficient extraction of the input information
and cannot be exploited, so that the performance is deteriorated: deep in the
MBL phase, input information is almost exclusively contained in the first
spin, and is erased after each partial trace. This effect can be seen in Fig.
6 (d).

In the spin-glass phase (region III) we find a mobility edge, but there are
arguments in favor of the presence of local conserved quantities: near-
integrable systems ($h\ll 1$) can present them [66, 67]; furthermore, as
conjectured in [68], emergent integrability can be extended to systems with
mobility edges, obtaining the conserved quantities from the localized
eigenstates.

##  III Dynamical response of the spin network reservoir

In the main text, we characterized the response of our system to the input by
showing the dynamics of the observables $\braket{\sigma^{z}_{i}}$. For the
sake of succinctness, we showed in Fig. 2 the dynamics of the observables that
correspond to the ergodic region (Fig. 2 (a)) and the transition from the
ergodic phase to MBL paramagnetic (Fig. 2 (b)). Here we provide an additional
characterization with examples in the remaining dynamical phases.

Fig. 6 presents the dynamics of the observables $\braket{\sigma^{z}_{i}}$ as
in Fig. 2 of the main text and over a longer time interval. The different
lines in Fig. 6 correspond to the value of the observable for each spin as a
function of time. The blue line represents the dynamics of the first spin,
where we inject the input every $\Delta t$. Figs. 6 (a) and (b), correspond to
the plots as in Fig. 2 of the main text, i.e. we have the parameters $W=0$ and
$h/J_{s}=10$ for (a) and $W/J_{s}=10$ and $h/J_{s}=1$ for (b), but with twice
the time window in the $x$-axis. Fig. 6 (c) represents a point in the spin-
glass (region III in Fig. 1 of the main text), with parameters $W=0$ and
$h/J_{s}=0.01$. Fig. 6 (d) represents a point in the MBL paramagnetic (region
I), with parameters $W/J_{s}=100$ and $h/J_{s}=1$.

Figure 6: Dynamics of observables $\braket{\sigma^{z}_{i}}$ with a binary
input ($s_{k}=\\{0,1\\}$) for $N=10$ spins. Parameters are $J_{s}\Delta t=10$,
(a) $W=0$, and $h/J_{s}=10$ (region IV); (b) $W/J_{s}=10$ and $h/J_{s}=1$
(transition between regions I and II); (c) $W=0$ and $h/J_{s}=0.01$ (region
III); (d) $W/J_{s}=100$ and $h/J_{s}=1$ (region I). Input is fed to the first
spin (blue line). The initial condition for all the realizations is a random
density matrix.

We note that Fig. 6 (a) is a paradigmatic example of a driven system, showing
a clear correlation between the dynamics of the observables and the input.
However, Figs. 6 (b), (c) and (d) present to us cases where the driving of the
input has a lesser influence than in Fig. 6 (a). The fact that some
observables in Fig. 2 (b) seem to follow better the trend of the first qubit
is a nontrivial point. After analyzing several examples of different
reservoirs, with couplings between spins both with random strengths and
any/only positive/only negative signs (not shown), our conclusion is that it
does not appear to be a local effect but rather a global one, determined by
the whole network structure. On the other hand, Figs. (c) and (d) correspond
to the extreme cases where there is no apparent correlation between
observables and input. In both cases, we argue that the presence of local
integrals of motion is the reason for this loss of correlation between the
dynamics of the observables and the input. In (d), we can identify that
localization freezes even the dynamics of the first spin after input
injections. On the other hand, Fig. 6 (c) also shows a case where the rest of
the spins do not respond to the input injection. In contrast, the dynamics of
the whole network depends on the initial condition (choosing different initial
conditions, not shown, leads to completely different trajectories of the
observables). Indeed, local conserved quantities involving the first spin
hinder the spreading of information, which explains the lack of response from
the other spins.

##  IV On the convergence properties of the different dynamical phases

In this Section, we tackle the dependence on $\Delta t$ and the number of
input injections for the convergence property. The current discussion relates
to the results presented in Fig. 3 of the main text. The evolution of Eq. (4)
in Sec. I is evaluated exactly for different initial conditions. First, our
numerical results indicate that this completely positive and trace-preserving
map is strictly contractive (for the considered random initial conditions),
i.e. $||\rho_{A}(k\Delta t)-\rho_{B}(k\Delta t)||<||\rho_{A}[(k-1)\Delta
t]-\rho_{B}[(k-1)\Delta t]||$. The only source of contractiveness leading to a
convergence of the dynamics after consecutive input injections is the partial
trace over the first spin, since the unitary dynamics leaves invariant the
distance between density matrices. Then, the more input injections we make
(with the associated partial trace), the smaller is the distance between the
two initial conditions. Interestingly, as discussed in Sec. I and illustrated
in Fig. 7 below, we find that the dynamical phase plays a determinant role in
the convergence.

Figure 7: Convergence of the system for two different (random) initial
conditions over repeated input injection for various periods of injection
$\Delta t$. Parameters are (a) $W=0$ and $h/J_{s}=0.01$ (region III); (b)
$W=0$ and $h/J_{s}=10$ (region IV); (c) $W/J_{s}=100$ and $h/J_{s}=1$ (region
I). Lines are averaged over 100 realizations. The number of spins is $N=10$.

Figs. 7 (a), (b) and (c) represent the distance between two random initial
conditions over repeated input injection. Figures 7 (a) and (c) correspond to
the spin-glass and MBL paramagnetic phases respectively, while Fig. 7 (b)
belongs to the ergodic phase. In the three cases, we present the distance for
different periods of injection $\Delta t$ such that we are able to show the
influence of this parameter in each phase. While for short $\Delta t$
($J_{s}\Delta t\leq 0.1$) the results are very similar for all dynamical
regimes (the dynamical map is basically the identity), we find a large
difference between the behavior of the convergence for the localized regimes
and the ergodic phase for large $\Delta t$. Figures 7 (a) and (c) present
curves with a small slope with respect to the number of consecutive injected
inputs. In addition, these two cases present a minor dependence on $\Delta t$.
In contrast, the value of $\Delta t$ has a strong influence over the
convergence property in the ergodic phase, being 200 consecutive input
injections more than enough to find a clear convergence when $J_{s}\Delta
t\geq 10$. This is consistent with the physical behavior characterizing the
ergodic phase: after the relaxation time all the spins of the network achieve
coherent coupling. If $\Delta t$ is smaller than this relaxation time, the
convergence is hindered because entanglement does not significantly build-up
and information about the initial conditions does not have enough time to flow
back to the first spin –and to be erased after each input injection. On the
other hand, looking at the localized phases, an extremely slow reduction of
the distance between states is reported, making the network in this phase
unsuitable for QRC. We remark that reducing at the minimum the injection time
interval is a main requirement in order to process data in RC, being
processing speed actually a major performance aspect. We have therefore chosen
$J_{s}\Delta t=10$ as the time interval for our simulations of the main text
because it represents a short time scale that still allows relaxation to the
thermal state in the ergodic phase, ensuring the validity of Eqs. (9) above.

##  V Notes on physical implementation

Several comments about the physical implementation of spin-based QRC are in
order. First, we describe the construction of the output layer. The $4^{N}$
observables of a system of $N$ qubits can be defined as $x_{j}(k\Delta
t)=\text{Tr}[B_{j}\rho(k\Delta t)]=\braket{B_{j}(k\Delta t)}$, where $B_{j}$
is the tensor product of $N$ Pauli operators taken from
$\\{I,\sigma^{x},\sigma^{y},\sigma^{z}\\}$ (for each qubit). The output layer
will be made of some of the observables of the system, such as the projections
$\braket{\sigma^{z}_{i}}$ of each spin over the $z$-axis or the spins
correlations $\braket{\sigma^{z}_{i}\sigma^{z}_{j}}$. Obtaining these
observables is a fundamental question in QRC: we want to harness the large
Hilbert space that a quantum system offers, but the processing of information
in RC is thought for an online mode. The present approach is based on an
ensemble of quantum systems as in [10], where one can disregard the
measurement back-action, allowing the online mode computation. Nuclear
magnetic resonance experiments, where localized phases can be observed, [69,
70] are suitable for this approach [21]. On the other hand, several platforms,
ranging from trapped-ion quantum simulators [43, 56, 57], to optical lattices
[58, 59] to superconducting circuits [60, 61], or photonic simulators [62, 63]
are already mature to establish the potential of QRC towards applications. In
this case we require the repetition of the experiment to obtain the expected
value of observables. But this can be seen as an opportunity to exploit a
large number of degrees of freedom. As it is done in ion-trap experimental
papers [71, 56], expected values of correlations of Pauli matrices can be
computed from projective measurements of the individual spins. For example,
consider that we want to compute the expected values
$\braket{\sigma^{z}_{1}}$, $\braket{\sigma^{z}_{3}}$ and
$\braket{\sigma^{z}_{1}\sigma^{z}_{3}}$. The expected value of the spin
projections are computed with the average of the projective measurements over
each spin. For the first spin, one can write this like
$\braket{\sigma^{z}_{1}}=\frac{1}{M}\sum^{M}_{i}z_{1i}$, where $M$ is the
number of experimental realizations and $z_{1i}$ accounts for the projection
obtained over the first spin for measurement $i$ (a projective measurement
produces either $z_{1i}=1$ or $z_{1i}=-1$). One can do the same for the third
spin: $\braket{\sigma^{z}_{3}}=\frac{1}{M}\sum^{M}_{i}z_{3i}$. Finally, to
compute the correlation $\braket{\sigma^{z}_{1}\sigma^{z}_{3}}$ one just needs
to use the product of the projections and make the average:
$\braket{\sigma^{z}_{1}\sigma^{z}_{3}}=\frac{1}{M}\sum^{M}_{i}z_{1i}z_{3i}$.
Thus, it is possible to obtain the correlations in the $z$-axis with the
measurements of the spin projections in the $z$-axis, being able to access
nontrivial degrees of freedom of our system. We note that one has to repeat
more experiments to obtain the projections and correlations in other
directions, but there are several papers that propose strategies to surpass
this overhead. In particular, there is a recent paper by Huang et al. [72]
where they theoretically and numerically demonstrate that one can efficiently
estimate a large number of observables (correlations with 2,3,… spins) from a
reasonable number of projective measurements, based on the formalism of shadow
tomography.

Another fundamental point in QRC is the characteristics of the quantum model
simulated by the physical platform. Different experimental platforms will
provide different network structures, and then, different couplings in the
Hamiltonian. However, the phenomenon of dynamical phase transition in the
transverse-field Ising model with long-range interactions of the couplings
seems rather universal, either from an order or glassy ground state towards
thermalization [56, 39, 40], or from the ergodic phase to the MBL phase [43,
42]. We have checked that our results are consistent with the convergence
property of other long-range configurations. However, we have not explored the
performance of systems without any kind of disorder. We foresee that reducing
the number of degrees of freedom with extra symmetries would decrease the
performance of the system. On the other hand, topologies with short-range
interactions have not been explored either.

##  VI Reservoir Computing Training Methodology

In RC, training of the output weights w is usually done using a linear
regression for vector $\textbf{y}=X\textbf{w}$ with respect to
$\bar{\textbf{y}}$, where $X$ is a $L\times(O+1)$ matrix that collects the
reservoir states at different times, $L$ is the length of the input sequence
and $O$ is the number of output variables, that we remind correspond to
expectation values of observables. A constant bias term $x_{k,O+1}=1$ is added
for an optimal training [4]. Sequences of $5000$ inputs are used for Fig. 4 of
the main text: $1000$ to wash out the initial conditions, $2000$ for training,
and $2000$ for prediction, that is to assess the performance in the task. The
linear regression is chosen as the training method to offload and simplify the
work of the output layer, harnessing the nonlinear processing of information
from the reservoir.

##  VII NARMA task

The nonlinear auto-regressive moving average (NARMA) task belongs to a large
family of problems that researchers in the Machine Learning community
(especially in the Recurrent Neural Networks community) try to solve. These
kinds of tasks consist in approximating a nonlinear and temporal function of
the input sequence. The NARMA task in particular consists in approximating a
quadratic nonlinear function of the input sequence up to a maximum fixed
delay. The general NARMA$n$ model is defined as:

$\bar{y}_{k}=0.3\bar{y}_{k-1}+0.05\bar{y}_{k-1}\left(\sum^{n}_{j=1}\bar{y}_{k-j}\right)+1.5s_{k-n}s_{k-1}+0.1,$
(10)

where $n$ is the maximum delay, $\bar{y}_{k}$ is the target and $s_{k}$ is a
random input uniformly distributed. In the main text, we fix the maximum delay
of Eq. (10) to $n=10$, which requires a significant memory capacity, and
restrict the input to the interval $[0,0.2]$ to prevent divergences [73].

To give some intuition about what it means to realize the NARMA task, we show
in Fig. 8 some of the results that we obtain for a single realization of the
learning algorithm. Figure 8 (a) contains a temporal window of the target and
output sequences. The target sequence (black line) represents the trajectory
of Eq. 10. The color lines represent the trajectories that the quantum
reservoir computer produces as an output after the training process. The
selection of parameters corresponds to the three different regions that Fig. 4
(a) of the main text contain: ergodic phase ($h=10$), transition ($h=0.1$) and
spin-glass phase ($h=0.01$). As we observe in this figure, the realization of
the transition region is the best one approximating the target function, while
the realization of the spin-glass phase does not follow the target. Trying to
get a better insight, we represent in Fig. 8 (b) the error between the target
trajectory and the outputs of the reservoir in the same temporal window. For
this time interval, it is clear that $h=0.01$ is the worst case scenario, but
it is less straightforward to assess if $h=0.1$ is better than $h=10$. Then,
we need to account for the differences between target and reservoir output
trajectories during a larger time interval. In Fig. 8 (c) we show a histogram
of all the differences obtained between target and outputs during the test
sequence (2000 time steps). We observed that the three examples have the
maximum of the distribution close to the origin. However, the transition
example ($h=0.1$) has the highest maximum with the shortest tail, representing
the case with the smallest errors.

Figure 8: Single realizations of the NARMA10 task for different values of $h$.
Figure (a) represents a temporal window of the target and output trajectories.
Figure (b) corresponds to the error between target and output reservoirs for
the same temporal window, and Fig. (c) is a histogram of all the errors during
the test. The total number of time steps was 5000: 1000 for washing out the
initial conditions, 2000 for the training and 2000 for the test. Parameters
are $N=10$, $W=0$ and $J_{s}\Delta t=10$. The initial condition for all the
realizations is the maximal coherent state
$\rho_{0}=\frac{1}{2^{N}}\sum\ket{i}\bra{j}$.

##  VIII On the information processing capacity

It is impractical to cover the performance of the quantum spin network for all
conceivable benchmark tasks. In this spirit, it was recently shown that one
can provide a task-independent characterization of the computational
capabilities of a dynamical system acting as a reservoir computer by measuring
its Information Processing Capacity (IPC) [49]. The authors demonstrated that
the total computational capacity of a dynamical system is bounded by the
number of linearly independent variables that we use for the output. In
addition, the computed IPC can only saturate the bound when the dynamical
system possesses fading memory, which is strongly connected to the convergence
property (or echo state property) and under some conditions (compactness of
the input and reservoir state spaces, contractivity and continuity of the
reservoir dynamics [48]) both of them are fulfilled simultaneously.

In this context, the computational capabilities of a reservoir are identified
as the ability of performing both temporal and nonlinear tasks. Then,
computing all of these capabilities would imply evaluating all the possible
linear and nonlinear functions of the input sequence that our system can
approximate. The original theory developed in [49] simplifies this problem by
choosing a set of orthogonal functions. The target functions are defined as
the product of Legendre polynomials of a given degree $d_{i}$, such that the
sum of the degrees of all the multiplied polynomials add up to a given degree
$d\geq 1$ of nonlinearity. Then, the total capacity can be divided into linear
($d=1$) and nonlinear contributions ($d>1$). Besides, the polynomials will be
functions of the input at different times in the past. In this way, our
targets not only account for nonlinear functions but also temporal maps.

The capacity coefficient to quantify how well a system reproduces a target
function can be defined as in [49]:

$C_{L}(X,\textbf{y})=1-\frac{\text{min}_{\textbf{w}}\text{MSE}_{L}(\textbf{y},\bar{\textbf{y}})}{\braket{\bar{\textbf{y}}^{2}}_{L}},$
(11)

where $X$ is the matrix of the reservoir variables at different times, y and
$\bar{\textbf{y}}$ are the prediction and target sequences respectively and w
is the vector of weights of the output layer. The mean square error MSE is the
cost function
$\text{MSE}_{L}(\textbf{y},\bar{\textbf{y}})=\frac{1}{L}\sum^{L}_{k=1}(y_{k}-\bar{y}_{k})^{2}$
and the bracket $\braket{}_{L}$ denotes the temporal average for sequences of
length $L$.

The target function for a given degree $d$ is:

$\bar{y}_{k}=\prod_{i}\mathcal{P}_{d_{i}}[\tilde{s}_{k-i}],\quad\sum_{i}d_{i}=d,$
(12)

where $\tilde{s}_{k}$ accounts for the input and $\mathcal{P}_{d_{i}}$ are the
Legendre polynomials of a given degree $d_{i}$. Due to the choice of
polynomials and to preserve their orthogonality, the input must be uniformly
distributed over the interval $[-1,1]$, rescaling the inputs introduced in the
spin system as $s_{k}=(1+\tilde{s}_{k})/2\in[0,1]$. The numerical simulations
have been done with input streams of length $L=10^{5}$ and maximum degree of
nonlinearity $d_{\text{max}}=6$, being sufficient to guarantee the saturation
of the total capacity in our case. An initial transient of $10^{4}$ inputs was
taken into account to wash out the initial conditions of the system. Since
individual capacities may be slightly overestimated in the numerical analysis,
a threshold is set to truncate the smaller contributions. These thresholds
have been numerically obtained and depend on the input sequence length $L$,
the number of observables in the output layer $O$ and the degree of
nonlinearity of the task (see [49] for more details).

Figure 9: Information processing capacity versus $h$ (a) and $W$ (b). We
represent the average value of the total capacity $C_{T}$ as well as the
different contributions by the maximum degree of nonlinearity of the
polynomials. Solid lines represent the average values while the shadows are
the standard deviation over 10 realizations. Both figures are computed for
$N=10$ and $J_{s}\Delta t=10$, using an output layer of $O=75$ observables
(see main text for details). We took $W=0$ for (a) (transition between regions
III and IV) and $h/J_{s}=10$ for (b) (transition between regions IV and I).
The grey-dashed line corresponds to the intermediate value $\braket{r}=0.46$
in Fig. 1 in the main text, used to guide the eye. The initial condition for
all the realizations is the maximal coherent state
$\rho_{0}=\frac{1}{2^{N}}\sum\ket{i}\bra{j}$.

The results for the normalized IPC are shown in Fig. 9. We represent the total
capacity of the system, denoted by $C_{T}$, as well as the linear ($d=1$) and
nonlinear ($d>1$) memory contributions to the IPC. We find that the localized
regimes (regions I and III) show the smallest values of $C_{T}$, while the
ergodic regime (region IV) shows a stable total capacity close to the maximum
value. We also find that linear and nonlinear memory present a trade-off at
the transitions between the dynamical phases. The vertical dashed lines in
Fig. 9, for an intermediate value of $\braket{r}$, interestingly mark the
building up of the total memory capacity for increasing magnetic field
strength ($h$) and decreasing onsite disorder strength ($W$), respectively.
The low IPC values at the localized regimes can be explained with the
convergence property: the convergence is so slow that we can assume the loss
of the fading memory property for finite input sequences, which is a necessary
condition to saturate the total capacity. This slow convergence rate is a
consequence of the weak contractiveness of the dynamical map in Eq. (4): in
the MBL phase the creation of entanglement between the first qubit and the
rest of the network is very weak (it only grows logarithmically in time), so
that the system cannot efficiently dissipate information about inputs injected
in the distant past.

## References

  * Jaeger [2021] H. Jaeger, Toward a generalized theory comprising digital, neuromorphic, and unconventional computing, Neuromorphic Computing and Engineering (2021).
  * Jaeger [2001] H. Jaeger, The “echo state” approach to analysing and training recurrent neural networks-with an erratum note, Bonn, Germany: German National Research Center for Information Technology GMD Technical Report 148, 13 (2001).
  * Maass  _et al._ [2002] W. Maass, T. Natschläger, and H. Markram, Real-time computing without stable states: A new framework for neural computation based on perturbations, Neural Computation 14, 2531 (2002).
  * Lukoševičius and Jaeger [2009] M. Lukoševičius and H. Jaeger, Reservoir computing approaches to recurrent neural network training, Computer Science Review 3, 127 (2009).
  * Appeltant  _et al._ [2011] L. Appeltant, M. C. Soriano, G. Van der Sande, J. Danckaert, S. Massar, J. Dambre, B. Schrauwen, C. R. Mirasso, and I. Fischer, Information processing using a single dynamical node as complex system, Nature Communications 2, 1 (2011).
  * Van der Sande  _et al._ [2017] G. Van der Sande, D. Brunner, and M. C. Soriano, Advances in photonic reservoir computing, Nanophotonics 6, 561 (2017).
  * Torrejon  _et al._ [2017] J. Torrejon, M. Riou, F. A. Araujo, S. Tsunegi, G. Khalsa, D. Querlioz, P. Bortolotti, V. Cros, K. Yakushiji, A. Fukushima, _et al._ , Neuromorphic computing with nanoscale spintronic oscillators, Nature 547, 428 (2017).
  * Tanaka  _et al._ [2019] G. Tanaka, T. Yamane, J. B. Héroux, R. Nakane, N. Kanazawa, S. Takeda, H. Numata, D. Nakano, and A. Hirose, Recent advances in physical reservoir computing: A review, Neural Networks (2019).
  * Arute  _et al._ [2019] F. Arute, K. Arya, R. Babbush, _et al._ , Quantum supremacy using a programmable superconducting processor, Nature 574, 505 (2019).
  * Fujii and Nakajima [2017] K. Fujii and K. Nakajima, Harnessing disordered-ensemble quantum dynamics for machine learning, Physical Review Applied 8, 024030 (2017).
  * Mujal  _et al._ [2021] P. Mujal, R. Martínez-Peña, J. Nokkala, J. García-Beni, G. L. Giorgi, M. C. Soriano, and R. Zambrini, Opportunities in quantum reservoir computing and extreme learning machines, arXiv preprint arXiv:2102.11831 (2021).
  * Nakajima  _et al._ [2019] K. Nakajima, K. Fujii, M. Negoro, K. Mitarai, and M. Kitagawa, Boosting computational power through spatial multiplexing in quantum reservoir computing, Physical Review Applied 11, 034021 (2019).
  * Chen and Nurdin [2019] J. Chen and H. I. Nurdin, Learning nonlinear input–output maps with dissipative quantum systems, Quantum Information Processing 18, 198 (2019).
  * Tran and Nakajima [2020] Q. H. Tran and K. Nakajima, Higher-order quantum reservoir computing, arXiv preprint arXiv:2006.08999 (2020).
  * Martínez-Peña  _et al._ [2020] R. Martínez-Peña, J. Nokkala, G. L. Giorgi, R. Zambrini, and M. C. Soriano, Information processing capacity of spin-based quantum reservoir computing systems, Cognitive Computation , 1 (2020).
  * Ghosh  _et al._ [2019a] S. Ghosh, A. Opala, M. Matuszewski, T. Paterek, and T. C. Liew, Quantum reservoir processing, npj Quantum Information 5, 35 (2019a).
  * Ghosh  _et al._ [2019b] S. Ghosh, T. Paterek, and T. C. H. Liew, Quantum neuromorphic platform for quantum state preparation, Phys. Rev. Lett. 123, 260404 (2019b).
  * Nokkala  _et al._ [2021] J. Nokkala, R. Martínez-Peña, G. L. Giorgi, V. Parigi, M. C. Soriano, and R. Zambrini, Gaussian states of continuous-variable quantum systems provide universal and versatile reservoir computing, Communications Physics 4, 53 (2021).
  * Govia  _et al._ [2021] L. C. G. Govia, G. J. Ribeill, G. E. Rowlands, H. K. Krovi, and T. A. Ohki, Quantum reservoir computing with a single nonlinear oscillator, Phys. Rev. Research 3, 013077 (2021).
  * Ghosh  _et al._ [2021] S. Ghosh, T. Krisnanda, T. Paterek, and T. C. H. Liew, Realising and compressing quantum circuits with quantum reservoir computing, Communications Physics 4, 105 (2021).
  * Negoro  _et al._ [2018] M. Negoro, K. Mitarai, K. Fujii, K. Nakajima, and M. Kitagawa, Machine learning with controllable quantum dynamics of a nuclear spin ensemble in a solid, arXiv preprint arXiv:1806.10910 (2018).
  * Chen  _et al._ [2020] J. Chen, H. I. Nurdin, and N. Yamamoto, Temporal information processing on noisy quantum computers, Phys. Rev. Applied 14, 024065 (2020).
  * Deutsch [1991] J. M. Deutsch, Quantum statistical mechanics in a closed system, Physical Review A 43, 2046 (1991).
  * Srednicki [1999] M. Srednicki, The approach to thermal equilibrium in quantized chaotic systems, Journal of Physics A: Mathematical and General 32, 1163 (1999).
  * D’Alessio  _et al._ [2016] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics, Advances in Physics 65, 239 (2016).
  * Abanin  _et al._ [2019] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Colloquium: Many-body localization, thermalization, and entanglement, Reviews of Modern Physics 91, 021001 (2019).
  * Žunkovič  _et al._ [2018] B. Žunkovič, M. Heyl, M. Knap, and A. Silva, Dynamical quantum phase transitions in spin chains with long-range interactions: Merging different concepts of nonequilibrium criticality, Physical review letters 120, 130601 (2018).
  * Huse  _et al._ [2013] D. A. Huse, R. Nandkishore, V. Oganesyan, A. Pal, and S. L. Sondhi, Localization-protected quantum order, Physical Review B 88, 014206 (2013).
  * Ponte  _et al._ [2015] P. Ponte, A. Chandran, Z. Papić, and D. A. Abanin, Periodically driven ergodic and many-body localized quantum systems, Annals of Physics 353, 196 (2015).
  * Lazarides  _et al._ [2015] A. Lazarides, A. Das, and R. Moessner, Fate of many-body localization under periodic driving, Physical review letters 115, 030402 (2015).
  * Tangpanitanon  _et al._ [2020] J. Tangpanitanon, S. Thanasilp, N. Dangniam, M.-A. Lemonde, and D. G. Angelakis, Expressibility and trainability of parametrized analog quantum systems for machine learning applications, Phys. Rev. Research 2, 043364 (2020).
  * Altshuler  _et al._ [2010] B. Altshuler, H. Krovi, and J. Roland, Anderson localization makes adiabatic quantum optimization fail, Proceedings of the National Academy of Sciences 107, 12446 (2010).
  * Laumann  _et al._ [2015] C. R. Laumann, R. Moessner, A. Scardicchio, and S. L. Sondhi, Quantum annealing: The fastest route to quantum computation?, The European Physical Journal Special Topics 224, 75 (2015).
  * Keating  _et al._ [2007] J. Keating, N. Linden, J. Matthews, and A. Winter, Localization and its consequences for quantum walk algorithms and quantum communication, Physical Review A 76, 012315 (2007).
  * Schreiber  _et al._ [2011] A. Schreiber, K. Cassemiro, V. Potoček, A. Gábris, I. Jex, and C. Silberhorn, Decoherence and disorder in quantum walks: from ballistic spread to localization, Physical review letters 106, 180403 (2011).
  * Oganesyan and Huse [2007] V. Oganesyan and D. A. Huse, Localization of interacting fermions at high temperature, Physical Review B 75, 155111 (2007).
  * [37] See Supplemental Material for additional notes on the dynamics, the convergence properties of the spin network, a note on the experimental implementations, the NARMA task, and the IPC of the system, including Refs. [64, 65, 66, 67, 68, 69, 70, 71, 72, 73].
  * Atas  _et al._ [2013] Y. Atas, E. Bogomolny, O. Giraud, and G. Roux, Distribution of the ratio of consecutive level spacings in random matrix ensembles, Physical Review Letters 110, 084101 (2013).
  * Baldwin  _et al._ [2017] C. Baldwin, C. Laumann, A. Pal, and A. Scardicchio, Clustering of nonergodic eigenstates in quantum spin glasses, Physical Review Letters 118, 127201 (2017).
  * Mukherjee  _et al._ [2018] S. Mukherjee, S. Nag, and A. Garg, Many-body localization-delocalization transition in the quantum sherrington-kirkpatrick model, Physical Review B 97, 144202 (2018).
  * Rademaker and Abanin [2020] L. Rademaker and D. A. Abanin, Slow nonthermalizing dynamics in a quantum spin glass, Phys. Rev. Lett. 125, 260405 (2020).
  * Maksymov and Burin [2020] A. O. Maksymov and A. L. Burin, Many-body localization in spin chains with long-range transverse interactions: Scaling of critical disorder with system size, Physical Review B 101, 024201 (2020).
  * Smith  _et al._ [2016] J. Smith, A. Lee, P. Richerme, B. Neyenhuis, P. W. Hess, P. Hauke, M. Heyl, D. A. Huse, and C. Monroe, Many-body localization in a quantum simulator with programmable random disorder, Nature Physics 12, 907 (2016).
  * Moudgalya  _et al._ [2020] S. Moudgalya, D. A. Huse, and V. Khemani, Perturbative instability towards delocalization at phase transitions between MBL phases, arXiv preprint arXiv:2008.09113 (2020).
  * Sahay  _et al._ [2021] R. Sahay, F. Machado, B. Ye, C. R. Laumann, and N. Y. Yao, Emergent ergodicity at the transition between many-body localized phases, Phys. Rev. Lett. 126, 100604 (2021).
  * Khemani  _et al._ [2017] V. Khemani, S. P. Lim, D. N. Sheng, and D. A. Huse, Critical properties of the many-body localization transition, Phys. Rev. X 7, 021013 (2017).
  * Panda  _et al._ [2020] R. K. Panda, A. Scardicchio, M. Schulz, S. R. Taylor, and M. Žnidarič, Can we study the many-body localisation transition?, EPL (Europhysics Letters) 128, 67003 (2020).
  * Grigoryeva and Ortega [2018] L. Grigoryeva and J.-P. Ortega, Echo state networks are universal, Neural Networks 108, 495 (2018).
  * Dambre  _et al._ [2012] J. Dambre, D. Verstraeten, B. Schrauwen, and S. Massar, Information processing capacity of dynamical systems, Scientific reports 2, 1 (2012).
  * Atiya and Parlos [2000] A. F. Atiya and A. G. Parlos, New results on recurrent network training: unifying the algorithms and accelerating convergence, IEEE transactions on neural networks 11, 697 (2000).
  * Carroll [2020] T. L. Carroll, Do reservoir computers work best at the edge of chaos?, Chaos 30, 121109 (2020).
  * [52] The localized regime shows an extremely much slower convergence than in the ergodic phase. For practical purposes, this is definitely an obstacle for implementing QRC even if strictly speaking convergence would asymptotically be achieved since the map presented in Eq. (2) of the manuscript is contractive due to the partial trace.
  * Swingle [2018] B. Swingle, Unscrambling the physics of out-of-time-order correlators, Nature Physics 14, 988 (2018).
  * Heyl  _et al._ [2013] M. Heyl, A. Polkovnikov, and S. Kehrein, Dynamical quantum phase transitions in the transverse-field ising model, Phys. Rev. Lett. 110, 135704 (2013).
  * Eddins  _et al._ [2019] A. Eddins, J. Kreikebaum, D. Toyli, E. Levenson-Falk, A. Dove, W. Livingston, B. Levitan, L. Govia, A. Clerk, and I. Siddiqi, High-efficiency measurement of an artificial atom embedded in a parametric amplifier, Physical Review X 9, 011004 (2019).
  * Zhang  _et al._ [2017] J. Zhang, G. Pagano, P. W. Hess, A. Kyprianidis, P. Becker, H. Kaplan, A. V. Gorshkov, Z.-X. Gong, and C. Monroe, Observation of a many-body dynamical phase transition with a 53-qubit quantum simulator, Nature 551, 601 (2017).
  * Kaplan  _et al._ [2020] H. B. Kaplan, L. Guo, W. L. Tan, A. De, F. Marquardt, G. Pagano, and C. Monroe, Many-body dephasing in a trapped-ion quantum simulator, Phys. Rev. Lett. 125, 120605 (2020).
  * Choi  _et al._ [2016] J.-y. Choi, S. Hild, J. Zeiher, P. Schauß, A. Rubio-Abadal, T. Yefsah, V. Khemani, D. A. Huse, I. Bloch, and C. Gross, Exploring the many-body localization transition in two dimensions, Science 352, 1547 (2016).
  * Kaufman  _et al._ [2016] A. M. Kaufman, M. E. Tai, A. Lukin, M. Rispoli, R. Schittko, P. M. Preiss, and M. Greiner, Quantum thermalization through entanglement in an isolated many-body system, Science 353, 794 (2016).
  * Houck  _et al._ [2012] A. A. Houck, H. E. Türeci, and J. Koch, On-chip quantum simulation with superconducting circuits, Nature Physics 8, 292 (2012).
  * Xu  _et al._ [2018] K. Xu, J.-J. Chen, Y. Zeng, Y.-R. Zhang, C. Song, W. Liu, Q. Guo, P. Zhang, D. Xu, H. Deng, K. Huang, H. Wang, X. Zhu, D. Zheng, and H. Fan, Emulating many-body localization with a superconducting quantum processor, Phys. Rev. Lett. 120, 050507 (2018).
  * Pierangeli  _et al._ [2019] D. Pierangeli, G. Marcucci, and C. Conti, Large-scale photonic ising machine by spatial light modulation, Physical review letters 122, 213902 (2019).
  * Pierangeli  _et al._ [2021] D. Pierangeli, M. Rafayelyan, C. Conti, and S. Gigan, Scalable spin-glass optical simulator, Phys. Rev. Applied 15, 034087 (2021).
  * Weinberg and Bukov [2017] P. Weinberg and M. Bukov, QuSpin: a Python Package for Dynamics and Exact Diagonalisation of Quantum Many Body Systems part I: spin chains, SciPost Phys. 2, 003 (2017).
  * Johansson  _et al._ [2012] J. R. Johansson, P. D. Nation, and F. Nori, Qutip: An open-source python framework for the dynamics of open quantum systems, Computer Physics Communications 183, 1760 (2012).
  * Ilievski  _et al._ [2015] E. Ilievski, J. De Nardis, B. Wouters, J.-S. Caux, F. H. Essler, and T. Prosen, Complete generalized Gibbs ensembles in an interacting theory, Physical Review Letters 115, 157201 (2015).
  * Vosk and Altman [2014] R. Vosk and E. Altman, Dynamical quantum phase transitions in random spin chains, Physical Review Letters 112, 217204 (2014).
  * Geraedts  _et al._ [2017] S. D. Geraedts, R. N. Bhatt, and R. Nandkishore, Emergent local integrals of motion without a complete set of localized eigenstates, Physical Review B 95, 064204 (2017).
  * Zobov and Lundin [2017] V. E. Zobov and A. A. Lundin, On the effect of an inhomogeneous magnetic field and many-body localization on the increase in the second moment of multiple-quantum nmr with time, JETP Letters 105, 514 (2017).
  * Wei  _et al._ [2018] K. X. Wei, C. Ramanathan, and P. Cappellaro, Exploring localization in nuclear spin chains, Physical Review Letters 120, 070501 (2018).
  * Richerme  _et al._ [2014] P. Richerme, Z.-X. Gong, A. Lee, C. Senko, J. Smith, M. Foss-Feig, S. Michalakis, A. V. Gorshkov, and C. Monroe, Non-local propagation of correlations in quantum systems with long-range interactions, Nature 511, 198 (2014).
  * Huang  _et al._ [2020] H.-Y. Huang, R. Kueng, and J. Preskill, Predicting many properties of a quantum system from very few measurements, Nature Physics 16, 1050 (2020).
  * Kubota  _et al._ [2019] T. Kubota, K. Nakajima, and H. Takahashi, Dynamical anatomy of narma10 benchmark task, arXiv preprint arXiv:1906.04608 (2019).