# Parameter estimation by learning quantum correlations in continuous photon-
counting data using neural networks

Enrico Rinaldi  Quantinuum K.K., Otemachi Financial City Grand Cube 3F, 1-9-2
Otemachi, Chiyoda-ku, Tokyo, Japan  Interdisciplinary Theoretical and
Mathematical Sciences (iTHEMS) Program, RIKEN, Wakoshi, Saitama 351-0198,
Japan  Center for Quantum Computing, RIKEN, Wakoshi, Saitama 351-0198, Japan
Theoretical Quantum Physics Laboratory, Cluster for Pioneering Research,
RIKEN, Wakoshi, Saitama 351-0198, Japan  Manuel González Lastre  Departamento
de Física Teórica de la Materia Condensada and Condensed Matter Physics Center
(IFIMAC), Universidad Autónoma de Madrid, 28049 Madrid, Spain  Instituto
Nicolás Cabrera, Universidad Autónoma de Madrid, 28049 Madrid, Spain  Sergio
García Herreros  Departamento de Física Teórica de la Materia Condensada and
Condensed Matter Physics Center (IFIMAC), Universidad Autónoma de Madrid,
28049 Madrid, Spain  Instituto Nicolás Cabrera, Universidad Autónoma de
Madrid, 28049 Madrid, Spain  Shahnawaz Ahmed  Department of Microtechnology
and Nanoscience, Chalmers University of Technology, 412 96 Gothenburg, Sweden
Maryam Khanahmadi  Department of Microtechnology and Nanoscience, Chalmers
University of Technology, 412 96 Gothenburg, Sweden  Franco Nori  Center for
Quantum Computing, RIKEN, Wakoshi, Saitama 351-0198, Japan  Theoretical
Quantum Physics Laboratory, Cluster for Pioneering Research, RIKEN, Wakoshi,
Saitama 351-0198, Japan  Physics Department, University of Michigan, Ann
Arbor, MI 48109-1040, USA  Carlos Sánchez Muñoz<EMAIL_ADDRESS>Departamento de Física Teórica de la Materia Condensada and Condensed Matter
Physics Center (IFIMAC), Universidad Autónoma de Madrid, 28049 Madrid, Spain
Instituto Nicolás Cabrera, Universidad Autónoma de Madrid, 28049 Madrid, Spain

###### Abstract

We present an inference method utilizing artificial neural networks for
parameter estimation of a quantum probe monitored through a single continuous
measurement. Unlike existing approaches focusing on the diffusive signals
generated by continuous weak measurements, our method harnesses quantum
correlations in discrete photon-counting data characterized by quantum jumps.
We benchmark the precision of this method against Bayesian inference, which is
optimal in the sense of information retrieval. By using numerical experiments
on a two-level quantum system, we demonstrate that our approach can achieve a
similar optimal performance as Bayesian inference, while drastically reducing
computational costs. Additionally, the method exhibits robustness against the
presence of imperfections in both measurement and training data. This approach
offers a promising and computationally efficient tool for quantum parameter
estimation with photon-counting data, relevant for applications such as
quantum sensing or quantum imaging, as well as robust calibration tasks in
laboratory-based settings.

##  I Introduction

Quantum parameter estimation refers to the fundamental problem of inferring
the value of unknown physical quantities within a quantum system [1]. It has
important practical applications in device characterization [2, 3], quantum
feedback and control [4, 5, 6, 7, 8, 9] and quantum communications [10, 11],
and most importantly, it is at the core of the field of quantum metrology [12,
13, 14, 15, 16, 17, 18, 19, 20], whereby quantum correlations are exploited to
estimate physical quantities with high precision.

The most conventional approach to quantum parameter estimation is based on the
following steps: i) prepare a quantum system in a given state; ii) allow it to
evolve through a unitary transformation that encodes the unknown parameters;
iii) perform a measurement on the system; iv) estimate the parameter from the
measurement results; v) repeat the process to gather statistical information
and increase the precision of the estimation [14, 16]. Beyond this standard
prepare-evolve-measure approach, there are relevant situations where, instead,
one performs a single-shot continuous quantum measurement on a quantum probe
[1]. In this case, one is interested in obtaining an estimation that will be
continuously updated as data is gathered over the time evolution of a single
quantum trajectory. Protocols of quantum parameter estimation from continuous
measurements have been developed and studied over the years [21, 22, 23, 24,
25, 26, 22, 27, 28, 29, 30, 31, 32, 33, 34]; such parameter estimation
techniques are particularly relevant for applications such as atomic
magnetometry [35, 36], quantum metrology with open quantum systems [27, 37,
38] or estimation of classical time-dependent signals [39]. Previous works
have established computational approaches to calculate the probability density
$P(\theta|D)$ [22, 28, 29, 30, 31] of unknown $n$-dimensional parameters
$\theta=[\theta_{1},\theta_{2},\ldots,\theta_{n}]$ conditioned on the
continuously recorded data, $D$, which could correspond, for instance, to
photon-counting records [29, 30], continuous weak measurements [40, 41, 42],
or homodyne measurements [31]. These calculations rely on the application of
Bayes’ theorem on the likelihood functions $P(D|\theta)$, obtained as the
trace of un-normalized conditional density matrices with the evolution
conditioned on the measured data record [43].

This Bayesian approach provides, in principle, the most accurate description
of the knowledge that can be gained about $\theta$ given the measurement data
$D$. However, the application of this technique faces several difficulties.
Firstly, it requires exact modelling of the open quantum system and any source
of error, which can become challenging for noisy, complex quantum systems
[44]. Secondly, the method can be computationally demanding. For instance,
when tackling the problem of multi-parameter estimation, where one wishes to
estimate several parameters simultaneously, Bayesian inference calls for the
use of stochastic methods to sample the posterior probability distribution
[28]. Generally, the Bayesian parameter estimation method will be
computationally time-consuming, even in simple systems. This hampers the
prospects for real-time estimation and for the integration of the inference
process in the actual measurement device taking the data, which would be
desirable for reduced latency and energy consumption [45, 46].

A novel approach to this problem that has gained attention in recent years is
the application of machine-learning methods in which a neural network (NN) is
trained to aid in the inference process. Crucially, this method allows to
circumvent the requirement for detailed knowledge of the underlying physical
system, as long as a reasonable amount of training data is available. The use
of NNs have been demonstrated for a variety of tasks in quantum metrology,
including parameter estimation in the presence of unavoidable noise [47, 48,
49, 50, 51], calibration [52, 53], phase estimation in interferometry [52, 53,
54], magnetometry [55, 51], control [50], and tomography [56, 57]. In the
context of continuously monitored systems, NNs have been employed for the
reconstruction of quantum dynamics [58, 2] and parameter estimation under
continuous weak dispersive measurements [59, 47, 48, 55], which are
characterized by a continuous, diffusive evolution of the quantum system [1,
60].

Here, we assess the potential of deep learning to tackle the parameter
estimation task from continuous measurements of the photon-counting type.
These measurements are described by point processes that, rather than yielding
a diffusive evolution of the system, are characterized by sudden quantum jumps
upon the acquisition of a discrete signal, such as the detection of a single
photon [1]. This type of measurements, ubiquitous in the field of quantum
optics and cavity-QED [61, 62, 63, 64, 65, 66, 67, 68], provide insights into
the particle-like nature of quantum systems and can produce strongly
correlated signals that directly violate classical inequalities, as in the
case of the renowned effect of photon antibunching [68]. Our work addresses
the question of whether NNs can effectively process and exploit the quantum
correlations present in this type of discrete signals.

Our proposed NN architectures approach this task as a regression problem,
taking the time delays between individual photodetection events as an input
and providing an estimate $\hat{\theta}$ for the desired parameters of
interest as the output, as sketched in Fig. 1. Our main finding is that the
NNs can take advantage of the quantum correlations within the signal,
providing higher precision than using an equivalent signal that lacks these
correlations, e.g., the integrated intensity value. By computing the root-
mean-squared error (RMSE) of the estimations, we find that, after training,
the NN can reach the limit of precision established by Bayesian inference,
while requiring dramatically fewer computational resources. We also
demonstrate that parameter estimation with NNs could be more resilient to
noise than a simple Bayesian approach with no detailed knowledge of the
underlying noise sources incorporated into the model, since it is possible for
NNs to indirectly learn the noise model during the training phase. In contrast
to a full Bayesian posterior distribution, the architectures studied in this
work provide a single point estimate of the parameters, without including a
confidence interval. Nevertheless, we evaluate their efficiency as estimators
by conducting statistical analyses of their performance across numerous data
samples.

Figure 1: Quantum parameter estimation strategies in open quantum systems.
Parameters are encoded in the dynamics of an open quantum system: here, these
are the frequency detuning between the qubit transition frequency and the
drive frequency, $\Delta=\omega_{q}-\omega_{L}$, and the Rabi frequency
$\Omega$ of the coherent drive, proportional to the dipole moment of the qubit
and the amplitude of the electromagnetic field driving it. Photon-counting
measurements of the quantum light radiated yield a record of the times of
photodetection. The unknown parameters can be reconstructed through Bayesian
parameter estimation (left). Here, we show the posterior probability
distribution $P(\theta|D)$ for the data obtained from a single quantum
trajectory, from which an estimator $\hat{\theta}$ is obtained by computing
the mean value of the distribution. The ground truth value for the parameters
is contained within the posterior distribution in a region of high
probability. An alternative approach is based on the use of neural networks
(right). The neural network learns the inverse process that maps observed data
to the parameters $x\to\text{NN}\to\hat{\theta}$ through an initial training
phase that employs data generated with known parameters. Once trained,
inference for new data is much faster using the neural network than repeating
the Bayesian parameter estimation with the new data. Figure 2: Classical and
quantum data. (a) Dynamics of an individual Monte Carlo quantum jump
trajectory for $\Delta=0.8\gamma$. Top panel shows the posterior probability
of the parameter $\Delta$ for the photon-counting data $D(t)$ acquired up to
time $t$ from a certain trajectory, starting from a flat prior. Bottom panel
shows the corresponding evolution of the population of the qubit for the same
trajectory, with red grid lines marking the times when a quantum jump occurs,
leading to the collapse of the population to zero and subsequently restarting
the system from the $\mathinner{|{0}\rangle}$ state. The list of time delays
between jumps constitutes the quantum data $D$ from which we infer the unknown
parameters $\theta$. (b) Waiting-time distribution $w(\tau)$ of the measured
time delays $\tau_{i}$ between photo-detection events during an individual
trajectory, which constitute the quantum signal. The plot shows the
distribution for several $\Delta$ parameters of the system. Clear patters
emerge which are $\Delta$-dependent. A classical signal corresponding to the
integrated intensity is defined as the mean of the measured delays,
$\bar{\tau}=\sum_{i=1}^{N}\tau_{i}/N$. Since all $\tau_{i}$ are independent
random variables coming from the same underlying distribution $w(\tau)$,
$\mathrm{E}[\bar{\tau}]=\mathrm{E}[\tau_{i}]=\mathrm{E}[\tau]$, shown as a
black-dashed line. The red vertical gridline marks the value
$\Delta=0.8\gamma$ used as ground truth in panel (a). The horizontal grid line
marks the corresponding expected value $\mathrm{E}[\tau](\Delta=0.8\gamma)$
(c) Heatmap of normalized histograms of $\tau$ obtained from $10^{4}$
trajectories with $N=48$ clicks each, simulated with $\Delta=0.8\gamma$. The
oscillatory patterns of the waiting-time distribution are clearly reproduced
in the histograms, which represent the quantum data. (d) Histogram of all the
values $\bar{\tau}$ obtained for the $10^{4}$ trajectories used in panel (c).
This represents the classical data, given by the mean of delays measured in a
given trajectory. The resulting Gaussian distribution is centered around
$\mathrm{E}[\tau]$, marked by the dashed red gridline. We fix $\Omega=\gamma$
in all panels. All units are defined in terms of $\gamma$.

##  II Results

###  II.1 System and data generation

The data we use to train the deep neural networks and perform parameter
inference consists of simulated measurements of single photo-detection times
from the light emitted by a driven two-level system (TLS). This system
represents a paradigmatic example of the types of quantum parameter estimation
problems studied extensively in the literature [21, 22, 27]. It is a
prototypical model that can describe various quantum emitters, such as quantum
dots [69], molecules [70], or color centres [71]. We consider a TLS with
states $\\{|0\rangle,|1\rangle\\}$, lowering operator
$\hat{\sigma}=|0\rangle\langle 1|$ and transition frequency $\omega_{q}$,
continuously excited by a coherent drive with frequency $\omega_{L}$ and a
Rabi frequency $\Omega$. This quantum probe is coupled to the environment,
giving rise to spontaneous emission with a rate $\gamma$. In the frame
rotating at the drive frequency, the dynamics of such a system is described by
the master equation (we set $\hbar=1$ herein):

$\partial_{t}\hat{\rho}=-i\left[\Delta\hat{\sigma}^{\dagger}\hat{\sigma}+\Omega(\hat{\sigma}+\hat{\sigma}^{\dagger}),\hat{\rho}\right]+\frac{\gamma}{2}D(\hat{\sigma})[\hat{\rho}],$
(1)

