Dissipation as a resource for Quantum Reservoir Computing

Dissipation induced by interactions with an external environment typically hinders the performance of quantum computation, but in some cases can be turned out as a useful resource. We show the potential enhancement induced by dissipation in the ﬁeld of quantum reservoir computing introducing tunable local losses in spin network models. Our approach based on continuous dissipation is able not only to re-produce the dynamics of previous proposals of quantum reservoir computing, based on discontinuous erasing maps but also to enhance their performance. Control of the damping rates is shown to boost popular machine learning temporal tasks as the capability to linearly and non-linearly process the input history and to forecast chaotic series. Finally, we formally prove that, under non-restrictive conditions, our dissipative models form a universal class for reservoir computing. It means that considering our approach, it is possible to approximate any fading memory map with arbitrary precision.


Introduction
Quantum science holds promise to revolutionize the technological perspectives in fields such as computing [1,2], communication [3], sensing [4], or cryptography [5].In particular, the introduction of quantum logic allows the development of algorithms that outperform the corresponding classical ones [6,7,8,9,10], with broad and interdisciplinary applications in chemistry [11], finance [12,13] and machine learning [14].Still, the experimental implementation of these algorithms is a challenging task with state-of-the-art technology [15,16].Currently, available devices mostly rely on a few (up to hundreds) of noisy qubits, an obstacle that must be overcome to achieve a quantum advantage in "real world" computational problems or quantum error correction.
A common problem is that dissipation, caused by the interaction between qubits and the external environment, produces decoherence phenomena that hinder fragile quantum resources.Dissipation, however, can be turned into a positive computational resource in some cases as for dissipation engineering, which exploits the systemenvironment interaction as an integral part of the computation process [17] and for quantum memories [18].Applications range from quantum control [19] to non-equilibrium quantum thermodynamics [20,21], quantum biology [22] or quantum synchronization [23,24].The aim of this work is to establish the beneficial effect of tunable dissipation on quantum machine learning and, in particular, on the field of quantum reservoir computing (QRC) [25].
Reservoir computing (RC) is a machine learning method rooted in Recurrent Neural Networks and is particularly suited for time series processing [26].Classical RC was initially proposed as Liquid State Machines [27] and Echo State Networks [28] and later generalized to include a broad spectrum of physical implementations [29].Two important advantages of RC are the simplicity of the training process, which requires only linear regression optimization rather than the gradient descent procedures of traditional neural networks, and the ability to multitask.RC is now wellestablished as a powerful analog neuromorphic method in classical machine learning.Successful experimental realizations of temporal processing with RC and applications have been reported in the last years, including electronic, spintronic, and photonic systems [30,31,32,33].
In spite of the variety of dynamical systems that can serve RC purposes, there are some key features needed for temporal series processing in this architecture (formally defined in Sec. 3).In a nutshell, the system is required to store some past input information (fading memory) and to forget its initial conditions (echo state property).For quantum RC, it is known that the presence of dissipation is a necessary feature [34,36] and, consequently, all the proposed QRC frameworks contain some dissipation.
The pioneering work of Fuji and Nakajima (FN) is based on a spin-network QRC scheme [34] and a discrete erase-and-write map.This alternates a unitary evolution followed by an instantaneous input encoding (resetting one qubit).The erasure can be realized via a local dissipation restricted to the input node(s) and modeled by a Lindblad master equation (Sec. 2 and Appendix A).An alternative paradigm for QRC that we propose is based on the presence of engineered losses for each of the spin nodes, giving rise to continuous dissipation (CD), acting at all times, and on a continuous input drive obtained by tuning an external magnetic field.We will show that the degree of adaptability of such a tailored dissipation allows the system to optimize its performance according to the task faced, thus outperforming the FN model, where the control over dissipation is very low.Furthermore, we will demonstrate that the CD model achieves universal QRC, which means that any generic task to be solved with RC can be arbitrarily well approximated by considering only these kinds of models.
The paper is structured as follows: in Sec. 2, we review the FN model and introduce the CD one, discussing its universality in Sec.3; in Sec.4.1, we analyze the memory properties of such models while in Sec.4.2 we test them in timeseries forecasting tasks; in Sec. 5, we study the effects of limiting the number of measurement samples on the performance of the systems, and we present possible experimental platforms for implementing the CD model; finally, discussion and conclusions are given in Sec. 6.

