# Physical reservoir computing using finitely-sampled quantum systems

Saeed A. Khan  Fangjun Hu  Gerasimos Angelatos  Hakan E. Türeci  Department of
Electrical and Computer Engineering, Princeton University, Princeton, NJ 08544

(August 27, 2024)

###### Abstract

The paradigm of reservoir computing exploits the nonlinear dynamics of a
physical reservoir to perform complex time-series processing tasks such as
speech recognition and forecasting. Unlike other machine-learning approaches,
reservoir computing relaxes the need for optimization of intra-network
parameters, and is thus particularly attractive for near-term hardware-
efficient quantum implementations. However, the complete description of
practical quantum reservoir computers requires accounting for their placement
in a quantum measurement chain, and its conditional evolution under
measurement. Consequently, training and inference has to be performed using
finite samples from obtained measurement records. Here we describe a framework
for reservoir computing with nonlinear quantum reservoirs under continuous
heterodyne measurement. Using an efficient truncated-cumulants representation
of the complete measurement chain enables us to sample stochastic measurement
trajectories from reservoirs of several coupled nonlinear bosonic modes under
strong excitation. This description also offers a mathematical basis to
directly compare the computational capabilities of a given physical reservoir
operated across classical and quantum regimes. Applying this framework to the
classification of quantum states of systems that are part of the same
measurement chain as the quantum reservoir computer, we assess and explain
measurement-contingent advantages and disadvantages of reservoir processing in
quantum regimes. Our results also identify the vicinity of bifurcation points
as presenting optimal nonlinear processing regimes of an oscillator-based
quantum reservoir. The considered models are directly realizable in modern
circuit QED experiments, while the framework is applicable to more general
quantum nonlinear reservoirs.

##  I Introduction

Reservoir computing is a machine-learning paradigm that emphasizes learning
efficiency [1, 2, 3, 4]: a linear combination of the accessible degrees of
freedom of a physical ‘reservoir’ is the only quantity trained for a given
task. Crucially, foregoing optimization of internal parameters does not
necessarily hinder the computational capacity of reservoir computers, which in
recent years have enabled efficient approaches to forecasting [5, 6, 7],
inference [8, 9], control [10], and similar resource-intensive signal
processing tasks [11, 12, 13, 14, 15]. By relaxing control over microscopic
parameters of the learning machine, reservoir computing embraces the
fundamental notion of computation with very general dynamical systems:
physical or artificial dynamical systems stimulated by incoming time-dependent
signals $\mathbf{u}(t)$ perform computation by transforming an input stream
into an output stream, effectively computing a function
$\mathcal{F}\\{\mathbf{u}(t)\\}$ [see Fig. 1(a)]. If the dynamical system is
sufficiently complex and high-dimensional, the expressive power of even a
linear classifier acting on the output state can enable approximation of a
large set of functions $\mathcal{F}\\{\cdot\\}$ [16, 17, 18]. In this vein,
diverse dynamical systems ranging from early mathematical black-box reservoirs
[1] to physical systems across optical [19, 20, 21, 22, 23, 24], electrical
[25, 26, 27, 28] and recently superconducting [29, 30] platforms have all been
explored as reservoir computers [31, 32, 33].

This flexibility has naturally led to the examination of quantum systems as
potential physical substrates for reservoir computing [34, 35, 36, 37, 38,
39], with the promise of exponentially-scaling reservoir dimensionality. Early
implementations of quantum reservoirs [see Fig. 1(a)] consider computation
with explicitly quantum components (e.g. interacting spins, bosonic or
fermionic degrees of freedom) and encode the input $\mathbf{u}(t)$ in the
Liouvillian governing state dynamics, or even in the quantum reservoir state
directly as a ‘quantum’ input. Reservoir output is often constructed assuming
access to expectation values of reservoir dynamical variables that correspond
loosely to measurable observables (e.g. occupation numbers via photodetectors,
quadratures via homodyne measurement).

While such schemes have been used to theoretically show certain performance
advantages compared to classical systems [40], critical foundational questions
have yet to be addressed. The experimentally-accessible degrees of freedom of
a quantum system are constrained given a specific input and output scheme,
which must be taken into account. Importantly, the role of the quantum theory
of measurement in sampling outputs from such a system for machine-learning
remains unexplored. Furthermore, a quantum reservoir has yet to be treated
from the perspective of a physically-motivated classical limit, the need for
which is twofold. First, this provides the basis to unequivocally identify any
possible advantages afforded by quantum reservoir processing. Secondly, the
energy expended during a given computational task is a fundamental concern in
physical computing; when considering computation with decreasing signal
powers, the transition from classical to quantum regimes of reservoir
operation must naturally be addressed.

This paper introduces a framework for quantum reservoir computing that allows
the same physical reservoir to be controllably operated across quantum or
classical regimes. To satisfy this requirement, we present an analysis built
upon an intrinsically quantum-mechanical description of the complete
measurement chain, including not only the physical reservoir but also inputs
and measured outputs necessary to define a complete computational framework.
The classical limit is then identified where the state of the measurement
chain during computation is well-approximated as a product of coherent states.
This unified approach allows us to define performance metrics (such as the
fidelity of executing a particular task, or the time-to-solution) that account
for limitations due to a prescribed input-output scheme and measurement noise,
and can be compared across classical or quantum regimes of operation of the
same physical reservoir. Such a comparative metric provides the theoretical
basis to identify task-specific and measurement-contingent advantages (or
disadvantages, as we shall see) of operating physical reservoirs in the
quantum regime.

The quantum-mechanical treatment of the complete measurement chain we require
demands the specification of details often simplified in earlier studies. A
representative layout of the quantum measurement chain we wish to describe is
shown in Fig. 1(b), depicted here in the circuit QED (cQED) architecture [41,
42] (many components are generalizable to other platforms as well). The
quantum reservoir computer (QRC) is modelled as a $K$-node quantum nonlinear
device, whose driven-dissipative dynamics are governed by a Liouvillian
superoperator $\mathcal{L}_{\rm QRC}$. Inputs to the QRC must be transmitted
via classical control fields used to interact with the quantum measurement
chain; while classical inputs can be directly encoded into these control
fields (shown here at a single frequency $\omega_{d}$), the question of
quantum inputs (also referred to as ’quantum data’ in earlier work [36]) is
more involved. We define quantum inputs as states of quantum systems
(described by a Liouvillian superoperator $\mathcal{L}_{\rm S}$) generated in
the same cryogenic environment as the QRC. The flow of signals from this
system to the reservoir is accurately modeled via a general coupling
superoperator $\mathcal{L}_{c}$, which typically includes directional
elements. This description of the composite measurement chain accurately
accounts for correlations between the quantum system and the QRC, including
the extreme case where the system and the quantum reservoir are entangled, a
scenario with no analog in classical reservoir computing. Evolution of the QRC
under classical or quantum inputs realizes a high-dimensional nonlinear
mapping of the input to the reservoir’s state space.

Information encoded in the QRC state is extracted from the measurement chain
via continuous heterodyne monitoring of the QRC nodes, which requires us to
address limitations on measured outputs imposed by quantum theory. These
measurements do not directly yield ensemble-averaged expectation values of
reservoir state variables, but noisy output field measurement records that
include quantum fluctuations non-linearly transformed by the quantum reservoir
itself, and are dependent upon the conditional quantum state of the
measurement chain. All this is accounted for by formulating and solving for
the conditional evolution of all dynamical elements in the quantum measurement
chain under the prescribed measurement scheme (often called quantum
trajectories [43] subject to measurement backaction). The final step in our
quantum reservoir computing framework is an output layer applied to these
measured outputs, performing linear operations that are trained on finite
measurement records and not directly on expectation values. The latter in
principle assumes access to an infinite number of measurement records,
ignoring a key resource cost of quantum reservoir computing, and for higher-
order expectation values can introduce classical nonlinear operations, whose
computational influence should be distinguished from that of the nonlinear
quantum reservoir alone. This output layer can be implemented in software or,
for real-time processing, via a customized closed-loop signal processing
hardware such as FPGAs at room temperature.

Starting from the formal mathematical description of the complete measurement
chain via a stochastic master equation (SME) [44, 45], we develop an efficient
computational framework to simulate its use for quantum reservoir computing.
By considering a representative machine-learning task of quantum state
classification with small-scale QRCs, we demonstrate computation with quantum
reservoirs on quantum inputs across classical and quantum regimes of reservoir
operation. Our principal result is that advantages of operating reservoirs in
the quantum regime are task-dependent: once limitations due to measurement are
taken into account, not all state classification tasks immediately benefit
from quantum reservoir processing. Specific tasks, which we identify, can be
carried out with reduced measurement resources in quantum regimes of reservoir
operation.

Before presenting the complete analysis of this quantum reservoir computing
framework that constitutes the main body of this paper, in the remaining
subsections of this introduction we preview the technical aspects that make
this analysis possible, and a summary of its main results.

###  I.1 Technical approach

The evolution of the conditional quantum state $\hat{\rho}^{c}$ of the
composite measurement chain under continuous measurement shown in Fig. 1(b) is
formally described by the SME

$\displaystyle d\hat{\rho}^{c}=\mathcal{L}_{\rm
S}\hat{\rho}^{c}\leavevmode\nobreak\
dt+\mathcal{L}_{c}\hat{\rho}^{c}\leavevmode\nobreak\ dt+\mathcal{L}_{\rm
QRC}\hat{\rho}^{c}\leavevmode\nobreak\ dt+\mathcal{S}_{\rm
meas}(dW)\hat{\rho}^{c}$ (1)

where $\mathcal{S}_{\rm meas}(dW)$ is the stochastic measurement
superoperator. For extended multi-mode and generally nonlinear networks of
interest, such as the one presented by the QRC setup of interest in Fig. 1(b),
Eq. (1) is typically unfeasible to simulate, limited by growing Hilbert space
constraints. To accurately yet efficiently account for dynamics of the full
measurement chain, a truncated-cumulants representation of its conditional
quantum state is developed; the choice of representation is strongly dependent
on the measurement scheme and Liouvillian dynamics in Eq. (1). This
description extends earlier work on continuously-monitored linear modes [46,
47] to nonlinear multimode quantum systems, yielding a set of stochastic
equations for measured observables that scales much more favourably with
increasing system complexity. The accuracy of this representation requires
specific conditions on the nonlinearity of the quantum reservoir and its
operating conditions which, as discussed below, can be readily satisfied
within the cQED architecture. Additionally, this framework admits a well-
defined classical limit, where the joint quantum state of the measurement
chain is well-approximated as a product of coherent states. Defined in this
way, the classical limit can be analyzed without having to redefine the
quantum-classical boundary (colloquially, the Heisenberg-von-Neumann cut)
separating the measured system from the classical observer, which we always
place at the termination of the measurement chain.

Figure 1: (a) General quantum reservoir computing frameworks. (b) Schematic
representation of quantum reservoir computing framework under continuous
measurement, describing a $K$-node quantum reservoir computer (QRC) as part of
a measurement chain for processing signals encoded in classical control fields
such as the single coherent drive (frequency $\omega_{d}$), or generated by
quantum systems in the same quantum environment. The dashed boundary separates
classical controls at room temperature (RT) from quantum devices at cryogenic
temperatures (mK). The quantum system is governed by Liouvillian
$\mathcal{L}_{\rm S}$, and is coupled to the QRC via $\mathcal{L}_{c}$. The
QRC nodes undergo heterodyne measurement at the completion of the measurement
chain, leading to its conditional evolution under measurement superoperators
$\mathcal{S}_{\rm meas}(dW)$. (c) The QRC is used to classify qubit pointer
states generated in a cavity (task I), and states with equal phase and mean
amplitude, but distinct variances (task II).

###  I.2 Machine-learning tasks using quantum inputs

As a quantum information processing task, we consider the problem of
classifying $C$ states of a quantum system (governed by $\mathcal{L}_{\rm S}$)
embedded in the same measurement chain as the QRC. These states
$\hat{\rho}^{(\sigma)}_{\rm S}(t)$ (indexed by $\sigma=1,\ldots,C$) are
dynamically generated by inputs to the measurement chain, and also evolve
under the backaction of quantum measurement of the QRC nodes. As such, they
are defined via the conditional quantum state of the measurement chain by
tracing out the QRC, $\hat{\rho}^{(\sigma)}_{\rm S}(t)={\rm tr}_{\rm
QRC}\\{\hat{\rho}^{c(\sigma)}(t)\\}$. Note that this is simply a formal
definition: our computational approach does not trace out any degrees of
freedom of the measurement chain. For concreteness, we consider the
classification of system states that (I) differ in their mean amplitudes and
phases, and (II) differ only in their variances [see Fig. 1(c)]. Task I is
closely related to the classification of pointer states for dispersive qubit
readout in cQED [48], while Task II can be applied to distinguishing entangled
and non-entangled quantum states [49]. In both cases, suitable physical QRCs
can process these quantum inputs to produce distinct measured outputs, whose
distributions include non-linearly transformed quantum fluctuations due to the
action of the nonlinear QRC itself. To successfully implement the function
$\mathcal{F}\\{\hat{\rho}_{\rm S}^{(\sigma)}(t)\\}=\sigma$ for quantum state
classification, as our framework does, the linear output layer must be trained
to take these non-classical fluctuations into account.

###  I.3 Performance across classical and quantum operating regimes

The quantum reservoir computing framework presented here is constructed to
assess computation fidelity as a function of quantum reservoir size $K$, and
more importantly to explore the role of quantum resources such as entanglement
[50] amongst the reservoir nodes. However, the generation and preservation of
entanglement itself requires additional resources and operating constraints
[51] on the QRC. As a first study, we find that comparing QRC performance
across the transition from a classical regime of QRC operation to a quantum
regime without necessarily demanding entanglement within the quantum reservoir
already yields important insights into the role of quantum reservoir
processing. Furthermore, this analysis provides the basis for future work to
precisely identify any further advantages of reservoir entanglement.

We show that the input power to the measurement chain provides a natural,
physically-motivated parameter to explore this classical-to-quantum
transition: with decreasing input power, maintaining computational accuracy
for quantum state classification requires increasing the QRC nonlinearity,
which we show corresponds to the quantum regime of QRC operation. Importantly,
across this transition, QRC performance on Tasks I and II shows important
qualitative differences. Task I, which requires the QRC to distinguish the
amplitudes and phases of quantum states, fares better in the classical regime.
Task II, on the other hand, requires the QRC to distinguish information
encoded in quantum fluctuations. We find that this task can be performed with
reduced measurement resources in increasingly quantum regimes (for details see
Sec. IV.4, V.4 respectively).

###  I.4 Nonlinear processing near classical bifurcation points

While the two quantum state classification tasks probe distinct nonlinear
processing capabilities of the QRC, we find that useful regimes of reservoir
operation are unified by their proximity to classical bifurcation points. Near
these bifurcation phases, we find an enhancement of the nonlinear response of
the QRC to its inputs; this suggests a general signature to identify parameter
regimes of optimal reservoir performance for a broader class of computational
tasks.

###  I.5 Organization of the paper

The remainder of this paper is organized as follows. In Sec. II, we introduce
the specific cQED realization of the $K$-node QRC and the measurement chain.
In Sec. III, we use a minimal version of this measurement chain, a QRC
comprised of a single nonlinear mode, to develop our theoretical approach to
the simulation of its conditional dynamics, and discuss its regimes of
validity. Sec. IV then considers the use of the QRC for classifying pointer
states, and Sec. V for the classification of general Gaussian states with
distinct variances (Tasks I and II respectively). In both cases, we provide a
characterization of the distinct quantum measurement chains and their measured
outputs, a detailed study of QRC parameters that enable successful
classification, and a comparison of QRC performance across classical and
quantum operating regimes.

##  II Physical Model for a Quantum Reservoir Computer

To describe the QRC, we adopt a model that is realizable within the cQED
architecture, based on $K$ coupled Kerr nonlinear modes furnished by Josephson
junctions [52]. The nonlinear modes have frequencies $\\{\omega_{k}\\}$,
nonlinearity strengths $\\{\Lambda_{k}\\}$, with linear coupling
$\\{\text{\it{g}}_{ij}\\}$ between modes $i$ and $j$, yielding the Liouvillian
governing QRC dynamics,

