Edge of Many-Body Quantum Chaos in Quantum Reservoir Computing

Introduction. —Quantum mechanics harbors immense potential to revolutionize modern technology [1]. Breakthroughs in quantum computing, communication, and sensing have vividly showcased the transformative power of quantum resources [2, 3, 4, 5]. These advances have, in turn, inspired extensive research efforts in quantum machine learning, with various approaches—notably variational quantum algorithms (VQAs)—being actively explored [6, 7, 8, 9, 10, 11].

More recently, the power of quantum has also integrated with reservoir computing (RC), a machine learning framework renowned for its efficacy in time-series analysis [12, 13, 14, 15, 16]. True to its name, the most distinctive feature of RC is the use of an input-driven dynamical system, known as a “reservoir”, to nonlinearly map input data onto a feature space. Notably, training in RC is limited to the post-processing stage, and the reservoir itself remains fixed. Extending this idea into the quantum domain, Fujii and Nakajima introduced quantum reservoir computing (QRC) [17]. QRC harnesses natural quantum dynamics as a reservoir, generating rich, nonlinear mappings in an exponentially large Hilbert space without encountering the trainability bottlenecks prevalent in VQAs. This capability has sparked significant theoretical progress [18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36], as well as promising experimental demonstrations [37, 38, 39, 40, 41, 42, 43, 44].

Given the optimization-free nature of RC, researchers naturally ask: What characterizes an effective reservoir system? While this question remains largely unanswered for QRC, several key criteria have been identified for classical RC. The most renowned among them is the concept of the “edge of chaos” [45, 46, 47, 48, 49]. This principle posits that the boundary between chaotic and nonchaotic dynamical regimes provides a performance sweet spot for RC systems, as empirically verified across diverse platforms [50, 51, 52]. By analogy, exploring a quantum version of the edge of chaos offers a promising route to addressing the fundamental question posed above.

In classical systems, chaos is naturally framed via phase space trajectories [53]. In quantum mechanics, however, the Heisenberg uncertainty principle precludes such a direct picture, and one must instead rely on more subtle diagnostics. Early investigations of quantum chaos therefore focused on single-particle systems whose classical limits are chaotic, leading to the conjecture (with some analytical supports) that their spectral statistics are governed by random matrix theory (RMT) [54, 55, 56, 57, 58]. In recent years, attention has shifted to many-body settings, spurred in large part by the exactly solvable Sachdev-Ye-Kitaev (SYK) model [59, 60, 61, 62]. Various diagnostics, including level statistics [63, 64, 65, 66], out-of-time-order correlators [67, 68, 69], and spectral form factors [70, 71, 72, 73] have enriched our understanding of chaos in interacting quantum systems. Experimental realizations of the SYK physics have also been proposed on several platforms [74, 75, 76, 77, 78]. Notably, the SYK model exhibits maximally quantum chaotic behavior, saturating the so-called bound on chaos [79, 69, 80], and has thus stood as a canonical model for the investigation of many-body quantum chaos. These insights motivate the exploration of QRC built upon the SYK model, aiming to address the open question of reservoir effectiveness from a quantum chaotic standpoint.

In this Letter, we explicitly demonstrate key characteristics of effective quantum reservoir systems. In particular, we examine the QRC performance near two distinct chaotic boundaries of the SYK model. The first is a temporal boundary defined by the Thouless time, beyond which the system’s dynamics conforms to RMT. The second is a parametric boundary marking the transition from quantum chaotic to integrable (nonchaotic) behavior, controlled by the strength of an additional noninteracting term. Using representative machine learning tasks that demand both memory retention and nonlinear transformation, we demonstrate that QRC performance is maximized near the onset of quantum chaotic regime in both temporal and parametric domains. These observations establish the “edge of many-body quantum chaos” as a foundational design principle for QRC systems.

Model. —Figure 1 schematically illustrates the QRC architecture implemented with the SYK model for fermions. The quantum reservoir system, whose dynamics is harnessed for information processing, is governed by the Hamiltonian

