Hybrid quantum-classical reservoir computing of thermal convection flow

We simulate the nonlinear chaotic dynamics of Lorenz-type models for a classical two-dimensional thermal convection flow with 3 and 8 degrees of freedom by a hybrid quantum--classical reservoir computing model. The high-dimensional quantum reservoir dynamics are established by universal quantum gates that rotate and entangle the individual qubits of the tensor product quantum state. A comparison of the quantum reservoir computing model with its classical counterpart shows that the same prediction and reconstruction capabilities of classical reservoirs with thousands of perceptrons can be obtained by a few strongly entangled qubits. We demonstrate that the mean squared error between model output and ground truth in the test phase of the quantum reservoir computing algorithm increases when the reservoir is decomposed into separable subsets of qubits. Furthermore, the quantum reservoir computing model is implemented on a real noisy IBM quantum computer for up to 7 qubits. Our work thus opens the door to model the dynamics of classical complex systems in a high-dimensional phase space effectively with an algorithm that requires a small number of qubits.


I. INTRODUCTION
Quantum computing (QC) and machine learning (ML) have changed our ways to process data fundamentally in the last years [1][2][3][4].Quantum algorithms accelerated the data search [5] or improved the sampling of probability distributions [6,7].These quantum advantages found already their way to various applications [8][9][10], even though we are still in the era of noisy intermediate scale quantum (NISQ) devices that suffer from decoherence and are limited to qubit numbers ∼ 10 2 with resulting shallow quantum circuit depths [11].
Meanwhile, ML algorithms in the form of deep convolutional neural networks extract features effectively and classify big data bases [12][13][14][15].Quantum machine learning ports such methods to a quantum computer [16][17][18] with the prospect that particularly high-dimensional problems can be solved much faster than with their classical counterparts.This expectation arises from two facts, (i) the data space dimension grows exponentially as 2 n with the number of qubits n, the smallest unit of information in QC; (ii) the entanglement of qubits creates highly correlated tensor product states that can represent complex features in the data effectively [19].Thus, for example, quantum support vector machines are expected to have the potential to determine nonlinear decision boundaries of classification problems in high-dimensional quantum enhanced feature Hilbert spaces more efficiently [20][21][22].
Recurrent neural networks (RNN) are specific ML algorithms with internal feedback loops which predict the time evolution of dynamical systems without knowing the underlying nonlinear ordinary or partial differential equations; they can be implemented either as gated RNNs in the form of long short-term memory networks [23] or as reservoir computing models (RCM) [24][25][26][27][28].As a consequence, RNNs have been used for the description of chaotic dynamics, fluid mechanical problems, and even turbulence [29][30][31][32].RCMs were also applied to represent low-dimensional chaotic models, one-dimensional Kuramoto-Sivashinsky equations [33,34], or even turbulent Rayleigh-Bénard convection [35][36][37][38].At the center of the RCM is the reservoir, a randomly initialized and fixed high-dimensionl network of perceptrons which is represented by an adjacency matrix.This specific implementation of an RNN requires only an optimization of the output layer, which maps the reservoir state back to the data space, and avoids costly backpropagation as required in most other ML algorithms [14].
In this work, we combine quantum algorithms with reservoir computing to a gate-based quantum reservoir computing model (QRCM) for a universal quantum computer to predict and reconstruct the dynamics of a thermal convection flow in the weakly nonlinear regime.The algorithm is of hybrid quantum-classical nature since the optimization of the output map is done by a classical ridge regression.The quantum reservoir is composed of a sequence of elementary single and two-qubit quantum gates which form a complex quantum circuit.Following the axioms of quantum mechanics, an elementary quantum gate performs a unitary transformation to a singleor two-qubit state.As a consequence, a highly entangled multi-quibit state will result.
Our first contribution is to demonstrate the feasibility of such a hybrid QRCM to describe the classical chaotic dynamics of a thermal convection flow on an actual NISQ device.The description of the thermal convection flow is based on Lorenz-type Galerkin models with N dof ≤ 8 degrees of freedom [39][40][41][42][43].This class of models is directly derived from the Boussinesq equations of two-dimensional thermal convection between two impermeable parallel plates, heated uniformly from below and cooled from above with free-slip boundary conditions for the velocity field [44,45].Here, we explore QRCMs in two different modes of operation [46]: Secondly, we directly compare the results of the QRCM to its classical counterpart for the same flow.We identify hyperparameters in both approaches that can be related to each other.Note that classical and quantum reservoir computing models differ essentially, which is primarily a consequence of the linearity and unitarity of the quantum dynamics [19].
We demonstrate finally that a systematic reduction of the degree of entanglement in the quantum reservoir by a stepwise transition from a fully to a weakly entan-gled configuration reduces the performance of the present QRCM algorithm.More detailed, this is done by the division of an n-qubit reservoir state into blocks of entangled p-qubit states, so called p-blocks [48].The strong encoding capabilities of fully entangled quantum reservoirs are demonstrated in the present flow case by runs with qubit numbers n < N dof .Also, we show for the open-loop scenario, that the number of operations of the QRCM circuit can be scaled with O(n) < O(2 n ) (where 2 n is the reservoir size).
The research on quantum reservoir computing models proceeds along two major frameworks [49,50].(1) The dynamics of an interacting boson or fermion many-particle quantum system is investigated in the analog framework, which is characterized by a Hamiltonian subject to a unitary time evolution.These systems have been established in the form of spin ensembles [51][52][53][54], circuit quantum electrodynamics [55], arrays of Rydberg atoms [56], or networks of linear quantum optical oscillators [57].In refs.[58,59], the phase transition from a thermalized to a localized many-particle quantum reservoir was studied in respect to the echo state property.The latter describes the ability of the (quantum) reservoir to forget its initial conditions.It is shown that thermalized quantum reservoirs close to the phase transition boundary, for which all spins or oscillators are still strongly entangled, show the best performance for nonlinear learning tasks.A further way to establish a quantum reservoir is by a single nonlinear oscillators [60].
(2) The digital gate-based framework, which sets the stage for the present work, uses circuits composed of universal quantum gates to build a quantum reservoir on NISQ devices [61][62][63].These configurations have been applied for the one-step prediction of nonlinear autoregressive moving-average (NARMA) time series or solutions of the nonlinear Mackey-Glass time-delay differential equation.Here we extend the applications to classical nonlinear dynamical systems with up to 8 degrees of freedom.Furthermore, we apply a reservoir update that blends a linear and a nonlinear activation term, as frequently done in classical reservoir models.
Our work opens the door for the application of quantum machine learning as a reduced-order dynamical model of a higher-dimensional classical complex dynamical nonlinear system.The study thus adds a further proof-of-concept for the potential use of quantum algorithms in studying turbulent flows.
The outline is as follows.In section II, we present the thermal convection flow model; technical details are collected in appendix A. Section III is dedicated to the closed loop scenario.In section IV the complexity of the quantum machine learning task is enhanced to the 8-dimensional model for which we apply an open-loop QRCM.We summarize our work and give a brief outlook in section V. Appendices B, C, and D provide additional material on n-qubit quantum states, the classical reservoir computing framework, and benchmarks of the QRCM with different leaking rates.