Quantum reservoirs
According to the RC theory, a general model of a reservoir must satisfy some necessary conditions to operate.One is the echo state property, which consists of the disappearance of the dependence of the initial condition on the reservoir dynamics over time [58].This is a necessary feature because otherwise, the training phase would also have to consider the initial state choice, making the whole procedure inefficient.Furthermore, a proper RC system needs fading memory to process the information of a time series without requiring an infinite amount of physical resources [59].Therefore, in the quantum case, dissipation is a crucial feature to satisfy these two conditions [36].Finally, an RC model must also be able to discriminate between different input sequences in order to adapt its behavior to the particular problem of interest [60].We will formalize and contextualize all these features for our case study in Appendix B. In this section, we will consider two different quantum reservoir computers focusing on the role played by dissipation in their functioning and performance.In all the cases, we will consider spin network-like reservoir systems, whose output layer is given by a linear combination of local and/or global observables.Moreover, as we will show, all the models presented are designed to process discrete-time signals.
Let us start our discussion by recalling the first model of QRC introduced by Fuji and Nakajima (FN) in Ref. [34] and based on a transverse-field Ising model characterized by the Hamiltonian where the i, j label the sites of the network, σ a i (a = x, y, z) are the Pauli matrices acting on the i-th site, h is the value of the homogeneous magnetic field and J ij is the spin-spin coupling following a uniform distribution in a pre-determined interval [−J s , J s ].We will consider a real input sequence of length M , {s k } M k=1 , and rescaled, so that s k ∈ [0, 1] ∀k.The FN updating rule of the reservoir is obtained by feeding the input to the state of one qubit of the network, for the sake of the definiteness we say the first one.In particular, the state of the first qubit is prepared in an input-dependent coherent superposition ρ (see [61] for a discussion of this and different encoding effects).Then, the system unitarily evolves for a certain interval of time ∆t, only following the dynamics generated by the Hamiltonian H.The complete update rule is: where Tr (1) {•} is the partial trace over the first qubit and e −iH∆t is the time evolution operator, assumed as unitary between the input injections.Still, the map (2) exhibits dissipation and decoherence and this occurs instantaneously at each input injection (as modeled by the partial trace).
A natural way to experimentally realize such injection on a NISQ device consists of realizing a measurement on the first spin and, subsequently, setting its state with a quantum gate conditioned by its outcome.In digital implementations like the IBM quantum computer, the reset of qubits after a measure is a recently implemented feature [62] although it is rather slow and then susceptible to uncontrollable decoherence [63].
Let us instead consider an alternative QRC approach characterized by interactions with an external environment under Markovian conditions.The most general time evolution of a density matrix is described by the Gorini-Kossakowski-Sudarshan-Lindblad (GKLS) Master Equation [64,65,66]: (3) where {γ i } are the decay rates of the qubits in each external environment and the operators {L i }, called jump operators, identify the environment action on the qubits.In Eq. (3), we identify unitary U and dissipative D superoperators: with U[ρ] = −i[H, ρ] and D[ρ] equal to the remaining (Lindbladian) term.We will consider dissipation leading to local losses (i.e.independent losses at each of the N reservoir nodes) modeled by N jump operators We notice that this model adequately accounts for Markovian dissipation in independent baths whenever the interaction between reservoir nodes is weak (for a more accurate discussion see [67]).Furthermore, protocols to engineer a Lindbladian with these characteristics are also known [17].In Appendix A, we show how dissipation acting locally on the first oscillator leads to an evolution that is able to reproduce the FN map (2).
Let us now introduce a different model of quantum reservoir computing characterized by continuous dissipation (CD) where the input is injected into the system through temporal driving, varying a Hamiltonian parameter.We consider, in particular, a variable magnetic field modulated in the x-direction into the Hamiltonian (1): where the spin-spin coupling coefficients are randomly chosen from a uniform distribution in [−J s , J s ] as in Eq. ( 1).The time-dependent h ′ (t) encodes the input time series and coherently modifies the evolution of the reservoir.For each input s k ∈ [0, 1], driving persists for a certain interval ∆t according to the assignment rule ).As dissipation is required for QRC [36], the unitary dynamics generated by the Hamiltonian (4) will not be sufficient for our purposes.The simplest kind of dissipation we can introduce consists of adding local, uniform losses to the reservoir nodes (γ i = γ ∀i).We will see that the decay rate γ strongly affects the memory and computational properties of the system.During each time interval, the reservoir dynamics is governed by where the local (L) dissipator is given by D and where Therefore, in this case, U and D L are functions of the input and γ respectively and the updating rule of the reservoir is: where the specific choice of continuousdissipation and input-driven dynamics is distinguished by the bar in L. A schematic representation of this CD model is shown in Fig. 1.

Universality
As anticipated in Sec. 2, a reservoir computing model, to properly work, must fulfill a set of necessary features, namely the echo state property [58], fading memory [59] and input separability [60].These properties clearly determine the class of problems that this machine-learning paradigm is designed to solve.A reservoir computing model is universal if it exhibits these properties.Then, in principle, its outputs can reproduce any fading memory map arbitrarily well [27,68,43,41].Indeed, as in other contexts, universality refers to the capability of a class of systems (reservoir computers here) to approximate any map in a much larger class, with arbitrary precision.As one of our main results, in Appendix B, we prove that the proposed CD model has the universality property.Previous results of universal QRC in other settings have been reported in Refs.[36,41].Beyond the specific model proposed here, our proof can be applied to a wider class of quantum reservoirs that evolve according to a master equation.In short, we showed that if a generic input-dependent generator L(s k ) (i) admits only one stationary state, (ii) is a continuous function of s k and (iii) takes different stationary states for different inputs, then the corresponding quantum reservoir model, for a value of ∆t sufficiently long, is universal.Interestingly, these general conditions can be used for finding new quantum reservoir models even beyond spin network implementations.

