Reduced-order modeling of two-dimensional turbulent Rayleigh-B\'enard flow by hybrid quantum-classical reservoir computing

Two hybrid quantum-classical reservoir computing models are presented to reproduce low-order statistical properties of a two-dimensional turbulent Rayleigh-B\'enard convection flow at a Rayleigh number Ra=1e+5 and a Prandtl number Pr=10. These properties comprise the mean vertical profiles of the root mean square velocity and temperature and the turbulent convective heat flux. Both quantum algorithms differ by the arrangement of the circuit layers of the quantum reservoir, in particular the entanglement layers. The second of the two quantum circuit architectures, denoted as H2, enables a complete execution of the reservoir update inside the quantum circuit without the usage of external memory. Their performance is compared with that of a classical reservoir computing model. Therefore, all three models have to learn the nonlinear and chaotic dynamics of the turbulent flow at hand in a lower-dimensional latent data space which is spanned by the time-dependent expansion coefficients of the 16 most energetic Proper Orthogonal Decomposition (POD) modes. These training data are generated by a POD snapshot analysis from direct numerical simulations of the original turbulent flow. All reservoir computing models are operated in the reconstruction mode. We analyse different measures of the reconstruction error in dependence on the hyperparameters which are specific for the quantum cases or shared with the classical counterpart, such as the reservoir size and the leaking rate. We show that both quantum algorithms are able to reconstruct the essential statistical properties of the turbulent convection flow successfully with similar performance compared to the classical reservoir network. Most importantly, the quantum reservoirs are by a factor of 4 to 8 smaller in comparison to the classical case.

In classical reservoir computing, the reservoir is the central building block of the neural network architecture.The reservoir is a sparse random network of neurons that substitutes the batch of successively connected hidden layers of deep convolutional neural networks of other machine learning algorithms [24,25].It introduces a short-term memory to process sequential data.This is the subject of the present investigation.Here, we will substitute the high-dimensional classical reservoir network by a small parametric quantum circuit in which n qubits span a 2 n -dimensional complex quantum state space for a highly entangled reservoir state to save memory and computational costs.Quantum reservoir computing can be implemented in two different ways which we describe in brief in the following.
The dynamics of an interacting many-particle quantum system -the quantum reservoir-is investigated in a so-called analogue framework.It is characterized by a Hamiltonian H subject to an unitary time evolution promoted by U (t).The time evolution of the density matrix ρ(t), which describes the quantum reservoir state, follows then to The operator H is the (many-particle) Hamiltonian and U † (t) the adjoint of U (t).These systems have been implemented in the form of spin ensembles [26][27][28][29], circuit quantum electrodynamics [30], arrays of Rydberg atoms [31], single oscillators, and networks of oscillators [32,33].They establish a closed quantum system in the ideal case that follows an ideal unitary time evolution after the input state is prepared.The pure state density matrix ρ(t) is given by the outer product of the (many-particle) state vector |Ψ(t)⟩ with itself ρ(t) = |Ψ(t)⟩⟨Ψ(t)|.
Beside the analogue framework , the digital gate-based framework uses parametric circuits composed of universal quantum gates.They are composed to a quantum reservoir on noisy intermediate-scale quantum devices in this case [34,35].The reservoir state is obtained by a repeated measurement of the equivalently prepared quantum system and gives the probabilities p j (t) for j = 1, . . ., 2 n of the n-qubit quantum state |Ψ(t)⟩ to collapse on the j-th eigenvector of the standard observable in a quantum computer, the Pauli-Z matrix [36].These probabilities correspond to the diagonal elements of the density matrix ρ and are summarized to the vector p(t) = (p 1 , . . ., p N ) T = (ρ 11 , . . ., ρ N N ) T with N = 2 n .They can be red out by a measurement.Again, we assumed that the initial state ρ(0) is a pure state.In ref. [37], the digital gate-based approach was realized in the form of an open quantum system which implies that parts of the short-term memory of the hybrid systems are kept outside the quantum reservoir.This method allowed to model the dynamics of low-dimensional nonlinear systems, such as the Lorenz 63 [38] and its 8-dimensional Lorenz-type model extension [39] on an actual IBM quantum computer.This algorithm will be denoted to as hybrid algorithm 1, in short H1, in the following.
In the present work, we advance our investigation on quantum reservoir computing with a proof-of-concept application to a realistic complex fluid mechanical problem.We seek to show that it is possible to achieve results for the following low-dimensional simulations which are comparable to classical reservoir computing.To this end, we present a full data processing pipeline for a twodimensional turbulent Rayleigh-Bénard convection flow [40] which contains a quantum computing module -the quantum reservoir.This flow is a paradigm for turbulence that is driven by buoyancy forces in many geophysical and astrophysical processes [41,42].The hybrid quantum-classical machine learning algorithm will serve as a data-driven reduced-order model of the turbulent flow without knowledge of the underlying mathematical equations of motion.
The hybrid nature of the quantum machine learning model includes a reduction of the high-dimensional turbulence data to a low-dimensional latent space.This is done by a classical snapshot-based Proper Orthogonal Decomposition (POD) [43].Similar to classical machine learning algorithms, this encoding/decoding step is necessary for a fully turbulent flow since the dimension of the classical input data is high; the actual number of degrees of freedom is here Ñdof = 3×384×96 = 110592, see also [44,45].The quantum machine learning algorithm thus operates in a latent data space of N dof = 16 in the present case and is able to reproduce relevant large-scale features and low-order statistics of the turbulent flow, such as the vertical profile of the mean convective heat flux across the convection layer.Particularly, the latter point is of particular interest in the present application.
We also extend our previous study with an improved hybrid quantum-classical reservoir computing model (RCM) which integrates more parts of the algorithm into the quantum computing part in comparison to the previous algorithm H1 from ref. [37].The present work is a first step away from the traditional von Neumann architecture, in which computation and memory are located in distinct components.The new algorithm will be denoted to as H2 in the following.It will be compared in detail with a classical RCM, in short C, and H1, our previous approach.
The hybrid nature of our algorithm implies additionally that the optimization of the reservoir output layer is performed classically by a direct solution of the minimization task.The full data processing pipeline of the algorithm comprising a combined POD-RCM model is sketched in Fig. 1.The figure sketches the classical and quantum reservoir in panels (b) and (c).The quantum reservoir builds on a low-qubit-number parametric quantum circuit which spans a high-dimensional reservoir state space based on a highly entangled n-qubit quantum state.
The paper is organized as follows.In Sec.II, we present the turbulent flow and in Sec.III the reduction to the low-dimensional latent space in which the quantum reservoir operates.Sec.IV follows with a compact presentation of the algorithms C, H1, and H2.Sec.V discusses the results in dependence on hyperparameters of all three reservoir computing models.Moreover, we also compare different error measures.One is adapted to the specific fluid mechanical application.A summary and an outlook are given in the last section.