$\displaystyle\mathcal{L}_{\rm QRC}\hat{\rho}^{c}=-i\Big{[}$
$\displaystyle\sum_{k}\omega_{k}\hat{b}_{k}^{\dagger}\hat{b}_{k}^{\
}-\sum_{k}\frac{\Lambda_{k}}{2}\hat{b}_{k}^{\dagger}\hat{b}_{k}^{\dagger}\hat{b}_{k}^{\
}\hat{b}_{k}^{\ },\hat{\rho}^{c}\Big{]}$ $\displaystyle-i\Big{[}$
$\displaystyle\sum_{ij}\text{\it{g}}_{ij}(\hat{b}_{i}^{\
}\hat{b}_{j}^{\dagger}+\hat{b}_{i}^{\dagger}\hat{b}_{j}^{\
}),\hat{\rho}^{c}\Big{]}+\sum_{k}\gamma_{k}\mathcal{D}[\hat{b}_{k}^{\
}]\hat{\rho}^{c}.$ (2)

The linear damping is described by the standard dissipative superoperator
$\mathcal{D}[\hat{o}]=\hat{o}\hat{\rho}\hat{o}^{\dagger}-\frac{1}{2}\\{\hat{o}^{\dagger}\hat{o},\hat{\rho}\\}$
, with rates $\\{\gamma_{k}\\}$. The decay channels for each node are
continuously monitored with unit efficiency via a heterodyne measurement
scheme, described by the stochastic measurement superoperator

$\displaystyle\mathcal{S}_{\rm meas}(dW)\hat{\rho}^{c}=$
$\displaystyle\sum_{k=1}^{K}\sqrt{\frac{\gamma_{k}}{2}}\left(\hat{b}_{k}^{\
}\hat{\rho}^{c}+\hat{\rho}^{c}\hat{b}_{k}^{\dagger}-\langle\hat{b}_{k}^{\
}+\hat{b}_{k}^{\dagger}\rangle^{c}\right)dW_{k}^{X}(t)+$
$\displaystyle\sum_{k=1}^{K}\sqrt{\frac{\gamma_{k}}{2}}\left(-i\hat{b}_{k}^{\
}\hat{\rho}^{c}+i\hat{\rho}^{c}\hat{b}_{k}^{\dagger}-\langle-i\hat{b}_{k}^{\
}+i\hat{b}_{k}^{\dagger}\rangle^{c}\right)dW_{k}^{P}(t),$ (3)

where $dW_{k}^{X,P}$ are independent Wiener increments, and
$\langle\hat{o}\rangle^{c}$ indicates the expectation value of an arbitrary
operator $\hat{o}$ with respect to the conditional quantum state,
$\langle\hat{o}\rangle^{c}={\rm tr}\\{\hat{\rho}^{c}\hat{o}\\}$. The results
of this continuous measurement are heterodyne currents $J_{k}^{X,P}(t)$ read
out from each reservoir node,

$\displaystyle J_{k}^{X}(t)\leavevmode\nobreak\ dt$
$\displaystyle=\sqrt{\frac{\gamma_{k}}{2}}\langle\hat{b}_{k}^{\
}+\hat{b}_{k}^{\dagger}\rangle^{c}\leavevmode\nobreak\ dt+dW_{k}^{X}(t),$ (4)
$\displaystyle J_{k}^{P}(t)\leavevmode\nobreak\ dt$
$\displaystyle=\sqrt{\frac{\gamma_{k}}{2}}\langle-i\hat{b}_{k}^{\
}+i\hat{b}_{k}^{\dagger}\rangle^{c}\leavevmode\nobreak\ dt+dW_{k}^{P}(t).$ (5)

As described by Eq. (1), the QRC is just one part of the measurement chain.
The signal-generating quantum system is described by the superoperator
$\mathcal{L}_{\rm S}$, and will be task-dependent; we will consider specific
realizations throughout this paper. The coupling superoperator
$\mathcal{L}_{c}$ completes the measurement chain by allowing the measured
quantum system and the QRC to interact. It can once again take an arbitrary
form, allowing for both coherent and dissipative, reciprocal and non-
reciprocal, and even nonlinear couplings. Depending on the type of coupling
interaction, $\mathcal{L}_{c}$ may also introduce local terms affecting the
measured quantum system and the QRC; we will make these clear when they arise.

Having defined every element of the measurement chain in Fig. 1(b), we will
now analyze the nature of its quantum dynamics, introducing an approach to
studying this dynamics that is scalable to extended, multimode quantum
systems.

##  III Quantum dynamics of the single-node Kerr QRC

To analyze the nonlinear processing capabilities of this QRC, it will prove
useful to begin with a simple minimal realization - a single Kerr node - and a
single coherent drive incident on the QRC directly, via a separate input port
with loss rate $\Gamma$. This situation describes QRC processing of classical
signals encoded in the generally time-dependent coherent drive amplitude and
phase. The resulting system schematic is shown in Fig. 2. We begin by
analyzing the unconditional dynamics of the resulting system, described by the
quantum master equation

$\displaystyle d\hat{\rho}\equiv\mathcal{L}\hat{\rho}\leavevmode\nobreak\
dt=\mathcal{L}_{\rm QRC}\hat{\rho}\leavevmode\nobreak\
dt-i[\hat{\mathcal{H}}_{d},\hat{\rho}]\leavevmode\nobreak\ dt$ (6)

where the drive Hamiltonian in the frame rotating at the drive frequency
$\omega_{d}$ is then

$\displaystyle\hat{\mathcal{H}}_{d}=\eta(\hat{b}+\hat{b}^{\dagger}).$ (7)

The Liouvillian $\mathcal{L}_{\rm QRC}$ describing the single-node QRC takes
the form

$\displaystyle\mathcal{L}_{\rm
QRC}\hat{\rho}=-i\left[-\Delta\hat{b}^{\dagger}\hat{b}-\frac{\Lambda}{2}\hat{b}^{\dagger}\hat{b}^{\dagger}\hat{b}\hat{b},\hat{\rho}\right]+(\gamma+\Gamma)\mathcal{D}[\hat{b}]\hat{\rho}$
(8)

where $\Delta=\omega_{d}-\omega$ is the detuning of the drive frequency from
the bare QRC node frequency $\omega$. For simplicity of notation, in what
follows we introduce for convenience the total linear damping
$\widetilde{\gamma}=\gamma+\Gamma$ of the nonlinear mode.

Figure 2: Schematic showing the fundamental node of the proposed QRC as part
of a simplified measurement chain: a ‘system’ described by a single frequency
coherent drive incident on the QRC node via an input channel that is not
measured, and continuous weak measurement of the QRC node output. The QRC node
can be realized in cQED as a tunable-frequency Kerr oscillator using a
capacitively shunted SQUID loop [53].

###  III.1 Equations of motion approach based on truncation of higher-order
cumulants

The quantum state described by an arbitrary density matrix $\hat{\rho}$ can
always be written in terms of normal-ordered cumulants, defined for a single
quantum mode $\hat{b}$ as

$C_{b^{\dagger
p}b^{q}}\equiv\left.\frac{\partial^{p+q}}{\partial(iz^{\ast})^{p}\partial(iz)^{q}}\ln\mathrm{tr}\\{\hat{\rho}e^{iz^{\ast}\hat{b}^{\dagger}}e^{iz\hat{b}}\\}\right|_{z=z^{\ast}=0}.$
(9)

The term in curly braces may immediately be recognized as the characteristic
function of $\hat{\rho}$ [54], associated with normal-ordered phase-space
distributions such as the $P$-representation; this connection is elaborated
upon in Appendix A. The above definition generalizes straightforwardly to
states $\hat{\rho}$ describing multiple bosonic modes.

Cumulants are defined by their order $p+q$, and provide an equivalent set of
variables to describe the dynamics of a quantum state. Specifically, we can
write down the equation of motion for the first-order cumulant of the single-
node QRC ($C_{b}\equiv\langle\hat{b}\rangle={\rm
tr}\\{\mathcal{L}\hat{\rho}\hat{b}\\}$),

$\displaystyle\langle\dot{\hat{b}}\rangle$
$\displaystyle=\left(i\Delta-\frac{\widetilde{\gamma}}{2}\right)\langle\hat{b}\rangle+i\Lambda\langle\hat{b}^{\dagger}\rangle\langle\hat{b}\rangle\langle\hat{b}\rangle-i\eta$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
+i\Lambda\left(C_{bb}\langle\hat{b}^{\dagger}\rangle+2C_{b^{\dagger}b}\langle\hat{b}\rangle+C_{b^{\dagger}bb}\right)$
(10)

where
$C_{bb}=\langle\hat{b}\hat{b}\rangle-\langle\hat{b}\rangle^{2},C_{b^{\dagger}b}=\langle\hat{b}^{\dagger}\hat{b}\rangle-\langle\hat{b}^{\dagger}\rangle\langle\hat{b}\rangle$,
and
$C_{b^{\dagger}bb}=\langle\hat{b}^{\dagger}\hat{b}\hat{b}\rangle-2\langle\hat{b}^{\dagger}\hat{b}\rangle\langle\hat{b}\rangle-\langle\hat{b}\hat{b}\rangle\langle\hat{b}^{\dagger}\rangle+2\langle\hat{b}^{\dagger}\rangle\langle\hat{b}\rangle^{2}$.
Expressions for higher-order cumulants become increasingly more unwieldy, but
can be systematically obtained, as discussed in Appendix C. We see immediately
that the first-order cumulant $\langle\hat{b}\rangle$ couples to cumulants of
second and higher order. We can similarly obtain equations of motion for
unique second-order cumulants $C_{b^{\dagger}b},C_{bb}$,

$\displaystyle\dot{C}_{b^{\dagger}b}$
$\displaystyle=-\widetilde{\gamma}C_{b^{\dagger}b}-i\Lambda(C_{bb}\langle\hat{b}^{\dagger}\rangle^{2}-C_{bb}^{*}\langle\hat{b}\rangle^{2})$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ -i\Lambda(C_{b^{\dagger}bb}\langle\hat{b}\rangle-
C_{b^{\dagger}bb}^{*}\langle\hat{b}^{\dagger}\rangle),$ (11a)
$\displaystyle\dot{C}_{bb}$
$\displaystyle=\left(i2\Delta-\widetilde{\gamma}+i\Lambda(1+4|\langle\hat{b}\rangle|^{2}\\!\\!+6C_{b^{\dagger}b})\right)C_{bb}$
$\displaystyle+i\Lambda$
$\displaystyle\langle\hat{b}\rangle^{2}(1+2C_{b^{\dagger}b})+i2\Lambda(C_{b^{\dagger}bbb}+2\langle\hat{b}\rangle
C_{b^{\dagger}bb}+C_{bbb}\langle\hat{b}^{\dagger}\rangle).$ (11b)

Clearly, second-order cumulants couple to cumulants of third and fourth order.
This is a general feature of nonlinear quantum systems: moments and thus
cumulants of a certain order can couple to those of higher-order, leading to
an infinite hierarchy of equations that do not form a closed set.

However, an infinite set of normal-ordered cumulants is not necessarily
required to describe all multimode quantum states. In particular, a multimode
quantum system in a product of coherent states is described entirely by its
nonzero first-order cumulants; all cumulants of order $p+q>1$ vanish (see
Appendix B for derivations). Multimode quantum states that are defined
entirely by their first and second-order cumulants admit Gaussian phase-space
representations, and are labelled Gaussian states. States with nonzero
cumulants of third or higher-order are thus by definition non-Gaussian states
[55, 56].

Our numerical approach leverages this efficient representation of specific
multimode quantum states in terms of cumulants. In particular, we consider an
ansatz wherein the quantum state of the complete measurement is described
entirely by cumulants up to a finite order $p+q\leq n_{\rm trunc}$; all
cumulants of order $p+q>n_{\rm trunc}$ are thus set to zero, truncating the
hierarchy and yielding a closed set of equations for the retained nonzero
cumulants. In this paper, we choose $n_{\rm trunc}=2$ for a quantum
measurement chain defined entirely by its first and second-order cumulants,
although the truncation can similarly be carried out at higher order. The
resulting Truncated Equations Of Motion (TEOMs), and their generalization to
conditional dynamics of multimode quantum systems introduced in Sec. III.3,
form the basis of our computational approach in this paper.

Since we are specifically interested in nonlinear quantum systems for
reservoir computing, which can generate higher-order cumulants in dynamics,
one must ask when such an ansatz may hold. In the two following subsections,
we identify system-specific drive and parameter regimes where this ansatz is
valid, and verify our conclusions via numerical comparisons with exact results
and (S)ME simulations.

###  III.2 Application to unconditional QRC node dynamics

We will use the single coherently-driven QRC node modeled as a Kerr oscillator
to benchmark the TEOMs approach. In doing so, we make extensive use of the
exact quantum steady-state solution of the coherently-driven Kerr oscillator,
obtained via the complex-$P$ representation [57]. This provides access to
steady-state cumulants of arbitrary order (see Appendix E for details), which
can be used to identify regions of parameter space where a truncation of
cumulants beyond second-order can be justified.

One such regime should be furnished by the classical limit of the single Kerr
QRC node, where the QRC is always in a coherent state and thus described
entirely by its first-order cumulants. To identify this classical limit, and
thus explore deviations from it due to quantum effects, it proves convenient
to work with specific scaled quantities: after introducing dimensionless time
$t^{\prime}=\widetilde{\gamma}t$ and energy scales
$(\Delta^{\prime},\Lambda^{\prime},\eta^{\prime})=(\Delta,\Lambda,\eta)/\widetilde{\gamma}$,
we scale
$\langle\hat{b}\rangle^{\prime}\to\sqrt{\Lambda^{\prime}}\langle\hat{b}\rangle$,
and impose our truncated ansatz, following which Eq. (10) becomes

$\displaystyle\frac{d}{dt^{\prime}}\langle\dot{\hat{b}}\rangle^{\prime}$
$\displaystyle=\left(i\Delta^{\prime}-\frac{1}{2}\right)\langle\hat{b}\rangle^{\prime}+i\langle\hat{b}^{\dagger}\rangle^{\prime}\langle\hat{b}\rangle^{\prime}\langle\hat{b}\rangle^{\prime}-i\mathcal{N}$
$\displaystyle+i\Lambda^{\prime}\left(C_{bb}\langle\hat{b}^{\dagger}\rangle^{\prime}+2C_{b^{\dagger}b}\langle\hat{b}\rangle^{\prime}\right).$
(12)

The first line in Eq. (12) makes no reference to second-order cumulants of the
QRC node state; it thus describes dynamics when the QRC is in a coherent state
(for which such cumulants vanish), corresponding to the classical limit.
Dynamics in this limit depend on the drive and the nonlinearity via a single
dimensionless parameter [58, 59],

$\displaystyle\mathcal{N}=\eta^{\prime}\sqrt{\Lambda^{\prime}}=\frac{\eta}{\widetilde{\gamma}}\sqrt{\frac{\Lambda}{\widetilde{\gamma}}}.$
(13)

The second line of Eq. (12) does however depend on cumulants; these terms
describe quantum dynamics that lead to a deviation from purely coherent state
evolution. In these scaled units, the TEOMs for the second-order cumulants
obtained from Eq. (11a), (11b) take the form

$\displaystyle\dot{C}_{b^{\dagger}b}$
$\displaystyle=-C_{b^{\dagger}b}-i(C_{bb}(\langle\hat{b}^{\dagger}\rangle^{\prime})^{2}-C_{bb}^{*}(\langle\hat{b}\rangle^{\prime})^{2}),$
(14a) $\displaystyle\dot{C}_{bb}$
$\displaystyle=\left(i2\Delta^{\prime}-1+i(\Lambda^{\prime}+4|\langle\hat{b}\rangle^{\prime}|^{2}\\!\\!+6\Lambda^{\prime}C_{b^{\dagger}b})\right)C_{bb}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
+i(\langle\hat{b}\rangle^{\prime})^{2}(1+2C_{b^{\dagger}b}).$ (14b)

Consider now the transformation where $\Lambda^{\prime}\to 0$, while
$\mathcal{N}$ is held fixed (by increasing the drive strength $\eta^{\prime}$
simultaneously). Terms in the first line of Eq. (12) describing classical
dynamics remain unchanged, while those in the second line describing quantum
deviations become smaller, provided second-order cumulants do not grow with
$\Lambda^{\prime}$. From Eqs. (14a), (14b), it is clear that
$C_{b^{\dagger}b}$ is independent of $\Lambda^{\prime}$, while $C_{bb}$ has a
dependence that is negligible as $\Lambda^{\prime}\to 0$ for
$|\langle\hat{b}\rangle^{\prime}|\neq 0$. Hence by keeping $\mathcal{N}$ fixed
while decreasing the nonlinearity strength, the influence of second-order
cumulants on dynamics of the QRC node amplitude is suppressed, describing the
classical limit of QRC operation. Physically, this transformation leads to an
increase in the unscaled drive $\eta^{\prime}$ and QRC node amplitude
$|\langle\hat{b}\rangle|$, and hence occupation number, decreasing the
relative impact of quantum fluctuations in agreement with conventional notions
of classicality. Conversely, increasing nonlinearity strength for fixed
$\mathcal{N}$ should enhance the impact of quantum fluctuations, marked by a
systematic increase in second-order cumulants, and taking the Kerr QRC
gradually towards a quantum regime of operation.

