# Quantum neuromorphic approach to efficient sensing of gravity-induced
entanglement

Tanjung Krisnanda  School of Physical and Mathematical Sciences, Nanyang
Technological University, 637371 Singapore, Singapore  Tomasz Paterek
Institute of Theoretical Physics and Astrophysics, Faculty of Mathematics,
Physics and Informatics, University of Gdańsk, 80-308 Gdańsk, Poland  Mauro
Paternostro  School of Mathematics and Physics, Queen’s University, Belfast
BT7 1NN, United Kingdom  Timothy C. H. Liew  School of Physical and
Mathematical Sciences, Nanyang Technological University, 637371 Singapore,
Singapore  MajuLab, International Joint Research Unit UMI 3654, CNRS,
Université Côte d’Azur, Sorbonne Université, National University of Singapore,
Nanyang Technological University, Singapore

###### Abstract

The detection of entanglement provides a definitive proof of quantumness. Its
ascertainment might be challenging for hot or macroscopic objects, where
entanglement is typically weak, but nevertheless present. Here we propose a
platform for measuring entanglement by connecting the objects of interest to
an uncontrolled quantum network, whose emission (readout) is trained to learn
and sense the entanglement of the former. First, we demonstrate the platform
and its features with generic quantum systems. As the network effectively
learns to recognise quantum states, it is possible to sense the amount of
entanglement after training with only non-entangled states. Furthermore, by
taking into account measurement errors, we demonstrate entanglement sensing
with precision that scales beyond the standard quantum limit and outperforms
measurements performed directly on the objects. Finally, we utilise our
platform for sensing gravity-induced entanglement between two masses and
predict an improvement of two orders of magnitude in the precision of
entanglement estimation compared to existing techniques.

##  I introduction

Following the success of the use of neural networks across different fields of
science [1, 2, 3, 4] for detecting patterns in data, proposals have been set
forth in the quantum regime [5, 6]. In this direction, a particular quantum
neural network architecture has emerged – termed _quantum reservoir
processing_ , in analogy to classical reservoir computing [7]. In such
architecture, the quantum network serving as a processor is composed of
randomly interacting quantum systems (the _nodes_), not requiring precise
control. The function of this kind of network is learned in a training
procedure that measures the system individually and fixes only a single output
layer, which makes it experimentally friendly. This architecture has been
proposed for executing classical tasks [5, 8, 9] (showing performance
advantage over classical networks) and quantum tasks such as state
characterisation [10, 11], quantum state preparation [12, 13], gate
compression [14], and quantum metrology [15] (see Refs. [16, 6] for reviews).
Remarkably, for characterisation and metrological tasks, it is not necessary
to perform correlation measurements, and it suffices to measure only local
observables such as average occupation numbers or intensities of the network
nodes. The platform is versatile and it holds the potential to directly
estimate important quantities such as quantum entanglement, which is the focus
of our study.

Entanglement is a special type of correlation between two or more objects, the
presence of which witnesses their quantum nature [17]. In experiments
involving objects that cannot be accessed directly, their quantum character
could be revealed by using such inaccessible systems as mediators between two
accessible probes. The revelation of an entanglement gain between the probes
then provides proof of another quantum signature – known as quantum discord –
of the mediators [18]. This experimental scheme has been put forward as a
proposal to probe quantum signatures of gravity through the observation of
gravity-induced entanglement between masses [19, 20, 21] (see also Refs. [22,
23, 24, 25, 26, 27, 28] for recent developments and discussion). This
motivates the general framework presented in this paper, which is aimed at
sensing (possibly weak) entanglement and its application to gravity-induced
entanglement.

To date, there are essentially two main schemes proposed for the observation
of gravity-induced entanglement, which suffer of different practical
difficulties. The Bose _et al._ -Marletto-Vedral (BMV) scenario [19, 20]
requires preparation of a macroscopic superposition of each of two nearby
massive bodies, whose later dynamics might showcase gravitational
entanglement. The challenging state-preparation stage is bypassed in the
proposal of Ref. [21], which resorts to continuous-variable (CV) entanglement
between masses that begin in natural and easy-to-arrange Gaussian states. In
this scheme, the entanglement detection remains as a demanding step. We show
that a relatively simple neural-network architecture is sufficient to achieve
a two-order-of-magnitude improvement in the precision with which entanglement
of massive systems can be estimated, when compared to state-of-the-art values
[29].

Specifically, we utilise a reservoir quantum network (QN) for precise
entanglement sensing. In particular, quantum objects whose entanglement we
want to scrutinise (the input) are put in contact with a QN. The observables
from the QN – which can be as simple as the mean excitation numbers of the
nodes – are post-processed through a single output layer. This layer is
trained so that the final output estimates quantum entanglement of the input
objects. Our general platform is particularly useful in situations where the
input is not accessible for direct measurements, the latter are complicated
(this is particularly the case for those that necessitate conditional or
correlation measurements), or in cases where the input is less resilient to
measurement errors than the QN. First, we will introduce the general framework
with generic quantum systems. We show that a QN can learn from a random set of
non-entangled input states and nevertheless is able to estimate the amount of
entanglement at the testing stage. For a more realistic scenario considering
measurement errors, we show that the entanglement precision scales better than
$\Delta_{E}\propto 1/\sqrt{\mathcal{N}}$, where $\mathcal{N}$ is the number of
measured observables. We shall refer to $1/\sqrt{\mathcal{N}}$ scaling as the
standard quantum limit (SQL). Finally, we demonstrate an explicit application
of our framework and its features to the recent endeavour whose goal is to
reveal quantum features of gravity by measuring gravity-induced entanglement
between masses. In particular, we show that measurements on cavity modes,
which have interacted with the masses, can be post-processed to estimate the
gravity-induced entanglement. Importantly, our approach offers better
sensitivity compared to direct measurements on the masses.

##  II The general framework

Our thought platform is depicted in Fig. 1. Consider that quantum objects,
whose entanglement is to be estimated, serve as the input. They come in
contact with a processor, namely, a QN composed of quantum nodes. Note that
the QN nodes require minimal control, e.g., they can be randomly interacting
with each other, the input, and environment. We also allow that they are
pumped by external coherent sources (e.g., in an optical system, lasers). The
purpose of the contact is a flow of information from the input to the QN. By
retrieving the observables from the QN and processing them via a single output
layer, one obtains a final output signal. We will show that by training a set
of weights and biases in the single output layer, the output signal estimates
entanglement of the input objects.

Figure 1: Illustration of a quantum neuromorphic platform for entanglement
sensing. It involves to-be-measured input objects that are connected to a
quantum network composed of uncontrolled nodes having random interactions with
the input $K_{nm}$ and between themselves $J_{mm^{\prime}}$. Both the input
and the nodes interact with their environments, denoted by $\gamma_{n}$ and
$\kappa_{m}$, respectively. The observables from the QN are processed by a
trained output layer, producing a signal that estimates input entanglement.

Let us consider generic quantum systems and their dynamics with which we
demonstrate the general framework described above. In what follows, we
consider continuous-variable systems (bosons). Additionally, the platform also
works for generic discrete systems (e.g., qubits), see Appendix A for details,
as well as hybrid discrete-continuous systems, see Appendix B. We begin by
modelling the dynamics of the input $\rho_{\text{in}}$ and the QN
$\rho_{\text{qn}}$ for a time $\tau$, after which the observables of the QN
are recorded. The coherent part of the dynamics is described by the following
Hamiltonian, written in a frame rotating with the pump frequency

$\displaystyle H$
$\displaystyle=\sum_{n=1}^{2}\hbar\Delta_{n}\hat{a}_{n}^{\dagger}\hat{a}_{n}+\sum_{m=1}^{M}\hbar(\Lambda_{m}\hat{b}_{m}^{\dagger}\hat{b}_{m}+P_{m}(\hat{b}_{m}+\hat{b}_{m}^{\dagger}))$
(1) $\displaystyle+\sum_{\llbracket n,m\rrbracket}\hbar
K_{nm}\mathcal{F}(\hat{a}_{n},\hat{b}_{m})+\sum_{\llbracket
m,m^{\prime}\rrbracket}\hbar
J_{mm^{\prime}}\mathcal{F}(\hat{b}_{m},\hat{b}_{m^{\prime}}),$

where $\hat{a}_{n}$ ($\hat{b}_{m}$) denotes the annihilation operator for the
$n$th input object ($m$th QN node). The detunings of all local frequencies
$\\{\omega_{n},\Omega_{m}\\}$ with respect to the frequency of the pump
$\Theta_{p}$ are denoted by $\Delta_{n}=\omega_{n}-\Theta_{p}$ and
$\Lambda_{m}=\Omega_{m}-\Theta_{p}$. The contact between the input and QN is
represented by the couplings $K_{nm}$, whereas the interactions within the QN
are denoted by $J_{mm^{\prime}}$. For simplicity, we take the operator
function to represent interactions that are ample in nature, i.e.,
$\mathcal{F}(\hat{X},\hat{Y})\equiv\hat{X}\hat{Y}^{\dagger}+\hat{Y}\hat{X}^{\dagger}$.
Each QN node may be coherently driven with strength $P_{m}$. The bracket
$\llbracket\cdot,\cdot\rrbracket$ denotes a particular configuration of the
couplings, e.g., all-to-all.

We note that the simulation of the system can be made efficient when dealing
with Gaussian states [30]. These tools are applicable as the generic dynamics
we consider here preserves Gaussianity, i.e., it involves a Hamiltonian that
is at most quadratic in operators (Eq. (1)) and Gaussian dissipative processes
(see below). In this case, complete description of the system is contained in
a covariance matrix (CM) $\bm{V}$ with elements $V_{ij}\equiv\langle
u_{i}u_{j}+u_{j}u_{i}\rangle/2-\langle u_{i}\rangle\langle u_{j}\rangle$,
where the vector
${\bm{u}}\equiv[\hat{q}_{1},\hat{r}_{1},\hat{q}_{2},\hat{r}_{2},\hat{x}_{1},\hat{p}_{1},\cdots,\hat{x}_{M},\hat{p}_{M}]^{T}$
is composed of dimensionless position and momentum quadratures (of the input
and QN nodes, respectively) that are expressed as
$\hat{q}_{n}=(\hat{a}_{n}+\hat{a}_{n}^{\dagger})/\sqrt{2}$,
$\hat{r}_{n}=(\hat{a}_{n}-\hat{a}_{n}^{\dagger})/(i\sqrt{2})$,
$\hat{x}_{m}=(\hat{b}_{m}+\hat{b}_{m}^{\dagger})/\sqrt{2}$, and
$\hat{p}_{m}=(\hat{b}_{m}-\hat{b}_{m}^{\dagger})/(i\sqrt{2})$. One can obtain
the dynamics of the quadratures in the Heisenberg picture from the Hamiltonian
of Eq. (1), which with added noise terms gives rise to a set of Langevin
equations (LEs) that can be written in a matrix form:
$\dot{\bm{u}}(t)={\bm{A}}{\bm{u}}(t)+{\bm{h}}(t)$. The drift matrix ${\bm{A}}$
contains the parameters
$\\{\Delta_{n},\Lambda_{m},K_{nm},J_{mm^{\prime}},\gamma_{n},\kappa_{m}\\}$
and the vector ${\bm{h}}(t)$ incorporates the pump and noise terms, see
Appendix C for details. The noise terms are of uncoloured Gaussian type, and
written as $\sqrt{2\gamma_{n}}\hat{a}_{n}^{\text{in}}$ and
$\sqrt{2\kappa_{m}}\hat{b}_{m}^{\text{in}}$, where
$\langle\hat{a}_{n}^{\text{in}}(t)\hat{a}_{n^{\prime}}^{\text{in},\dagger}(t^{\prime})\rangle=\delta_{nn^{\prime}}\delta(t-t^{\prime})$
and
$\langle\hat{b}_{m}^{\text{in}}(t)\hat{b}_{m^{\prime}}^{\text{in},\dagger}(t^{\prime})\rangle=\delta_{mm^{\prime}}\delta(t-t^{\prime})$
[31].

The solution of the LEs is given by

${\bm{u}}(t)={\bm{W}}_{+}(t){\bm{u}}(0)+{\bm{W}}_{+}(t)\int_{0}^{t}dt^{\prime}{\bm{W}}_{-}(t^{\prime}){\bm{h}}(t^{\prime}),$
(2)