II. THERMAL CONVECTION FLOW
We start with a compact description of the flow.The thermal convection flow consists of a two-dimensional fluid layer which is heated uniformly from below with a temperature T bot and cooled from above with T top thus giving ∆T = T bot − T top > 0. The convection flow domain is A = [0, Γ] × [0, 1].The velocity u(x, t) = (u x (x, z, t), u z (x, z, t)) and (total) temperature T (x, z, t) are coupled by the balances of mass, momentum, and energy.The fluid is incompressible and the mass density ρ depends linearly on θ in the buoyancy term only.This is known as the Boussinesq approximation in thermal convection [45].The total temperature is decomposed into T (x, z, t) = 1 − z + θ(x, z, t) where T eq (z) = 1 − z is the static equilibrium profile and θ(x, z, t) is the temperature deviation.The non-dimensional equations are then given by where we use the kinematic pressure p.The Rayleigh number Ra and the Prandtl number σ are the two parameters that characterize the strength of the thermal driving via the temperature difference ∆T and the ratio of momentum to temperature diffusion, respectively.The boundary conditions in x-direction are periodic.At z = 0, 1, one takes They correspond to isothermal, impermeable, free-slip walls.Incompressibility and two-dimensionality allow to reduce the velocity vector field further to a scalar stream function ζ(x, z, t) by This ansatz satisfies (1) automatically and the equations of motion ( 2)-( 3) are now given by with boundary conditions in the vertical direction Equations ( 6) and ( 7) are then reduced by an expansion into trigonometric Fourier modes which satisfy the boundary conditions for the stream function and temperature and encode the spatial structure of the thermal convection flow, see appendix A for further technical details.A subsequent truncation to N and M real time-dependent amplitudes is done for the stream function, {A 1 (τ ), . . ., A N (τ )}, and the temperature, {B 1 (τ ), . . ., B M (τ )}, respectively.This step leads to a class of low-dimensional Lorenztype Galerkin models of the thermal convection flow starting with the original three-dimensional Lorenz 63 model [39] for N = 1 and M = 2 (where N dof = N + M ).The resulting coupled nonlinear system of ordinary differential equations is given by for i, j = 1 . . .N and k, l = 1 . . .M .Here, σ is again the Prandtl number, r the relative Rayleigh number, and b an aspect ratio parameter, see appendix A. Furthermore, F i and G k are quadratic nonlinear functions of the amplitudes A i (τ ) and B k (τ ).We will consider two implementations, the Lorenz 63 model (L63) [39] with N = 1 and M = 2 and an extended 8-dimensional model [42] with N = M = 4 that introduces shear in the flow and causes tilted convection rolls and shearing motion.It thus displays a more complex fluid motion further away from the primary instability point at r = 1 or Ra c = 27π 4 /4 [44].
Figures 2(a) shows two instances of the temperature and velocity fields with the counter-rotating circulation rolls that cause a rise of warm and a descent of cold fluid.These two flow states corresponds to trajectory points of L63 in each of the two butterfly-like wings in Fig. 2(c).

