{'fraction_non_alphanumeric': 0.04642929395104188, 'fraction_numerical': 0.015375278171151122, 'mean_word_length': 4.2019205472244145, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Higher-Order Quantum Reservoir Computing

Quantum reservoir computing (QRC) is an emerging paradigm for harnessing the natural dynamics of quantum systems as computational resources that can be used for temporal machine learning (ML) tasks. In the current setup, QRC is difficult to deal with high-dimensional data and has a major drawback of scalability in physical implementations. We propose higher-order QRC, a hybrid quantum-classical framework consisting of multiple but small quantum systems that are mutually communicated via classical connections like linear feedback. By utilizing the advantages of both classical and quantum techniques, our framework enables an efficient implementation to boost the scalability and performance of QRC. We demonstrate the effectiveness of our framework in emulating large-scale nonlinear dynamical systems, including complex spatiotemporal chaos, which outperforms many of the existing ML techniques in certain situations.


Introduction
It is postulated that quantum computers may outperform classical computers when it comes to ML tasks due to the superior ability of quantum mechanics to generate counter-intuitive patterns [1]. Quantum machine learning (QML) is an active interdisciplinary research area proposed from this motivation to improve existing ML methods through the advantages of quantum mechanics [2]. For the foreseeable future, QML algorithms are expected to run on noisy intermediate-scale quantum (NISQ) devices, which includes a few tens of qubits and supports only non-error corrected computations [3]. Existing quantum techniques are utilized to their fullest extent in ML tasks on these devices via hybrid quantum-classical methods that combine classical learning regimes with the advantages of quantum systems [4,5].
Reservoir computing (RC) [6][7][8][9] is a framework that originated from recurrent neural networks to efficiently solve temporal ML tasks. Conventional RC consists of a randomly connected network called a reservoir and a trainable readout part for pattern analysis from output states of the reservoir. The input stream is fed into the reservoir, which functions like a nonlinear processing unit to project low-dimensional input into a high-dimensional dynamical system. Since the training in RC is simple and extremely fast, it is highly suitable and amendable for hardware implementation in a wide variety of physical systems [10][11][12][13][14][15][16]. Quantum reservoir computing (QRC) is a variation of RC, where the reservoir is implemented as a quantum many-body system such as a set of interacting qubits [17,18] or a set of fermions [19,20] driven by a Hamiltonian dynamics. The random connections in the classical reservoir are replaced by basic quantum tunneling. The input stream then drives the transition state, which is evolved through a unitary operator based on the dynamics of the system. Thereafter, quantum measurements are performed to obtain signals that can be considered as reservoir states for training [ Fig. 1(a)]. Recently, the nuclear-magnetic-resonance (NMR) spin-ensemble in a molecular solid [21] and the quantum circuits on superconducting quantum processors [22] have been reported as the physical implementations of QRC. However, aspects such as the difficulty in using multidimensional inputs or the operating limitations when increasing the number of qubits remain uncertain. Figure 1: Schematic of the QRC and the higher-order quantum reservoir (HQR). (a) The input stream drives the transition of the state of qubits, and the quantum system continues evolving itself for a specific time interval to spread the information of the input. Measurements are performed on some qubits, which correspond to signals of readout nodes. After collecting a sufficient amount of input-output signal pairs, the readout weights are trained to emulate the target signals. (b) The HQR consists of multiple different quantum reservoirs that are mutually connected via classical connections to enable multidimensional inputs and enhance both quantum and classical advantages.
In this paper, we propose a general hybrid quantum-classical framework for RC, namely, higher-order quantum reservoir computing (HQRC). We aim to utilize ensemble small quantum systems as a big reservoir and enhance its computational power through random connections between sub-systems [ Fig. 1(b)]. Our scheme is called "higher-order" in terms of different quantum systems placed in an ensemble reservoir, where each system can be considered as a node in the reservoir. Multidimensional inputs can be easily fed into the system by injecting each element to each sub-system. The input of each sub-system is a linear combination of the common input streams and the signals obtained from other systems. This scheme enables to equip a massive amount of computational nodes with the controllable linear feedback, which is expected to increase the expressive power of the system.

Contribution.
We provide a detailed procedure for the temporal processing in HQRC with theoretical explanations of parameter design, which are useful for evaluating a quantum reservoir (QR) as a learning system. We numerically verify that the computational power of QRC, like the memory capacity can be enhanced through higher-order settings. We further demonstrate the effectiveness of our proposal in emulating nonlinear dynamical systems, including high-dimensional spatiotemporal chaos. Interestingly, the experimental results indicate that our approach performs comparatively with classical approaches and even outperforms them in certain situations, such as in the limitation of the training data. From the physical implementation perspective, our framework also paves several opportunities for effective design of recently proposed experimental platforms for QRC [21,22].

Preliminaries
We describe some standard notions that are necessary in our analysis. First, we define M m×n and R m×n as the sets of m×n-dimensional matrix with complex elements and real elements, respectively. For U ∈ M m×n , we denote U † ∈ M n×m as the conjugate transpose of U . A complex square matrix U is called a unitary matrix if U U † = U † U is the identity matrix, and a Hermitian matrix if U = U † .
Pure and mixed states. A pure state |ψ of a spin is a two-dimensional complex vector spanned by the eigenstates {|0 , |1 } of the Pauli operator σ z = 1 0 0 −1 . We denote ψ| = (|ψ ) † ; therefore, . Given a quantum system S comprising N qubits, the Hilbert space H S of S is a tensor product space of two-dimensional individual spin spaces. A pure quantum state is represented by a normalized 2 N -dimensional vector |Ψ . Generally, the quantum state can be considered a statistical mixture of pure states, which is described by a density matrix ρ ∈ M 2 N ×2 N . For a pure state |Ψ , the density matrix is defined as ρ = |Ψ Ψ|.
If we consider a quantum state of the system that is sampled from a set of pure states {|Ψ k } with a probability distribution {p k }, the density matrix is given by ρ = k p k |Ψ k Ψ k |. The density matrix ρ is Hermitian and positive semidefinite; moreover, Tr(ρ) = 1, where Tr(X) denotes the trace of square matrix X.
Partial trace. Given a composite system S of two subsystems S 1 and S 2 , the Hilbert space of S is the tensor product H S1 ⊗ H S2 . From the density matrix ρ of S, we can recover the marginal "reduced density matrix" for a subsystem via the partial trace operation Tr 1 (ρ), which is defined as the linear extension of the mapping Tr 1 : Observables and time evolution. Observables of quantum system S correspond to Hermitian matrices O ∈ M 2 N ×2 N . If we perform measurements on these observables at state ρ, it means that S is interacting with measuring apparatus in the presence of the surrounding environment. The outcomes of the measurement are recorded on the dial on the measuring apparatus, and the statistics of measurement outcomes is determined by the expected value of O as O = Tr[ρO]. If S is a closed system, the time evolution is generated by a Hamiltonian H via the Schrödinger equation where H ∈ M 2 N ×2 N is the Hermitian matrix that defines the system dynamics. For a time τ , the evolution from t to t + τ is given by |Ψ(t + τ ) = U τ |Ψ(t) , or by an expression in terms of density matrix ρ(t + τ ) = U τ ρ(t)U † τ , where U τ = e −iHτ is a unitary matrix. Most generally, the time evolution law is described by a map called a CPTP map L : ρ → ρ with the following properties: linear, trace preserving, Hermiticity preserving, and completely positive. 1

Quantum reservoir computing
We briefly explain the background of QRC. For simplicity, we consider a one-dimensional input and output case. We are given an input sequence u = {u 1 , . . . , u L } and the corresponding target sequencê y = {ŷ 1 , . . . ,ŷ L }, where u k is a continuous variable in [0, 1]. QRC emulates a nonlinear function Q to produce the output y k = Q(w, ρ (0) , {u l } k l=1 ) (k = 1, . . . , L). Here, ρ (0) is the initial state of the quantum system, and w is the parameter that needs to be optimized. A temporal learning task consists of three phases: a washout phase, a training phase, and an evaluation phase. In the washout phase, the system evolves for the first T transient steps to washout the initial conditions from the dynamics. The training phase to optimize w is performed with training data ({u k } L1 k=T , {y k } L1 k=T ), where 1 ≤ T < L 1 < L, such that the mean-square error between y k andŷ k over k = T, . . . , L 1 becomes minimum. The trained parameter w is used to generate outputs in the evaluation phase.
For an N -qubits system, at time t = (k − 1)τ , the input u k ∈ [0, 1] is fed to the system by setting the density matrix of the first spin to u k = (1 − u k ) |0 0| + u k |1 1| ∈ M 2×2 . Therefore, the density matrix ρ ∈ M 2 N ×2 N of the entire system is mapped by a CPTP map where Tr 1 denotes a partial trace with respect to the first qubit. After the input is set, the system continues evolving itself during time interval τ . The dynamics are governed by the Schrödinger equation and the information of the input sequence encoded in the first spin spreads through the system. It follows that the state of the system before the next input u k+1 is where ρ (k) = ρ(kτ ) is the density matrix at t = kτ . It is noted that the size of ρ k grows exponentially with the number of qubits. We obtain partial information regarding ρ (k) by measuring local observables on qubits. The signals for training are obtained from average values of measurement results on each qubit. More precisely, if we employ the ordered basis {O j } in the operator space, then the observed signals at time t are the first N out elements s j (t) = Tr[ρ(t)O j ], where the selection of observables depends on the physical implementation of the system. For example, if we consider the spin receiving the input is an ancilla qubit that we cannot perform the measurement on it, the observed operators act only on the other spins. Since this does not change our main points, for the sake of convenience, we consider a situation in which we can perform measurements on all qubits.
In practical applications, the temporal multiplexing scheme is introduced to improve the performance in extracting dynamics. Here, the signals are measured not only at time kτ but also at each of the subdivided V time intervals during the evolution in interval τ to construct V virtual nodes. The density matrix is then updated by where U (τ /V ) = e −iH(τ /V ) . Therefore, we can obtain temporal signals from N out V nodes, and then the learning procedure is straightforward as we parameterize the linear readout function as Here, x k0 = 1.0 are introduced as constant bias terms, and w = [w 0 w 1 . . . w NoutV ] T represents the readout weight parameters. If we denote K as the number of time steps used in the training phase, w is optimized via the linear regression, or the Ridge regression in the matrix form w = (X X + βI) −1 X ŷ. Here,ŷ = [ŷ 1 . . .ŷ K ] is the target sequence, X = (x ki ) ∈ R K×(NoutV +1) is the training data matrix and β is the parameter serves as the positive constant shifting the diagonals introduced to avoid the problem of the near-singular moment matrix.
We employ the fully connected transverse field Ising model, which is the standard workhorse to build the QR. The Hamiltonian is given by where all the spins interact with each other in x-direction and are coupled to an external magnetic field in z-direction.
Here, σ γ j (γ ∈ {x, y, z}) is the operator measuring the spin j along the γ direction, which can be described as an N -tensor product of 2 × 2-matrices as where 4 Higher-order quantum reservoir computing 4

.1 Model
We propose an effective and practical design to enhance the computational ability of the QR. Our higher-order quantum reservoir (HQR) consists of an ensemble of N qr QRs, such that the lth system Q l has the Hamiltonian H l with N l qubits (which are known as true nodes) and V l virtual nodes. We denote HQR-n as the HQR comprising n QRs. For one-dimensional input, these QR systems are driven by a common input stream, while they can receive different input streams for the multidimensional input setting in general. Here, we consider the setting of one-dimensional input for simplicity. 2 Algorithm 1 presents the temporal processing steps of the HQR. We denote ρ (k) l as the density matrix of Q l at time t = kτ , where ρ (0) l is the initialized density matrix. The reservoir states of Q l at time t = kτ are represented by a V l N l -dimensional vector z kl , which is initialized at t = 0 as the zero vector. At t = (k − 1)τ , the input u kl , which is injected into Q l , is the linear combination between the u k and the linear-scaled reservoir states from other QRs (step 3). The connection coefficients vectors c lm ∈ R (VmNm)×1 from Q m to Q l are randomly generated and fixed, such that elements in c lm are nonnegative and m =l c∈c lm c = α (0 ≤ α ≤ 1), where α is defined as the connection strength parameter. After injecting u kl into Q l , ρ l is transformed by the CPTP map, T u kl , and is then consequently evolved in each τ /V time. The training is performed with the temporal reservoir states of the entire system, which are rewritten with bias terms as the matrix X ∈ R K×( Nqr l=1 N l V l +1) .

Algorithm 1 Temporal processing of higher-order quantum reservoir
Require: The input stream {u k }, N qr QRs Q l with corresponding configurations (H l , N l , V l , ρ (0) l ), the time interval τ , and the nonnegative connection coefficient vectors c lm ∈ R (VmNm)×1 (m = l) such that m =l c∈c lm c = α ≤ 1. z kl ∈ R (V l N l )×1 are denoted as the linear scaled reservoir states of Q l at t = kτ by scaling z → (z + 1)/2 such that their elements are in [0, 1]. 1: for input u k do 2: for l = 1, . . . , N qr do 3: T u kl is a CPTP map defined in Eq. (1) 5: for j = 1, . . . , N l do 8: Measure the average spin values in the z-direction 15: end for 16: return Reservoir temporal states: z = {z 1 , z 2 , . . .}

Properties of higher-order quantum reservoir dynamics
In the classical regime, the reservoir is required to satisfy an asymptotic stability property termed the echo state property (ESP), which ensures that the computations are performed independently with the initial state of the reservoir [24]. ESP is a similar concept with the fading memory property, which states that the reservoir must produce the close outputs if the corresponding inputs are close in recent times [25]. We explore the asymptotic behavior and memory capacity of HQR systems.
Asymptotic stability. A "robust" HQR must produce trajectories that are robust to small perturbations to the system-that is, the computations performed by an HQR system for the same input are independent of its initial density matrix. We define this as quantum echo state property (QESP).
Definition 1 An HQR system whose reservoir dynamics are governed by Algorithm 1 is said to satisfy the quantum echo state property (QESP) when for each initial density matrix ρ (0) ,ρ (0) , and for any input sequence u L = {u l } L l=1 of length L, it holds that ρ (L) −ρ (L) p → 0 as L → ∞. Here, · p denotes the Schatten p-norm for p ≥ 1, defined as From Eq. (1) and Algorithm 1, we can rewrite the evolution from ). Since a CPTP map is a contraction map, the density matrices satisfy decreasing system distinguishability-that is, However, it is not trivial to theoretically explore QESP in higher-order settings [26,27]. Instead, we consider an empirical perspective to the study of the asymptotic stability. We extend the algorithm in [28] to define the QESP index to evaluate the average deviation of observed signals generated from random initial states to reference signals starting from a fixed state [29]. The first T time-steps are discarded as the washout phase, and the QESP index is averaged over P randomly generated initial density matrices. Figure 2(a) presents the average QESP indexes along the variation of J and τ over 10 random trials for the HQR-5 of six qubits, V = 1, and the connection strength α = 0.5 for the range of dynamics in [T, T + 1000] (T = 9000 time steps). The stabilizing effect is empirically confirmed by increasing Jτ , which increases the contractivity of L k . The empirical QESP validity is characterized by a border of a sharp stable-unstable transition. Memory capacity (MC). The property of the HQR in storing information of recent inputs to the current states is commonly measured by memory capacity MC = ∞ d=0 MF(d), where MF(d) is the memory function that evaluates the capacity to reconstruct the previous d steps of the input [30]. This implies that if we set the input as a random sequence {u k } in [0, 1], then the forecast and target are cov(x, y) and σ(x) express the covariance and the standard deviation, respectively. A dynamical system with high MC is superior for temporal forecasting tasks that need to utilize historical patterns. As shown in Fig. 2(b), for each J, there exists an optimal τ to maximize the MC, which can be understood as a trade-off between the influence of the input and the effect of the relaxation dynamics in the QR. Between two consecutive inputs u k and u k+1 , the QR continues evolving for a time interval τ as ξ(τ ) = e −iHτ ξ(0)e iHτ , where we denote ξ(0) as the density matrix after setting u k . The auto correlation function ξ(0)ξ(τ ) reduces to the spin glass order parameter with the exponential decay coefficient λ = ∆J, where ∆ is the smallest nonzero eigengap of H/J [31] (in our Hamiltonian, changing J does not affect the value of ∆). If Jτ is too small, the dynamics approach an identity map, therefore the signals are almost linear without prediction power. If Jτ is too large, the dynamics relax exponentially fast between inputs, thereby decreasing the effects of past inputs into the system.
We further look into the dynamics of the HQR-5 via the bifurcation diagrams where we consider τ is the driving parameter, J = 1.0, V = 1. In Fig. 2(d), the values of four representative signals in one QR (excluding the signals from the input qubit) are overlayed for the range of dynamics in [T, T + 1000]. We consider the diagram at α = 0.5 as an example to explore three intriguing regimes: (A) the dynamics approach fixed points for long enough transient time T , (C) the dynamics exhibits two separated regions in responding to the external inputs, and (B) the transition region between (A) and (C). We note that the left part of (A) appears as a broad band even if the external input is withdrawn (α = 1.0), or setting the constant inputs without the feedback (α = 0.0) [ Fig. 2(e)]. Therefore, this part remains as the effect of the initial state during the short transient time T . As we increase T , the region of τ in (B) will extend to the left. While setting τ in (B) leads to the optimal trade-off for the effect of the past input and the relaxation dynamics, which gives the maximum MC, in practical applications we should set a slightly longer τ to avoid the effect of the initial state.