where ${\bm{W}}_{\pm}(t)=\exp{(\pm{\bm{A}}t)}$. This further gives the
dynamical equation for the CM:
$\dot{\bm{V}}={\bm{A}}{\bm{V}}(t)+{\bm{V}}(t){\bm{A}}^{T}+{\bm{D}}$, where
${\bm{D}}=\mbox{diag}[\gamma_{1},\gamma_{1},\gamma_{2},\gamma_{2},\kappa_{1},\kappa_{1},\cdots,\kappa_{M},\kappa_{M}]$.
The observables of the QN
$\langle\hat{O}_{mk}\rangle=\mbox{tr}(\rho(\tau)\hat{O}_{mk})$ ($k$ labels
observables from the same $m$th node) at time $\tau$ can be obtained from
${\bm{u}}(\tau)$ and ${\bm{V}}(\tau)$. For more detailed expressions, see
Appendix C. Here we consider local observables, for simplicity. We will see
that it is sufficient to work with average occupation numbers (intensities) as
the observables, although any additional variables that can be measured can
further help. The observables define an output layer upon which a training
procedure is used to find a linear combination of the observables that will
define the final system output.

The training is performed with ridge regression using a random set of input
CMs $\\{\bm{V}_{\text{in},l}\\}_{l=1}^{N_{\text{tr}}}$ as follows. Each of the
input CMs will be in contact with the QN and produce a set of $N_{\text{ob}}$
observables at time $\tau$, recorded as a vector $\bm{v}_{l}$. The observables
are used to first estimate the input state (its _unique_ elements), from which
entanglement is calculated. In the present case, each element of the CM
(labelled $s$) is estimated linearly as $f_{s}=\bm{\beta}_{s}[1;\bm{v}_{l}]$,
where $\bm{\beta}_{s}=[\beta_{0},\beta_{1},\cdots,\beta_{N_{\text{ob}}}]$
contains the coefficients to be obtained with ridge regression. In particular,
$\bm{\beta}_{s}=(\bm{X}^{T}\bm{X}+\lambda\mathbb{1})^{-1}\bm{X}^{T}\bm{Y}_{s}$,
where
$\bm{X}=[1,\bm{v}_{1}^{T};1,\bm{v}_{2}^{T};\cdots;1,\bm{v}_{N_{\text{tr}}}^{T}]$
contains all the observables in the training set, $\bm{Y}_{s}$ contains the
target $s$th element, and $\lambda$ is the ridge parameter. This allows us to
obtain an estimated input CM $\tilde{\bm{V}}_{\text{in}}$ from the trained
output layer $\\{\bm{\beta}_{s}\\}$, given measured QN observables.
Consequently, the estimated entanglement is computed using the logarithmic
negativity $E=L_{e}(\tilde{\bm{V}}_{\text{in}})$ [32]. In what follows, we
define the entanglement estimation error as

$\Delta_{E}=\sqrt{\sum_{l^{\prime}=1}^{N_{\text{te}}}\frac{(E_{\text{est},l^{\prime}}-E_{\text{in},l^{\prime}})^{2}}{N_{\text{te}}}},$
(3)

where $N_{\text{te}}$ is the number of random input CMs in the testing set.

##  III Entanglement estimation

Here we present the performance of entanglement estimation. In simulations,
the parameters are taken as random
$\\{\Delta_{n},\Lambda_{m},K_{nm},J_{mm^{\prime}},P_{m},10\gamma_{n},10\kappa_{m}\\}\in[0,1]\Gamma$,
where $\Gamma$ is an overall strength in units of frequency, and evolution
time $\tau=\pi/2\Gamma$. One set of random parameters will be taken to define
one particular QN. When assessing the performance of the scheme, we will
average over different parameter choices, to provide a general assessment of
the architecture rather than any specific parameter choice. Indeed, one
advantage of our scheme is that the considered systems do not need precise
control of their parameters.

Figure 2(a) shows the entanglement estimation error against the number of QN
nodes. The sudden shift shown by the arrow indicates $\Delta_{E}\sim 10^{-10}$
is obtained for QNs having at least $4$ nodes. This can be understood as
follows. Recall that the number of independent parameters required to fully
characterise an $N$-mode Gaussian state is $2N(2N+1)/2$. This suggests that to
faithfully estimate the state of a two-mode Gaussian input, one requires at
least $10$ observables from the QN. This is fulfilled by having at least $4$
QN nodes as each node is itself in a Gaussian state and hence requires three
independent real parameters (e.g., we take two diagonal and one off-diagonal
entries from the local CM) to be determined. The inset shows the entanglement
profile of the input CM ${\bm{V}}_{\text{in}}$ used in both training and
testing, with $N_{\text{tr}}=50$ and $N_{\text{te}}=100$, respectively. See
Appendix D for the generation of random input CMs. A closer look at the
comparison between the estimated and input entanglement during testing is
plotted in Figs. 2(b) and (c) for the case where the QN is composed of $3$ and
$4$ nodes, respectively. It can be seen that the latter offers minute errors.

Figure 2: Performance of entanglement sensing for generic dynamics. (a)
Estimation error vs number of nodes used in the QN. The inset shows the
profile of entanglement of the randomly generated input CMs in both training
and testing. Panels (b) and (c) present explicit comparison between the
estimated and input entanglement during testing where the QN is composed of
$3$ and $4$ nodes, respectively.

We have also simulated the case where we record one observable (the mean
excitation $\langle\hat{b}_{m}^{\dagger}\hat{b}_{m}\rangle$) from each QN
node. In this case, one requires either the addition of two-photon pump, i.e.,
$\sum_{m}P_{m}^{\prime}(\hat{b}^{2}_{m}+\hat{b}_{m}^{\dagger 2})$ with random
strengths (relatively weaker) $P_{m}^{\prime}\in[0,1]\Gamma/10$ or the
presence of ultra-strong coupling
$\mathcal{F}(\hat{b}_{m},\hat{b}_{m^{\prime}})=(\hat{b}_{m}+\hat{b}_{m}^{\dagger})(\hat{b}_{m^{\prime}}+\hat{b}_{m^{\prime}}^{\dagger})$.
The reason for this is that simpler interactions or drives in Eq. (1) are not
sufficient for complex information transfer during the dynamics, which would
allow the mean excitation of the QN nodes to completely recover information
regarding the input objects. We found that the shift to low estimation error
requires at least 10 QN nodes, again consistent with the number of independent
parameters of the input CM, see Appendix E for details.

As the scheme estimates the CM of the input objects before computing
entanglement, it opens up the possibility to use a training set consisting of
separable input CMs without affecting its entanglement-testing capabilities.
We used the same setup as in Fig. 2(a) based on 4 QN nodes and performed
training using only non-entangled input CMs. The testing was performed with
entangled input CMs, finding a profile similar to the inset in Fig. 2(a).
Indeed, the comparison between the estimated and input entanglement is similar
to the one in Fig. 2(c) (see Appendix F for details). We note that although
each input CM in the training set is not entangled, they are still correlated.
Similarly, learning from separable input objects is also possible for discrete
systems (cf. Appendix F).

##  IV Scaling beyond the SQL

For a more realistic model, we incorporate measurement errors of the
observables from the QN. The observables now read
$\langle\hat{O}_{mk}\rangle\rightarrow\langle\hat{O}_{mk}\rangle+\varepsilon_{mk}$,
where $\\{\varepsilon_{mk}\\}$ are generated from a normal distribution with
zero mean and standard deviation $\zeta/2$. In what follows, we take
$\zeta=10^{-3}$. Similar to Fig. 2(a), we present the estimation errors in
Fig. 3(a). The dots indicate the scaling of the error with respect to the
number of nodes $M$. From Fig. 3(a), one can see the signature of the shift
previously observed in Fig. 2(a). In particular, the scaling of the estimation
error becomes clearer for $M\geq 4$ in Fig. 3(a). It can be seen that
$\Delta_{E}$ can exhibit scaling beyond the dashed curve, i.e., beyond the SQL
$\propto 1/\sqrt{\mathcal{N}}=1/\sqrt{3M}\propto 1/\sqrt{M}$.

Another alternative to obtain independent observables from the QN is through
time-multiplexing. For instance, we consider a single observable from each QN
node, i.e., the mean excitation
$\langle\hat{b}_{m}^{\dagger}\hat{b}_{m}\rangle$ and measure it at
$\mathcal{T}$ different times. This gives a total of
$N_{\text{ob}}=M\mathcal{T}$ observables. We demonstrate the case for
$\mathcal{T}=3$, i.e., at $\tau=\\{1,2,3\\}\pi/2\Gamma$ in Fig. 3(b). In this
case, we have added random two-photon pump $P_{m}^{\prime}\in[0,1]\Gamma/10$
(see Appendix E for the case with ultra-strong coupling). One can see similar
scaling as in panel (a).

Figure 3: Scaling beyond the SQL. (a) Estimation error plotted against the
number of nodes $M$. Panel (b), taking only mean excitation from each QN node
as the observable with time-multiplexing performed $\mathcal{T}=3$ times. The
scaling for the SQL is given by the dashed curve in each panel.

##  V Gravity-induced entanglement

We present an application of the entanglement sensing scheme to estimate
gravity-induced entanglement (GIE) generated between masses. Consider two
identical spherical objects, each with mass $m$, trapped in a 1D harmonic
potential. This configuration has been theoretically predicted to generate
entanglement between the masses through gravitational interactions [21]. Here,
each mass is probed by a cavity mode, see Fig. 4(a). The probes are turned on
by the pump on the respective cavities, $\mathcal{E}_{a}$ (left) and
$\mathcal{E}_{b}$ (right). This way, the observables from the cavity modes can
be processed through an output layer, which then produces an estimate of the
GIE.

Figure 4: Sensing gravity-induced entanglement (GIE) between masses. (a) The
setup of two trapped masses interacting gravitationally. After a time
$\tau_{0}$, two cavity modes are turned on, whose observables are processed to
estimate the GIE. (b) The estimated GIE for different initial evolution of the
masses $\tau_{0}$ vs testing instances. One observes that the standard
deviation $\delta_{E}<10^{-4}$. The parameters used in simulations are $m=1$
kg made of Osmium with density $22.59$ g/cm3, $\omega=0.1$ Hz with
$\omega/\gamma\gg 1$, $r_{0}=1.73$, $\tau=1$ $\mu$s, $L_{a(b)}=25$ mm, laser
wavelength $1064$ nm, $P_{a(b)}=50$ mW,
$\\{\kappa_{a},\kappa_{b},\Delta_{a},\Delta_{b}\\}=2.36\times 10^{5}$ Hz,
$L\simeq 2R$ with $R$ the radius of the mass, and $\zeta/2=10^{-2}$. See
Appendix G for references of the parameters.

First, we consider the dynamics without the probes, in which the Hamiltonian
reads

$H_{0}=\frac{\hbar\omega}{2}(\hat{p}_{A}^{2}+\hat{x}_{A}^{2}+\hat{p}_{B}^{2}+\hat{x}_{B}^{2})-\frac{\hbar
Gm}{\omega L^{3}}(\hat{x}_{A}-\hat{x}_{B})^{2},$ (4)

where $\hat{x}_{A(B)}$ denotes the dimensionless displacement of mass $A(B)$,
$\omega$ the frequency of the trapping potentials, and $L$ the equilibrium
distance between the masses. We have used
$\hat{x}_{A(B)}=x_{A(B)}\sqrt{m\omega/\hbar}$ and
$\hat{p}_{A(B)}=p_{A(B)}/\sqrt{\hbar m\omega}$, where $x_{A(B)}$ and
$p_{A(B)}$ are the displacement and momentum operators, respectively. The
gravitational interaction is expanded from
$-Gm^{2}/(L-(\hat{x}_{A}-\hat{x}_{B})\sqrt{\hbar/m\omega})$ up to a quadratic
term, $(\hat{x}_{A}-\hat{x}_{B})^{2}$, which is necessary for entanglement
generation as it contains non-local coupling $\propto\hat{x}_{A}\hat{x}_{B}$
acting on both masses. We have neglected the constant and linear term
$\propto(\hat{x}_{A}-\hat{x}_{B})$ as the former is simply an energy offset
and the latter a bi-local operator (cannot create entanglement) that
constitutes to shifting the equilibrium position of the masses. One can
construct a set of LEs from Eq. (4) with the addition of damping $\gamma$ and
Brownian-like noises $\hat{\xi}_{A(B)}$ affecting the masses (see Appendix G
for details). As we deal with Gaussianity-preserving dynamics, we use the
tools for CV systems. This includes the description of the system within a CM
and its evolution to $\bm{V}(\tau_{0})$ from which properties of the system
can be calculated (see Appendix G).