Computational benchmark tasks
We now proceed to numerically evaluate the computing performance of the FN and CD models.We will tackle two families of benchmark tasks that usually appear in the RC literature, namely memory and forecasting tasks.We set the reservoir to N = 5 spins, which implies a density matrix with a number of elements ∼ 4 5 .This reservoir size allows us to achieve reasonably good performance in all considered tasks, while still not requiring excessive numerical resources and being accessible on a standard desktop computer.The effect of increasing or decreasing the reservoir size is discussed in Appendix E.
For the readout layer, we choose a linear combination of the expectation values of Pauli strings with length one (⟨σ a i ⟩) and two (⟨σ a i σ a j ⟩), with a = x, y, z, and 1 ≤ i, j ≤ N , i ̸ = j.Altogether, the output layer is made of 45 nodes which are the parameters determined during the training phase.This set of observables is adequate for solving all the tasks we have tackled with high accuracy while utilizing a large number of degrees of freedom that scales polynomially with the number of qubits.
In all the cases studied, we have subjected the systems to a washout phase before carrying out a training and an evaluation phase.It consists in letting the reservoir evolve for a certain number of time steps in order to guarantee the ESP [55].We have found that a number of 1000 time steps is always a sufficient value for the washout.The subsequent 1000 points of the dynamics, were used to train the free weights of the output layer (training phase) [30], while the next 1000 points were used to evaluate the performance (evaluation or test phase).We have numerically verified that the data sets employed in all the tasks considered are sufficiently long to avoid overfitting and to collect meaningful statistics to properly evaluate performance.More details about how we trained the reservoirs of interest can be found in Appendix D.

Memory tasks
In this section, we will present the performance of tasks related to the capacity of the systems to linearly and nonlinearly process the memory of previous inputs, starting with a linear memory test: the short-term memory task (STM) [69].Following the standard procedure, the input originates from a uniform random distribution in the interval [0, 1] and the expected target at each time step (y k ) is to reproduce the previous input for a given delay τ y k = s k−τ .
Indeed, the STM task measures the ability of a system to store information about the input received a certain number τ of time steps in the past, which is an indicator of linear memory.The metric chosen to evaluate all the presented memory tasks is the capacity: where y and ȳ are respectively the time series of the targets and the predictions, cov(•) is the covariance and σ(•) is the standard deviation.The coefficient C ranges between C = 0 (complete mismatch of the predictions) and C = 1 (perfect accuracy).
Our analysis encompasses a coarse-grained exploration to optimize the hyperparameters of the models for all the tasks discussed in this paper.Working in the units of J s , the degrees of freedom h, ∆t, and γ are varied by orders of magnitude, taking values from the following set: {0.01, 0.1, 1, 10}.For each possible combination of these hyperparameters, we simulated 100 different random pairs of coupling sets {J ij } and input sequences and took the metric average over them as a representative value.Finally, the combination with the average of the maximum performance in the given task was considered optimal.
The results of the STM task, specifying the values of the hyperparameters, are shown in Fig. 2 (a).In general, for different orders of delay, the optimal set of free parameters varies for all the models, which makes it evident that the learning capability of the reservoir is strongly influenced by the choice of such hyperparameters.We found that FN reservoir memory is maximized for h = 0.1 and ∆t = 10 when τ < τ * 1 = 5 while in the complementary case ∆t and h change to 1.For our proposal (CD model ( 6)) four different optimal sets have been found.Indicating the free hyperparameters through a triple of the form (h, ∆t, γ), for the regions of delay τ < τ 1 = 2, τ 1 ≤ τ < τ 2 = 6, τ 2 ≤ τ < τ 3 = 16 and τ 3 ≤ τ ≤ 20 the optimal sets are respectively (0.1, 10, 1), (0.01, 10, 0.1), (0.01, 10, 0.01) and (0.01, 1, 0.1).For all the values of the delay, the CD dynamical map (6) is able to reach better performances than the FN model.This suggests that the tunable damping rate introduced in the dynamics, which is an additional free hyperparameter with no equivalent in the FN map, allows us to have more control over the linear dependence from the past injected input.
The generality of this result is also addressed in Appendix E where we report the STM achieved at delay 10 for a smaller (N = 3, 4) and a larger (N = 6, 7) number of qubits.Our results indicate that the advantage in the performance of the CD model is sustained.
Another well-established and challenging benchmark test studied in the context of RC is the nonlinear auto-regressive moving average (NARMA) task, first introduced in the context of RNNs [70].This task, beyond requiring linear memory, also adds the requirement of non-linear memory, and its formulation depends on a given order of delay n.At any time step, the NARMA(n) target is As usual, to prevent divergences until an order of delay equal to 20, it is necessary to define the target of the task rescaling the input in the interval [0, 0.02].The NARMA task performances, for both systems, are shown in Fig. 2 (b).Also in this case the results show that the CD model performs better than the FN one, broadening our conclusions about the STM task, also for a target that includes some nonlinearity.We conclude our analysis of the memory properties considering the parity-check task [71].In this case, we work with a binary random input sequence s k ∈ {0, 1}.The desired output for a delay τ is given by This task is strictly non-linear and, in fact, for a delay equal to τ , it has the degree of nonlinearity of a monomial of degree τ of the previous inputs.In order to reach performances not negligible for τ > 1, we have introduced in both the QRC models the time multiplexing, considering V = 15 virtual nodes.It means that, for each input injection, during the evolution time ∆t, the observables in the readout have been collected not only in the final time but in V equidistant times of the dynamics [34,38].The optimal performances for this task are shown in Fig. 2 (c).Even for this fully nonlinear task, the CD model outperforms the FN one.
Therefore, from the results obtained, we can conclude that the CD model, with the introduction of the new tunable dissipation γ, reaches better performances for all the memory tasks considered.The proposed approach of continuous dissipation allows indeed to achieve better control of the dynamical system response.We also remark that the optimal results for the FN model are in agreement with the ones theorized in [72], where the role of the QRC performance in connection with dynamical phase transitions was discussed.However, for the CD case, the optimal hyperparameters cannot be related to the same physical phenomena.In fact, the CD model does not present any dynamical phases and, as expected, we have found different optimal values of h and ∆t for the tasks considered.In this case, the physical interpretation is related to the interplay between the unitary and dissipative parts of the system evolution.

