Analytical evidence of nonlinearity in qubits and continuous-variable quantum reservoir computing

The natural dynamics of complex networks can be harnessed for information processing purposes. A paradigmatic example are artificial neural networks used for machine learning. In this context, quantum reservoir computing (QRC) constitutes a natural extension of the use of classical recurrent neural networks using quantum resources for temporal information processing. Here, we explore the fundamental properties of QRC systems based on qubits and continuous variables. We provide analytical results that illustrate how nonlinearity enters the input–output map in these QRC implementations. We find that the input encoding through state initialization can serve to control the type of nonlinearity as well as the dependence on the history of the input sequences to be processed.

Recently, also the capability of quantum networks to process information has begun to be explored. By combining the properties of neural networks with the unique features arising in the quantum realm, quantum neural networks are expected to offer several advantages with respect to their classical counterparts, such as higher effective dimension, exponentially increased memory capacity, and faster training and learning [38,39]. In this context, a first proposal for extending reservoir computing (RC) to the quantum domain has also been put forward recently, based on a network of qubits [40]. RC is a three-layer (recurrent) neural network singularly suited to solve temporal tasks [41]. Several realizations of classical RC [41][42][43] were already established in the last years in photonics, spintronics, mechanics and biological systems [44][45][46][47][48][49][50][51][52][53].
A good performance in RC is known to be achieved by exploiting a high-dimensional physical system, internal memory, and nonlinearity [41,54]. As for the system size, large reservoir networks can be considered in classical systems or, as a promising alternative, in quantum ones. In fact for quantum networks, the large Hilbert space displayed even with a reduced number of nodes is one of the main motivations to extend RC into the quantum regime [40,[55][56][57]. As for the internal memory, it refers to the ability of the reservoir to store the information of previous inputs up to sometime in the past and it is generally related to the extension and dynamics of the reservoir computer. As for the nonlinearity, this is a property needed in the full input-output transformation and it depends on the specific algorithm definition as well as on the dynamical properties of the reservoir. Regarding physical implementations, most classical realizations of RC are based on nonlinear systems [44][45][46] with a few exceptions introducing nonlinearity through input encoding or detection [58][59][60][61][62]. In quantum systems, despite the fact that the state dynamics are linear, efficient RC has been reported and even in the absence of any nonlinearity in the reservoir [56]. Indeed quantum bosonic networks restricted to Gaussian states, where also observables' dynamics are governed by a fully linear (and finite) set of equations, are a promising platform for quantum RC. The origin of the nonlinearity in quantum reservoir computing (QRC) in this linear system can actually be traced back to the input encoding as first shown in [56] and recently reviewed in [63]. In fact, nonlinear encoding has been emphasized as an important tool in the field of quantum machine learning, outsourcing data processing and gaining power [64][65][66], as well as providing universal approximation properties in different settings [67,68].
In this work, we specifically address how nonlinearity enters the input-output map in quantum RC implementations both with spin and bosonic networks, providing insightful analytical results that also allow us to establish new nonlinear processing strategies through input encoding.

Classical and quantum RC
We start by defining RC, a machine learning technique belonging to the family of recurrent neural networks that exploits the natural dynamics of input-driven complex systems for information processing [41]. This algorithm can be divided into three steps associated with the relative system layers: (a) inject an input into the dynamical system, here a complex quantum network; (b) let the dynamical system (also known as a reservoir) evolve under its natural dynamics; (c) extract information from the reservoir, using all or some of its degrees of freedom via an output layer.
The main advantage of RC in the machine learning context with respect to other neural-network algorithms is its minimal requirements for training [69,70]. Indeed, only the weights of the output connections are optimized toward the target function via e.g. linear regression, while the reservoir layer does not need any fine tuning of parameters. RC is especially suited for temporal data processing, as the input temporal dependence can be rooted in the complex dynamical behavior of the reservoir. Despite the simplicity, classical RC has been successful in numerous practical applications as diverse as the recognition of human actions in video streams, arrhythmia detection in human heartbeats or chaotic time-series forecasting [71][72][73].
QRC is the natural extension of the classical setting toward implementing quantum dynamical systems as reservoirs. It has already been successfully demonstrated for time-series processing [40,74] and quantum state measurement [75], and the processing capacity has been quantified in different implementations [56,76]. The main advantage of quantum systems is that they exhibit a large number of degrees of freedom that can be exploited for information processing.
Since the first proposal of networks of quantum spins as reservoir substrates [40], there have been a variety of works exploring the possibilities that the spin-based implementation can offer [40,55,74,[76][77][78][79][80][81][82][83]. These spin-based systems stand as the most analyzed platform for QRC at the moment, but also other finitedimensional systems have been considered, like fermionic networks [84,85], and infinite-dimensional ones, like bosonic reservoirs [56,86,87], as recently reviewed in reference [88]. In all these applications, computational performance is enhanced from having a nonlinearity in the input-output map. In the following sections, we focus on this aspect considering both qubits (section 2) and bosonic continuous variables (section 3) networks and we provide an analytical treatment.

Nonlinearity of the input-output map QRC with qubits
In the present section, we deal with a quantum reservoir consisting of a complex network of N randomly coupled qubits in a finite-dimensional Hilbert space of dimension 2 N . This quantum reservoir computer was introduced in reference [40] and its potential was further explored in references [55,89], while the performance in terms of its dynamical phases was established in reference [76]. Figure 1 presents a schematic representation of this spin-based QRC approach. We first examine the presence of nonlinear dependence on the input-output QRC map and derive an analytical form of the observables as a function of the state of the reservoir. We focus on the factors that influence the precise shape of the nonlinearities like the form of input encoding, the input history or the dynamical regime. To this end, we derive a general expression for the system observables, which is independent of the particular Hamiltonian that describes the unitary evolution of the system between input injections. In order to exemplify the effect of working in different dynamical regimes, we numerically simulate Figure 1. Schematic representation of a RC system with a network of quantum spins. The input information is injected into the system through the state of one or more of the spins. Then, the system evolves under the Hamiltonian dynamics and finally some of the observables are read out.
the system of reference [76], which presents a dynamical phase transition from a many-body localized (MBL) to an ergodic phase.
In the following, for illustration purposes and for simplicity, we look at the case of a discrete scalar input that is given to the system at each time step k, so that the input is a series of random real numbers, {s k }, produced from a uniform distribution in the interval [0, 1]. In a first approach, we treat the case of the input being encoded into a single qubit and then we move to other configurations.

Input encoding in one qubit
As a starting point, we consider the scenario in which the input is introduced into the reservoir system through the reinitialization of the state of a single qubit for different types of states.
In the first place, we analyze a widely used input encoding in QRC as the single-qubit pure state [40,55,76,77] |ψ with a normalized input, such that s k ∈ [0, 1] (the effect of a phase factor will be addressed in section 2.1.2). As remarked in reference [63], this input choice implies a linear encoding in the z direction, ψ k |σ z |ψ k = 1 − 2s k , a nonlinear one in the x direction ψ k |σ x |ψ k = √ s k (1 − s k ), and that no information is codified in the y direction, ψ k |σ y |ψ k = 0. Below, we show how this dependence is transferred to the output observables of the whole reservoir system, which are eventually used to obtain the output when performing a given task [88].
In terms of a density matrix, the input state is written as: where r k = √ s k (1 − s k ). We arbitrarily chose the first qubit to feed the input to the reservoir network in an instantaneous way between steps of unitary evolution of duration Δt under the HamiltonianĤ, which governs the coherent evolution of the qubit network. Therefore, the map that describes this process for the state of the whole system at each time step k is the following: whereρ k−1 = Tr 1 (ρ k−1 ) can be represented as a matrix of dimension 2 N−1 × 2 N−1 and Tr 1 is the partial trace performed with respect to the first qubit. The term inside the parenthesis in equation (3) describes the state of the system after the input injection and before the time evolution that can be decomposed as where 0 is the zero matrix. The state at time step k displays three contributions: (i) a s k -independent term; (ii) a term proportional to s k ; and (iii) a nonlinear term in s k through r k = √ s k (1 − s k ). After evolving during Δt, the state of the system takes the following form: where ρ (0) k , ρ (1) k , and ρ (nl) k depend on the previous state ρ k−1 , the particular form of the HamiltonianĤ and the values of its internal parameters, and the time interval Δt. Besides that, from equation (5), we can immediately write the functional form of the expected value of any observableÔ as The insight gained from equation (6) is that the functional form that relates the current input and all possible output expectation values is the same, irrespective of the choice of observables. For instance, singlequbit observables or pairwise correlations similarly follow this functional relationship. Furthermore, in the full non-unitary evolution of the state (see equation (5)), we can identify the effects of both the input encoding (determining the specific functional dependence on s k ) and Hamiltonian choices. More specifically, for the observables, in equation (6), the relative importance between the three terms contributing to the sum can vary from a linear trend in s k to have a dominant contribution coming from the last nonlinear term. This will depend on the observable in consideration and, more importantly, on the particular structure of the reservoir and on the dynamical evolution of the system, both established by the Hamiltonian. Remarkably, this fact brings out the importance of the reservoir's design and working regime. To illustrate this, in the following lines, we concentrate on the system Hamiltonian (and regimes) studied in reference [76], which has demonstrated a high performance in realizing nonlinear demanding tasks as the NARMA task.
The physical system receiving attention is made of N qubits (or spins) that will evolve in between input injections under the unitary evolution given by the transverse-field Ising Hamiltonian: The spin-spin couplings J ij are randomly generated from a uniform distribution in the interval The functional dependence on s k in equation (6) can be numerically illustrated in the following way. First, we evolve the dynamical system up to a time step k − 1. Then, we compute the expected value of observables at time step k for a sampling set of the input s k in the interval [0, 1], in such a way that we can plot Ô k versus s k . Note that this is done for a unique input history {s k−1 , s k−2 , . . .} such that the terms ρ (0) k , ρ (1) k , and ρ (nl) k remain constant for all the considered s k values. In other words, by fixing the history up to the time step k − 1, we can sample the input s k and identify the observables' dependence on it.
In figure 2, we show the response of two arbitrarily selected single-qubit observables. The first row corresponds to the expected value of the spin projections in the x-axis ( σ x 1 k and σ x 4 k in (a) and (b), respectively) and the second row contains the spin projections in the z-axis ( σ z 1 k and σ z 4 k in (c) and (d), respectively). We recall that the input is introduced through qubit 1. In all panels, blue dots represent the response of the observables in the ergodic regime, while red dots correspond to the MBL phase. At a first sight, we can distinguish that observables related to the x component have a nonlinear shape (produced by the combination of s k and r k ) and the observables related with the z direction show a linear dependence. In particular, σ x 1 k and σ z 1 k in the MBL phase (red dots) of panels (a) and (c) show the same function as the input encoding itself: σ z 1 k = 1 − 2s k and σ x 1 k = r k , while as expected figures 2(b) and (d) show a negligible response of the observables (red dots) for a qubit different from the input one. Only in the thermal phase a clear (linear or nonlinear) response arises (blue dots).
As we mentioned earlier, the response is determined by the terms ρ (1) k and ρ (nl) k , which depend on the Hamiltonian of the system. In the quantum ergodic regime, information spreads through the reservoir network and all local observables can respond to input injection [76]. However, dynamical regimes of localization prevent a good response of all the observables to the input, due to the presence of an extensive number of local conserved quantities that hinder the spreading of information from the first qubit to the rest of the network.
The functional behavior presented in figure 2 is not exclusive of the single-qubit observables. ). Again, observables in the ergodic phase show a dependence on the input while the MBL phase suppresses it since the input information remains localized in the first qubit.

A mixed state as a linear input encoding
We now explore an alternative input encoding in order to see if the nonlinear dependence on the input is always present. In particular, we choose a mixed state (equivalent to set r k = 0 in equation (2)) as in references  [ 63, 74, 78-80, 89, 90] so that with the information encoded only in the z direction, ψ k |σ x |ψ k = ψ k |σ y |ψ k = 0. Then, equation (6) is simplified to and all the observables indeed remain linear in s k .

Input encoding in all directions
When the encoding of the input is done not only in the relative amplitude but also in the relative phase between basis states with φ k ∈ [0, 2π], there is information in all directions. In this case, we deal with the superposition The dependence is linear in the z direction as before ψ k |σ z |ψ k = 1 − 2s k , and nonlinear contributions are found in the other two directions with additional sinusoidal factors [63] of the relative phase, ψ k |σ x |ψ k = 2r k cos(φ k ), and ψ k |σ y |ψ k = 2r k sin(φ k ).
For example, if the input is a single scalar, as considered before, we could reintroduce it into the system through the phase as φ k = 2πs k . Now, the state of the input qubit reads: and equation (5) is rewritten for the reservoir state as: From this last expression, it follows that the general form of the expectation values of the reservoirs' observables has explicit contributions with sinusoidal nonlinear factors: We notice that this further nonlinearity arises because the input is encoded also in the phase, producing an input dependent rotation.

Dependence of the observables on the input history
The expression in equation (5), where the state of the reservoir is decomposed in three contributions, is valid for any time step k. Thus, we can replace ρ k−1 insideρ k−1 = Tr 1 (ρ k−1 ) in equation (3) to relate the previous input, s k−1 , to the state ρ k : where For any observable, from equation (13), we derive the following expression showing the dependencies on two consecutive inputs: Equations (6) and (14) both imply that the particular input history at the time step k when we evaluate the observables is relevant for the nonlinear input dependence. In addition, the functional form of the observables at a given time is restricted to combinations of s k and r k with themselves at previous time steps (s k−1 or r k−1 ).
As an example, for s k−1 = s k in equation (14), we see the effect of reintroducing the same input [91,92] in two consecutive time steps on the functional form of the observables. This would increase the number of different nonlinear terms compared to equation (6), where the same input is only introduced once. In return, this would involve using more system memory for processing the same input.
Just to give a taste of these implications, we represent in figure 4 the response of a local observable at two different time steps. Panel (a) shows the dependence of the observable σ x 4 on the input at time steps k (red dots) and k − 1 (blue dots). The lines are generated with the same input history {s 1 , s 2 , . . .} up to k − 1, where we sample the interval [0, 1]. In the case of time step k (red dots), s k−1 is fixed at an arbitrary random value. We can clearly see in figure 4(a) that the precise shape of the input dependence of σ x 4 is determined by the input history, even with only one time step of difference. In figure 4(b), we represent all the possible dependencies of this observable respect to inputs s k and s k−1 (for the selected parameters). The axes x and y contain the values of the input while the z axis displays the magnitude of σ x 4 (the color of the surface goes from the maximum value in yellow to the lowest value in dark blue), illustrating the rich complex nonlinear response of the observable for different s k and s k−1 combinations and reflecting the nontrivial capabilities of spin-based quantum reservoir computers [55].

Beyond one-qubit encoding
Among other possible extensions of the input encoding, here we examine the effect of encoding the same input in several (N E ) qubits, assuming the possibility of accessing a subset of the qubits in order to reinitialize their state recurrently. With respect to the previously considered one-qubit encoding, the information is more evenly and redundantly distributed into the system. Firstly, we consider the encoding in a product of identical states of the form: Within this framework, the input-output mapping is rewritten as where in the present caseρ k−1 = Tr 1,...,N E (ρ k−1 ) corresponds to a matrix of dimension 2 N−N E × 2 N−N E . Following the same procedure to decompose the state as before in the one-qubit case, remembering that r 2n k = [s k (1 − s k )] n , each of the matrix elements of the state (15) is a linear combination of powers of s k multiplied either by r k or by 1. Then, the state ρ k can be formally written as where the explicit form of the matrices ρ (N E ) , ρ (j,0) and ρ (j,1) is omitted. Thus, the expectation value of any observableÔ is given by: As we have seen in section 2.1.1, in the case of mixed-state encoding over a qubit, the nonlinear input dependence is lost. Interestingly, extending such an encoding to multiple qubits, the nonlinearity is immediately restored. Indeed when the input is introduced via multiple qubits and we use a product of identical mixed states (1 − s k )|0 0| + s k |1 1|, the general form of the expectation value of the observables is a polynomial whose degree is determined by the number of qubits used to codify the input, N E :

Nonlinearity of the input-output map QRC with continuous variables
Besides qubits or other finite-dimensional systems, RC can also be implemented by harnessing infinitedimensional quantum systems such as the bosonic harmonic oscillator. We will focus on the framework introduced in [56], where a complex network of interacting oscillators is considered and focusing on Gaussian states. A small subset of them-the ancillae-are subject to periodic state resets conditioned on the input whereas the rest play the role of the reservoir. The output is some trained function of the observables of the reservoir oscillators, while the effect of measurements on the dynamics is disregarded as in the previous section. Let us consider N identical unit mass oscillators with frequencies ω 0 and work in such units that = 1 and k B = 1. Letx = {q 1 ,p 1 ,q 2 ,p 2 , . . . ,q N ,p N } be the 2N vector of position and momentum operators of the network oscillators. Then the covariance matrix σ(x) and first moments vector x are defined as [56] where i, j = 1, . . . , 2N. Throughout this section, we indicate the elements of the matrix and of the vector with superindices (i) and (j) whereas the subindex k is reserved for time steps.
In the case of a single ancilla in a Gaussian state, the state is completely determined by [93,94] σ(x A ) = n th + 1 2 (y + z c )/ω 0 z s z s (y − z c )ω 0 , x A = |α| cos(arg(α)) 2/ω 0 |α| sin(arg(α)) √ 2ω 0 , (22) wherex A = {q A ,p A } is the vector of positionq A and momentump A operators of the ancilla, whereas σ(x A ) and x A are the corresponding covariance matrix and first moments vector, respectively. The covariance matrix depends on the amount of thermal excitations n th as well as the magnitude r and phase ϕ of squeezing [95] which have been collected into terms y = cosh(2r), z c = cos(ϕ) sinh(2r), and z s = sin(ϕ) sinh(2r). The first moments vector in turn is a function of the magnitude |α| and phase arg(α) of displacement. Although equation (22) describe a displaced squeezed thermal state, an arbitrary single mode Gaussian state can be recovered by a judicious choice of the parameters n th , r, ϕ, |α| and arg(α). These are also the parameters that can be separately conditioned on the input s k ; this is a quite natural choice as thermal excitations, squeezing and displacement correspond to different state preparation methods [93,94]. As before, the input is a real uniformly-distributed random number, s k ∈ [0, 1]. Let the time between state resets be Δt and the network HamiltonianĤ be quadratic in the elements ofx, with the form:Ĥ where the momentum and position of the oscillators are, respectively, in the vectorsp T = {p 1 , . . . ,p N }, andq T = {q 1 , . . . ,q N }; Δ ω = ω 0 ; and L is a Laplacian matrix with components L (ij) = δ (ij) l g (il) − (1 − δ (ij) )g (ij) determined by spring-like links of strength g (ij) . Then, given the current network state, the next state is completely determined through the action of a 2N × 2N symplectic matrix S = S(Ĥ, Δt) [56], which is a real matrix that preserves the canonical commutation relations. The form of the vector of operators at some time step k is indicated byx k . In practicex k =x R k ⊕x A k+1 , wherex R k are the reservoir operators at time step k andx A k+1 are the ancilla operators conditioned on the input s k+1 . Then, the operators at the following time step have the formx k+1 = Sx k . Consequently, it follows from the definitions given by equation (21) To better understand how the reservoir oscillators evolve, we divide the symplectic matrix S to blocks as where A is 2(N − 1) × 2(N − 1) and D is 2 × 2. We can assume a single ancilla to be the Nth oscillator without a loss of generality. Then the update map reads which governs how the next form of reservoir covariances and first moments depend on the previous form and the ancilla state conditioned on the input. Provided that the reservoir state is also Gaussian, all possible reservoir observables are either elements of the covariance matrix or the first moments vector or a function of these elements. By extension, the output is also a function of the elements of these quantities.
One may now ask how the reservoir observables depend on the input s k . To this end, one may solve the recurrence relations of equation (25) as in [56]. Here we present the solution for the first moments vector; the solution for the covariance matrix is similar. Given the initial state x R 0 , the solution reads where only x A k depends on the input s k . Since both the reservoir first moments and reservoir covariances are linear in the ancilla first moments and covariances, respectively, it is clear that any nonlinearity in these observables is due to the nonlinearity in the ancilla observables of equation (22) only. By direct observation of this equation we may then state that there can be a nonlinear input dependence in covariances only if r or ϕ (or both) are functions of the input, and a nonlinear input dependence in first moments only if arg(α) is a function of the input. Otherwise, there is only a linear input dependence in these reservoir observables. Depending on the input encoding, different linear and nonlinear dependencies on the input sequence history, i.e. memory, can arise in this system as reported in reference [56]. We note in passing that displaced squeezed states have been found to be useful for other forms of quantum machine learning as well since, e.g., the Hilbert-Schmidt distance between two such states is nonlinear in both displacement and squeezing [96][97][98].
It should be pointed out that since each term in the sum in equation (26) is a function of the input at only a single time step, there are no cross terms between different time steps as in, e.g., s k s k−1 . This is detrimental to the performance in tasks that require such terms. As pointed out in [99], this drawback can be avoided by forming the output from reservoir observables of the form x R,(i)xR,(j) +x R,(j)xR,(i) , which introduces such cross terms as long as the reservoir first moments depend on the input. This can be seen from equation (21) by solving for such terms in the left equation.
At this point, we may briefly compare the nonlinearity found in the case of oscillators in Gaussian states interacting according to a quadratic Hamiltonian to that of the spins analyzed in the previous section. In the latter case and with encoding to pure states, there is a richer nonlinear input history dependence since the observables in equation (6) depend on a sum of products of linear and nonlinear functions of the most recent input and the previous state, which is further shown to lead to products of these functions for different time steps with themselves in equation (13). In contrast, the first moments and covariances of the oscillators are linear combinations of such terms and exhibit a limited range of nonlinear input history dependencies.
The reservoir response can be checked as in figure 2, i.e. a single Hamiltonian is generated with arbitrary random fixed links g (ij) along with the input sequence up to the step k − 1. The final element s k is then varied. Specifically, we consider a small network of N = 3 identical unit mass oscillators interacting with springlike couplings, one of them being the ancilla. The bare frequency of all oscillators is ω 0 = 0.25 and interaction strength between some oscillator i and some oscillator j are g (ij) = g (ji) ∈ [0, 0.2], distributed independently and uniformly. The initial state is the ground state of the network Hamiltonian. RC requires in particular that the reservoir state becomes independent of the initial conditions at the limit of many inputs. Looking at equation (26), it is clear that this requires that ρ(A) < 1, where ρ(A) is the spectral radius of matrix A. This ensures that lim m→∞ A m = 0. In fact, it can be shown [56] that this condition is both necessary and sufficient to ensure that the reservoir observables become well-approximated by continuous functions of only a finite number of past inputs. This in turn guarantees, e.g., that the reservoir state never diverges as long as the input is bounded and the spectral radius condition is satisfied. In practice, the condition can be satisfied for a given network Hamiltonian H by tuning the time between inputs Δt.
In figure 5, we show all together both the spectral radius of matrix A and reservoir response of one of the quantum harmonic oscillators as a function of Δt. The response is checked considering a single covariance of the reservoir corresponding to the position of the first reservoir oscillator. The input has been injected using squeezed vacuum states for the ancilla where r k = 1, ϕ k = 2πs k . The diagonally hatched regions indicate values that are finite but higher than the plot range and horizontally hatched regions indicate regimes where reservoir dynamics diverges. There is a clear connection between the hatched regions and the spectral radius ρ(A). Indeed, they indicate regions where ρ(A) is large. Near the boundaries where the spectral radius approaches values that are too high, the reservoir seems to become unresponsive to the input. This might be due to a relatively strong dependency on multiple earlier inputs as more and more terms in the sum of equation (26) become significant. In the particular case at hand, these appear to be approximately the boundaries defined by ρ(A) 0.983 as indicated by the shaded regions under the curve. The regions include ρ(A) 1, indicated by the dashed horizontal line, where RC becomes impractical. In the regime useful for RC where the spectral radius is sufficiently small, the response is nonlinear in s k as hinted by the observable responses in the bottom part of figure 5. This is a direct consequence of the ancilla covariances being trigonometric functions of the input. In particular, when the spectral radius is relatively small only a few of the terms of equation (26) survive, which can be expected to lead to a stronger response to the most recent input. This can be observed, e.g., near the region Δtω 0 ∈ [2, 2.5].
In figure 6, we have fixed Δt to correspond to the global minimum of the spectral radius in figure 5 and consider the response of both reservoir oscillators to the input for multiple different input encodings. In panel (a) we have set |α k | = s k , arg(α k ) = 0 and consider all reservoir first moments. The response is linear, as expected. In panels (b) and (c) the ancilla is prepared in squeezed vacuum states whereas the diagonal elements of the covariance matrix σ (ii) (x R k ) = x 2 k (i) − ( x k (i) ) 2 are considered. In (b) we have set r k = s k , ϕ k = 0 which leads to a nonlinear response since this makes the ancilla covariances hyperbolic functions of the input; notice that logarithmic scaling is used. In (c) phase encoding is used with r k = 1, ϕ k = 2πs k , as in figure 5, and the response is similar. Finally in (d) the input is encoded into phase of displacement according to |α k | = 0, arg(α k ) = 2πs k but reservoir second moments x 2 k (i) are shown with logarithmic scaling. The observed shape originates from the dependency of x 2 k (i) on ( x k (i) ) 2 where each x k (i) is a trigonometric function of s k . Finally, it should be pointed out that the analysis presented above for the single ancilla case generalizes directly to the case of multiple ancillae; in particular the forms of equation (24) through (26) still apply. The nonlinearities analyzed and exemplified here can be harnessed for solving typical benchmark tasks with high performance, as was done in [99]. Here it was also found that high performance comes mostly from choosing an encoding compatible with the requirements of the task at hand, leaving considerable freedom over the form of the Hamiltonian (here actually restricted to linear modes dynamics). This further increases the experimental Figure 6. Response of the harmonic oscillator reservoir to the input at time step k. Here observables from all reservoir oscillators (first and second) are shown for a random but fixed network of N = 3 oscillators with the third oscillator as ancilla. In panel (a) input is encoded to coherent states of the ancilla through the magnitude of displacement, which induces a linear response in reservoir first moments. Panels (b) and (c) focus on encoding the input to squeezed vacuum states, which affects reservoir covariances. The encoding is to the magnitude and phase of squeezing, respectively. In panel (d) the encoding is again to coherent states via the phase of displacement but reservoir second moments are shown. Notice that in panels (b) and (d) logarithmic scaling is used. The elementsx (i) with i = 1, . . . , 4 correspond to the observablesq 1 ,p 1 ,q 2 , andp 2 , respectively. convenience besides using only Gaussian states. Allowing for correlations between the ancillae might introduce new ways to encode the input, however the limitations of the nonlinearity in the reservoir first moments and covariances still hold, i.e. they cannot depend on products of the input at different timesteps.

Conclusions
Quantum complex networks can provide the necessary ingredients for a good performance in RC. High dimensionality is guaranteed by making use of a quantum system by an exponential increase of the Hilbert space with the number of nodes or constituents. Furthermore, several proposals in different systems have numerically shown the possibility to introduce nonlinearities and to have the necessary memory to perform standard machine learning benchmark tasks.
In this work, we have complemented the previous studies by showing in an analytical and explicit way the input-output map nonlinear dependence for quantum networks of qubits and oscillators, respectively. The knowledge of the functional form of the observables of both systems depending on the input has clarified the role of all factors contributing to the presence of nonlinearity and it has also made visible the dependence on past inputs, which is strongly related to the memory of the system.
In the case of the network of qubits, we have presented a formalism dependent on the specific map adopted but independent of the particular Hamiltonian to derive the explicit functional dependence of the quantum state of the system and its observables on the input. In general, the reinitialization of the state of one or more qubits to encode the input leads to the presence of nonlinearities in all types of observable quantities. The form of the nonlinear contributions, e.g. polynomial or sinusoidal, can be tuned by proper choice of the quantum state used for encoding the input. As an exception, with a single-qubit encoding and using a mixed state there are only linear dependences, but, remarkably, the nonlinearity is recovered with a multiqubit encoding. This makes particularly evident the substantial effect the form of the input encoding can entail. We have also shown that the form of the Hamiltonian and, in particular, being in different dynamical phases, plays a relevant role to determine the response of the diverse output observables to the input. Furthermore, it allows us to tailor the optimal input strategy depending on the specific task at hand and the required degrees of nonlinearity.
In the continuous-variable case, for a bosonic network of harmonic oscillators considering only Gaussian states, the importance of the input encoding becomes apparent in our analysis. We show how encoding the input in quantum properties as the magnitude and phase of squeezing explicitly determines the expressions for the covariances and first moments of the state of the system. The input-induced nonlinearities in these output quantities arise from their functional dependence in the ancilla used for codifying the input. The presence of nonlinear terms in the covariances is determined by the squeezing, whereas in the first moments it comes from encoding in the displacement of the ancilla state. In contrast to the qubit network, we have found that the dependence in the input history consists in a sum of functions without crossed terms between consecutive inputs. This limitation, however, could be overcome with further post-processing of the output. Beyond the general functional form, we have illustrated the influence of the dynamics by means of numerical calculations. For a good performance, the system is required to work in a regime in which it forgets its initial state. The condition to forget the initial state and to have a convergent dynamics depends on the Hamiltonian of the system and the amount of time the system is allowed to evolve before the next input injection. Therefore, the role of the system's dynamical properties in establishing the dependence on the history of the system of the output observables is again found to be decisive for a proper operation of quantum reservoir computers.
The analysis presented in this work could be extended in several promising directions and assist in the search for the best design of QRC platforms. In particular, the form of encoding can include spatial and temporal multiplexing strategies, temporal parameter modulation, or can even be extended to quantum inputs [88]; physical reservoirs of qudits systems, including more complex network interactions or operating beyond continuous-variable Gaussian states [57,86] will also provide further forms on nonlinearity; finally, also the observable choice of the output layer plays a key role. Models with increasing complexity can entail an improved performance of QRC in tasks requiring nonlinear memory, but generally, this can only be assessed numerically. A major advantage of the platforms considered here, both with qubits and continuous variables, is the possibility to achieve analytical insight, being therefore convenient models to test novel strategies in QRC.