At time $\tau_{0}$, the probes are turned on, where the Hamiltonian (in a
rotating frame with the frequency of the lasers) now reads

$\displaystyle H$
$\displaystyle=H_{0}+\hbar\Delta_{0a}\hat{a}^{\dagger}\hat{a}+\hbar\Delta_{0b}\hat{b}^{\dagger}\hat{b}+i\hbar\mathcal{E}_{a}(\hat{a}^{\dagger}-\hat{a})$
(5) $\displaystyle+i\hbar\mathcal{E}_{b}(\hat{b}^{\dagger}-\hat{b})-\hbar
G_{0a}\hat{a}^{\dagger}\hat{a}\hat{x}_{A}+\hbar
G_{0b}\hat{b}^{\dagger}\hat{b}\hat{x}_{B},$

where $\hat{j}=\hat{a},\hat{b}$ denotes the annihilation operator of the left
and right cavity mode, $\Delta_{0j}=\omega_{j}-\omega_{lj}$ the cavity-laser
detuning, $\mathcal{E}_{j}=\sqrt{2P_{j}\kappa_{j}/\hbar\omega_{lj}}$ the
driving strength of the cavity, $P_{j}$ the laser power with frequency
$\omega_{lj}$, $\kappa_{j}=\pi c/2F_{j}L_{j}$ the cavity decay rate with
finesse $F_{j}$ and length $L_{j}$,
$G_{0j}=(\omega_{j}/L_{j})\sqrt{\hbar/m\omega}$ the optomechanical coupling
strength. From Eq. (5), one can construct a set of linearised LEs, which are
then used to evolve the CM $\bm{V}(\tau_{0})$ to $\bm{V}(\tau_{0}+\tau)$ at
which the observables from the cavity modes are recorded.

In what follows, we take into account the features shown previously for
entanglement sensing using generic systems. As the task is estimating
entanglement of a two-mode CM (of the masses), at least 10 observables are
required for recording. This is taken from 10 independent CM elements of the
joint cavity modes. From the central limit theorem it follows that
$\zeta/2\propto 1/\sqrt{N_{\text{rep}}}$, where $N_{\text{rep}}$ is the number
of repetitions that an element is measured. To make a comparison with
entanglement measurement in Ref. [29] whereby $N_{\text{rep}}=10^{4}$, we
shall assume error statistics with $\zeta\sim 2/\sqrt{N_{\text{rep}}}=2\times
10^{-2}$. As the initial CM at $t=0$, we use squeezed (local) thermal state
for the masses
$\mbox{diag}[e^{2r_{0}},e^{-2r_{0}},e^{2r_{0}},e^{-2r_{0}}](1+2\bar{n})/2$
with $r_{0}$ being the squeezing strength and $\bar{n}$ the mean thermal
phonon number, and vacuum for the cavity modes. The training is performed
using random separable input states $\bm{V}_{\text{in}}(\tau_{0})$, which are
generated using random $\bar{n}>0$. This is such that entanglement does not
yet grow for initial thermal states within $\tau_{0}$. On the other hand,
testing is performed with $\bar{n}=0$. For better precision, one can use time-
multiplexing during the dynamics with the probes at
$\\{1,2,\cdots,\mathcal{T}\\}\tau$.

We present the estimated GIE for different initial accumulation time
$\tau_{0}$ in Fig. 4(b). We have taken $\mathcal{T}=4$ (see Appendix H for the
scaling of standard deviation $\delta_{E}$ against $\mathcal{T}$),
$N_{\text{tr}}=50$, and $N_{\text{te}}=100$. The standard deviations of the
GIE in Fig. 4(b) follow $\delta_{E}<10^{-4}$, which is two orders of magnitude
better than the experimentally achieved $\sim 10^{-2}$ in Ref. [29]. We also
computed the estimated GIE from direct measurements, which is done by adding
measurement errors directly to the elements of $\bm{V}(\tau_{0})$. In this
case, standard deviation $\delta_{E}\sim 10^{-4}$ is only possible if the
system permits three orders of magnitude weaker measurement error strength
$\zeta/2=10^{-5}$. This demonstrates the efficiency of our method, which
requires less number of single-shot measurements $N_{\text{rep}}$ to obtain
precision comparable to measurements directly on the masses, i.e., with noisy
$\bm{V}(\tau_{0})$.

Figure 4(b) shows that our method is able to estimate GIE efficiently for
$\tau_{0}=0.5$ s. We note that this is shorter than the coherence times
resulting from thermal photons from environment and collisions with air
molecules (both in the range of about $5$ s) if the experiments were conducted
on Earth with liquid Helium in ultrahigh vacuum [21].

##  VI DISCUSSION

We have shown that a simple neural network (quantum reservoir processor) can
be used for efficient estimation of quantum entanglement. Our main motivation
for development of such a method is provided by present efforts to design
experiments capable of detection of gravity-induced entanglement. The
introduced method shows that the entanglement precision can be improved by two
orders of magnitude from what was achieved in Ref. [29].

The entanglement sensing step is crucial for masses initialised in natural
Gaussian states and any improvement on it relaxes other requirements of the
setup. The most direct one is the requirement on coherence times: since
smaller values of entanglement become detectable, the system can be measured
earlier. With entanglement estimation accuracy on the order $10^{-4}$
detection of GIE could be performed within decoherence times available on
Earth, whereas accuracy $10^{-2}$ would rather require an experiment in space
[21]. Moreover, in order to understand how other experimental parameters can
be changed, let us recall that the _figure of merit_ for entanglement
generated via gravity between trapped masses $m$ separated by a distance $L$
is given by $2Gm/\omega^{2}L^{3}$, where $\omega$ characterises the trapping
potential or spread of the initial wave function of each mass [21]. Therefore,
better entanglement precision also translates to smaller masses in the
experiment that could be placed further apart.

The method presented in this paper also holds potential for other settings
where one estimates entanglement of the easily accessed probes with precision
advantage and reveals quantumness of a macroscopic mediating object. In
particular, this includes an extension of Refs. [33, 34, 35, 36] towards
showing quantum properties of photosynthetic bacteria [37] or that of a
macroscopic mechanical membrane in the membrane-in-the-middle optomechanics
setting [38, 39, 18]. Additionally, we note that our scheme can work not only
for CV or discrete systems, but also hybrid configurations such as discrete
systems as input and CV systems as the QN or vice versa (see Appendix B).

###### Acknowledgements.

We thank Sanjib Ghosh, Kevin Dini, and Yvonne Gao for stimulating discussion.
T.K. and T.C.H.L. acknowledge the support by the Singapore Ministry of
Education under its AcRF Tier 2 grant MOE2019-T2-1-004. T.P. is supported by
the Polish National Agency for Academic Exchange NAWA Project No.
PPN/PPO/2018/1/00007/U/00001. MP acknowledges the support by the European
Union’s Horizon 2020 FET-Open project TEQ (766900), the Leverhulme Trust
Research Project Grant UltraQuTe (grant RGP-2018-266), the Royal Society
Wolfson Fellowship (RSWF/R3/183013), the UK EPSRC (EP/T028424/1), and the
Department for the Economy Northern Ireland under the US-Ireland R&D
Partnership Programme (USI 175).

Author contributions: TK, TP, and TCHL conceived the initial project
direction; TK carried out all the calculations and derivations under the
supervision of TP, MP, and TCHL; MP and TP assisted in designing the setup for
gravity-induced entanglement; TK wrote the paper with contributions from TP,
MP, and TCHL. All authors discussed the results and revised the paper.

Competing interests: The authors declare no competing interests.

Data availability: All data needed to evaluate the conclusions in the paper
are present in the paper and/or the Appendix.

##  Appendix A Generic discrete systems: Entanglement estimation and scaling

Here, we consider that all quantum systems, i.e., the input objects and the QN
nodes, are qubits. Each qubit has two energy levels, the ground state
$|g\rangle$ and excited state $|e\rangle$. Let us take the generic Hamiltonian
in Eq. (1) in the main text, where now $\hat{a}_{n}$ ($\hat{b}_{m}$) denotes
the lowering operator $|g\rangle\langle e|$ for the $n$th input qubit ($m$th
QN node).

In addition, the input and QN nodes may interact with their environment,
adding an incoherent element to the dynamics. We consider a simple dissipative
process such that the dynamics of the whole system is described within the
Lindblad master equation

$\dot{\rho}=-\frac{i}{\hbar}[H,\rho]+\sum_{n}\frac{\gamma_{n}}{2}\mathcal{L}(\rho,\hat{a}_{n})+\sum_{m}\frac{\kappa_{m}}{2}\mathcal{L}(\rho,\hat{b}_{m}),$
(6)

where $\mathcal{L}(\rho,\hat{X})\equiv
2\hat{X}\rho\hat{X}^{\dagger}-\\{\hat{X}^{\dagger}\hat{X},\rho\\}$ and the QN
nodes are initialised in their ground state $|g\rangle$. The dissipation rate
of the $n$th input and $m$th QN node are denoted by $\gamma_{n}$ and
$\kappa_{m}$, respectively. These processes are not essential for our scheme,
but are included to show robustness in their presence. After a time $\tau$,
the observables $\langle\hat{O}_{mk}\rangle=\mbox{tr}(\rho(\tau)\hat{O}_{mk})$
are recorded as a vector $\bm{v}$ and sent to a trained output layer (the
training with ridge regression is described in the main text). Note that the
index $k$ denotes different observables from the same $m$th QN node. The
trained output layer is used to estimate the unique elements of the input
state, giving us $\tilde{\rho}_{\text{in}}$, from which the entanglement is
quantified using negativity $E=N_{e}(\tilde{\rho}_{\text{in}})$ [32].

In simulations, the system parameters are randomised in the same way as that
described in Section III in the main text. The procedure to generate random
input states for training $\\{\rho_{\text{in},l}\\}_{l=1}^{N_{\text{tr}}}$ and
testing $\\{\rho_{\text{in},l^{\prime}}\\}_{l^{\prime}=1}^{N_{\text{te}}}$ are
described below in Section D.

We tested the scheme to estimate entanglement of two-qubit input states. The
estimation error is plotted in Fig. 5(a) against the number of qubits used in
the QN. Here we recorded 3 observables from each qubit in the QN at
$\tau=\pi/2\Gamma$, i.e., $\langle\sigma_{x}\rangle$,
$\langle\sigma_{y}\rangle$, and $\langle\sigma_{z}\rangle$, where
$\sigma_{\\{x,y,z\\}}$ stand for the Pauli matrices. Similar shift is seen
where $\Delta_{E}\sim 10^{-11}$ is obtained for a QN with at least 5 qubits.
This is because to fully characterise an $N$-qubit input state, one requires
$2^{2N}-1$ parameters. This way, to estimate entanglement of a two-qubit input
state ideally, at least 5 qubits are needed in the QN, corresponding to a
total of 15 observables. The direct comparison between the estimated and input
entanglement can be seen in Figs. 5(b) and (c) when the QN is composed of 4
and 5 qubits, respectively. We note that, in principle, if one were to record
one observable from each qubit in the QN, it would require at least 15 qubits,
which is too demanding to simulate on classical computers. In this case, we
show below that time-multiplexing is of help.

Figure 5: (a) Estimation error vs number of qubits used in the QN. The inset
shows the profile of entanglement of the randomly generated input states in
both training and testing. Panels (b) and (c) present explicit comparison
between the estimated and input entanglement during testing where the QN is
composed of $3$ and $4$ qubits, respectively.

We also performed simulations by taking into account measurement errors with
$\zeta=10^{-3}$. We present the estimation errors in Fig. 6(a) against the
number of qubits used in the QN. One can also utilise time-multiplexing with
only measurements of $\langle\hat{b}^{\dagger}_{m}\hat{b}_{m}\rangle$ from
each QN node. An example of this is plotted in Fig. 6(b), where the
measurements are performed three times at $\tau=\\{1,2,3\\}\pi/2\Gamma$ on
each QN node. One can see that both panels in Fig. 6 show error scaling beyond
the SQL (dashed curves) $\propto 1/\sqrt{\mathcal{N}}\propto
1/\sqrt{3Q}\propto 1/\sqrt{Q}$.