Time series forecasting
Apart from memory tasks, RC is especially suited for chaotic time series forecasting.A popular benchmark application is the prediction of the well-known Mackey-Glass (MG) dynamical evolution [73].The target obeys the following differential equation: being in the chaotic regime for τ = 17, as reported in previous works on time series forecasting [74,75].
During the training phase, we have injected, as input sequence into the system, the numerical solutions of Eq. (8) sampled with a time resolution t r = 3 (see [76] for more details) and with input values rescaled in the interval [0, 1].The readout weights were trained to solve the one-step-ahead prediction task: During the evaluation phase, the systems have been left to evolve in an autonomous way taking, at each time step, their previous prediction as a current input.
Here, we can distinguish between the shortterm (weather-like) forecast of the dynamical trajectories and the long-term (climate-like) forecast of the chaotic attractor [77].We first show in Fig. 3 (a) that both the CD and FN models are able to reproduce the shape of the chaotic attractor of the Mackey-Glass time series (climate-like forecast).In order to perform a quantitative analysis, we have also tested the capabilities of the systems to predict the oscillations of the target trajectory from the start of the autonomous phase (weatherlike forecast).The metrics used to evaluate the performances for this task was the mean squared error: M where the index i iterates the times and M is the number of points analyzed.Setting M = 150, we were able to see differences as plotted in Fig. 3 (b).The average MSE values for the CD and FN models over these 150 samples are 8 • 10 −3 and 5 • 10 −2 , respectively.This last result allows us to conclude that also for the forecasting task explored here, the introduction of local losses makes it possible to achieve higher performance.

Experimental feasibility
In this section, we discuss relevant aspects of the experimental implementation of the CD model on real hardware.First, we consider the case where access to a number of measurement samples is limited, going beyond ideal conditions, [55,56,57], and address the achievable performance.Then, we will present some potential experimental platforms where the proposed QRC model can be implemented.

Measurement effects
Access to the information contained in the dynamics of a quantum system is generally limited by the stochastic nature of quantum mechanics.In the benchmark tasks presented so far, our analysis has been performed in the ideal case where the readout layer is built up from the exact expectation values of selected observables.In a real experiment, however, only finite measurement samples can be accessed.We will now show how the performance of the two QRC models under consideration (FN and CD) is affected by this fundamental limitation.In particular, denoting by N s the number of samples (i.e. the size of the available ensemble), we will construct the output layers of the models by approximating each observable of interest ⟨ Ô⟩ with the mean of its measurement: where we use O i to indicate a generic value of the measurement.We follow the same approach as in [55]: considering that N s >> 1, we can apply the central limit theorem so that each ⟨ Ô⟩ Ns is randomly drawn from a Gaussian distribution whose mean is the ideal expectation ⟨ Ô⟩, while its maximum standard deviation is σ max (N s ) = 1/ √ N s .The effect of this finite-size statistic is assessed by varying the order of magnitude of N s in the case of the STM task presented above.For each combination of delay value τ and number of samples N s , we have optimized and tested the models in the same way as done in Sec.4.1, keeping the same dataset sizes and exploring the hyperparameter space with the same criteria.From Fig. 4 we can clearly see that the prediction ac- curacy for both models increases with N s as a consequence of the better precision of the observable estimates.We have considered both a short (τ = 2) and a longer (τ = 10) delay.For these two STM tasks, we observe that for a relatively small sample size, the FN model achieves higher capacities than the CD model, although the performance is not optimal for either.We also observe that the FN model can saturate its best performance for smaller samples (e.g.N s ∼ 10 7 for τ = 10) but in general, the CD model can outperform the FN model.Actually, even for a large delay, while the STM of the FN model saturates to rather poor values (for N s ≳ 10 7 ), the CD achieves a very good performance.
It is important to note that the relationship between the number of measurements and the performance found is specific to the model being considered, and can be assessed precisely only in some special cases [56].While the precision of observable expectation values is expected to improve with a general trend of approximately 1/ √ N , predicting how this translates into the performance in different tasks is more complex and cannot be done in general.As expected, the performance for averages obtained considering a sufficiently broad sample saturates to the ideal value, where averages are estimated over an infinite ensemble (represented by the rightmost values in Fig. 3).
In a real experimental implementation, due to the back-action effect of projective measurements after each sample is taken, the reservoir can no longer be used and a new experiment must be started, as analyzed in Ref. [55].The fading memory property, which will be discussed in detail in section B, can significantly reduce the experimental time required.In fact, it is not necessary to inject the whole time series in each experiment, but only a shorter number of points equal to the time width that the reservoir is able to remember.In addition, the experimental time can be further reduced by using a weak measurement protocol [55].The relationship between continuous dissipation models and weak measurement schemes is an issue to be addressed in the future.