III. CLOSED-LOOP SCENARIO FOR THREE-DIMENSIONAL LORENZ MODEL A. Quantum and classical reservoirs
The design of our time-discrete and gate-based QRCM builds on a n-qubit tensor product state at time t.In appendix B, we provide a compact primer on qubits, tensor product spaces, and entangled or separable states.The n-qubit state in Dirac notation [19] is given by with N res = 2 n .Here |k is the standard basis of the nqubit quantum register.The measured probabilities p k are given by with a t k from eq. ( 11).Here, we have 2 n probabilities; 2 n − 1 of them are linear independent since they have to sum up to 1.The reservoir state evolves from time t to t + ∆t with a fixed time step width ∆t as follows.First, the dynamical part is updated by three blocks of unitary linear transformations with random rotation angles β = (β 1 , ..., β n ), reservoir state probability amplitudes p t = (p t 1 , ..., p t 2 n ), and the past system state vector x t = (x t 1 , ..., x t N +M ), the latter of which summarizes (A 1 , ..., B M ).We simplify the notation in (13) by switching from t + ∆t to t + 1 (or later t+m with m ∈ N).The initial n-qubit state vector |0 ⊗n implies that every qubit is in the basis state |0 .With eq. ( 12) for the probability amplitudes p t+1 k which are obtained from | ψt+1 , the RCM update step outside the quantum reservoir is given by the following iteration The update rule thus contains two terms, a first linear memory term and a second nonlinear activation term.The nonlinearity is connected to the classical data loading as will be discussed in the next subsection.Equation ( 14) contains a leaking rate 0 ≤ ε ≤ 1 that blends both terms.In the classical reservoir computing model, the update of the reservoir state ψ t c would be given by which reveals a similar structure to (14).Here, W r is the reservoir matrix and W in the input matrix.More details are provided in appendix C. A leaking rate of ε = 1 implies that only a nonlinear activation by the hyperbolic tangent is present -a mode in which several, but not all classical reservoir computing models are operated [24,65].Inubushi and Yoshimura [66] term this split the memory-nonlinearity trade-off since the nonlinear activation term typically will degrade the memory of the system.Differently to the analog framework of quantum reservoir computing [51][52][53][54][55][56][57][58][59] that processes the state without external memory in the reservoir, we add external memory by the first term in (14).An improved performance of the present hybrid quantum-classical reservoir computing model for ε < 1 in comparison to one with ε = 1 is demonstrated in appendix D for two common benchmark cases.

B. Classical data loading and reservoir state evolution
We now specify the loading procedure of the classical data into the quantum reservoir.The unitary transformations of eq. ( 13) consist of single-qubit rotation gates R Y and subsequent two-qubit controlled NOT (in short CNOT) gates.The R Y -gate is defined by shows the corresponding circuit diagram of the quantum reservoir which consists of three circuit blocks as lined out in eq. ( 13).The first block of unitary transformations U (4πp t ) loads the reservoir state probability amplitudes of the previous time step t (indicated as the blue box).This is done by rotation gates R Y (4πp t k ).In the circuit diagram of Fig. 2, these 2 9 = 512 probabilities with 0 ≤ p t k ≤ 1 are summarized to a vector to keep the notation less crowded.For this as for all the following blocks, the combination of R Y and CNOT gates is continued until the last qubit is reached.There, the CNOT is applied to the previous qubit and if not yet finished, the constructor starts at the upper qubit again.
The application of an R Y rotation gate, which is parametrized by the input value, is a nonlinear operation in terms of the amplitudes [62,67].It can be considered as an analogy to the nonlinear activation in a classical RCM, e.g., by tanh(•), see also table I where we summarize hyperparameters of the classical and quantum RCMs.Note also that all 2 9 amplitudes p t k are loaded into the reservoir in this case.
Similarly, the degrees of freedom x t i at time t are loaded into the quantum reservoir in the second block U (4πx t ) before the model advances further (indicated as the yellow box) to the last block (indicated as the gray box).This third block U (β) performs additional rotations by angles β i .It stands for a unitary evolution step of the full reservoir state which enhances the entanglement and randomization.The rotation angles β i are sampled initially in a reproducible way from a uniform distribution between 0 and 2π, which corresponds to a random seed initialization of a classical reservoir.
The loading of the full reservoir state, which becomes exponentially more costly, was necessary to obtain the reported prediction horizons for the closed-loop scenario.This is discussed in the next subsection.