where $\Delta\equiv\omega_{q}-\omega_{L}$ and $D(\hat{x})[\hat{\rho}]\equiv
2\hat{x}\hat{\rho}\hat{x}^{\dagger}-\left\\{\hat{x}^{\dagger}\hat{x},\hat{\rho}\right\\}$
follows the definition of the Lindblad superoperator [72]. The interplay
between continuous driving and losses eventually brings the system into a
steady state.

The estimation is performed using the data $D$ generated by the continuous
photon-counting measurement of the radiation emitted by the system during a
single run of the experiment, considering an initial state
$|\psi(0)\rangle=|0\rangle$ and a total evolution time, $T$, much longer than
the relaxation time towards the steady state. We simulate this process using
the Monte Carlo method of quantum trajectories [73, 1, 43] implemented in the
QuTiP library [74, 75]. Under continuous photon-counting measurements, the
system evolves stochastically under the following rules:

  * •

at every differential time step, $dt$, the system can experience a quantum
jump $|\psi(t+dt)\rangle\propto\hat{\sigma}|\psi(t)\rangle$ with a probability
$p(t)=dt\gamma\langle\hat{\sigma}^{\dagger}\hat{\sigma}\rangle(t)$, which
corresponds to the detection of a photon emitted by the system;

  * •

otherwise, the system follows a non-Hermitian evolution
$|\psi(t+dt)\rangle\propto[\mathbf{I}-i(\hat{H}-i\hat{\sigma}^{\dagger}\hat{\sigma}/2)dt]|\psi(t)\rangle$.
The wavefunction is normalized in each time step.

Each stochastic simulation of a trajectory corresponds to an individual
realization of the experiment, generating a data vector by storing the
specific set of times $(t_{1},\ldots,t_{N})$ when a photodetection leading to
a quantum jump occurs during the course of a total evolution time $T$, with
$N$ being the total number of jumps recorded. For convenience, we define our
data as the set of time delays between these jumps,
$D=[\tau_{1},\ldots,\tau_{N}]$, where $\tau_{i}=t_{i}-t_{i-1}$, $i=1,\ldots N$
and $t_{0}=0$. An important feature of this particular system is that each
time interval $\tau_{i}$ is an independent random variable, given that the
system collapses to the ground state $|0\rangle$ after each emission [29].

Each trajectory is simulated for an evolution time $T$ chosen to generate a
data record $D$ of precisely $N=48$ jumps. This number is chosen to meet a
compromise; it is sufficiently large to ensure that the system is in a steady
state, while also allowing for inference with a limited number of photons in a
relatively short timescale spanning just a few tens of emitter lifetimes,
which is typically within the range of nanoseconds for solid-state quantum
emitters [62]. It is important to note that fixing $N$ requires that each
trajectory will have a different evolution time $T$ [76, 77]. More detailed
information about the unravelling of the dynamics in quantum trajectories can
be found in the Supplemental Material 111See Supplemental Material for further
details on the theory of quantum jump trajectories, the choice of different
estimators from posterior distributions, and the Fisher information and the
CRB.

After obtaining the data $D$ from a single trajectory or experiment, our goal
is to estimate the unknown parameters $\theta$ that govern the dynamics of the
system. In this study, we specifically focus on the parameters $\Delta$ and/or
$\Omega$, considering that $\gamma$ is known. While this might be an idealized
setting, it is instructive enough to benchmark the power of quantum
correlations and neural networks for the task of parameter estimation. We are
interested in benchmarking the NNs against the ultimate precision limit set by
the Bayesian inference process. This process yields a posterior probability
distribution for the unknown parameters conditioned on the measured data $D$:
following Bayes’ rule we have $P(\theta|D)={P(D|\theta)}/{\int d\theta
P(D|\theta)}$. A comprehensive explanation of Bayesian inference from
continuous measurements is provided in Ref. [27] and, within the context of
this work, in Ref. [78]. An illustrative individual trajectory, generated
through our simulation approach, is depicted in Fig. 2(a), which shows the
evolution of the posterior probability distribution for a single parameter,
namely $\theta=[\Delta]$, in the top panel, and the TLS population in the
bottom panel. Red grid lines mark the times when a quantum jump takes place.
One can see how the regions of high probability in the posterior distribution
concentrate more around the ground truth as time increases and more data is
collected.

Once the posterior distribution is obtained, we build an estimator by taking
the mean, $\hat{\theta}=E[\theta|D]=\int d\theta P(\theta|D)\theta$. While
reducing the full posterior to a single estimation $\hat{\theta}$ implies a
loss of information and involves a subtle choice of the best estimator, this
step allows us to straightforwardly compare the result with the estimation
provided by the NNs. Here we also note that we use a flat (uniform) prior on
the parameters of interest, while in general the Bayesian framework allows us
to introduce our knowledge of the system using arbitrary priors.

A central question of our study is whether NNs can leverage the quantum
correlations in the data. The emission from a TLS is strongly correlated; the
fact that the emitter can only store one excitation at a time gives rise to
the well-known non-classical phenomenon of antibunching [68], i.e., the
suppressed probability of detecting two photons simultaneously. These quantum
correlations can be evidenced through several observables, such as the
waiting-time distribution $w(\tau)$ (the probability of two consecutive
emissions to be separated by a time delay $\tau$), which presents a vanishing
value at zero delay, $w(\tau=0)=0$, and an oscillatory behaviour that depends
strongly on the system parameters. These features, shown in Fig. 2(b), are
translated directly into the histograms of time delays observed in each
individual trajectory, as we show in Fig. 2(c). We will refer to a signal
containing these time delays as a quantum signal.

To assess the information contained in these correlations, we compare it to an
equivalent classical signal where the only relevant value is the integrated
intensity, i.e. the total number of counts detected per unit time, $I=N/T$. We
can interpret such a classical signal as measurements lacking the necessary
time resolution to detect individual photons, providing an integrated
intensity without any access to the internal temporal correlations between
photon emissions. In fact, we choose the inverse of the intensity as the
classical signal, since it conveniently represents that we have reduced the
list of measured delays, that constitute the quantum signal, by taking their
average, $\bar{\tau}\equiv\sum_{i=1}^{N}\tau_{i}/N=I^{-1}$. Since $\bar{\tau}$
is a sum or independent random variables $\tau_{i}$ that all follow the
probability distribution $w(\tau_{i})$, the central limit theorem ensures
that, when $N$ is large, $\bar{\tau}$ is itself a random variable that follows
instead a normal distribution, with $\text{E}[\bar{\tau}]=\text{E}[\tau_{i}]$.
This normal distribution stands in contrast to the more complicated
probability distribution $w(\tau_{i})$ followed by the individual random
variables $\tau_{i}$. In the steady state, this distribution is centered
around $\text{E}[\tau_{i}]=\text{E}[\tau]\equiv\int d\tau
w(\tau)\tau=(\langle\hat{\sigma}^{\dagger}\hat{\sigma}\rangle_{\mathrm{ss}}\gamma)^{-1}$
[see Fig. 2(d)].

Access to $\text{E}[\tau]$ instead of $w(\tau)$ represents an obvious loss of
information, where any contributions from quantum correlations is completely
discarded. Thanks to this, the RMSE of the estimation made by using the
classical signal serves as a reference for comparison when doing estimations
with the full quantum signal: whenever the obtained RMSE is lower than that of
the classical signal, we understand that the estimation method is taking some
advantage from quantum correlations. The process of parameter estimation from
this classical signal, described in more detail the Methods section, also
follows a Bayesian approach in which the chosen estimator $\hat{\theta}$ is
the mean of the posterior $P(\theta|\bar{\tau})$.

###  II.2 Training

We considered different neural-network architectures with the common structure
of inputs and outputs that we illustrate in Fig. 1. The input consists of a
vector $x=D=[\tau_{1},\tau_{2},\ldots,\tau_{N}]$ of time delays between
photodetection events, while the output is a vector of size $n$ corresponding
to the estimations of the $n$-dimensional estimation problem,
$y=[\hat{\theta}_{1},\hat{\theta}_{1},\ldots,\hat{\theta}_{n}]$. Here, we
treat both one-dimensional estimation problems, with $\theta=[\Delta]$ in the
1D case, and $\theta=[\Delta,\Omega]$ in the 2D case. We fixed the size of the
input vector to $N=48$, indicating that our networks perform estimations based
on a single experiment run once 48 photons have been detected.

Estimation via NNs requires a training stage akin to a calibration process
(learning phase). In this work, the input data used for training,
$x_{\mathrm{train}}$, consists of a list of time delays obtained from
trajectories simulated in QuTiP, with the corresponding target data
$y_{\mathrm{train}}$ being the ground truth parameters used to run the
simulations. Our best results are achieved using two types of architectures:

  1. 1.

_RNN_ : A recurrent neural network (RNN) with two Long Short-Term Memory
(LSTM) layers. These networks can learn correlations in sequential data,
rendering them particularly promising for exploiting temporal quantum
correlations. It is worth noting that there are no such correlations between
time delays in our present system, since the $\tau_{i}$ are independent due to
the collapse following a quantum jump. Nonetheless, these networks still
deliver a very good performance. More information on these networks and their
application in parameter estimation from continuous signals can be found in
[55].

  2. 2.

_Hist-Dense_ : A fully connected neural network with a first custom layer
which incorporates the physical knowledge of the absence of correlations
between time delays by performing a histogram of the input data. The geometry
of the bins of the histogram is defined by the number of bins and their range.
We set these as non-trainable hyperparameters fixed at 700 total bins and a
range given by $\tau_{\mathrm{min}}=0$ and $\tau_{\mathrm{max}}=100/\gamma$.

Figure 3:  Performance of different estimation strategies. (a) Scatter plot of
estimations $\hat{\Delta}$ across a range of $\Delta$ values from classical
(left) and quantum (right) versions of the photon-counting signal, computed
for $10^{4}$ trajectories for each $\Delta$ value. All the different
estimators used for quantum data (Bayesian, NN) give a similar result, so only
one is shown (NN-Hist). The blue shade of the points represents the underlying
probability density function obtained through kernel density estimation. The
dashed red line represents the ground truth. Orange points represent the mean
of the the estimations. (b) RMSE of the estimators for on classical (red) and
quantum (yellow, blue, green) versions of the signal. Estimations on quantum
signals clearly outperform the classical estimation, signaling the useful
information encoded in quantum correlations. The predictions made by NNs
(blue, green) are on par with Bayesian inference (yellow), which is the
optimal process in terms of information retrieval. Dashed lines mark the lower
limit set by the biased Cramér-Rao bound for the classical and quantum
signals. (c) Bias of the predictions for the same cases as in (b).  Figure 4:
Multidimensional parameter estimation for $\theta=[\Delta,\Omega]$. The qubit-
drive detuning $\Delta$ and the Rabi frequency $\Omega$ are estimated
simultaneously. (a) RMSE of the estimations made through Bayesian inference.
(b) RMSE of the estimations made by the NN. (c) Contours corresponding to
Bayesian (red) and NN (blue) are overlaid for the sake of comparison.

Further details on the networks and the training process can be found in the
Methods section. In all the numerical calculations, the unit of time is set by
$1/\gamma$ (which we assume to be a known parameter) by fixing $\gamma=1$.
After training, we benchmark the estimations made by the NNs against the
estimation made through Bayesian inference and, in the single-parameter
estimation case, also against the estimation done on the classical version of
the signal. Details of this validation process are in the Methods section.

###  II.3 Single-parameter estimation (1D)

First, we focus on the simplest situation in which a single parameter is
estimated: the detuning $\theta=[\Delta]$. Figure 3(a) depicts a scatter plot
of the predictions made by each method for different values of $\Delta$. These
plots show that the predictions made using the quantum signal—regardless of
whether Bayesian inference or NNs are used—are more densely concentrated
around the ground truth than the predictions obtained using the classical
version of the signal. This already suggests that quantum correlations in the
data can provide an advantage.

This advantage is clearly quantified in Fig. 3(b), which compares the RMSE of
the predictions obtained by each of the methods computed on a different set of
validation trajectories (see Methods for details on the validation data). The
estimation from classical data (red) always has a higher RMSE than estimation
methods that have access to the quantum correlations in the data. Figure 3(b)
also shows that the NN architectures used in this work learn to extract
information encoded in these correlations, with the resulting RMSE reaching
values very close to the ultimate limit set by Bayesian inference. The fact
that neural networks perform as well as the Bayesian approach by leveraging on
quantum correlations in a completely model-agnostic manner is the main result
of our work.

We note, however, that none of the approaches studied here provided an
estimator that is unbiased for all $\Delta$, as we show in Fig. 3(c). There,
we plot the bias of several estimators, defined as
$\text{bias}(\hat{\theta})\equiv\text{E}(\hat{\theta}-\theta)$. This bias is
not problematic _per se_ , since it is known that biased estimators can be
preferable to unbiased ones in some cases [79]. However, this bias needs to be
taken into account at the time of comparing the obtained RMSE with the lower
bound set by the Fisher information and the Cramér-Rao bound (CRB) [80, 81].
As we discuss in more detail in Ref. [78], one needs to consider the _biased_
CRBs [82]. We calculated the Fisher information and the biased CRBs for the
photon-counting measurements considered here [28, 27, 29], incorporating the
calculated biases shown in Fig. 3(c). The resulting lower limits for the RMSE
are shown as dashed lines in Fig. 3(b). We show that, for the quantum signal,
both the Bayesian and NN estimators saturate the CRB for a wide range of
values of $\Delta$.

