Also at ]GRID - Grupo de Investigación en Ciencia de Datos
Valencian International University - VIU, Carrer Pintor Sorolla 21, 46002
València, Spain.

# Recovering discrete delayed fractional equations from trajectories

J. Alberto Conejero<EMAIL_ADDRESS>Instituto Universitario de Matemática
Pura y Aplicada
Universitat Politècnica de València, 46022 València, Spain.  Òscar Garibo-i-
Orts  [  Instituto Universitario de Matemática Pura y Aplicada
Universitat Politècnica de València, 46022 València, Spain.  Carlos Lizama
Departamento de Matemática y Ciencia de la Computación
Universidad de Santiago de Chile, Las Sophoras 173, Estación Central,
Santiago, Chile.

(October 2022)

###### Abstract

We show how machine learning methods can unveil the fractional and delayed
nature of discrete dynamical systems. In particular, we study the case of the
fractional delayed logistic map. We show that given a trajectory, we can
detect if it has some delay effect or not, and also to characterize the
fractional component of the underlying generation model.

Fractional dynamical systems; delayed discrete fractional systems; chaotic
systems; machine learning; recurrent neural networks.

##  I Introduction

Discrete Fractional Calculus (DFC) permits a description of dynamical systems
of real-world processes with discrete nature where there exist long-term
connections between different time steps. Such connections are introduced as
discrete-memory effects in the model, and DFC has been shown a reliable
approach in order to capturing them. Recent examples can be found in the study
of chaos [45], variable-order discrete-time recurrent neural networks [16],
and fuzzy systems on discrete time domains [17]. As a matter of fact, there
have been found interesting applications in different settings such as image
encryption [3, 28, 42, 44] or earthquakes intensity modeling [22, 8], and
[48].

The origins of DFC are dated to 1956, when Kutter [25] mentioned for the first
time differences of fractional order. In 1974, Diaz and Osler [9] introduced a
discrete fractional difference operator defined as an infinite series. In
1988, Grey and Zhang [14] developed a fractional calculus for the discrete
$\nabla$ (backward) difference operator. Miller and Ross [32] defined a
fractional sum via the solution of a linear difference equation. Their
definition is the discrete analog of the Riemann-Liouville fractional
integral, which can be obtained via the solution of a linear differential
equation. In 2007, Atici and Eloe [2] introduced the Riemann-Liouville like
fractional difference by using the definition of a fractional sum of Miller
and Ross, and developed some of its properties that allow one to obtain
solutions of certain fractional difference equations. The presence of chaos in
these models was first studied for the logistic map [45]; see also [43] for
sine maps.

Later, Wu and Baleanu also studied the chaos in a delayed version of the
logistic map. They started from an expression where the forward Euler operator
$\Delta$ was equal to the nonlinear right term of the logistic, that is
$\Delta x(n):=x(n+1)-x(n)$, and then they replaced the Euler operator with the
left Caputo discrete difference operator $\Delta^{\alpha}$, obtaining

$x(n)=x(0)+\frac{\mu}{\Gamma(\nu)}\sum_{j=1}^{n}\frac{\Gamma(n-j+\nu)}{\Gamma(n-j+1)}x(j-1)(1-x(j-1)),$
(1)

where $\mu$ is a parameter and $\nu$ is a scaling factor. Such a relationship
can also be obtained by convolution, using the Cesàro numbers of order $\nu$,
$k^{\nu}(j)=\frac{\Gamma(\nu+j)}{\Gamma(\nu)\Gamma(j+1)}$ with
$j\in\mathbb{N}_{0}$, as a memory kernel in terms of the scaling factor $\nu$
[7]. It is worth mentioning that new fractional difference equations have been
introduced with a general logarithm function on time scales as a kernel
function. Such equations are known as Caputo–Hadamard fractional difference
equations [47, 39].

After their first work on discrete fractional chaos, Wu and Baleanu later
introduced a two-dimensional dynamical system including a delay term in [46],
that reads as

$\begin{array}[]{ll}x(n)&=\displaystyle
x(0)+\frac{\mu}{\Gamma(\nu)}\sum_{j=1}^{n}\frac{\Gamma(n-j+\nu)}{\Gamma(n-j+1)}x(j-1)(1-y(j-1))\\\
y(n)&=x(n-1).\end{array}$ (2)

This delayed fractional discrete dynamical system also presents chaos as can
be noticed in the Feigenbaum diagrams computed in terms of the parameter
$\mu$. As a matter of fact, in Figures 1 2, we show the Feigenbaum diagrams
computed for the initial conditions $x(0)=y(0)=0.3$ and $x(0)=y(0)=0.8$.

Figure 1: Feigenbaum plots for the dynamical system given by (1) for
$x(0)=y(0)=0.3$ and $\nu=0.01$ (top left), $\nu=0.2$ (top right), $\nu=0.6$
(bottom left), and $\nu=1$ (bottom right), For each value $\mu$, we compute
200 terms of the sequence, and we plot the last 100 values.

Figure 2: Feigenbaum plots for the dynamical system given by (1) for
$x(0)=y(0)=0.8$ and $\nu=0.01$ (top left), $\nu=0.2$ (top right), $\nu=0.6$
(bottom left), and $\nu=1$ (bottom right), For each value $\mu$, we compute
200 terms of the sequence, and we plot the last 100 values.

The emergence of machine learning as a research tool has been extended to
almost all research fields. It has shown its potential in Physics being used
for studying anomalous diffusion noisy trajectories, which are characterized
as the ones whose variance of the mean square displacement grows with respect
to the time in terms of $t^{\alpha}$, with $\alpha$ being the fractional
diffusion exponent [34]. Examples of models generating those types of
trajectories are Annealed Transient Time Motion (ATTM) [31], Continuous Time
Random Walk (CTRW) [38], Fractional Brownian Motion (FBM) [30, 19], Lévy Walks
(LW) [33, 21], and Scaled Brownian Motion (SBM) [26].

Several models have been developed in the frame of the Andi Challenge for
inferring the fractional exponent $\alpha$ and for determining the generating
model of 1D, 2D, and 3D noisy trajectories, see for instance [1, 11, 12]. We
refer the reader for a full comparison between these models to [33].

Machine learning and artificial intelligence methods have been also
successfully incorporated in the study of formal mathematical problems, as is
the case of finding the solution of nonlinear models [36, 37] or the recent
success of discovering new multiplication algorithms [10].

Fractional models are suitable for representing dynamical systems with a
memory effect, either in space or time, due to the non-local character of
fractional operators, as it can be shown in the following examples [4, 6, 18,
40]. There is not a universal definition of fractional derivative [41] and in
some cases, it is not straightforward to know which definition would be better
to choose [35].

Recently, we have studied how machine learning models can be used in
connection with fractional models in order to infer the $\mu$ parameter and
the scaling factor $\nu$ of the fractional version of the logistic equation
[5]. In this work, we study up to which point we can predict the $\mu$ and
$\nu$ parameters from short trajectories generated by the discrete fractional
delayed dynamical system given by (2). We also analyze if we can discriminate
if trajectories contained a delayed component by comparing trajectories
generated by the aforementioned systems described in (2) and (1).

For this purpose, we have used a model that combines convolutional and
recurrent neural networks, that presented successful results for inferring the
exponent $\alpha$ of anomalous diffusion trajectories [11] and for predicting
the number of new COVID19 cases [29]. We will show in which cases we can
clearly determine from a given trajectory (1) which were the parameters used
for generating them, and (2) if we can affirm that it contains a delay
component or not.

In Section II, we briefly outline the architecture of the machine learning
model used. We also show how we construct the training, validation, and test
data sets. We introduce the results in Section III, and we draw some
conclusions in Section IV.

##  II Methodology

Figure 3: Machine learning architecture used for inferring the generating
$\mu$ and $\nu$ parameters of the delayed logistic model trajectories.

We start by describing our network architecture we depict in Figure 3 and that
consists of 3 parts:

  1. 1.

First, a trajectory is processed by two convolutional layers. They are
intended for feature extraction since they retain spatial structure inherent
to trajectories by setting a patch, which behaves as a sliding window, that
goes along the complete trajectory, connecting each patch of size kernel size,
to a single neuron. By doing so, different neurons specialize in different
regions in the trajectory [23, 24].

The first convolutional layer consists of 32 filters, and the second one of
64, with a sliding window (kernel) of size 5, following the doubling rule [23,
24]. The different number of filters returns features at different scales
while retaining the spatial structure.

  2. 2.

Secondly, the features feed 3 stacked bidirectional Long Short Term Memory
(LSTM) of 64 inputs and 32 units each one [15, 27] using Adam optimizer [20].
This type of recurrent neural network is able to learn sequential dependencies
from the extracted features. For this purpose, LSTM layers include gated cells
to avoid short-term memory or gradient vanishing, allowing the layer to learn
from the complete trajectory, even from positions at the beginning of it. A
cell state channel is selectively updated to allow getting rid of irrelevant
information while retaining the important pieces.

In other words, they are able to retain memory about early positions in the
trajectory, thus considering the complete trajectory while adjusting to the
target ($\mu$ and $\nu$). We also include a dropout layer after each of these
LSTM of the $10\%$ neurons to avoid over-fitting.

  3. 3.

Lastly, the output of the LSTM blocks is passed to two dense layers of 20 and
1 unit. This last choice depends if we want to predict a single parameter or
both of them at the same time. The activation function of this last dense
layer depends on the problem to be addressed: It is a linear function for
regression of the $\mu$ and $\nu$ parameters and a sigmoid function for
determining whether a trajectory comes from a delayed model or not, that is,
for the classification problem.

In order to train this model, we built two separate data sets for training and
validating while training, the previous model. Once the model was trained, we
also constructed an additional data set for testing. These data sets are
generated according to the following indications:

  * •

$\mu\in[0.0,2.0]$ with increments of $0.001$.

  * •

$\nu\in[0.01,1]$ with increments of $0.01$.

  * •

trajectory length $N$, with $N\in[10,50]$ randomly selected.

  * •

$u_{0}\in[0,1]$ randomly chosen with a resolution of $10^{-2}$.

We point out that in order to capture the chaotic dynamics that appear in some
regions of the $\mu$ parameter (see Figures 1 and 2), the resolution of $\mu$
is increased by 10 times respect to the one of $\nu$. We summarize the data
generation procedure in the algorithm II.

Data sets creation algorithm

1:for <every $\mu$> do

2: for <every $\nu$> do

3: <i = 0>

4: while i $<$ 5 do

5: <Select random index>

6: <Select random $x(0)$>

7: <Select random trajectory length between 10 and 50

8: <Generate trajectory>

9: if trajectory length $>$ 9 then

10: if index $\leq$ 0.20 then

11: <Save test trajectory>

12: else

13: if index $\leq$ 0.35 then

14: <Save validation trajectory>

15: else

16: <Save train trajectory>

17: end if

18: end if

19: else

20: <pass>

21: end if

22: <i += 1>

23: end while

24: end for

25:end for

As described in II, for every pair of ($\mu$, $\nu$) values, we generate 5
trajectories or random length (between 10 and 50), each one with a random
initial condition chosen in each one of these 5 ranges [0.0, 0.2], (0.2, 0.4],
(0.4, 0.6], (0.6, 0.8] and (0.8, 1.0]. We also have set some restrictions on
the admissible values in a trajectory: no value in the trajectory can be
smaller than $-1$ nor greater than $3$. If one of these conditions is
satisfied, we remove the last point added to the trajectory and save the
trajectory provided that it has a length greater than 9. Once a trajectory is
generated, it is randomly assigned to the train, validation, or test data
sets, resulting in a train-validation-test split with 1,110,266 trajectories
in the train data set (65%), 255,667 trajectories in the validation data set
(65%), and 341,374 in the test data set (20%).

Since the first layer in our neural network architecture is a convolutional
layer, all trajectories have to be padded (if needed) with 0s to reach the
maximum length (50). In other words, 0s are inserted to the left of the
trajectory to make all of them have a length equal to 50, see [13, Ch. 5 & Ch.
9]. We also set a patience value equal to 20 to allow the training to end
before the selected maximum number of epochs (200) if the training process
does not improve the validation mean average error (MAE) after 20 consecutive
epochs, thus saving training time. Finally, and for the sake of completeness,
we have used a computer with 16 cores configured with 128 GB RAM and Nvidia
RTX 3090 GPU with 22 GB RAM, running Ubuntu 20.10. The complete training
process took less than 3 hours, running up to 24 epochs.

##  III Results

In Figure 4, we depict the MAE distribution for $\mu$ and $\nu$, comparing the
ground truth (x-axis) and the predicted values (y-axis). The more values we
have in the diagonal, the better the models can infer these values. The
inference of the $\mu$ and $\nu$ parameters is very consistent along the
entire range of possible values of each parameter. However, the MAE results
for the inference of $\nu$ are slightly worse than the results for $\mu$. We
also notice that the trajectory length affects the accuracy of the
predictions. Since we are dealing with time series and recurrent neural
networks, it is expected that long trajectories will return lower MAE values.

Figure 4: Truth vs predicted values of $\mu$ and $\nu$ in the validation data
set (top) and for trajectories of length between 40 and 50 (bottom).

In Table 1, we show the MAE results for different lengths bins. MAEs both for
$\mu$ and $\nu$ are in the same order of magnitude ($10^{-2}$). This is the
order taken for the $\nu$ discretization. This also supports the use of a
thinner discretization for $\mu$ ($10^{-3}$) respect to the one used for
$\nu$.

MAEs for short trajectories are in fact considerably larger than for the rest
of lengthsfor instance, more than 20000 thousand for length equal to 10 but
around 6000 for lengths between 40 and 50. As it was also noticed in [5] the
MAE slightly increases for trajectories with lengths between $40$ and $50$. It
is worth to mention, that as the length increases, the number of trajectories
in each bin decreases due to the bounds set for stopping the generation of
trajectories (no term can be between smaller than -1 nor greater than 3).

Length | MAE ($\mu$) | MAE ($\nu$)
---|---|---
10-19 | 0.0284 | 0.0318
20-29 | 0.0209 | 0.0229
30-39 | 0.0188 | 0.0217
40-50 | 0.0216 | 0.0227
All | 0.0233 | 0.0259
Table 1:  $\mu$ and $\nu$ MAE’s in the test data set.

In order to get deeper insights into the predictions obtained for the
parameter $\nu$, we explore in more detail the ones that result in high values
of MAEs. In Figure 5, we represent the frequency of trajectories with MAEs
higher than 0.05 (left) in terms of the initial conditions $x(0)$, trajectory
length, and the $\mu$ and $\nu$ parameters. We appreciate than in most cases,
the main causes are short trajectory lengths, values of $\mu$ near to $0$ and
an initial condition $x(0)$ close to 0.0 or to 1.0. On Figure 5 (right), we
show the same histograms but just for trajectories of length greater than 15.
Here, we see that despite the trajectory length disappears, we still
appreciate the same conclusions for the values of $x(0)$ and $\mu$. Besides,
we also notice that whereas trajectories with values of $\mu$ between 0.4 and
1.1 return low MAEs for $\nu$, values of $\mu$ close to 0 result in high MAE
values for $\nu$.

Figure 5: Histograms of $x(0)$, trajectory lengths, $\mu$, and $\nu$ values
for trajectories with $\nu$ MAE higher than 0.05 (left) and for trajectories
with $\nu$ MAE higher than 0.05 and length greater than 15 (right).

We also confront the values of the initial conditions with the values of $\mu$
and $\nu$. In Figure 6, we represent the average MAE for pairs of initial
conditions and values of the $\mu$ (left) and $\nu$ (right) parameters. The
dark horizontal lines represent initial conditions that provide systematically
higher values of MAE. On the one hand, for predictions of $\mu$, it seems that
the initial condition is much more relevant for values of $\mu$ smaller than
1.66 than for greater ones. On the other hand, for the predictions of $\nu$,
we can find initial conditions returning high values of MAE, mainly for values
of $\nu$ smaller than 0.4.

Figure 6: $\mu$ (left) and $\nu$ (right) values vs $x(0)$ by $\mu$ and $\nu$
MAE respectively.

We also study the MAE distribution for every $\mu$ and $\nu$ truth value. In
Figure 7, we show the $Q_{1}$, $Q_{2}$, and $Q_{3}$ quartiles of these error
distributions for each value of $\mu$ (left) and $\nu$ (right). On the one
hand, it can be clearly seen that MAE error has greater variance for $\mu$
than for $\nu$, in particular for values of $\mu$ greater than 1.25, which is
consistent with what we observed in Figure 6. On the other hand, MAEs are
considerably higher for $\nu$ rather than for $\mu$, as we have previously
deduced from the Feigenbaum diagrams in Figures 1 and 1, and from the
comparison of predicted and truth values of $\mu$ and $\nu$ in Figure 4.

Figure 7: Quartiles Q1 (red), Q2 (blue), and Q3 (green) of the MAE
distribution on the evaluation data for $\mu$ (left) and $\nu$ (right).

Due to the fractional nature of equation (2), the model has a memory component
that is strongly dependent on the initial condition $x(0)$. In Figure 8, we
show boxplots of the MAE distribution for each initial condition $x(0)$. The
initial condition is more influential in predicting $\mu$ rather than in
predicting $\nu$, since boxes (green) and whiskers (grey) are wider if we
compare with the corresponding ones for $\nu$. Again, we see that the results
are worse for initial conditions close to $0$ and to $1$, due in part to the
form of the nonlinear logistic terms in (2).

Figure 8: Box and whiskers plots for the MAE distribution obtained from
predictions of $\mu$ (left) and $\nu$ (right) on the evaluation data set in
terms of the initial condition $x(0)$. The boxes are painted in green and the
whiskers on grey.

Finally, we wonder if machine learning methods are able to determine which
trajectories are generated from (1) and which ones are given by (2). In other
words, we want to see if they can determine if a delayed is incorporated to
the model or not. To do so, we have used the data set described in this work
jointly with the data set used in [5]. This second data set was also generated
following the same procedure described in Algorithm II but with the following
specifications.

  * •

$\mu\in[2,3.2]$ discretized with a step size of $10^{-3}$,

  * •

$\nu\in[0.01,1]$ discretized with a step size of $10^{-2}$,

  * •

trajectory lengths between 10 and 50,

  * •

and initial conditions $x(0)\in[0,1]$ randomly chosen with a resolution of
$10^{-2}$.

Since the range of values for $\mu$ is bigger in the present work than in the
other, we have randomly sampled the train, validation and test data sets in
order to obtain 3 perfectly balanced pairs of sets, consisting of 618,199
trajectories for each training data set, 142,800 for each validation data set,
190,334 for each the testing data set. For this classification task we have
used the architecture shown in Figure 3 but changing the activation function
in the last dense layer to a sigmoid function. The training finishes after 43
epochs and 3 hours and 30 minutes because the 20 epochs of patience are
reached. The accuracy in the validation data set is of a 99.38%, which is very
close to the area under the curve (AUC) value obtained for the test data set
of 99.375%, see Figure 9 where we show the receiver operating characteristic
(ROC) curve and the AUC value.

Figure 9: Receiver Operating Characteristic (ROC) curve on the test data set
for the classification problem of fractional logistic trajectories with and
without delay.

##  IV Conclusions

In this work, we have considered a delayed version of the fractional logistic
equation introduced by Wu and Baleanu in [46]. We have seen that a combination
of convolutional and recurrent neural networks succeeds in training a model
that given an new trajectory produced by this model is able to infer the $\mu$
and $\nu$ parameters of the model. We have seen that the MAE of the
predictions falls almost in the order of magnitude of the discretization used.
It would be interesting to check up to which point the results can be improved
using other deep learning models, training with larger data sets, and with
data sets of longer trajectories.

In our results, we have seen the importance of the initial condition, as the
main contributor to the memory effect of the generating model. We have seen,
that the initial conditions near $0$ and $1$ tend to accumulate predictions
with higher errors in comparison with other values. However, we have detected
some bins of initial conditions that produce higher errors than other
neighbour bins, and this happens almost independently of the values of the
parameters $\mu$ and $\nu$, see Figure 6. We have also seen that the predictor
of the parameters provides pretty accurate predictions for trajectory lengths
greater than 15 and, in general, the accuracy of the predictions increases as
long as the trajectory length increases.

Finally, we have seen that these methods almost perfectly detect if a
trajectory comes from a particular model with or without delay. It would be
interesting to study whether this approach is able to determine the effect of
the delayed term and of the fractional component in other fractional dynamical
systems.

## Acknowledgements

JAC acknowledges funding from grant PID2021-124618NB-C21 funded by MCIN/AEI/
10.13039/501100011033 and by “ERDF A way of making Europe”, by the “European
Union”.

## References

  * [1] A. Argun, G. Volpe, and S. Bo.  Classification, inference and segmentation of anomalous diffusion with recurrent neural networks.  J. Phys. A: Math. Theor., 54(29):294003, 2021.
  * [2] F. M. Atici and P. W. Eloe.  A transform method in discrete fractional calculus.  Int. J. Difference Equ., 2:165 – 176, 2007.
  * [3] Y.-R. Bai, D. Baleanu, and G.-C. Wu.  A novel shuffling technique based on fractional chaotic maps.  Optik, 168:553–562, 2018.
  * [4] D. Baleanu, Z. Güvenç, Burhanettin, J. A. Tenreiro Machado, et al.  New trends in nanotechnology and fractional calculus applications, volume 10.  Springer, 2010.
  * [5] J Alberto Conejero, Òscar Garibo-i Orts, and Carlos Lizama.  Inferring the fractional nature of wu baleanu trajectories.  Nonlinear Dynamics, pages 1–11, 2023.
  * [6] J.A. Conejero, J. Franceschi, and E. Picó-Marco.  Fractional vs. ordinary control systems: What does the fractional derivative provide?  Mathematics, 10(15):2719, 2022.
  * [7] J.A. Conejero, C. Lizama, A. Mira-Iglesias, and C. Rodero.  Visibility graphs of fractional Wu–Baleanu time series.  J. Differ. Equ. Appl., 25(9-10):1321–1331, 2019.
  * [8] L. Cristofaro, R. Garra, E. Scalas, and I. Spassiani.  A fractional approach to study the pure-temporal Epidemic Type Aftershock Sequence (ETAS) process for earthquakes modeling.  arXiv preprint arXiv:2301.02586, 2022.
  * [9] J. B. Diaz and T. J. Osler.  Differences of fractional order.  Math. Comp., 28:185–202, 1974.
  * [10] A. Fawzi, M. Balog, A. Huang, T. Hubert, B. Romera-Paredes, M. Barekatain, A. Novikov, F.J.R. Ruiz, J. Schrittwieser, G. Swirszcz, et al.  Discovering faster matrix multiplication algorithms with reinforcement learning.  Nature, 610(7930):47–53, 2022.
  * [11] Ò. Garibo-i Orts, A. Baeza-Bosca, M.A. Garcia-March, and J.A. Conejero.  Efficient recurrent neural network methods for anomalously diffusing single particle short and noisy trajectories.  J. Phys. A: Math. Theor., 54(50):504002, 2021.
  * [12] A. Gentili and G. Volpe.  Characterization of anomalous diffusion classical statistics powered by deep learning (CONDOR).  J. Phys. A: Math. Theor., 54(31):314003, 2021.
  * [13] I. Goodfellow, Y. Bengio, and A. Courville.  Deep Learning.  MIT Press, 2016.
  * [14] H.L. Gray and N.F. Zhang.  On a new definition of the fractional difference.  Math. Comp., 50:513–529, 1988.
  * [15] S. Hochreiter and J. Schmidhuber.  Long short-term memory.  Neural Comput., 9(8):1735–1780, 1997.
  * [16] L.-L. Huang, J.H. Park, G.-C. Wu, and Z.-W. Mo.  Variable-order fractional discrete-time recurrent neural networks.  J. Comput. Appl. Math., 370:112633, 2020.
  * [17] L.-L. Huang, G.-C. Wu, D. Baleanu, and H.-Y. Wang.  Discrete fractional calculus for interval–valued systems.  Fuzzy Sets Syst., 404:141–158, 2021.
  * [18] E. İlhan and İO. Kıymaz.  A generalization of truncated M-fractional derivative and applications to fractional differential equations.  Applied Math. Nonlinear Sci., 5(1):171–188, 2020.
  * [19] J.H. Jeon and R. Metzler.  Fractional Brownian motion and motion governed by the fractional Langevin equation in confined geometries.  Phys. Rev. E, 81:021103, 2010.
  * [20] D.P. Kingma and J. Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [21] J. Klafter and G. Zumofen.  Lévy statistics in a Hamiltonian system.  Phys. Rev. E, 49:4873–4877, 1994.
  * [22] H. Kong, G. Yang, and C. Luo.  Modeling aftershocks by fractional calculus: Exact discretization versus approximate discretization.  Fractals, 29(08):2140038, 2021.
  * [23] A. Krizhevsky, V. Nair, and G. Hinton.  CIFAR-10 (Canadian Institute for Advanced Research).  URL http://www. cs. toronto. edu/kriz/cifar. html, 5(4):1, 2010\.
  * [24] A. Krizhevsky, I. Sutskever, and G.E. Hinton.  ImageNet classification with deep convolutional neural networks.  In Proceedings of the 25th International Conference on Neural Information Processing Systems - Volume 1, page 1097–1105, Red Hook, NY, USA, 2012. Curran Associates Inc.
  * [25] B. Kuttner.  On differences of fractional order.  Proc. London Math. Soc., 3:453–466, 1957.
  * [26] S.C. Lim and S.V. Muniandy.  Self-similar Gaussian processes for modeling anomalous diffusion.  Phys. Rev. E, 66(2 Pt 1):021114, 2002.
  * [27] Z.C. Lipton.  A critical review of recurrent neural networks for sequence learning.  arXiv, abs/1506.00019, 2015.
  * [28] F. Liu, M. Meerschaert, R. McGough, P. Zhuang, and Q. Liu.  Numerical methods for solving the multi-term time-fractional wave-diffusion equation.  Fract. Calc. Appl. Anal., 16(1):9–25, 2013.
  * [29] M.A. Lozano, Ò. Garibo-i Orts, E. Piñol, M. Rebollo, K. Polotskaya, M.A. Garcia-March, J.A. Conejero, F. Escolano, and N. Oliver.  Open data science to fight COVID-19: Winning the 500k Xprize pandemic response challenge.  In Joint European Conference on Machine Learning and Knowledge Discovery in Databases, pages 384–399, 2021.
  * [30] B.B. Mandelbrot and J.W. Van Ness.  Fractional Brownian motions, fractional noises and applications.  SIAM Rev., 10(4):422–437, 1968.
  * [31] P. Massignan, C. Manzo, J.A. Torreño-Pina, M.F. Garcia-Parajo, M. Lewenstein, and G.J. Lapeyre.  Nonergodic subdiffusion from Brownian motion in an inhomogeneous medium.  Phys. Rev. Lett., 112(15):150603, 2014.
  * [32] K. S. Miller and B. Ross.  Fractional difference calculus.  Univalent functions, fractional calculus, and their applications,(Köriyama, 1988), Horwood, Chichester., pages 139–152, 1989.
  * [33] G. Muñoz-Gil, G. Volpe, M.A. Garcia-March, E. Aghion, A. Argun, C.B. Hong, T. Bland, S. Bo, J.A. Conejero, N. Firbas, and Òscar et al. Garibo-i Orts.  Objective comparison of methods to decode anomalous diffusion.  Nature Commun., 12(1), 2021.
  * [34] G. Muñoz-Gil, G. Volpe, M.A. García-March, R. Metzler, M. Lewenstein, and C. Manzo.  The Anomalous Diffusion challenge: objective comparison of methods to decode anomalous diffusion.  In G.Volpe, J.B. Pereira, D.Brunner, and A. Ozcan, editors, Emerging Topics in Artificial Intelligence (ETAI) 2021, volume 11804, page 1180416\. Int. Soc. Opt. Photonics, SPIE, 2021.
  * [35] M. Ortigueira and J.A.T. Machado.  Which derivative?  Fract. Fractional, 1(1):3, 2017.
  * [36] M.A.Z. Raja, M. Umar, Z. Sabir, J.A. Khan, and D. Baleanu.  A new stochastic computing paradigm for the dynamics of nonlinear singular heat conduction model of the human head.  Eur. Phys. J. Plus, 133(9):1–21, 2018.
  * [37] Z. Sabir, M.A.Z. Raja, D. Baleanu, K. Cengiz, and M. Shoaib.  Design of Gudermannian Neuroswarming to solve the singular Emden–Fowler nonlinear model numerically.  Nonlinear Dynam., 106(4):3199–3214, 2021.
  * [38] H. Scher and E.W. Montroll.  Anomalous transit-time dispersion in amorphous solids.  Phys. Rev. B, 12:2455–2477, 1975.
  * [39] T.T. Song, G.-C. Wu, J.-L. Wei, and et al.  Hadamard fractional calculus on time scales.  Fractals, 30(07):1–14, 2022.
  * [40] H. G. Sun, Y. Zhang, D. Baleanu, W. Chen, and Y. Q. Chen.  A new collection of real world applications of fractional calculus in science and engineering.  Commun. Nonlinear Sci. Numer. Simulat., 64:213–231, 2018.
  * [41] D. Valério, M.D. Ortigueira, and A.M. Lopes.  How many fractional derivatives are there?  Mathematics, 10(5):737, 2022.
  * [42] G.-C. Wu, D. Baleanu, and Z.-X. Lin.  Image encryption technique based on fractional chaotic time series.  J. Vib. Control, 22(8):2092–2099, 2016.
  * [43] G.-C. Wu, D. Baleanu, and S.-D. Zeng.  Discrete chaos in fractional sine and standard maps.  Phys. Lett. A, 378(5-6):484–487, 2014.
  * [44] G.-C. Wu, Z.-G. Deng, D. Baleanu, and D.-Q. Zeng.  New variable-order fractional chaotic systems for fast image encryption.  Chaos, 29(8):083103, 2019.
  * [45] G.C. Wu and D. Baleanu.  Discrete fractional logistic map and its chaos.  Nonlinear Dynam., 75(1-2):283–287, 2014.
  * [46] G.C. Wu and D. Baleanu.  Discrete chaos in fractional delayed logistic maps.  Nonlinear Dynam., 80(4):1697–1703, 2015.
  * [47] G.C. Wu, T.T. Song, and S.Q. Wang.  Caputo–hadamard fractional differential equations on time scales: Numerical scheme, asymptotic stability, and chaos.  Chaos, 32(9):093143, 2022.
  * [48] M. Ye, J. Li, X. Han, and H. Jiang.  Global analysis of a fractional-order reaction-diffusion rumor-spreading model in multi-lingual environment.  Eur. Phys. J. Plus, 138(1):102, 2023.


# Neural networks for Bayesian quantum many-body magnetometry

Yue Ban1, Jorge Casanova2,3,4 and Ricardo Puebla5 1 TECNALIA, Basque Research
and Technology Alliance (BRTA), 48160 Derio, Spain  2 Department of Physical
Chemistry, University of the Basque Country UPV/EHU, Apartado 644, 48080
Bilbao, Spain  3 EHU Quantum Center, University of the Basque Country UPV/EHU,
Leioa, Spain  4 IKERBASQUE, Basque Foundation for Science, Plaza Euskadi 5,
48009 Bilbao, Spain  5 Departamento de Física, Universidad Carlos III de
Madrid, Avda. de la Universidad 30, 28911 Leganés, Spain

###### Abstract

Entangled quantum many-body systems can be used as sensors that enable the
estimation of parameters with a precision larger than that achievable with
ensembles of individual quantum detectors. Typically, the parameter estimation
strategy requires the microscopic modelling of the quantum many-body system,
as well as a an accurate description of its dynamics. This entails a
complexity that can hinder the applicability of Bayesian inference techniques.
In this work we show how to circumvent these issues by using neural networks
that faithfully reproduce the dynamics of quantum many-body sensors, thus
allowing for an efficient Bayesian analysis. We exemplify with an XXZ model
driven by magnetic fields, and show that our method is capable to yield an
estimation of field parameters beyond the standard quantum limit scaling. Our
work paves the way for the practical use of quantum many-body systems as
black-box sensors exploiting quantum resources to improve precision
estimation.

##  1 Introduction

The key objective of quantum metrology [1, 2, 3] consists in enhancing the
precision of measured physical parameters exploiting quantum resources. In
this context, the quantum Fisher information plays a central role as it sets a
bound on the maximum achievable precision of the target parameter [4, 5]. In
particular, it has been shown that a quantum entangled many-body system acting
as a quantum sensor offers advantages with respect to individual quantum
registers: An entangled state of $N$ quantum particles can reach the so-called
Heisenberg limit where the precision scales as $1/N$ [3, 1, 2, 6] rather than
$1/\sqrt{N}$ that corresponds to the standard quantum limit. Note, however,
that decoherence effects can push the precision scaling back to the standard
quantum limit [7]. Remarkably, for specific scenarios and interactions even
better precision scalings have been identified [8, 9, 10]. Hence, devising
optimal quantum metrological protocols to boost parameter precision is
inherently linked with quantum many-body physics [3, 11, 12, 13, 14]. This is
particularly relevant for critical systems, i.e. systems featuring a
continuous phase transition, where the estimation benefits from the high
sensitivity of the system close to its critical point [15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25].

In another vein, machine learning (ML) [26] has become an effective and
increasingly popular tool to address complex problems in distinct fields of
modern quantum technologies [27, 28, 29]. Assisted by ML, several protocols
have been put forward for quantum phase estimation [30, 31, 32], parameter
estimation [33, 34, 35, 36, 37, 38], calibration of quantum sensors [39, 40],
quantum optimal control [41, 42, 43, 44, 45, 46, 47], and even to design
quantum experiments [48, 49, 50, 51], among many others. For example, in the
context of quantum sensing, neural networks (NNs) allow for the use of quantum
sensors even in regimes where the output response becomes complex and even in
the presence of large shot noise, thus extending its working regime [52, 53].
Employing a NN that replaces the computation of a physical observable offers a
drastic reduction in the amount of sample statistic [54].

On the other hand, Bayesian inference schemes [55, 56] enable the accurate
estimation of external parameters from a complex sensor response in quantum
systems [57]. The combination of ML techniques for Bayesian analysis have led
to powerful inference protocols [36]. However, in the context of many-body
system, Bayesian inference poses a challenge as only very specific many-body
systems admit a closed form solution. This is so because the computation of
the dynamics of a quantum many-body system can be a daunting task even for few
spins, and totally unfeasible for larger systems.

In this paper, we design a quantum many-body magnetometer assisted by NNs that
enable Bayesian analysis. We illustrate our protocol for a XXZ model that acts
as a magnetometer, and show that target magnetic fields can be detected
efficiently with high precision (i.e. beyond the standard quantum limit
scaling) via Bayesian inference. Playing the role in constructing the
likelihood, NNs reproduce the microscopic dynamics with high accuracy, even
with the consideration of shot-noise measurements and decoherence.
Importantly, the training of the NNs could be performed solely based on
experimental data, thus opening the door to employ experimental quantum many-
body systems as sensors that are beyond the classical simulation capabilities.
Our results demonstrate that the super-Heisenberg scaling reported in [18] is
recovered when replacing the exact dynamics by the output provided by a
suitably trained NN. In Sec. 2, we introduce our protocol of Bayesian
inference which utilizes NNs to do parameter estimation. In Sec. 3, we
elaborate the approximated dynamics with the application of NNs for magnetic
fields with one or two vector components. In Sec. 4, we analyze the precision
scaling law with the system size. Finally, we present the conclusions in Sec.
5.

##  2 Bayesian inference assisted by neural networks

In this Section, we introduce Bayesian inference of parameters assisted by a
NN. Bayesian inference is a method of statistical inference in which Bayes’
theorem is used to update the probability distributions over a set of $M$
target parameters $\Theta={\theta_{1},\theta_{2},\ldots,\theta_{M}}$ that we
aim to estimate from a certain prior hypothesis $P(\Theta)$ given a set of
experimental observations D. That is, one can apply Bayes’ theorem to find the
posterior distribution

$\displaystyle P(\Theta|\textbf{D})\propto P(\textbf{D}|\Theta)P(\Theta),$ (1)

where $P(\textbf{D}|\Theta)$ denotes the likelihood that such observations D
correspond to a set of parameters $\Theta$. Let us stress that D denotes the
experimental data from which we aim to obtain the information of the
parameters $\Theta$. For a particular case in which the observations D are the
result of measurements on a qubit, the outcomes can be either $0$ or $1$ for
the qubit being in a suitably chosen basis $|0\rangle$ or $|1\rangle$, i.e.
$\textbf{D}=\left\\{X_{j}\right\\}$ with $j=1,\ldots,N_{T}$ where $X_{j}$
denotes the number of successes obtained after measuring $N_{m}$ times at time
$t_{j}$. That is, the values $X_{j}$ from $N_{m}$ measurements at each time
$t_{j}$ form the data D. In this manner, the likelihood can be written as

$\displaystyle
P(\textbf{D}|\Theta)=\Pi_{j=1}^{N_{T}}f(X_{j},N_{m},p(\Theta)).$ (2)

The function $f$ denotes the probability of having observed exactly $X_{j}$
success outcomes from $N_{m}$ trials with success probability $p(\Theta)$,
i.e. Bernoulli distribution, where
$f(X_{j},N_{m},p)=N_{m}!/(X_{j}!(N_{m}-X_{j})!)p^{X_{j}}(1-p)^{N_{m}-X_{j}}\propto
p^{X_{j}}(1-p)^{N_{m}-X_{j}}$. Note that the success probability depends on
the target parameters $\Theta$, which in turn depends on the specific
microscopic Hamiltonian model $\hat{H}(\Theta)$ that governs the dynamics of
the system. Without lack of generality, let us assume that we monitor an
observable $\hat{A}$, which corresponds to the probability of finding the
state of a qubit in the upper state $|1\rangle$ 111Trivially, one could
consider other observables by taking into account its corresponding
probability $p(\Theta)$ in the likelihood.

$\displaystyle\langle\hat{A}(t;\Theta)\rangle={\rm{Tr}}\left[\hat{A}\hat{U}(t)\hat{\rho}_{0}\hat{U}^{\dagger}(t)\right],$
(3)

where $\hat{U}(t)=e^{-i\hat{H}(\Theta)t}$ is the time evolution operator of
the whole system under the target parameters $\Theta$ and $\hat{\rho}_{0}$
denotes the initial state ($\hbar=1$). In this manner, the likelihood can be
written as

$\displaystyle P(\textbf{D}|\Theta)=\Pi_{j=1}^{N_{T}}f(X_{j},N_{m},\langle
A(t_{j};\Theta)\rangle).$ (4)

from where the posterior follows (cf. Eq. (1)), and thus the estimation of the
target parameters.

At this stage it is important to stress that the computation of
$\langle\hat{A}(t;\Theta)\rangle$ poses a challenge to perform Bayesian
inference on generic quantum many-body systems. Indeed, Bayesian inference
requires the computation of $\langle\hat{A}(t;\Theta)\rangle$ for each of the
candidate values $\Theta$. Since integrable quantum many-body systems are
scarce, a closed-form solution for $\langle\hat{A}(t;\Theta)\rangle$ is not
available. Thus, obtaining Eq. (4) leads to a large bottleneck, if not
completely unfeasible even for moderate system sizes, which hinders the
application of estimation protocols for general quantum many-body systems.
Moreover, note that $\langle\hat{A}(t;\Theta)\rangle$ requires a microscopic
modelling $\hat{H}(\Theta)$ that accurately describes the system including all
sort of potential experimental imperfections, as otherwise the target
parameters may be poorly inferred.

Figure 1:  Sketch of the estimation strategy for quantum many-body systems.
The calibration stage consists in training a NN that faithfully reproduces the
dynamics of a chosen observable $\hat{A}$ for different candidate target
parameters $\Theta$. This allows to use the quantum many-body system as a
block-box sensor, lifting the need for an accurate microscopic model and the
computation of its dynamics, which becomes unfeasible even for a few number of
particles. Upon calibration, a set of observations on the quantum many-body
system ${\bf D}$ together with the NN enables an efficient computation of the
likelihood, and so the estimation of target parameters via Bayesian inference.
This hybrid scheme can improve the precision with respect to the standard
quantum limit.

Fortunately, NNs can help to simplify the computation and reproduce the
microscopic dynamics with high accuracy. Let us denote $F(\Theta)$ the action
of a NN which takes a parameter array $\Theta$ of dimension $M$ as input, so
that $F(\Theta)$ outputs the an array ${\bf Y}$ of dimension $N_{T}$
corresponding to the expectation value of the chosen observable $\hat{A}$ at
times $t_{j}$ ($j=1,...,N_{T}$). The suitably trained NN establishes the
relation between the inputs $\Theta=\\{\theta_{1},...,\theta_{M}\\}$ and the
outputs ${\bf Y}$, i.e. we aim at finding

$\displaystyle{\bf Y}_{j=1,\ldots,N_{T}}\equiv
F(\Theta)_{j=1,\ldots,N_{T}}\approx\langle\hat{A}(t_{j=1,\ldots,N_{T}};\Theta)\rangle$
(5)

within a tolerant error. As a consequence, the likelihood can be rewritten as

$\displaystyle
P(\textbf{D}|\Theta)=\Pi_{j=1}^{N_{T}}f(X_{j},N_{m},F(\Theta)_{j,\ldots,N_{T}}).$
(6)

Once $F(\Theta)$ is established by the NN, its evaluation is highly efficient,
which significantly decreases the computational cost. As customary in hybrid
quantum-classical algorithms, the observations ${\bf D}$ stem from
measurements on a quantum system, while the data processing is done in
classical computer. In Fig. 1, we show a sketch of our hybrid estimation
scheme, where NNs play the role of imitating the microscopic dynamics, that
enables a fast computation of the likelihood and thus to do parameter
estimation. Note that the NN is able to reproduce the microscopic dynamics of
the many-body system for a chosen observable, taken into account all
experimental imperfections and for any system size. Lifting the requirement of
possessing an accurate microscopic description $\hat{H}(\Theta)$, this hybrid
scheme comes at the cost of a calibration step, where the many-body system
response to different candidate target parameters $\Theta$ must be scanned to
then train the NN. Such calibration stage may be repeated at later times to
ensure a proper functioning of the NN. Upon the calibration, the considered
quantum many-body system can be used as black-box sensor allowing for an
efficient Bayesian parameter estimation, while still profiting from quantum
resources to boost the precision and overcoming the unfeasible computation of
the dynamics for moderate and large system sizes. Note, however, that an
improvement in the precision thanks to quantum resources depends on different
factors [3], such as the initial state $\hat{\rho}_{0}$, measured observable,
specific quantum many-body system $\hat{H}(\Theta)$, experimental
imperfections and decoherence.

Figure 2:  (a-d) Training results of the NN to establish the relation
$F_{1}(g_{x})_{j=1,...,N_{T}}\approx\langle\hat{A}_{1}(t_{j=1,...,N_{T}};g_{x})\rangle$
in the presence of 1D magnetic field using a XXZ model with $N=8$ spins. The
NN is trained based on the data directly derived from the exact expectation
value of the observable $\hat{A}_{1}$: (a) Cost function or mean squared
error; (b) gradient of the training set at different epoches; at the 35th
epoch (c) histogram and (d) Regression of the outputs of the NN trained by the
data derived without shot noise with respect to the ideal targets
$\textbf{T}=\\{\langle\hat{A}_{1}(t_{1};g_{x})\rangle,...,\langle\hat{A}_{1}(t_{101};g_{x})\rangle\\}$.
(e) Regression of the outputs of the NN trained based on the data obtained
from the numerical simulation including finite number of measurements, here
$N_{m}=100$, with respect to the ideal targets T. In both (d) and (e), the fit
lines (solid black) given in the equations shown below each plot, almost
overlap with the line $\textbf{Y}=\textbf{T}$, while $R$ is the correlation
coefficient of the outputs and the targets. (f) Comparison of the outputs from
the NNs trained by the data derived without shot noise (dashed red) and
including shot noise $N_{m}=100$ (solid blue) in accordance with the numerical
ideal results $\langle\hat{A}_{1}(t_{j=1,...,N_{T}};g_{x})\rangle$ from the
Hamiltonian $\hat{H}_{1}$ (dotted black), where $g_{x}=0.1253$.

##  3 Simulation of quantum many-body dynamics by neural networks

In order to illustrate the capability of NNs to reproduce the dynamics of a
chosen observable for a quantum many-body system, we consider a XXZ spin-1/2
chain in the presence of an external field that we aim to estimate. We stress
that the specific many-body system is irrelevant for the proposed hybrid
estimation protocol. However, we choose the XXZ spin-1/2 chain since this
model has been shown to display a precision scaling beyond the quantum
standard limit, see Ref. [18]. This quantum many-body system constitutes
therefore a good example to test whether the proposed estimation strategy can
still render a quantum enhancement, which we will discuss in Sec. 4. In the
following we show two examples to test the performance of NNs to reproduce the
quantum many-body dynamics of interest, namely, for a one and two dimensional
external field. In addition, in A we show how the NNs are able to accurately
reproduce the dynamics under potential experimental imperfections, such as
when decoherence effects are present.

###  3.1 1D external magnetic field

We first start from the ferromagnetic XXZ spin-1/2 for $N$ particles with open
boundary conditions as a quantum many-body system in the presence of a 1D
external field. The Hamiltonian of the system read as

$\displaystyle\hat{H}_{1}(g_{x})=-\sum_{i=1}^{N-1}(\hat{\sigma}_{i}^{x}\hat{\sigma}_{i+1}^{x}+\hat{\sigma}_{i}^{y}\hat{\sigma}_{i+1}^{y}+J_{z}\sigma_{i}^{z}\sigma_{i+1}^{z})+g_{x}\sum_{i=1}^{N}\hat{\sigma}_{i}^{x}.$
(7)

This model displays a critical point at $g_{x}=0$ [58, 18]. In the following
we consider the initial state $\hat{\rho}_{0}$ as the ground state without the
external field $g_{x}=0$, i.e. at the critical point. Moreover, we take
$J_{z}=0$ as considered in Ref. [18]. The strength of the external field can
be estimated by measuring the response delivered by the system, in particular,
we choose a local observable $\hat{A}_{1}=(\hat{\sigma}_{N/2}^{x}+1)/2$, so
that $\Theta=\\{g_{x}\\}$. Our aim is to build a NN that encodes the relation

$\displaystyle
F_{1}(g_{x})_{j=1,...,N_{T}}\approx\langle\hat{A}_{1}(t_{j=1,...,N_{T}};g_{x})\rangle$
(8)

from a finite number of measurements, i.e. without having access to the exact
expectation value $\langle\hat{A}_{1}(t_{j=1,...,N_{T}};g_{x})\rangle$ as it
would be the case in any experimental realization. This corresponds to the
calibration stage in Fig. 1.

Yet, before proceeding further, it is convenient to first consider an ideal
scenario where our NN can be trained using the exact expectation value
$\langle\hat{A}_{1}(t_{j=1,...,N_{T}};g_{x})\rangle$. This will allow us to
quantify the deviations from this ideal scenario when training the NN with
shot noise measurements. To generate the dataset, we derive the expectation
value $\langle\hat{A}_{1}(t_{j=1,...,N_{T}};g_{x})\rangle$ as the target data
during the time interval $0\leq t\leq t_{f}=5$ at $N_{T}=101$ time instants,
i.e., $t\in\\{0,0.05,...,5\\}$. As the inputs of the NN, the magnetic field is
in the range $0\leq g_{x}\leq 0.5$ with $51$ examples such that
$g_{x}\in\\{0,0.01,...,0.5\\}$. The whole dataset is then composed by $n=51$
examples, each of which has the output data
$\textbf{Y}=\\{y_{1},y_{2},...,y_{101}\\}$ approaching the targets
$\textbf{T}=\\{a_{1},...,a_{101}\\}=\\{\langle\hat{A}_{1}(t_{1};g_{x})\rangle,...,\langle\hat{A}_{1}(t_{101};g_{x})\rangle\\}$
–a string of $101$ numbers– while the input data is $\Theta=\\{g_{x}\\}$. The
training/validation/test datasets are built with the data randomly chosen
$70\%$ / $15\%$ / $15\%$ from the total dataset. We introduce the cost
function as the mean squared error,

$\displaystyle
C=\sum_{i=1}^{n}\sum_{j=1}^{N_{T}}\frac{1}{nN_{T}}(y_{i}^{j}-a_{i}^{j})^{2}.$
(9)

The NN is adjusted by minimizing the cost function for a training set with $n$
examples employing gradient descent, that ultimately leads to
$F_{1}(g_{x})_{j=1,...,N_{T}}\approx\langle\hat{A}_{1}(t_{j=1,...,N_{T}};g_{x})\rangle$.
We take a XXZ system with $N=8$ spins as an example and measure the
expectation value
$\langle\hat{A}_{1}\rangle=(\langle\hat{\sigma}_{4}^{x}\rangle+1)/2$. Tuning
the parameters of this NN with $1$ input neuron and $101$ output neurons, we
find that four hidden layers with $6,12,25,50$ neurons with a hyperbolic-
tangent activation function [26] give the highest accuracy. The training of
such a NN stops when $C<10^{-5}$ at the 35th epoch leading to the results
plotted in Fig. 2 (a)-(d). The regression of the outputs of the NN with
respect to the targets T derived from numerical simulation is shown in Fig. 2
(d). As an illustration to corroborate the good performance of the NN in
reproducing the dynamics of the system, in Fig. 2 (f), we show the exact
expectation value $\langle\hat{A}_{1}(t_{j=1,...,N_{T}};g_{x})\rangle$ (dotted
black) for a randomly chosen example $g_{x}=0.1253$ –that does not belong to
the training/validation/test datasets– and the output of the NN (dashed red),
which overlap to a very good approximation.

Figure 3:  (a-c) Training results of the NN to establish the relation
$F_{2}(g_{x},g_{y})_{j=1,...,N_{T}}\approx\langle\hat{A}_{2}(t_{j=1,...,N_{T}};g_{x},g_{y})\rangle$
in the presence of 2D magnetic field in $x-y$ plane and $N=6$ with a finite
number of measurements, $N_{m}=100$: (a) Cost function or mean squared error;
(b) gradient of the training set; (c) Regression of the outputs of such a NN
with respect to the exact expectation values (targets)
$\textbf{T}=\langle\hat{A}_{2}(t_{j=1,...,N_{T}};g_{x},g_{y})\rangle$. The fit
line (solid black), expressed in the equations shown below, approximately
matches the ideal relation $\textbf{Y}=\textbf{T}$, while $R$ is the
correlation coefficient of the outputs and the targets. (d) Comparison of the
outputs from the NNs trained by the data with shot-noise measurements
$N_{m}=100$ (solid blue) with the exact expectation value
$\langle\hat{A}_{2}\rangle$ obtained from the Hamiltonian $\hat{H}_{2}$
(dashed red), when $g_{x}=0.0998$, $g_{y}=0.1236$.

We now turn our attention to a more realistic scenario, where the calibration
of the NN is performed using shot noise measurements. In order to simulate an
experimental acquisition, we numerically draw a set of $N_{m}$ binary
measurements $x_{m}^{j}$ with $m=1,\ldots,N_{m}$ per time instant $t_{j}$,
where $x_{m}^{j}\in\\{0,1\\}$. The outcomes are drawn from the probability of
measuring the observable $\hat{A}_{1}$ at time time $t_{j}$ and the
corresponding external field $g_{x}$, i.e.
$\langle\hat{A}_{1}(t_{j};g_{x})\rangle$. That is, $x_{m}^{j}\sim
B(1,\langle\hat{A}_{1}(t_{j};g_{x})\rangle)$ where $B(1,p)$ denotes the
Bernoulli distribution. We then use $\Theta=\\{g_{x}\\}$ as the input of the
NN, while the corresponding set of averaged values
$d_{j}=\sum_{m=1}^{N_{m}}x_{m}^{j}/N_{m}$ as the target dataset to train the
NN, i.e. ${\bf T}=\\{d_{1},\ldots,d_{N_{T}}\\}$. As in the previous case, we
consider $N_{T}=101$ time instances in the interval $0\leq t\leq 5$, so that
the number of input and output neurons of the NN are, as before, $1$ and
$101$, respectively. Training the NN with shot-noise measurements introduces
statistical fluctuations. On the one hand, for a large number of measurements
$N_{m}\ll 1$ the shot-noise averaged value converges to the exact expectation
value $d_{j}\approx\langle\hat{A}_{1}(t_{j};g_{x})\rangle$. A large number of
measurements, however, requires long experimental runs and thus leads to a
large overhead in experimental resources. On the other hand, a reduced number
of measurements may speed up the calibration stage at the expense of
introducing further errors in the approximation of the dynamics by the NN, and
thus in the estimation precision. In the following we set a reasonable number
of measurements per time instant, $N_{m}=100$, for which the training provides
good results. We have tested cases with other $N_{m}$; in particular for
$N_{m}=50$ one still finds a good accuracy of the NN, while for $N_{m}=20$
measurements, shot noise dominates and it deteriorates the performance of the
NN. Moreover, in order to mitigate the impact of statistical fluctuations in
the NN, we perform $10$ repetitions per calibrated value of $g_{x}$, so the
whole dataset contains $n=51\times 10$ examples from which $70\%$, $15\%$,
$15\%$ represent the training, validation and test sets. The resulting
regression of the NN considering a finite number of measurements is shown in
Fig. 2 (e). The NN trained with data including shot noise provides a high
accuracy regression with respect to the ideal targets
$\langle\hat{A}_{1}(t_{j=1,...,N_{T}};g_{x})\rangle$. Taking a randomly chosen
value $g_{x}=0.1253$, the outputs from this NN (solid blue) coincide well with
the exact numerical result calculated from the Hamiltonian $\hat{H}_{1}$
(dotted black) (cf. Fig. 2 (f)). The NN is therefore capable of establishing
the relation given in Eq. (8), so that the NN can be used to compute the
likelihood, cf. Eq. (6), and thus to perform Bayesian inference.

###  3.2 2D external magnetic field

We consider now an extension of the previous case, namely, a ferromagnetic XXZ
spin-1/2 chain under a 2D magnetic field in $x-y$ plane we want to estimate.
The Hamiltonian of the system can be written as

$\displaystyle\hat{H}_{2}(g_{x},g_{y})=-\sum_{i=1}^{N-1}(\hat{\sigma}_{i}^{x}\hat{\sigma}_{i+1}^{x}+\hat{\sigma}_{i}^{y}\hat{\sigma}_{i+1}^{y}+J_{z}\hat{\sigma}_{i}^{z}\hat{\sigma}_{i+1}^{z})+g_{x}\sum_{i=1}^{N}\hat{\sigma}_{i}^{x}+g_{y}\sum_{i=1}^{N}\hat{\sigma}_{i}^{y}.$
(10)

Again, the initial state $\hat{\rho}_{0}$ is taken as the ground state for
$g_{x}=g_{y}=0$, and we consider $J_{z}=0$. Letting the system evolve from the
initial state for some time $0\leq t\leq 5$, we measure the local observable
$\hat{A}_{2}=(\hat{\sigma}^{y}_{N/2}+1)/2$ for different values of $g_{x}$ and
$g_{y}$ in the range $0\leq g_{x}\leq 0.5$ and $0\leq g_{y}\leq 0.5$. For the
considered examples, the local magnetization along the $y$ axis provides a
better sensitivity than along $x$ or $z$ axes.

For the calibration of the NN we proceed as before, including shot-noise
measurements. Recall, however, that now the inputs $\Theta=\\{g_{x},g_{y}\\}$
require $2$ neurons for the input layer of the NN. We train the NN so that the
outputs ${\bf Y}$ reproduce the targets that are the averaged value shot-noise
measurements of $\hat{A}_{2}$, $\textbf{T}=\\{d_{1},\ldots,d_{N_{T}}\\}$. As
before, we consider $N_{T}=101$, while the total number of datasets becomes
$n=51^{2}\times 10$ due to the estimated parameters in two dimension and the
$10$ repetitions for mitigating the impact of shot-noise fluctuations. In this
case, the trained NNs allow us to find the relation

$\displaystyle F_{2}(g_{x},g_{y})_{j=1,...,N_{T}}\approx\langle
A_{2}(t_{j=1,...,N_{T}};g_{x},g_{y})\rangle.$ (11)

See Fig. 3 for the NN training results to reproduce the dynamics of the system
with $N=6$ spins for different values of $g_{x}$ and $g_{y}$.

##  4 Estimation precision

The quantum Fisher information (QFI) [4, 5] provides the ultimate bound to the
achievable precision of a parameter $\theta$ encoded in a quantum state
$|\Psi_{\theta}\rangle$ over all possible quantum measurements. The QFI is
defined as

$\displaystyle
I_{\theta}^{\rm{QFI}}=4\left(\left\langle\partial_{\theta}\Psi|\partial_{\theta}\Psi\right\rangle-|\langle\Psi|\partial_{\theta}\Psi\rangle|^{2}\right),$
(12)

where $|\partial_{\theta}\Psi\rangle$ denotes the derivative of the state with
respect to the parameter $\theta$. The minimum variance for an unbiased
estimator of $\theta$ for a single measurement is set by the Cramer-Rao bound,

$\displaystyle(\Delta\theta)^{2}\geq(\Delta\theta^{\rm{QFI}})^{2}=\frac{1}{I_{\theta}^{\rm
QFI}}.$ (13)

The scaling with the number of probes or particles $N$ is determined by the
QFI. Recall that standard quantum limit fixes the precision scaling as
$\Delta\theta\propto N^{-1/2}$, while the Heisenberg limit refers to a better
scaling of the form $\Delta\theta\propto N^{-1}$.

In Sec. 3 we have shown how the NN allows us to faithfully reproduce the
dynamics of the quantum many-body system upon a calibration stage using finite
number measurements. Hence, we can compute the likelihood from a set of
observations ${\bf D}$, and so the posterior distribution over the relevant
parameters, $\Theta=\\{g_{x}\\}$ or $\Theta=\\{g_{x},g_{y}\\}$ for the one or
two-dimensional external field, respectively. The posterior distribution
follows from the Bayes’ theorem, $P(\Theta|{\bf D})\propto P({\bf
D}|\Theta)P(\Theta)$, and we take the prior distribution as flat or
uninformative in the considered range for $g_{x}$ and $g_{y}$ (cf. Sec. 3).
The estimated parameter and its variance can be computed from its marginal
posterior distribution, i.e. for $\theta_{j}$ we have as
$\theta_{j}^{\rm{est}}=\int d\theta_{j}\theta_{j}P(\theta_{j}|\textbf{D})$ and
$(\Delta\theta_{j}^{\rm{est}})^{2}=\int
d\theta_{j}(\theta_{j}-\theta_{j}^{\rm{est}})^{2}P(\theta_{j}|\textbf{D})$,
where the marginal distribution is given by
$P(\theta_{j}|\textbf{D})=\int\Pi_{i\neq j}d\theta_{i}P(\Theta|\textbf{D})$.

Figure 4:  (a) Posterior probability distributions of the one-dimensional
external field strength $g_{x}$ for a $N=4$ (solid, blue), $N=6$ (dotted,
black), $N=8$ (red, dashed) and $N=10$ (dot-dashed, magenta) computing
$P(\textbf{D}|g_{x})$ assisted by the NN. The target value corresponds to
$g_{x}=0.1$. (b-d) The scaling of estimation precision (squared, solid blue).
The red dashed lines provide a guide for the eye with various scalings.

Figure 5:  (a) Posterior probability distributions for $N=6$ in terms of
$g_{x}=0.1$ and $g_{y}=0.1$ using the likelihood $P(\textbf{D}|g_{x},g_{y})$
assisted by the NN. (b-d) The scaling of estimation precision of
$\overline{\Delta g_{x}^{\rm{est}}}$ (squared, solid blue) and
$\overline{\Delta g_{y}^{\rm{est}}}$ (diamond, solid magenta) for different
$g_{x}$ and $g_{y}$, where the value with each $N$ is derived from Bayesian
inference with the likelihood assisted by the NN. The red dashed lines
indicate the standard quantum and Heisenberg limit scalings, i.e. $N^{-0.5}$
and $N^{-1}$.

Employing the trained NNs, we compute the resulting mean and standard
deviations for different target values $g_{x}$ (and $g_{y}$ for the two-
dimensional case) simulating an experimental acquisition ${\bf D}$. This
experimental acquisition is simulated in the same manner as when generating
the datasets for the calibration step, cf. Sec. 3, albeit with an arbitrary
number of measurements $N_{p}$ per time instant. Without lack of generality,
we set $N_{p}=N_{m}=100$. Note that the precision scaling with the number of
particles $N$ is not significantly altered by considering different number of
measurements $N_{p}$. This results in $g_{x}^{\rm est}$ and $\Delta g_{x}^{\rm
est}$ (and $g_{y}^{\rm est}$ and $\Delta g_{y}^{\rm est}$). Moreover, in order
to compute a statistically significant precision, we average the resulting
standard deviation for $10$ different observations under the same target
values. In this manner we find $\overline{\Delta g_{x}^{\rm est}}$ (and
$\overline{\Delta g_{y}^{\rm est}}$ for the two-dimensional case).

In Fig. 4 we summarize the results to determine the strength of a one-
dimensional external field $g_{x}$ using the NNs trained with the datasets
with $N_{m}=100$ measurements. In Fig. 4(a) we show the posterior distribution
for an increasing number of spins for $g_{x}=0.1$, while the other panels of
Fig. 4 show the scaling of the averaged standard deviation for three different
values of $g_{x}$, i.e. $0.05$, $0.1$ and $0.15$. The specific values for
$g_{x}^{\rm est}$ and $\Delta g_{x}^{\rm est}$, as well as their averaged
values over $10$ repetitions can be found in B. Remarkably, the behavior of
the standard deviation for small sizes indicates a Heisenberg ($\Delta
g_{x}\propto N^{-1}$) or even a super-Heisenberg scaling ($\Delta g_{x}\propto
N^{-\alpha}$ with $\alpha>1$), although the scaling breaks down for some $N$.
Note that the larger $g_{x}$, the smaller the size where the scaling breaks.
These issues have been previously discussed [18] for the same model. Indeed,
as $g_{x}$ gets larger the system is pushed out of the critical region
($g_{x}=0$) where the sensitivity is higher, reducing the precision as more
probes are included to finally saturate (see [18] for more details). In
addition, it is worth noting that fixing a maximum time $t_{f}$ for the
dynamical protocol ($0\leq t\leq t_{f}=5$) also limits the precision for
magnetic-field strengths such that $t_{f}g_{x}\ll 1$.

For the estimation of the strength of a two-dimensional magnetic field we
proceed in a similar manner. In Fig. 5(a) we show the bi-dimensional posterior
distribution $P(g_{x},g_{y}|{\bf D})$ for $N=6$ spins with $g_{x}=g_{y}=0.1$.
In the other panels of Fig. 5, the behavior of the averaged standard deviation
with an increasing number of spins is shown. As for the one-dimensional case,
the model presents an improvement in the precision for small sizes and
strengths, which is captured by our inference protocol. The numerical values
can also be found in B.

Besides these expected issues, it is important to stress that this hybrid
metrological scheme incorporating the NN to perform Bayesian inference is able
to capture the quantum-enhanced scaling in the precision as long as the
quantum many-body system used as a sensor allows for it.

##  5 Conclusion

In this article we propose a protocol to design a quantum many-body
magnetometer using Bayesian inference techniques based on the ability of
neural networks to faithfully reproduce the dynamics of the quantum many-body
system. Bayesian inference techniques require the computation of the
microscopic dynamics of the quantum many-body system under different
parameters. For a general quantum many-body system, the computation of the
dynamics soon becomes impractical even for small system sizes. To overcome
this obstacle, we make use of a suitably trained neural network that is able
to reproduce its dynamics. The neural network enables then an efficient
Bayesian inference of external magnetic fields. In this manner, the neural
network circumvents the need for a microscopic model description of the system
at the expense of requiring a calibration step to establish a relation between
dynamics and external magnetic fields. The calibration of the NNs could be
performed solely based on experimental data, thus opening the door to employ
experimental quantum many-body systems as sensors that are beyond the
classical simulation capabilities. We illustrate this procedure in a XXZ model
that allows for Heisenberg or even super-Heisenberg scaling in the precision
in a parameter regime. Our numerical results show that the external field
strength can be obtained with high accuracy even for a finite number of
measurements. This protocol showcases the advantage of Bayesian inference
making use neural networks for parameter estimation, and thus paves the way to
use quantum many-body systems as sensors under realistic conditions while
benefiting from quantum resources to improve their precision.

## Acknowledgments

J. C. acknowledges the Ramón y Cajal (RYC2018-025197-I) research fellowship,
the financial support from Spanish Government via EUR2020-112117 and Nanoscale
NMR and complex systems (PID2021-126694NB-C21) projects, the EU FET Open Grant
Quromorphic (828826), and the Basque Government grant IT1470-22. This work was
supported in part by the Basque Government through the ELKARTEK Program under
Grant KK-2022/00062, “QFIRST - Dispositivos en Tecnologías Cuánticas, and
under Grant KK-2022/00041, “BRTA QUANTUM: Hacia una especialización armonizada
en tecnologías cuánticas en BRTA”.

##  Appendix A Neural network reproducing the dynamics including decoherence

Figure 6:  (a-c) Training results of the NN to establish the relation
$F_{1}(g_{x})_{j=1,...,N_{T}}\approx\langle\hat{A}_{1}(t_{j=1,...,N_{T}};g_{x})\rangle$
for 1D magnetic field and $N=4$ with a finite number of measurements
$N_{m}=100$ and decoherence with rate $\gamma=1/10$: (a) mean squared error;
(b) gradient of the training set; (c) regression of the outputs of such a NN
with respect to the ideal targets
$\textbf{T}=\langle\hat{A}_{1}(t_{j=1,...,N_{T}};g_{x})\rangle$. The fit line
(solid black) demonstrates that $\textbf{Y}=\textbf{T}$ to a good
approximation, while $R$ is the correlation coefficient of the outputs and the
targets. (d) Comparison of the outputs from the NNs trained by the data with
$N_{m}=100$ (solid blue) in accordance with the exact numerical calculation
$\langle\hat{A}_{1}(t;g_{x})\rangle$ (dashed red) from the Lindblad master
equation, Eq. (14) when $g_{x}=0.2053$.

As an example to stress the generality of the neural network to reproduce the
dynamics of the quantum many-body system with imperfections, we illustrate the
performance when including decoherence. In this case, the considered that
dynamics of the system is described by means of a Lindblad master equation
[59, 60]

$\displaystyle\frac{d\hat{\rho}}{dt}=-i[\hat{H},\hat{\rho}]+\sum_{k}\gamma_{k}\left(\hat{L}_{k}\hat{\rho}\hat{L}_{k}^{\dagger}-\frac{1}{2}\left\\{\hat{L}_{k}^{\dagger}\hat{L}_{k},\hat{\rho}\right\\}\right).$
(14)

For simplicity, we consider equal decoherence rates $\gamma=\gamma_{k}$ for
all qubits and the spin operator $\hat{L}_{k}=\hat{\sigma}_{-}$. Due to the
large computational cost to simulate Eq. (14), we access system with $N=4$ and
$6$ spins for the 1D magnetic field, setting a moderate decoherence rate
$\gamma=1/10$. In both cases the trained NN reproduces the exact dynamics to a
very good approximation. The training process of such a NN is equivalent to
one described in Sec. 3. In Fig. 6, we show the training results of
performance, gradients and regression of the outputs of the NN with respect to
the exact numerical results including decoherence. Inputting one randomly
chosen example $g_{x}=0.2053$, which does not belong to any dataset, we find
the outputs from the NN that corresponds to the targets
$\langle\hat{A}_{1}(t_{j=1,...,N_{T}};g_{x})\rangle$ to a very good
approximation. Note that the precision in estimating $g_{x}$ may be spoiled
depending on the decoherence rate $\gamma$, but the NN still allows for an
efficient Bayesian estimation. This hybrid scheme can achieve quantum-enhanced
scaling in the precision as long as the quantum many-body system used as a
sensor allows for it.

##  Appendix B Numerical results

In the following Tab. 1, 2 and 3 we provide the numerical values that are
plotted in Figs. 4 and 5 for the estimated magnetic-field strengths $g_{x}$,
$g_{y}$, their corresponding standard deviations as well as their averaged
values over $10$ repetitions.

Target $g_{x}$ | $N$ | $g_{x}^{\rm{est}}$ | $\Delta g_{x}^{\rm{est}}$
---|---|---|---
0.1 | 4 | 0.09729 | 0.004343
| 6 | 0.09749 | 0.002931
| 8 | 0.09817 | 0.001949
| 10 | 0.09824 | 0.002105
Table 1:  Estimated magnetic field strength $g_{x}$ and its standard deviation $\Delta g_{x}^{\rm{est}}$ for $N=4,6,8,10$ spins obtained from Bayesian inference employing the trained NN with a finite-number of measurements. The simulated observations ${\bf D}$ have been obtained for a target value $g_{x}=0.1$. Target $g_{x}$ | $N$ | $\overline{g_{x}^{\rm{est}}}$ | $\overline{\Delta g_{x}^{\rm{est}}}$
---|---|---|---
0.05 | 4 | 0.04788 | 0.004513
| 6 | 0.04872 | 0.002251
| 8 | 0.04769 | 0.002376
| 10 | 0.04931 | 0.002239
0.1 | 4 | 0.1005 | 0.004513
| 6 | 0.1005 | 0.002850
| 8 | 0.1006 | 0.002062
| 10 | 0.1003 | 0.002238
0.15 | 4 | 0.1529 | 0.004112
| 6 | 0.1511 | 0.002566
| 8 | 0.1515 | 0.002444
| 10 | 0.1498 | 0.002973
0.2 | 4 | 0.2016 | 0.003225
| 6 | 0.2001 | 0.002625
| 8 | 0.2013 | 0.003672
| 10 | 0.2015 | 0.004255
Table 2:  Averaged values for the estimated magnetic-field strength $\overline{g_{x}^{\rm{est}}}$ and for its standard deviation $\overline{\Delta g_{x}^{\rm{est}}}$ for $10$ repetitions and for different system sizes. Target $g_{x}$, $g_{y}$ | $N$ | $\overline{g_{x}^{\rm{est}}}$ | $\overline{g_{y}^{\rm{est}}}$ | $\overline{\Delta g_{x}^{\rm{est}}}$ | $\overline{\Delta g_{y}^{\rm{est}}}$
---|---|---|---|---|---
0.05, 0.05 | 4 | 0.04153 | 0.04909 | 0.01859 | 0.003323
| 6 | 0.04752 | 0.04999 | 0.004651 | 0.002452
| 8 | 0.04798 | 0.04881 | 0.002235 | 0.002088
| 10 | 0.04923 | 0.04983 | 0.001535 | 0.002327
0.075, 0.075 | 4 | 0.06909 | 0.07339 | 0.01062 | 0.003817
| 6 | 0.07590 | 0.07662 | 0.002804 | 0.002991
| 8 | 0.07505 | 0.07477 | 0.002509 | 0.003322
| 10 | 0.07343 | 0.07687 | 0.003050 | 0.003502
0.1, 0.1 | 4 | 0.09779 | 0.09799 | 0.004168 | 0.003942
| 6 | 0.1002 | 0.09752 | 0.003011 | 0.003913
| 8 | 0.1023 | 0.09630 | 0.003542 | 0.004176
| 10 | 0.1022 | 0.09813 | 0.004438 | 0.003917
Table 3:  Averaged values for the estimated magnetic-field strengths
$\overline{g_{x}^{\rm{est}}}$, $\overline{g_{y}^{\rm{est}}}$ and for their
standard deviations $\overline{\Delta g_{x}^{\rm{est}}}$, $\overline{\Delta
g_{y}^{\rm{est}}}$.

## References

## References

  * [1] Giovannetti V, Lloyd S and Maccone L 2004 Science 306 1330–1336 (Preprint https://www.science.org/doi/pdf/10.1126/science.1104149) URL https://www.science.org/doi/abs/10.1126/science.1104149
  * [2] Giovannetti V, Lloyd S and Maccone L 2006 Phys. Rev. Lett. 96(1) 010401 URL https://link.aps.org/doi/10.1103/PhysRevLett.96.010401
  * [3] Giovannetti V, Lloyd S and Maccone L 2011 Nat. Phot. 5 222–229 ISSN 1749-4893 URL https://doi.org/10.1038/nphoton.2011.35
  * [4] Braunstein S L and Caves C M 1994 Phys. Rev. Lett. 72(22) 3439–3443 URL https://link.aps.org/doi/10.1103/PhysRevLett.72.3439
  * [5] Paris M G A 2009 Int. J. Quant. Inf. 07 125–137 (Preprint https://doi.org/10.1142/S0219749909004839) URL https://doi.org/10.1142/S0219749909004839
  * [6] Leibfried D, Barrett M D, Schaetz T, Britton J, Chiaverini J, Itano W M, Jost J D, Langer C and Wineland D J 2004 Science 304 1476–1478 (Preprint https://www.science.org/doi/pdf/10.1126/science.1097576) URL https://www.science.org/doi/abs/10.1126/science.1097576
  * [7] Huelga S F, Macchiavello C, Pellizzari T, Ekert A K, Plenio M B and Cirac J I 1997 Phys. Rev. Lett. 79(20) 3865–3868 URL https://link.aps.org/doi/10.1103/PhysRevLett.79.3865
  * [8] Boixo S, Flammia S T, Caves C M and Geremia J 2007 Phys. Rev. Lett. 98(9) 090401 URL https://link.aps.org/doi/10.1103/PhysRevLett.98.090401
  * [9] Roy S M and Braunstein S L 2008 Phys. Rev. Lett. 100(22) 220501 URL https://link.aps.org/doi/10.1103/PhysRevLett.100.220501
  * [10] Pang S and Jordan A N 2017 Nat. Commun. 8 14695 ISSN 2041-1723 URL https://doi.org/10.1038/ncomms14695
  * [11] Jones G S, Bose S and Bayat A 2020 arXiv:2003.02308 URL https://arxiv.org/abs/2003.02308
  * [12] Montenegro V, Jones G S, Bose S and Bayat A 2022 Phys. Rev. Lett. 129(12) 120503 URL https://link.aps.org/doi/10.1103/PhysRevLett.129.120503
  * [13] Mishra U and Bayat A 2022 Sci. Rep. 12 14760 ISSN 2045-2322 URL https://doi.org/10.1038/s41598-022-17381-y
  * [14] Yang J, Pang S, Chen Z, Jordan A N and del Campo A 2022 Phys. Rev. Lett. 128(16) 160505 URL https://link.aps.org/doi/10.1103/PhysRevLett.128.160505
  * [15] Zanardi P, Paris M G A and Campos Venuti L 2008 Phys. Rev. A 78 042105 URL https://link.aps.org/doi/10.1103/PhysRevA.78.042105
  * [16] Boyajian W L, Skotiniotis M, Dür W and Kraus B 2016 Phys. Rev. A 94(6) 062326 URL https://link.aps.org/doi/10.1103/PhysRevA.94.062326
  * [17] Fernández-Lorenzo S and Porras D 2017 Phys. Rev. A 96(1) 013817 URL https://link.aps.org/doi/10.1103/PhysRevA.96.013817
  * [18] Rams M M, Sierant P, Dutta O, Horodecki P and Zakrzewski J 2018 Phys. Rev. X 8(2) 021022 URL https://link.aps.org/doi/10.1103/PhysRevX.8.021022
  * [19] Garbe L, Bina M, Keller A, Paris M G A and Felicetti S 2020 Phys. Rev. Lett. 124(12) 120504 URL https://link.aps.org/doi/10.1103/PhysRevLett.124.120504
  * [20] Ivanov P A 2020 Phys. Rev. A 102(5) 052611 URL https://link.aps.org/doi/10.1103/PhysRevA.102.052611
  * [21] Chu Y, Zhang S, Yu B and Cai J 2021 Phys. Rev. Lett. 126(1) 010502 URL https://link.aps.org/doi/10.1103/PhysRevLett.126.010502
  * [22] Gietka K, Ruks L and Busch T 2022 Quantum 6 700 ISSN 2521-327X URL https://doi.org/10.22331/q-2022-04-27-700
  * [23] Garbe L, Abah O, Felicetti S and Puebla R 2022 Quantum Sci. Technol. 7 035010 URL https://dx.doi.org/10.1088/2058-9565/ac6ca5
  * [24] Garbe L, Abah O, Felicetti S and Puebla R 2022 Phys. Rev. Res. 4(4) 043061 URL https://link.aps.org/doi/10.1103/PhysRevResearch.4.043061
  * [25] Ilias T, Yang D, Huelga S F and Plenio M B 2022 PRX Quantum 3(1) 010354 URL https://link.aps.org/doi/10.1103/PRXQuantum.3.010354
  * [26] Bishop C M 2006 Pattern Recognition and Machine Learning (Springer)
  * [27] Carleo G, Cirac I, Cranmer K, Daudet L, Schuld M, Tishby N, Vogt-Maranto L and Zdeborová L 2019 Rev. Mod. Phys. 91(4) 045002 URL https://link.aps.org/doi/10.1103/RevModPhys.91.045002
  * [28] Torlai G and Melko R G 2020 Annual Review of Condensed Matter Physics 11 325–344 (Preprint https://doi.org/10.1146/annurev-conmatphys-031119-050651) URL https://doi.org/10.1146/annurev-conmatphys-031119-050651
  * [29] Krenn M, Landgraf J, Foesel T and Marquardt F 2022 arXiv:2208.03836 URL https://arxiv.org/abs/2208.03836
  * [30] Lumino A, Polino E, Rab A S, Milani G, Spagnolo N, Wiebe N and Sciarrino F 2018 Phys. Rev. Appl. 10(4) 044033 URL https://link.aps.org/doi/10.1103/PhysRevApplied.10.044033
  * [31] Xiao T, Huang J, Fan J and Zeng G 2019 Sci. Rep. 9 12410 ISSN 2045-2322 URL https://doi.org/10.1038/s41598-019-48551-0
  * [32] Palittapongarnpim P and Sanders B C 2019 Phys. Rev. A 100(1) 012106 URL https://link.aps.org/doi/10.1103/PhysRevA.100.012106
  * [33] Xu H, Li J, Liu L, Wang Y, Yuan H and Wang X 2019 npj Quantum Information 5 82 ISSN 2056-6387 URL https://doi.org/10.1038/s41534-019-0198-z
  * [34] Peng Y and Fan H 2020 Phys. Rev. A 101(2) 022107 URL https://link.aps.org/doi/10.1103/PhysRevA.101.022107
  * [35] Schuff J, Fiderer L J and Braun D 2020 New J. Phys. 22 035001 URL https://dx.doi.org/10.1088/1367-2630/ab6f1f
  * [36] Fiderer L J, Schuff J and Braun D 2021 PRX Quantum 2(2) 020303 URL https://link.aps.org/doi/10.1103/PRXQuantum.2.020303
  * [37] Xiao T, Fan J and Zeng G 2022 npj Quantum Information 8 2 ISSN 2056-6387 URL https://doi.org/10.1038/s41534-021-00513-z
  * [38] Madsen C N, Valdetaro L and Mølmer K 2021 Phys. Rev. A 104(5) 052621 URL https://link.aps.org/doi/10.1103/PhysRevA.104.052621
  * [39] Cimini V, Gianani I, Spagnolo N, Leccese F, Sciarrino F and Barbieri M 2019 Phys. Rev. Lett. 123(23) 230502 URL https://link.aps.org/doi/10.1103/PhysRevLett.123.230502
  * [40] Nolan S, Smerzi A and Pezzè L 2021 npj Quantum Information 7 169 ISSN 2056-6387 URL https://doi.org/10.1038/s41534-021-00497-w
  * [41] Niu M Y, Boixo S, Smelyanskiy V N and Neven H 2019 npj Quantum Information 5 33 ISSN 2056-6387 URL https://doi.org/10.1038/s41534-019-0141-3
  * [42] Bukov M, Day A G R, Sels D, Weinberg P, Polkovnikov A and Mehta P 2018 Phys. Rev. X 8(3) 031086 URL https://link.aps.org/doi/10.1103/PhysRevX.8.031086
  * [43] Brown J, Sgroi P, Giannelli L, Paraoanu G S, Paladino E, Falci G, Paternostro M and Ferraro A 2021 New J. Phys. 23 093035 URL https://dx.doi.org/10.1088/1367-2630/ac2393
  * [44] Borah S, Sarma B, Kewming M, Milburn G J and Twamley J 2021 Phys. Rev. Lett. 127(19) 190403 URL https://link.aps.org/doi/10.1103/PhysRevLett.127.190403
  * [45] Porotti R, Essig A, Huard B and Marquardt F 2022 Quantum 6 747 ISSN 2521-327X URL https://doi.org/10.22331/q-2022-06-28-747
  * [46] Ding Y, Ban Y, Martín-Guerrero J D, Solano E, Casanova J and Chen X 2021 Phys. Rev. A103(4) L040401 URL https://doi.org/10.1103/PhysRevA.103.L040401
  * [47] Ai M Z, Ding Y, Ban Y, Martín-Guerrero J D, Casanova J, Cui J M, Huang Y F, Chen X, Li C F and Guo G C 2022 Sci. China Phys. Mech. Astron.  65 250312 URL https://doi.org/10.1007/s11433-021-1841-2
  * [48] Krenn M, Malik M, Fickler R, Lapkiewicz R and Zeilinger A 2016 Phys. Rev. Lett. 116(9) 090405 URL https://link.aps.org/doi/10.1103/PhysRevLett.116.090405
  * [49] Melnikov A A, Nautrup H P, Krenn M, Dunjko V, Tiersch M, Zeilinger A and Briegel H J 2018 PNAS 115 1221–1226 (Preprint https://www.pnas.org/doi/pdf/10.1073/pnas.1714936115) URL https://www.pnas.org/doi/abs/10.1073/pnas.1714936115
  * [50] Krenn M, Erhard M and Zeilinger A 2020 Nat. Rev. Phys. 2 649–661 ISSN 2522-5820 URL https://doi.org/10.1038/s42254-020-0230-4
  * [51] Cervera-Lierta A, Krenn M and Aspuru-Guzik A 2022 Quantum 6 836 ISSN 2521-327X URL https://doi.org/10.22331/q-2022-10-13-836
  * [52] Ban Y, Echanobe J, Ding Y, Puebla R and Casanova J 2021 Quantum Sci. Technol. 6 045012 URL https://dx.doi.org/10.1088/2058-9565/ac16ed
  * [53] Chen Y, Ban Y, He R, Cui J M, Huang Y F, Li C F, Guo G C and Casanova J 2022 arXiv:2203.05849 URL https://arxiv.org/abs/2203.05849
  * [54] Torlai G, Mazzola G, Carleo G and Mezzacapo A 2020 Phys. Rev. Res. 2(2) 022060 URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.022060
  * [55] Von der Linden W D V and U V T 2014 Bayesian Probability Theory (Cambridge, UK: Cambridge University Press)
  * [56] Gelman A C J B and B R D 2004 Bayesian Data Analysis 2nd ed (New York, US: Chapman & Hall/CRC)
  * [57] Puebla R, Ban Y, Haase J, Plenio M, Paternostro M and Casanova J 2021 Phys. Rev. Appl. 16(2) 024044 URL https://link.aps.org/doi/10.1103/PhysRevApplied.16.024044
  * [58] Affleck I and Oshikawa M 1999 Phys. Rev. B 60(2) 1038–1056 URL https://link.aps.org/doi/10.1103/PhysRevB.60.1038
  * [59] Lindblad G 1976 Comm. Math. Phys. 48 119–130 ISSN 1432-0916 URL https://doi.org/10.1007/BF01608499
  * [60] Breuer H P and Petruccione F 2002 The theory of open quantum systems (Oxford, UK: Oxford University Press)


* (96) T. D. Lee and M. Nauenberg, _Degenerate systems and mass singularities_ , _Phys. Rev._ 133 (1964) B1549.
  * (97) I. Balitsky and A. Tarasov, _Power corrections to TMD factorization for Z-boson production_ , _JHEP_ 05 (2018) 150 [1712.09389].
  * (98) M. A. Ebert, I. Moult, I. W. Stewart, F. J. Tackmann, G. Vita and H. X. Zhu, _Subleading power rapidity divergences and power corrections for q T_, _JHEP_ 04 (2019) 123 [1812.08189].
  * (99) L. Cieri, C. Oleari and M. Rocco, _Higher-order power corrections in a transverse-momentum cut for colour-singlet production at NLO_ , _Eur. Phys. J. C_ 79 (2019) 852 [1906.09044].
  * (100) M. A. Ebert, J. K. L. Michel, I. W. Stewart and F. J. Tackmann, _Drell-Yan $q_{T}$ resummation of fiducial power corrections at N3LL_, _JHEP_ 04 (2021) 102 [2006.11382].
  * (101) V. Moos and A. Vladimirov, _Calculation of transverse momentum dependent distributions beyond the leading power_ , _JHEP_ 12 (2020) 145 [2008.01744].
  * (102) C. Oleari and M. Rocco, _Power corrections in a transverse-momentum cut for vector-boson production at NNLO: the $qg$-initiated real-virtual contribution_, _Eur. Phys. J. C_ 81 (2021) 183 [2012.10538].
  * (103) M. A. Ebert, A. Gao and I. W. Stewart, _Factorization for azimuthal asymmetries in SIDIS at next-to-leading power_ , _JHEP_ 06 (2022) 007 [2112.07680].
  * (104) A. Gehrmann-De Ridder and E. W. N. Glover, _A Complete O (alpha alpha-s) calculation of the photon + 1 jet rate in e+ e- annihilation_ , _Nucl. Phys. B_ 517 (1998) 269 [hep-ph/9707224].
  * (105) M. Ritzmann and W. J. Waalewijn, _Fragmentation in Jets at NNLO_ , _Phys. Rev. D_ 90 (2014) 054029 [1407.3272].
  * (106) J. M. Campbell and E. W. N. Glover, _Double unresolved approximations to multiparton scattering amplitudes_ , _Nucl. Phys. B_ 527 (1998) 264 [hep-ph/9710255].
  * (107) S. Catani and M. Grazzini, _Collinear factorization and splitting functions for next-to-next-to-leading order QCD calculations_ , _Phys. Lett. B_ 446 (1999) 143 [hep-ph/9810389].
  * (108) T. Gehrmann, A. von Manteuffel and T.-Z. Yang, _Renormalization of twist-two operators in QCD and its application to singlet splitting functions_ , in _16th DESY Workshop on Elementary Particle Physics: Loops and Legs in Quantum Field Theory 2022_ , 7, 2022, 2207.10108.