where $N$ denotes the system size, and $c^{\dagger}_{i}$ ($c_{i}$) is the spinless fermionic creation (annihilation) operator at site $i$. These operators satisfy the anticommutation relations $\{c_{i},c_{j}^{\dagger}\}=\delta_{i,j}$ and $\{c_{i},c_{j}\}=\{c_{i}^{\dagger},c_{j}^{\dagger}\}=0$. In Eq. (1), the first term describes the SYK₄ interaction: the all-to-all couplings $J_{ijkl}$ are drawn from a zero-mean complex Gaussian distribution of variance $\overline{|J_{ijkl}|^{2}}=J_{4}^{2}/N^{3}$ ($\overline{\bullet}$ denotes the ensemble average). The second term is the noninteracting SYK₂ contribution, with coupling constants $\kappa_{ij}$ sampled from a zero-mean complex Gaussian distribution of variance $\overline{|\kappa_{ij}|^{2}}=\kappa_{2}^{2}/(2N)$. The couplings satisfy the symmetry relations $\kappa_{ij}=\kappa_{ji}^{*}$, $J_{ijkl}=J_{klij}^{*}$, and $J_{ijkl}=-J_{jikl}=-J_{ijlk}$. Both coupling terms conserve the total fermion number $N_{p}\equiv\sum_{i}\langle c_{i}^{\dagger}c_{i}\rangle$. Unless otherwise specified, we set $N=8$, $J_{4}=1$, and $N_{p}=N/2$ (note that the QRC input protocol introduced below does not conserve $N_{p}$).

To diagnose the chaotic nature of quantum many-body systems, we first examine the level statistics [63, 64, 65, 66]. Let the ordered eigenvalues of $\mathcal{H}$ be $E_{n}$; the spacings between consecutive energy levels are given by $s_{n}=E_{n+1}-E_{n}$. In an integrable system, the spacing distribution follows Poisson statistics, whereas in a quantum chaotic system, it is well described by Wigner-Dyson (WD) distribution. A convenient metric to characterize these distributions is the level spacing ratio, defined as $r_{n}=\min(s_{n}/s_{n+1},s_{n+1}/s_{n})$. Its mean value $\langle r\rangle$ assumes $\approx 0.3862$ for Poisson distribution and $\approx 0.5996$ for WD statistics [63, 64]. Figure 2(a) presents the distribution of the level spacing ratio $P(r)$ for the SYK₄ model [$(J_{4},\kappa_{2})=(1,0)$] and the SYK₂ model [$(J_{4},\kappa_{2})=(0,1)$], in excellent agreement with the WD and Poisson predictions, respectively. Notably, when both terms are present, the system can be tuned between these extremes: as the ratio of the typical coupling strengths $\kappa_{2}/J_{4}$ increases, Fig. 2(b) displays a gradual shift of $\langle r\rangle$ from the WD to the Poisson value. Despite the smearing effects associated with finite system sizes, this result strongly suggests the emergence of the quantum chaotic to integrable transition [81, 82, 83, 84, 85] (see the Supplemental Material [86]).

Framework of the QRC. —We now describe the QRC framework for information processing; all procedures employed herein follow the conventional protocols established in earlier investigations [17, 18]. The QRC architecture comprises four principal stages: signal encoding, reservoir dynamics, read-out, and post-processing, through which the input at the $k$-th step $u^{(k)}$ is processed in an attempt to produce its corresponding target value $\overline{y}^{(k)}$ (Fig. 1). Consider a discretized input sequence $\{u^{(k)}\}\equiv\{u^{(0)},u^{(1)},\dots\}$ with $u^{(k)}\in[0,1]$. Each input $u^{(k)}$ is encoded into the reservoir by replacing the quantum state at site $1$ with $\rho_{\mathrm{in}}^{(k)}\equiv|\psi_{\mathrm{in}}^{(k)}\rangle\langle\psi_{\mathrm{in}}^{(k)}|$, where $|\psi_{\mathrm{in}}^{(k)}\rangle\equiv\sqrt{1-u^{(k)}}|0\rangle+\sqrt{u^{(k)}}|1\rangle$. Subsequently, the system evolves unitarily under the Hamiltonian in Eq. (1). The encoding operation is repeated at intervals $\Delta t_{\mathrm{in}}$, which determines the intrinsic time-scale of the QRC system. The overall update rule is summarized by the completely positive trace-preserving map:

where $\mathrm{Tr}_{1}\{\bullet\}$ denotes the partial trace over site $1$; see the Supplemental Material [86] for details on the observed dynamics. During the evolution, the expectation value of the single-site particle number is measured at discrete times $t_{v}=v\Delta t_{\mathrm{in}}/V$, where $0<v\leq V$ and $V$ is the number of virtual computational nodes (here, $V=10$). At the $k$-th step, the outcomes $\mathrm{Tr}\{\rho^{(k)}(t_{v})c_{i}^{\dagger}c_{i}\}$ are assembled into a one-dimensional state vector $\bm{\mathrm{x}}^{(k)}$; including an additional constant term, $\bm{\mathrm{x}}^{(k)}$ has $NV+1$ components in total. The output layer then applies a linear transformation $y^{(k)}=W_{\mathrm{out}}\bm{\mathrm{x}}^{(k)}$, where the output weight $W_{\mathrm{out}}$ is trained by minimizing the mean squared error between the computed output $y^{(k)}$ and the target output $\bar{y}^{(k)}$. In our simulations, a total of $10000$ inputs are processed. Starting from a random initial state, the first $4000$ inputs are discarded to wash out the effects of initial conditions; for the convergence property, see the Supplemental Material [86]. The subsequent $3000$ inputs are employed for training $W_{\mathrm{out}}$, and the final $3000$ are reserved for evaluation. Performance results are averaged over $500$ independent random realizations.

We evaluate the reservoir performance on two complementary tasks: the short-term memory (STM) task, which probes linear memory capacity, and the nonlinear auto-regressive moving average (NARMA) task, which evaluates combined memory and nonlinear processing capabilities [14, 17].

In the STM task, the target output at delay $d$ is defined as the input value $d$ steps earlier in a random input sequence $\{u^{(k)}\}$, i.e., $\bar{y}^{(k)}_{d}=u^{(k-d)}$. During the testing phase, we collect the outputs $\bm{\mathrm{y}}=\{y^{(k)}\}$ and targets $\bar{\bm{\mathrm{y}}}_{d}=\{\bar{y}^{(k)}_{d}\}$, and compute the determination coefficient $R^{2}_{d}\equiv\mathrm{cov}^{2}(\bm{\mathrm{y}},\bar{\bm{\mathrm{y}}}_{d})/[\sigma^{2}(\bm{\mathrm{y}})\sigma^{2}(\bar{\bm{\mathrm{y}}}_{d})]$, where $\mathrm{cov}$ and $\sigma^{2}$ denote the covariance and the variance, respectively. $R^{2}_{d}$ is bounded between $0$ and $1$; values closer to $1$ indicate better memory performance.

In the NARMA task, the target sequence at a given order $n$ is defined by the recursive relation: $\bar{y}^{(k+1)}=\alpha\bar{y}^{(k)}+\beta\bar{y}^{(k)}\sum_{j=0}^{n-1}\bar{y}^{(k-j)}+\gamma u^{(k-n+1)}u^{(k)}+\delta$, where ($\alpha$, $\beta$, $\gamma$, $\delta$) are set to ($0.3$, $0.05$, $1.5$, $0.1$) [17, 20]. The input stream is sampled from the range $[0,0.2]$ to prevent divergence in $\{\bar{y}^{(k)}\}$, which is then rescaled to $[0,1]$ when encoded to the quantum reservoir [17]. Performance is evaluated in terms of the deviation between the output and the target, which is quantified by the normalized mean-squared error (NMSE) $\|\bm{\mathrm{y}}-\bar{\bm{\mathrm{y}}}\|^{2}/\|\bar{\bm{\mathrm{y}}}\|^{2}$ with $\|\bullet\|$ being the Euclidean norm. A lower NMSE indicates superior performance in accurately reproducing the target nonlinear dynamics.