Platforms
Various driven-dissipative systems have been experimentally realized in numerous platforms such as cavity QED and cold atoms [78,79,80], circuit QED [81], and trapped ions [82,83,84], due to the rapid progress of experimental techniques in the last 20 years.In addition, they are prominent platforms for quantum simulation and quantum computation [17], especially relevant for quantum computing systems in the NISQ era.
In particular, several experimental platforms capable of implementing dissipative Ising models, which include our CD model, have been proposed in the literature.We mention significant results obtained with Rydberg atoms [85,86,87,88], ions traps [88,89,90,91] and cold atoms [92,93].For all these physical systems, the local independent dissipation present in Eq. 5 can be easily implemented using standard techniques.Following the pioneering work of Ref. [17], the desired dissipation can be achieved by coupling each qubit to an ancillary one that presents spontaneous emission through the well-known strategies of Ref. [94].Moreover, this model of dissipation can be easily implemented by means of optical pumping [95].
Beyond the analog system approach to QRC [96], we also mention that quantum circuits can efficiently simulate generic Markovian dynamics, like that of Eq. 5.In particular, the strategy consists of implementing the collision model algorithm [97] as it has been experimentally demonstrated on the IBM platform [63].

Discussion and Conclusions
In this paper, we have shown that tunable losses in external environments can be turned into a crucial factor for quantum reservoir computing, to tailor and optimize the memory capabilities of the system.In our analysis, we have compared the performances of continuous dissipation (CD) QRC maps with alternative ones based on a discontinuous erase-and-write, as proposed by Fuji and Nakajima (FN) [34].The CD reservoir is based on an input-dependent generator of the dynamics L(s k ) and on the presence of Markovian dissipation modeled through a local master equation in the Lindblad form, with uniform and tunable decay rates.
Through a set of non-linear memory and forecasting tasks commonly employed to benchmark time series processing, we have shown that the degree of tunability of the losses represents a powerful tool that allows the system to be reconfigured with respect to the specific problem under consideration.Indeed, for each task, we have found a range of values of the losses for which the performance indicators exceed those obtained with the FN model, both for memory and temporal prediction tasks.This improvement was also shown to be robust in the realistic case of a finite number of measurement samples.
As a further major result, we have analytically shown that the CD model fulfills the three necessary conditions for time series processing as a reservoir computer (the echo state, fading memory and input separability properties) proving that it forms a class of universality, approximating any fading memory map with arbitrary precision.Finally, our proof has shown that this last result is a more general feature of Markovian dynamics in open quantum systems.In fact, universality is achieved if the generator of the dynamics has the following general, mild properties: for each input injection at sufficiently long ∆t, it must admit only one stationary state; it must be a continuous function of the input; finally, for different inputs, the corresponding stationary states must be in turn different.Some considerations can be added about the generality and interpretation of our results.Looking at the FN and CD CPTP maps, we can say that the FN model of Eq. ( 2) can be described as a sequential map composed by a dissipation (Φ D ) and a Hamiltonian evolution (Φ U ): A slightly different approach was taken in Refs.[44,45] in the attempt to model the decoherent noise in a quantum circuit, where the roles of Φ U and Φ D were reversed: Φ = Φ D •Φ U .The Markovian continuousdissipation map of Eq. ( 6) goes beyond these proposals as dissipation and driving act in a continuous and, in general, non-factorizable way.
Interestingly, our map can approximate all these factorized ones [34,44,45] by a proper choice of the dissipation rates and their time dependence.Indeed, for the erase and write map of Eq. ( 2) this is shown in Appendix A. On the other hand, the model used in Refs.[44,45] assumes a noise modeled by a Markovian channel Φ D so that also Φ is Markovian, and in general can be written in the GKLS form (3). In conclusion, we remark that the major generality and tunability obtained with our formalism has been revealed to be useful in practical applications as shown in the performance improvements found in Sec. 4.
From a more physical point of view, we can relate the substantial improvements found in the CD model with the form of dissipation and driving (continuous information erasure and injection) in the reservoir.Indeed, an essential ingredient for RC to work is to find a balance between information spreading across the reservoir and dissipation, which is essential for the fading memory property.In the FN model, the input is injected locally and needs to spread during the evolution [72], while dissipation is operated by the partial trace and not optimized.On the contrary, in the CD case, information is already injected over the entire reservoir, and at the same time, the degree of dissipation can be fine-tuned by the parameter γ.These two limiting features of the FN map are at the origin of the higher performance of the CD one, even varying the dimensions as shown in Appendix E.
While a clear improvement is already obtained here assuming the simplest form of dissipation, the idea of dissipation engineering can be explored also considering non-local losses [24], as well as non-Markovian dissipation [98], opening the way to study a wider range of quantum reservoir computers.In conclusion, in this work, we have established a new paradigm of QRC based on Markovian dynamics that is suitable for generalization to more complex dissipation engineering techniques. is also acknowledged.GLG is funded by the Spanish Ministerio de Educación y Formación Profesional/Ministerio de Universidades and co-funded by the University of the Balearic Islands through the Beatriz Galindo program (BG20/00085).AS acknowledges funding from the CSIC hub on AI through the scholarship JAEIntroAIHUB2-19.The project that gave rise to these results received the support of a fellowship from the "la Caixa" Foundation (ID 100010434).The fellowship code is LCF/BQ/DI23/11990081. RMP acknowledges the QCDI project funded by the Spanish Government and part of this work was also funded by MICINN/AEI/FEDER and the University of the Balearic Islands through a pre-doctoral fellowship (Grant No. MDM-2017-0711-18-1).

A Fuji-Nakajima map and Lindblad dynamics
In the following, we want to frame the FN model into our Markovian strategy to obtain dissipation by designing an equivalent input injection strategy.The idea is to replace the role played by the partial trace by introducing a strong loss on the first qubit and then preparing its state by applying the rotation operator to it.More precisely the network's damping rates will be and the input s k will be converted into an angle of rotation θ k according to θ k = arccos √ 1 − s k .Accounting for the two steps, we can determine the conditions for the entire process, at each input injection, to be equivalent to the updating step of Eq. (2), i.e., R (1)  y (θ k+1 )e