To explore this transition, we start with the classical limit of the QRC as
determined by terms in the first line of Eq. (12), for which the steady-state
phase diagram can be easily computed (see Appendix G). This phase diagram,
shown in $(\Delta,\mathcal{N})$-space in the center panel of Fig. 3(a), is the
well-known result for a coherently-driven classical Kerr oscillator [58, 59]:
the orange region depicts the classical bistability, which emerges for
sufficiently negative detuning given our choice of the sign of the Kerr
nonlinearity.

Figure 3: Classical and quantum dynamics of a single Kerr QRC node. In this
plot, colored solid curves are exact complex-$P$ results, colored dots are
TEOMs simulations in the long-time limit, and dashed orange curves are
classical results. Other curves are specified when relevant. (a) Center panel:
Classical phase diagram of the single coherently-driven Kerr QRC node. Left
panel: QRC amplitude as a function of $\mathcal{N}$ for fixed detuning
$\Delta/\gamma=-1.0$ (across indicated vertical cross-section of the phase
diagram). Right panel: QRC amplitude as a function of detuning $\Delta$ for
fixed $\mathcal{N}=0.385$ (across indicated horizontal cross-section of the
phase diagram). Dashed black lines in both plots indicate the response for a
linear oscillator. (b) Steady-state scaled amplitude
$|\langle\hat{b}^{\prime}\rangle|$ for nonlinearity strengths
$\Lambda/\gamma\in[0.005,0.02,0.05]$, as a function of $\mathcal{N}$ for fixed
$\Delta/\gamma=-1.0$ (region within the gray box in (a), left panel). (c) Same
as (b) but as a function of detuning $\Delta$ for fixed $\mathcal{N}=0.385$
(region within the gray box in (a), right panel). (d) Steady-state second-
order cumulants ($C_{b^{\dagger}b}$, purple, and $|C_{bb}|$, green) and higher
(third, fourth) order cumulants (dashed gray) across the cross-section in (b);
nonlinearity increases from left to right. (e) Same as (d) but now for the
cross-section in (c).

This bistability is a feature specific to nonlinear systems; in its vicinity,
we can expect an enhanced nonlinear response of the QRC observables to QRC
inputs. This is explicitly observed in the scaled field amplitude
$|\langle\hat{b}\rangle^{\prime}|$, plotted in Fig. 3(b) as a function of
$\mathcal{N}$ for fixed $\Delta/\gamma=-1.0$ (that is, across the indicated
vertical dashed line of the classical phase diagram). The solid blue curve is
the exact complex-$P$ result for $\Lambda/\gamma=0.005$, while blue dots are
obtained by evolving the TEOMs, Eqs. (10), (11a), (11b), to their steady state
numerically. For small $\mathcal{N}$, the response follows that of a linear
oscillator (dashed black line), but becomes nonlinear with increasing
$\mathcal{N}$. Also shown in dashed orange is the response of the classical
nonlinear oscillator (explicitly dropping second-order cumulants), which does
not depend on the actual nonlinearity strength but only on $\mathcal{N}$.

The influence of nonlinearity is similarly evident when plotting
$|\langle\hat{b}^{\prime}\rangle|$ as a function of $\Delta$, for fixed
$\mathcal{N}=0.385$ in Fig. 3(c) (across the horizontal dashed line of the
classical phase diagram). The black dashed curve indicates the standard
Lorentzian response of a linear oscillator, centered around the bare frequency
$\omega_{b}$. The response of the nonlinear QRC node is clearly that of a
Lorentzian deformed towards lower frequencies due to the Kerr-induced
frequency shift.

In both results, we find excellent agreement between the exact quantum
results, the TEOMs, and classical nonlinear dynamics. This is indicative of
operating parameters where not only is a truncated description using first and
second-order cumulants sufficient, but the influence of second-order cumulants
on first-order cumulants is negligible as well. Consequently, we expect the
QRC node here to function as a classical reservoir node, still allowing for
nonlinear processing of information encoded in the drive amplitude, phase,
and/or frequency.

If we now keep $\mathcal{N}$ fixed but increase the nonlinearity strength
$\Lambda$, the QRC can be moved to a quantum regime of operation. We plot the
QRC response as a function of $\mathcal{N}$ for $\Delta/\gamma=-1.0$ in Fig.
3(b), and as a function of $\Delta$ for $\mathcal{N}=0.385$ in Fig. 3(c), for
three different nonlinearity strengths $\Lambda/\gamma\in[0.005,0.02,0.05]$.
In both cases, we see that with increasing nonlinearity, the response deviates
from that of the classical nonlinear oscillator (which recall is essentially
equivalent to the quantum result for weak $\Lambda/\gamma=0.005$). From Eq.
(12), it is clear that this deviation is due to the coupling of first-order
moments to second-order cumulants; this is precisely the effect accounted for
by the TEOMs, which thus agree with the full quantum result.

Also shown in Figs. 3(d) and (e) are steady-state cumulants plotted across the
same cross-sections as Figs. 3(b) and (c) respectively, with panels from left
to right indicating stronger nonlinearities. Second-order cumulants are
plotted in green and purple (lines are exact complex-$P$ results, dots are
TEOMs); they clearly display non-monotonic behaviour as a function of
$\mathcal{N}$ and $\Delta$, acquiring large magnitudes for specific operating
parameters in the vicinity of the classical bistability. Physically, second-
order cumulants are related to the maximum and minimum (dimensionless)
quadrature variances of the QRC nodes:
$\frac{1}{2}+C_{b^{\dagger}b}\pm|C_{bb}|$ respectively, and their magnitudes
are thus measures of amplification and squeezing of quantum fluctuations due
to the underlying Kerr nonlinearity.

Third and fourth order cumulants are shown in the same plots in dashed gray
(exact complex-$P$ only). There are two features of note. Firstly, the overall
magnitude of these higher-order cumulants relative to second-order cumulants
appears to increase for stronger nonlinearities. This is an indicator of the
emergence of non-Gaussian features with increasing $\Lambda$, as expected, and
generally leads to greater deviation between TEOMs and exact quantum results
for first and second-order cumulants. Secondly, even for a fixed nonlinearity
strength, the magnitude of cumulants depends on operating parameters
($\mathcal{N}$ and $\Delta$) that lead to operation in the vicinity of the
classical bistability. Hence, while results obtained from TEOMs for such
operating parameters can deviate from the exact quantum solution, away from
these regions the agreement improves, even for strong nonlinearities.

###  III.3 Conditional QRC node dynamics

We have so far used steady-state, unconditional quantities to analyze the
nonlinear processing capabilities of a single Kerr-based QRC node, and
identified parameter regimes where an approach based on truncated EOMs is
valid. However, to allow for processing of time-dependent information by the
QRC, and to analyze the output from the QRC obtained as individual measurement
records, an approach is needed that is able to accurately capture the
conditional, dynamical evolution of the QRC node. This is described by the SME

$\displaystyle d\hat{\rho}^{c}=\mathcal{L}\hat{\rho}^{c}\leavevmode\nobreak\
dt+\mathcal{S}_{\rm meas}(dW)\hat{\rho}^{c}$ (15)

where $\mathcal{S}_{\rm meas}(dW)$ is defined in Eq. (3), and our current
system assumes the $K=1$ case. Expectation values with respect to the
conditional quantum state under measurement $\hat{\rho}^{c}$ are thus also
conditional, which we indicate by the superscript $c$. For example, the
equation of motion for $\langle\hat{b}\rangle^{c}$ is now given by ${\rm
tr}\\{(\mathcal{L}\leavevmode\nobreak\ dt+\mathcal{S}_{\rm
meas}(dW))\hat{\rho}^{c}\hat{b}\\}$, which takes the form (see Appendix D for
details)