Figure 5: Parameter estimation with noisy data. (a) RMSE of NN predictor
trained and evaluated with data containing jitter noise given by a Gaussian
distribution with standard deviation $\sigma_{\tau}$ (green circles). This is
compared to the RMSE of Bayesian inference lacking a description of the noise
(yellow crosses). Results shown for a fixed $\Delta=0.27\gamma$. (b) Same as
in panel (a), shown versus $\Delta$ and with a fixed jitter noise
$\sigma_{\tau}=0.76\gamma$. RMSE of Bayesian predictor in the absence of noise
shown for comparison (gray crosses). (c) RMSE of NN predictor trained with
noisy target training data $y_{\mathrm{train}}$, with noise following a
Gaussian distribution with standard deviation $\sigma_{y}$. We observe that
the RMSE of the NN predictor (green circles) can remain smaller than the
standard deviation of the training data $\sigma_{y}$ (gray dot-dashed line).
The horizontal grid line marks the RMSE of Bayesian inference with noiseless
data for comparison . Results shown for a fixed $\Delta=0.27\gamma$. (d) Sames
as in panel (e), shown as a function of $\Delta$. Different colors represent
different values of $\sigma_{y}$, which are also represented as horizontal
grid lines, confirming that the RMSE of NN predictor remains smaller than
$\sigma_{y}$ for all values of $\Delta$.

###  II.4 Multi-parameter estimation (2D)

We also benchmarked our results for the 2D parameter-estimation case, where
both $\Delta$ and $\Omega$ are unknown and estimated simultaneously,
$\theta=[\Delta,\Omega]$.

Our results for the 2D estimation case are shown in Fig. 4. In panel (a), we
depict the RMSE of Bayesian estimation from the quantum signal in a region of
the plane $[\Delta,\Omega]$. The resulting map shows a nontrivial structure
that reflects the underlying physical complexity of the system and that is
consistent with previous reports in the literature [28, 29]. The same physical
complexity is captured by the NNs, as shown in Fig. 4(b), which depicts the
RMSE obtained by the estimation of a _Hist-Dense_ architecture and shows that
the RMSE reflects the same underlying structure as the one obtained with
Bayesian inference, yielding comparably small values overall. Both results can
be more clearly compared in Fig. 4(c), depicting overlaying contour lines of
the two plots in panels (a) and (b).

##  III Discussion

Our results clearly establish that NNs are able to exploit the information
hidden in quantum correlations in order to perform parameter estimation. Given
their comparable precision to Bayesian inference, a careful evaluation of the
benefits of opting for NNs as an alternative becomes important.

###  III.1 Computation efficiency

A first aspect in which NNs can bring an important advantage is the efficiency
of the computation, which is particularly relevant for the case in which one
wishes to achieve real-time estimation. For the 2D case of simultaneous
estimation of two parameters, the posterior for a single trajectory for
Bayesian estimation can be computed using the UltraNest package [83] (see
Methods) in a typical wall time of 10 seconds on a MacBook Pro 2023 with an
Apple M2 Max Chip. On the other hand, an estimation using the NNs built in
this work, while being less informative (unlike the Bayesian posterior they
provide a single estimation without a confidence interval) can be computed in
a wall time of less than 1 ms, on the same device. The efficiency of the
required calculation is such that trained NNs can be converted into compact
and efficient TensorFlow Lite models, with a size of approximately 90kB. These
models are suitable for deployment and execution on microcontrollers, which
would minimize resource demands and enable data collection and estimation on
the same device, leading to a significant reduction in latency [45, 46]. We
can expect that, during data collection of photon-counting measurements,
different chunks of measured time-delays are fed into the trained NN models
for inference in real time. Such a situation would not only allow for precise
and robust real time inference, but could also inform about time-drifting
parameters by looking at the temporal sequence of predictions.

###  III.2 The challenge of noise

Another important challenge faced by the Bayesian inference approach is the
requirement of a precise modelling of the underlying physical system. This
become particularly problematic in the presence of noise, which not only
requires a perfect characterization and modelling of the different sources of
imperfections, but also makes the subsequent calculation of posteriors much
more involved. In contrast, the model-agnostic approach of using NNs allows
for a robust prediction without the need of any modelling of the sources of
noise, only by training the network directly on noisy data. In the particular
case of photon-counting with single-photon detectors, numerous sources of
noise can be attributed just to the detection process, such as the efficiency,
bandwidth, dead time, dark counts or jitter noise [67]. All these sources of
noise have non-trivial consequences on the statistics of the photodetection
events [84] and on the quantum trajectories associated to these imperfect
measurements [85]. Hence, the necessity of a perfect modelling and
characterization of these processes seriously hampers the utilization of
Bayesian inference.

###  III.3 Noise in $x$: time jitter

We illustrate the potential advantage brought in by the estimation with NNs by
considering the particular case of time jitter [84], which we describe by
adding a noise term to the values of time delays $\tau$ used both in the
training and estimation of the NNs for the 1D estimation case
$\theta=[\Delta]$. We consider noise following a normal distribution with
standard deviation $\sigma_{\tau}$, so that $x\rightarrow
x+x_{\mathrm{noise}}(\sigma_{\tau})$, with
$x_{\mathrm{noise}}(\sigma_{\tau})\sim\mathcal{N}(0,\sigma_{\tau})$. In order
to assess the resilience of the NN predictors to this type of noise, we
defined a list of values of $\sigma_{\tau}$ in the range
$\sigma_{\tau}=[0,\gamma^{-1}]$. For each value in this list, we sampled a
population of $x_{\mathrm{noise}}$ and trained a model with the corresponding
noisy data.

The resulting RMSE of the predictions returned by each of these models is
computed using the same validation trajectories used in Fig. 3, but containing
an added Gaussian noise with the same $\sigma_{\tau}$ as the data used to
train the model. The resulting RMSEs are shown in Fig. 5(a), compared to the
RMSE of the prediction via an imperfect Bayesian inference process in which
the noise is not explicitly modeled. We observe that, in the case of NNs, the
results are more insensitive to the jitter noise and the RMSE of the NN
predictor is always lower than that of the Bayesian predictor. This result
remains true for all values of $\Delta$, as shown in Fig. 5(b). We emphasize
that this did not require any modelling of the imperfections, as would have
been required for a correct process of Bayesian inference. The latter would
presumably recover the values reported by the NN, but requiring a perfect
characterization of the sources of noise, which may be challenging in some
scenarios.

###  III.4 Noise in $y_{\mathrm{train}}$: imperfect calibration

Another relevant question is the effect of noise in the training target data
$y_{\mathrm{train}}$. In particular, we consider the case in which the target
data used for training differs from the underlying ground truth used to
generate the quantum trajectories by a Gaussian noise term with standard
deviation $\sigma_{y}$, so that $y_{\mathrm{train}}\rightarrow
y_{\mathrm{train}}+y_{\mathrm{noise}}(\sigma_{y})$, with
$y_{\mathrm{noise}}\sim\mathcal{N}(0,\sigma_{y})$. This could reflect, for
instance, the finite precision of a calibration device employed to
experimentally measure the set of target values used to train the network.

The question we address is whether a NN will be able to perform parameter
estimation with a lower RMSE than the standard deviation $\sigma_{y}$ of the
noisy target data used for training. A positive answer to this question would
indicate that the precision of the NN estimator can outperform that of the
calibration method or device used to generate the training data. Figures 5(c)
and 5(d) illustrate that this scenario is indeed possible. Similarly to the
process described for Figure 5(a), we defined a series of values of
$\sigma_{y}$ and, for each value, we sampled a population of
$y_{\mathrm{noise}}$ and trained a model with noisy target data. The
validation of the trained models is done using the same set of trajectories as
in Fig. 3. Figure 5(c) shows that, once the standard deviation of the noisy
training target data $\sigma_{y}$ (dot-dashed line) becomes greater than the
RMSE achievable in the noiseless case (horizontal grid line), the predictions
made by the NNs trained with this noisy data remain robust and approximately
equal to the ideal RMSE. As $\sigma_{y}$ increases, the RMSE of the NN
predictor eventually increases as well, but it remains much smaller than
$\sigma_{y}$. This observation holds for all values of $\Delta$, as shown in
Fig. 5(d). This result highlights the potential of the NNs to improve sub-
optimal estimations used to feed the training process.

###  III.5 Conclusions and outlook

We have established the potential of NNs to perform parameter estimation from
the data generated by continuous measurements of the photon-counting type. We
have carefully benchmarked their precision against the ultimate limit set by
Bayesian inference, showing that NN approaches can reach comparable
performance in a much more computationally efficient way, while also being
more robust against noise.

These results pave the way for future research that can explore more complex
quantum optical systems and measurement schemes. For example, future studies
could investigate the radiation from cavity-QED systems with non-trivial
quantum statistics [66, 86, 87, 61] or several quantum emitters, as well as
multiple simultaneous detection channels, such as those provided by arrays of
single-photon avalanche diodes for generalized Hanbury-Brown and Twiss
measurements [64, 62], streak cameras [63], or single-photon fibre bundle
cameras [65, 88].

One aspect we leave for future work is the exploration of NN methods that
incorporate mechanisms to do uncertainty quantification (UQ) on the
predictions (compute prediction intervals) [89]. For example, Bayesian Neural
Networks [90] train a set of weights distributions and the predictions from
those networks are obtained by sampling the posterior of the weights after
training. Other approaches to UQ include bootstrapping and ensemble methods
which are easily implemented and only require extra training time, maintaining
all the benefits of NN estimators that we have shown in our work.

Another aspect deserving further exploration is leveraging domain knowledge
and incorporating physical models into the learning process [2]. We note that
the _Hist-Dense_ architecture used in this study could already be considered a
physics-informed approach, taking advantage of a known feature about the
model: the absence of information contained in the particular ordering of
delays in the input data $x$, given that the source is a TLS that resets its
state after every emission [29].

Moreover, it is worth noting that protocols of parameter estimation from
photon-counting data, as the one described here, may find applications in
fields such as fluorescence quantum microscopy [91, 92, 65, 88, 93] or quantum
imaging [94, 95]. In these areas, experimental quantum signals beyond the
capabilities of classical simulations—due to the number of emitters and
detectors involved—are already within reach, calling for efficient methods of
analysis for the optimal image reconstruction [96].

##  IV Methods

###  IV.1 Photodetection time data

In practice, rather than working with the absolute times of detection,
$t_{i}$, we find it more convenient to consider $D$ to be equivalently
composed of a list of time delays $\tau_{i}=t_{i}-t_{i-1}$, with $i=1,\ldots
N$ and $t_{0}=0$. For simplicity of calculations, we consider trajectories
with different values of total evolution time $T$ but a fixed number of total
detections $N$, and assume that the measurement is immediately terminated
after the last detection, i.e., $t=t_{N}$. This means that $t-t_{N}$ will
always be zero, becoming an irrelevant time interval that we do not need to
include in our data, which will then only consist of $N$ values,
$D=[\tau_{1},\ldots\tau_{N}$]. This format provides a set of data records with
a fixed number of photodetection events or “clicks” that can be easily used as
inputs for the NNs that we use in this work.

###  IV.2 Bayesian Inference

_Bayes’ Rule.—_ We consider that, at the beginning of our estimation protocol,
our knowledge of the unknown set of parameters is characterized by a constant
prior probability distribution $P_{0}(\theta)$ over an arbitrarily large but
finite support (we define an initial support that is big enough not to affect
the final estimation, and alternative prior distributions can also be used).
Our goal is to update this knowledge upon the new evidence provided by some
measurement data $D$, i.e., to compute the posterior probability
$P(\theta|D)$. According to Bayes’ rule, this is given by

$P(\theta|D)=\frac{P(D|\theta)}{P(D)}P_{0}(\theta)=\frac{P(D|\theta)}{\int
d\theta P(D|\theta)P_{0}(\theta)}P_{0}(\theta),$ (2)

where the integral is defined over the domain of possible values of $\theta$,
i.e., the support of $P_{0}(\theta)$. If we have no prior knowledge of
$\theta$, as we assume here, $P_{0}(\theta)$ is constant over its support and
can be taken out of the integral in the denominator in Eq. (2), leaving us
with a simple relation between $P(\theta|D)$ and $P(D|\theta)$:

$P(\theta|D)=\frac{P(D|\theta)}{\int d\theta P(D|\theta)}.$ (3)

Thus, the required step to compute the posterior $P(\theta|D)$ is the
calculation of the likelihood $P(D|\theta)$, which in general can be computed
as the norm of a conditional density matrix evolving according to the list of
quantum-jump times in data $D$, as explained in detail in Ref. [78].