Figure 6: (a) Estimation error plotted against the number of qubits $Q$ in QN.
Panel (b), taking only mean excitation from each QN node as the observable
with time-multiplexing performed $\mathcal{T}=3$ times. The scaling for SQL is
given by the dashed curve in each panel.

##  Appendix B Hybrid systems

Here we show that the general scheme introduced in the main text is not
limited to particular quantum systems (only CV or discrete systems). In what
follows, we demonstrate this with a simple hybrid system: discrete input with
CV QN. It is important to note that as this involves discrete systems, the
dynamics will not preserve Gaussianity of the CV systems in general.
Covariance matrix does not fully describe the involved CV systems. Therefore,
we describe all quantum systems by their density matrices (truncated dimension
for CV systems at $d=20$).

Let us demonstrate sensing entanglement of two-qubit input with a single
bosonic mode as the QN. In particular, take the Hamiltonian as

$H=\hbar\sum_{n=1}^{2}\Delta_{n}\hat{a}_{n}^{\dagger}\hat{a}_{n}+\hbar\Lambda\hat{b}^{\dagger}\hat{b}+\hbar\sum_{n=1}^{2}K_{n}(\hat{a}_{n}\hat{b}^{\dagger}+\hat{b}\hat{a}_{n}^{\dagger})+\hbar
P(\hat{b}+\hat{b}^{\dagger}),$ (7)

where $\hat{a}_{n}$ denotes the lowering operator ($|g\rangle\langle e|$) for
the $n$th input qubit and $\hat{b}$ is the bosonic lowering operator for the
single QN node. For simplicity, we take the evolution as unitary, i.e.,
$\rho(t)=\hat{U}\rho(0)\hat{U}^{\dagger}$ with $\hat{U}=\exp(-iHt/\hbar)$. The
initial state is random $\rho_{\text{in}}$ for the qubits and vacuum
$|0\rangle$ for the QN node.

The parameters are randomised as
$\\{\Delta_{n},\Lambda,K_{n},2P\\}\in[1,2]\Gamma$. As the readout, we take the
mean excitation of the QN node $\langle\hat{b}^{\dagger}\hat{b}\rangle$ with
time-multiplexing at $\mathcal{T}$ different times, i.e.,
$\tau=\\{1,2,\cdots,\mathcal{T}\\}\pi/10\Gamma$. We tested this architecture
for $10$ different realisations of the parameters, where in each we performed
training with $N_{\text{tr}}=50$ and testing with $N_{\text{te}}=100$. Our
simulations show a transition to low entanglement estimation error
$\Delta_{E}\sim 10^{-8}$ for $\mathcal{T}\geq 15$. Again, this is because it
requires at least 15 different parameters to characterise two-qubit input
states.

Note that one can also apply the scheme to estimate entanglement of CV input
systems (with truncated dimension) using discrete systems as the QN. In
general, the following requirements set the guidelines for a particular setup
to be viable:

  1. 1.

The number of independent observables from the QN nodes has to be at least
equal to the number of independent parameters required to characterise the
state of the input objects.

  2. 2.

A dynamics ensuring that sufficient information about the inputs are carried
forward to the QN observables. This requires the essential interactions
between the input and QN nodes as well as within the QN nodes.

  3. 3.

The independent QN observables may be obtained from different QN nodes or/and
time-multiplexing (measurement of observables at different times). Normally,
for simpler QN observables such as mean excitations, relatively richer
dynamics is necessary. For CV systems, the latter can be achieved by, e.g.,
adding two-photon pumping for the QN nodes, having ultra-strong coupling
between the involved quantum systems, or even nonlinearity.

##  Appendix C Generic CV systems: Details

From the Hamiltonian of Eq. (1) in the main text, a set of LEs is obtained
from the equations of motion in Heisenberg picture and the addition of noise
terms:

$\displaystyle\dot{\hat{a}}_{n}$ $\displaystyle=$
$\displaystyle-(\gamma_{n}+i\Delta_{n})\hat{a}_{n}-i\sum_{m}K_{nm}\hat{b}_{m}+\sqrt{2\gamma_{n}}\;\hat{a}_{n}^{\text{in}},$
$\displaystyle\dot{\hat{b}}_{m}$ $\displaystyle=$
$\displaystyle-(\kappa_{m}+i\Lambda_{m})\hat{b}_{m}-i\sum_{n}K_{nm}\hat{a}_{n}-i\sum_{m^{\prime}}J_{mm^{\prime}}\hat{b}_{m^{\prime}}-iP_{m}+\sqrt{2\kappa_{m}}\;\hat{b}_{m}^{\text{in}},$
(8)

where $\hat{a}_{n}^{\text{in}}$ and $\hat{b}_{m}^{\text{in}}$ are zero mean
Gaussian noise operators with correlation functions
$\langle\hat{a}_{n}^{\text{in}}(t)\hat{a}_{n^{\prime}}^{\text{in},\dagger}(t^{\prime})\rangle=\delta_{nn^{\prime}}\delta(t-t^{\prime})$
and
$\langle\hat{b}_{m}^{\text{in}}(t)\hat{b}_{m^{\prime}}^{\text{in},\dagger}(t^{\prime})\rangle=\delta_{mm^{\prime}}\delta(t-t^{\prime})$
[31]. This allows us to write the LEs in terms of dimensionless position and
momentum quadratures

$\displaystyle\dot{\hat{q}}_{n}$ $\displaystyle=$
$\displaystyle-\gamma_{n}\hat{q}_{n}+\Delta_{n}\hat{r}_{n}+\sum_{m}K_{nm}\hat{p}_{m}+\sqrt{2\gamma_{n}}\;\hat{q}_{n}^{\text{in}}$
$\displaystyle\dot{\hat{r}}_{n}$ $\displaystyle=$
$\displaystyle-\gamma_{n}\hat{r}_{n}-\Delta_{n}\hat{q}_{n}-\sum_{m}K_{nm}\hat{x}_{m}+\sqrt{2\gamma_{n}}\;\hat{r}_{n}^{\text{in}}$
$\displaystyle\dot{\hat{x}}_{m}$ $\displaystyle=$
$\displaystyle-\kappa_{m}\hat{x}_{m}+\Lambda_{m}\hat{p}_{m}+\sum_{n}K_{nm}\hat{r}_{n}+\sum_{m^{\prime}}J_{mm^{\prime}}\hat{p}_{m^{\prime}}+\sqrt{2\kappa_{m}}\;\hat{x}_{m}^{\text{in}}$
$\displaystyle\dot{\hat{p}}_{m}$ $\displaystyle=$
$\displaystyle-\kappa_{m}\hat{p}_{m}-\Lambda_{m}\hat{x}_{m}-\sum_{n}K_{nm}\hat{q}_{n}-\sum_{m^{\prime}}J_{mm^{\prime}}\hat{x}_{m^{\prime}}-\sqrt{2}P_{m}+\sqrt{2\kappa_{m}}\;\hat{p}_{m}^{\text{in}}.$
(9)

We have used the following quadrature relations:

$\displaystyle\hat{q}_{n}$ $\displaystyle=$
$\displaystyle\frac{\hat{a}_{n}+\hat{a}_{n}^{\dagger}}{\sqrt{2}},\>\>\hat{r}_{n}=\frac{\hat{a}_{n}-\hat{a}_{n}^{\dagger}}{i\sqrt{2}},\>\>\hat{x}_{m}=\frac{\hat{b}_{m}+\hat{b}_{m}^{\dagger}}{\sqrt{2}},\>\>\hat{p}_{m}=\frac{\hat{b}_{m}-\hat{b}_{m}^{\dagger}}{i\sqrt{2}},$
$\displaystyle\hat{q}_{n}^{\text{in}}$ $\displaystyle=$
$\displaystyle\frac{\hat{a}_{n}^{\text{in}}+(\hat{a}_{n}^{\text{in}})^{\dagger}}{\sqrt{2}},\>\>\hat{r}_{n}^{\text{in}}=\frac{\hat{a}_{n}^{\text{in}}-(\hat{a}_{n}^{\text{in}})^{\dagger}}{i\sqrt{2}}\>\>\hat{x}_{m}^{\text{in}}=\frac{\hat{b}_{m}^{\text{in}}+(\hat{b}_{m}^{\text{in}})^{\dagger}}{\sqrt{2}},\>\>\hat{p}_{m}^{\text{in}}=\frac{\hat{b}_{m}^{\text{in}}-(\hat{b}_{m}^{\text{in}})^{\dagger}}{i\sqrt{2}}.$
(10)

The LEs in Eq. (9) can be written in a matrix form

$\dot{\bm{u}}(t)={\bm{A}}{\bm{u}}(t)+{\bm{h}}(t),$ (11)

where the vector

$\displaystyle\bm{u}=[\hat{q}_{1},\hat{r}_{1},\hat{q}_{2},\hat{r}_{2},\hat{x}_{1},\hat{p}_{1},\hat{x}_{2},\hat{p}_{2},\cdots,\hat{x}_{M},\hat{p}_{M}]^{T}$
(12)

contains the quadratures, the drift matrix reads

$\bm{A}=\left[\begin{array}[]{ccccccccccc}-\gamma_{1}&\Delta_{1}&0&0&0&K_{11}&0&K_{12}&\cdots&0&K_{1M}\\\
-\Delta_{1}&-\gamma_{1}&0&0&-K_{11}&0&-K_{12}&0&\cdots&-K_{1M}&0\\\
0&0&-\gamma_{2}&\Delta_{2}&0&K_{21}&0&K_{22}&\cdots&0&K_{2M}\\\
0&0&-\Delta_{2}&-\gamma_{2}&-K_{21}&0&-K_{22}&0&\cdots&-K_{2M}&0\\\
0&K_{11}&0&K_{21}&-\kappa_{1}&\Lambda_{1}&0&J_{11}&\cdots&0&J_{1M}\\\
-K_{11}&0&-K_{21}&0&-\Lambda_{1}&-\kappa_{1}&-J_{11}&0&\cdots&-J_{1M}&0\\\
0&K_{12}&0&K_{22}&0&J_{11}&-\kappa_{2}&\Lambda_{2}&\cdots&0&J_{2M}\\\
-K_{12}&0&-K_{22}&0&-J_{11}&0&-\Lambda_{2}&-\kappa_{2}&\cdots&-J_{2M}&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\
0&K_{1M}&0&K_{2M}&0&J_{1M}&0&J_{2M}&\cdots&-\kappa_{M}&\Lambda_{M}\\\
-K_{1M}&0&-K_{2M}&0&-J_{1M}&0&-J_{2M}&0&\cdots&-\Lambda_{M}&-\kappa_{M}\end{array}\right],$
(13)

and the vector

$\bm{h}=\left[\begin{array}[]{c}0\\\ 0\\\ 0\\\ 0\\\ 0\\\ -\sqrt{2}P_{1}\\\
0\\\ -\sqrt{2}P_{2}\\\ \vdots\\\ 0\\\
-\sqrt{2}P_{M}\end{array}\right]+\left[\begin{array}[]{c}\sqrt{2\gamma_{1}}\hat{q}_{1}^{\text{in}}\\\
\sqrt{2\gamma_{1}}\hat{r}_{1}^{\text{in}}\\\
\sqrt{2\gamma_{2}}\hat{q}_{2}^{\text{in}}\\\
\sqrt{2\gamma_{2}}\hat{r}_{2}^{\text{in}}\\\
\sqrt{2\kappa_{1}}\hat{x}_{1}^{\text{in}}\\\
\sqrt{2\kappa_{1}}\hat{p}_{1}^{\text{in}}\\\
\sqrt{2\kappa_{2}}\hat{x}_{2}^{\text{in}}\\\
\sqrt{2\kappa_{2}}\hat{p}_{2}^{\text{in}}\\\ \vdots\\\
\sqrt{2\kappa_{M}}\hat{x}_{M}^{\text{in}}\\\
\sqrt{2\kappa_{M}}\hat{p}_{M}^{\text{in}}\end{array}\right]$ (14)

contains the pump and noise terms. The solution to Eq. (11) is given by

${\bm{u}}(t)={\bm{W}}_{+}(t){\bm{u}}(0)+{\bm{W}}_{+}(t)\int_{0}^{t}dt^{\prime}{\bm{W}}_{-}(t^{\prime}){\bm{h}}(t^{\prime}),$
(15)

where ${\bm{W}}_{\pm}(t)=\exp{(\pm{\bm{A}}t)}$. One can then form the CM at
time $t$ and show that it follows