B Proof of universality
In this section, we will prove that the CD model proposed in Eq. (6) forms a universality class for reservoir computing.In the context of QRC, universality corresponds to an approximation property and, in particular, is referred to the capability to arbitrarily well approximate any fading memory map [27,68,43,41].It will be shown that our proposal satisfies both the echo state property (ESP) and the fading memory property (FMP).Finally, we will prove separability verifying that the hypotheses of the Stone-Weierstrass theorem are fulfilled, then proving the universality.

B.1 Echo state property
The ESP refers to the capability of a reservoir to forget the initial conditions in the limit of an infinite input sequence [58].It has been recently shown that this is equivalent to the usual definition of a unique solution in the RC literature, under some mild conditions (Theorem 1 in [99]).In our case of study, it means that given a left infinite and time-ordered input sequence {• • • , s −1 , s 0 } (where the 0-th time step corresponds to the last input) and given two different initial quantum states ρ 1 and, ρ 2 , whose dynamics is driven by this sequence, the evolved states will become indistinguishable: where the limit is considered to be pointwise and ∥ • ∥ indicates any matrix norm (since all norms are equivalent in finite-dimensional spaces).
A sufficient condition for this relation to hold is strict contractivity (Theorem 2.2 in Ref. [100]): for all s k and any pair of states ρ 1 and ρ 2 , where 0 ≤ r < 1.In particular, this can be ensured for Markovian maps that provide a unique stationary state for a generic input s k , a property that our model possesses according to Theorem 2 [101], setting a long enough value of ∆t.More precisely, for a generic generator L(s k ), we can always define a finite mixing time, denoted as ∆t mix , such that if ∆t ≥ ∆t mix then the map is strictly contractive and the ESP holds.
As a further observation, we estimated the ∆t mix trend as the system dimensions increase.Our numerical analysis suggests that the ∆t mix upper bound has an asymptotic behavior ∼ N γ (see Appendix C for more details).Concluding, the mixing time of all the proposed maps e L(s k )∆t in addition to being finite, according to what we have observed numerically, scales efficiently with the number of qubits.

B.2 Fading memory
Considering the set of left infinite sequences of inputs belonging to the interval [0, 1], K − ([0, 1]), the fading memory property is a condition of continuity of functionals defined on this set with a given norm.For a sequence s ∈ K − ([0, 1]), this norm is defined starting from a generic null sequence, say w = {w k } k≥0 such that lim k→∞ w k = 0, in the following way: where s k and w −k are respectively elements of s and w and Z − is the set of negative integers.By definition, we say that a map has fading memory if, for any null sequence w, it is continuous in (K − ([0, 1]), ∥ • ∥ w ).In order to prove that the model of Eq. (6) has fading memory, it is sufficient to prove that e L(s k )∆t , with ∆t > ∆t mix to ensure the ESP, is a continuous function of s k according to [36] (Lemma 3).If L(s k ) is continuous then the continuity of its exponential will be a direct consequence.Let s k ,u k ∈ [0, 1], ∥ • ∥ 2 be the Hilbert-Schmidt norm and ρ be a unit density matrix belonging to C 2 N ×2 N .Then where the first inequality is a consequence of the following property of the commutator: which can be proven expanding ρ as a linear combination of Pauli strings and using their commutation rules.The last term of Eq. ( 11) implies continuity because it shows that an arbitrarily small distance between the two functions, let us say ϵ, is reached when the distance of the two inputs is 2N .We then conclude that the family, say L, of functionals defined by Eq. (6), with the output layer of expected values, belongs to the set of fading memory functionals C(K − ([0, 1]), ∥ • ∥ w ), which map input sequences into real expected values of observables.

B.3 Stone-Weierstrass theorem
The universality condition is obtained if L is also a dense subset of C(K − ([0, 1]), ∥ • ∥ w ).It implies that any fading memory map f ∈ C(K − ([0, 1])) can be arbitrarily well approximated by an element of L. Quantitatively for any f ∈ C(K − ([0, 1])) and for any ϵ > 0 exist l ∈ L such that the following condition holds: A sufficient condition to obtain this result is given by the well-known Stone-Weierstrass theorem: Let E be a compact metric space and C(E) be the set of real-valued continuous functions defined on E. If a subalgebra A of C(E) contains the constant functions and separates points of E, then A is dense in C(E).