_Likelihood computation.—_ The calculation of the likelihood is particularly
simple for the dynamics of a TLS, since each time interval is independent of
the others, given that the system is completely reset to the ground state
$|0\rangle$ after each emission [29]. Consequently, the probability for the
set of data $D=[\tau_{1},\ldots,\tau_{N}]$ is simply given by the product of
the probabilities of each interval

$P(D|\theta)=\prod_{i=1}^{N}w(\tau_{i};\theta),$ (4)

where the waiting-time distribution $w(\tau;\theta)$ gives the probability of
two successive clicks being separated by a time interval $\tau$ and we made
explicit its dependence on the system parameters $\theta$. In our system,
governed by the master equation (1), the waiting-time distribution has the
following analytical form,

$w(\tau;\theta)=\frac{8\gamma\Omega^{2}}{R}\exp[-\gamma\tau/2]\\\
\times\sum_{\eta=-1,1}\eta\cosh\left(\frac{\tau\sqrt{\gamma^{2}-4\left(\Delta^{2}+4\Omega^{2}\right)+\eta
R}}{2\sqrt{2}}\right)$ (5)

with
$R=\sqrt{\left[\gamma^{2}+4\left(\Delta^{2}+4\Omega^{2}\right)\right]^{2}-64\gamma^{2}\Omega^{2}}$.

_Posterior computation.—_ Once the method to compute the likelihood
$P(D|\theta)$ is established, the only obstacle to using Bayes rule in Eq. (2)
is the calculation of the normalization factor in the denominator. This factor
is called the evidence or the marginal likelihood, as it marginalizes the
likelihood over all possible values of the parameters. While this integral can
be computed analytically for simple likelihood distributions, or numerically
efficiently for a low-dimensional parameter space, many sampling techniques
are typically used in Bayesian inference to bypass its direct computation.
Markov Chain Monte Carlo (MCMC) [97] is a stochastic sampling method that
allows us to compute samples from the posterior $P(\theta|D)$ as draws from
multiple Markov chains evolving in Monte Carlo time according to a stationary
distribution proportional to the posterior. Another method is Nested Sampling
(NS) [98], which is able to provide an estimate for the marginal likelihood
itself by converting the original multi-dimensional integral over the
parameter space to a one-dimensional integral over the inverse volume of the
prior support. NS methods have been shown to be effective also in the case of
multi-modal posteriors (related to parameter degeneracy) and when the
posterior distribution has heavy or light tails.

For the multi-dimensional case $n=2$, we compute the posterior probability
distributions $P(\theta|D)$ using the UltraNest package [83] in Python.
UltraNest contains several advanced algorithms to improve the efficiency and
correctness of nested sampling, including bootstrapping uncertainties on the
marginal likelihood and handling vectorized likelihood functions over
parameter sets (together with additional MPI support).

###  IV.3 Parameter estimation from classical signals

We consider that the classical version of the measured quantum signal is given
by the sample mean of the $N$ observed delays
$\bar{\tau}=\sum_{i=1}^{N}\tau_{i}/N$. Being a sum of many independent random
variables $\tau_{i}$, $\bar{\tau}$ is itself a random variable that, in the
limit of large $N$, will follow a Gaussian distribution
$\bar{\tau}\sim\mathcal{N}(\mu,\sigma)$, where $\mu=\text{E}[\tau_{i}]$ and
$\sigma^{2}=\text{Var}[\tau_{i}]/N$. Since the independent random variables
follow the waiting-time distribution given by Eq. (5), $\tau_{i}\sim
w(\tau;\theta)$, $\mu$ and $\sigma$ can be directly obtained:

$\displaystyle\mu(\theta)$
$\displaystyle=\left[\gamma\langle\hat{\sigma}^{\dagger}\hat{\sigma}\rangle_{\text{ss}}(\theta)\right]^{-1}=\frac{\gamma^{2}+4\Delta^{2}+8\Omega^{2}}{4\gamma\Omega^{2}},$
(6a) $\displaystyle\sigma(\theta)^{2}$
$\displaystyle=\frac{(\gamma^{2}+4\Delta^{2})^{2}-8(\gamma^{2}-12\Delta^{2})\Omega^{2}+64\Omega^{4}}{N16\gamma^{2}\Omega^{4}},$
(6b)

were we assume that the emission is recorded when the driven emitter has
reached a steady state and we made explicit the dependence of these quantities
on the unknown parameters $\theta$. Within a Bayesian framework, this means
that the likelihood of the classical data is

$P(\bar{\tau}|\theta)=\frac{1}{\sigma\sqrt{2\pi}}\exp\left[-\frac{1}{2}\left(\frac{\bar{\tau}-\mu}{\sigma}\right)^{2}\right],$
(7)

and the posterior is given by Eq. (3), assuming a constant prior. From here,
we obtain a single estimation by taking the mean of the posterior, i.e.,
$\hat{\theta}=\text{E}[\theta|D]=\int d\theta\theta P(\theta|D)$. We make this
choice in order to be consistent with the estimator used for Bayesian
inference with quantum data, and because it shows overall lower RMSE and bias
than the maximum-likelihood estimator, see Ref. [78].

_Recurrent Neural Network (RNN)_
---
Layer | Output shape | Activation | # Parameters
LSTM | 17 | ReLU | 1292
LSTM | 17 | ReLU | 2380
Dense | 1 | Linear | 18
Trainable params. | 3,690 |  |
Epochs | 1200 |  |
_Hist-Dense_
Layer | Output shape | Activation | # Parameters
Histogram | 700 | - | 0
Dense | 100 | ReLU | 70100
Dense | 50 | ReLU | 5050
Dense | 30 | ReLU | 1530
Dense | 1 | Linear | 31
Trainable params. | 76,711 |  |
Epochs | 1200 |  |
_Hist-Dense 2D_
Layer | Output shape | Activation | # Parameters
Histogram | 700 | - | 0
Dense | 100 | ReLU | 70100
Dense | 50 | ReLU | 5050
Dense | 30 | ReLU | 1530
Dense | 20 | ReLU | 620
Dense | 10 | ReLU | 210
Dense | 2 | Linear | 22
Trainable params. | 77,532 |  |
Epochs | 1200 |  |
Table 1: Summary of the NN architectures used for parameter estimation. All
networks are trained with the Adam optimizer, with a learning rate $0.001$,
using a batch size of 12800 and with a mean squared logarithmic error (MSLE)
loss function.

###  IV.4 Training and design of the neural networks

For the 1D estimation case, the networks are trained with $4\times 10^{6}$
trajectories in which the first 48 photodetection times are stored. From
these, 80% are taken as training data, and 20% as validation data during the
training process. Each trajectory is simulated by fixing $\Omega=\gamma$, and
choosing randomly a value of the detuning within the domain
$\Delta\in[0,5\gamma]$, which should be large enough to cover the expected
range of possible values of $\Delta$. We choose this domain to be the same as
the support of the prior in the Bayesian inference, so that both methods rely
on the same underlying assumptions about the probability distribution of
$\Delta$. All the networks are trained using the Adam optimizer. A summary of
these architectures—including number of trainable parameters and details about
the batch size and number of epochs used in training—is provided in Table 1.
All the networks are defined and trained using the Tensorflow library [99].

For the 2D parameter estimation, we train the network with a new set of
$4\times 10^{6}$ trajectories generated for random values of $\Delta$ and
$\Omega$ in the range $\Delta\in[0.,3\gamma]$ and
$\Omega\in[0.25\gamma,5\gamma]$ (lower values of $\Omega$ are not used in
order to avoid longer simulations times associated with the lower pumping
rate).

###  IV.5 Validation process

####  IV.5.1 1D case

In order to validate the performance of the NNs, we generate a set of
validation trajectories for a uniform grid of 40 values of $\Delta$ (with
$\Omega=\gamma$) in the range $\Delta\in[0,2.1\gamma]$. For each value of
$\Delta$, we generated $10^{4}$ trajectories containing 48 photodetection
events each. For each trajectory, we obtained an estimation of the parameter
for each of the methods, taking the mean of the posterior as the estimation in
the Bayesian case. We then computed the MSE from the $10^{4}$ different
predictions with each of the methods for each value of $\Delta$.

####  IV.5.2 2D case

Validation is done for a grid of $(\Delta,\Omega)$ pairs built from a grid of
40 values of $\Delta$ in the range $\Delta\in[0.,2.1\gamma]$ and a grid of 40
values of $\Omega$ in the range $\Omega\in[0.25\gamma,2.1\gamma]$, spanning a
square grid of 1600 points. For each of these points, we simulate $10^{4}$
validation trajectories with 48 clicks each.

##  V Data availability

All the data necessary to reproduce the results in this work, including
simulated quantum trajectories for training of the NNs, validation
trajectories, trained models and Bayesian inference calculations in 2D via
nested sampling are openly available on Zenodo [100].

##  VI Code availability

The codes used to generate the results in this work, including the
implementation of the proposed NN architectures as TensorFlow models, are
openly available on Zenodo and GitHub [101].