$\frac{d{\bm{V}}}{dt}={\bm{A}}{\bm{V}}(t)+{\bm{V}}(t){\bm{A}}^{T}+{\bm{D}},$
(16)

where
${\bm{D}}=\mbox{diag}[\gamma_{1},\gamma_{1},\gamma_{2},\gamma_{2},\kappa_{1},\kappa_{1},\kappa_{2},\kappa_{2},\cdots,\kappa_{M},\kappa_{M}]$.

In simulations, the initial CM of the system is taken as

$\bm{V}(0)=\left[\begin{array}[]{cc}\bm{V}_{\text{in}}(0)&\bm{0}\\\
\bm{0}&\bm{V}_{\text{qn}}(0)\end{array}\right],$ (17)

where $\bm{V}_{\text{in}}(0)$ is the initial CM for the input objects, whereas
that for the QN nodes is initiated with vacuum
$\bm{V}_{\text{qn}}(0)=\openone/2$.

##  Appendix D The generation of random input states for generic dynamics

For CV systems, the random two-mode input CMs are generated dynamically as
follows. We take a Hamiltonian of the form

$H=\hbar\Delta_{1}\hat{a}_{1}^{\dagger}\hat{a}_{1}+\hbar\Delta_{2}\hat{a}_{2}^{\dagger}\hat{a}_{2}+\hbar
K(\hat{a}_{1}\hat{a}_{2}^{\dagger}+\hat{a}_{2}\hat{a}_{1}^{\dagger})+\hbar
P^{\prime}_{1}(\hat{a}_{1}\hat{a}_{1}+\hat{a}_{1}^{\dagger}\hat{a}_{1}^{\dagger})+\hbar
P^{\prime}_{2}(\hat{a}_{2}\hat{a}_{2}+\hat{a}_{2}^{\dagger}\hat{a}_{2}^{\dagger}),$
(18)

where the two modes are coupled and pumped (two-photon drives). The drives
have terms similar to single mode squeezing operations. As the initial CM, we
take vacuum $\openone/2$. The parameters are taken as random, i.e.,
$(\Delta_{1},\Delta_{2},K,P^{\prime}_{1},P^{\prime}_{1})\in(1,1,1,0.3,0.3)[0,1]\Gamma$,
and the evolution time $\tau_{0}=\pi/2\Gamma$. The random input CMs used for
CV systems in the main text are sampled from $\bm{V}(\tau_{0})$. The
entanglement profile resulting from this distribution is plotted in the inset
of Fig. 2(a) in the main text. For the scheme where learning is done with non-
entangled states, the evolution time is taken as $\tau_{0}=1/2\Gamma$ with
thermal states (CM is $3\openone/2$) as the initial condition.

For discrete systems (qubits), the random input states are sampled as follows.

$\displaystyle Z/10$ $\displaystyle=$ $\displaystyle
2(\upsilon_{1}+i\upsilon_{2})-(1+i)\mathcal{J}+h.c.,$
$\displaystyle\rho_{\text{in}}$ $\displaystyle=$ $\displaystyle
ZZ^{\dagger}/\text{tr}(ZZ^{\dagger}),$ (19)

where $\upsilon_{1,2}$ is a random $4\times 4$ matrix whose elements are
sampled from standard normal distribution and $\mathcal{J}$ is a $4\times 4$
matrix of ones. This sampling results in entanglement profile shown in the
inset of Fig. 5(a). For generating classically correlated states, one can
simply destroy the entanglement in $\rho_{\text{in}}$ by projective
measurements on one input object, i.e., $\rho_{\text{in,
sep}}=\Pi_{1}\rho_{\text{in}}\Pi_{1}+\Pi_{2}\rho_{\text{in}}\Pi_{2}$ with
$\\{\Pi_{1},\Pi_{2}\\}$ as random projection operators and
$\Pi_{1}+\Pi_{2}=\openone$.

##  Appendix E Entanglement estimation for CV systems using mean excitations

In the main text, we have demonstrated that by taking 3 observables from each
QN node, the shift to low estimation error requires at least 4 QN nodes, see
Fig. 2(a). Here, we present the case where we utilise the mean excitation from
each QN node instead. In this case, we present two options both of which
require additional (necessary) ingredient. First, one can add two-photon pump
to each QN node (the two-photon pump can be relatively weaker in strength than
the single-photon pump). This is carried out by adding
$\sum_{m}P_{m}^{\prime}(\hat{b}_{m}\hat{b}_{m}+\hat{b}_{m}^{\dagger}\hat{b}_{m}^{\dagger})$
to the Hamiltonian of Eq. (1) in the main text. For simulations, we take
$P_{m}^{\prime}\in[0,1]\Gamma/10$. We present the estimation error in Fig.
7(a), where the shift to low estimation error is achieved for a QN having at
least 10 nodes. For option two, the interactions between the QN nodes are
taken following the ultra-strong coupling type. In this case, one simply
replaces the operator function in Eq. (1) with
$\mathcal{F}(\hat{b}_{m},\hat{b}_{m^{\prime}})=(\hat{b}_{m}+\hat{b}_{m}^{\dagger})(\hat{b}_{m^{\prime}}+\hat{b}_{m^{\prime}}^{\dagger})$.
Similarly, the estimation error is plotted in Fig. 7(b).

Similar to Fig. 3(b) in the main text, we present the estimation error in Fig.
7(c) for CV systems using ultra-strong coupling. The time-multiplexing is
$\mathcal{T}=3$ and the strength of measurement errors is $\zeta=10^{-3}$.

Figure 7: Performance of entanglement sensing for CV systems using mean
excitations $\\{\langle\hat{b}_{m}^{\dagger}\hat{b}_{m}\rangle\\}$: (a) with
the addition of two-photon pump or (b) the use of ultra-strong coupling.
Measurement errors are not taken into account in panels (a) and (b). (c) The
estimation error as in panel (b) with time-multiplexing at
$\tau=\\{1,2,3\\}\pi/2\Gamma$, taking into account measurement errors with
$\zeta=10^{-3}$. The dashed curve represents the SQL scaling.

We appreciate that the function of the quantum reservoir is to map the input
state (or its parameters) to the measured local observables. The mapping is
unknown, as the parameters defining the reservoir are allowed to be random.
However, one can assume that the measured local observables attained by the
mapping must contain sufficient information to be a representation of the
initial input state. Apparently, without two-photon pumping or ultra-strong
coupling information is not well spread, i.e., many distinct states are mapped
to the same QN state, and therefore, information is lost. The two-photon
pumping creates and destroy two excitations on QN nodes and ultra-strong
coupling ensures creation/annihilation of two excitations in addition to the
normal hopping type coupling between the QN nodes. Both processes result in
more states being populated, explore higher-dimensional subspace of the QN
Hilbert space, which makes the spread of information more complex.

##  Appendix F Learning with non-entangled states: Results

As described in the main text (Fig. 2) and Section A above, entanglement
sensing may be performed using generic CV or discrete systems. There, the
training and testing both use random entangled input states, with entanglement
profile given in the inset of Fig. 2(a) and Fig. 5(a). We present similar
analysis in Fig. 8 where the training only utilises separable input states.
For CV systems, panel (a1) shows the entanglement profile of the input states
used in testing and (b1) the comparison between estimated and input
entanglement for a QN composed of 4 nodes. Similarly, the case for discrete
systems are presented in panels (a2) and (b2) using a QN composed of 5 qubits.
One can see that panels (b1) and (b2) are similar to Fig. 2(c) in the main
text and Fig. 5(c), respectively.

Figure 8: Learning with non-entangled states. For CV systems: (a1)
entanglement profile of the input states in testing and (b1) estimated vs
input entanglement in testing using 4 QN nodes. The corresponding plots for
discrete systems are presented in panels (a2) and (b2) using 5 qubits in the
QN. All input states for training have zero entanglement $E_{\text{in}}=0$.

##  Appendix G Gravity-induced entanglement: Details

For the case without probes, the equations of motion in Heisenberg picture
read

$\displaystyle\dot{\hat{x}}_{A}$ $\displaystyle=$
$\displaystyle\omega\hat{p}_{A},\>\>\dot{\hat{x}}_{B}=\omega\hat{p}_{B},$
$\displaystyle\dot{\hat{p}}_{A}$ $\displaystyle=$
$\displaystyle-\omega(1-\eta)\hat{x}_{A}-\omega\eta\hat{x}_{B}-\gamma\hat{p}_{A}+\hat{\xi}_{A},$
$\displaystyle\dot{\hat{p}}_{B}$ $\displaystyle=$
$\displaystyle-\omega(1-\eta)\hat{x}_{B}-\omega\eta\hat{x}_{A}-\gamma\hat{p}_{B}+\hat{\xi}_{B},$
(20)

where $\eta\equiv 2Gm/\omega^{2}L^{3}$, $\gamma$ represents the damping of
each mass, and $\hat{\xi}_{A(B)}$ the Brownian-like noise for the masses. We
assume high mechanical quality factor, $\omega/\gamma\gg 1$, where the noises
can be treated as uncoloured noise with correlation function
$\langle\hat{\xi}_{j}(t)\hat{\xi}_{j}(t^{\prime})+\hat{\xi}_{j}(t^{\prime})\hat{\xi}_{j}(t)\rangle/2=\gamma(2\bar{n}+1)\delta(t-t^{\prime})$,
$j=\\{A,B\\}$ [40, 41]. The thermal phonon number $\bar{n}$ is related to the
temperature of the environment $T$ as
$\bar{n}=1/(\exp(\hbar\omega/k_{B}T)-1)$.

The matrix form of the LEs in (20) is written as
$\dot{\bm{u}}(t)={\bm{A}}{\bm{u}}(t)+{\bm{h}}(t)$, where

$\bm{u}=[\hat{x}_{A},\hat{p}_{A},\hat{x}_{B},\hat{p}_{B}]^{T},$ (21)
$\bm{A}=\left[\begin{array}[]{cccc}0&\omega&0&0\\\
-\omega(1-\eta)&-\gamma&-\omega\eta&0\\\ 0&0&0&\omega\\\
-\omega\eta&0&-\omega(1-\eta)&-\gamma\end{array}\right],$ (22)

and

$\bm{h}=[0,\hat{\xi}_{A},0,\hat{\xi}_{B}]^{T}.$ (23)

The solution to the quadratures $\bm{u}(t)$ and CM $\bm{V}(t)$ are obtained as
in Eqs. (15) and (16), where
$\bm{D}=\mbox{diag}[0,\gamma(2\bar{n}+1),0,\gamma(2\bar{n}+1)]$.

As the initial state for each mass, we take thermal squeezed state, i.e.,
$\bm{V}(0)=\mbox{diag}[e^{2r_{0}},e^{-2r_{0}},e^{2r_{0}},e^{-2r_{0}}](1+2\bar{n})/2$.
The system is evolved for a time $\tau_{0}$, after which one obtains
$\bm{u}(\tau_{0})$ and $\bm{V}(\tau_{0})$.

When the probes are turned on at $\tau_{0}$, the new dynamics is now described
as follows. From the Hamiltonian of Eq. (5) in the main text, the new LEs read

$\displaystyle\dot{\hat{x}}_{A}$ $\displaystyle=$
$\displaystyle\omega\hat{p}_{A},\>\>\dot{\hat{x}}_{B}=\omega\hat{p}_{B},$
$\displaystyle\dot{\hat{p}}_{A}$ $\displaystyle=$
$\displaystyle-\omega(1-\eta)\hat{x}_{A}-\omega\eta\hat{x}_{B}-\gamma\hat{p}_{A}+G_{0a}\hat{a}^{\dagger}\hat{a}+\hat{\xi}_{A},$
$\displaystyle\dot{\hat{p}}_{B}$ $\displaystyle=$
$\displaystyle-\omega(1-\eta)\hat{x}_{B}-\omega\eta\hat{x}_{A}-\gamma\hat{p}_{B}-G_{0b}\hat{b}^{\dagger}\hat{b}+\hat{\xi}_{B},$
$\displaystyle\dot{\hat{a}}$ $\displaystyle=$
$\displaystyle-(\kappa_{a}+i\Delta_{0a})\hat{a}+iG_{0a}\hat{a}\hat{x}_{A}+\mathcal{E}_{a}+\sqrt{2\kappa_{a}}\;\hat{a}^{\text{in}}$
$\displaystyle\dot{\hat{b}}$ $\displaystyle=$
$\displaystyle-(\kappa_{b}+i\Delta_{0b})\hat{b}-iG_{0b}\hat{b}\hat{x}_{B}+\mathcal{E}_{b}+\sqrt{2\kappa_{b}}\;\hat{b}^{\text{in}},$
(24)