C. Quantum reservoir readout and classical optimization
A projection-valued measure in the standard basis of the Pauli-Z operator provides the probabilities p t k from K 2 n independent circuit simulations, known as shots.These probabilities are mapped to the updated dynamical system state by the output matrix, with the optimized weights which are summarized in the matrix W out * ∈ R (N +M )×Nres .We note once more that the output matrix is optimized by a classical algorithm similar to the classical RCM case.This optimization seeks a minimum of the cost function C(W out ) which is given in appendix C. Panel (d) of Fig. 2 and table I compare the classical and quantum RCM with the numerical simulation of the equations of motion obtained by a 4th-order Runge-Kutta method.The integration time is rescaled by the largest Lyapunov exponent of the system, λ 1 = 0.9056, which quantifies the deterministic chaos of the model [64].
The training phase comprises N train = 2000 time steps, both for the classical and quantum case.The first 50 time steps out of N train are used for the washout of the initial reservoir state.For times t ≥ 0 the reservoirs are exposed to unseen test data predicting the dynamics autonomously.It is seen that the prediction horizon of the QRCM with 9 qubits is about 1.5λ 1 t in this example.This result remains nearly the same for different reservoir seeds.Either the approximation of the L63 model by reservoir dynamics or the additional white noise in the Qiskit simulator will cause a switch of the trajectory into the other wing of the butterfly-like Lorenz attractor.Note that the classical RCM (CRCM) prediction will also deviate from the ground truth at about 3λ 1 t which is not shown here.The leaking rates in this example are given in table I.
Two hyperparameters are varied, the leaking rate ε and the number of qubits n that determines the reservoir dimension N res .We identify a minimum of the cost function in the form of a mean squared error (MSE) around ε = 0.025.This is the statistical minimum while the single best representation shown in Fig. 2(d) has ε = 0.05, as stated in table I.The larger the number of qubits the smaller MSE, although the improvements in the Qiskit simulations remain small (and thus the difference of the displayed to the optimal case).Note also that each data point for the MSE in Fig. 2 memory-dominated blended with a small nonlinear contribution [66], as we detailed already in Sec.III A.
The Tikhonov regularization parameter γ which is added to the cost function to avoid overfitting was set to γ = 0 in the present as well as in the NISQ device runs.The noise in the ideal quantum simulator and the decoherence in the NISQ device have a regularizing effect, see discussions of this aspect in ref. [65] for classical and in [63] for quantum devices.We will come back to this hyperparameter in the next section.All details on the classical reservoir computing model, the hyperparameters, and the cost function are in appendix C in order to keep the work self-contained.
The quantum reservoir readout requires K projective measurements of the n-qubit reservoir state on identically prepared quantum systems which might be costly in comparison to the CRCM case.However, a readout of the data has to be done in the classical case as well.Throughout this study, we used K = 2 10+n shots to reduce the measurement noise.A further increase of the number of shots did not reduce the MSE.Since the NISQ devices allow a maximum number of K = 8192 shots only, batch jobs with a pre-defined Qiskit function (combine counts) were used.
In the closed-loop scenario, we conduct these measurements after each step to monitor the time evolution at equidistant instants from t to t m = t+m (m > 1, m ∈ N).This is the same procedure as in the classical algorithm such that both can be compared to each other.With an increasing number of qubits the state vector grows exponentially.Larger qubit numbers can also require a number of shots K for the quantum simulator or quantum device that goes beyond our presently suggested one.In the open-loop scenario, which will be discussed in the next section IV, we advance the system from t to t + 1 only to obtain the results.This step always closes with a readout in the form of a measurement.