$\displaystyle d\langle\hat{b}\rangle^{c}={\rm
tr}\\{(\mathcal{L}\hat{\rho}^{c}\hat{b}\\}\leavevmode\nobreak\ dt+$
$\displaystyle\sqrt{\frac{\gamma}{2}}\left(C^{c}_{b^{\dagger}b}+C^{c}_{bb}\right)dW^{X}(t)+i\sqrt{\frac{\gamma}{2}}\left(C^{c}_{b^{\dagger}b}-C^{c}_{bb}\right)dW^{P}(t).$
(16)

The first line $\propto\mathcal{L}$ simply includes terms from Eq. (10), with
expectation values replaced by their conditional counterparts. Terms in the
second line are due to the measurement superoperator, and render the equation
of motion for $\langle\hat{b}\rangle^{c}$ stochastic, conditioned on the
trajectory-specific realizations of $dW^{X,P}(t)$. Dropping cumulants higher
than second order in Eq. (16), arising in the first line, yields the
stochastic truncated EOMs (STEOMs) for $\langle\hat{b}\rangle^{c}$.

We can similarly obtain truncated equations of motion for the conditional
cumulants $C^{c}_{bb}$, $C^{c}_{b^{\dagger}b}$ under heterodyne measurement:

$\displaystyle dC^{c}_{b^{\dagger}b}$ $\displaystyle=\left[{\rm
tr}\\{\mathcal{L}\hat{\rho}^{c}\hat{b}^{\dagger}\hat{b}\\}\\!-\\!\langle\hat{b}^{\dagger}\rangle^{c}{\rm
tr}\\{\mathcal{L}\hat{\rho}^{c}\hat{b}\\}\\!-\\!\langle\hat{b}\rangle^{c}{\rm
tr}\\{\mathcal{L}\hat{\rho}^{c}\hat{b}^{\dagger}\\}\right]\\!dt$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
-{\gamma}\left[(C^{c}_{b^{\dagger}b})^{2}+C^{c}_{bb}(C^{c}_{bb})^{*}\right]\leavevmode\nobreak\
dt,$ (17a) $\displaystyle dC^{c}_{bb}$ $\displaystyle=\left[{\rm
tr}\\{\mathcal{L}\hat{\rho}^{c}\hat{b}\hat{b}\\}-2\langle\hat{b}\rangle^{c}{\rm
tr}\\{\mathcal{L}\hat{\rho}^{c}\hat{b}\\}\right]dt$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
-2{\gamma}C^{c}_{bb}C^{c}_{b^{\dagger}b}\leavevmode\nobreak\ dt.$ (17b)

Again, terms $\propto\mathcal{L}$ are as found in Eqs. (11a), (11b) post-
truncation, with expectation values replaced by their conditional
counterparts. The second line of each equation describes the evolution due to
measurement, which at first glance appears deterministic: note that Wiener
increments $dW^{X,P}(t)$ make no appearance. For linear quantum systems under
continuous weak measurement, this is in fact the case: second-order cumulants
form a closed set described by the above equations, and no stochastic terms
arise [46, 47, 60, 61]. However, for nonlinear quantum systems of interest
here, second-order cumulants can couple to the stochastic first-order moments
(here, via terms $\propto\mathcal{L}$), rendering the conditional evolution of
second-order cumulants generally stochastic as well.

Figure 4: Benchmarking STEOMs against full SME simulations. Simulated
parameters are as marked by the teal ‘X’ in the phase diagram of Fig. 3(a):
$\mathcal{N}=0.385$, $\Delta/\gamma=-1.0$, for (a) $\Lambda/\gamma=0.005$ and
(b) $\Lambda/\gamma=0.02$. Left panel shows real and imaginary parts of the
conditional expectation value $\langle\hat{b}\rangle$, and center panel shows
conditional second-order cumulants $C_{bb},C_{b^{\dagger}b}$; colored (white)
curves are SME (STEOMs) simulations. Right panel shows measured quadratures
$I^{X}(t)$, $I^{P}(t)$: the SME result is plotted in solid black, while the
STEOMs result is plotted in yellow, with a lighter shade indicating earlier
times. (c) Measured quadratures $\\{I^{X}(t),I^{P}(t)\\}$ at $\gamma t=20.0$
obtained using SME (left panel) and STEOMs (right panel), for
$\mathcal{N}=0.385$, $\Lambda/\gamma=0.02$, and $\Delta/\gamma=-1.0$. (d) Same
as (c) but for $\Delta/\gamma=0.0$.

Eqs. (16), (17a), (17b) define the STEOMs for the single Kerr QRC node. To
assess their validity, we compare their simulation results against integration
of the SME, Eq. (15). For concreteness, we consider the operating point marked
by the teal ‘X’ in the phase diagram of Fig. 3(a): $\mathcal{N}=0.385$,
$\Delta/\gamma=-1.0$. Note that the same noise realizations $dW^{X,P}(t)$ are
used for both simulations to ensure the compared trajectories are conditioned
on the same measurement record. Results are included in Fig. 4(a), (b) for
$\Lambda/\gamma\in[0.005,0.02]$ respectively, showing first-order moments
(left panel) and second-order cumulants (center panel), with STEOMs results in
white, and SME results obtained using QuTiP [62] in color. Both methods agree
very well, especially for first-order cumulants. The increasing fluctuations
in the second-order cumulants are a direct signature of the nonlinearity of
the quantum system under study.

However, in a real experiment, such information about the QRC state must be
extracted from measurement records. For the case of heterodyne measurement,
the obtained records $J^{X,P}(t)$ are defined in Eqs. (4), (5). These records
are typically processed to reduce noise, most commonly via temporal filtering;
the processed records define measured quadratures

$\displaystyle I^{X,P}(t)$
$\displaystyle=\frac{1}{t}\int_{0}^{t}d\tau\leavevmode\nobreak\ J^{X,P}(\tau)$
(18)

which here describe processing of single-shot readout records. Measured
quadratures obtained using a single set of measurement records are shown in
the measured phase space given by $(I^{X}(t),I^{P}(t))$ in the right panel of
Fig. 4(a), (b) $\Lambda/\gamma\in[0.005,0.02]$ as before; SME results are
shown in solid black, with the STEOMs results plotted on top, using a
colorscale going from black to yellow indicating increasing time. Both show
excellent agreement; this is unsurprising since $J^{X,P}(t)$ yield information
about the underlying conditional dynamics of $\langle\hat{b}\rangle^{c}$,
which also agree very well between the two methods.

Note that the measured quadratures are themselves stochastic quantities, whose
statistics are correlated with the underlying measured system state. This can
be seen via the measured quadrature distributions in Fig. 4(c), obtained by
simulating several measurement records (here $100$ in total) and plotting
$\\{I^{X}(t),I^{P}(t)\\}$ at $\gamma t=20.0$, for $\Delta/\gamma=-1.0$ as
before. At this operating point near the classical bistable region, the Kerr
nonlinearity amplifies the magnitude of cumulants
$C_{bb}\equiv|C_{bb}|e^{i\phi_{bb}}$ [see Fig. 3 (d)] leading to squeezing of
the internal QRC field along the axis determined by $\phi_{bb}/2$ (and
amplification along the orthogonal quadrature). This internal QRC squeezing
manifests as squeezing of the measured quadrature distribution; the squeezing
axis is unchanged since the temporal filter defining measured quadratures via
Eq. (18) is quadrature-agnostic, and thus only reduces the overall noise power
while preserving its relative strength amongst measured quadratures. On the
other hand, for $\Delta/\gamma=0.0$, further away from the bistable region,
$C_{bb}$ is much smaller in magnitude and the measured quadrature distribution
in Fig. 4(d) exhibits no squeezing. Excellent agreement between SME (left) and
STEOMs (right) is observed for both cases.

###  III.4 Validity and scalability of (S)TEOMs

The key observations from our analysis of Kerr QRC operating regimes and the
generation of higher-order cumulants (summarized in Fig. 3) are twofold.
Primarily, the magnitude of nonzero higher-order cumulants increases with the
strength of the Kerr nonlinearity relative to the QRC node damping rate,
$\Lambda/\widetilde{\gamma}$, which coincides with the QRC under coherent
driving transitioning from classical to quantum regimes. However, operating
regimes in the vicinity of the classical bistability also lead to an
enhancement in higher-order cumulants, regardless of nonlinearity strength.

To be confident of the validity of the truncated cumulants approach, we
therefore analyze QRCs with nonlinearity strengths
$\Lambda/\widetilde{\gamma}\leq 0.02$, and operate near but not within the
classically-bistable region. This truncated cumulants approach provides a
highly efficient mathematical description when applied to multimode quantum
systems. For a measurement chain comprising an $N$-mode quantum system coupled
to a $K$-node QRC, the composite system of $N_{\rm T}=N+K$ total quantum modes
is described by $2N_{\rm T}^{2}+3N_{\rm T}$ unknowns using (S)TEOMs, scaling
quadratically with $N_{\rm T}$ instead of the exponential growth in required
Hilbert space size for full (S)ME simulations. This description also places no
constraints on modal occupation numbers, enabling our exploration of a well-
defined classical limit of the coherently-driven Kerr QRC. Our (S)TEOMs
approach is built to be a scalable theoretical framework via an efficient
computer-algebra implementation, which allows the calculation of equations of
motion of retained cumulants up to order $n_{\rm trunc}$ (here, second-order)
for nonlinear quantum systems comprising arbitrary numbers of bosonic modes
under continuous measurement. Supplementary benchmarking simulations for
systems with more than one mode are included in Appendix F.

In the remainder of this paper, we employ (S)TEOMs to simulate dynamics of
measurement chains comprising measured quantum systems coupled to multimode
nonlinear QRCs. We find that such QRCs are able to successfully perform the
quantum state classification tasks we have set out as objectives in the
introduction, while operating within the aforementioned constraints.

##  IV Cavity measurement: classifying qubit pointer states

We are now in a position to apply our quantum reservoir computing framework to
the task of quantum state classification, beginning with Task I as identified
in Fig. 1(c): classifying Gaussian states with mean values that differ in
amplitude or phase. A useful example of this task is the classification of
pointer states generated during cavity-mediated dispersive-qubit measurements
[48, 63], when qubit evolution during measurement is negligible.
Classification for nontrivial qubit evolution during measurement has been
considered using classical reservoirs in Ref. [30]. We begin with a
description of all elements of the measurement chain included in Eq. (1) for
this task, depicted schematically in Fig. 5(a).

###  IV.1 Hardware measurement chain

####  IV.1.1 Single cavity mode and nonreciprocal coupling to QRC

Our measured quantum system is described by the Liouvillian superoperator
$\mathcal{L}_{\rm S}$ defined as

$\displaystyle\mathcal{L}_{\rm
S}\hat{\rho}^{c}=-i[-\Delta^{(\sigma)}_{a}\hat{a}_{1}^{\dagger}\hat{a}_{1}+\eta(\hat{a}_{1}+\hat{a}_{1}^{\dagger}),\hat{\rho}^{c}]+\kappa\mathcal{D}[\hat{a}_{1}]\hat{\rho}^{c}$
(19)

which defines a single cavity mode $\hat{a}_{1}$, with frequency
$\omega_{a_{1}}$, damping rate $\kappa$, and coherently driven at frequency
$\omega_{d}$ with drive strength $\eta$. The superoperator is written in the
frame rotating at the drive frequency, rendering the drive term time-
independent and introducing the cavity mode detunings

$\displaystyle\Delta^{(\sigma)}_{a}=\omega_{d}-(\omega_{a_{1}}+\delta^{(\sigma)}).$
(20)

Here frequency shifts $\delta^{(\sigma)}$ indexed by $\sigma$ model the case
of dispersive qubit readout: for the measurement of two qubits
$\hat{\sigma}_{1,2}^{z}$ dispersively coupled to the cavity with strengths
$\chi_{1,2}$ respectively, the dispersive shifts corresponding to the four
joint qubit $\hat{z}$-basis eigenstates are given by
$\delta^{(\sigma)}\in\\{\chi_{1}-\chi_{2},\chi_{1}+\chi_{2},-\chi_{1}-\chi_{2},-\chi_{1}+\chi_{2}\\}$.
Under coherent driving, the cavity state evolves to a distinct amplitude and
phase for each dispersive shift, serving as a pointer to the corresponding
qubit state. The computational task set for the QRC is then to distinguish
these pointer states. For concreteness, we choose the drive to be resonant
with the bare measured cavity frequency, $\omega_{d}=\omega_{a}$, and choose
cavity detunings that map to effective dispersive shifts
$\chi_{1}=2\kappa,\chi_{2}=0.5\kappa$.

The description of this quantum system alone is insufficient to define its
conditional dynamics when embedded within a measurement chain, which includes
the measured QRC and a coupling element. The specific form of the coupling is
important as well: ideally, we require the information from the measured
quantum system to be seen by the QRC, but for the QRC to not strongly couple
to, and thus renormalize the properties of, the measured quantum system.

These requirements can be met by ensuring that $\mathcal{L}_{c}$ describes a
non-reciprocal coupling between the measured quantum system and the QRC. Here
we couple the single cavity mode $\hat{a}_{1}$ to a single node $\hat{b}_{1}$
of the QRC via a directional amplifier interaction, realized by balancing a
coherent QND interaction with its appropriately chosen dissipative counterpart
[64]. This is defined by the specific coupling superoperator $\mathcal{L}_{c}$
defined as

$\displaystyle\mathcal{L}_{c}\hat{\rho}=-i[-\text{\it{g}}_{c}\hat{P}_{1}\hat{X}_{a},\hat{\rho}]+\Gamma_{c}\mathcal{D}[\hat{X}_{a}+i\hat{P}_{1}]\hat{\rho}$
(21)

where $\hat{X}_{a}=\frac{1}{\sqrt{2}}(\hat{a}_{1}+\hat{a}_{1}^{\dagger})$,
$\hat{P}_{a}=-\frac{i}{\sqrt{2}}(\hat{a}_{1}-\hat{a}_{1}^{\dagger})$ are the
canonically conjugate quadratures of the measured cavity mode, and
$\hat{X}_{1},\hat{P}_{1}$ are analogously-defined quadratures of the QRC
$\hat{b}_{1}$ mode. The form of the interaction, realizable in cQED using
multiple parametric pumps [65, 66], ensures that the quadratures
$\hat{P}_{a}$, $\hat{X}_{1}$ undergo coupled evolution, as defined by the
equations of motion

$\displaystyle\frac{d}{dt}\langle\hat{P}_{a}\rangle^{c}\\!$
$\displaystyle=\\!{\rm tr}\\{[\mathcal{L}_{\rm S}+\mathcal{S}_{\rm
meas}(dW)]\hat{\rho}^{c}\hat{P}_{a}\\}-\left[\text{\it{g}}_{c}-\Gamma_{c}\right]\langle\hat{P}_{1}\rangle^{c},$
(22a) $\displaystyle\frac{d}{dt}\langle\hat{X}_{1}\rangle^{c}\\!$
$\displaystyle=\\!{\rm tr}\\{[\mathcal{L}_{\rm QRC}+\mathcal{S}_{\rm
meas}(dW)]\hat{\rho}^{c}\hat{X}_{1}\\}-\left[\text{\it{g}}_{c}+\Gamma_{c}\right]\langle\hat{X}_{a}\rangle^{c}\\!\\!.$
(22b)

The first term in each equation describes the uncoupled dynamics of the
corresponding system, while the second term includes the conditional evolution
due to the measurement superoperator. The third term is the coupling
contribution: under the specific choice of coupling strength and phase

$\displaystyle\text{\it{g}}_{c}=\Gamma_{c}$ (23)

it is clear from Eq. (22b) that $\mathcal{L}_{c}$ realizes a directional
amplifier: the cavity $\hat{X}_{a}$ quadrature is amplified by a factor
$\propto\Gamma_{c}$ and impinges on the QRC, while the cavity itself is not
driven by the QRC in the opposite direction.

Even if the cavity-QRC coupling is nonreciprocal, measurement of the QRC nodes
downstream leads to measurement-conditioned evolution of the cavity
$\propto\mathcal{S}_{\rm meas}(dW)$ [67]. This is evident from the conditional
equation of motion for the cavity $\langle\hat{X}_{a}\rangle^{c}$ quadrature
(which drives the QRC),

$\displaystyle\frac{d}{dt}$ $\displaystyle\langle\hat{X}_{a}\rangle^{c}={\rm
tr}\\{\mathcal{L}_{\rm S}\hat{\rho}^{c}\hat{X}_{a}\\}$ $\displaystyle+$
$\displaystyle\sum_{k}\frac{\sqrt{\gamma_{k}}}{2}\left[C^{c}_{ab_{k}}+(C^{c}_{ab_{k}})^{*}+C^{c}_{a^{\dagger}b_{k}}+(C^{c}_{a^{\dagger}b_{k}})^{*}\right]dW^{X}_{k}(t)$
$\displaystyle-i$
$\displaystyle\sum_{k}\frac{\sqrt{\gamma_{k}}}{2}\left[C^{c}_{ab_{k}}-(C^{c}_{ab_{k}})^{*}+C^{c}_{a^{\dagger}b_{k}}-(C^{c}_{a^{\dagger}b_{k}})^{*}\right]dW^{P}_{k}(t)$
(24)

where the noted stochastic terms are visible in the second and third lines.
These terms depend on joint cumulants of the cavity-QRC quantum state and are
thus influenced by the QRC evolution; in particular, they can be nonzero even
if Eq. (23) is satisfied. This principle of measurement-conditioned evolution
has been analyzed in cQED experiments, for example in quantum trajectories of
dispersively measured qubits nonreciprocally coupled to a quantum-limited
phase-sensitive amplifier [68], and is accurately captured by the joint
description of the measurement chain we employ here. In the present case,
typical conditional cavity dynamics are illustrated for
$\langle\hat{X}_{a}\rangle^{c}$ in Fig. 5(b), color-coded for the four states
defined by $\Delta_{a}^{(\sigma)}$.

####  IV.1.2 QRC

Our choice of coupling leads to the cavity $\hat{X}_{a}$ quadrature driving
the QRC via its $\hat{X}_{1}$ quadrature, whose resulting conditional
evolution is given by:

$\displaystyle\frac{d}{dt}\langle\hat{X}_{1}\rangle^{c}$ $\displaystyle={\rm
tr}\\{\mathcal{L}_{\rm
QRC}\hat{\rho}^{c}\hat{X}_{1}\\}-2\Gamma_{c}\langle\hat{X}_{a}\rangle^{c}$
$\displaystyle+$
$\displaystyle\sum_{k}\frac{\sqrt{\gamma_{k}}}{2}\left[C^{c}_{b_{1}b_{k}}+(C^{c}_{b_{1}b_{k}})^{*}+2C^{c}_{b_{1}^{\dagger}b_{k}}\right]dW^{X}_{k}(t)$
$\displaystyle-i$
$\displaystyle\sum_{k}\frac{\sqrt{\gamma_{k}}}{2}\left[C^{c}_{b_{1}b_{k}}-(C^{c}_{b_{1}b_{k}})^{*}\right]dW^{P}_{k}(t)$
(25)

The QRC nodes thus evolve under the cavity signal, but also the QRC’s
intrinsic network dynamics governed by $\mathcal{L}_{\rm QRC}$, and
conditional evolution that depends on cumulants of the joint QRC state.
Typical conditional dynamics of the $\langle\hat{X}_{1}\rangle^{c}$ quadrature
of QRC node $\hat{b}_{1}$ for each generated cavity state is shown in Fig.
5(c). For this illustration, we have chosen two-node QRC parameters:
$\gamma_{1}=\gamma_{2}=\kappa$, $\Lambda_{1}=\Lambda_{2}=0.005\gamma_{1}$,
$\Delta_{1}=\Delta_{2}=0.0$, and $\text{\it{g}}_{12}=\gamma_{1}$.

Nonlinear processing of the cavity signal is encoded in this time evolution of
the multimode QRC state. However, processed information about this state must
be extracted from noisy measurement records $J_{k}^{X,P}(t)$ defined in Eqs.
(4), (5) obtained from the measurement chain, which we discuss next.

Figure 5: (a) Measurement chain using a 2-node QRC for the classification of
states of a single cavity mode. The cavity mode is coupled to the QRC using a
quantum-limited directional amplifier. Each node of the QRC is continuously
monitored using separate output channels via heterodyne measurement.
Conditional dynamics of (b) the cavity quadrature
$\langle\hat{X}_{a}\rangle^{c}$, (c) the QRC node quadrature
$\langle\hat{X}_{1}\rangle^{c}$ it drives, and (d) measurement records
$J_{1}^{X}(t)$, are shown for one initialization of the measurement chain per
state, with colors blue, orange, green, and red corresponding to cavity mode
detunings $\Delta_{a}^{(\sigma)}\in\\{2.5,1.5,-1.5,-2.5\\}\kappa$
respectively. (e) Dynamics of QRC outputs projected onto the measured subspace
$\\{I_{1}^{\varphi},I_{2}^{\varphi}\\}$ as a function of time, with panels (1)
through (4) corresponding to $\kappa t\in[1.0,2.0,4.0,8.0]$ respectively.
$100$ trajectories are shown for each state $\sigma$. (f) Classification
accuracy as a function of time, computed for a test set of size $Q_{\rm
test}=200$ measured outputs per state. Black and gray curves correspond to a
typical nonlinear QRC and a linear QRC respectively.

###  IV.2 Training and classification using finitely-sampled measurement
distributions

Measurement records $\\{J_{k}^{X,P}(t)\\}$ contain experimentally-accessible
information about the multimode QRC state for a given initialization of the
measurement chain. We consider the case of single-shot reservoir readout,
where the output $\mathbf{x}(t)$ of the QRC is constructed as the vector of
all measured quadratures obtained via linear boxcar filtering of single-shot
heterodyne measurement records [69],

$\displaystyle\mathbf{x}(t)\equiv\begin{pmatrix}I_{1}^{X}(t)\\\
I_{1}^{P}(t)\\\ I_{2}^{X}(t)\\\ I_{2}^{P}(t)\end{pmatrix},\leavevmode\nobreak\
I_{k}^{X,P}(t)=\frac{1}{t}\int_{0}^{t}d\tau\leavevmode\nobreak\
J_{k}^{X,P}(\tau).$ (26)

Measured quadratures are themselves stochastic, so that the QRC output
provides a distinct stochastic trajectory for repeated runs of the measurement
chain; QRC output corresponding to $\Delta_{a}^{(\sigma)}$ is thus labelled
$\mathbf{x}^{(\sigma)}_{(q)}$, indexed by both $\sigma$ and a trajectory label
$q$.

Having defined each element of the measurement chain and the measured QRC
output, we can finally describe its operation for reservoir computing. We
employ an output layer configured to yield a predicted class label
$\sigma^{p}$ using the measured QRC output,
$\sigma^{p}=f_{N}\\{\mathbf{W}_{\rm
O}\mathbf{x}^{(\sigma)}_{(q)}+\mathbf{b}\\}$, where $\mathbf{W}_{\rm O}$
defines a matrix of time-independent trainable weights, $\mathbf{b}$ a vector
of trainable biases, and where $f_{N}\\{\cdot\\}$ is a fixed (that is,
untrained) normalizing function. Perfect classification yields
$\sigma^{p}=\sigma$ for all $\sigma=1,\ldots,C$, irrespective of $q$.

Formally, training is performed by minimizing a cost function, here the
squared error over a training set whose output is known, and importantly is a
convex minimization problem that is guaranteed to converge; full details of
this training are included in Appendix J. Here we instead present a visual
representation of training and classification, shedding light into how the QRC
dynamics enable classification. Classification using linear weights
$\mathbf{W}_{\rm O}$ is equivalent to constructing linear decision boundaries
to separate reservoir outputs $\mathbf{x}^{(\sigma)}_{(q)}$ for different
$\sigma$ in the four-dimensional measured phase space spanned by
$\mathbf{x}=(I_{1}^{X},I_{1}^{P},I_{2}^{X},I_{2}^{P})$, as shown in Appendix
J. While this four-dimensional space is difficult to visualize, we find that
for the present task two optimally chosen quadratures prove to be sufficient
for successful classification. We thus rewrite the matrix of trainable weights
as $\mathbf{W}_{\rm O}\mathbf{x}=\widetilde{\mathbf{W}}_{\rm
O}\mathcal{M}^{\varphi}\mathbf{x}$, where $\mathcal{M}^{\varphi}$ acts to
project the measured QRC output onto a two-dimensional measured subspace
$\\{I_{1}^{\varphi},I_{2}^{\varphi}\\}$,

$\displaystyle\begin{pmatrix}I_{1}^{\varphi}\\\
I_{2}^{\varphi}\end{pmatrix}\equiv\mathcal{M}^{\varphi}\mathbf{x}\\!=\\!\begin{pmatrix}\cos\varphi_{1}&\sin\varphi_{1}&0&0\\\
0&0&\cos\varphi_{2}&\sin\varphi_{2}\end{pmatrix}\\!\\!\\!\begin{pmatrix}I_{1}^{X}\\\
I_{1}^{P}\\\ I_{2}^{X}\\\ I_{2}^{P}\end{pmatrix}.$ (27)

We note that all QRC node quadratures are still measured in each measurement
run;. The projection phases are simply applied as part of the output layer,
and are learned during training.

In Fig. 5(d), we project QRC outputs $\\{\mathbf{x}^{(\sigma)}_{(q)}\\}$ onto
the measured subspace to visualize their evolution. Panels (1) through (4)
show these outputs as a function of time, $\kappa t\in[1.0,2.0,4.0,8.0]$
respectively. We see clearly that for short times (1), all four states are
indistinguishable and the measured quadrature distributions overlap. This is
not surprising, as the measurement chain is initialized to vacuum for all
states. As the cavity is populated, distributions for states with distinct
phases first separate (2) while the QRC response is still effectively linear.
By (3), the cavity signal corresponding to higher amplitude states (orange,
green) has increased sufficiently to drive the QRC into a nonlinear regime.
This is clear from the rotation of these distributions relative to those
corresponding to the low amplitude states, as well as their visible
distortion; the latter is a signature of amplified and squeezed QRC node
fluctuations due to the Kerr nonlinearity, also seen in Fig. 3(c). The effect
of the nonlinearity clearly separates all four distributions into
(correspondingly shaded) distinct phase space regions.

It is then clear that linear decision boundaries (dashed black) can be drawn
in the measured subspace to separate distributions for all $C$ states.
Training using a set of QRC output trajectories
$\\{\mathbf{x}^{(\sigma)}_{(q)}\\}$ for $Q_{\rm Train}=100$ trajectories per
state ($q=1,\ldots,Q_{\rm Train}$) simply learns the parameters
$\\{\widetilde{\mathbf{W}}_{\rm O},\mathbf{\varphi}\\}$ uniquely specifying
the optimal decision boundaries. Note that these decision boundaries, and
therefore the training protocol to learn them, must take into account both the
mean and variance of the finitely-sampled QRC output distributions. These
quantities are affected by factors such as measurement time and filtering;
therefore, training on measured output distributions is tailored to the finite
experimental resources available for output processing. This is in sharp
contrast to training on unconditional expectation values, which in principle
assumes access to infinitely-many measurement records (for further discussion,
see Appendix K). Finally, we emphasize that the linear output layer employed
here serves two important, yet simple purposes: (i) filtering out measurement
noise obscuring the underlying QRC dynamics, and (ii) drawing linear decision
boundaries that optimally separate regions of the measured phase space into
classes. The nontrivial processing step that enables classification - the
evolution of measured outputs into linearly-separable regions of the measured
subspace - is performed by the nonlinear QRC.

Once the QRC is trained, we can quantify its performance by predicting state
labels $\sigma^{p}$ for a distinct test set of reservoir outputs
$\\{\mathbf{x}^{(\sigma)}_{(q)}\\}$ for $q=1,\ldots,Q_{\rm Test}=200$. Using
the trained output layer, we determine $\sigma^{p}$ for each
$\mathbf{x}^{(\sigma)}_{(q)}$ and calculate the classification accuracy as a
function of time, defined as the fraction of $C\times Q_{\rm Test}$ total
trajectories that were correctly classified. The classification accuracy is
plotted in Fig. 5(f), and is characterized by two metrics: the maximum
classification accuracy $\mathcal{C}_{\rm max}$ achieved within $\kappa
t=10.0$, and the time $t_{\rm max}$ taken to reach this maximum. We note that
the specific QRC shown here achieves perfect classification ($\mathcal{C}_{\rm
max}=1$) before the measurement chain dynamics have settled into steady-state
[see Figs. 5(b), (c)], indicating the ability of the QRC to process the
information encoded in time-evolving quantum signals.

Also shown in gray is the classification accuracy for a linear QRC,
$\bar{\Lambda}=0$, which saturates to $0.5$. The typical dynamics in the
measured phase space are analogous to that shown in panel (2) of Fig. 5(e);
the linear QRC is able to separate the pairs of states with different phases
(blue-orange from red-green) but not distinguish states within this subset.
Having used pecisely the same form of output layer for both linear and
nonlinear QRCs, we are able to unequivocally identify the QRC nonlinearity as
being critical for the present classification task.

It is important to note, however, that a completely arbitrary nonlinear QRC
will not succeed at this classification task; internal QRC parameters
introduced in Eq. (2) must satisfy some task-dependent constraints, which we
explore next.

###  IV.3 Classification performance as a function of hyperparameters

Hyperparameters.$-$ For large classical reservoirs, it is conventional - and
practical - to describe reservoir parameters in terms of average
hyperparameters instead of individual microscopic values. Here we adopt a
similar convention despite the small size of the two-node QRC, parameterizing
it as

$\displaystyle\Lambda_{k}$
$\displaystyle=\bar{\Lambda}+[-\epsilon,\epsilon]\cdot\bar{\Lambda},\leavevmode\nobreak\
\gamma_{k}=\bar{\gamma}=\kappa\leavevmode\nobreak\ \forall\leavevmode\nobreak\
k,$ $\displaystyle\Delta_{k}$
$\displaystyle=\bar{\Delta}+[-\epsilon,\epsilon]\cdot\bar{\gamma},\leavevmode\nobreak\
\text{\it{g}}_{12}=\bar{\text{\it{g}}}$ (28)

where $[x,y]$ denotes random samples from a uniform distribution with
probability density between $x$ and $y$; we set $\epsilon=0.1$ to observe the
effect of random variations in QRC nonlinearities and frequencies that may
occur in experiment. For convenience, the QRC damping rates are all set equal
to the measured quantum system decay rate $\kappa$.

We analyze QRC classification performance as a function of the above
hyperparameters in terms of the maximum classification accuracy
$\mathcal{C}_{\rm max}$, and the time taken to reach this maximum, $t_{\rm
max}$. In doing so, we ensure that all other aspects of the measurement chain,
training of the output layer, and testing remain unchanged from Sec. IV.2, so
that the particular role of QRC hyperparameters can be highlighted. Our key
findings are recounted in a summary at the end of this section.

####  IV.3.1 Dependence on nonlinearity $\bar{\Lambda}$

We begin by simulating QRCs with different values of the nonlinearity
hyperparameter $\bar{\Lambda}$, while fixing $\bar{\text{\it{g}}}/\gamma=1.0$,
$\bar{\Delta}/\bar{\gamma}=0.0$, and calculating the maximum classification
accuracy $\mathcal{C}_{\rm max}$. The results are shown in Fig. 6(a). Each red
dot corresponds to a different random QRC, as defined by Eq. (28). The solid
black curve is the mean over these QRCs, while the gray shaded region marks
the range between worst and best performing QRCs. 5 random QRCs are simulated
for each set of hyperparameters, amounting to over 200 distinct QRCs shown in
Fig. 6(a).

For weak nonlinearities, $\mathcal{C}_{\rm max}$ approaches $0.5$, the limit
of a linear QRC as seen earlier. With increasing nonlinearity however, the
QRCs are able to perform the classification task, and $\mathcal{C}_{\rm max}$
approaches unity. The time $t_{\rm max}$ to reach $\mathcal{C}_{\rm max}$ can
further distinguish operating hyperparameters: we plot $\kappa t_{\rm max}$ in
Fig. 6(a) (right-hand axis) in green, for parameters where $\mathcal{C}_{\rm
max}\geq 0.99$ (yellow shaded region); green dots again indicate different
random QRCs with the same hyperparameter values, while the solid green curve
indicates the average, and the (small) green shaded region marks the range
between the worst and best performing QRCs. We see that there exists a subset
of optimal nonlinearity strengths for which the QRC correctly classifies all
states in the shortest amount of time.

This calculation of $\mathcal{C}_{\rm max}$ and $t_{\rm max}$ obtained from
simulation of the complete measurement chain captures the full complexity of
QRC evolution under time-dependent cavity signals with distinct amplitudes and
phases. However, we find that a heuristic understanding of QRC performance
characteristics can be developed by focusing on its steady-state response to a
time-independent drive signal. For the pointer state classification task, we
can define an effective drive strength $\eta_{\rm eff}^{(\sigma)}$ as the
long-time limit of the cavity signal corresponding to $\Delta_{a}^{(\sigma)}$
incident on the QRC,

$\displaystyle\eta_{\rm eff}^{(\sigma)}$
$\displaystyle=-2\Gamma_{c}\\!\lim_{t\to\infty}\\!\\{\sqrt{2}\langle\hat{X}_{a}^{(\sigma)}(t)\rangle\\}=-2\Gamma_{c}\eta\left[\\!\frac{\Delta_{a}^{(\sigma)}}{(\Delta_{a}^{(\sigma)})^{2}+\frac{\kappa^{2}}{4}}\\!\right]\\!,$
(29)

where the second equality is obtained from the steady-state solution of the
unconditional version of Eq. (24). As visible in cavity signals in Fig. 5(b),
for $|\Delta_{a}^{(\sigma)}|\geq\kappa$, $\sigma=1,4$ (blue, red) correspond
to larger $|\Delta_{a}^{(\sigma)}|$ and yield weaker effective drives than
$\sigma=2,3$ (orange, green).

We now analyze the steady-state QRC response in a regime enabling perfect
classification, $\bar{\Lambda}/\bar{\gamma}=0.005$ in Fig. 6(a) (and setting
$\epsilon=0$). To do so, we plot the steady-state values of
$|\langle\hat{b}_{1,2}\rangle|$ (top panel) and ${\rm
arg}\langle\hat{b}_{1,2}\rangle$ (bottom panel) under a constant drive
strength $\eta_{\rm eff}$ in Fig. 6(b), with the particular effective drives
$\eta_{\rm eff}^{(\sigma)}$ indicated by correspondingly colored lines.
Clearly, the QRC amplitude and phase response to higher amplitude cavity
states $\sigma=2,3$ (orange, green) is non-linearly related to its response to
lower amplitude states $\sigma=1,4$ (blue, red). This leads to the nonlinear
amplitude response and phase rotation observed in the measured quadrature
distributions for higher amplitude states in Fig. 5(e), which enables
successful classification of the pointer states.

Figure 6: (a) $\mathcal{C}_{\rm max}$ (left-hand axis) vs. nonlinearity
hyperparameter $\bar{\Lambda}$, for fixed coupling
$\bar{\text{\it{g}}}/\bar{\gamma}=1.0$ and detuning
$\bar{\Delta}/\bar{\gamma}=0.0$. For $\mathcal{C}_{\rm max}\geq 0.99$ (yellow
shaded region), $\kappa t_{\rm max}$ is plotted along the right-hand axis
(green). Dots indicate different random QRCs, solid curves are averages over
these QRCs, and shaded regions mark the range between worst and best
performing QRCs (see text for more details). (b) Steady-state response of
$|\langle\hat{b}_{1,2}\rangle|$, ${\rm arg}\langle\hat{b}_{1,2}\rangle$ as a
function of constant drive strength $\eta_{\rm eff}$. (c) Classical phase
diagram of the two-node QRC in $(\bar{\text{\it{g}}},\mathcal{N}_{\rm eff})$
space, for fixed detuning $\bar{\Delta}/\bar{\gamma}=0.0$. (d) Steady-state
response of $|\langle\hat{b}_{1,2}\rangle|$ as a function of increasing
nonlinearity $\bar{\Lambda}$ across the vertical dashed line
($\bar{\text{\it{g}}}/\bar{\gamma}=1.0$) in (c).

We find that the nonlinear QRC response observed in Fig. 6(b) is typical when
the two-node QRC is successful at the pointer state classification task; it is
therefore helpful to locate such regimes in terms of general QRC
hyperparameters. The classical phase diagram of the two-node QRC ends up
providing the pathway to identifying such regimes. Furthermore, this phase
diagram depends on the effective QRC drive and nonlinearity strengths only via
the effective parameter (for details see Appendix G)

$\displaystyle\mathcal{N}_{\rm eff}=\frac{\eta_{\rm
eff}}{\bar{\gamma}}\sqrt{\frac{\bar{\Lambda}}{\bar{\gamma}}}.$ (30)

To construct the QRC phase diagram informed by the pointer state
classification task, we choose $\eta_{\rm eff}\equiv{\rm
max}_{\sigma}\\{|\eta_{\rm eff}^{(\sigma)}|\\}$, corresponding to higher
amplitude cavity states that drive the QRC into nonlinear processing regimes.

We plot the resulting phase diagram in Fig. 6(c) as a function of
$\mathcal{N}_{\rm eff}$ and coupling $\bar{\text{\it{g}}}$, for fixed
$\bar{\Delta}/\bar{\gamma}=0.0$. The orange region indicates the classical
bistability, analogous to that found for the single QRC node. The vertical
dashed line indicates $\bar{\text{\it{g}}}/\bar{\gamma}=1.0$, namely the
cross-section along which Fig. 6(a) is plotted. Also shown in Fig. 6(d) is the
steady-state response of the QRC node amplitudes and phases along this cross-
section; the solid blue line in (c) and yellow shaded region in (d) mark
$\mathcal{N}_{\rm eff}$ values for which $\mathcal{C}_{\rm max}\geq 0.99$. For
weak nonlinearities where classification is unsuccessful, the mode amplitudes
and phases are effectively constant, indicating the linear regime of operation
where QRC states generated by different $\eta_{\rm eff}$ cannot be separated.
For increasing $\bar{\Lambda}$, approaching the classical bistability results
in a nonlinear QRC response enabling successful classification. When the
nonlinearity is further increased however, $|\langle\hat{b}_{1,2}\rangle|$
begins to decrease; this is a result of the increasing Kerr shift moving the
node resonances further from the steady-state input signals. The reduced QRC
amplitude is responsible for the increase in $t_{\rm max}$ seen in Fig. 6(a).

####  IV.3.2 Dependence on coupling $\bar{\text{\it{g}}}$ and detuning
$\bar{\Delta}$

We can similarly analyze the performance of QRCs as a function of coupling and
detuning hyperparameters, for fixed nonlinearity (here,
$\bar{\Lambda}/\bar{\gamma}=0.002$, which fixes $\mathcal{N}_{\rm eff}=0.8$).
In cQED realizations of the proposed QRCs, these hyperparameters can be
rendered in-situ tunable (via flux-modulated SQUIDs and parametric couplers
respectively), thus providing experimental control over QRC performance.

We begin by calculating $\mathcal{C}_{\rm max}$ as a function of the coupling
hyperparameter $\bar{\text{\it{g}}}$ in Fig. 7(a), fixing the detuning
$\bar{\Delta}/\bar{\gamma}=0.0$. This cross-section of
$(\bar{\Delta},\bar{\text{\it{g}}})$ space is indicated in the center panel
classical phase diagram with a vertical dashed line, with the corresponding
QRC steady-state amplitude and phase response plotted in the left panel of
Fig. 7(c). For weak couplings, the $\hat{b}_{2}$ is correspondingly weakly
excited, since the input couples only to $\hat{b}_{1}$. Only the $\hat{b}_{1}$
is effectively available for processing; this is insufficient for the task at
hand and the QRC does not attain $\mathcal{C}_{\rm max}=1$ within the
measurement time.

With increasing coupling, the QRC achieves perfect classification;
furthermore, an optimal range of couplings exists where $t_{\rm max}$ (green,
right-hand axis) is shortest. From the left panel of Fig. 7(c), we see that
here the QRC response for both nodes is strongest. Curiously, the performance
dips again for strong couplings: here the dynamics due to the linear coupling
dominate over the nonlinearity, as the QRC polariton modes are further split
and respond weakly to the signal emanating from the cavity, resulting in the
decreased response visible in the steady state plots.

Figure 7: $\mathcal{C}_{\rm max}$ (left-hand axis) for fixed nonlinearity
$\bar{\Lambda}/\bar{\gamma}=0.002$ (and hence fixed $\mathcal{N}_{\rm
eff}=0.8$) vs. (a) coupling hyperparameter $\bar{\text{\it{g}}}$ for fixed
detuning $\bar{\Delta}/\bar{\gamma}=0.0$, and vs. (b) average detuning
hyperparameter $\bar{\Delta}$ for fixed coupling
$\bar{\text{\it{g}}}/\bar{\gamma}=2.0$. For both plots, where
$\mathcal{C}_{\rm max}\geq 0.99$ (yellow shaded region), $t_{\rm max}$ is
plotted along the right-hand axis (green). Dots indicate different random
QRCs, solid curves are averages over these QRCs, and shaded regions mark the
range between worst and best performing QRCs (see text for more details). (c)
Classical phase diagram of the two-node QRC in
$(\bar{\Delta},\bar{\text{\it{g}}})$ space for fixed $\mathcal{N}_{\rm
eff}=0.8$. Cross-sections of the steady-state response of QRC nodes are shown
in (d) for fixed detuning $\bar{\Delta}/\bar{\gamma}=0.0$ as a a function of
coupling $\bar{\text{\it{g}}}$, and in (e) for fixed coupling
$\bar{\text{\it{g}}}/\bar{\gamma}=2.0$ as a a function of detuning
$\bar{\Delta}$.

To demonstrate the advantage of in-situ control, we now consider QRCs with
sub-optimal coupling $\bar{\text{\it{g}}}/\bar{\gamma}=2.0$ (where
$\mathcal{C}_{\rm max}\simeq 0.8$, see Fig. 7(a)). Figure 7(b) shows the
resulting $\mathcal{C}_{\rm max}$ as a function of the detuning hyperparameter
$\bar{\Delta}$, a range indicated via the horizontal dashed line in the phase
of diagram of Fig. 7(c). The corresponding steady-state QRC amplitude and
phase responses are shown in the right panel. We clearly see that
$\mathcal{C}_{\rm max}$ approaches unity once more for specific positive and
negative detunings where the QRC response is enhanced.

####  IV.3.3 Dependence on damping $\bar{\gamma}$

Thus far, the QRC node damping rate has appeared as the natural scaling
parameter for hyperparameters in our analysis. Hence if $\bar{\gamma}$ is
varied but the remaining QRC hyperparameters are suitably adjusted so that
$\\{\mathcal{N}_{\rm
eff},\bar{\text{\it{g}}}/\bar{\gamma},\bar{\Delta}/\bar{\gamma}\\}$ are
unchanged, the QRC will operate in a fixed regime of the classical phase
diagram; provided this regime enables nonlinear processing as previously
analyzed, successful classification will be possible.

However, two aspects of QRC dynamics controlled by $\bar{\gamma}$ are not
captured by its steady-state phase diagram: the response time of QRC nodes to
time-dependent input signals, as well as the magnitude of QRC response for
fixed input signal strengths. To isolate the influence of $\bar{\gamma}$ on
QRC performance via these specific factors, we must ensure that the QRC
operates in a fixed nonlinear regime of the classical phase diagram even as
$\bar{\gamma}$ is varied.

This analysis is carried out in Appendix H, with our main results summarized
here. We find that for $\bar{\gamma}\ll\kappa$, classification is still
possible in nonlinear QRC regimes; however, as the QRC responds much slower
than the evolution timescale of cavity signals (set by $\kappa$), $t_{\rm
max}$ for classification increases. As $\bar{\gamma}$ increases the QRC
response time is reduced, leading to a decrease in $t_{\rm max}$ until a
minimum is attained around $\bar{\gamma}\simeq 2\kappa$. For larger
$\bar{\gamma}$, while the QRC response time is still fast relative to the
quantum system evolution timescale, the QRC node amplitude is suppressed such
that the amplitude of extracted output signals from the measurement chain is
reduced relative to the measurement noise, thereby increasing $t_{\rm max}$
once more.

####  IV.3.4 Summary

The results of this section indicate that a range of hyperparameter values
$\\{\mathcal{N}_{\rm
eff},\bar{\text{\it{g}}}/\bar{\gamma},\bar{\Delta}/\bar{\gamma}\\}$ exist for
which the QRC is able to perform nonlinear processing of input signals, and
consequently succeed at the task of pointer state classification considered
here. Importantly, regions of classical bistability - key signatures of
classical nonlinearity - are strong (although not exclusive) indicators of
such parameter regimes. The mechanism of classification is the nonlinear
response of QRC node amplitudes with respect to the amplitude and phase of
incident signals from the cavity mode.

###  IV.4 Dependence on input signal strength: classical vs. quantum regimes
of operation

Our analysis of QRC performance as a function of hyperparameters highlights
the importance of operating in dynamical regimes that enable nonlinear
processing of input signals. Importantly, the effective nonlinearity parameter
$\mathcal{N}_{\rm eff}$ that is central to such regimes depends not only on
the QRC nonlinearity $\bar{\Lambda}$, but also on the effective drive
$\eta_{\rm eff}$. The latter is ultimately related to the strength of the
drive $\eta$ incident on the measurement chain, which in turn can be
constrained by various physical considerations. For composite measurement
chains considered here, the incident field initially drives the measured
quantum system upstream of the QRC. Hence, its amplitude is dictated by the
desired response of this quantum system, namely the amplitude required to
generate quantum states of interest, as well as to avoid possible non-
idealities that may appear under strong driving. Secondly, to decrease the
energy of physical computation, one must naturally consider a reduction in the
input power to the measurement chain containing the QRC.

With this view, we consider separate instances of pointer state classification
tasks for the same values of detunings $\Delta_{a}^{(\sigma)}$ defining the
pointer states as considered earlier, but for different values of drive input
to the measurement chain, $\eta/\bar{\gamma}\in\\{26.0,15.0,9.0\\}$. These
input drive values translate to effective QRC drives $\eta_{\rm
eff}/\bar{\gamma}\in\\{31.2,18.0,10.8\\}$. We then determine classification
metrics $\mathcal{C}_{\rm max}$ and $\kappa t_{\rm max}$ as a function of QRC
nonlinearity strength $\bar{\Lambda}$ as before, holding the remaining QRC
hyperparameters fixed, $\bar{\text{\it{g}}}/\bar{\gamma}=1.0$,
$\bar{\Delta}/\bar{\gamma}=0.0$.

Figure 8: $\mathcal{C}_{\rm max}$ (top panel) and $t_{\rm max}$ (bottom panel)
as a function of average nonlinearity hyperparameter $\bar{\Lambda}$, for
instances of the pointer state classification task corresponding to three
different effective QRC drives $\eta_{\rm eff}/\bar{\gamma}$. Dots indicate
different random QRCs, solid curves are averages over these QRCs, and shaded
regions mark the range between worst and best performing QRCs. Crosses
indicate QRCs with fixed $\mathcal{N}_{\rm eff}\simeq 1.4$ for the three
considered task instances; for weaker $\eta_{\rm eff}/\bar{\gamma}$,
corresponding QRCs have stronger nonlinearity and operate in more quantum
regimes (see text for more details).

The achieved $\mathcal{C}_{\rm max}$ is plotted in the top panel of Fig. 8 for
the aforementioned $\eta_{\rm eff}/\bar{\gamma}$ values (blue, purple, and
peach respectively). Immediately, we note that for weaker input drives, a
larger QRC nonlinearity $\bar{\Lambda}$ is required for successful pointer
state classification. When interpreted in the context of the effective
nonlinearity $\mathcal{N}_{\rm eff}$, this result makes sense: the latter
depends on the product of the nonlinearity $\bar{\Lambda}$ of the QRC, and the
effective drive $\eta_{\rm eff}$ it sees. To operate QRCs in a nonlinear
processing regime given by a fixed value of $\mathcal{N}_{\rm eff}$ (for
example, $\mathcal{N}_{\rm eff}\simeq 1.4$ marked by colored crosses in Fig.
8), a weaker $\eta_{\rm eff}$ will require stronger $\bar{\Lambda}$.

This intuitive result encourages a complementary perspective on reservoir
processing across varying input drive strengths, enabled by our full quantum
treatment of the measurement chain using truncated cumulants. As shown in Sec.
III, increasing $\bar{\Lambda}$ for fixed $\mathcal{N}_{\rm eff}$ controllably
transitions Kerr QRCs from classical to quantum regimes of operation.
Therefore, considering instances of the pointer state classification task
across drive strengths allows us to analyze and compare QRC performance on the
task throughout this transition.

In this context, extrapolating the results of Fig. 8 to weaker nonlinearities
leads to an immediate observation: even a QRC operating in the classical
nonlinear regime ($\bar{\Lambda}\to 0$ for fixed $\mathcal{N}_{\rm eff}$) can
perform the pointer state classification task we consider here, albeit
requiring increasing input strengths. In this limit, the quantum state of the
measurement chain is described as a product of coherent states, defined
entirely by its first-order cumulants. Fortunately, the information that
distinguishes qubit pointer states is also encoded in their first-order
cumulants (amplitudes and phases) and not their second-order cumulants
(quantum fluctuations). As a result, the nonlinear mapping of this information
to measured outputs by the nonlinear QRC is described primarily by dynamics of
first-order cumulants, and is therefore possible in the classical limit of QRC
operation.

Fig. 8 indicates that for weaker incident drive strengths, classification is
enabled by QRCs with stronger nonlinearities, operating in more quantum
regimes; however, not all classification metrics are unchanged. In particular,
the time $t_{\rm max}$ to reach maximum classification accuracy (plotted in in
the lower panel of Fig. 8) increases for weaker input drives. This can be
understood as follows. Decreasing the input drive strength reduces the
amplitude (and hence first-order cumulants) of the pointer states to be
classified. This reduced amplitude can be compensated by increasing the QRC
nonlinearity to perform the desired nonlinear mapping to measured outputs;
however, the amplitude of QRC outputs is still correspondingly reduced. In
terms of the measured subspace $\\{I_{1}^{\varphi},I_{2}^{\varphi}\\}$ shown
in Fig. 5(e), the QRC output distributions to be separated have centroids that
are closer to each other. The spread of these distributions, on the other
hand, is determined by measurement noise (independent of drive strength) and
QRC-mediated amplified or squeezed quantum fluctuations (mostly unchanged for
fixed $\mathcal{N}_{\rm eff}$); it primarily decreases with integration time
$t$ due to filtering. Hence a decrease in input strengths reduces the relative
separation of the distributions for a given measurement time, making them
harder to separate, and demanding longer integration times $t_{\rm max}$ to
enable classification.

The above results have an important consequence for the analysis of quantum
reservoir processing for general computational tasks. Information processing
with reduced input signals necessitates strongly-nonlinear QRCs, and thus
provides a physically-motivated probe of QRC operation in quantum regimes.
However, by using a framework that accounts for the complete measurement
chain, we explicitly observe that the decrease in input power can lead to an
increase in required measurement resources for computation; in this particular
case, the measurement time. Such costs must be accounted for in any analysis
of the possible advantages of reservoir computing in the quantum regime.

Finally, we note that processing with low input signals is commonplace in cQED
experiments. Constraints on measurement time due to qubit decay and the use of
noisy HEMT amplifiers [70] in cQED experiments means that the reduced output
signal is overcome by including quantum-limited amplifiers [71, 72, 73] as
part of the measurement chain. Such amplifiers can be included as part of our
QRC framework for operation with low signal powers, as additional quantum
measurement resources. The key considerations in doing so are discussed in
Appendix L (albeit applied to our second computational task, which is
considered next).

##  V Amplifier measurement: classifying Gaussian states with equal means and
distinct variances

We will now apply our quantum reservoir computing framework to Task II: the
classification of Gaussian states that differ only in their second-order
moments, or variances. To generate such states, we consider a minimal system
based on a two-mode amplifier; as a shorthand, we refer to this task as that
of amplifier state classification from here on.

###  V.1 Hardware measurement chain

The quantum measurement chain is depicted schematically in Fig. 9(a). The
quantum system is coupled to a single-node QRC, which we will see proves
sufficient for the classification task at hand. We again begin with a
description of each element of this measurement chain, depicted in Fig. 9(a).

####  V.1.1 Two-mode amplifier and nonreciprocal coupling to QRC

Our measured quantum system is a general two-mode amplifier, described by the
Liouvillian superoperator $\mathcal{L}_{\rm S}$ defined as

$\displaystyle\mathcal{L}_{\rm
S}\hat{\rho}=\\!-i[\hat{\mathcal{H}}_{0}+\hat{\mathcal{H}}_{\rm
S}^{(\sigma)}+\eta^{(\sigma)}(-i\hat{a}_{1}+i\hat{a}_{1}^{\dagger}),\hat{\rho}]\\!+\\!\sum_{n}\kappa_{n}\mathcal{D}[\hat{a}_{n}]\hat{\rho},$
(31)

where the system Hamiltonians take the form:

$\displaystyle\hat{\mathcal{H}}_{0}$
$\displaystyle=-\sum_{n={1,2}}\Delta_{an}\hat{a}_{n}^{\dagger}\hat{a}_{n},$
$\displaystyle\hat{\mathcal{H}}_{\rm S}^{(\sigma)}$
$\displaystyle=\frac{1}{2}G_{1}^{(\sigma)}e^{-i\pi/2}\hat{a}_{1}^{2}+G_{12}^{(\sigma)}\hat{a}_{1}\hat{a}_{2}+h.c.,$
(32)

where $\Delta_{an}=\omega_{d}-\omega_{an}$. The two-mode system is defined by
a degenerate parametric amplifier interaction $\propto G_{1}$ acting on mode
$\hat{a}_{1}$, and a non-degenerate interaction $\propto G_{12}$. Note that
these interaction strengths can also be set by tunable parametric drives in
cQED implementations. Defining for convenience
$\hat{X}_{a_{1}}=\frac{1}{\sqrt{2}}(\hat{a}_{1}+\hat{a}_{1}^{\dagger})$,
$\hat{P}_{a_{1}}=-\frac{i}{\sqrt{2}}(\hat{a}_{1}-\hat{a}_{1}^{\dagger})$ as
the canonically conjugate quadratures of the $\hat{a}_{1}$ mode, note that our
choice of drive Hamiltonian $\propto\hat{P}_{a_{1}}$ drives the amplifier
$\hat{X}_{a_{1}}$ quadrature only. Distinct states of this system are thus
defined by the drive strength and values of the two interaction parameters,
$(\eta^{(\sigma)},G_{1}^{(\sigma)},G_{12}^{(\sigma)})$.

As before, understanding the dynamics of this two-mode amplifier requires
specifying how it is embedded within the measurement chain including the QRC.
We choose to couple the system $\hat{a}_{1}$ mode to the $\hat{b}_{1}$ node of
the QRC; distinct from Sec. IV, however, we now employ a non-reciprocal
hopping interaction (modeling a standard circulator), defined by the coupling
superoperator $\mathcal{L}_{c}$,

$\displaystyle\mathcal{L}_{c}\hat{\rho}=-i\left[\frac{i}{2}\text{\it{g}}_{c}\hat{a}_{1}^{\dagger}\hat{b}_{1}+h.c.,\hat{\rho}\right]+\Gamma_{c}\mathcal{D}[\hat{a}_{1}+\hat{b}_{1}]\hat{\rho}.$
(33)

To understand the role of this coupling, we can once again analyze the
conditional dynamics of the $\hat{X}_{a_{1}}$ quadrature that couples to the
QRC:

$\displaystyle\frac{d}{dt}$
$\displaystyle\langle\hat{X}_{a_{1}}\rangle^{c}={\rm tr}\\{\mathcal{L}_{\rm
S}\hat{\rho}^{c}\hat{X}_{a_{1}}\\}+\frac{\text{\it{g}}_{c}-\Gamma_{c}}{2}\langle\hat{X}_{1}\rangle^{c}-\frac{\Gamma_{c}}{2}\langle\hat{X}_{a_{1}}\rangle^{c}$
$\displaystyle+$
$\displaystyle\frac{\sqrt{\gamma_{1}}}{2}\left[C^{c}_{a_{1}b_{1}}+(C^{c}_{a_{1}b_{1}})^{*}+C^{c}_{a_{1}^{\dagger}b_{1}}+(C^{c}_{a_{1}^{\dagger}b_{1}})^{*}\right]dW^{X}_{1}(t)$
$\displaystyle-i$
$\displaystyle\frac{\sqrt{\gamma_{1}}}{2}\left[C^{c}_{a_{1}b_{1}}-(C^{c}_{a_{1}b_{1}})^{*}+C^{c}_{a_{1}^{\dagger}b_{1}}-(C^{c}_{a_{1}^{\dagger}b_{1}})^{*}\right]dW^{P}_{1}(t)$
(34)

Setting $\text{\it{g}}_{c}=\Gamma_{c}$ once again ensures that the amplifier
does not see the QRC field, while the QRC is driven by the field from mode
$\hat{a}_{1}$. However, in addition to this coupling term, the form of
$\mathcal{L}_{c}$ introduces local damping $\propto\Gamma_{c}$ of the
$\hat{a}_{1}$ mode, which modifies the system’s dynamics. For convenience, we
choose damping parameters such that the total damping rate of system modes
$\hat{a}_{1}$ and $\hat{a}_{2}$ is equal after taking this additional damping
term into account; more precisely, we set
$\kappa_{1}+\Gamma_{c}=\kappa_{2}\equiv\kappa$. For simplicity, we further set
$\kappa_{1}=\Gamma_{c}=\kappa/2$, and $\omega_{d}=\omega_{a1}=\omega_{a2}$ so
that the detunings $\Delta_{a1}=\Delta_{a2}=0$.

We are now in a position to define our classification task: we wish to use a
QRC to classify $C=2$ different quantum states of the two-mode amplifier,
defined by specific values of the drive and amplifier interaction strengths
$(\eta^{(\sigma)},G_{1}^{(\sigma)},G_{12}^{(\sigma)})\in\\{(5.0,0.3,0.0),(8.0,0.0,0.3)\\}\kappa$.
Leaving details for Appendix I, these parameter choices generate amplifier
states with
$\langle\hat{X}_{a_{1}}^{(1)}\rangle=\langle\hat{X}_{a_{1}}^{(2)}\rangle\neq
0$, and
$\langle\hat{P}_{a_{1}}^{(1)}\rangle=\langle\hat{P}_{a_{1}}^{(2)}\rangle=0$,
namely with equal first-order moments for the mode $\hat{a}_{1}$ that the QRC
is directly coupled to. This is visible in the conditional dynamics of
$\langle\hat{X}_{a_{1}}\rangle^{c}$ plotted in Fig. 9(b), with blue and orange
colors corresponding to the two states. However the two states are distinct,
as is clear from their different second-order cumulants, whose unconditional
steady-state values are shown in the bar plots in Fig. 9(b).

Figure 9: (a) Measurement chain using a single-node QRC for the classification
of Gaussian states generated by a two-mode amplifier. A single mode
$\hat{a}_{1}$ of the system is coupled to the QRC using a directional hopping
interaction (a circulator). The single QRC node is continuously monitored via
heterodyne measurement. (b) Conditional dynamics of the $\hat{a}_{1}$ mode’s
$\langle\hat{X}_{a}\rangle^{c}$ quadrature, and steady-state unconditional
cumulants. (c) Same as (b), for the $\hat{b}_{1}$ QRC node’s
$\langle\hat{X}_{1}\rangle^{c}$ quadrature and unconditional cumulants.
Conditional dynamics are shown for one initialization of the measurement chain
per state. Colors blue and orange correspond to parameter sets
$(\eta^{(\sigma)},G_{1}^{(\sigma)},G_{12}^{(\sigma)})\in\\{(5.0,0.3,0.0),(8.0,0.0,0.3)\\}\kappa$
respectively. (d) Measured phase space quadratures
$\\{I_{1}^{X}(t_{f}),I_{1}^{P}(t_{f})\\}$ for (1) $N_{\rm S}=10$, and (2)
$N_{\rm S}=150$, for a linear QRC, $\bar{\Lambda}=0.0$. (e) Same as (d) but
for a nonlinear QRC (see text for parameters). (f) Classification accuracy as
a function of $N_{\rm S}$, computed for a test set of size $Q_{\rm
Test}=3000$; gray and black curves correspond to the linear and nonlinear QRCs
of (d) and (e) respectively.

####  V.1.2 QRC

Under directional coupling $\text{\it{g}}_{c}=\Gamma_{c}$, the amplifier
$\hat{a}_{1}$ mode drives the QRC via the $\hat{b}_{1}$ mode. For specificity,
we can look at the conditional evolution of the $\hat{X}_{1}$ quadrature of
the QRC:

$\displaystyle\frac{d}{dt}\langle\hat{X}_{1}\rangle^{c}$ $\displaystyle={\rm
tr}\\{\mathcal{L}_{\rm
QRC}\hat{\rho}^{c}\hat{X}_{1}\\}-\Gamma_{c}\langle\hat{X}_{a1}\rangle^{c}-\frac{1}{2}\Gamma_{c}\langle\hat{X}_{1}\rangle^{c}$
$\displaystyle+$
$\displaystyle\frac{\sqrt{\gamma_{1}}}{2}\left[C^{c}_{b_{1}b_{1}}+(C^{c}_{b_{1}b_{1}})^{*}+2C^{c}_{b_{1}^{\dagger}b_{1}}\right]dW^{X}_{1}(t)$
$\displaystyle-i$
$\displaystyle\frac{\sqrt{\gamma_{1}}}{2}\left[C^{c}_{b_{1}b_{1}}-(C^{c}_{b_{1}b_{1}})^{*}\right]dW^{P}_{1}(t)$
(35)

The QRC nodes thus evolve under the amplifier signal and also experiences
local damping due to the coupling superoperator, in addition to its internal
dynamics governed by $\mathcal{L}_{\rm QRC}$, and conditional evolution that
depends on cumulants of the QRC state. Typical conditional dynamics of
$\langle\hat{X}_{1}\rangle^{c}$ are shown in Fig. 5(c) for the two amplifier
states. The specific single-node QRC parameters we have chosen here are
$\gamma_{1}=\kappa$, $\Lambda_{1}=0.0027\gamma_{1}$, and
$\Delta_{1}=-\gamma_{1}$.

The QRC is driven to distinct quantum states under the two amplifier state
inputs, as is clear from the unconditional second-order QRC cumulants shown in
the bar plots in Fig. 5(c). What is not immediately clear is that the
resulting QRC states can also differ in their first-order moments, even though
the two amplifier states do not; this has implications for state
classification, as we analyze next.

###  V.2 Training and classification using finitely-sampled measurement
distributions

Measurement records $\\{J_{1}^{X,P}(t)\\}$ obtained from the measurement chain
contain information about QRC dynamics that will enable state classification.
In using these measurement records to construct the measured QRC output, we
will consider a slightly different scheme that is also experimentally
relevant, to highlight the flexibility afforded by our description.

In contrast to the pointer state classification task, the measured QRC output
$\mathbf{x}(t)$ here is constructed as the vector of measured quadratures
obtained using both temporal filtering and ensemble-averaging of heterodyne
measurement records, calculated over $N_{\rm S}$ distinct runs of the
measurement chain, or ‘shots’,

$\displaystyle\mathbf{x}(t)\equiv\begin{pmatrix}I_{1}^{X}(t)\\\
I_{1}^{P}(t)\end{pmatrix},\leavevmode\nobreak\
I_{1}^{X,P}(t)=\Big{<}\frac{1}{t}\int_{t_{0}}^{t}d\tau\leavevmode\nobreak\
J_{1}^{X,P}(\tau)\Big{>}_{N_{\rm S}}.$ (36)

This more general form of readout (which reduces to the single-shot case when
$N_{\rm S}=1$) can be used to reduce noise power in readout using multiple
measurement records obtained for the same inputs to the measurement chain.
Measured quadratures are themselves stochastic, so that the QRC output
provides a distinct stochastic trajectory for repeated runs of the measurement
chain. QRC output corresponding to different parameters
$\\{\eta^{(\sigma)},G_{1}^{(\sigma)},G_{12}^{(\sigma)}\\}$ is thus labelled
$\mathbf{x}^{(\sigma)}_{(q)}$, where $q$ is the trajectory label. Note that
each trajectory $q$ for a given $\sigma$ is constructed using $N_{\rm S}$
measurement records, each obtained from a distinct run of the measurement
chain for the same input parameters
$\\{\eta^{(\sigma)},G_{1}^{(\sigma)},G_{12}^{(\sigma)}\\}$.

Secondly, instead of integrating the entire obtained measurement records, we
introduce a waiting time $t_{0}$. By choosing $t_{0}$ to be much longer than
the amplifier response time for the two states (specifically,
$t_{0}=35.0/\kappa\gg\frac{2}{\kappa-2G_{1}^{(1)}}\simeq\frac{2\kappa}{\kappa^{2}-4(G_{12}^{(2)})^{2}}$),
we ensure that initial transients have settled, so that any information about
the states contained in their slightly-distinct transient dynamics is lost and
cannot influence classification.

The classification task is then set up as in Sec. IV. We employ an output
layer to predict a class label $\sigma^{p}$ using the measured QRC output,
$\sigma^{p}=f_{N}\\{\mathbf{W}_{\rm
O}\mathbf{x}^{(\sigma)}_{(q)}(t_{f})+\mathbf{b}\\}$. Here, $\mathbf{W}_{\rm
O}$ defines a matrix of time-independent trainable weights, $\mathbf{b}$ a
vector of trainable biases, while $f_{N}\\{\cdot\\}$ is the same fixed (that
is, untrained) normalizing function from earlier. Perfect classification
yields $\sigma^{p}=\sigma$ for all $\sigma=1,\ldots,C$. Note that here we
consider the measured phase space as being spanned by measured quadratures
evaluated at the end of the measurement, namely
$\mathbf{x}(t_{f})=\\{I_{1}^{X}(t_{f}),I_{1}^{P}(t_{f})\\}$ with
$t_{f}=120.0/\kappa$. This corresponds to classification using a fixed
measurement time; we will instead analyze performance as a function of the
other measurement resource, the number of shots $N_{\rm S}$.

Beginning with the case of a complete linear QRC ($\bar{\Lambda}=0$), the QRC
outputs $\\{\mathbf{x}^{(\sigma)}_{(q)}\\}$ are shown in measured phase space
in Fig. 9(d), for $q=1,\ldots,1500$ per state $\sigma$ (blue and orange). Each
panel shows outputs for the indicated value of $N_{\rm S}$. We see (1) that
for $N_{\rm S}=10$, the distributions corresponding to the two states have
distinct profiles but are overlapping, and hence not separable by a linear
decision boundary (black dashed line). The spread of the distributions can be
decreased by averaging over an increasing number of shots; this is seen in
(2), where $N_{\rm S}=150$ (the phase space area shown is the same as in (1),
allowing direct comparison of the distribution area). However, the centers of
the measured distributions are unmoved, so that they still overlap and their
linear separation is not possible.

The observed distributions are very different for an appropriately-chosen
nonlinear QRC. The measured phase space is shown in Fig. 9(e), again for
increasing values of $N_{\rm S}$. For (1) $N_{\rm S}=10$, we see that the
distributions are significantly distorted in comparison to the linear QRC,
showing amplified and squeezed axes due to the Kerr nonlinearity. However the
distributions still overlap significantly and are not linearly separable. With
increasing (2) $N_{\rm S}=150$, the width of the distributions decreases, and
a relative displacement of their centers (means) ensures that they can be
separated by the dashed black linear decision boundary, allowing for perfect
classification.

Note that both the linear and nonlinear QRC’s are driven to distinct quantum
states under input from the two states to be classified; this is clear from
the distinct distributions observed in the measured phase space. However, that
this difference manifests in a change to first order moments that are linearly
related to the measured quadratures is a signature of nonlinear quantum
systems, and not specific to the Kerr model we have chosen here: first-order
moments can couple directly to higher-order moments, and can thus be modified
by differences in these higher-order moments. Precisely this effect was
observed in the steady-state response of a single driven Kerr QRC node in
Figs. 3(b), (c).

As before, an arbitrary nonlinear QRC will not generally carry out this
classification task optimally. To understand the required operating regimes,
it proves useful to analyze the classification performance more
quantitatively. Using a training set $\\{\mathbf{x}^{(\sigma)}_{(q)}\\}$ for
$q=1,\ldots,Q_{\rm Train}=1500$ QRC outputs per state, we learn the optimal
weights ${\mathbf{W}}_{\rm O}$ and biases $\mathbf{b}$ defining the linear
decision boundary separating measured output distributions. Then, the trained
QRC is used to predict state labels $\sigma^{p}$ for a distinct test set
$\\{\mathbf{x}^{(\sigma)}_{(q)}\\}$ for $q=1,\ldots,Q_{\rm Test}=3000$, and
the classification accuracy as a function of $N_{\rm S}$ is calculated. The
result is plotted in Fig. 9(f), in gray and black for the linear and the
specific nonlinear QRC respectively, and is characterized by two metrics: the
maximum classification accuracy $\mathcal{C}_{\rm max}$ as before, and the
number of shots used in constructing measured quadratures to reach this
maximum, $N_{\rm S}^{\rm max}$. In the following section, we will analyze
these performance metrics as a function of QRC hyperparameters to determine
constraints that enable the amplifier state classification task.

Figure 10: (a) $\mathcal{C}_{\rm max}$ (left-hand axis) as a function of
nonlinearity $\bar{\Lambda}$, for fixed detuning
$\bar{\Delta}/\bar{\gamma}=-1.0$. Where $\mathcal{C}_{\rm max}\geq 0.99$
(yellow shaded region), $N_{\rm S}^{\rm max}$ is plotted along the right-hand
axis (green). Inset indicates the classical phase diagram. (b) Unconditional
QRC node amplitude difference
$\mathcal{B}^{12}=|\langle\hat{b}_{1}^{(1)}\rangle-\langle\hat{b}_{1}^{(2)}\rangle|$
calculated using TEOMs as a function of nonlinearity $\bar{\Lambda}$. Lower
panel shows unconditional QRC cumulants corresponding to the states
$\hat{\rho}_{\rm QRC}^{(1),(2)}$ (blue and orange respectively) in a region of
small $\mathcal{B}^{12}$ (left) and large $\mathcal{B}^{12}$ (right). Note the
difference in magnitude of cumulants. (c), (d) Same as (a), (b) respectively,
but now as a function of detuning $\bar{\Delta}$, for fixed nonlinearity
$\bar{\Lambda}/\bar{\gamma}=0.0027$, $\mathcal{N}_{\rm eff}\simeq 0.354$.

###  V.3 Classification performance as a function of hyperparameters

Hyperparameters.$-$ For a single-node QRC, the parameter space to vary is
reduced relative to the task considered in Sec. IV, which employed a two-node
QRC. Regardless, for consistency with our prior analysis we employ a similar
notation, setting $\Lambda_{1}\equiv\bar{\Lambda}$,
$\Delta_{1}\equiv\bar{\Delta}$, and setting the damping rate equal to the
(common) decay rate $\kappa$ of amplifier modes,
$\gamma_{1}\equiv\bar{\gamma}=\kappa$. Once again, all remaining aspects of
the measurement chain and training of the output layer remain unchanged as
hyperparameters are varied and QRC performance observed. Our key findings are
summarized at the end of this section.

####  V.3.1 Dependence on nonlinearity $\bar{\Lambda}$

We begin by fixing the QRC detuning $\bar{\Delta}/\bar{\gamma}=-1.0$, and
calculating $\mathcal{C}_{\rm max}$ for varying QRC nonlinearity
$\bar{\Lambda}$; the results are shown in Fig. 10(a). Clearly, an optimal
range of nonlinearity strengths exists where $\mathcal{C}_{\rm max}\to 1$. For
the orange shaded region where $\mathcal{C}_{\rm max}\geq 0.99$, we also plot
$N_{\rm S}^{\rm max}$ along the right-hand axis (green curve). A nonlinearity
strength exists for which $N_{\rm S}^{\rm max}$ is minimized, thereby reducing
the measurement resources required for classification.

The binary nature of the current task allows us to extract further intuition
about the observed QRC performance. As discussed earlier, the ability to
linearly separate integrated output records processed by the QRC corresponding
to the two amplifier states depends on the measured noisy distributions
acquiring a sufficient displacement. The unconditional QRC state variable that
is a measure of this displacement is simply the difference in QRC node
amplitude for the two states,
$\mathcal{B}^{12}\equiv|\langle\hat{b}_{1}^{(1)}\rangle-\langle\hat{b}_{1}^{(2)}\rangle|$.
We plot this quantity calculated using TEOMs in the top panel of Fig. 10(b)
for the same parameter regime as (a), together with the yellow shaded region
where $\mathcal{C}_{\rm max}\geq 0.99$. The lower panel shows unconditional
cumulants of the QRC state corresponding to the two amplifier states (blue and
orange respectively). Outside the yellow shaded region, the QRC state is
driven to distinct cumulants for the two amplifier states, but the magnitude
of cumulants, and the resulting displacement of QRC node amplitudes
$\mathcal{B}^{12}$, is small. However, within the yellow shaded region, the
distinct cumulants are much larger, leading to larger value of
$\mathcal{B}^{12}$. In contrast, for a linear QRC (dashed gray line),
$\mathcal{B}^{12}$ is always zero.

We emphasize here that a nonzero value of this displacement is by itself
insufficient for successful discrimination between the amplifier states. The
displacement is contained within noisy measurement records, and classification
amounts to being able to distinguish their filtered, ensemble-averaged
distributions in measured phase space. These factors are all accounted for by
our training on measured quadrature distributions, which yield all our results
for $\mathcal{C}_{\rm max}$ and $N_{\rm S}^{\rm max}$.

Operating regimes where $\mathcal{B}^{12}$ becomes nonzero can once again be
analyzed in the context of the classical phase diagram of the single-node QRC.
To this end, we introduce an effective drive $\eta_{\rm eff}$ experienced by
the QRC due to its coupling to the amplifier mode $\hat{a}_{1}$. Unlike the
QND interaction employed in Sec. IV which coupled only the cavity
$\hat{X}_{a}$ quadrature to the QRC, the directional hopping interaction here
ensures that mode $\hat{b}_{1}$ sees both quadratures of mode $\hat{a}_{1}$.
The effective drive $\eta_{\rm eff}^{(\sigma)}$ seen by the QRC for each
amplifier state then takes the form

$\displaystyle\eta_{\rm
eff}^{(\sigma)}=-\frac{\Gamma_{c}}{\sqrt{2}}\\!\left(\langle\hat{X}_{a_{1}}^{(\sigma)}\rangle+i\langle\hat{P}_{a_{1}}^{(\sigma)}\rangle\right)\\!=\frac{2\Gamma_{c}\kappa\cdot\eta^{(\sigma)}}{4(G^{(\sigma)}_{12})^{2}+\kappa(2G_{1}^{(\sigma)}\\!\\!-\\!\kappa)}$
(37)

where the second equality follows from the steady-state solutions for
$\langle\hat{X}_{a_{1}}^{(\sigma)}\rangle$,
$\langle\hat{P}_{a_{1}}^{(\sigma)}\rangle$ (see Appendix I). The current task
is specifically chosen so that $\eta_{\rm eff}^{(1)}=\eta_{\rm eff}^{(2)}$; we
hence set $\eta_{\rm eff}\equiv|\eta_{\rm eff}^{(1,2)}|$ to once again define
an effective nonlinearity $\mathcal{N}_{\rm eff}$ that defines the classical
phase diagram,

$\displaystyle\mathcal{N}_{\rm eff}=\frac{\eta_{\rm
eff}}{\bar{\gamma}+\Gamma_{c}}\sqrt{\frac{\bar{\Lambda}}{\bar{\gamma}+\Gamma_{c}}}$
(38)

similar to Eq. (30), but now with $\bar{\gamma}\to\bar{\gamma}+\Gamma_{c}$ to
account for the modified decay rate of the QRC due to the different form of
coupling $\mathcal{L}_{c}$. The resulting phase diagram is plotted in the
inset of Fig. 10(a). The dashed vertical line marks the cross-section along
which Fig. 10(a) is calculated, while the solid blue line coincides with the
yellow shaded region where $\mathcal{C}_{\rm max}\geq 0.99$. Clearly, this
region places the QRC in the vicinity of the classical bistability, where
cumulants are amplified, precisely as observed in Sec. III for a single-node
QRC under coherent driving.

####  V.3.2 Dependence on detuning $\bar{\Delta}$

Next we explore the classification performance as a function of detuning
$\bar{\Delta}$ for fixed nonlinearity strength
$\bar{\Lambda}/\bar{\gamma}=0.0027$, which sets $\mathcal{N}_{\rm eff}\simeq
0.354$ via Eq. (38). These results are shown in Fig. 10(c). Once more, we find
that an optimal range of detuning values exists for which $\mathcal{C}_{\rm
max}$ approaches unity and $N_{\rm S}^{\rm max}$ is minimized. Analysis of
$\mathcal{B}^{12}$ as a function of detuning values in Fig. 10(d) shows how
the QRC node amplitude displacement peaks around a detuning value
$\bar{\Delta}/\bar{\gamma}\simeq-1.0$ for the specific $\mathcal{N}_{\rm eff}$
value considered. This operating regime corresponds to $\mathcal{C}_{\rm
max}\geq 0.99$ (yellow shaded region), and once again places the QRC in the
vicinity of the classical bistability (solid blue line, Fig. 10(c) inset).

####  V.3.3 Dependence on damping $\bar{\gamma}$

The QRC node damping rate appears as the natural scaling parameter for QRC
hyperparameters, just as in the pointer state classification task. Hence if
$\bar{\gamma}$ is varied, the remaining QRC hyperparameters
$\\{\mathcal{N}_{\rm eff},\bar{\Delta}/\bar{\gamma}\\}$ can be suitably
adjusted to ensure that the QRC operates in a nonlinear processing regime,
enabling successful classification.

Beyond parameterizing such nonlinear processing regimes, $\bar{\gamma}$ also
influences the QRC node response time and the amplitude of QRC output fields
for fixed input signal strengths. Recall that our setup of the amplifier state
classification task - using measurement records only beyond a time $t_{0}$
once transients have settled - effectively analyzes the QRC’s processing of
steady-state signals from the two-mode measured quantum system. Hence
considerations of the QRC responding promptly to a dynamically-evolving signal
as in the pointer state classification task do not apply directly. Regardless,
the QRC damping rate $\bar{\gamma}$ can influence the time $t_{0}$ if it is
slower than the amplifier response rate for the states under consideration;
more precisely, in this case we would choose $t_{0}\gg 1/\bar{\gamma}$. This
in principle increases the total required measurement time. On the other hand,
larger $\bar{\gamma}$ for fixed strength of incident drives will suppress the
amplitude of QRC output, which is in general detrimental to classification.
These considerations imply that $\bar{\gamma}$ should ideally be chosen to be
on the order of the response rate of the amplifier (or the measured quantum
system in general).

####  V.3.4 Summary

The results of this section identify a range of QRC hyperparameters
$\\{\mathcal{N}_{\rm eff},\bar{\Delta}/\bar{\gamma}\\}$ where the QRC is able
to successfully classify amplifier states with equal means, finding a strong
connection to operating points in the vicinity of the classical bistability of
the single coherently-driven QRC node, similar to the pointer state
classification task. Importantly, the mechanism of classification is distinct:
the input-dependent enhancement of second-order cumulants of the QRC state,
whose nonlinearity leads to this difference manifesting via the QRC node
amplitudes.

###  V.4 Classical vs. quantum regimes of operation

As with the task of pointer state classification, we now analyze QRC
performance for amplifier state classification across decreasing input drives
to the measurement chain. To this end, we consider different instances of this
task, defined by the same two-mode amplifier interaction strengths
$G_{1}^{(\sigma)},G_{12}^{(\sigma)}$, but different values of drive input to
the measurement chain,
$(\eta^{(1)}/\bar{\gamma},\eta^{(2)}/\bar{\gamma})\in\\{(2.0,3.2),(5.0,8.0),(10.8,12.0)\\}$.
Via Eq. (38), these instances yield effective QRC drive strengths $\eta_{\rm
eff}/\bar{\gamma}\in\\{20.0,12.5,5.0\\}$, for which we calculate the
classification accuracy metrics $\mathcal{C}_{\rm max}$ and $N_{\rm S}^{\rm
max}$ as a function of nonlinearity $\bar{\Lambda}$. We have fixed
$\bar{\Delta}/\bar{\gamma}=0.0$ and $\bar{\gamma}=\kappa$ as before.

We plot $\mathcal{C}_{\rm max}$ in the lower panel of Fig. 11 for the
corresponding values of $\eta_{\rm eff}/\bar{\gamma}$ (solid blue, purple, and
peach respectively), with the yellow shaded region indicating nonlinearity
strengths for which $\mathcal{C}_{\rm max}\geq 0.99$. Similar to the pointer
state classification task, we see that for lower $\eta_{\rm eff}$, the
nonlinearity strength required for successful classification increases. This
is captured by the form of effective nonlinearity $\mathcal{N}_{\rm eff}$;
colored crosses indicate QRCs with the same value of $\mathcal{N}_{\rm
eff}\simeq 0.354$ for the different instances of this classification task,
which all succeed at classification but require stronger $\bar{\Lambda}$ for
weaker $\eta_{\rm eff}$.

Decreasing the input power therefore requires QRCs with stronger nonlinearity,
providing a natural probe of the quantum regime of QRC operation. To analyze
QRC performance across this transition, we also plot $N_{\rm S}^{\rm max}$ in
the lower panel of Fig. 11 (corresponding lighter colors, right-hand axis).
Interestingly, we find that for decreasing $\bar{\Lambda}$ indicating the
classical limit of QRC operation, the ampifier state classification task
becomes more difficult, as the number of shots $N_{\rm S}^{\rm max}$ required
to obtain $\mathcal{C}_{\rm max}\geq 0.99$ increases. Conversely, for
operation in more quantum regimes, $N_{\rm S}^{\rm max}$ decreases, so that
the task has a reduced resource cost.

Figure 11: $\mathcal{C}_{\rm max}$ (solid curves) and $N_{\rm S}^{\rm max}$
(lighter curves, right-hand axis) as a function of average nonlinearity
hyperparameter $\bar{\Lambda}$ for instances of the amplifier state
classification task with different effective drive amplitudes $\eta_{\rm eff}$
as indicated. Crosses indicate QRCs with fixed $\mathcal{N}_{\rm eff}\simeq
0.354$ for the three considered task instances; for weaker $\eta_{\rm
eff}/\bar{\gamma}$, corresponding QRCs have stronger nonlinearity and operate
in more quantum regimes (see text for more details). Top panel plots the
unconditional QRC node amplitude difference
$\mathcal{B}^{12}=|\langle\hat{b}_{1}^{(1)}\rangle-\langle\hat{b}_{1}^{(2)}\rangle|$,
which increases for QRCs operating in more quantum regimes. The dashed gray
line indicating the (vanishing) value of $\mathcal{B}^{12}$ for linear QRCs.

To understand this observation further, we use the TEOMs to calculate the
difference in the unconditional QRC node amplitude for the two states,
$\mathcal{B}^{12}=|\langle\hat{b}_{1}^{(1)}\rangle-\langle\hat{b}_{1}^{(2)}\rangle|$.
The result is plotted in the top panel of Fig. 11 for each considered instance
of the classification task. We clearly see that $\mathcal{B}^{12}$ increases
as we move towards QRCs operating deeper in the quantum regime (blue to purple
to peach). For conditional evolution, this increase manifests in a larger
relative displacement of measured quadrature distributions, resulting in a
reduction in ensemble-averaging required for successful classification.
Conversely, extrapolating the results to the classical limit
($\bar{\Lambda}\to 0$ for fixed $\mathcal{N}_{\rm eff}$), we see that
$\mathcal{B}^{12}\to 0$, so that the QRC is unable to distinguish the
amplifier states in the classical limit.

We thus observe a qualitative difference in QRC performance across classical
and quantum regimes when compared to the pointer state classification task.
This is because unlike the latter, the information that distinguishes
amplifier states is encoded in their second-order cumulants. In the classical
limit governed primarily by first-order cumulants of the measurement chain,
the influence of these second-order cumulants is suppressed; this directly
reduces $\mathcal{B}^{12}$, rendering the classification task more difficult
(requiring larger $N_{\rm S}^{\rm max}$). Conversely, in more quantum regimes
the influence of second-order cumulants on QRC node amplitudes is enhanced,
reducing $N_{\rm S}^{\rm max}$. These results indicate that the nature of
input encoding into QRC evolution can play an important role in extracting
possible processing advantages in quantum operating regimes.

However, one aspect of probing the quantum regime is shared across both
classification tasks: the absolute amplitude of measured QRC outputs is
reduced for weaker input drives in this regime. For the amplifier state
classification task, the relative difference in QRC node amplitudes
$\mathcal{B}^{12}$ is still larger in this regime, which is why the required
$N_{\rm S}^{\rm max}$ decreases in our simulations of the complete measurement
chain. In a real cQED experiment, however, the reduced amplitude can lead to
the desired signal being swamped by noisy (classical) HEMT amplifiers that are
also part of the measurement chain. As discussed earlier, this limitation is
overcome by including quantum-limited amplifiers as part of the measurement
chain, in this case to amplify the QRC output at cryogenic temperatures. We
discuss the inclusion of such an amplifier for the current task in Appendix L.
We find that even when using phase-preserving quantum-limited amplifiers that
add unavoidable amplified quantum noise to the measured quadrature
distributions, the associated amplification of the difference in QRC node
amplitudes $\mathcal{B}^{12}$ ensures that the benefit of operating the QRC in
more quantum regimes is preserved.

##  VI Discussion and Outlook

In this paper we present and analyze a framework for quantum reservoir
computing that we believe is timely within the current landscape of quantum
computation. Research efforts geared towards the realization of a universal
quantum computer to outperform classical computational paradigms have
accelerated in recent years [74, 75]; however, several important challenges
remain in scaling nascent quantum computers to larger and larger numbers of
quantum degrees of freedom, in particular maintaining coherence times,
correcting for the inevitable errors that occur during computation, and
managing growing hardware complexity [76]. Nevertheless, as part of these and
related efforts, modern experiments have continued to evolve in scope and
complexity, realizing extended multimode nonlinear quantum systems with an
unprecedented level of control. This naturally raises the question of whether
such near-term quantum devices can already be used for nontrivial computation.
This paper lays the theoretical foundations to assess the potential
computational capabilities of physical quantum systems as they are currently
realized, with constrained input and output control and subject to specific
measurement schemes. The machine-learning paradigm of reservoir computing
using physical systems, with relaxed requirements over internal device
optimization, is ideally suited to harness these capabilities.

To this end, our strategy is to develop a framework that (1) is built upon a
quantum-mechanical description of a physical quantum systems as part of its
complete quantum measurement chain, outputs from which are sampled via a
prescribed measurement scheme, and (2) has a physically-motivated classical
limit of operation. Such a framework allows us to determine performance
metrics, such as fidelity and time-to-solution, of performing a computational
task subject to operational constraints on inputs and outputs, and using only
accessible information obtained via measurement, therefore accounting for
experimental resource costs such as measurement time and repeated runs. Then,
comparing these metrics across well-defined classical and quantum regimes of
reservoir operation pinpoints advantages of quantum reservoir processing for a
given computational task.

We believe such an approach can be applied broadly in seeking advantages or
disadvantages to computation using quantum devices across very general
computational tasks. For a concrete analysis, in this paper we realize this
framework in the cQED architecture, and apply it to two quantum state
classification tasks: states with distinct means (task I: pointer state
classification), and states with equal means but distinct variances (task II:
amplifier state classification). This requires developing an efficient
computational approach to simulate the conditional dynamics of nonlinear
quantum measurement chains incorporating Josephson-junctions, signal-
generating quantum systems, with the entire chain evolving under continuous
heterodyne measurement. We present such an approach based on a set of
Stochastic Truncated Equations Of Motion (STEOMs) for measured observables and
quantum cumulants. Different from earlier quantum reservoir computing
approaches, this allows us to consider training of a linear output layer
directly on measured outputs, which in addition to measurement noise include
quantum fluctuations non-linearly transformed by the action of the quantum
reservoir itself.

The imposition of a trainable linear output layer applied to obtained
measurement records is an important choice: it allows any useful nonlinear
processing that enables quantum state classification to be attributed to
nonlinearity in the quantum measurement chain, and not to post-processing
operations. In our realization, this nonlinearity is provided by the QRC, and
is in fact critical: considering linear QRCs under the same operating
conditions and output processing does not yield perfect classification. We
note that the required nonlinearity could in principle be provided by other
elements of a measurement chain, for example via alternate measurement schemes
such as photodetection, or distinct input encoding schemes. This intuition is
supported by work on classical reservoirs, which have demonstrated the power
of computation with entirely linear reservoirs but employing non-linear
readout layers [4, 77]. This flexibility speaks to the ability to implement
physical reservoir computing across a range of quantum device platforms,
tailored to specific available nonlinearities and measurement schemes. A
comparison of the relative performance of different quantum reservoir
computing implementations is an interesting avenue for future exploration.

The ability to analyze QRC performance across classical-to-quantum transitions
within this framework provides a systematic approach to identifying the role
of reservoir processing in regimes where the operation is governed by quantum
statistical features. Carrying out this analysis, we find that pointer states,
distinguished by their first-order cumulants, can be successfully classified
even by QRCs operating in the classical limit. In addition, from the point of