## References

  * Wiseman and Milburn [2010] H. M. Wiseman and G. J. Milburn,  _Quantum Measurement and Control_ (Cambridge Univ. Press, 2010).
  * Genois  _et al._ [2021] É. Genois, J. A. Gross, A. Di Paolo, N. J. Stevenson, G. Koolstra, A. Hashim, I. Siddiqi, and A. Blais, Quantum-Tailored Machine-Learning Characterization of a Superconducting Qubit, PRX Quantum 2, 040355 (2021).
  * Carrasco  _et al._ [2021] J. Carrasco, A. Elben, C. Kokail, B. Kraus, and P. Zoller, Theoretical and Experimental Perspectives of Quantum Verification, PRX Quantum 2, 010102 (2021).
  * Zhang  _et al._ [2017] J. Zhang, Y.-x. Liu, R.-B. Wu, K. Jacobs, and F. Nori, Quantum feedback: Theory, experiments, and applications, Physics Reports Quantum Feedback: Theory, Experiments, and Applications, 679, 1 (2017).
  * Wiseman and Milburn [1993a] H. M. Wiseman and G. J. Milburn, Interpretation of quantum jump and diffusion processes illustrated on the Bloch sphere, Physical Review A 47, 1652 (1993a).
  * Ashhab and Nori [2010] S. Ashhab and F. Nori, Control-free control: Manipulating a quantum system using only a limited set of measurements, Physical Review A 82, 062103 (2010).
  * Sayrin  _et al._ [2011] C. Sayrin, I. Dotsenko, X. Zhou, B. Peaudecerf, T. Rybarczyk, S. Gleyzes, P. Rouchon, M. Mirrahimi, H. Amini, M. Brune, J.-M. Raimond, and S. Haroche, Real-time quantum feedback prepares and stabilizes photon number states, Nature 477, 73 (2011).
  * Yonezawa  _et al._ [2012] H. Yonezawa, D. Nakane, T. A. Wheatley, K. Iwasawa, S. Takeda, H. Arao, K. Ohki, K. Tsumura, D. W. Berry, T. C. Ralph, H. M. Wiseman, E. H. Huntington, and A. Furusawa, Quantum-Enhanced Optical-Phase Tracking, Science 337, 1514 (2012).
  * Cui  _et al._ [2013] J.-M. Cui, F.-W. Sun, X.-D. Chen, Z.-J. Gong, and G.-C. Guo, Quantum Statistical Imaging of Particles without Restriction of the Diffraction Limit, Physical Review Letters 110, 153901 (2013).
  * Leverrier [2015] A. Leverrier, Composable Security Proof for Continuous-Variable Quantum Key Distribution with Coherent States, Physical Review Letters 114, 070501 (2015).
  * Liang  _et al._ [2022] K. Liang, G. Chai, Z. Cao, Y. Yuan, X. Chen, Y. Lu, and J. Peng, Bayesian Parameter Estimation for Continuous-Variable Quantum Key Distribution, Physical Review Applied 18, 054077 (2022).
  * Giovannetti  _et al._ [2006] V. Giovannetti, S. Lloyd, and L. Maccone, Quantum Metrology, Physical Review Letters 96, 010401 (2006).
  * Giovannetti  _et al._ [2011] V. Giovannetti, S. Lloyd, and L. Maccone, Advances in quantum metrology, Nature Photonics 5, 222 (2011).
  * Degen  _et al._ [2017] C. L. Degen, F. Reinhard, and P. Cappellaro, Quantum sensing, Reviews of Modern Physics 89, 035002 (2017).
  * Pirandola  _et al._ [2018] S. Pirandola, B. R. Bardhan, T. Gehring, C. Weedbrook, and S. Lloyd, Advances in photonic quantum sensing, Nature Photonics 12, 724 (2018).
  * Pezzè  _et al._ [2018] L. Pezzè, A. Smerzi, M. K. Oberthaler, R. Schmied, and P. Treutlein, Quantum metrology with nonclassical states of atomic ensembles, Reviews of Modern Physics 90, 035005 (2018).
  * Ma  _et al._ [2011] J. Ma, X. Wang, C. P. Sun, and F. Nori, Quantum spin squeezing, Physics Reports 509, 89 (2011).
  * Liu  _et al._ [2016] Z.-P. Liu, J. Zhang, Ş. K. Özdemir, B. Peng, H. Jing, X.-Y. Lü, C.-W. Li, L. Yang, F. Nori, and Y.-x. Liu, Metrology with PT-Symmetric Cavities: Enhanced Sensitivity near the PT-Phase Transition, Physical Review Letters 117, 110802 (2016).
  * Xu  _et al._ [2022] K. Xu, Y.-R. Zhang, Z.-H. Sun, H. Li, P. Song, Z. Xiang, K. Huang, H. Li, Y.-H. Shi, C.-T. Chen, X. Song, D. Zheng, F. Nori, H. Wang, and H. Fan, Metrological Characterization of Non-Gaussian Entangled States of Superconducting Qubits, Physical Review Letters 128, 150501 (2022).
  * Lee  _et al._ [2023] K.-Y. Lee, J.-D. Lin, A. Miranowicz, F. Nori, H.-Y. Ku, and Y.-N. Chen, Steering-enhanced quantum metrology using superpositions of noisy phase shifts, Physical Review Research 5, 013103 (2023).
  * Mabuchi [1996] H. Mabuchi, Dynamical identification of open quantum systems, Quantum and Semiclassical Optics: Journal of the European Optical Society Part B 8, 1103 (1996).
  * Gambetta and Wiseman [2001] J. Gambetta and H. M. Wiseman, State and dynamical parameter estimation for open quantum systems, Physical Review A 64, 14 (2001).
  * Verstraete  _et al._ [2001] F. Verstraete, A. C. Doherty, and H. Mabuchi, Sensitivity optimization in quantum parameter estimation, Physical Review A 64, 032111 (2001).
  * Chase and Geremia [2009] B. A. Chase and J. M. Geremia, Single-shot parameter estimation via continuous quantum measurement, Physical Review A 79, 022314 (2009).
  * Guţă [2011] M. Guţă, Fisher information and asymptotic normality in system identification for quantum Markov chains, Physical Review A 83, 062324 (2011).
  * Ralph  _et al._ [2011] J. F. Ralph, K. Jacobs, and C. D. Hill, Frequency tracking and parameter estimation for robust quantum state estimation, Physical Review A 84, 052119 (2011).
  * Gammelmark and Mølmer [2014] S. Gammelmark and K. Mølmer, Fisher Information and the Quantum Cramér-Rao Sensitivity Limit of Continuous Measurements, Physical Review Letters 112, 170401 (2014).
  * Gammelmark and Mølmer [2013] S. Gammelmark and K. Mølmer, Bayesian parameter inference from continuously monitored quantum systems, Physical Review A 87, 032115 (2013).
  * Kiilerich and Mølmer [2014] A. H. Kiilerich and K. Mølmer, Estimation of atomic interaction parameters by photon counting, Physical Review A 89, 052110 (2014).
  * Kiilerich and Mølmer [2015] A. H. Kiilerich and K. Mølmer, Parameter estimation by multichannel photon counting, Physical Review A 91, 12119 (2015).
  * Kiilerich and Mølmer [2016] A. H. Kiilerich and K. Mølmer, Bayesian parameter estimation by continuous homodyne detection, Physical Review A 94, 032103 (2016).
  * Cortez  _et al._ [2017] L. Cortez, A. Chantasri, L. P. García-Pintos, J. Dressel, and A. N. Jordan, Rapid estimation of drifting parameters in continuously measured quantum systems, Physical Review A 95, 012314 (2017).
  * Ralph  _et al._ [2017] J. F. Ralph, S. Maskell, and K. Jacobs, Multiparameter estimation along quantum trajectories with sequential Monte Carlo methods, Physical Review A 96, 052306 (2017).
  * Albarelli  _et al._ [2018] F. Albarelli, M. A. C. Rossi, D. Tamascelli, and M. G. Genoni, Restoring Heisenberg scaling in noisy quantum metrology by monitoring the environment, Quantum 2, 110 (2018).
  * Kominis  _et al._ [2003] I. K. Kominis, T. W. Kornack, J. C. Allred, and M. V. Romalis, A subfemtotesla multichannel atomic magnetometer, Nature 422, 596 (2003).
  * Geremia  _et al._ [2003] JM. Geremia, J. K. Stockton, A. C. Doherty, and H. Mabuchi, Quantum Kalman Filtering and the Heisenberg Limit in Atomic Magnetometry, Physical Review Letters 91, 250801 (2003).
  * Macieszczak  _et al._ [2016] K. Macieszczak, M. Guţă, I. Lesanovsky, and J. P. Garrahan, Dynamical phase transitions as a resource for quantum enhanced metrology, Physical Review A 93, 022103 (2016).
  * Ilias  _et al._ [2022] T. Ilias, D. Yang, S. F. Huelga, and M. B. Plenio, Criticality-Enhanced Quantum Sensing via Continuous Measurement, PRX Quantum 3, 010354 (2022).
  * Tsang  _et al._ [2011] M. Tsang, H. M. Wiseman, and C. M. Caves, Fundamental Quantum Limit to Waveform Estimation, Physical Review Letters 106, 090401 (2011).
  * Smith  _et al._ [2006] G. A. Smith, A. Silberfarb, I. H. Deutsch, and P. S. Jessen, Efficient Quantum-State Estimation by Continuous Weak Measurement and Dynamical Control, Physical Review Letters 97, 180403 (2006).
  * Ashhab  _et al._ [2009] S. Ashhab, J. Q. You, and F. Nori, The information about the state of a qubit gained by a weakly coupled detector, New Journal of Physics 11, 083017 (2009).
  * Gong  _et al._ [2018] B. Gong, Y. Yang, and W. Cui, Quantum parameter estimation via dispersive measurement in circuit QED, Quantum Information Processing 17, 301 (2018).
  * Wiseman and Milburn [1993b] H. M. Wiseman and G. J. Milburn, Quantum theory of field-quadrature measurements, Physical Review A 47, 642 (1993b).
  * Chen  _et al._ [2022] Y. Chen, Y. Ban, R. He, J.-M. Cui, Y.-F. Huang, C.-F. Li, G.-C. Guo, and J. Casanova, A neural network assisted 171Yb+ quantum magnetometer, npj Quantum Information 8, 1 (2022).
  * Verhelst and Murmann [2020] M. Verhelst and B. Murmann, Machine Learning at the Edge, in  _NANO-CHIPS 2030: On-Chip AI for an Efficient Data-Driven World_, The Frontiers Collection, edited by B. Murmann and B. Hoefflinger (Springer International Publishing, Cham, 2020) pp. 293–322.
  * Homrighausen  _et al._ [2023] J. Homrighausen, L. Horsthemke, J. Pogorzelski, S. Trinschek, P. Glösekötter, and M. Gregor, Edge-Machine-Learning-Assisted Robust Magnetometer Based on Randomly Oriented NV-Ensembles in Diamond, Sensors 23, 1119 (2023).
  * Liu  _et al._ [2019] W. Liu, J. Huang, Y. Li, H. Li, C. Fang, Y. Yu, and G. Zeng, Parameter estimation via weak measurement with machine learning, Journal of Physics B: Atomic, Molecular and Optical Physics 52, 045504 (2019).
  * Ban  _et al._ [2021] Y. Ban, J. Echanobe, Y. Ding, R. Puebla, and J. Casanova, Neural-network-based parameter estimation for quantum detection, Quantum Science and Technology 6, 045012 (2021).
  * Nolan  _et al._ [2021] S. P. Nolan, L. Pezzè, and A. Smerzi, Frequentist parameter estimation with supervised learning, AVS Quantum Science 3 (2021).
  * Xiao  _et al._ [2022] T. Xiao, J. Fan, and G. Zeng, Parameter estimation in quantum sensing based on deep reinforcement learning, npj Quantum Information 8, 1 (2022).
  * Ban  _et al._ [2022] Y. Ban, J. Casanova, and R. Puebla, Neural networks for Bayesian quantum many-body magnetometry 10.48550/arXiv.2212.12058 (2022).
  * Cimini  _et al._ [2019] V. Cimini, I. Gianani, N. Spagnolo, F. Leccese, F. Sciarrino, and M. Barbieri, Calibration of Quantum Sensors by Neural Networks, Physical Review Letters 123, 230502 (2019).
  * Cimini  _et al._ [2021] V. Cimini, E. Polino, M. Valeri, I. Gianani, N. Spagnolo, G. Corrielli, A. Crespi, R. Osellame, M. Barbieri, and F. Sciarrino, Calibration of Multiparameter Sensors via Machine Learning at the Single-Photon Level, Physical Review Applied 15, 044003 (2021).
  * Cimini  _et al._ [2023] V. Cimini, M. Valeri, E. Polino, S. Piacentini, F. Ceccarelli, G. Corrielli, N. Spagnolo, R. Osellame, and F. Sciarrino, Deep reinforcement learning for quantum multiparameter estimation, Advanced Photonics 5, 016005 (2023).
  * Khanahmadi and Mølmer [2021] M. Khanahmadi and K. Mølmer, Time-dependent atomic magnetometry with a recurrent neural network, Physical Review A 103, 032406 (2021).
  * Ahmed  _et al._ [2021a] S. Ahmed, C. Sánchez Muñoz, F. Nori, and A. F. Kockum, Quantum State Tomography with Conditional Generative Adversarial Networks, Physical Review Letters 127, 140502 (2021a).
  * Ahmed  _et al._ [2021b] S. Ahmed, C. Sánchez Muñoz, F. Nori, and A. F. Kockum, Classification and reconstruction of optical quantum states with deep neural networks, Physical Review Research 3, 033278 (2021b).
  * Flurin  _et al._ [2020] E. Flurin, L. S. Martin, S. Hacohen-Gourgy, and I. Siddiqi, Using a Recurrent Neural Network to Reconstruct Quantum Dynamics of a Superconducting Qubit from Physical Observations, Physical Review X 10, 011006 (2020).
  * Greplova  _et al._ [2017] E. Greplova, C. K. Andersen, and K. Mølmer, Quantum parameter estimation with a neural network 10.48550/arXiv.1711.05238 (2017).
  * Bartolo  _et al._ [2017] N. Bartolo, F. Minganti, J. Lolli, and C. Ciuti, Homodyne versus photon-counting quantum trajectories for dissipative Kerr resonators with two-photon driving, The European Physical Journal Special Topics 226, 2705 (2017).
  * Zubizarreta Casalengua  _et al._ [2020] E. Zubizarreta Casalengua, J. C. López Carreño, F. P. Laussy, and E. del Valle, Conventional and Unconventional Photon Statistics, Laser & Photonics Reviews 14, 1900279 (2020).
  * De Santis  _et al._ [2017] L. De Santis, C. Antón, B. Reznychenko, N. Somaschi, G. Coppola, J. Senellart, C. Gómez, A. Lemaître, I. Sagnes, A. G. White, L. Lanco, A. Auffèves, and P. Senellart, A solid-state single-photon filter, Nature Nanotechnology 12, 663 (2017).
  * Silva  _et al._ [2016] B. Silva, C. Sánchez Muñoz, D. Ballarini, A. González-Tudela, M. de Giorgi, G. Gigli, K. West, L. Pfeiffer, E. del Valle, D. Sanvitto, and F. P. Laussy, The colored Hanbury Brown–Twiss effect, Scientific Reports 6, 37980 (2016).
  * Rundquist  _et al._ [2014] A. Rundquist, M. Bajcsy, A. Majumdar, T. Sarmiento, K. Fischer, K. G. Lagoudakis, S. Buckley, A. Y. Piggott, and J. Vučković, Nonclassical higher-order photon correlations with a quantum dot strongly coupled to a photonic-crystal nanocavity, Physical Review A 90, 023846 (2014).
  * Israel  _et al._ [2017] Y. Israel, R. Tenne, D. Oron, and Y. Silberberg, Quantum correlation enhanced super-resolution localization microscopy enabled by a fibre bundle camera, Nature Communications 8, 1 (2017).
  * Lambert  _et al._ [2010] N. Lambert, Y.-N. Chen, and F. Nori, Unified single-photon and single-electron counting statistics: From cavity QED to electron transport, Physical Review A 82, 063840 (2010).
  * Hadfield [2009] R. H. Hadfield, Single-photon detectors for optical quantum information applications, Nature Photonics 3, 696 (2009).
  * Kimble  _et al._ [1977] H. J. Kimble, M. Dagenais, and L. Mandel, Photon Antibunching in Resonance Fluorescence, Phys. Rev. Lett. 39, 691 (1977).
  * Lodahl  _et al._ [2015] P. Lodahl, S. Mahmoodian, and S. Stobbe, Interfacing single photons and single quantum dots with photonic nanostructures, Reviews of Modern Physics 87, 347 (2015).
  * Toninelli  _et al._ [2021] C. Toninelli, I. Gerhardt, A. S. Clark, A. Reserbat-Plantey, S. Götzinger, Z. Ristanović, M. Colautti, P. Lombardi, K. D. Major, I. Deperasińska, W. H. Pernice, F. H. L. Koppens, B. Kozankiewicz, A. Gourdon, V. Sandoghdar, and M. Orrit, Single organic molecules for photonic quantum technologies, Nature Materials 20, 1615 (2021).
  * Sipahigil  _et al._ [2016] A. Sipahigil, R. E. Evans, D. D. Sukachev, M. J. Burek, J. Borregaard, M. K. Bhaskar, C. T. Nguyen, J. L. Pacheco, H. A. Atikian, C. Meuwly, R. M. Camacho, F. Jelezko, E. Bielejec, H. Park, M. Lončar, and M. D. Lukin, An integrated diamond nanophotonics platform for quantum-optical networks, Science 354, 847 (2016).
  * Breuer and Petruccione [2002] H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford Univ. Press, 2002).
  * Plenio and Knight [1998] M. B. Plenio and P. L. Knight, The quantum-jump approach to dissipative dynamics in quantum optics, Reviews of Modern Physics 70, 101 (1998).
  * Johansson  _et al._ [2012] J. R. Johansson, P. D. Nation, and F. Nori, QuTiP: An open-source Python framework for the dynamics of open quantum systems, Computer Physics Communications 183, 1760 (2012).
  * Johansson  _et al._ [2013] J. R. Johansson, P. D. Nation, and F. Nori, QuTiP 2: A Python framework for the dynamics of open quantum systems, Computer Physics Communications 184, 1234 (2013).
  * Gneiting  _et al._ [2021] C. Gneiting, A. V. Rozhkov, and F. Nori, Jump-time unraveling of Markovian open quantum systems, Physical Review A 104, 062212 (2021).
  * Gneiting  _et al._ [2022] C. Gneiting, A. Koottandavida, A. V. Rozhkov, and F. Nori, Unraveling the topology of dissipative quantum systems, Physical Review Research 4, 023036 (2022).
  * Note [1] See Supplemental Material for further details on the theory of quantum jump trajectories, the choice of different estimators from posterior distributions, and the Fisher information and the CRB.
  * Jaynes [2003] E. T. Jaynes,  _Probability theory: The logic of science_ (Cambridge University Press, Cambridge, 2003).
  * Paris [2009] M. G. Paris, Quantum estimation for quantum technology, International Journal of Quantum Information 7, 125 (2009).
  * Cramér [1999] H. Cramér,  _Mathematical Methods of Statistics_ (Princeton University Press, 1999).
  * Matson and Haji [2006] C. L. Matson and A. Haji, Biased Cramér-Rao lower bound calculations for inequality-constrained estimators, Journal of the Optical Society of America. A, Optics, Image Science, and Vision 23, 2702 (2006).
  * Buchner [2021] J. Buchner, UltraNest - a robust, general purpose Bayesian inference engine, The Journal of Open Source Software 6, 3001 (2021).
  * López Carreño  _et al._ [2022] J. C. López Carreño, E. Zubizarreta Casalengua, B. Silva, E. Del Valle, and F. P. Laussy, Loss of antibunching, Physical Review A 105, 023724 (2022).
  * Warszawski  _et al._ [2002] P. Warszawski, H. M. Wiseman, and H. Mabuchi, Quantum trajectories for realistic detection, Physical Review A 65, 023802 (2002).
  * Muñoz  _et al._ [2014] C. S. Muñoz, E. del Valle, A. G. Tudela, K. Müller, S. Lichtmannecker, M. Kaniber, C. Tejedor, J. J. Finley, and F. P. Laussy, Emitters of N-photon bundles, Nature Photonics 8, 550 (2014).
  * Hamsen  _et al._ [2017] C. Hamsen, K. N. Tolazzi, T. Wilk, and G. Rempe, Two-Photon Blockade in an Atom-Driven Cavity QED System, Physical Review Letters 118, 133604 (2017).
  * Tenne  _et al._ [2019] R. Tenne, U. Rossman, B. Rephael, Y. Israel, A. Krupinski-Ptaszek, R. Lapkiewicz, Y. Silberberg, and D. Oron, Super-resolution enhancement by quantum image scanning microscopy, Nature Photonics 13, 116 (2019).
  * Khosravi  _et al._ [2011] A. R. Khosravi, S. Nahavandi, D. Creighton, and A. F. Atiya, Comprehensive review of neural network-based prediction intervals and new advances, IEEE Transactions on Neural Networks 22, 1341 (2011).
  * Wang and Yeung [2020] H. Wang and D. Yeung, A survey on Bayesian deep learning, ACM Computing Surveys 53, 1 (2020).
  * Delaubert  _et al._ [2008] V. Delaubert, N. Treps, C. Fabre, H. A. Bachor, and P. Réfrégier, Quantum limits in image processing, EPL (Europhysics Letters) 81, 44001 (2008).
  * Speirits  _et al._ [2017] F. C. Speirits, M. Sonnleitner, and S. M. Barnett, From retrodiction to Bayesian quantum imaging, Journal of Optics 19, 044001 (2017).
  * He  _et al._ [2023] Z. He, Y. Zhang, X. Tong, L. Li, and L. V. Wang, Quantum microscopy of cells at the Heisenberg limit, Nature Communications 14, 2441 (2023).
  * Albarelli  _et al._ [2020] F. Albarelli, M. Barbieri, M. G. Genoni, and I. Gianani, A perspective on multiparameter quantum metrology: From theoretical tools to applications in quantum imaging, Physics Letters A 384, 126311 (2020).
  * Wolley  _et al._ [2022] O. Wolley, T. Gregory, S. Beer, T. Higuchi, M. Padgett, and H. K. Photonics, Quantum imaging with a photon counting camera, Scientific Reports 2022 12:1 12, 1 (2022).
  * Altmann  _et al._ [2018] Y. Altmann, S. McLaughlin, M. J. Padgett, V. K. Goyal, A. O. Hero, and D. Faccio, Quantum-inspired computational imaging, Science 361, eaat2298 (2018).
  * Speagle [2020] J. S. Speagle, A Conceptual Introduction to Markov Chain Monte Carlo Methods 10.48550/arXiv.1909.12313 (2020).
  * Buchner [2023] J. Buchner, Nested sampling methods, Statistics Surveys 17, 169 (2023).
  * Abadi  _et al._ [2016] M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean, M. Devin, S. Ghemawat, G. Irving, M. Isard, M. Kudlur, J. Levenberg, R. Monga, S. Moore, D. G. Murray, B. Steiner, P. Tucker, V. Vasudevan, P. Warden, M. Wicke, Y. Yu, and X. Zheng, TensorFlow: A system for large-scale machine learning 10.48550/arXiv.1605.08695 (2016).
  * Sánchez Muñoz and Rinaldi [2023a] C. Sánchez Muñoz and E. Rinaldi, Dataset: Parameter estimation by learning quantum correlations in continuous photon-counting data using neural networks (2023a).
  * Sánchez Muñoz and Rinaldi [2023b] C. Sánchez Muñoz and E. Rinaldi, Code: Parameter estimation by learning quantum correlations in continuous photon-counting data using neural networks (2023b).