IV. OPEN-LOOP SCENARIO FOR 8-DIMENSIONAL LORENZ-TYPE MODEL
A. Quantum reservoir with one continually available degree of freedom We proceed from the standard L63 model to an extension with 8 degrees of freedom, which is listed and explained further in appendix A. As shown by Gluhovsky et al. [42], this extension can be decomposed into subgroups, so-called gyrostats.The model conserves total energy and vorticity.They are given by with a kinetic and potential energy term and with the vorticity ω = −∇ 2 ζ and the convection domain size A. The open-loop scenario of the QRCM implies that a subset of the degrees of freedom will remain continually available in the reconstruction phase after the training phase.In this subsection, we will take N in = 1 which will be A 4 (t).The leading Lyapunov exponent was computed numerically by a method proposed in [68] and turned out to be λ 1 = 0.825.Figure 3 displays the results for the 8-dimensional Lorenz-type model which receives the time series A 4 (t) to reconstruct the remaining 7 degrees of freedom of the thermal convection flow model.Panel (a) compares the times series of the ground truth (GT) with the results of a CRCM, and a QRCM which was run on n = 7 qubits on an ideal Qiskit simulator.We see that the data remain closely together for the displayed interval of more than 16 Lyapunov time units.The QRCM runs again through a training phase of N train = 2000 integration time steps with a leaking rate ε = 0.05 after an initial washout of 50 time steps.Figure 3 the lower-dimensional Lorenz 63 model and lead to tilted convection cells, as can be seen in the panels.
The dimension of the quantum reservoir is N res = 128, while the one of the CRCM in Fig. 3 is N res = 512.We now compare the mean squared error (MSE) as a function of the dimension of the CRCM N res and the regularization parameter γ in Fig. 4. We also show the behavior of the QRCM N res = 128 for comparison.The MSE is given as see also eq. ( 16) and T test = 2000.Subscript tg stands for target and denotes the ground truth (GT) which is obtained by time integration of the model eqns.( 9) and (10).The MSE of the CRCM is large for all reservoir dimensions when γ is very small, but improves as γ rises.
The minimal MSE for the CRCM is ≈ 2.8×10 −3 at γ = 1 and N res = 512 while the minimal value of the QRCM is ≈ 1.36 × 10 −3 .The MSE of the QRCM remains practically unchanged for γ 10 −3 and starts to grow then moderately.As we discussed in the last section already, the noise of the quantum algorithms seems to provide already enough regularization.We can conclude that the MSE of the QRCM is found to be fairly close to that of the CRCM with a reservoir dimension that is reduced by a factor of 4.

B. Implementation on an actual quantum device
The 8-dimensional convection flow model is finally implemented on an actual noisy quantum device.Figure 5(a) displays the quantum reservoir for the implementation.The circuit depth on real devices is still rather lim- ited by the decoherence of the elementary quantum gates that are installed in the form of microwave-controlled superconducting SQUIDs.The figure shows that we had to reduce the original three-block-structure of the quantum reservoir, which was used for the results that are shown in Figs. 2 and 3, to one block.Two further steps were necessary: (1) Instead of one continually available variable in the reconstruction phase, we provide now A 4 (t) and B 3 (t).The reason is that the shallower circuit was found to be too noisy for the reconstruction of 7 degrees of freedom from one continually available degree of freedom.(2) Only 14 out of the 128 components of the reservoir state measurement vector p t are fed back into reservoir together with the two degrees of freedom.This reduces the cost of loading data into a quantum register significantly.The total qubit number was limited to n = 7.The studies were conducted on two devices, ibmq ehningen, a 27-qubit quantum computer in Germany, and ibm perth, a 7-qubit machine.Figure 5(c) displays the arrangement of the 7 qubits on ibm perth which was mainly used for the results in Fig. 5(b).The quantum computer ibmq ehningen was used for reservoir computing runs with ε = 1, which gave a reduced network performance.The additional option to take the best callibrated 7 qubits out of 27 did not lead to significant improvements.Entanglement operations, e.g. by C-NOT gates, are only possible for qubits which are connected by the bars in Fig. 5(c).No error correction was performed.
We backed up this investigation by two runs on the Qiskit simulator with the same configuration.One is the ideal simulator that has been used already before.The other simulation was done on a noisy Qiskit simulator for which you can prescribe the probabilities of measurement errors, here p m = 0.05, gate errors, here p = 0.1, and qubit resets, here p r = 0.03.Values have been chosen such that they come close to those on real devices.All environments are compared in Fig. 5(b).The data from the noisy Qiskit simulator and the real quantum device do partly deviate, but are found to follow the overall trend fairly well.This proves the concept of a hybrid QRCM for a classical dynamical system on a NISQ device.