NARMA tasks
We first demonstrate the performance of the HQRC in the NARMA benchmark, which is commonly used for evaluating the computational capability for temporal processing with long time dependence. The NARMA system is formulated as the nth-order nonlinear dynamical system, which has the following form: y k = κy k−1 + ηy k , where y k andŷ k are the prediction and the target of system at time step k. We present two types of connections in HQR: the mutual connection where any two QRs are mutually connected, and the forward connection where each QR only connects with the next QR in a forward direction. Figure 3(a) illustrates the averaged NMSE over 10 random trials for the HQRs comprising N qr = 1-5 QRs with five qubits, V = 20, and α = 0.0-0.9. We set J and τ in the stable region of the QESP index such as J = 1.0 and τ = 2.0. These values are close to the stable-unstable border in Fig. 2(a) to increase the MC, but at the beginning of (C) as shown in Fig. 2(d) to avoid the effect of the initial state. Since higher-order NARMA tasks require high memory to predict, increasing α will boost the performance for both types of connections. Figure 3(b) illustrates the typical outputs for the HQR-5 in the evaluation phase. The outputs with higher α fit well to the targets for all tasks, even for a difficult task like NARMA20.

Emulating chaotic systems
Next, we demonstrate the performance of the HQRC on emulating chaotic systems, particularly for high-dimensional input that the normal QRC is unable to implement. In this task, the system learns the input of the next step. After learning the readout weights to fit the targets in the stage called teacher forcing stage, the external input streams are withdrawn and the output is fed into the input so that the system is closed. Thereafter, the system can automatically evolve and replicate the dynamical evolution of the target system but within a relatively short prediction horizon. We define  [35], to evaluate the prediction performance, we compute the valid prediction time VPT = Λ −1 1 argmax t f {NRMSE(y t ) ≤ ε, ∀t ≤ t f }, which is the largest time t f (normalized with respect to Λ 1 of the chaotic system) where the NRMSE error is smaller than ε (ε = 0.5 in our experiments). Large VPT means long prediction horizon in the performance of the model. We employ two chaotic systems: the Lorenz attractor and the Kuramoto-Sivashinsky equation (KSE) [32,33] with spatiotemporally chaotic solutions. The Lorenz attractor is given by three ordinary differential equations: dx/dt = a(y − x), dy/dt = x(b − z) − y, dz/dt = xy − cz, where (a, b, c) = (10, 28, 8/3). The time series is obtained by the fourth-order Runge-Kutta method with step size ∆t = 0.01, and the model attempts to predict x(t + ∆t) from x(t). The first K = 10 4 steps are used for training, and 10 3 steps are predicted iteratively. The input signals are min-max scaled to be in the range [0.0, 1.0]. Our HQR-5 of six qubits, with τ = 4.0, J = 2.0, and V ∈ {10, 15, 20}, corresponds with 6 × 5 × V computational nodes. For chaos emulating tasks, the reservoir needs to learn an arbitrarily good approximation of the chaotic recurrence. The non-linear transformation is required, but adding it will degrade the MC [34], thus, we need careful considerations on this memory nonlinear trade-off of the reservoir [29]. In this context, α adjusts this trade-off, which is tuned for the best performance in our experiments. We adopt the framework in [35] to generate the simulation data and to implement the common ML models such as long shortterm memory (LSTM) [36], gated recurrent unit (GRU) [37], and echo state network (ESN) [7]-for example, LSTM-n for the LSTM with n hidden units, and ESN-n for the ESN with n nodes, where n ∈ {80, 100, 120, 150, 500, 1000, 1500, 3000}. Here, we set the number of layers in the LSTM and GRU equal to one. The box plots in Fig. 4(a) demonstrate the distribution of the VPT over 100 random tests for the HQR model, and the two best results in each of the other models. HQR-5 shows superior performance, as its VPTs are highest, while other models cannot capture long time steps in the prediction as shown in the typical examples placed in the right plots of Fig. 4(a). As demonstrated in Fig. 4(b), for different K, the HQR are still better even with a small number of the computational nodes (V = 15 corresponds with 450 nodes). Here, the best results of the ESN, LSTM, and GRU models are selected to plot in Fig. 4(b).
The KSE is the partial differential equation ∂u ∂t + u ∂u ∂x + ∂ 2 u ∂x 2 + ∂ 4 u ∂x 4 = 0 for the scalar function u(x, t) in the interval x ∈ [0, L) with periodic boundary conditions u(x, t) = u(x + L, t). The KSE is integrated on a spatial domain of M = 64 uniform subintervals with ∆t = 0.25, L = 22, thereby yielding a simulation of M -dimensional time series {u k }. We employ the parallel architecture in [38] to build 32 HQRs, where each HQR predicts a spatially 2-dimensional local region. We divide the input u k into 32-local groups u k comprises the ith and (i + 1)th element in u k . The ith HQR is the HQR-10 with six qubits to receive 10-dimensional inputs to predict u (i) k+1 . The input for each HQR is the concatenated vector of u (j) 15, 20, 25}, and tune α for the best performance. The first K = 10 4 steps are used for training, and 400 steps are predicted iteratively. We employ the same parallel architecture for ESN-n, LSTM-n, and GRU-n models, where n ∈ {80, 100, 120, 150, 500, 1000} is the number of nodes or hidden units in each group. Figure 4(c) illustrates contours and error plots for the typical forecasts of HQR and ESN, which demonstrates the ability of the HQR to emulate the spatiotemporal chaos in approximately two Lyapunov time Λ −1 1 . Figure 4(d) illustrates the evolution of NRMSE averaged over 10 random tests for each model. The HQR outperforms other models even with a small number of computational nodes. This may be mainly due to the exponential numbers of degrees of freedoms behind the quantum measurements, which leads to the quantum computational supremacy region.