B.4 Input separability
Now we prove the separability of L. It means that given two sequences s, s ∈ K − ([0, 1]) such that s ̸ = s exists at least one element of L able to separate them.In our context it is translated to the condition that the corresponding density matrices at the 0-th time step must be different: ∥ρ 0 − ρ0 ∥ ̸ = 0.It is useful at this point to show an important property of Eq. (6): Lemma 1: Given two different inputs s k , u k then the corresponding stationary states of L(s k ) and L(u k ) are different. Proof: We firstly give a necessary condition for a stationary state ρ ss which corresponds to a generator L(s k ).It is helpful to expand ρ ss in the basis of Pauli strings: where the indexes {i n } label single Pauli matrices: i n = x, y, z or the identity matrix, say σ 0 = I, for the Hilbert space of a single qubit.A necessary condition can be found projecting the definition of ρ ss (i.e.L(s k )ρ ss = 0) on a given Pauli string; we will consider σ z i : Observing that we can write the Lindbladian of Eq. (6) as a sum of single qubit dissipators: its action on the single qubit Hilbert space: implies the following expression for Eq. ( 12): where α z i = Tr{σ z i ρ ss }, α y,x i,j = Tr σ y i σ x j ρ ss , α y i = Tr{σ y i ρ ss }.Considering now another generic input u k , such that u k ̸ = s k , if ρ ss is a common stationary state of both, the following relation must be satisfied: Equation (14) gives necessary conditions for the coefficients of ρ ss and, among these, it gives these relations: Tr{σ a i ρ ss } = α a i = 0, Tr σ a i σ x j ρ ss = α ax ij = 0 with a = y, z and 1 ≤ i, j ≤ N .Then Eq. (13) becomes γ 2 N = 0 which is not satisfied because we are assuming that γ ̸ = 0 in order to guarantee the ESP.As a consequence, we can conclude that Lemma 1 holds.■ Another useful statement to arrive at the input separability is the following: Lemma 2: Let s k and ρ ss (s k ) be a generic input and its corresponding unique stationary state, let P be the set of the density matrices of the considered N qubits system and let d : P × P → R + be the distance induced by the Hilbert-Schmidt product, then the following function g : R + × P → R + , (t, ρ) → d e L(s k )t ρ, ρ ss (s k ) is bounded, its maximum d (s k ) max always exists and it is strictly decreasing in t when t > ∆t mix .

Proof:
We first observe that for any fixed value of the time g is a continuous function on P.This is because it is a composition of continuous functions: the action of e L(s k )t composed to the action of the distance from a fixed point.As a consequence because P is a compact set and because the continuous image of a compact set is a compact set we arrive at the following relation: We observe that for the boundaries of the set of the density matrices d Considering now the action of g for a fixed density matrix in P, we know that g(t)| ρ = d(e L(s k )t ρ, ρ ss (s k )) is strictly decreasing for all ρ ̸ = ρ ss (s k ) if t > ∆t mix as a consequence of the already mentioned strict contractivity of e L(s k )t .Finally the same property must hold for d We can now return to the problem of the separability in which we have the two generic different sequences s and s.Considering the smallest index J such as s −J ̸ = s−J , applying Lemma 1 we know that two different stationary states exist: ρ ss (s −J ) and ρ ss (s −J ).We define two open balls of radius r respectively centered on the two states: B r (ρ ss (s −J )) and B r (ρ ss (s −J )) such that B r (ρ ss (s −J )) ∩ B r (ρ ss (s −J )) = ∅.The last condition is easily satisfied if r < d(ρ ss (s −J ), ρ ss (s −J ))/2.In order to find an element of L which separates the sequences, applying Lemma 2, we have to set ∆t into a sufficiently long value so that: d With this condition, it is ensured that after the applications of the two inputs, the resulting states of the reservoirs will be different regardless of their states at the time-step −J.After this injection by hypothesis, the subsequent inputs will be the same, and necessary the corresponding states of the reservoir at the time 0 will be different because e L(s k )t is a full rank linear operator.We can in this way conclude that the input separability is satisfied.

B.5 Polynomial algebra and final considerations about the universality
Another hypothesis that must be satisfied in order to have universality is the presence of a subalgebra in L. Since we are working on a sys-tem that gives as output a linear combination of observables of the density matrix we have to look for a polynomial algebra.It can be obtained by adding to the family the model of Eq. (6) in spatial multiplexing.It means that we can consider V different and independent states {ρ (1) , • • • , ρ (V ) } whose dynamics will be lead by generators of the form of Eq. (6) with, in general, different value of the characteristic hyperparameters: { L1 , • • • , LV }.We can write the total updating rule of the total reservoir in the following way: Considering the set of the polynomial outputs for all the single states {P 1 , • • • , P V } we achieve the algebra considering as output for the reservoir a linear combination of the polynomials: These newly added reservoirs satisfy the fading memory condition according to [36] (Lemma 5).
The only condition that remains to be proven in order to assert the Stone-Weierstrass theorem is the presence of constant functions in the family L but it is trivially satisfied due to the fact that we are working with polynomial inputs.We can now conclude that L is a universal class for reservoir computing.
Finally, we notice that this proof does not hold only for Eq.(6) but for a more general class of reservoirs working with a master equation.The conditions required for the generator L are the following: (i) it admits only one stationary state for each input s k , (ii) it must be a continuous function of s k , (iii) given two different generic inputs, say s k and u k , the stationary states of L(s k ) and L(u k ) are also different.