C. Stepwise reduction of reservoir entanglement and quantum advantage
Finally, we investigated if a simulation of the Lorenztype model with the simpler quantum reservoir than from Fig. 5(a) is successful when the corresponding quantum circuit is decomposed into several p-qubit-blocks which are disentangled.If p = n the circuit is fully entangled, for p = 1 the n-qubit quantum state is fully separable; see appendix B for the definitions of both possible multiqubit quantum states.The decomposition is illustrated in Figs.6(b,c).The rational behind this analysis, which we did with the ideal Qiskit simulator for n = 8, is that a p-blocked structure might be simulated efficiently on a classical computer loosing the quantum advantage [48].
In Fig. 6(a), we summarize the MSE in a diagram for circuits with 3 ≤ n ≤ 8 and possible block size 2 ≤ p ≤ 8.For example, n = 4 and p = 3 imply that a single qubit remains which is disentangled from the 3-qubitblock.In general, the number of blocks of size p follows by n p = n/p .Each possible p-blocked reservoir configurations at n qubits was trained and then run for 100 different realizations to gather statistics.Since the reduced structure of Fig. 5(a) is used, there is no U (β) block.To generate the 100 different realizations, we took 14 values p i randomly out of the 128 possible probabilities as input.The block diagram shows that the MSE decreases when the block size is increased.We can conclude from this analysis that the entanglement of the qubits in the reservoir is essential for the performance of the QRCM.This is different compared to the classical reservoir which is a sparse network for which 20% of the network nodes are actively only, but the number of perceptrons is by 2 to 3 order of magnitude larger in comparison to the number of qubits of the quantum reservoir.
Finally, we estimate the number of operations for the CRCM and our QRCM.While most tasks have equal computational effort, our comparison is related to the second term in ( 14) which provides the nonlinear activation.In the classical case, this requires mainly O((N c res ) 2 ) operations for the first term in the tanh(•) activation; see ref. [36] for an analysis of the occupation of the reservoir matrix.Superscripts c and q of N res stand for classical and quantum reservoir dimension, respectively.The second term with N c res × N in operations is subdominant when N in N c res .In the quantum case, we have O(ξn) gate operations for a shallow circuit as the one given in Fig. 5(a) that spans a reservoir of dimension N q res = 2 n .Here, the prefactor of ξ 3 is the approximate amount of operations for a single qubit.This has to be multi- plied with the number of shots, which was set here to K = 2 10+n .We compare the amount of operations for N = N c res = N q res with N being sufficiently large.Thus follows the inequality or, in terms of states N = 2 n , for having less reservoir operations of the quantum case.
In order to show a quantum advantage for this framework rigorously, one needs to prove that the formula K = 2 10+n is still appropriate and that ξ can still be choosen constant for increasing qubit number.Since the QRCM requires typically less nodes than the CRCM, i.e., N q res N c res as seen in Fig. 4, we expect that the QRCM might be able to outperform its classical counterpart, at least for the class of problems discussed here.It is clear that future investigations have to show if this is the case.