##  VII Acknowledgments

Authors thank Yexiong Zeng, Andrey Kardashin, Clemens Gneiting, Yanming Che,
Fabrizio Minganti, Anton Frisk Kockum, Yuta Kikuchi, and Marcello Benedetti
for critical reading and feedback. E. R. was supported by Nippon Telegraph and
Telephone Corporation (NTT) Research during the early stages of this work. S.
A. and M. K. acknowledge support from the Knut and Alice Wallenberg Foundation
through the Wallenberg Centre for Quantum Technology (WACQT). F. N. is
supported in part by Nippon Telegraph and Telephone Corporation (NTT)
Research, the Japan Science and Technology Agency (JST) [via the Quantum Leap
Flagship Program (Q-LEAP) and the Moonshot R&D Grant No. JPMJMS2061], the
Asian Office of Aerospace Research and Development (AOARD) (via Grant No.
FA2386-20-1-4069), and the Foundational Questions Institute Fund (FQXi) via
Grant No. FQXiIAF19-06. C.S.M. acknowledges that the project that gave rise to
these results received the support of a fellowship from “la Caixa” Foundation
(ID 100010434) and from the European Union’s Horizon 2020 research and
innovation programme under the Marie Skłodowska-Curie Grant Agreement
No.847648, with fellowship code LCF/BQ/PI20/11760026, and financial support
from the MCINN project PID2021-126964OB-I00 (QENIGMA) and the Proyecto
Sinérgico CAM 2020 Y2020/TCS- 6545 (NanoQuCo-CM).

Supplemental Material

##  Appendix A Unravelling dissipative dynamics into quantum trajectories

To better understand the origin of this decomposition of the dynamics into
trajectories experiencing quantum jumps, we may write the master equation in
Eq. (1) in superoperator form as
$\partial_{t}\hat{\rho}=\mathcal{L}_{\theta}\hat{\rho}$, where
$\mathcal{L}_{\theta}$ is the Liouvillian superoperator, and the subscript
$\theta$ makes explicit its dependence on the unknown set of parameters we
want to estimate, e.g., $\theta=[\Delta,\Omega]$, which are encoded into the
state through the Liouvillian evolution. We can separate the Liouvillian in
two terms, $\mathcal{L}_{\theta}=\mathcal{J}+(\mathcal{L}-\mathcal{J})$, where
$\mathcal{J}$ represents action of a quantum jump onto the system [43]:

$\mathcal{J}\equiv\gamma\hat{\sigma}\hat{\rho}\hat{\sigma}^{\dagger}.$ (S1)

The unravelling into different trajectories becomes apparent by writing the
density matrix at any time $T$ as a mixture of non-normalized conditional
density matrices using a generalized Dyson expansion of the solution of the
master equation,

$\rho(T,\theta)=\sum_{N=0}^{\infty}\int_{0}^{T}dt_{N}\int_{0}^{t_{N}}dt_{N-1}\cdots\int_{0}^{t_{2}}dt_{1}\\\
\times\tilde{\rho}(t_{1},\cdots,t_{N};T;\theta),$ (S2)

each of these matrices being defined as

$\tilde{\rho}(t_{1},\cdots,t_{N};T;\theta)=\mathcal{S}_{\theta}(T-t_{N})\mathcal{J}\\\
\times\mathcal{S}_{\theta}(t_{N}-t_{N-1})\cdots\mathcal{J}\mathcal{S}_{\theta}(t_{1})\rho(0)\,,$
(S3)

where
$\mathcal{S}_{\theta}(\tau)\equiv\exp[(\mathcal{L}_{\theta}-\mathcal{J})\tau]$.
This equation represents a conditional, un-normalized density matrix
characterized by periods of evolution without a quantum jump generated by the
operator $S_{\theta}$, and quantum jumps taking place at times $t_{i}$ through
the operator $\mathcal{J}$. Equation (S2) represents the unraveling of the
evolution of $\rho$ in the time interval $[0,T]$ into all of possible
individual realizations, i.e., the quantum trajectories. A trajectory is fully
characterized by its duration, $T$, and the number of photons $N$ emitted
during that time interval as well as their corresponding times of emission. In
other words, the data array that contains all the photo-emission times,
$D=[t_{1},\ldots,t_{N}]$ is enough to tag each individual quantum trajectory.
Equation (S2) tells us that each of these realizations, tagged by the data
array $D$, has a probability of occurring given by