C Mixing time scaling
In Sec.B.1, we have defined the mixing time ∆t mix as the minimum interval of time ∆t required of e L(s k )∆t to fulfill the Echo State Property.From a computational point of view, we are interested in estimating how ∆t mix scales with the system size.We will now show the strategy used to achieve it through numerical simulations.We, firstly, recall that e L(s k )∆t can always be expanded in terms of its dual basis except for a countable number of points in the hyperparameter space that have zero probability to occur [103].Let {λ i } be the set of eigenvalues of L(s k ) and let {|r i ⟩⟩} and {|l i ⟩⟩} be the corresponding set of orthogonal and normalized right and left eigenvectors respect the Hilbert-Schmidt product.As a consequence of the uniqueness of the stationary state of L(s k ) (Theorem 2 of Ref. [101]), when its action is restricted to the set of traceless Hermitian matrices, which we will denote as H 0 , we can conclude that the real part of all its eigenvalues will be strictly less than zero.For the sake of definition, we will order the eigenvalues such that Re{λ 0 } = 0 > Re{λ We have numerically computed a maximal value of τ spanning the system dimension from 3 to 6 spins.For each case, the magnetic field h took the values from the set {0.01, 0.05, 0.1, 0.5, 1,5, 10}, the damping rate γ from {0.01, 0.1, 1, 10} while the input has been fixed to be s k = 0.For all the possible combinations we have simulated 100 realizations of the Hamiltonian couplings {J i,j } and we have taken the average value of τ as a representative.In our analysis, we have selected the max average for each N, which computationally corresponds to the worst case, in order to estimate an upper bound.As shown in Fig. 6, our numerics suggest that the max order of magnitude of ∆t mix scales proportionally to N γ .As a result, we have numerically observed that the minimum value of ∆t required to reach the ESP, for the map of our model generated by L(s k ), scales in an efficient way with respect to the number of qubits.

D Training details
As explained in Sec. 5, we benchmarked the considered quantum reservoirs by computing their performance for different network realizations.In the case of memory tasks, we also varied the random inputs injected.It is important to note that training a reservoir for a specific task depends on the reservoir evolution rule.Therefore, we performed independent training for each different network coupling sampled.
Going more into details, we indicate with {⟨ Ô⟩ k,j } L,M k=1,j=1 the set of the expectation values of the reservoir observables involved in the training phase, where L is the number of training points and M the number of chosen observables.Considering the free weights that determine the reservoir output at each time {w j } M j=1 , to be determined by the training, and set of targets which defines the task of interest {y k } L k=1 , the implemented training consisted on minimizing the following euclidean distance: We have numerically performed this optimization by using the LAPACK library [104], whose strategy consists of making use of the singular value decomposition.

E Performances for different reservoir sizes
While in the main text we considered five qubit reservoirs, here we explore the effect of increasing or decreasing the size of the reservoir.Following the same optimization procedure described in Sec. 4 regarding the choice of optimal hyperparameters, in Fig. 7 we show the optimal STM task performances of the FN model and the CD model and for a delay τ = 10.
The dimension of the output layer changes with the system size.Thus, as expected, the performance of both the FN and CD models improves with the number of qubits (more degrees of freedom are exploited for solving the task).Interestingly, the CD model outperforms the FN model for all considered values of N .For the task at hand, it is reasonable to expect that for a sufficiently large reservoir, both QRCs are able to achieve optimal performance.However, we see that the CD model requires a smaller number of qubits to saturate.In fact, N = 6 qubits are already sufficient to reach a capacity value C ≃ 1, while the FN would require N > 7, implying the need for more resources to achieve the same performance.

Figure 1 :
Figure 1: Schematic of the CD model proposed in Eq. (6): spin network (with Hamiltonian H ′ given in Eq. (4)) with uniform losses at rate γ governed by D L .The discrete input (bottom right) is injected uniformly into the network nodes through a time-dependent magnetic field (bottom left).The output layer is constructed by measuring a portion of the observables of the density matrix of the reservoir.

Figure 2 :
Figure 2: Performance evaluation of different memory tasks for the FN model described in Eq. (2) (red) and for CD in Eq. (6) (blue).Capacity for STM task (a), NARMA task (b), and parity check (c).Optimized hyperparameters as specified in the main text for the STM task, while for NARMA and parity check the hyperparameters are also optimized but not reported, being not relevant for the conclusions.A number of 1000 time steps have always been used for the washout, training, and evaluation phases.In all the plots the shadow regions cover one standard deviation taken as statistical error over the 100 random realizations of spin coupling values for the Hamiltonians of the systems {J ij }.

Figure 3 :
Figure 3: Results of the Mackey-Glass time-series prediction of the systems.In all the cases we have set 1000 points for the washout and training phase.The optimal free hyperparameters found, and fixed in the following plots, are h = 1, ∆t = 10 for FN model and h = 0.1, ∆t = 0.1, γ = 10 for the CD model.(a) Chaotic attractor autonomous reproductions (CD and FN) in the phase space y(t)-y(t − 6) compared to the target attractor (MG) for one exemplary realization of the systems couplings {J ij }. 2000 points have been used for the plots.(b) Predictions of 150 points in the autonomous phase contrasted to target values of the series averaged over 100 random realizations of the coupling values.The shadows cover one standard deviation taken as a statistical error.

Figure 4 :
Figure 4: STM task performance as a function of the number of samples N S .Two delay values are considered, τ = 2 and τ = 10, while varying the magnitude of N S from 10 4 to 10 12 .The right part of the plot shows the ideal performance.

Figure 5 :
Figure 5: Hilbert-Schmidt distance of states evolving according to the left-hand side or the right-hand side of Eq. (9), say ρ A and ρ B , as a function of the amount of dissipation.Values are averaged over 100 random pairs of density matrices and inputs, respectively ρ k and s k in the formula.

Figure 6 :
Figure 6: Max average values of τ in γ units as a function of the number of spins.The purple points indicate the averages while the purple vertical lines cover one standard deviation taken as statistical error over the {J i,j } realizations.The linear fit result is represented by the red line.

Figure 7 :
Figure 7: STM task optimal performances of the FN model and the CD one as a function of the number of qubits N , fixing a delay τ = 10.The optimal values have been computed in the same way as in Sec. 4. We averaged over 100 network realizations and input series for the cases of N ≤ 5 while we took a statistics of 10 realizations for the complementary ones.The error bars refer to one standard deviation taken as a statistical error.