Conclusion and discussion
In this paper, we propose the HQRC framework, which offers an effective potential means for using quantum dynamics in machine learning tasks. It is suggested that QRC has higher expressive power than classical RC, even for the same number of computational nodes; however, technical scalability is a major drawback. Our framework introduces an implementation to solve this problem. Here, local operations are performed on each QR and the results of those operations are communicated in a classical manner. The experimental results on emulating nonlinear systems, including highdimensional spatiotemporal chaos, demonstrate that our framework can boost the computational power and scalability of QRC through both quantum and classical advantages. Since only time evolution, according to the Hamiltonian as an interaction between nodes, is permitted in QRC, the design for an arbitrary nonlinear function remains as a future challenge. A complete exploration of capabilities of other quantum systems as reservoirs and determining the optimal setting for the classical communications can be a possible direction for the future research.

Broader Impact
As a positive impact in the foreseeable future, our framework can be a typical case for quantum machine learning in the NISQ era for utilizing the noise-robust property of the classical system while exploiting the quantum computational supremacy region. The flexibility of selecting the physical system as a computational resource, the low operation cost when running the hybrid quantum-classical system consisting of multiple small quantum systems, and exhibiting high computational power can enable our research to be utilized in real-world applications. Currently, this theoretical work does not present any foreseeable negative societal consequence. We believe that the experimental efforts in implementing physical quantum reservoir systems, such as the NMR ensemble system [21], the random and noisy quantum circuits on super conducting quantum processors [22], and the photonics system, can have the most benefit from our work. This supplementary material provides a detailed description of the calculations, the experiments introduced in the main text, and the additional figures. The equation, figure, and table numbers in this section are prefixed with S (e.g., Eq. (S1) or Fig. S1, Table S1), while numbers without the prefix (e.g., Eq. (1) or Fig. 1, Table 1) refer to items in the main text.