A. Model equations and parameters
We consider a two-dimensional Rayleigh-Bénard system where a fluid is enclosed between two impermeable plates with constant temperature difference ∆T = T bot − T top > 0 [40].The Boussinesq equations connect the incompressible velocity field u = (u x , u z ) with the temperature and pressure fields, T and P.They are given by Here, α, g, ν, and κ are the thermal expansion coefficient, the acceleration due to gravity, the kinematic viscosity, and the thermal diffusivity, respectively.We set P = P/ρ 0 in eq. ( 3).Furthermore, T 0 and ρ 0 are constant reference values of the temperature and mass density, respectively.These equations stand for the differential balances of mass density (2), momentum density (3), and energy density (4) of a fluid parcel.In the Boussinesq approximation, it is assumed that the fluid is incompressible (or divergence-free) and that small deviations of the density of the fluid from the reference values depend linearly on temperature deviations [40].This leads to the last term in the momentum equation (3) which couples the temperature field with the vertical velocity component.
The system is made dimensionless by the choice of the free-fall velocity scale U f = √ αg∆T H, the free-fall time scale T f = H/U f , and ∆T .Here, H is the height of the convection layer, the characteristic spatial scale in the setting.In this way all material parameters and scales can be summarized in dimensionless parameters that determine the operating point of the turbulent flow.These parameters are the Rayleigh and the Prandtl numbers, Ra = αg∆T H 3 /(νκ) and Pr = ν/κ.They take values of Ra = 10 5 and Pr = 10 in the present proof-of-concept study.Alternatively, one can use U diff = κ/H as a characteristic velocity [37,38].This does not affect the physical outcome.Figure 2 illustrates the configuration which we want to investigate in the following by the hybrid quantum-classical algorithm.
We conduct direct numerical simulations using the spectral element solver Nek5000 [46] to solve the Rayleigh-Bénard system (2)-(4) in a domain A = L × H with aspect ratio Γ = L/H = 2 √ 2. Dimensionless co-ordinates are thus x ∈ [0, 2 √ 2] and z ∈ [0, 1].Dirichlet boundary conditions are imposed for the temperature field at top and bottom, T (z = 0) = 1 and T (z = 1) = 0. Furthermore, we choose free-slip boundary conditions for the velocity field in z-direction, u z (z = 0, 1) = 0 and ∂u x /∂z = 0 at z = 0, 1. Periodic boundaries for all fields are taken in the horizontal x-direction.The chosen boundary conditions, aspect ratio and Prandtl number correspond to a popular standard case for the Lorenz systems [37], but at a higher Ra and thus fully turbulent in contrast to our previous work.