where $\hat{a}^{\text{in}}$ and $\hat{b}^{\text{in}}$ are Gaussian noise
operators with
$\langle\hat{a}^{\text{in}}(t)\hat{a}^{\text{in},\dagger}(t^{\prime})\rangle=\delta(t-t^{\prime})$
and
$\langle\hat{b}^{\text{in}}(t)\hat{b}^{\text{in},\dagger}(t^{\prime})\rangle=\delta(t-t^{\prime})$
[31].

The linearised version of the LEs is obtained through the following
transformations

$\displaystyle\hat{x}_{A}$ $\displaystyle\rightarrow
x_{As}+\delta\hat{x}_{A},\>\hat{p}_{A}\rightarrow p_{As}+\delta\hat{p}_{A},$
$\displaystyle\hat{x}_{B}$ $\displaystyle\rightarrow
x_{Bs}+\delta\hat{x}_{B},\>\hat{p}_{B}\rightarrow p_{Bs}+\delta\hat{p}_{B},$
$\displaystyle\hat{a}$ $\displaystyle\rightarrow\alpha_{s}+\delta\hat{a},$
$\displaystyle\hat{b}$ $\displaystyle\rightarrow\beta_{s}+\delta\hat{b},$

and by ignoring any nonlinear term such as
$\delta\hat{a}^{\dagger}\delta\hat{a}$, $\delta\hat{a}\delta\hat{x}_{A}$,
$\delta\hat{b}^{\dagger}\delta\hat{b}$, and $\delta\hat{b}\delta\hat{x}_{B}$
in the fluctuation operators. In what follows, we will consider a much shorter
evolution time such that one may neglect the contribution from the
gravitational coupling ($\eta\approx 0$). In particular, we have

$\displaystyle\delta\dot{\hat{x}}_{A}$ $\displaystyle=$
$\displaystyle\omega\delta\hat{p}_{A},\>\>\delta\dot{\hat{x}}_{B}=\omega\delta\hat{p}_{B},$
$\displaystyle\delta\dot{\hat{p}}_{A}$ $\displaystyle=$
$\displaystyle-\omega\delta\hat{x}_{A}-\gamma\delta\hat{p}_{A}+G_{a}\delta\hat{x}_{a}+\hat{\xi}_{A},$
$\displaystyle\delta\dot{\hat{p}}_{B}$ $\displaystyle=$
$\displaystyle-\omega\delta\hat{x}_{B}-\gamma\delta\hat{p}_{B}-G_{b}\delta\hat{x}_{b}+\hat{\xi}_{B},$
$\displaystyle\delta\dot{\hat{x}}_{a}$ $\displaystyle=$
$\displaystyle-\kappa_{a}\delta\hat{x}_{a}+\Delta_{a}\delta\hat{p}_{a}+\sqrt{2\kappa_{a}}\;\hat{x}_{a}^{\text{in}},$
$\displaystyle\delta\dot{\hat{p}}_{a}$ $\displaystyle=$
$\displaystyle-\kappa_{a}\delta\hat{p}_{a}-\Delta_{a}\delta\hat{x}_{a}+G_{a}\delta\hat{x}_{A}+\sqrt{2\kappa_{a}}\;\hat{p}_{a}^{\text{in}}$
$\displaystyle\delta\dot{\hat{x}}_{b}$ $\displaystyle=$
$\displaystyle-\kappa_{b}\delta\hat{x}_{b}+\Delta_{b}\delta\hat{p}_{b}+\sqrt{2\kappa_{b}}\;\hat{x}_{b}^{\text{in}},$
$\displaystyle\delta\dot{\hat{p}}_{b}$ $\displaystyle=$
$\displaystyle-\kappa_{b}\delta\hat{p}_{b}-\Delta_{b}\delta\hat{x}_{b}-G_{b}\delta\hat{x}_{B}+\sqrt{2\kappa_{b}}\;\hat{p}_{b}^{\text{in}},$
(25)

where

$\displaystyle p_{As}$ $\displaystyle=$ $\displaystyle 0,\>\>p_{Bs}=0,$
$\displaystyle x_{As}$ $\displaystyle=$
$\displaystyle\frac{G_{0a}|\alpha_{s}|^{2}}{\omega},\>\>x_{Bs}=-\frac{G_{0b}|\beta_{s}|^{2}}{\omega},$
$\displaystyle\alpha_{s}$ $\displaystyle=$
$\displaystyle\frac{|\mathcal{E}_{a}|}{\sqrt{\kappa_{a}^{2}+\Delta_{a}^{2}}},\>\>\beta_{s}=\frac{|\mathcal{E}_{b}|}{\sqrt{\kappa_{b}^{2}+\Delta_{b}^{2}}}.$
(26)

In these equations, we have introduced the quantities
$\Delta_{a}=\Delta_{0a}-G_{0a}x_{As}$, $\Delta_{b}=\Delta_{0b}+G_{0b}x_{Bs}$,
$G_{a}=G_{0a}\alpha_{s}\sqrt{2}$, and $G_{b}=G_{0b}\beta_{s}\sqrt{2}$. Note
that $\alpha_{s}$ and $\beta_{s}$ have been assumed real, which can be done by
tuning the phase of the laser $\mathcal{E}_{a}$ and $\mathcal{E}_{b}$,
respectively. We have also used the following quadrature relations:

$\displaystyle\hat{x}_{a}$ $\displaystyle=$
$\displaystyle\frac{\hat{a}+\hat{a}^{\dagger}}{\sqrt{2}},\>\>\hat{p}_{a}=\frac{\hat{a}-\hat{a}^{\dagger}}{i\sqrt{2}},\>\>\hat{x}_{b}=\frac{\hat{b}+\hat{b}^{\dagger}}{\sqrt{2}},\>\>\hat{p}_{b}=\frac{\hat{b}-\hat{b}^{\dagger}}{i\sqrt{2}},$
$\displaystyle\hat{x}_{a}^{\text{in}}$ $\displaystyle=$
$\displaystyle\frac{\hat{a}^{\text{in}}+(\hat{a}^{\text{in}})^{\dagger}}{\sqrt{2}},\>\>\hat{p}_{a}^{\text{in}}=\frac{\hat{a}^{\text{in}}-(\hat{a}^{\text{in}})^{\dagger}}{i\sqrt{2}}\>\>\hat{x}_{b}^{\text{in}}=\frac{\hat{b}^{\text{in}}+(\hat{b}^{\text{in}})^{\dagger}}{\sqrt{2}},\>\>\hat{p}_{b}^{\text{in}}=\frac{\hat{b}^{\text{in}}-(\hat{b}^{\text{in}})^{\dagger}}{i\sqrt{2}}.$
(27)

One can write the LEs for the fluctuation of the quadratures in (25) as
$\dot{\bm{u}}(t)={\bm{A}}{\bm{u}}(t)+{\bm{h}}(t)$, where now

$\bm{u}=[\delta\hat{x}_{A},\delta\hat{p}_{A},\delta\hat{x}_{B},\delta\hat{p}_{B},\delta\hat{x}_{a},\delta\hat{p}_{a},\delta\hat{x}_{b},\delta\hat{p}_{b}]^{T},$
(28) $\bm{A}=\left[\begin{array}[]{cccccccc}0&\omega&0&0&0&0&0&0\\\
-\omega&-\gamma&0&0&G_{a}&0&0&0\\\ 0&0&0&\omega&0&0&0&0\\\
0&0&-\omega&-\gamma&0&0&-G_{b}&0\\\ 0&0&0&0&-\kappa_{a}&\Delta_{a}&0&0\\\
G_{a}&0&0&0&-\Delta_{a}&-\kappa_{a}&0&0\\\
0&0&0&0&0&0&-\kappa_{b}&\Delta_{b}\\\
0&0&-G_{b}&0&0&0&-\Delta_{b}&-\kappa_{b}\\\ \end{array}\right],$ (29)

and

$\bm{h}=[0,\hat{\xi}_{A},0,\hat{\xi}_{B},\hat{x}_{a}^{\text{in}},\hat{p}_{a}^{\text{in}},\hat{x}_{b}^{\text{in}},\hat{p}_{b}^{\text{in}}]^{T}.$
(30)

The solution $\bm{u}(t)$ and $\bm{V}(t)$ are obtained from Eqs. (15) and (16),
where
$\bm{D}=\mbox{diag}[0,\gamma(2\bar{n}+1),0,\gamma(2\bar{n}+1),\kappa_{a},\kappa_{a},\kappa_{b},\kappa_{b}]$.
As the initial CM, we take $\bm{V}(\tau_{0})$ for the masses and vacuum
$\openone/2$ for the cavity modes.

The parameters chosen in the main text (see the caption of Fig. 4) are
motivated as follows. Mechanical mirrors of mass $m\sim 1$ kg and frequency
$\omega\sim 0.1$ Hz have been cooled down near their ground state [42], see
also Ref. [43]. The initial state for each mass (squeezed thermal) can be
prepared by appropriate optical driving [44, 45], and the strength
$r_{0}=1.73$ is motivated by the squeezing of light mode [46] and advances in
the state transfer in optomechanics [38]. Cavity length ($25$ mm) and laser
wavelength $1064$ nm are typical in optomechanics [47], see also Ref. [38] for
cavity finesse up to $10^{5}$. For example, a cavity finesse $F=8\times
10^{4}$ gives a decay rate $\Gamma=\pi c/2FL_{a(b)}\approx 2.36\times 10^{5}$
Hz. This is used in simulations as a basis for the random cavity decay rates
and effective detunings, i.e.,
$\\{\kappa_{a},\kappa_{b},\Delta_{a},\Delta_{b}\\}\in[1,2]\Gamma.$

##  Appendix H Gravity-induced entanglement: Error scaling

Figure 9 presents the entanglement estimation error against time-multiplexing
instances $\mathcal{T}$ for the case of Fig. 4(b) with $\tau_{0}=5$ s. The
estimation error is taken as standard deviation,

$\delta_{E}=\sqrt{\sum_{l^{\prime}}^{N_{\text{te}}}\frac{(E_{\text{est},l^{\prime}}-E_{\text{in}})^{2}}{N_{\text{te}}-1}}.$
(31)

One can see that the scaling of the estimation error is beyond the SQL (dashed
curve). In fact, it follows a Heisenberg-like scaling with
$\delta_{E}\propto\mathcal{T}^{-1}$ (dashed-dotted curve). We note that for
$\mathcal{T}=4$, the estimation error is comfortably below $10^{-4}$, which is
two orders of magnitude lower than what was experimentally achieved in Ref.
[29].

Figure 9: Error scaling. The standard deviation for the setup in Fig. 4(b) in
the main text for $\tau_{0}=5$ s vs instances of time-multiplexing
$\mathcal{T}$.