ALGORITHMS FOR EXPLORING ASYMPTOTIC STABILITY
Algorithm S1 Calculate the quantum echo state property (QESP) index Require: Higher-order quantum reservoir system Q of Nqr QR systems Q l (l = 1, . . . , Nqr), number of washout time steps T , number of evaluation time steps E, number of trials P . Denote F (ρ (0) , {ui} k i=1 ) = [x k1 , . . . , x kN total ] as the reservoir states of the system (N total = l N l V l ) at time t = kτ for input sequence {ui} k i=1 and the initial density matrix ρ (0) . 1: Initialize the density matrices ρ Initialize the density matrices σ end for 9: ∆p = δp(k) k Average the displacement along the trajectory 10: end for 11: ∆ ← ∆p p Average over the trials 12: return The QESP index ∆

MEMORY CAPACITY
The property of a reservoir system in storing information of recent inputs to the current states is commonly measured by memory capacity MC = ∞ d=0 MF(d), where MF(d) is the memory function that evaluates the capacity of system to reconstruct the previous d steps of the input. To calculate MF(d), the system learns to reproduce the previously injected random input of d time steps by using the current states. This implies that if we set the input as a random sequence u k in [0, 1], then the forecast and target are y k andŷ k = u k−d , respectively. MF(d) is then defined as where cov(x, y) and σ(x) express the covariance and the standard deviation, respectively. We note that a dynamical system with high memory capacity is superior for the temporal forecasting tasks that need to utilize historical patterns. In Fig. S1, MF(d) is plotted as a function of d for the HQR system, which contains N qr QR systems. Hereafter, we denote HQR-n for the HQR system consisting of n QR systems. We consider the coupling magnitude J = 1.0 in the Hamiltonian, the number of virtual nodes V = 1 in each QR system, the number of QR systems N qr = 1, 5, and the connection strength α = 0.0, 0.9 with different settings of the time interval τ to inject the input. Here, each QR has five qubits and the number of time steps used in the washout, training, and evaluation phases are 1000, 3000, and 1000, respectively. The memory function is evaluated on 10 samples of HQR with respect to the random coupling. It is evident that increasing N qr will increase the value of MF(d) for any τ , thereby indicating the relatively large value of MF(d) in the region of the large delay d. For the HQR-5 system, the connection strength α = 0 implies that we have a spatial multiplexing setup where all QR systems are disjoint. By introducing the feedback, for example, α = 0.9, the values of MF(d) are higher in the region of the small delay, d. Figure S2 shows MC as functions of τ for the HQR along with different α, N qr [ Fig. S2(a)], J [ Fig. S2(b)], and V [ Fig. S2(c)]. MC is calculated for delay d = 0, 1, . . . , 200, and is averaged over 10 different runs with random trials of Hamiltonian coefficients and connection coefficients in the HQR. We observe from Fig. S2(a) and Fig. S2(c) that increasing N qr and V will lead to an increase in the MC of the system, while the MC becomes saturated around V = 15 [ Fig. S2(c)].
It is interesting to see the optimal values of J, τ in Fig. S2(b) to maximize the memory capacity. These optimal values are close to the border of stable-unstable transition for the QESP index presented in Fig. 2(a) in the main text. The theoretical reasons for this behavior can be understood as a trade-off between the influence of the input and the effect of the relaxation dynamics in the QR. Between two consecutive inputs u k and u k+1 , the QR continues evolving for a time interval τ as ξ(τ ) = e −iHτ ξ(0)e iHτ , where we denote ξ(0) as the density matrix after setting u k , and H is the Hamiltonian of the QR. The auto correlation function G(τ ) = ξ(0)ξ(τ ) at the absolute zero temperature can be written as the ground state (the state of the lowest possible energy) average of the operator ξ(τ ) as the following where denotes the element-wise product, and g f t , g i t , g o t ∈ R d h are the forget, input, and output gates signals, respectively. q t = [h t−1 ,ô t ] ∈ R d h +do is the concatenated vector of the hidden state h t−1 and the observable input The activation functions σ f , σ i , and σ h are sigmoid functions. Similar to ESN, the hidden-to-output mapping is set to o t+1 = W o,h h t , where W o,h ∈ R d h ×do . In our study, we use back-propagation through time (BPTT) algorithm to train LSTM and its variations. We referenced Ref. [2] for a more detailed explanation of BPTT and the implementation of BPTT in the software framework with hyper-parameter settings used in our study.