Temporal edge of many-body quantum chaos. —Our investigation begins with an examination of the QRC performance of the SYK₄ and SYK₂ models, in which clear fingerprints of many-body quantum chaos manifest themselves in the temporal domain. As a time-domain probe of the energy-spectrum correlations characteristic of quantum chaos, we analyze the spectral form factor (SFF) $K(t)=\overline{\sum_{m,n}e^{i(E_{m}-E_{n})t}}/\mathcal{N}^{2}$, where $\mathcal{N}$ denotes the Hilbert space dimension [70, 71, 72, 73]. In many-body chaotic systems, the SFF generally exhibits a “ramp-plateau” profile governed by two time-scales: the Thouless time $t_{\mathrm{Th}}$, after which the evolution shows a universal “ramp” and follows RMT predictions, and the plateau time $t_{\mathrm{p}}$, beyond which $K(t)$ saturates at a plateau. Figures 3(a) and 3(b) illustrate the SFF for the SYK₄ and SYK₂ limits, respectively. Notably, in the SYK₄ case, the SFF exhibits the pronounced ramp beginning around $t_{\mathrm{Th}}\sim 4.5$ and the subsequent plateau beyond $t_{\mathrm{p}}\sim 80$, providing a clear signature of many-body quantum chaos. While $K(t)$ in the SYK₂ model also features a dip and a plateau, they occur on a much shorter timescale and with a different system-size scaling of $t_{\mathrm{p}}$ compared to the chaotic system (see the Supplemental Material for details [86]).

Turning to computational benchmarks, Figs. 3(c) and 3(d) plot the QRC performance of the SYK₄ and SYK₂ models as functions of the input interval $\Delta t_{\mathrm{in}}$, which defines the intrinsic time-scale governing the QRC dynamics. The upper panels display the STM performance $R^{2}_{d=n-1}$ at delay $d=n-1$, and the lower panels show the NMSE on the NARMA-$n$ tasks of orders $n=2$, $3$, $5$, and $7$. The chosen STM delay corresponds to the oldest input term explicitly appearing in the recursive definition of the NARMA-$n$ dynamics; additional past inputs are nevertheless implicitly encoded within the target sequence $\{\bar{y}^{k}\}$ (see the Supplemental Material for the full delay dependence of $R^{2}_{d}$ [86]). To establish a baseline for comparison, we also analyze the Haar QRC model, in which the evolution operator $e^{-i\mathcal{H}t}$ is replaced by a unitary $U$ drawn from the Haar measure [87, 88]. Its performance, shown as horizontal dashed lines, deteriorates monotonically with increasing $n$ as the task complexity grows.

A striking pattern emerges in the QRC performance for the SYK₄ model [Fig. 3(c)]. When $\Delta t_{\mathrm{in}}$ is very small, the time evolution of the quantum reservoir does not yield an effective feature map, resulting in poor performance on both the STM and NARMA tasks. However, as $\Delta t_{\mathrm{in}}$ increases toward the Thouless time $t_{\mathrm{Th}}$, the QRC performance improves dramatically and reaches its maximum just before $t_{\mathrm{Th}}$. Beyond this point, both metrics level off at the baseline of the Haar QRC model for each order $n$, in accordance with the onset of universal RMT descriptions typical of quantum chaotic systems. In other words, the Thouless time effectively marks the temporal boundary where fully quantum chaotic signatures begin to dominate the QRC performance. The observed peak, situated just prior to $t_{\mathrm{Th}}$, thus explicitly underscores the performance enhancement at the “edge of many-body quantum chaos” in the temporal domain.

In contrast, the QRC with the integrable SYK₂ model behaves quite differently [Fig. 3(d)]. As $\Delta t_{\mathrm{in}}$ increases, the NARMA performance improves only modestly from its initial low values and quickly saturates. No pronounced performance peak emerges therein. This behavior starkly contrasts with that of the SYK₄ model, indicating that the performance enhancements observed in Fig. 3(c) are fundamentally attributable to the latter’s intrinsic many-body quantum chaotic nature. Notably, despite its poor NARMA performance, the STM metric for the SYK₂ model remains close to $1$ for sufficiently large values of $\Delta t_{\mathrm{in}}$, surpassing both the SYK₄ model and the Haar QRC reference. This exemplifies the well-established trade-off in RC between memory capacity and nonlinear transformation capability [89, 90, 91]: QRC with the integrable SYK₂ model retains strong memory but lacks nonlinear richness required for complex tasks. These model-dependent behaviors motivate our exploration of interpolations between the SYK₂ and SYK₄ models.