$P(D|\theta)=\text{Tr}\left[\tilde{\rho}(t_{1},\ldots,t_{N};T,\theta\right].$
(S4)

This equation gives us the likelihood of the data. From this, we can
reconstruct the posterior probability distribution of the unknown parameters
$\theta$ by direct application of Bayes’ rule.

The unravelling of the master equation discussed previously establishes a
connection to measurement theory and the concept of monitored realizations.
When knowledge is acquired about the system, the state of the system is
consequently updated, in accordance with the conditional density matrices
presented in Eq. (S3). It is important to note that there are multiple ways to
unravel the dynamics, as the choice of the jump operator $\mathcal{J}$ can
vary while still yielding the same master equation that we presented in (1)
[60]. Each different choice of the jump operator corresponds to a distinct
unraveling and measurement strategy. Nevertheless, by averaging the ensemble
of trajectories from each type of unraveling, we obtain the same density-
matrix evolution, which does not assume that any knowledge has been retrieved.

Notable examples of different unravellings in quantum optics include those
arising from photon-counting measurements, which exhibit abrupt jumps in
trajectories, and homodyne measurements, which result in a diffusive
evolution. These two measurement strategies access the particle-like and wave-
like nature of the electromagnetic field, respectively. Prior works utilizing
neural networks in continuously monitored systems have primarily focused on
the latter type of measurement, i.e. in continuous, diffusive signals [55,
44]. In contrast, our work focuses on unravellings with abrupt jumps,
corresponding to a “click” in the detector, producing discrete signals
consisting of these precise photodetection times.

A single trajectory describing the evolution of an individual run of the
experiment in which the system is continuously monitored cannot be predicted
deterministically, since its evolution is stochastic. This evolution preserves
the purity of the state, allowing us to work in a framework that only uses
wavefunctions. For a trajectory describing the specific case of continuous
monitoring through photon-counting measurements, the wavefunction is governed
by a nonlinear stochastic Schrödinger equation [1]

$d|\psi(t)\rangle=dN(t)\left(\frac{\hat{\sigma}}{\sqrt{\langle\hat{\sigma}^{\dagger}\hat{\sigma}\rangle(t)}}-\mathbf{I}\right)|\psi(t)\rangle\\\
+dt\left(\frac{\langle\hat{\sigma}^{\dagger}\hat{\sigma}\rangle(t)}{2}-\frac{\hat{\sigma}^{\dagger}\hat{\sigma}}{2}-i\hat{H}\right)|\psi(t)\rangle,$
(S5)

where $dN(t)^{2}=dN(t)$ and
$E[dN(t)]=dt\langle\psi(t)|\hat{\sigma}^{\dagger}\hat{\sigma}|\psi(t)\rangle$.
Notice that, following the discrete nature of photon-counting measurements,
$dN(t)$ is a discrete random variable that can be either 0 or 1, and that
$E[dN(t)]$ corresponds to the probability of detecting a photon in the time
window $(t,t+dt)$. Solving the previous equation by discretizing time leads to
the usual recipe of Monte Carlo evolution described in the main text [73].

##  Appendix B Choosing an estimator from the posterior

One of the purposes of this work is to benchmark the precision of the single-
point estimates produced by the proposed NN architectures against the
estimation provided by Bayesian inference. Making an equivalent comparison
between the two requires to build a single point estimator summarizing the
information contained in the full Bayesian posterior distribution (which
inevitably entails some loss of information). Moreover, we diminish the role
of the prior in the Bayesian inference procedure by focusing on a uniform
prior, since we do not equip the NN architectures with a mechanism similar to
a prior distribution on parameters. It is well known that, when the data is
scarce, the importance of a prior in a Bayesian protocol can be substantial.

Choosing the best estimator can be a subtle problem, since several single-
point estimators can be built from the posterior. The most relevant ones are
arguably the _maximum a posteriori_ (MAP) estimator, defined as

$\hat{\theta}_{\text{MAP}}=\operatorname*{arg\,max}_{\theta}P(\theta|D),$ (S6)

and the _mean estimator_ , defined as the mean of the posterior

$\hat{\theta}_{\text{MEAN}}=\int d\theta\,\theta P(\theta|D).$ (S7)

Notice that the MAP estimator is equivalent to the _maximum-likelihood
estimator_ if there is no previous knowledge about the parameters, i.e., if
the prior distribution $P_{0}(\theta)$ is flat.

In all the comparisons shown in this work, we have consistently selected the
mean estimator $\hat{\theta}_{\text{MEAN}}$ as the one representing the
performance of Bayesian inference, both for classical and quantum signals. We
have done this without explicitly using the subscript when referring to
$\hat{\theta}$, for the sake of notational simplicity. In this section, we
justify our choice by showing in Fig. S1 the RMSE and bias of both MAP and
mean estimators in the 1D estimation case, $\theta=[\Delta]$. These results
show that, for most values of $\Delta$, mean estimators show a lower RMSE and
bias than the MAP estimator, which justifies our choice of the mean estimator
as the best representation of Bayesian inference. Notably, we find that the
RMSE of the mean estimator coincides with the RMSE naturally reached by the NN
estimators after their training.

Figure S1: Performance of different estimators built from the Bayesian
posterior distribution $P(\theta|D)$. (a) RMSE of estimators obtained as the
mean of the posterior distribution (squares) and as the maximum a posteriori
probability (MAP) estimate (crosses). For most values of $\Delta$, the
posterior mean estimate yields a lower RMSE. (b) Bias of the same estimators
shown in (a).

##  Appendix C Fisher information

Fisher information is ubiquitous as a tool to establish the ultimate limits of
precision of parameter estimation approaches in quantum metrology [80]. The
Cramér-Rao bound (CRB) sets a limit for the variance of an _unbiased_
estimator given by [81]:

$\text{var}\hat{\theta}\geq\frac{1}{\eta F(\theta)}$ (S8)

where $\eta$ is the number of measurements and $F(\theta)$ is the Fisher
information,

$F(\theta)=\int dx\,P(x|\theta)\left(\frac{\partial\ln
P(x|\theta)}{\partial\theta}\right)^{2},$ (S9)

where we defined $x$ as the value of a measurement outcome and $P(x|\theta)$
as the probability to measure such a value given a specific value of the
parameter $\theta$.

In the context of quantum measurements, the function $P(x|\theta)$ is
determined by both the quantum state in which the measurement is performed and
by the choice of measurement, established by a particular positive operator-
valued measure (POVM), $\Pi$. For a given quantum state, the Fisher
information is thus bounded by the _quantum Fisher information_ (QFI), $H$,
which is the Fisher information corresponding to the optimal POVM:

$H=\max_{\Pi}F.$ (S10)

The notions of Fisher information and quantum Fisher information have both
been introduced and used in the context of continuous measurements in open
quantum systems that we discuss in our work [28, 27, 29]. In this case, the
quantum state for which the quantum Fisher information is defined is the
combined quantum state of system and environment [28], and the measurement
strategy that saturates the Fisher information can, in general, be a very
complicated measurement defined in this combined space. For a particular
measurement strategy, $F$ can be calculated by computing $P(x|\theta)$ from
Eq. (S4). The integral in Eq. (S9) can then be approximated by sampling
measurement outcomes, e.g., performing multiple simulations of quantum jump
trajectories [28].

Figure S2: Limits set by the Cramér-Rao bound. RMSE and bias of different
estimators compared to the lower limit set by the biased CRB and QCRB, using
the bias of each estimator.

We have shown that all the estimators used in this work show some degrees of
bias. In that case, the CRB defined by Eq. (S8) does not apply. Instead, we
need to use a generalized _biased_ CRB, which is given by [82]

$\text{var}\hat{\theta}\geq\left(1+\frac{d\,\text{bias}(\hat{\theta})}{d\theta}\right)^{2}\frac{1}{\eta
F(\theta)}.$ (S11)

It is clear from this equation that estimators with a negative derivative of
its bias will yield lower variance than that predicted by the unbiased CRB, as
it is our case. For reference, see the bias curves in Fig. 3(c).

In order to assess how far our estimators are from the lower limit set by the
CRB, we computed the Fisher information for photon counting, $F$, by
calculating the likelihood functions $P(D|\theta)$ and their derivatives over
samples of quantum trajectories, as discussed in more detail in Ref. [28]. We
also computed the quantum Fisher information, $H$, following the approach in
Ref. [27], by which $H$ is given by

$H(\theta)=4T\partial_{\theta_{1}}\partial_{\theta_{2}}\text{Re}[\lambda_{s}(\theta_{1},\theta_{2})]_{\theta_{1}=\theta_{2}=\theta},$
(S12)

where we set the time as the average duration of the trajectory
$T=N/(\gamma\langle\hat{\sigma}^{\dagger}\hat{\sigma}\rangle_{\mathrm{ss}})$,
and $\lambda_{s}(\theta_{1},\theta_{2})$ is the smallest eigenvalue of the
generalized Liouvillian $\mathcal{L}_{\theta_{1},\theta_{2}}$ which, in our
case, is defined as

$\mathcal{L}_{\theta_{1},\theta_{2}}[\hat{\rho}]\equiv-i\left(\hat{H}(\theta_{1})\hat{\rho}-\hat{\rho}\hat{H}(\theta_{2})\right)+\frac{\gamma}{2}D(\hat{\sigma})[\hat{\rho}],$
(S13)

with
$H(\theta)=\Delta\hat{\sigma}^{\dagger}\hat{\sigma}+\Omega(\hat{\sigma}+\hat{\sigma}^{\dagger})$
being the system Hamiltonian. Substituting $F$ by $H$ in Eq. (S11) we obtain a
_quantum_ Cramér-Rao bound (QCRB) for an estimator with a given bias.

Figure S2 shows the RMSE of our estimators compared with the limits
established by the biased CRB and QCRB, using the respective bias of each
estimator. We observe that the RMSE of both Bayesian and NN estimators for the
quantum data saturates the CRB for a large region of values of $\Delta$. The
lower limit set by the QCRB is smaller and not saturated by our estimators.
This is expected, since photon counting does not necessarily represent the
most optimal measurement.

##  Appendix D Deploying NN predictors on edge devices with TensorFlow Lite

The largest models trained in this work have a size of $\sim 1$MB. While this
figure is by no means large, we can further reduce the size by converting
these models into compact and efficient data structures using the TensorFlow
Lite library. This opens up the opportunity to deploy the models into
microcontrollers or other edge devices, making it possible to perform
inference on the edge in real time. Figure S3 depicts the RMSE of a NN
predictor that has been converted into a TensorFlow Lite model weighting only
93kB, showing that the performance of the model remains almost unchanged after
the compression, thanks to optimization techniques such as post-training
quantization and weights pruning.

Figure S3: Inference at the edge. RMSE of the predictions made by a trained
model that has been compressed into a TensorFlow Lite model (red crosses). The
RMSE of the Bayesian inference and of the original NN predictor are shown for
comparison.


$(ii)$ below, one can replace (5.24), (5.27), and (5.28) respectively by

$\mathbb{E}[I_{\ell}]\leq C^{*}a(R_{q_{m}})R_{q_{m}}^{d-1},\ \
\ell\in\\{1,2,\dots,i\\},$

and

$v_{m+1}^{2i}\mathbb{E}[C_{m}]\,{\bf
1}\big{\\{}|{\mathcal{Y}}\cap{\mathcal{Y}}^{\prime}|=0\big{\\}}=o\big{(}a(R_{q_{m}})R_{q_{m}}^{d-1}\big{)},\
\ \ m\to\infty,$

and

$v_{m+1}^{2i}\mathbb{E}[D_{m}]\,{\bf
1}\big{\\{}|{\mathcal{Y}}\cap{\mathcal{Y}}^{\prime}|=0\big{\\}}\leq
C^{*}a(R_{q_{m}})R_{q_{m}}^{d-1}.$

Since the rest of the argument is totally the same as that for (5.11), we now
conclude that

$\sup_{m\geq 1}\big{(}a(R_{q_{m}})R_{q_{m}}^{d-1}\big{)}^{-1}A_{m}<\infty,\ \
\text{ and }\ \ \sup_{m\geq
1}\big{(}a(R_{q_{m}})R_{q_{m}}^{d-1}\big{)}^{-1}B_{m}<\infty,$

as desired.

Finally, for (5.41) we have

$\displaystyle\big{(}a(R_{q_{m}})R_{q_{m}}^{d-1}\big{)}^{-1}\,\mathbb{E}\big{[}V_{m,M}(t)\big{]}$
$\displaystyle=\sum_{i=M+1}^{\infty}i^{k+1}\binom{v_{m+1}}{i}\big{(}a(R_{q_{m}})R_{q_{m}}^{d-1}\big{)}^{-1}\,$
$\displaystyle\quad\times\mathbb{P}\Big{(}\check{C}\big{(}\\{X_{1},\dots,X_{i}\\},t\big{)}\text{
is connected},\,\mathcal{M}(X_{1},\dots,X_{i})\geq R_{q_{m}}\Big{)}.$

Proceeding as in (5.29) and using Lemma 5.11 $(ii)$, it turns out that

$\big{(}a(R_{q_{m}})R_{q_{m}}^{d-1}\big{)}^{-1}\,\mathbb{E}\big{[}V_{m,M}(t)\big{]}\leq
C^{*}\sum_{i=M+1}^{\infty}i^{k-1}\big{(}\lambda(1+\delta)e\omega_{d}\big{)}^{i}<\infty.$

Therefore, combining Lemma 5.11 $(i)$ and the dominated convergence theorem
yields (5.41). ∎

The result below is analogous to Lemma 5.6 when the density has an
exponentially decaying tail.

######  Lemma 5.11.

$(i)$ Under the assumptions of Theorem 3.4 $(ii)$, for every $i\geq k+2$ and
$t\in[0,1]$,

$\displaystyle
v_{m+1}^{i}\big{(}a(R_{q_{m}})R_{q_{m}}^{d-1}\big{)}^{-1}\,\mathbb{P}\Big{(}\check{C}\big{(}\\{X_{1},\dots,X_{i}\\},t\big{)}\emph{
is connected},\,$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\mathcal{M}(X_{1},\dots,X_{i})\geq
R_{q_{m}}\Big{)}\to\lambda^{i}\zeta_{i}^{\prime}(t),\ \ m\to\infty,$

where $\zeta_{i}^{\prime}(t)$ is defined in (5.42).
$(ii)$ Moreover, let $\delta>0$ be a constant so small that
$\lambda(1+\delta)e\omega_{d}<1$. Then, for each $i\geq k+2$ and $t\in[0,1]$,
there exists $N\in\mathbb{N}$ such that for all $m\geq N$,

(5.55) $\displaystyle
v_{m+1}^{i}\big{(}a(R_{q_{m}})R_{q_{m}}^{d-1}\big{)}^{-1}\,\mathbb{P}\Big{(}\check{C}\big{(}\\{X_{1},\dots,X_{i}\\},t\big{)}\emph{
is connected},\,$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\mathcal{M}(X_{1},\dots,X_{i})\geq
R_{q_{m}}\Big{)}\leq C^{*}\lambda^{i}(1+\delta)^{i}i^{i-2}\omega_{d}^{i-1},$

where $C^{*}>0$ is a constant independent of $i$ and $m$.

###### Proof.

We first show (5.55). Performing the same change of variables as in (5.43) and
(5.44), we have

(5.56) $\displaystyle
v_{m+1}^{i}\big{(}a(R_{q_{m}})R_{q_{m}}^{d-1}\big{)}^{-1}\,\mathbb{P}\Big{(}\check{C}\big{(}\\{X_{1},\dots,X_{i}\\},t\big{)}\emph{
is connected},\,\mathcal{M}(X_{1},\dots,X_{i})\geq R_{q_{m}}\Big{)}$
$\displaystyle=v_{m+1}^{i}f(R_{q_{m}})^{i}\int_{S^{d-1}}J(\theta)\int_{0}^{\infty}\Big{(}1+\frac{a(R_{q_{m}})}{R_{q_{m}}}\rho\Big{)}^{d-1}$
$\displaystyle\qquad\qquad\times\int_{(\mathbb{R}^{d})^{i-1}}{\bf
1}\Big{\\{}\check{C}\big{(}\\{0,{\bf y}\\},t\big{)}\text{ is
connected}\Big{\\}}\,\frac{f\big{(}R_{q_{m}}+a(R_{q_{m}})\rho\big{)}}{f(R_{q_{m}})}$
$\displaystyle\qquad\qquad\times\prod_{\ell=1}^{i-1}\frac{f\big{(}\big{\|}\big{(}R_{q_{m}}+a(R_{q_{m}})\rho\big{)}\theta+y_{\ell}\big{\|}\big{)}}{f(R_{q_{m}})}\,{\bf
1}\Big{\\{}\big{\|}\big{(}R_{q_{m}}+a(R_{q_{m}})\rho\big{)}\theta+y_{\ell}\big{\|}\geq
R_{q_{m}}\Big{\\}}\,d{\bf y}\,d\rho\,d\theta.$

By Lemma 5.9, there exists $N_{1}\in\mathbb{N}$ so that for all $m\geq N_{1}$,
we have $v_{m+1}f(R_{q_{m}})\leq\lambda(1+\delta)$. Employing the bounds
derived in the proof of (5.39), there exists $N_{2}\in\mathbb{N}$ such that
for all $m\geq N:=N_{1}\vee N_{2}$, one can bound (5.56) by

$\displaystyle\big{(}\lambda(1+\delta)\big{)}^{i}\int_{S^{d-1}}J(\theta)\int_{0}^{\infty}2(1\vee\rho)^{d-1}\int_{(\mathbb{R}^{d})^{i-1}}{\bf
1}\Big{\\{}\check{C}\big{(}\\{0,{\bf y}\\},t\big{)}\text{ is
connected}\Big{\\}}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\times\sum_{\ell=0}^{\infty}{\bf
1}\big{\\{}0<\rho\leq\epsilon^{-1}e^{(\ell+1)\epsilon}\big{\\}}e^{-\ell}\,d{\bf
y}\,d\rho\,d\theta$
$\displaystyle=C^{*}\big{(}\lambda(1+\delta)\big{)}^{i}\int_{(\mathbb{R}^{d})^{i-1}}{\bf
1}\Big{\\{}\check{C}\big{(}\\{0,{\bf y}\\},t\big{)}\text{ is
connected}\Big{\\}}d{\bf y}$ $\displaystyle\leq
C^{*}\big{(}\lambda(1+\delta)\big{)}^{i}i^{i-2}\omega_{d}^{i-1},$

where $\epsilon\in(0,1/d)$ is determined in (5.49) and $C^{*}>0$ is a constant
independent of $i$ and $m$.

Finally, combining the convergences (5.45), (5.48), and (5.54), together with
the dominated convergence theorem, can show that (5.56) tends to
$\lambda^{i}\zeta_{i}^{\prime}(t)$ as $m\to\infty$, as required. ∎

The last two lemmas below will be used for the proof of Theorem 3.4 $(i)$. The
proof of these lemmas are considerably similar to those of Lemmas 5.7 and 5.8,
so we skip their proofs.

######  Lemma 5.12.

Under the assumptions of Theorem 3.4 $(i)$, for every $t\in[0,1]$, as
$m\to\infty$,

$\displaystyle\big{(}v_{m+1}^{k+2}a(R_{q_{m}})R_{q_{m}}^{d-1}f(R_{q_{m}})^{k+2}\big{)}^{-1}\mathbb{E}\big{[}S_{m}^{\uparrow}(t)\big{]}$
$\displaystyle\to\frac{\xi_{k}^{(k+2,1,+)}(t;0)}{(k+2)!},$
$\displaystyle\big{(}v_{m}^{k+2}a(R_{p_{m}})R_{p_{m}}^{d-1}f(R_{p_{m}})^{k+2}\big{)}^{-1}\mathbb{E}\big{[}S_{m}^{\downarrow}(t)\big{]}$
$\displaystyle\to\frac{\xi_{k}^{(k+2,1,+)}(t;0)}{(k+2)!},$

and also,

$\displaystyle\big{(}v_{m+1}^{k+3}a(R_{q_{m}})R_{q_{m}}^{d-1}f(R_{q_{m}})^{k+3}\big{)}^{-1}\mathbb{E}\big{[}W_{m}(t)\big{]}$
$\displaystyle\qquad\qquad\to\frac{1}{(k+3)!}\,\int_{S^{d-1}}\int_{0}^{\infty}\int_{(\mathbb{R}^{d})^{k+2}}{\bf
1}\Big{\\{}\check{C}\big{(}\\{0,{\bf y}\\},t\big{)}\text{ is
connected}\Big{\\}}$ $\displaystyle\qquad\qquad\qquad\qquad\qquad\times
e^{-\rho(k+3)-c^{-1}\sum_{\ell=1}^{k+2}\langle\theta,y_{\ell}\rangle}\prod_{\ell=1}^{k+2}{\bf
1}\big{\\{}\rho+c^{-1}\langle\theta,y_{\ell}\rangle\geq 0\big{\\}}\,d{\bf
y}\,d\rho\,J(\theta)\,d\theta,$

where ${\bf y}=(y_{1},\dots,y_{k+2})\in(\mathbb{R}^{d})^{k+2}$. Moreover,

$\displaystyle\sup_{m\geq
1}\big{(}v_{m+1}^{k+2}a(R_{q_{m}})R_{q_{m}}^{d-1}f(R_{q_{m}})^{k+2}\big{)}^{-1}\text{Var}\big{(}S_{m}^{\uparrow}(t)\big{)}<\infty,$
$\displaystyle\sup_{m\geq
1}\big{(}v_{m}^{k+2}a(R_{p_{m}})R_{p_{m}}^{d-1}f(R_{p_{m}})^{k+2}\big{)}^{-1}\text{Var}\big{(}S_{m}^{\downarrow}(t)\big{)}<\infty.$

######  Lemma 5.13.

Under the assumptions of Theorem 3.4 $(i)$, for every $t\in[0,1]$ we have, as
$m\to\infty$,

$\displaystyle\big{(}v_{m+1}^{k+2}a(R_{q_{m}})R_{q_{m}}^{d-1}f(R_{q_{m}})^{k+2}\big{)}^{-1}\bigg{\\{}\sum_{\begin{subarray}{c}{\mathcal{Y}}\subset{\mathcal{X}}_{v_{m+1}},\\\
|{\mathcal{Y}}|=k+2\end{subarray}}h_{t}^{(k+2,1,+)}({\mathcal{Y}})\,{\bf
1}\big{\\{}\mathcal{M}({\mathcal{Y}})\geq
R_{q_{m}}\big{\\}}-S_{m}^{\uparrow}(t)\bigg{\\}}\to 0,\ \ \text{a.s.},$
$\displaystyle\big{(}v_{m}^{k+2}a(R_{p_{m}})R_{p_{m}}^{d-1}f(R_{p_{m}})^{k+2}\big{)}^{-1}\bigg{\\{}\sum_{\begin{subarray}{c}{\mathcal{Y}}\subset{\mathcal{X}}_{v_{m}},\\\
|{\mathcal{Y}}|=k+2\end{subarray}}h_{t}^{(k+2,1,+)}({\mathcal{Y}})\,{\bf
1}\big{\\{}\mathcal{M}({\mathcal{Y}})\geq
R_{p_{m}}\big{\\}}-S_{m}^{\downarrow}(t)\bigg{\\}}\to 0,\ \ \text{a.s.},$

and further,

$\displaystyle\big{(}v_{m+1}^{k+3}a(R_{q_{m}})R_{q_{m}}^{d-1}f(R_{q_{m}})^{k+3}\big{)}^{-1}$
$\displaystyle\qquad\qquad\qquad\times\bigg{\\{}\sum_{\begin{subarray}{c}{\mathcal{Y}}\subset{\mathcal{X}}_{v_{m+1}},\\\
|{\mathcal{Y}}|=k+3\end{subarray}}{\bf
1}\big{\\{}\check{C}({\mathcal{Y}},t)\text{ is
connected},\,\mathcal{M}({\mathcal{Y}})\geq
R_{q_{m}}\big{\\}}-W_{m}(t)\bigg{\\}}\to 0,\ \ \text{a.s.}$

## References

  * [1] R. J. Adler, O. Bobrowski, and S. Weinberger.  Crackle: The homology of noise.  Discrete & Computational Geometry, 52:680–704, 2014.
  * [2] S. Bachmann and M. Reitzner.  Concentration for Poisson $U$-statistics: subgraph counts in random geometric graphs.  Stochastic Processes and their Applications, 128:3327–3352, 2018\.
  * [3] G. Balkema and P. Embrechts.  Multivariate excess distributions.  www.math.ethz.ch/ embrecht/ftp/guuspe08Jun04.pdf, 2004.
  * [4] G. Balkema and P. Embrechts.  High Risk Scenarios and Extremes: A Geometric Approach.  European Mathematical Society, 2007.
  * [5] A. Björner.  Topological methods. In Handbook of Combinatorics.  Elsevier, Amsterdam, 1995.
  * [6] L. Decreusefond, M. Schulte, and C. Thäle.  Functional Poisson approximation in Kantorovich-Rubinstein distance with applications to $U$-statistics and stochastic geometry.  The Annals of Probability, 44:2147–2197, 2016.
  * [7] L. de Haan and A. Ferreira.  Extreme Value Theory: An Introduction.  Springer, New York, 2006.
  * [8] P. Embrechts, C. Klüppelberg, and T. Mikosch.  Modelling Extremal Events: for Insurance and Finance.  Springer, New York, 1997.
  * [9] R. Ghrist.  Elementary Applied Topology.  Createspace, 2014.
  * [10] A. Goel, K. T. Duy, and K. Tsunoda.  Strong law of large numbers for Betti numbers in the thermodynamic regime.  Journal of Statistical Physics, 174, 2019.
  * [11] Y. Hiraoka, T. Shirai, and K. D. Trinh.  Limit theorems for persistence diagrams.  The Annals of Applied Probability, 28:2740–2780, 2018.
  * [12] M. Kahle and E. Meckes.  Limit theorems for Betti numbers of random simplicial complexes.  Homology, Homotopy and Applications, 15:343–374, 2013.
  * [13] P. Niyogi, S. Smale, and S. Weinberger.  A topological view of unsupervised learning from noisy data.  SIAM Journal on Computing, 40:646–663, 2011.
  * [14] T. Owada.  Functional central limit theorem for subgraph counting processes.  Electronic Journal of Probability, 22, 2017.
  * [15] T. Owada.  Limit theorems for Betti numbers of extreme sample clouds with application to persistence barcodes.  The Annals of Applied Probability, 28:2814–2854, 2018.
  * [16] T. Owada and R. J. Adler.  Limit theorems for point processes under geometric constraints (and topological crackle).  The Annals of Probability, 45:2004–2055, 2017.
  * [17] T. Owada and O. Bobrowski.  Convergence of persistence diagrams for topological crackle.  Bernoulli, 26:2275–2310, 2020.
  * [18] M. Penrose.  Random Geometric Graphs, Oxford Studies in Probability 5.  Oxford University Press, Oxford, 2003.
  * [19] S. Resnick.  Extreme Values, Regular Variation and Point Processes.  Springer-Verlag, New York, 1987.
  * [20] S. Resnick.  Heavy-Tail Phenomena: Probabilistic and Statistical Modeling.  Springer, New York, 2007.
  * [21] M. Schulte and C. Thäle.  The scaling limit of Poisson-driven order statistics with applications in geometric probability.  Stochastic Processes and their Applications, 122:4096–4120, 2012\.
  * [22] A. M. Thomas and T. Owada.  Functional limit theorems for the Euler characteristic process in the critical regime.  Advances in Applied Probability, 53:57–80, 2021.
  * [23] A. M. Thomas and T. Owada.  Functional strong law of large numbers for Euler characteristic processes of extreme sample clouds.  To appear in _Extremes_ , 2021.
  * [24] D. Yogeshwaran, E. Subag, and R. J. Adler.  Random geometric complexes in the thermodynamic regime.  Probability Theory and Related Fields, 167:107–142, 2017.