Gated Recurrent Unit (GRU)
The gated recurrent unit (GRU) [5] is a variation of the LSTM without output gates. GRU uses update gate and reset gate instead to decide what information should be passed to the output; thus, information from long ago can be kept in the training process without excluding information that is irrelevant to the prediction. The recurrent mappings of the GRU are given by

EMULATING CHAOTIC SYSTEMS
In our HQR for the tasks of emulating chaos, we employ the simple form of the Hamiltonian: where g = 2.0, and h i,j is distributed randomly in [−1.0, 1.0].

The Lorenz attractor
The Lorenz attractor is given by three ordinary differential equations: dx/dt = a(y−x), dy/dt = x(b−z)−y, dz/dt = xy − cz, where (a, b, c) = (10, 28, 8/3). The time series is obtained by the fourth-order Runge-Kutta method with step size ∆t = 0.01, and the model attempts to predict x(t + ∆t) from x(t). The first K = 10 4 steps are used for training, and 10 3 steps are predicted iteratively. The input signals are min-max scaled to be in the range [0.0, 1.0]. Our HQR comprises five QR systems (HQR-5) of six qubits, with τ = 4.0, J = 2.0, and V ∈ {10, 15, 20}, corresponds with 6 × 5 × V computational nodes. We adopt the framework in Ref. [2] to generate the simulation data and implement the ESN, LSTM and the GRU-for example, LSTM-n(m) for the LSTM with n hidden units and m layers (if l > 1), and ESN-n for the ESN with n nodes, where n ∈ {80, 100, 120, 150, 500, 1000, 1500, 3000}. We note that in the main manuscript, we consider LSTM and GRU model with one layer to obtain the results in Fig. 4(a)(b). In the preparation of this supplemental material, we became aware that increasing the number of layers in LSTM and GRU leads to an improvement in the performance of these models. We present the results here with the number m of layers is set in {1, 2, 3}. This modification will be addressed in the revised version of our main manuscript. We present in the accompanied software package the full experimental results and the details for the parameter settings.
We evaluate the prediction accuracy using the normalized root mean square error where y t = (y ti ) i andŷ t = (ŷ ti ) i are the forecast and target vectors at time t, andσ i is the standard deviation of the target in time of each component i. To evaluate the prediction performance, we compute the valid prediction time (VPT): which is the largest time t f (normalized with respect to Λ 1 of the chaotic system) where the NRMSE error is smaller than ε (ε = 0.5 in our experiments). Fig. S4 demonstrates the box plots to depict the distribution of VPT over 100 random tests for HQR-5 model with different settings of the connection strength α. We can observe that large α will reduce the average value of VPT. We believe that this may be due to the trade-off of linear and nonlinear dynamics in the reservoir model [7,8]. If we increase the connection strength α in the HQR model, the reservoir becomes more linear and then extensive memory capacity can be obtained. However, for effective functioning of computing in order to address the linearly inseparable problem, a nonlinear transformation of the input signal is required to introduce into reservoir dynamics. Further, by introducing the nonlinearity in general signal-driven dynamical systems, the memory capacity of the reservoir will be degraded (Jaeger conjecture [9]). If the tasks require extensive memory like NARMA tasks, the linear property must be preferred. Therefore increasing connection strength in HQR leads to an improvement in the performance in NARMA tasks. For chaotic emulating tasks with output feedback, the reservoir needs to learn an arbitrarily good approximation of the chaotic recurrence. Therefore, a non-linear transformation is required in these tasks and, thus, we need careful consideration of the memory nonlinear trade-off of the HQR model. In this context, the connection strength α provides a good parameter to adjust this trade-off, which depends on the learning tasks. Figures S5-S7 depict box plots to indicate the distribution of VPT over 100 random tests for HQR, LSTM, GRU, and ESN model. Here, the number of time steps for training is set as K = 10 3 steps (Fig. S5), K = 10 4 steps (Fig. S6), and K = 10 5 steps (Fig. S7). We plot HQR-5 for V=10, 15, 20; and the three best results (with the highest average VPT) in each of the ESN, LSTM, and GRU model. As shown in Fig. S5(a), with only K = 10 3 time steps for training, HQR-5 shows superior performance as its VPTs are highest, while other models cannot capture long time steps in the prediction as shown in the typical examples [ Fig. S5(b)]. As demonstrated in Fig. S8, increasing K leads to improvements for LSTM and GRU models. Nevertheless, the HQR are still better even with a small number of the computational nodes (V = 15 corresponds with 450 nodes). Here, the best results of the ESN, LSTM, and GRU models are selected to plot in Fig. S8.

The Kuramoto-Sivashinsky equation
The Kuramoto-Sivashinsky equation (KSE) is the partial differential equation for the scalar function u(x, t) in the interval x ∈ [0, L) with periodic boundary conditions u(x, t) = u(x + L, t). The KSE is integrated on a spatial domain of M = 64 uniform subintervals with ∆t = 0.25, L = 22, thereby yielding a simulated data set with M -dimensional time series, which can be used as multi-dimensional input-target pairs. The simulation is performed using the fourth-order method for stiff PDEs [10]. We utilize the framework in Ref. [2] to simulate the solution of KSE up to T = 6 × 10 4 , which corresponds to 24 × 10 4 samples. The first 4 × 10 4 are excluded for initial transients, and the remaining data are divided into training and testing dataset of 10 5 samples in each dataset. For the HQR and the ESN model, we employ the augmentation technique proposed in Ref. [11]. Here, the hidden states are augmented such that the hidden states are squared in half of the computational nodes. In order to reproduce the results in the main text, we present the details for the parameter settings in the running scripts of the accompanied software package.