Parametric edge of many-body quantum chaos. —We here demonstrate the QRC performance across the chaotic-integrable transition by mixing the SYK₄ and SYK₂ Hamiltonians. To ensure a fair comparison as the coupling strength is varied, each Hamiltonian is rescaled by its spectral norm: $\mathcal{H}\rightarrow\mathcal{H}/||\mathcal{H}||$ [92]. (Note that this correspondingly changes the time-scale relative to Fig. 3.) Figures 4(a) and 4(b) display the STM and NARMA performances versus the ratio $\kappa_{2}/J_{4}$ for long ($\Delta t_{\mathrm{in}}=50$) and short ($\Delta t_{\mathrm{in}}=1$) input intervals, respectively. The vertical dashed lines serve as a guide of the eye for the onset of WD and Poisson statistics, while the horizontal dotted lines represent the performance of the Haar QRC model for reference. We note that these input intervals are chosen to probe two distinct time-scales: the longer input interval $\Delta t_{\mathrm{in}}=50$ evaluates QRC performance within the universal RMT regime, while the shorter interval $\Delta t_{\mathrm{in}}=1$ characterizes QRC behavior in the nonuniversal, short-time window outside the scope of RMT (see the Supplemental Material [86]).

The most pronounced trend is observed for $\Delta t_{\mathrm{in}}=50$ [Fig. 4(a)]. When $\kappa_{2}/J_{4}$ is varied, the NARMA performance exhibits a clear peak at a coupling strength just before the onset of the many-body quantum chaotic regime. This peak provides another manifestation of the “edge of many-body quantum chaos”, in this case in the parameter domain. By contrast, for the short input interval $\Delta t_{\mathrm{in}}=1$, varying the coupling strength yields no performance enhancement [Fig. 4(b)]. At such a short time-scale, the QRC does not reflect chaotic characteristics, so the performance changes only gradually, without any distinct peak at the chaos boundary.

We remark that the specifics of performance enhancement are task-dependent. For instance, in the higher-order tasks detailed in the Supplementary Material, both the position and strength of performance peaks vary with the order [86]. Nevertheless, the fundamental trend remains robust: the most pronounced improvements appear at the edge of many-body quantum chaos. Moreover, since the NARMA task captures nonlinear temporal dynamics that mirror those in real-world systems, the performance boosts shown on this task strongly suggest that similar improvements will extend to a broad class of practical temporal processing applications.

Discussions. —Recent advances in QRC have yielded related performance enhancements, and we believe that a concise comparison highlights the significance of our contributions. One of the authors and his collaborators identified a performance peak near the operational regime boundary in their feedback-driven QRC protocol [32]; however, this enhancement appears predominantly due to classical instabilities in feedback connections, rather than intrinsically quantum phenomena. In the purely quantum regime, Martínez-Peña et al. demonstrated enhanced performance near the ergodic-localization transition in the spin-based QRC [20, 23, 33]. While this finding highlights the importance of dynamical phases in QRC, the role of many-body quantum chaos, and its relation to RMT, remain largely unexplored. Recently, Llodrà et al. analyzed QRC performance in connection with quantum chaos at the quantum phase transition between superfluid and Mott insulator in a clean system, which is complementary to our study on random systems [34]. From a standpoint of temporal domain, investigations into quantum chaos remain sparse: Vetrano et al. explored the role of scrambling times, though within a different framework (quantum extreme learning machines) that does not consider memory of sequential inputs [88].

Building upon these studies, our work clarifies the key characteristics that endow QRC with high computational performance. We have extended the concept of the “edge of chaos” into the quantum realm, investigating how many-body quantum chaos influences QRC in both temporal and parametric domains. The SYK model—a canonical many-body quantum chaotic system—serves as an ideal testbed. Through a unified analysis, we have revealed two distinct optimal regimes: one is the “temporal edge”, which emerges before the Thouless time, and the other is the “parametric edge”, which appears near the quantum chaotic-to-integrable transition. These findings integrate existing knowledge and establish a general principle: optimal QRC performance occurs at the “edge of many-body quantum chaos”. Our study therefore provides clear design principles for engineering quantum reservoirs, unlocking advanced information processing with quantum many-body systems.

Acknowledgments. —The authors thank Shusei Wadashima for fruitful discussions. This work was supported by the JSPS KAKENHI (No. JP25H01247). K.K. was supported by the Program for Leading Graduate Schools (MERIT-WINGS) and JST BOOST (No. JPMJBS2418). Part of computation in this work has been done using the facilities of the Supercomputer Center, the Institute for Solid State Physics, the University of Tokyo.