V. SUMMARY AND OUTLOOK
The main objective of our present work was to show the feasibility of a hybrid quantum-classical reservoir computing model to predict and to reproduce the dynamical evolution of a classical, nonlinear thermal convection flow, on an actual quantum computer with up to 7 superconducting qubits.In a nutshell, quantum reservoir computing models are recurrent machine learning algorithms for which the reservoir state is built by a highly entangled tensor product quantum state that grows exponentially in dimension with the number of qubits.
Our work showed that a quantum reservoir has a qubit number that is by about 2 orders of magnitude smaller than that of the perceptrons in a classical one.On the one hand, we could thus take advantage of the data compression capabilities of quantum machine learning algorithms where the dimension of the data space grows exponentially with the number of qubits which is essential for the modeling of higher-dimensional nonlinear dynamical systems.On the other hand, it shows that a classical reservoir state which is caused by a sparsely occupied network matrix of dimension 10 3 can be substituted by a highly entangled quantum state that is caused by the application of unitary transformations.The qubit number was n < 10 in the present case.
The study can be extended into several directions.It is clear that the present thermal convection flow model is still very low-dimensional and thus far away from convective turbulence.Our efforts should be considered as one first step to model real fluid flows on a quantum computer.It provides a possible route beside other directions, such as quantum embeddings of nonlinear dynamical systems by the Koopman operator framework [69] or variational quantum algorithms for the direct solution of the equations of motion [70], see also ref. [71] for further directions such as lattice Boltzmann methods.Extensions to higher-dimensional models are currently still limited by the technological capabilities of quantum computers.As the technological progress in this field is very fast, it can be expected that Galerkin models with significantly more modes will be modeled on upcoming devices with chips with a higher noise-resilience and lower error rates at the gates.The model that we applied here can be systematically extended towards turbulent convection, as discussed in detail in refs.[43,72].A QRCM with n ∼ 10 might thus be able to run a two-dimensional turbulent convection flow usable as a subgrid-scale superparametrization in a global circulation model [73].
In the present work, we have not systematically optimized the circuit architectures for the different tasks.Further reductions of the number of gates caused always reduced prediction and reconstruction capabilities in the closed-and open-loop scenarios, respectively.A possible route of research would thus be to compose the different quantum reservoirs more systematically from first principles, e.g. in the form of a multilayer tensorial network that potentially improves the performance of quantum algorithms on NISQ devices, see e.g.ref. [74].
sorting the terms with respect to the wavenumbers, and truncating higher wavenumbers leads to closed systems of coupled nonlinear ordinary differential equations.In the present work, we consider Lorenz-type models up to order 8 (N = M = 4); higher-dimensional models have been investigated for example in refs.[42,43].In detail, we take the expansions In this appendix, we provide details about the reservoir computing approach, a recurrent supervised machine learning algorithm, which here is implemented in the form of an echo state network with a leaking rate ε.The training of the RCM proceeds as follows.The dynamical system state at time t, which is denoted more compactly as is mapped to the reservoir state ψ t via the randomly initialized input weight matrix W in ∈ R Nres×Nin .Here, N res N in is the reservoir dimension.The reservoir state ψ t is updated as follows [26,28,35], see eq. ( 15) in the main text, This update rule comprises external forcing by the inputs x t as well as a self-interaction with the reservoir state ψ t .The two terms on the right hand side of (C2) are combined by the leaking rate ε.The hyperbolic tangent tanh (•) is the nonlinear activation function of each reservoir node.The randomly initialized matrix W r represents the reservoir, a sparse random network of neurons [65].Thus the leaking rate ε ∈ (0, 1] moderates the linear and nonlinear contributions.The updated reservoir state ψ t+1 is mapped via the output matrix W out ∈ R Nin×Nres to form the reservoir output x t+1 ∈ R Nin The elements of W out have to be computed.Therefore, a set of T training data instances {x t+1 , x t+1 tg }, where t = −T, −T + 1, ..., −1, needs to be prepared.The target output x t+1 tg (also denoted as ground truth (GT)) represents the desired output that the RCM should produce for the given input x t .The resulting data pairs are assembled into a mean squared cost function C (W out ) with a Tikhonov regularization term which is given by and has to be minimized corresponding to W out * = arg min C(W out ). Superscript T denotes the transposed.The regularization parameter γ > 0 avoids overfitting [14].The optimized output matrix is given by where I is the identity matrix.U tg and R are matrices where the t-th column is the target output x t tg and reservoir state ψ t , respectively.The optimization of the output weights thus becomes computationally inexpensive.The hyperparameters of the classical RCM are N res , ε, γ, the reservoir density D, and the spectral radius ρ(W r ).
Once the output weights are optimized and the hyperparameters are tuned the RCM can run either in the prediction (closed-loop scenario) or reconstruction mode (open-loop scenario).Equation (C2) changes in the closed-loop scenario to ψ t+1 =(1 − ε)ψ t + ε tanh W in W out * ψ t + W r ψ t . (C5) Now the RCM can work independently of training input.The prediction for the dynamical system follows by x t+1 = W out * ψ t+1 .Equation (C2) remains the same in the open-loop scenario, except that the continually available input vector is very low-dimensional in this regime, see Fig. 1.The full state reconstruction follows again by x t+1 = W out * ψ t+1 .The latter case is also called onestep prediction since x t+1 is not used as a new input for x t+2 , differently to the former prediction mode.Here, α = 1, β = 2, γ = 1, and T = 2.The time τ is measured here in multiples of the time step width ∆τ = 0.1, i.e., τ = k∆τ with k ∈ N. Figure 7 compares the aforementioned benchmarks for two leaking rates which were processed with our QRCM, either with ε = 1 or ε = 0.2.All runs were done in Qiskit.The first 100 time steps are used for washout, the subsequent 400 time steps for training.We clearly observe a significant improvement of the performance of the hybrid quantum-classical reservoir computing model with a leaking rate of ε < 1.The QRCM prediction with ε = 0.2 follows the ground truth nearly perfectly.

FIG. 1 : 2 .
FIG. 1: Sketch of the two scenarios in which the reservoir computing model is run.(a) Closed-loop scenario for autonomous prediction of the dynamics.(b) Open-loop scenario for reconstruction of dynamics from continually available data.The matrix W out * stands for the classically optimized output layer.N dof is the number of degrees of freedom of the dynamical system, Nin is the number of continually available components of the system state vector.The dimensionality of the reservoir is Nres N dof .For (a), Nin = N dof and for (b), Nin N dof .