## References

  * Ching  _et al._ [2018] T. Ching, D. S. Himmelstein, B. K. Beaulieu-Jones, A. A. Kalinin, B. T. Do, G. P. Way, E. Ferrero, P.-M. Agapow, M. Zietz, M. M. Hoffman, _et al._ , Opportunities and obstacles for deep learning in biology and medicine, Journal of The Royal Society Interface 15, 20170387 (2018).
  * Topol [2019] E. J. Topol, High-performance medicine: the convergence of human and artificial intelligence, Nature Medicine 25, 44 (2019).
  * Hannun  _et al._ [2019] A. Y. Hannun, P. Rajpurkar, M. Haghpanahi, G. H. Tison, C. Bourn, M. P. Turakhia, and A. Y. Ng, Cardiologist-level arrhythmia detection and classification in ambulatory electrocardiograms using a deep neural network, Nature Medicine 25, 65 (2019).
  * Mehta  _et al._ [2019] P. Mehta, M. Bukov, C.-H. Wang, A. G. Day, C. Richardson, C. K. Fisher, and D. J. Schwab, A high-bias, low-variance introduction to machine learning for physicists, Physics Reports 810, 1 (2019).
  * Fujii and Nakajima [2017] K. Fujii and K. Nakajima, Harnessing disordered-ensemble quantum dynamics for machine learning, Physical Review Applied 8, 024030 (2017).
  * Ghosh  _et al._ [2021a] S. Ghosh, K. Nakajima, T. Krisnanda, K. Fujii, and T. C. H. Liew, Quantum Neuromorphic Computing with Reservoir Computing Networks, Advanced Quantum Technologies 4, 2100053 (2021a).
  * Montavon  _et al._ [2012] G. Montavon, G. Orr, and K.-R. Müller,  _Neural networks: tricks of the trade_ , Vol. 7700 (springer, 2012).
  * Govia  _et al._ [2021] L. C. G. Govia, G. J. Ribeill, G. E. Rowlands, H. K. Krovi, and T. A. Ohki, Quantum reservoir computing with a single nonlinear oscillator, Physical Review Research 3, 013077 (2021).
  * Xu  _et al._ [2021] H. Xu, T. Krisnanda, W. Verstraelen, T. C. H. Liew, and S. Ghosh, Superpolynomial quantum enhancement in polaritonic neuromorphic computing, Physical Review B 103, 195302 (2021).
  * Ghosh  _et al._ [2019a] S. Ghosh, A. Opala, M. Matuszewski, T. Paterek, and T. C. H. Liew, Quantum reservoir processing, npj Quantum Information 5, 35 (2019a).
  * Ghosh  _et al._ [2020] S. Ghosh, A. Opala, M. Matuszewski, T. Paterek, and T. C. H. Liew, Reconstructing Quantum States With Quantum Reservoir Networks, IEEE Transactions on Neural Networks and Learning Systems 32, 3148 (2020).
  * Ghosh  _et al._ [2019b] S. Ghosh, T. Paterek, and T. C. H. Liew, Quantum Neuromorphic Platform for Quantum State Preparation, Physical Review Letters 123, 260404 (2019b).
  * Krisnanda  _et al._ [2021] T. Krisnanda, S. Ghosh, T. Paterek, and T. C. H. Liew, Creating and concentrating quantum resource states in noisy environments using a quantum neural network, Neural Networks 136, 141 (2021).
  * Ghosh  _et al._ [2021b] S. Ghosh, T. Krisnanda, T. Paterek, and T. C. H. Liew, Realising and compressing quantum circuits with quantum reservoir computing, Communications Physics 4, 105 (2021b).
  * Krisnanda  _et al._ [2022] T. Krisnanda, S. Ghosh, T. Paterek, W. Laskowski, and T. C. H. Liew, Phase Measurement Beyond the Standard Quantum Limit Using a Quantum Neuromorphic Platform, Physical Review Applied 18, 034011 (2022).
  * Marković and Grollier [2020] D. Marković and J. Grollier, Quantum neuromorphic computing, Applied Physics Letters 117, 150501 (2020).
  * Horodecki  _et al._ [2009] R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki, Quantum entanglement, Reviews of Modern Physics 81, 865 (2009).
  * Krisnanda  _et al._ [2017] T. Krisnanda, M. Zuppardo, M. Paternostro, and T. Paterek, Revealing Nonclassicality of Inaccessible Objects, Physical Review Letters 119, 120402 (2017).
  * Bose  _et. al._ [2017] S. Bose  _et. al._ , Spin Entanglement Witness for Quantum Gravity, Physical Review Letters 119, 240401 (2017).
  * Marletto and Vedral [2017] C. Marletto and V. Vedral, Gravitationally Induced Entanglement between Two Massive Particles is Sufficient Evidence of Quantum Effects in Gravity, Physical Review Letters 119, 240402 (2017).
  * Krisnanda  _et al._ [2020] T. Krisnanda, G. Y. Tham, M. Paternostro, and T. Paterek, Observable quantum entanglement due to gravity, npj Quantum Information 6, 12 (2020).
  * Al Balushi  _et al._ [2018] A. Al Balushi, W. Cong, and R. B. Mann, Optomechanical quantum cavendish experiment, Physical Review A 98, 043811 (2018).
  * Belenchia  _et al._ [2018] A. Belenchia, R. M. Wald, F. Giacomini, E. Castro-Ruiz, Č. Brukner, and M. Aspelmeyer, Quantum superposition of massive objects and the quantization of gravity, Physical Review D 98, 126009 (2018).
  * Qvarfort  _et al._ [2020] S. Qvarfort, S. Bose, and A. Serafini, Creating esoscopic entanglement through central potential interactions, Journal of Physics B 53, 235501 (2020).
  * van de Kamp  _et al._ [2020] T. W. van de Kamp, R. J. Marshman, S. Bose, and A. Mazumdar, Quantum gravity witness via entanglement of masses: Casimir screening, Physical Review A 102, 062807 (2020).
  * Rijavec  _et al._ [2021] S. Rijavec, M. Carlesso, A. Bassi, V. Vedral, and C. Marletto, Decoherence effects in non-classicality tests of gravity, New Journal of Physics 23, 043040 (2021).
  * Margalit  _et al._ [2021] Y. Margalit, O. Dobkowski, Z. Zhou, O. Amit, Y. Japha, S. Moukouri, D. Rohrlich, A. Mazumdar, S. Bose, C. Henkel, _et al._ , Realization of a complete stern-gerlach interferometer: Toward a test of quantum gravity, Science advances 7, eabg2879 (2021).
  * Pedernales  _et al._ [2022] J. S. Pedernales, K. Streltsov, and M. B. Plenio, Enhancing gravitational interaction between quantum systems by a massive mediator, Physical Review Letters 128, 110401 (2022).
  * Palomaki  _et al._ [2013] T. Palomaki, J. Teufel, R. Simmonds, and K. W. Lehnert, Entangling mechanical motion with microwave fields, Science 342, 710 (2013).
  * Adesso  _et al._ [2014] G. Adesso, S. Ragy, and A. R. Lee, Continuous variable quantum information: Gaussian states and beyond, Open Systems & Information Dynamics 21, 1440001 (2014).
  * Walls and Milburn [2007] D. F. Walls and G. J. Milburn,  _Quantum optics_ (Springer Science & Business Media, 2007).
  * Vidal and Werner [2002] G. Vidal and R. F. Werner, Computable measure of entanglement, Physical Review A 65, 032314 (2002).
  * Lambert  _et al._ [2013] N. Lambert, Y.-N. Chen, Y.-C. Cheng, C.-M. Li, G.-Y. Chen, and F. Nori, Quantum biology, Nature Physics 9, 10 (2013).
  * Scholes  _et al._ [2017] G. D. Scholes, G. R. Fleming, L. X. Chen, A. Aspuru-Guzik, A. Buchleitner, D. F. Coker, G. S. Engel, R. Van Grondelle, A. Ishizaki, D. M. Jonas, _et al._ , Using coherence to enhance function in chemical and biophysical systems, Nature 543, 647 (2017).
  * Collini  _et al._ [2010] E. Collini, C. Y. Wong, K. E. Wilk, P. M. Curmi, P. Brumer, and G. D. Scholes, Coherently wired light-harvesting in photosynthetic marine algae at ambient temperature, Nature 463, 644 (2010).
  * Panitchayangkoon  _et al._ [2010] G. Panitchayangkoon, D. Hayes, K. A. Fransted, J. R. Caram, E. Harel, J. Wen, R. E. Blankenship, and G. S. Engel, Long-lived quantum coherence in photosynthetic complexes at physiological temperature, Proceedings of the National Academy of Sciences 107, 12766 (2010).
  * Krisnanda  _et al._ [2018] T. Krisnanda, C. Marletto, V. Vedral, M. Paternostro, and T. Paterek, Probing quantum features of photosynthetic organisms, npj Quantum Information 4, 60 (2018).
  * Aspelmeyer  _et al._ [2014] M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, Cavity optomechanics, Reviews of Modern Physics 86, 1391 (2014).
  * Paternostro  _et al._ [2007] M. Paternostro, D. Vitali, S. Gigan, M. Kim, C. Brukner, J. Eisert, and M. Aspelmeyer, Creating and probing multipartite macroscopic entanglement with light, Physical Review Letters 99, 250401 (2007).
  * Giovannetti and Vitali [2001] V. Giovannetti and D. Vitali, Phase-noise measurement in a cavity with a movable mirror undergoing quantum Brownian motion, Physical Review A 63, 023812 (2001).
  * Benguria and Kac [1981] R. Benguria and M. Kac, Quantum Langevin equation, Physical Review Letters 46, 1 (1981).
  * Abbott  _et al._ [2009] B. Abbott, R. Abbott, R. Adhikari, P. Ajith, B. Allen, G. Allen, R. Amin, S. Anderson, W. Anderson, M. Arain, _et al._ , Observation of a kilogram-scale oscillator near its quantum ground state, New Journal of Physics 11, 073032 (2009).
  * Whittle  _et al._ [2021] C. Whittle, E. D. Hall, S. Dwyer, N. Mavalvala, V. Sudhir, R. Abbott, A. Ananyeva, C. Austin, L. Barsotti, J. Betzwieser, _et al._ , Approaching the motional ground state of a 10-kg object, Science 372, 1333 (2021).
  * Vanner  _et al._ [2013] M. Vanner, J. Hofer, G. Cole, and M. Aspelmeyer, Cooling-by-measurement and mechanical state tomography via pulsed optomechanics, Nature Communications 4, 2295 (2013).
  * Rashid  _et al._ [2016] M. Rashid, T. Tufarelli, J. Bateman, J. Vovrosh, D. Hempston, M. Kim, and H. Ulbricht, Experimental Realization of a Thermal Squeezed State of Levitated Optomechanics, Physical Review Letters 117, 273601 (2016).
  * Vahlbruch  _et al._ [2016] H. Vahlbruch, M. Mehmet, K. Danzmann, and R. Schnabel, Detection of 15 dB Squeezed States of Light and their Application for the Absolute Calibration of Photoelectric Quantum Efficiency, Physical Review Letters 117, 110801 (2016).
  * Gröblacher  _et al._ [2009] S. Gröblacher, K. Hammerer, M. R. Vanner, and M. Aspelmeyer, Observation of strong coupling between a micromechanical resonator and an optical cavity field, Nature 460, 724 (2009).


* 37 Royston P, Parmar M. Restricted mean survival time: an alternative to the hazard ratio for the design and analysis of randomized trials with a time-to-event outcome. BMC Medical Research Methodology 2013; 152: 152.
  * 38 Cox D. Regression Models and Life-Tables. Journal of the Royal Statistical Society (Series B) 1972; 34: 187–220.
  * 39 Aalen O. A linear regression model for the analysis of life times. Statistics in Medicine 1989; 8: 907–925.
  * 40 Aalen O, Borgan Ø, Gjessing H. Survival and Event History Analysis: A Process Point of View. New York: Springer .  2008\.
  * 41 VanderWeele T. Concerning the consistency assumption in causal inference. Epidemiology 2009; 20: 880-883.
  * 42 Daniel R, Cousens S, De Stavola B, Kenward M, Sterne J. Methods for dealing with time-dependent confounding. Statistics in Medicine 2013; 32: 1584–1618.
  * 43 D’Agostino R, Lee ML, Belanger A, Cupples L, Anderson K, Kannel W. Relation of pooled logistic regression to time-dependent Cox regression analysis: The Framingham Heart Study. Statistics in Medicine 1990; 9: 1501–1515.
  * 44 Røysland K, Gran J, Ledergerber B, Wyl vV, Young J, Aalen O. Analyzing direct and indirect effects of treatment using dynamic path analysis applied to data from the Swiss HIV Cohort. Statistics in Medicine 2011; 30: 2947–2958.
  * 45 Richardson T, Rotnitzky A. Causal Etiology of the Research of James M. Robins. Statisticl Science 2014; 29: 459-484.
  * 46 Hernán M, Robins J. Causal Inference: What If. Boca Raton: Chapman and Hall/CRC .  2010\.
  * 47 Morris T, White I, Crowther M. Using simulation studies to evaluate statistical methods. Statistics in Medicine 2019; 38: 2074–2102.
  * 48 McKone E, Goss C, Aitken M. CFTR Genotype as a Predictor of Prognosis in Cystic Fibrosis. Chest 2006; 130: 1141–1147.
  * 49 Cole S, Edwards J, Hall H, et al. Incident AIDS or Death After Initiation of Human Immunodeficiency Virus Treatment Regimens Including Raltegravir or Efavirenz Among Adults in the United States. Clinical Infectious Diseases 2017; 64: 1591–1596.
  * 50 van der Laan M, Petersen M, Joffe M. History-adjusted marginal structural models and statically-optimal dynamic treatment regimens.. International Journal of Biostatistics 2005; 1: article 4.
  * 51 Bembom O, van der Laan M, Haight T, Tager I. Leisure-time Physical Activity and All-cause Mortality in an Elderly Cohort. Epidemiology 2009; 20: 424–430.
  * 52 Schaubel D, Wolfe R, Port F. A sequential stratification method for estimating the effect of a time-dependent experimental treatment in observational studies.. Biometrics 2006; 62: 910–917.
  * 53 Seaman S, Dukes R, Vansteelandt S. Adjusting for time‐varying confounders in survival analysis using structural nested cumulative survival time models. Biometrics 2019\.
  * 54 Sawicki G, Signorovitch J, Zhang J, Latremouille-Viau M, Wu E, Shi L. Reduced mortality in cystic fibrosis patients treated with tobramycin inhalation solution. Pediatric Pulmonology 2012; 47: 44–52.
  * 55 Cole S, Hernán M. Constructing inverse probability weights for marginal structural models. American Journal of Epidemiology 2008; 168: 656–664.