B. Turbulent fluctuations and heat transfer
A central physical question in turbulent convection flows is that on the turbulent transfer of heat and momentum across the convection layer in dependence on the two dimensionless parameters, the Rayleigh number Ra and the Prandtl number Pr.In response to both, these transfers can be summarized in two further dimensionless parameters, the Nusselt number Nu for the turbulent heat transfer and the Reynolds number Re for the The symbol ⟨•⟩ x,t stands for a combined average with respect to horizontal x-direction with x ∈ [0, Γ] and time t.The root mean square (rms) velocity is given by u rms = (⟨u where the average is now a combination of averages with respect to the simulation plane A and time t.The two terms in the definition of the Nusselt number stand for two heat currents across the layer, the convective and the diffusive one.Their sum is constant and equal to Nu for each height z ∈ [0, 1].However their contributions to the total turbulent heat flux differ with respect to z, caused for example by the boundary condition u z = 0 at the bottom and top.
It is exactly the mean profile of the convective heat flux ⟨u z T ⟩ x,t as a function of z, which we want to obtain as a reservoir computing model output.The same holds for the vertical profiles of the root mean square velocity and temperature.This is the low-order statistics of the turbulent flow which is generated by the different reservoir computing models without a knowledge of the nonlinear Boussinesq equations of motion.We will come back to these results in Sec.V D.

III. DATA REDUCTION TO LATENT SPACE
The numerical simulations are performed on a nonuniform grid of size of 384 × 96 points with a 2nd-order equidistant time stepping of ∆t = 5 • 10 −4 .For the analysis, the simulation data were interpolated to a uniform grid of size N x × N z = 128 × 32.The data set consists of a sequence of snapshots of the fields u(x, z, t m ) and T (x, z, t m ) with m = 1, 2, ..., 10 4 ; they are equidistant in time with τ = 0.25H/U f .The sequence covers the statistically stationary regime of the turbulent convection flow.The three turbulent fields possess an input vector of size 12288.In order to circumvent high computational effort for the machine learning, we add a pre-processing step.
Here, we apply a POD in the form of a snapshot method [43,47,48].It is a linear method, where the data reduction is realized by a truncation to a set of Galerkin modes.For this, we decompose the physical fields into time mean and fluctuations, Finally, we perform the snapshot POD to the fluctuating component fields g k into time dependent coefficients a i (t) and spatial modes Φ i (x, z), such that the truncation error is minimized.The degrees of freedom Ñdof can then be reduced, by taking only Ñdof ≪ N dof modes and coefficients with the most variance into account, with g k = {u ′ x , u ′ z , θ ′ } and thus k = 1, 2, 3. Figure 3(af) compares the reconstruction of the temperature and velocity fields from 3 and 16 modes with the original simulation data at a time instant.The time series of the expansion coefficients of three first POD modes, a 1 (t) to a 3 (t), will be fed into the recurrent network in the reconstruction phase after training.
In the following, we use the more compact notation a t = (a 1 (t), a 2 (t), ..., a M (t)) with M = N dof and the discrete snapshot time superscript t.This is our dynamical system state which has to be learned by the reservoir computing model.We choose the cutoff at N dof = 16 and capture 87% of the variance of the original fields as seen in Fig. 3(g).It also implies that the notation ρ(t) changes to ρ t and so on, see again eq.(1).

IV. RESERVOIR COMPUTING FRAMEWORKS
A reservoir computing model [24,25] is one realization of recurrent machine learning beside long shortterm memory networks or gated recurrent units [49,50].FIG.4: Quantum circuits of the two applied hybrid quantum classical algorithms which represent the quantum computing version of the reservoir.Left: hybrid quantum-classical reservoir computing algorithm H1, which represents our previous approach of ref. [37] with a new hyperparameter l to denote the circuit complexity (or depth).Right: hybrid quantumclassical reservoir computing algorithm H2 which we will in short denote to as the full quantum reservoir computing model.The new structure incorporates an efficient pseudo-initialization of the last reservoir time step as well as an efficient memory implementation by the leaking rate ε, which now scales the dynamical activation of the approach without classical postprocessing.In both cases, every qubit of the quantum register starts in the basis state |0⟩ to the left.Each quantum circuit is shown for 7 qubits and l = 3 layers of RY rotation gate at a discrete time step t.The CNOT gates between the RY gate connect neighboring qubits and generate entanglement.Since this pairwise 2-qubit operation is done successively over all qubits, one ends up with a fully entangled 7-qubit quantum state.The tilde symbol indicates that the rotation gate arguments are re-scaled to cover the respective range, as given in the text.Also, in the right circuit diagram init(ω) = RY (2 arccos(ω)).The last column of symbols (to the right) in both panels stands for the measurement to read out the reservoir state.
Its fundamental element is the reservoir state vector r t ∈ R Nres , which evolves from t to t + 1 by a certain update equation characterizing the approach.As indicated in the introduction, this work compares three RCMs, algorithm C a classical RCM [44,45,48,[51][52][53] with linear memory and nonlinear activation, see eq. ( 10), algorithm H1 a hybrid quantum-classical RCM [37] with classical linear memory and nonlinear quantum dynamics, see eq. ( 12), and algorithm H2 a full quantum RCM which induces memory by reduced gate parameterization, see eq. (14).Figure 4 illustrates the different quantum circuits of H1 (left) and H2 (right).All algorithms run here in the open-loop or reconstruction mode [54], where the reservoir receives a t in = (a t 1 , a t 2 , a t 3 ) at each time step and outputs the missing expansion coefficients at the next time step, ât+1 out = (â t+1 4 . . ., ât+1 N dof ) with N dof = 16 for our study.The hat symbol identifies always the network prediction in the following.Note that in this scenario the reservoir receives the coefficients corresponding to the large-scale fluid motion while returning the small-scale features of the higher modes.

A. Classical reservoir computing algorithm
The classical algorithm C is characterized by the following iterative equation for the reservoir state vector r C [45,48,55], where W res ∈ R Nres×Nres is the sparsely occupied random reservoir matrix and W in ∈ R Nres×Nin is the random input matrix.The scalar ε ∈ [0, 1] is called leaking rate, as it scales the influence of the first memory term and the second dynamical term.This discrete time stepping from snapshot t to t + 1 comprises external forcing by the inputs a t in as well as a self-interaction with the reservoir state r t C .The hyperbolic tangent tanh (•) is the nonlinear activation function, which is applied to the elements of its argument vector.The reservoir matrix W res is specified further by the reservoir density D, its sparsity, and the spectral radius ϱ(W res ), the magnitude of the maximum eigenvalue.The hyperparameters of the classical RCM are thereby N res , ε, D and ϱ(W res ).In all reservoir computing models, the optimization is done directly for the output layer only which is represented by the matrix W out * ∈ R Nout×Nres .The asterisk stands for the optimized matrix after training, see refs.[37,48] for the details.Also there, the Tikhonov parameter γ is explained further, which is used as a prefactor of an additional penalty term in the cost function.This hyperparameter is only applied in the classical RCM.
To give an explicit example for the one-step reconstruction mode in which the model will be used here: the classical RCM utilizes eq. ( 10) to calculate r t+1 C from r t C and a t in with its three components to obtain the estimate with ât+1 out = (â t+1 4 , . . ., ât+1 16 ) in the present turbulent convection flow study.

B. Hybrid quantum-classical reservoir computing algorithm 1
The algorithm H1 follows the classical reservoir iteration procedure (10) closely, but the update equation receives another dynamical part.We introduced H1 in [37].It is given by with where |e j ⟩ is the j-th basis vector of the n-qubit Hilbert space C ⊗n with a dimension N = 2 n .Furthermore, |e 1 ⟩ = |0⟩ ⊗n is the n-qubit basis vector for which every of the n individual qubits is in the base state |0⟩.The quantum circuit U (ζ, a t in , r t H1 ) is parameterized with random values ζ, the current input a t in and the last reservoir state r t H1 .We initialize r 0 H1 as a uniform-amplitude probability vector at all entries, which guarantees that r t H1 will remain a probability vector for all times t.
The structure of the circuit is a repeated pattern of R Y -gate and separating CNOT-gate layers, where the arguments of the R Y -gates are the only difference between the layers.The circuit always starts and ends with an R Y -gate layer, as shown in Fig. 4 Note that the loading of the reservoir state and the dynamical system modes by rotations into H1 and H2 introduces a nonlinearity in the (linear) quantum dynamics.
The CNOT-gate layer uses n CNOT gates such that every qubit is control and target qubit once.The specific arrangement is random but fixed.We always assure that it is impossible to separate any subgroup of qubits from the remaining one; see [37] on the importance of the entanglement for the reconstruction quality measured by the mean squared error.The sorting of the R Y gate arguments is such that the last layer is filled with random values, all other layers receive probabilities, and the first three entries in the second layer are always the inputs.There are too many possibilities to proof that this specific circuit construction is the optimum, but we tested many architectures and choose the described due to best overall performance.

C. Hybrid quantum-classical reservoir computing algorithm 2
The new hybrid algorithm H2 is modified such that the complete execution on a quantum computer is enabled.It is given by with An identity transformation is imposed by This is a consequence of the inclusion of the leaking rate ε in the arguments for H2, the major difference to H1.Furthermore, the approximate reservoir vector is where P t i,|1⟩ denotes the probability of measuring qubit i of the whole n-qubit register in basis state |1⟩ at reservoir time step t.In other words, the two probabilities, 1 − P t i,|1⟩ and P t i,|1⟩ are the diagonal elements of the 2 × 2 density matrix ρt i which is obtained by the following partial trace of the original density matrix of the n-qubit state, ρ t , which traces out all qubits except qubit i.It is given by In a quantum algorithm, this is realized by an individual measurement at qubit i of the whole quantum register only.We thus structure the circuit such that it is initialized by a completely separable approximation of the last reservoir state and thus ease the reservoir initialization at each step, see refs.[23,56] for alternative solutions to circumvent this bottleneck in analogue quantum reservoir computing.This combines two advantages: It is a minimal initialization for the quantum circuit with only n operations which contains the integral information on the reservoir.Combined with the adapted definition of the unitary U in eq. ( 14), we include the previously external memory in the quantum circuit iteration step.That is, the dynamics will correspond to an identity transformation, see eq. ( 15), once the leaking rate is ε = 0.
Here, the structure of the quantum circuit starts with the preparation of |r t H2 ⟩, which comes down to an R Ygate layer parameterized with the P t i,|1⟩ , as done in upcoming eq. ( 23).The following circuit layer has pairwise CNOT-gate layers, where every second CNOT layer is the inverse of the first CNOT layer, as indicated in the right panel of Fig. 4. Thereby, we satisfy the central condition of eq. ( 15), as R Y (0) = I 2 .Beside the input gates, all R Y gates are filled with random values ζ k .Further details of the circuit are adopted from H1.Note however, that the circuit ends with a CNOT layer for an even number of R Y layers.

D. Advantages of hybrid reservoir computing approaches
There are several aspects which motivate the presented hybrid quantum-classical approaches H1 and particularly H2.First, we want to investigate whether the nonlinear activation of the quantum circuits is preferable over tanh(•) of eq.(10).The inherent nonlinearity of product concatenations from the unitary gates to larger matrices approximates the inherent nonlinearity of the original flow problem, see the second term on the left hand side of eq. ( 3).
Secondly, we note that for online processing of the data, we can strongly reduce the required memory once the output of a time step is generated.For instance, one needs to save the used probabilities only, i.e., for H1 the r t H,i and for H2 the P t i,|1⟩ , to recompute the reservoir state at the subsequent time step again.This would not be possible for the classical reservoir.Especially for H2, this always means that for a reservoir of size N , one keeps log 2 (N ) values only to store sufficient information of the reservoir evolution.If necessary, it is then possible to reconstruct all 2 n reservoir states in parallel.
Thirdly, H2 is in contrast to H1 fully executable on a quantum computer, as the leakage rate of the reservoir step is included in the rotation gate parameters of the quantum circuit.Therefore, no external memory is required, cf.eq.(12a) for H1.Additionally, it might be possible to extend the present scheme H2 to a multistepping approach on the quantum computer if one implements the output weights as well.This could be done by encoding the output on an ancilla qubit to either measure it directly or insert it back for the next time step via a rotation controlled by this ancilla qubit.Future work in this direction might elaborate whether this approach can strongly reduce the computational complexity in comparison to the classical approach by reducing the number of sampling steps (or shots) for multiple time steps in the hybrid case.

V. COMPARISON OF THE MODELS A. Error quantification measures
In order to evaluate the reconstruction quality of the proposed RCMs, we need appropriate error measures.The first standard is the root mean squared error of the prediction.Since multiple modes have to be reconstructed, we take the normalized error with respect to each mode and combine the N out = 13 individual errors of the reconstructed modes to the normalized root mean squared error (NRMS).This results to where T = t b − t a is the length of the testing phase (measured in discrete time steps as discussed in Sec.III).
The maxima and minima are determined with respect to the time interval [t a , t b ].  18)-( 20) and three amplitude encoding methods, which are defined in Eqns.( 22) to ( 24) for both hybrid RCM cases, H1 and H2.They are termed R1, R2, and R3.Here, each data point is the median obtained from ten different random reservoir seeds ζ.The hyperparameters l and ε are optimal in each case.For the optimization, 2 ≤ l ≤ 7 and 0.05 ≤ ε ≤ 1.We trained for 5000 and evaluated on the subsequent 500 discrete time steps.
A second popular approach is the correlation error, also known as the coefficient of determination, or R2-score, which computes the correlation between the original and predicted modes [23].Here, we average the square of correlation, such that the correlation error is given by Here, σ(•) is the standard deviation and cov(•, •) the covariance of the arguments.Note that by the square of the correlation, we value anti-correlation as much as positive correlation, thus 0 ≤ E corr ≤ 1. Strongly correlated time series send E corr → 0. Both error measures are applicable to any dynamical system.However in the present work, we consider a turbulent flow; the RCM application is focused to reconstructed statistical properties, as motivated in the introduction and in Sec.II.Such properties can be the mean vertical profiles of the velocity components or the temperature field.Therefore, we define an additional measure which is directly related to the low-order statistical reconstruction results, the normalized average relative error (NARE), which has been used in classical RCM applications [48,57] and is given by the L 1 norm, with the normalization constant Here, ûz and T are the reconstructed flow values which are obtained by the Proper Orthogonal Expansion (POE) from the modes ât j .As the convective heat flux is prone to the error of two fields, it is a suited measure of the accordance of the inferred convection flow.We compare these three errors for the hybrid quantum RCMs H1 and H2 in Fig. 5 as a function of the reservoir size controlled by the qubit number n and for different amplitude encoding methods, which will be detailed in the next subsection.
First, it can be observed that E corr is relatively large with a minimum for 8 to 9 qubits.The reason is that the accurate reconstruction of the POD modes is difficult as frequent deviations of the time series ât i from the ground truth are inevitable in this higher-dimensional, turbulent flow problem.In contrast to low-dimensional dynamical systems, such as the Lorenz model, a reservoir computing model will not reconstruct the exact systems trajectory in the high-dimensional phase space with a sampling time step of 0.25t f , but generate a trajectory which gives the right low-order statistics.This is even the case for the classical RCMs [45,48,52,57].Thus it is not appropriate to optimize the network on the basis of E corr .Particularly weak correlations are not directly linked to the dynamical quality of the reconstruction.
Secondly, we observe from the figure that both, E NRMS and E NARE , grow eventually with a large qubit number n.However, it can be observed that particularly for 8, 9, and 10 qubits, the physical error E NARE improves by almost one order of magnitude while E NRMS remains relatively constant and insensitive.We thus conclude that while the accuracy of the reconstructed modes is difficult to improve further, the physical properties of the flow are more sensitive.Thereby, for most of the remaining analysis, we will continue our RCM evaluation with E NARE only, as it is the physically relevant measure for the fluid mechanical application of the quantum algorithm.

B. Different amplitude encodings of classical data
Besides the hyperparameters, which will be discussed further below, the quantum circuits need to encode the classical input data a t j .This can be realized in different ways.We discuss the following three amplitude encoding methods, The tilde symbol in the equations indicates again, that the input mode a t j needs to be re-scaled such that it only varies in the interval [−1, 1].Encoding R 1 ensures that ãt j is the component of the corresponding qubit, while encoding R 2 reveals ãt j after a measurement.Encoding R 3 is a natural encoding inside the R Y -gate, where we only re-scale the input to harness the largest but still unique range of the trigonometric functions.
Each approach induces a specific nonlinear characteristic and the superiority of each encoding may change for different learning tasks.We come back to Fig. 5 where we plotted all three error measures versus qubit number n for R 1 to R 3 .The error measure E NRMS shows an approximate independence on the specific encoding scheme for n ≤ 9 for both H1 and H2.Only for the larger qubit numbers R 3 performs best.In case of E NARE a local minimum can be observed for each encoding scheme.All three amplitude encodings have their optimum at a magnitude of approximately 10 −2 .This is obtained for R 1 and R 2 already at a smaller number of 8 qubits.Again, for the lower reservoir dimensions of n = 6 and n = 7 all errors are of the same order of magnitude.As already said, for the largest reservoir dimensions, the error increases, similar to E NRMS , for n ≥ 10.In conclusion, we do not fix the specific amplitude encoding for the following analysis, but use it as a further degree of freedom to be optimized.

C. Hyperparameter dependence
Table I summarizes all hyperparameters that appear either in the classical or in the hybrid quantum-classical RCMs.The hyperparameters, which exist in the quantum case only, are the circuit layer depth l and the type of amplitude encoding.The latter was already discussed in the past subsection.Joint parameters of the classical The number of used qubits follows by n = log 2 (Nres) in H1 and H2.The hyperparameters, which are not plotted, in the figures are always pre-optimized, that is, we take the optimal values such that the single-best representation is illustrated.
FIG. 6: Evaluation of the NARE as a function of circuit depth l for the hybrid quantum cases H1 and H2.Each data point is the median of ENARE, obtained from ten random reservoir seeds, using the optimal leaking rate ε and optimal input amplitude encoding.In nearly all cases, the optimal values are ε = 0.2 and R1 for the input encoding.
and quantum case, which will be studied in the following, are the reservoir size N res and the leaking rate ε.We also mention at this point that all studies for H1 and H2 are conducted with the IBM Qiskit statevector simulator where the measurements are not subject to shot noise [58].We train all cases for 5000 output time steps and validate the trained network on the subsequent 500 output time steps.

Quantum circuit layer depth
The quantum circuits of H1 and H2 can be parameterized by the amount of R Y -layers l and the positions of FIG.7: Comparison of the error ENARE for the convective heat flux ⟨uzT ⟩x,t versus the leaking rate ε.We show the results for the three reservoir computing approaches C, H1, and H2 and for different reservoir sizes 2 n .Displayed are the median values for 10 seeds with H1 and H2 and 100 seeds for C. The optimum of each curve is the single-best median for the respective approach at the given reservoir size.Note that the best classical reservoir performance is achieved at substantially larger reservoir sizes (n = 11) as compared to the quantum algorithms (n = 8 for H1 and n = 9 for H2).
the CNOT-gates.The later aspect showed no significant influence during our analysis, that is, the performance of both algorithms is relatively constant as long as all qubits are connected, as described in Sec.IV B. Meanwhile, the depth of the circuit is of interest, as it is directly proportional to the computational effort of the approach, i.e., the number of operations, as well as the realizability on real quantum computers.Figure 6 displays the dependence of E NARE on the layer depth for H1 and H2.We collect the results for n = 7, 8, 9, and 10.Except for n = 10, the overall trend of the error E NARE is a decrease with growing circuit depth l.For n = 10, we observe a clear advantage of H2 compared to H1.For n = 8 and n = 9, the cases H1 and H2 perform similarly well, almost at their optimal error measure.We finally mention that the median was taken over a rather small number of different reservoir realizations.

Reservoir size and leaking rate
In Fig. 7, we show the median of E NARE in dependence on the leaking rate ε for different reservoir sizes N res = 2 n .For all points, we averaged here over 100 seeds for C and 10 seeds for H1 and H2, while all other parameters are pre-optimized, that is, we choose the optimal spectral radius ρ(W res ) and the Tikhonov parameter γ in the classical RCM case, the optimal input encoding R i and the amount of layers l in the hybrid quantumclassical cases.We choose this optimum such that the single-best median is illustrated for the respective approach and reservoir size.We observe that H1 and H2 seem to outperform the classical approach for qubit numbers n < 10.The global optimum, i.e., the minimal amplitudes of E NARE , are obtained for the new architecture H2 at n = 9, though the other RCMs can perform simi-larly well if the reservoir is large enough.

D. Statistical analysis of the reconstructed fields
Of particular relevance in turbulent convection is the mean vertical convective heat flux profile, which is the one-point correlation of the vertical velocity component and the temperature field, ⟨u z T (z)⟩ x,t , see Sec.II B. It is a measure of the amount of heat transported by fluid motion from the bottom of the layer to the top of the layer.Such a vertical profile is more difficult to reconstruct, as it combines the statistics of two reconstructed fields.In addition, we will monitor root mean square profiles of the velocity components and the temperature.These are essential low-order statistical properties of the flow at hand.They are also important when the turbulence cannot be modeled down to the smallest physically relevant scale and has to be parametrized.This is the case for in subgrid-scale parametrizations in global circulation models, e.g. in atmospheric turbulence [41].
We illustrate in Fig. 8 the single best reconstruction of each RCM model combined with the mean statistical profiles of the flow.We evaluated the previous grid search over the hyperparameters for the single best results.Figure 8(a) displays the reconstructed time series of four POD expansion coefficients from C, H1, and H2 in comparison to the ground truth (GT).It is seen that the curves are not followed exactly, but that the overall trends and thus the low-order statistics are represented well.Note, that this is not only for the case for H1 and H2, but also for C. Illustrated are the best cases for the physical error measure E NARE .In other words, we optimized for the lower part of the figure.
The mean convective heat flux profile in the most right panel of Fig. 8(b) is the most sensitive, which is the rea- show vertical mean profiles of essential low-order statistical measures.These are the thermal variance ⟨θ ′2 ⟩x,t, the horizontal velocity fluctuations ⟨u ′2 x ⟩ 1/2 x,t , the vertical velocity fluctuations ⟨u ′2 z ⟩ 1/2 x,t , and the convective heat flux ⟨uzT ⟩x,t.
FIG. 9: POD reconstruction based on the latent space variables in Fig. 8.We show the horizontal and vertical velocity fluctuation fields, u ′ x in panels (a,d,g,j) and u ′ z in panels (b,e,h,k).Furthermore, the local convective heat flux uzT is displayed in panels (c,f,i,l).The fields of the classical reservoir algorithm reconstruction are shown in panels (d-f), while the quantum algorithm results for H1 and H2 are shown in (g-i) and (j-l), respectively.For comparison, the ground truth, i.e., the POD model, is shown in panels (a-c).It can be seen that both quantum reservoirs, as well as the classical counterpart, reproduce the spatial features of all three fields of the ground truth (GT) well.
son why we utilized it for the hyperparameter optimization.This statistical correlation is connected to the hot rising and the cold falling plumes which are visible in Fig. 3(e).Figure 9 displays finally the spatial reconstruction of velocity and heat flux fields for all three cases.We find that both quantum algorithms, as well as the classical counterpart, produce statistical profiles that are in good accordance with the ground truth, for both, the root mean square profiles of the velocity components and temperature as well as for the convective heat flux.This demonstrates the applicability of the hybrid quantumclassical reservoir computing model as a reduced-order model in combination with the encoder/decoder module in the form of POD/POE, respectively.

VI. FINAL DISCUSSION
Two central objectives can be given for the present proof-of-concept study.First, we wanted to extend the application of hybrid quantum-classical reservoir computing algorithms towards more complex classical dynamical systems.Starting with the well-known Lorenz 63 benchmark case and its extension to 8 degrees of freedom in [37], we increased here the complexity of the task to be learned further by proceeding to a turbulent convection flow at the same geometry and Prandtl number as in the Lorenz cases, but at a significantly higher Rayleigh number (the latter of which measures the driving by buoyancy forces of the flow).We integrated the quantum circuit therefore into a combined encoder/decoderreservoir computing pipeline which has to be used as well when classical machine learning is applied to turbulent flows [44], see again Fig. 1.Since the phase space of the Rayleigh-Bénard system is higher-dimensional and thus the dimensionality of the turbulent attractor, the reservoir computing algorithm is only able to predict loworder turbulence statistics rather than exactly following a specific dynamical systems trajectory for a longer time.But this is exactly the task that we had in mind, reproducing 2nd-order statistics, such as the convective heat flux, in a data-driven model without solving the full nonlinear partial differential equation system of Rayleigh-Bénard convection.
The second objective is related to the modified architecture of H2 in comparison to H1.We have compared both hybrid quantum-classical algorithms with respect to various hyperparameters and found that they mostly perform equally well for the reconstruction tasks.Nonetheless, the update of the circuit architecture H2 can be evaluated completely on a quantum computer, which enables further steps of the hybrid algorithm on the quantum device and avoids additional external memory as H1.Additionally, the simulation with the circuit architecture H2 can be realized more efficiently than the one for H1 once the layer depth is l > 3. The reason is that in H2 every operation that follows the input encoding can be summarized to one pre-computable matrix which acts on the time-dependent inputs.This is not possible for the architecture H1, as further subsequent circuit layers also have to be filled with the time-dependent components of the probability amplitude vector p t .We also investigated, if a random value encoding which is used in H2 works for the original hybrid algorithm H1, but it was found that this method strongly impairs the performance of H1.It can be concluded therefore that H2 is a more efficient implementation of the reduced-order model of the turbulent convection flow by means of a hybrid quantum-classical reservoir computing algorithm; it is comparable to the best-case scenarios of the classical reservoir computing approaches which however need at least twice as large reservoirs in the present case, see Fig. 7.In future applications, this could reduce the numerical effort for both, hyperparameter optimization and production runs.
The demonstration of the capabilities and potential of the present framework, but also its current limitations, has to our point of view its value for following studies in this subject in the future applied to realistic fluid flow problems.
A first open point for our future work is to solve the sampling problem.We used the ideal Qiskit statevector simulator for both algorithms, H1 and H2, which circumvents the crucial problem of approximating the necessary probabilities [58].A deeper analysis of this aspect shows that the computational overhead to approximate the probabilities, both, for the best cases of H1 and H2, is big.In detail, more than 2 20 samples (or shots) are necessary for comparable results.This seems to damp the prospects for an application on current noisy intermediate-scale quantum devices.However, a repetition of the hyperparameter grid search with sample-based probabilities and the additional implementation of weak measurements as in ref. [22] might ease this problem.Furthermore, it has to be evaluated if the hybrid reservoir computing approach can be further scaled up to more vigorous turbulence, i.e., flows at higher Rayleigh number.This would imply a higher dimension of the latent data space and possibly a different encoding/decoding scheme, see [44,45] for the classical cases.These investigations are going on and will be reported elsewhere.
other words, the qubit can be consequently found in infinitely many superposition states, all the points that fill the surface of the unit sphere.It is the building block of an n-qubit system, also denoted as an n-qubit quantum register.They are formed by successive tensor products of qubits.For example, a two-qubit state vector is the tensor product of two single-qubit vectors, The basis of this tensor product space is given by 4 vectors: An n-qubit state vector is called fully separable if it can be written as a tensor product, where |q i ⟩ are single qubit quantum states given by eq.(A1).It is called separable if a tensor product decomposition of |Ψ⟩ into blocks is possible with at least one multi-qubit quantum state |q i ⟩, that is not fully separable.Multi-qubit quantum states which are not separable are called entangled.An n-qubit quantum state is called fully entangled if no subspace of separable qubits exists.Fully entangled quantum states are characterized by correlations which do not classically exist.The entanglement is a unique property of quantum computing; it is supposed to be responsible for the quantum advantage of some quantum algorithms with respect to their classical counterparts, such as prime factorization [36], as single qubit operations act non-trivially on a large quantum state and produce global parallel processing by a single computational step.
The time evolution of a quantum state is described by a unitary transformation, |Ψ⟩(t) = U (t)|Ψ⟩(0) with U (t) −1 = U (t) † .(A5) Elementary unitary transformation are supplied by quantum gates, for example rotations of a single qubit.Note that they are reversible transformations.The R Y rotation gate is defined by eq.(13)   well as for the unitary evolution of the same.The quantum version of the reservoir is composed of exactly these gates.The readout of information is done by a measurement process which causes the collapse of the n-qubit quantum state.

FIG. 1 :
FIG. 1: Pipeline of data processing in the classical and hybrid quantum-classical reservoir computing models.(a) Sketch of the whole algorithm.Direct numerical simulation snapshots of the turbulent convection flow to the left are encoded in a low-dimensional latent data space by a Proper Orthogonal Decomposition (POD).A reservoir computing model operates in the latent space receiving the input sequences ai(tn) with i = 1, ..., Nin < Nout.Its output is a sequence aj(t n+k ) with j = 1, ..., N dof ; typically one takes k = 1.The output is finally decoded by a Proper Orthogonal Expansion (POE) to simulation data at later times and/or a subsequent statistical analysis.The reservoir takes two different forms in this work.(b) The classical reservoir, which we use for comparison, consists of a sparse random network of neurons that possesses a dynamic memory for temporal data processing.(c) The quantum reservoir is a digital gate-based quantum circuit.In this work, two different quantum algorithms are employed for the quantum reservoir.Both, classical and quantum reservoir contain input and output layers to the left and right, respectively.The reservoir architectures are further detailed in Sec.IV.

FIG. 2 :
FIG. 2: Two-dimensional Rayleigh-Bénard convection flow.The fluid with a kinematic viscosity ν and temperature diffusivity κ is confined between two impermeable boundary planes at the top (z = H) and bottom (z = 0).The bottom plane has a uniform temperature T bot ; the top one a uniform temperature Ttop < T bot .The fluid layer has a length L = 2 √ 2 and a height H = 1, such that the aspect ratio Γ = L/H = 2 √ 2 results.The figure displays a snapshot of the turbulent flow.Background contours stand for the temperature field T ; the vectors illustrate the 2-component velocity field (ux, uz).The turbulent heat transfer across the layer from the bottom to the top is manifest by the thin plume structures which are visible: hot lighter fluid (in red) rises to the top while cold heavier fluid (in blue) falls downwards.Since the velocity field is divergence-free a circulation roll pattern develops.

FIG. 3 :
FIG. 3: Snapshot of two-dimensional convection flow.Panel (a,c,e) show the temperature field T (x, t0), panels (b,d,f) the corresponding velocity vector field u(x, t0) = (ux(x, z, t0), uz(x, z, t0)).Panels (a,b) display the reconstruction by means of the first 3 POD modes, (c,d) of the first 16 POD modes, and (e,f) display the full snapshot.(g) Cumulative spectrum of the Proper Orthogonal Decomposition.The shaded region indicates the chosen cutoff mode at N dof = 16.The first three POD modes include 63% of the total energy, the first 16 modes include 87%.
(left).As a new hyperparameter, we introduce the depth l as the amount of R Y -gate layers inside the circuit.An R Y -gate layer applies an R Y -gate on every qubit, where the arguments are scaled versions of the aforementioned variables.The random values ζ vary in [0, 4π], the probabilities r t H1 are multiplied with π to vary in [0, π] and the inputs a t in are re-scaled to vary in [−π, π].For example, the random values ζ can build various unitary matrices, as the matrix definition of the R Y (φ) gate is given by

FIG. 5 :
FIG. 5: Comparison of three error measures, which are defined in eqns.(18)-(20) and three amplitude encoding methods, which are defined in Eqns.(22) to (24) for both hybrid RCM cases, H1 and H2.They are termed R1, R2, and R3.Here, each data point is the median obtained from ten different random reservoir seeds ζ.The hyperparameters l and ε are optimal in each case.For the optimization, 2 ≤ l ≤ 7 and 0.05 ≤ ε ≤ 1.We trained for 5000 and evaluated on the subsequent 500 discrete time steps.

FIG. 8 :
FIG. 8: Single-best reservoir computing results for Nin = 3.We compare the ground truth (GT) with a classical reservoir computing model (C) with Nres = 2048 and the hybrid quantum-classical quantum algorithms (H1, H2) with Nres = 512 and 6 layers in both cases.All networks output 500 time steps.Panels (a) display the time series of modes a5(t) to a8(t).Panels (b)show vertical mean profiles of essential low-order statistical measures.These are the thermal variance ⟨θ ′2 ⟩x,t, the horizontal velocity fluctuations ⟨u ′2x ⟩

TABLE I :
List of hyperparameters for the classical and both hybrid quantum-classical algorithms.
in Sec.IV B. A second central gate is the controlled NOT gate (in short CNOT) which connects two qubits.It represents a flip of the target qubit once the control qubit is in state |1⟩.The logical table of the two-qubit CNOT gate is shown in table II which can be transformed into a 4x4 unitary matrix.Rotation and CNOT gates are elementary gates which are composed to quantum circuits that are required for the input of the classical data into a quantum algorithm as

TABLE II :
Logical table of CNOT gate.