FIG. 2 :
FIG. 2: Quantum reservoir computing model (QRCM) for the Lorenz 63 dynamical system.(a) Two instantaneous convection flow states which display the velocity vector field (ux, uz) together with the total temperature field T as a colored background (blue for Ttop and red for T bot ).(b) Circuit diagram for the QRCM.Here (x1, x2, x3) = (A1, B1, B2).The three groups of unitary operations are indicated by differently colored boxes.The scaling of input parameters of the RY -rotation gates is specified in the text.(c) Trajectory plot of the Lorenz 63 system in the phase space which is spanned by one stream function mode A1(t) and two temperature modes, B1(t) and B2(t), i.e., N = 1 and M = 2.(d) Comparison of classical and quantum reservoir computing in the prediction phase (yellow background).Time is rescaled by the largest Lyapunov exponent λ1 = 0.9056, see e.g.ref. [64].The two flow configurations in (a) are also indicated in panels (c,d).(e) Mean squared error (MSE) as a function of the leaking rate ε and the number of qubits n.The model parameter are σ = 10, b = 8/3, and r = 28.All displayed QRCM runs are done with the Qiskit simulator.
FIG. 3: Comparison of the 8-dimensional Lorenz-type model, time-integrated system of equations as ground truth (GT), classical reservoir computing model (CRCM), and quantum reservoir computing model (QRCM).Panel (a) displays the time evolution of all variables.(b) Reconstruction of the flow and temperature fields at times t1 to t4 which are indicated in A1(t).The model parameters are σ = 10, b = 8/3, and r = 28.Here, N = 4 and M = 4. Mode A4 is the only input and always given accurately into each reservoir.

FIG. 4 :
FIG. 4: Comparison of classical and quantum reservoir computing models in the reconstruction phase for different regularization parameters γ which have been increased by factors of 10 from 10 −10 to 10.Each data point corresponds to the median of 30 different random reservoir realizations.The remaining hyperparameters remain fixed throughout this study.

FIG. 5 :
FIG. 5: Quantum reservoir computing model run of the 8-dimensional time-integrated Lorenz-type model an actual quantum device.(a) Sketch of the quantum reservoir which had to be reduced due to decoherence in comparison with the one that is displayed in Fig. 2. (b) Time series of the extended Lorenz model.We compare the ground truth (Model) with an ideal and noisy Qiskit simulator and the 7-qubit quantum computer (IBM Q).The number of training steps was again Ntrain = 2000 and the leaking rate ε = 0.2.The two degrees of freedom that are continually available in the reconstruction phase are indicated.(c) Connection of the 7 qubits on the ibm perth quantum computer.

FIG. 6 :
FIG. 6: Performance of the quantum reservoir computing model for different reservoir architectures.(a) Mean squared error on a logarithmic scale as a function of the total number of qubits and the size of the blocks of entangled qubits.The dark cells in the lower left stand for impossible decompositions.(b) Sketch of an example case.Fully entangled 4-qubit-quantum circuit which is the normal setting.(c) Two fully entangled 2-qubit blocks (p = 2) build the 4-qubit-quantum circuit.
Appendix D: NARMA-2 model and Mackey-Glass equation for different leaking rates In this appendix, we demonstrate the necessity of ε < 1 for two common reservoir computing benchmark cases.The first case is the NARMA model, an input-output model class with input u k with k ∈ N given by u k = 0.1[sin(2παk) sin(2πβk) sin(2πγk) + 1] .(D1) Here, α = 2.11/T , β = 3.73/T , γ = 4.11/T , and T = 100.The output y k is then given by the following iteration rule y k+1 = 0.4y k + 0.4y k y k−1 + 0.6u 3 k + 0.1 .(D2) We use y 0 = y 1 = 0.19.The recursive character of the discrete time series can be enhanced by adding further terms from the past.Here, we take a NARMA-2 model since y k+1 depends on y k and y k−1 .The second case is the Mackey-Glass equation, a nonlinear time-delay differential equation, which is given by dx(τ ) dτ = βα n x(τ − T ) α n + x(τ − T ) n − γx(τ ) .(D3)

TABLE I :
(e) is obtained as an average over 50 different random seeds of the quantum reservoir, i.e., 50 different random vectors with angles β i .A small leaking rate implies that the reservoir dynamics is Comparison of classical and quantum reservoir computing models.Different essential quantities including optimal hyperparameters for the Lorenz 63 model in the closed-loop scenario are listed.The spectral radius ρ(W r ) in the quantum case is always equal to 1 since unitary transformations are norm-preserving.The number of qubits is n.Two additional hyperparameters are used in the classical RCM: a reservoir density D = 0.2 which determines the percentage of active nodes in the matrix W r and an additional Tikhonov regularization term with a parameter γ = 10 −1 in the cost function C(W out ), see appendix C.