Causal inference in survival analysis using longitudinal observational data:
Sequential trials and marginal structural models

Supplementary Material

##  A1 Inverse probability weights

###  A1.1 MSM-IPTW

To estimate MSMs using IPTW, the weight at time $t$ for individual $i$ is the
inverse of their probability of their observed treatment pattern up time time
$t$ given their time-dependent covariate history 2, 1:

$W_{i}(t)=\prod_{k=0}^{\lfloor
t\rfloor}\frac{1}{\Pr(A_{k}=A_{k,i}|\bar{L}_{k,i},\bar{A}_{k-1,i},T\geq k)}$
(A1)

Some individuals can have very large weights, which can results in the
parameters of the MSM being estimated very imprecisely, and therefore
stabilized weights are typically used. The stabilized weight for individual
$i$ is:

$SW_{i}(t)=\prod_{k=0}^{\lfloor
t\rfloor}\frac{\Pr(A_{k}=A_{k,i}|\bar{A}_{k-1,i},T\geq
k)}{\Pr(A_{k}=A_{k,i}|\bar{L}_{k,i},\bar{A}_{k-1,i},T\geq k)}$ (A2)

The model in the numerator of the stabilized weights can also include
covariates measured at time 0, $\bar{L}_{0}$:

$SW_{i}(t)=\prod_{k=0}^{\lfloor
t\rfloor}\frac{\Pr(A_{k}=A_{k,i}|\bar{A}_{k-1,i},\bar{L}_{0,i},T\geq
k)}{\Pr(A_{k}=A_{k,i}|\bar{L}_{k,i},\bar{A}_{k-1,i},T\geq k)},$ (A3)

in which case the MSM must also condition on these variables, as noted in the
main text.

The models used in the weights can be estimated using logistic regression, and
they may be fitted using all visits combined or separately by visit.

###  A1.2 Sequential trials

For individuals in the trial starting at visit $k$ but who do not initiate
treatment at that time ($\bar{A}_{k}=0$) the weight at time $t-k$ from the
start of the trial is

$W^{0}_{i,k}(t-k)=\prod_{j=k+1}^{k+\lfloor
t\rfloor}\frac{1}{\Pr(A_{j}=0|\bar{L}_{j,i},\bar{A}_{(k,j-1),i}=0,T\geq k)}$
(A4)

where $\bar{A}_{(k,j-1)}=\\{A_{k},A_{k+1},\ldots,A_{j-1}\\}$ denotes the
treatment status from visit $k$ to visit $j-1$. For those who initiate
treatment at the start of trial $k$ ($\bar{A}_{k-1}=0,A_{k}=1$) the weight is

$W^{1}_{i,k}(t-k)=\prod_{j=k+1}^{k+\lfloor
t\rfloor}\frac{1}{\Pr(A_{j}=1|\bar{L}_{j,i},\bar{A}_{k-1,i}=0,\bar{A}_{(k,j-1),i}=1,T\geq
k)}$ (A5)

In the trial starting at visit $k$, the weights are equal to 1 for all
individuals between time $k$ and $k+1$, i.e. up to follow-up time 1 in trial
$k$.

The corresponding stabilized weights are

$SW^{0}_{i,k}(t-k)=\prod_{j=k+1}^{k+\lfloor
t\rfloor}\frac{\Pr(A_{j}=0|\bar{L}_{k,i},\bar{A}_{(k,j-1),i}=0,T\geq
k)}{\Pr(A_{j}=0|\bar{L}_{j,i},\bar{A}_{(k,j-1),i}=0,T\geq k)}$ (A6)

and

$SW^{1}_{i,k}(t-k)=\prod_{j=k+1}^{k+\lfloor
t\rfloor}\frac{\Pr(A_{j}=1|\bar{L}_{k,i},\bar{A}_{k-1,i}=0,\bar{A}_{(k,j-1),i}=1,T\geq
k)}{\Pr(A_{j}=1|\bar{L}_{j,i},\bar{A}_{k-1,i}=0,\bar{A}_{(k,j-1),i}=1,T\geq
k)}$ (A7)

where the probability in the numerator is conditional on the covariates
$\bar{L}_{k}$ observed at the start of trial $k$.

The models used in the weights can be estimated using logistic regression, and
they may be fitted using all visits combined or separately by visit, and
across all trials combined or separately by trial.

##  A2 Assumptions

The no interference assumption is that the counterfactual event time for a
given individual, $T^{\underline{a}_{0}}$, does not depend on the treatment
received by any other individuals. The positivity assumption is that each
individual has a strictly non-zero probability of receiving each given pattern
of treatments over time. Consistency means that an individual’s observed
outcome is equal to the counterfactual outcome when the assigned treatment
pattern is set to that which was actually received,
$T_{i}=T_{i}^{\underline{A}_{0,i}}$. The conditional exchangeability
assumption can be expressed formally as
$T^{\bar{A}_{k-1},\underline{a}_{k}}\perp\\!\\!\\!\perp
A_{k}|\bar{A}_{k-1},\bar{L}_{k},T\geq k$ for all feasible $\underline{a}_{k}$,
where $T^{\bar{A}_{k-1},\underline{a}_{k}}$ denotes the counterfactual event
time had an individual followed their observed treatment pattern up to time
$k-1$, $\bar{A}_{k-1}$, and had their treatments been set to
$\underline{a}_{k}$ from time $k$ onwards, given survival to time $k$. The
conditional exchangeability assumption means that among individuals who remain
at risk of the event at time $k$, the treatment $A_{k}$ received at time $k$
may depend on past treatment and covariates $\bar{A}_{k-1}$ and $\bar{L}_{k}$,
but that, conditional on these, it does not depend on the remaining lifetime
that would apply if all future treatments were set to any particular values
$\underline{a}_{k}$.

##  A3 Compatibility of MSMs used in the two approaches

To obtain expressions (28) and (29) in the main text we use the result that
the probability $\Pr(Y_{2}^{\bar{A}_{1}=a}=0|Y_{1}^{A_{0}=a}=0)$ can be
written

$\begin{split}\Pr(Y_{2}^{\bar{A}_{1}=a}=0|&Y_{1}^{A_{0}=a}=0)=\sum_{l_{0}}\Pr(Y_{2}^{\bar{A}_{1}=a}=0|Y_{1}^{A_{0}=a}=0,L_{0}=l_{0})\frac{\Pr(Y_{1}=0|A_{0}=a,L_{0}=l_{0})\Pr(L_{0}=l_{0})}{\sum_{l_{0}^{\prime}}\Pr(Y_{1}=0|A_{0}=a,L_{0}=l_{0}^{\prime})\Pr(L_{0}=l_{0}^{\prime})}\end{split}$
(A8)

and the first term in the sum can then be written

$\begin{split}\Pr(Y_{2}^{\bar{A}_{1}=a}=0|&Y_{1}^{A_{0}=a}=0,L_{0})=\sum_{l_{1}}\Pr(Y_{2}=0|Y_{1}=0,\bar{A}_{1}=a,L_{0},L_{1}=l_{1})\Pr(L_{1}=l_{1}|Y_{1}=0,A_{0}=a,L_{0}).\end{split}$
(A9)

Table A1: Simulation results. Summary of number of individuals observed at each time, number always treated and never treated, and corresponding percentages. Numbers and percentages are the mean across the 1000 simulations. |  | MSM-IPTW |  | Sequential trials
---|---|---|---|---
|  |  Time $k$ |  |  Time $k$
|  | 0 | 1 | 2 | 3 | 4 |  | 0 | 1 | 2 | 3 | 4
Scenario 1 |  |  |  |  |  |
Observed at time $k$ | N | 1000 | 829 | 694 | 588 | 502 |  | 2264 | 1414 | 873 | 509 | 238
| % | 100 | 83 | 69 | 59 | 50 |  | 100 | 62 | 39 | 22 | 10
Always treated from time 0 to $k$ | N | 279 | 237 | 204 | 177 | 155 |  | 643 | 514 | 397 | 282 | 155
| % | 28 | 29 | 29 | 30 | 31 |  | 28 | 36 | 46 | 55 | 65
Never treated from time 0 to $k$ | N | 721 | 424 | 249 | 144 | 82 |  | 1621 | 900 | 475 | 227 | 82
| % | 72 | 51 | 36 | 25 | 16 |  | 72 | 64 | 54 | 45 | 35
Scenario 2 |  |  |  |  |  |
Observed at time $k$ | N | 1000 | 821 | 675 | 556 | 458 |  | 3180 | 2178 | 1403 | 806 | 349
| % | 100 | 82 | 68 | 56 | 46 |  | 100 | 68 | 44 | 25 | 11
Always treated from time 0 to $k$ | N | 53 | 45 | 39 | 34 | 29 |  | 197 | 142 | 98 | 61 | 29
| % | 5 | 5 | 6 | 6 | 6 |  | 6 | 7 | 7 | 8 | 8
Never treated from time 0 to $k$ | N | 947 | 731 | 560 | 425 | 319 |  | 2983 | 2036 | 1305 | 745 | 319
| % | 95 | 89 | 83 | 76 | 70 |  | 94 | 93 | 93 | 92 | 92
Scenario 3 |  |  |  |  |  |
Observed at time $k$ | N | 1000 | 832 | 701 | 596 | 510 |  | 2083 | 1333 | 865 | 543 | 287
| % | 100 | 83 | 70 | 60 | 51 |  | 100 | 64 | 42 | 26 | 14
Always treated from time 0 to $k$ | N | 387 | 326 | 278 | 241 | 210 |  | 700 | 563 | 446 | 337 | 210
| % | 39 | 39 | 40 | 40 | 41 |  | 34 | 42 | 52 | 62 | 73
Never treated from time 0 to $k$ | N | 613 | 352 | 213 | 129 | 77 |  | 1383 | 770 | 419 | 206 | 77
| % | 61 | 42 | 30 | 22 | 15 |  | 66 | 58 | 48 | 38 | 27
Figure A1: Simulation results: bias in estimation of the risk difference using
the sequential trials analysis and the MSM-IPTW analysis. The black line shows
the bias at each time point and the grey area shows the Monte-Carlo 95% CI at
each time point. The MSM-IPTW results are from the MSM not conditional on
$L_{0}$. Figure A2: Simulation results: relative efficiency of the sequential
trials analysis compared with the MSM-IPTW analysis, defined as the ratio of
the empirical variances of the risk difference estimates at each time. The
MSM-IPTW results are from the MSM not conditional on $L_{0}$. Figure A3:
Simulation results using truncated weights: bias in estimation of the risk
difference using the sequential trials analysis and the MSM-IPTW analysis. The
black line shows the bias at each time point and the grey area shows the
Monte-Carlo 95% CI at each time point. The MSM-IPTW results are from the MSM
conditional on $L_{0}$.

.

Figure A4: Simulation results using truncated weights: relative efficiency of
the sequential trials analysis compared with the MSM-IPTW analysis, defined as
the ratio of the empirical variances of the risk difference estimates at each
time. The MSM-IPTW results are from the MSM conditional on $L_{0}$. Figure A5:
Plots showing the distribution of the weights used in the MSM-IPTW and
sequential trials analyses of the application using UK CF Registry data.

## References

  * Cole and Hernan 2008 Cole, S. and Hernan, M. (2008).  Constructing inverse probability weights for marginal structural models.  American Journal of Epidemiology 168, 656–664.
  * Daniel et al. 2013 Daniel, R., Cousens, S., De Stavola, B., Kenward, M., and Sterne, J. (2013).  Methods for dealing with time-dependent confounding.  Statistics in Medicine 32, 1584–1618.