# A greedy reconstruction algorithm for the identification of spin
distribution

S. Buchwald111Universität Konstanz, Universitätsstraße 10, 78464 Konstanz,
Germany, G. Ciaramella222MOX-Laboratory, Dipartimento di Matematica,
Politecnico di Milano, Piazza Leonardo da Vinci 32, 20133 Milano, Italy, J.
Salomon333INRIA Paris, ANGE team, Building A, office 323, 2 rue Simone Iff,
75589 Paris, France, D. Sugny444Laboratoire Interdisciplinaire Carnot de
Bourgogne (ICB), UMR 6303 CNRS-Université Bourgogne-Franche Comté, 9 Av. A.
Savary, BP 47 870, F-21078 Dijon Cedex, France<EMAIL_ADDRESS>
###### Abstract

We propose a greedy reconstruction algorithm to find the probability
distribution of a parameter characterizing an inhomogeneous spin ensemble. The
identification is based on the application of a number of constant control
processes during a given time for which the final ensemble magnetization
vector is measured. From these experimental data, we show that the
identifiability of a piecewise constant approximation of the probability
distribution is related to the invertibility of a matrix which depends on the
different control protocols applied to the system. The algorithm aims to
design specific controls which ensure that this matrix is as far as possible
from a singular matrix. Numerical simulations reveal the efficiency of this
algorithm in different examples. A systematic comparison with respect to
random constant pulses is done.

##  1 Introduction

The identification of parameters that characterize the dynamics of a quantum
system is a fundamental prerequisite for controlling its evolution [1, 2, 3,
4, 5, 6, 7, 8] and is of practical interest for realizing specific tasks in
quantum technologies [9]. This aspect is crucial in open-loop configurations
for which the control protocols are designed without any experimental feedback
from the system during the control process [1, 2, 10, 11, 5]. In the context
of quantum systems, the problem of identifying unknown parameters (or
functions) has been explored by a large number of studies and for a variety of
applications ranging from molecular physics [12, 13, 14] and magnetic
resonance [15, 16, 17, 18] to quantum information science [19, 20, 21, 22, 23,
24, 25, 26, 27, 28, 29] and open quantum systems [30, 31, 32, 33]. Some
mathematical results have also been established in this direction [34, 35, 36,
37, 38, 39, 40, 41, 42]. On the basis of different measurement processes and
specific control protocols, the goal of these works is generally to estimate
the value of one or several parameters of the system Hamiltonian. When
controlling an ensemble of identical quantum systems, such a parameter may
vary in a given range due to experimental limitations or uncertainties. A key
example comes from the spatial inhomogeneities of the external control [43,
44, 2, 45, 46, 47]. In this case, all the systems are not subjected exactly to
the same control. This aspect has to be taken into account in the modeling of
the dynamics and in the computation of the control procedure. Robust control
protocols against such inaccuracies have been developed recently to solve this
experimental issue [2, 45, 48, 49, 50, 51, 52, 53, 54, 55]. However, the
variation range of the unknown parameter is not the only crucial quantity, the
probability distribution of this parameter (i.e. the number of systems for
each value of the parameter) may play also a major role. It is generally
assumed that this probability distribution is flat or has a simple Gaussian or
Lorentzian form. In these cases, the probability distribution can be quite
easily characterized. However, the problem of identifying probability
distributions becomes much more difficult if these have complex structures,
with e.g. several peaks, or if no information is known about them. It is
therefore essential to be able to identify with a great precision these
unknown probability distributions.

This paper aims at taking a step toward the answer to this open question by
developing a numerical algorithm, called a Greedy Reconstruction Algorithm
(GRA). By definition, an algorithm is said to be greedy if it takes the best
choice available at each iterative step. Greedy algorithms find generally a
sub-optimal solution, but in a computational time which may be very small
compared to the one of a global optimization procedure. Such algorithms have
been recently applied to the identification of quantum systems [56, 57] and we
propose to adapt them to the reconstruction of probability distribution. For
the sake of clarity, we focus in this study on a specific example, although
our algorithm applies to a large variety of systems. We consider the case of a
spin ensemble in Nuclear Magnetic Resonance (NMR) [44, 58, 59, 60, 61]
subjected to an inhomogeneous radio-frequency magnetic field whose range of
variation is known, but not its probability distribution. The probability
distribution is approximated by a piecewise constant function taking at most
$\displaystyle K$ values. The algorithm then designs a series of
$\displaystyle K$ controls for GRA (or less for the optimized version) which
are, in a second step, applied to the spin ensemble. The $\displaystyle K$
measured ensemble magnetization vectors at the final time are then used to
identify the probability distribution. More precisely, the identification
process is related to the invertibility of a matrix which depends on the
different controls. The aim of the algorithm is therefore to design specific
control protocols which ensure that this matrix is as far as possible from
being singular. The precision of the identification process can be understood
from the eigenvalues and eigenvectors of the matrix. In the examples analyzed
in this study, constant controls will be sufficient to find the probability
distribution with a very good accuracy, but time-dependent controls could also
be used. We show that the optimization procedure of the algorithm has a unique
solution and a good convexity structure leading to fast convergence. We point
out that the controls only depend on the model system and not on the spin
distribution or on the available data. The numerical efficiency of the
algorithms is shown on different illustrative examples, namely a double-peak
distribution and a step one. A systematic comparison with random constant
pulses is also done.

The paper is organized as follows. Section 2 describes the model system.
Section 3 is dedicated to the theoretical framework, while the algorithm is
presented in Sec. 4. The efficiency of GRA is numerically demonstrated in Sec.
5 on two standard examples, a double-peak and a step probability
distributions. Conclusion and prospective views are given in Sec. 6. A
mathematical description of the method and the development of an optimized
greedy algorithm are reported respectively in Appendices A and B. Additional
numerical results are presented in the supplemental material [62].

##  2 The model system

To illustrate our study, we consider a basic control problem in NMR, i.e. a
spin ensemble subjected to an inhomogeneous radio-frequency magnetic field
[43, 44, 45, 61, 46]. In a given rotating frame, we assume that all the spins
have the same resonance offset $\displaystyle\omega$. Each isochromat is
characterized by a Bloch vector
$\displaystyle\textbf{M}(\alpha)=[M_{x},M_{y},M_{z}]^{\intercal}$ whose
dynamics are governed by the following equations of motion:

$\displaystyle\begin{cases}\dot{M}_{x}=-\omega
M_{y}+(1+\alpha)\omega_{y}M_{z},\\\ \dot{M}_{y}=\omega
M_{x}-(1+\alpha)\omega_{x}M_{z},\\\
\dot{M}_{z}=(1+\alpha)\omega_{x}M_{y}-(1+\alpha)\omega_{y}M_{z},\end{cases}$

where the coordinates of the Bloch vector satisfy $\displaystyle
M_{x}^{2}+M_{y}^{2}+M_{z}^{2}=M_{0}^{2}$, with $\displaystyle M_{0}$ the
equilibrium magnetization. $\displaystyle\omega_{x}$ and
$\displaystyle\omega_{y}$ are time-dependent controls that correspond to the
components of the magnetic field along the $\displaystyle x$\- and the
$\displaystyle y$\- directions. In this study, we assume that these controls
are constant in time. We show in Sec. 5 that this hypothesis is sufficient for
the different examples to identify the probability distributions. The
parameter $\displaystyle\alpha$ is used to model the control field
inhomogeneities which are of the order of few percents in standard experiments
[60]. The controls $\displaystyle\frac{\omega_{x}}{2\pi}$ and
$\displaystyle\frac{\omega_{y}}{2\pi}$ are expressed in Hz. We consider a
typical field amplitude $\displaystyle\omega_{0}$ that can be fixed, for
instance, to $\displaystyle\omega_{0}=2\pi\times 100$ Hz. We introduce
normalized coordinates as follows:

$\displaystyle
u_{x}=2\pi\frac{\omega_{x}}{\omega_{0}};~{}u_{y}=2\pi\frac{\omega_{y}}{\omega_{0}};~{}t^{\prime}=\frac{\omega_{0}}{2\pi}t;~{}\Delta=2\pi\frac{\omega}{\omega_{0}};\textbf{X}=\frac{\textbf{M}}{M_{0}}.$

We omit the ′ in the time below to simplify the notations. We deduce that the
differential system can be expressed in normalized units as:

$\begin{cases}\dot{x}=-\Delta y+(1+\alpha)u_{y}z\\\ \dot{y}=\Delta
x-(1+\alpha)u_{x}z\\\ \dot{z}=(1+\alpha)u_{x}y-(1+\alpha)u_{y}z\end{cases}$
(1)

with $\displaystyle x^{2}+y^{2}+z^{2}=1$. The initial state of the dynamics
for each spin is the thermal equilibrium point, i.e. the north pole of the
Bloch sphere, $\displaystyle\textbf{X}_{0}=(0,0,1)^{\intercal}$. We neglect
the relaxation effect and we consider a control time of the order of 100 ms.
This corresponds to a normalized time $\displaystyle t_{f}$ of the order of
10. In the numerical simulations, we add the constraints
$\displaystyle|u_{x}|\leq u_{m}$ and $\displaystyle|u_{y}|\leq u_{m}$ where
$\displaystyle u_{m}$ is the maximum amplitude of each component. In NMR, only
the first two coordinates of the magnetization vector can be directly
measured. We do not have accessed to the $\displaystyle z$\- component due to
the strong magnetic field applied along this direction. We denote by
$\displaystyle\textbf{Y}(t)=(x(t),y(t))^{\intercal}$ the projection of the
Bloch vector onto the first two coordinates. We point out that this aspect is
not a limiting point for the application of the identification process.

##  3 Identification of spin distribution

We consider an ensemble of $\displaystyle N$ spins whose dynamics are governed
by Eq. (1). We assume that the control amplitudes $\displaystyle(u_{x},u_{y})$
belong to the admissible set
$\displaystyle\mathcal{U}=\\{(u_{x},u_{y})\in\mathbb{R}^{2}\mid|u_{x}|\leq
u_{m},|u_{y}|\leq u_{m}\\}$. The objective of the control procedure is to
identify the probability distribution of the parameter $\displaystyle\alpha$.
To simplify the recognition process, we assume that the ensemble of spins can
be decomposed into a set of $\displaystyle K$ subgroups with the same value of
the parameter $\displaystyle\alpha_{\ell}$, $\displaystyle 1\leq\ell\leq K$.
However, the discrete probability distribution $\displaystyle P_{\star}$ for
$\displaystyle\alpha$, namely the probability of each possible outcome
$\displaystyle\alpha_{\ell}$, or in other words the number of elements
$\displaystyle N_{\ell}$ of each subgroup, is unknown. By definition, we have
$\displaystyle P_{\star}(\ell)=\frac{N_{\ell}}{N}$ and
$\displaystyle\sum_{\ell=1}^{K}P_{\star}(\ell)=1$.

The projected solution onto the first two coordinates at time $\displaystyle
t_{f}$ of Eq. (1) is denoted by
$\displaystyle\textbf{Y}_{\textbf{u},\alpha}(t_{f})$ where the dependance on u
and $\displaystyle\alpha$ has been explicitly mentioned. The corresponding
experimental realization of this controlled dynamic leads to
$\displaystyle\textbf{Y}^{\textrm{exp}}_{\textbf{u}}(t_{f})=(x^{\textrm{exp}}_{\textbf{u}}(t_{f}),y^{\textrm{exp}}_{\textbf{u}}(t_{f}))^{\intercal}$,
where $\displaystyle\textbf{Y}^{\textrm{exp}}_{\textbf{u}}(t_{f})$ can be
viewed as the average at time $\displaystyle t_{f}$ of the experimental
measures of all the spins of the set subjected to the control u. The
coordinates $\displaystyle x^{\textrm{exp}}_{\textbf{u}}$ and $\displaystyle
y^{\textrm{exp}}_{\textbf{u}}$ are the ones of this measured magnetization
vector.

The relation between the theoretical description of the dynamical system to
the experimental outcome can be expressed as:

$\textbf{Y}^{\textrm{exp}}_{\textbf{u}}(t_{f})=\sum_{\ell=1}^{K}P_{\star}(\ell)\textbf{Y}_{\textbf{u},\alpha_{\ell}}(t_{f}),$
(2)

in which the two sides of the equation crucially depend on the control u. A
specific control protocol is not sufficient to identify the probability
distribution $\displaystyle P_{\star}$ which generally requires the
implementation of $\displaystyle K$ control processes with $\displaystyle K$
different controls denoted $\displaystyle\textbf{u}_{k}$, $\displaystyle
k=1,\cdots,K$. Note that in the optimized version of the GRA presented in
Appendix B, the number of controls can be different from $\displaystyle K$.

On the basis of the experimental outputs, a straightforward way to determine
$\displaystyle P_{\star}$ is to solve the following minimization problem:

$\min_{P\in\mathbb{P}}\sum_{k=1}^{K}\|\textbf{Y}^{\textrm{exp}}_{\textbf{u}_{k}}(t_{f})-\sum_{\ell=1}^{K}P(\ell)\textbf{Y}_{\textbf{u}_{k},\alpha_{\ell}}(t_{f})\|^{2},$
(3)

where $\displaystyle\mathbb{P}$ is the set of all the possible probability
distributions $\displaystyle P$ that satisfy $\displaystyle P(\ell)\geq 0$ for
$\displaystyle 1\leq\ell\leq K$ and $\displaystyle\sum_{\ell=1}^{K}P(\ell)=1$.
Mathematically, we point out that $\displaystyle\mathbb{P}$ is a convex and
closed set. $\displaystyle\|\cdot\|$ denotes the standard Euclidean vector
norm. Note that Eq. (3) can be rewritten as:

$\min_{P\in\mathbb{P}}\sum_{k=1}^{K}\|\sum_{\ell=1}^{K}(P_{\star}(\ell)-P(\ell))\textbf{Y}_{\textbf{u}_{k},\alpha_{\ell}}(t_{f})\|^{2}.$
(4)

At this point, it is clear that a key ingredient of the accuracy of the
identification process relies on the choice of the controls
$\displaystyle\textbf{u}_{k}$.

To clarify this problem, we introduce a set
$\displaystyle\\{\varphi_{j}\\}_{j=1}^{K}$ of linearly independent functions
$\displaystyle\varphi_{j}:\\{1,\cdots,K\\}\to\mathbb{R}$ such that
$\displaystyle\mathbb{P}\subset{\rm span}(\\{\varphi_{j}\\}_{j=1}^{K})$, where
span denotes the vector space generated by the functions. Expressing
respectively $\displaystyle P_{\star}$ and $\displaystyle P$ as $\displaystyle
P_{\star}(\ell)=\sum_{j=1}^{K}\beta_{\star,j}\varphi_{j}(\ell)$ and
$\displaystyle P(\ell)=\sum_{j=1}^{K}\beta_{j}\varphi_{j}(\ell)$, the
minimization problem (4) becomes:

$\min_{\beta\in\widehat{\mathbb{R}}^{K}}\sum_{k=1}^{K}\|\sum_{\ell,j=1}^{K}(\beta_{\star,j}-\beta_{j})\varphi_{j}(\ell)\textbf{Y}_{\textbf{u}_{k},\alpha_{\ell}}(t_{f})\|^{2},$
(5)

where the vector $\displaystyle\beta=(\beta_{j})_{j=1}^{K}$ is taken in
$\displaystyle\widehat{\mathbb{R}}^{K}$, a subset of
$\displaystyle\mathbb{R}^{K}$, so that $\displaystyle
P=\sum_{j}\beta_{j}\varphi_{j}$ is a probability distribution. Equation (5)
can be rewritten in a compact form as follows:

$\min_{\beta\in\widehat{\mathbb{R}}^{K}}\langle\beta_{\star}-\beta|W|\beta_{\star}-\beta\rangle,$
(6)

where $\displaystyle W$ is a symmetric and positive semi-definite
$\displaystyle K\times K$\- matrix whose elements are defined as:

$W_{\ell,j}=\sum_{k}\langle\gamma_{\ell}(\textbf{u}_{k})|\gamma_{j}(\textbf{u}_{k})\rangle$
(7)

with

$\displaystyle\gamma_{j}(\textbf{u}_{k})=\sum_{\ell}\varphi_{j}(\ell)\textbf{Y}_{\textbf{u}_{k},\alpha_{\ell}}(t_{f}).$

Since the set of vectors $\displaystyle\beta$ is a convex subset of
$\displaystyle\mathbb{R}^{K}$, we deduce that the problem is uniquely solvable
if the matrix $\displaystyle W$ is positive definite, i.e. if $\displaystyle
W$ has a non-zero determinant. In the case $\displaystyle W$ has a non-trivial
kernel, infinitely many solutions may exist which lead to wrong probability
distributions different from the experimental one $\displaystyle P_{\star}$.
We stress that the non-triviality of the kernel depends completely on the
choice of the controls $\displaystyle\textbf{u}_{k}$.

We show in this study that GRA allows us to design a set of controls
$\displaystyle\textbf{u}_{k}$ so that the matrix $\displaystyle W$ is positive
definite with a trivial kernel. The algorithm is composed of two steps, namely
an offline and an online steps. In the first stage, GRA computes the controls
$\displaystyle\textbf{u}_{k}$. In this phase, only the theoretical model is
needed without any experimental input. The derived controls are used in a
second step in which the different magnetization vectors are measured and the
minimization problem (3) is solved. Note that the controls are the same for
any probability distribution to identify and only depend on the model system
under study. Finally, we point out that, while in a first algorithm we
consider that all control pulses have the same duration $\displaystyle t_{f}$,
in a second version described in Sec. 4.2, the duration of each pulse is
considered as a variable to be optimized together with its amplitude. The
generality of GRA allows one to tackle this situation in a straightforward
manner.

##  4 A greedy reconstruction algorithm

We present in this section the GRA in its classical form, an optimized
extension called optimized GRA (OGRA) is described in Appendix B. For
pedagogical purposes, we have limited the mathematical derivation of the
algorithm to its strict minimum. The interested reader can find mathematical
details about the algorithms in [56] and [57] for the standard and optimized
GRA, respectively.

###  4.1 Optimizing the control amplitudes for a fixed control time

The GRA computes the controls $\displaystyle\textbf{u}_{k}$ by solving a
sequence of fitting-step and discriminatory-step problems, in which the goal
of the first step is to identify a nontrivial kernel of a sub-matrix of
$\displaystyle W$, while the second phase designs a new control which is aimed
to correct this discrepancy and to eliminate the identified non-trivial
kernel. The explicit formulation of the algorithm is given in terms of the
function $\displaystyle\textbf{h}^{(k)}$ defined by:

$\textbf{h}^{(k)}(\beta,\textbf{u})=\sum_{\ell=1}^{K}\sum_{j=1}^{k}\beta_{j}\varphi_{j}(\ell)\textbf{Y}_{\textbf{u},\alpha_{\ell}}(t_{f}),$
(8)

for any $\displaystyle\beta$ in $\displaystyle\mathbb{R}^{k}$. GRA is
described below. Some mathematical statements of the different steps of the
algorithm are described in Appendix A. Its numerical implementation is
presented and discussed in Sec. 5.

Greedy Reconstruction Algorithm (GRA): Given a set of $\displaystyle K$
linearly independent functions
$\displaystyle(\varphi_{1},\ldots,\varphi_{K})$.
Solve the initialization problem

$\max_{\textbf{u}\in\mathcal{U}}\|\textbf{h}^{(1)}(1,\textbf{u})\|^{2},$ (9)

that gives the control $\displaystyle\textbf{u}_{1}$, and set $\displaystyle
k=1$.
While $\displaystyle k\leq K-1$

  1. 1.

Fitting step: Find $\displaystyle(\beta^{k}_{j})_{j=1,\dots,k}$ that solves
the problem

$\min_{\beta\in\mathbb{R}^{k}}\sum_{m=1}^{k}\|\textbf{h}^{(K)}(\textbf{e}_{k+1},\textbf{u}_{m})-\textbf{h}^{(k)}(\beta,\textbf{u}_{m})\|^{2},$
(10)

where $\displaystyle\textbf{e}_{k+1}$ is the $\displaystyle(k+1)$-th canonical
vector in $\displaystyle\mathbb{R}^{K}$.

  2. 2.

Discriminatory step: Find $\displaystyle\textbf{u}_{k+1}$ that solves the
problem

$\max_{\textbf{u}\in\mathcal{U}}\|\textbf{h}^{(K)}(\textbf{e}_{k+1},\textbf{u})-\textbf{h}^{(k)}(\beta_{k},\textbf{u})\|^{2}.$
(11)

  3. 3.

Update $\displaystyle k+1\rightarrow k$.

End while

The basic principles of GRA can be detailed by its two first iterations for
$\displaystyle K=2$. Using

$\displaystyle\textbf{h}^{(1)}(1,\textbf{u})=\sum_{\ell=1}^{K}\varphi_{1}(\ell)\textbf{Y}_{\textbf{u},\alpha_{\ell}}(t_{f}),$

the initialization problem can be expressed as:

$\displaystyle\max_{\textbf{u}\in\mathcal{U}}\|W_{11}(\textbf{u})\|^{2}.$

We deduce that the goal of this step is to maximize the modulus of this
$\displaystyle W$\- matrix element, so as to be as far as possible from a zero
of $\displaystyle W$. We then consider the first step of the algorithm with
$\displaystyle k=1$. We omit below for clarity the dependence on u of
$\displaystyle W$. By definition, we have:

$\displaystyle\begin{cases}W_{11}=\|\sum_{\ell=1}^{K}\varphi_{1}(\ell)\textbf{Y}_{\textbf{u},\alpha_{\ell}}(t_{f})\|^{2}\\\
W_{22}=\|\sum_{\ell=1}^{K}\varphi_{2}(\ell)\textbf{Y}_{\textbf{u},\alpha_{\ell}}(t_{f})\|^{2}\\\
W_{12}=W_{21}=\sum_{\ell,\ell^{\prime}=1}^{K}\varphi_{1}(\ell)\varphi_{2}(\ell^{\prime})\langle\textbf{Y}_{\textbf{u},\alpha_{\ell}}(t_{f})|\textbf{Y}_{\textbf{u},\alpha_{\ell^{\prime}}}(t_{f})\rangle\end{cases}$

and we deduce that the quantity to minimize in the fitting step can be written
as:

$\displaystyle\|\textbf{h}^{(2)}(e_{2},\textbf{u}_{1})-\textbf{h}^{(1)}(\beta,\textbf{u}_{1})\|^{2}=W_{11}\beta^{2}-2W_{12}\beta+W_{22},$

where $\displaystyle\beta$ is here a real number. The minimum is reached for
$\displaystyle\beta_{1}=W_{11}^{-1}W_{12}$ where $\displaystyle W$ is computed
for the control $\displaystyle\textbf{u}_{1}$. This value can be associated to
a vector $\displaystyle(\beta_{1},-1)^{\intercal}$ of the kernel of the
following $\displaystyle 2\times 2$\- submatrix of W:

$\begin{pmatrix}W_{11}&W_{12}\\\ W_{12}&W_{22}\end{pmatrix}$ (12)

The fitting step of GRA can thus be interpreted as a systematic way to find a
basis of the kernel of larger and larger sub-matrices of $\displaystyle W$.
Setting $\displaystyle\beta$ to $\displaystyle\beta_{1}$, the discriminatory
step consists in adding a new control $\displaystyle\textbf{u}_{2}$ to correct
this singularity, i.e. in selecting this control such that the corresponding
quantity is as far as possible from a zero. Mathematically, it can be shown
that this procedure has always a solution and that the new matrix
$\displaystyle W^{(2)}$ has a non-trivial $\displaystyle 2\times 2$ sub-matrix
(12) (see Appendix A for details).

###  4.2 Optimizing amplitude controls and time horizon

Until now, we have considered a fixed control time $\displaystyle t_{f}$.
However, it is also possible to consider controls with different control
times, up to a fixed boundary $\displaystyle t_{f}^{\max}>0$. In this case, we
also maximize with respect to time, meaning that the initialization and
discriminatory step problems at iteration $\displaystyle k$ would change to

$\max\limits_{\begin{subarray}{c}\textbf{u}\in\mathcal{U},\\\
t_{f}\in[0,t_{f}^{\max}]\end{subarray}}\|\textbf{h}^{(1)}(1,\textbf{u};t_{f})\|^{2}$
(13)

and

$\max\limits_{\begin{subarray}{c}\textbf{u}\in\mathcal{U},\\\
t_{f}\in[0,t_{f}^{\max}]\end{subarray}}\|\textbf{h}^{(K)}(\boldsymbol{e}^{k+1},\textbf{u};t_{f})-\textbf{h}^{(k)}(\beta^{k},\textbf{u};t_{f})\|^{2},$
(14)

respectively. In Eq. (13) and (14), the function
$\displaystyle\textbf{h}^{(k)}$ is still defined as in (8), only with the
control time as an additional variable. Similarly, one can adapt the
corresponding problems in OGRA, the optimized version described in Appendix B.
We denote by GRAt and OGRAt the two resulting algorithms.

##  5 Numerical results

###  5.1 The case of a double peak distribution

As a first illustrative example, we investigate in this paragraph the
identification of a symmetric double peak probability distribution, displayed
in Fig. 1. Similar results have been achieved for other smooth distributions
with one or several peaks. Numerical details are described in the supplemental
material [62].

In the numerical simulations, we consider a control time $\displaystyle
t_{f}=t_{f}^{max}=16$. The amplitude $\displaystyle u_{m}$ is equal to 10. The
normalized offset resonance is set for all the spins to
$\displaystyle\frac{\pi}{10}$, i.e. to 30 Hz. We also assume that
$\displaystyle\alpha\in[-0.2,0.2]$ and $\displaystyle K=30$. The
$\displaystyle K$ discrete values $\displaystyle\alpha_{\ell}$ of
$\displaystyle\alpha$ are regularly spaced in the interval of variation of
$\displaystyle\alpha$, i.e.
$\displaystyle\alpha_{\ell}=-0.2+0.4\frac{\ell-1}{K-1}$. Since the control
protocols are constant in time, Eq. (1) is solved numerically by directly
evaluating the exponential matrix corresponding to the exact solution. All
optimization problems are solved by a BFGS descent-direction method. We also
mention that the exact number of uncoupled spins in the ensemble is not
relevant for all the computations, since we are only interested in their
probability distribution. However, we use a total number of $\displaystyle
10^{5}$ spins in the numerical simulations.

For GRA and GRAt, we consider a random and orthonormal basis
$\displaystyle\\{\varphi_{k}\\}_{k=1}^{30}$. Note that any basis of this space
can be used in the respective algorithms. For OGRA and OGRAt, we extend the
basis from GRA by 30 randomly chosen probability distributions
$\displaystyle\\{\varphi_{k}\\}_{k=31}^{60}$. The tolerance used in the OGRA
and OGRAt (see Appendix B) is set to be
$\displaystyle\textnormal{tol}=10^{-14}$. The controls generated by the
algorithms and corresponding to the numerical results discussed below are
described in the supplemental material [62]. To test whether it is even
necessary to run the algorithm or if the same results could be achieved with
other control protocols, we also consider two sets of 30 random and constant
controls. For the first and second cases, we use respectively completely
random constant values in the set $\displaystyle\mathcal{U}$, with a control
time $\displaystyle t_{f}$ or with different and random control times in the
interval $\displaystyle[0,t_{f}]$. We denote by RCC and RCCt the two sets of
controls.

The robustness of the different control functions is evaluated by considering
a 30-dimensional hypercube centered in the global minimum $\displaystyle
P_{\star}$ of our identification problem, with a radius of $\displaystyle
100\|P_{\star}\|$, and we repeat the minimization process for 100
initialization vectors randomly chosen in this hypercube. We then compute the
minimum norm difference
$\displaystyle\frac{\|P_{\star}-P_{f}\|}{\|P_{\star}\|}$ over all optimization
runs, where $\displaystyle P_{f}$ denotes the solution given by the
optimization algorithm. We obtain the results reported in Tab. 1.

Control set | GRA | GRAt | OGRA | OGRAt | RCC | RCCt
---|---|---|---|---|---|---
Min. error | 0.0045 | 0.0098 | 0.0005 | 0.0009 | 0.4685 | 0.0841
Table 1: Minimum relative norm error for different control sets.

As can be seen in Tab. 1, the errors of OGRA and OGRAt are ten times smaller
than the ones of GRA and GRAt, which themselves are respectively 10 and 100
times smaller than the errors of both sets of random controls. Similar results
have been achieved for other smooth distributions, which show the efficiency
of the two proposed algorithms. Figures 1 and 2 display respectively the true
distribution and the minimal solution for all control sets for fixed and
variable control times.

Figure 1: Plot as a function of $\displaystyle\alpha$ of the true distribution
(red crosses) and approximated solutions with minimum error, computed by the
different optimization algorithms for the identification problem using control
sets with a fixed control time. In particular, controls generated by OGRA
(green, vertical lines) and GRA (black, circles) and RCC (blue, squares) are
plotted. Dimensionless units are used. Figure 2: Same as Fig. 1, but for
variable control times.

We observe that the solutions computed with controls generated by any
algorithm match the true distribution. On the other hand, RCC completely
fails, showing a third peak in the middle, while RCCt can at least identify
the two peaks of the distribution.

###  5.2 The case of a step distribution.

As a second illustrative example, we consider a non-continuous step
distribution, displayed in Fig. 3, in which only spins with a positive
parameter $\displaystyle\alpha$ can be observed in the sample. We repeat the
numerical simulations of Sec. 5.1 and we obtain the results reported in Tab.
2.

Control set | GRA | GRAt | OGRA | OGRAt | RCC | RCCt
---|---|---|---|---|---|---
Min. error | 0.0295 | 0.0181 | 0.0018 | 0.0021 | 0.4204 | 0.1943
Table 2: Minimum relative norm error for different control sets.

As can be seen in Tab. 2, the difference in magnitude of errors is similar to
the one for the double peak distribution. These results are displayed in Fig.
3 and 4.

Figure 3: Same as Fig. 1 but for a step distribution. Figure 4: Same as Fig.
3, but for variable control times.

We observe that OGRA and OGRAt are still able to identify the true
distribution, while GRA and GRAt already show small discrepancies. RCC
completely fails again, but also RCCt shows major visible differences in the
upper part of the step distribution. Arguments based on the properties of the
matrix $\displaystyle W$ are given in Sec. 5.3 to explain such numerical
observations.

###  5.3 Eigenvalues and eigenvectors of $\displaystyle W$

We explain qualitatively in this section the numerical results observed in
Sec. 5.1 and 5.2 through the properties of the matrix $\displaystyle W$, i.e.
its eigenvalues and eigenvectors. We present the spectra of the matrix
$\displaystyle W$ for different control sets in Fig. 5 and 6. A very large
difference is observed between the eigenvalues associated with the optimized
controls and the random ones. Note that this observation is the same if the
control time is also optimized. This difference is quantitatively measured by
the condition number of $\displaystyle W$, i.e. the ratio between the largest
and the smallest eigenvalues, which is given in Tab 3. As could be expected,
these results show that, using random controls, the matrix $\displaystyle W$
can be close to being singular. In the example of Fig. 6, while most of the
eigenvalues are larger than 1, four of them are smaller than $\displaystyle
10^{-10}$. Hence, the matrix $\displaystyle W$ has a very bad condition
number. We stress the very good result achieved by OGRA for which all the
eigenvalues have almost the same value. This analysis may also explain the
difference between a smooth and a non-continuous probability distributions. As
a matter of fact, numerical results reveal that random controls have more
difficulty identifying non-smooth probability distribution as illustrated in
Sec. 5.2. This aspect can be understood from the behavior of the eigenvectors.
Indeed, we observe numerically that the modes with a large number of
oscillations correspond to the smallest eigvenvalues. Such modes have to be
used to reconstruct probability distributions with rapid and abrupt
variations. For random controls, these eigenvectors lead to large errors and
to wrong probability distributions.

Figure 5: Spectrum of the matrix $\displaystyle W$ defined in (7) for controls generated by OGRA (green, vertical lines), GRA (black, circles) and for RCC (blue, squares). Dimensionless units are used. Figure 6: Same as Fig. 5, but for variable control times. Control set | GRA | GRAt | OGRA | OGRAt | RCC | RCCt
---|---|---|---|---|---|---
cond($\displaystyle W$) | 4.9$\displaystyle\cdot 10^{3}$ | 6.9$\displaystyle\cdot 10^{7}$ | 19.55 | 16.6178 | 4.3$\displaystyle\cdot 10^{9}$ | 1.42$\displaystyle\cdot 10^{17}$
Table 3: Condition number of the matrix $\displaystyle W$ for different
control sets.

##  6 Conclusions

We have introduced in this work a Greedy Reconstruction Algorithm with an
application to spin dynamics. The algorithm provides a systematic way to
identify the probability distribution of a parameter of the Hamiltonian system
varying in a given range. The efficiency of the identification process has
been illustrated in the case of a spin ensemble subjected to an inhomogeneous
radio-frequency magnetic field. After having described some mathematical
properties of the algorithm, numerical simulations have revealed the
efficiency of GRA and its quite large basin of convergence. We have shown that
GRA is able to identify non-trivial probability distributions with several
peaks or with a step variation. An optimized version of this algorithm can be
derived to further improve the identification process. We have limited the
study to constant controls, but similar results can be achieved with time-
dependent pulses. A quantitative comparison with random constant controls have
highlighted the non-trivial recognition process realized by the algorithms.
The numerical observations can also be partly explained by the computation of
the eigenvalues and eigenvectors of the matrix $\displaystyle W$.

This analysis paves the way for further investigations in magnetic resonance.
An interesting direction is the study of the sensitivity of the algorithm to
experimental imperfections or to the presence of noise. It could be also used
to identify probability distribution of other parameters, such as the
resonance offset. These greedy algorithms could also be transferred to other
domains such as quantum optics and atomic and molecular physics. Finally, we
hope that our method will be used in relevant experimental applications in
magnetic resonance in a near future.

Acknowledgment
The work of the first author was supported by the DFG via the collaborative
research center SFB1432, Project-ID 425217212. This research has been
partially supported by the ANR project “QUACO” ANR-17-CE40-0007-01. This
project has received funding from the European Union Horizon 2020 research and
innovation program under Marie-Sklodowska-Curie Grant No. 765267 (QUSCO).

##  Appendix A Mathematical description of GRA

We give in this section some mathematical details about GRA. Straightforward
computations show that the different steps of GRA can be expressed in matrix
form as follows:

  * •

The initialization problem (9) is equivalent to:

$\displaystyle\max_{\textbf{u}\in\mathcal{U}}[W(\textbf{u})]_{1,1}$

  * •

The fitting-step problem (10) is equivalent to:

$\displaystyle\min_{\beta\in\mathbb{R}^{k}}\langle\beta|W^{k}_{[1:k,1:k]}|\beta\rangle-2\langle
W^{k}_{[1:k,k+1]}|\beta\rangle$

where $\displaystyle W^{k}=\sum_{m=1}^{k}W(u_{m})$. $\displaystyle
W^{k}_{[1:k,1:k]}$ and $\displaystyle W^{k}_{[1:k,k+1]}$ denote respectively
the $\displaystyle k\times k$ upper-left block of $\displaystyle W^{k}$ and a
column vector containing the first $\displaystyle k$ components of the
$\displaystyle k+1$-th column of $\displaystyle W^{k}$.

  * •

The discriminatory-step problem (11) is equivalent to:

$\displaystyle\max_{u\in\mathcal{U}}\langle\textbf{v}|[W(\textbf{u})]_{[1:k+1,1:k+1]}|\textbf{v}\rangle$

where $\displaystyle\textbf{v}=(\beta_{k}^{\intercal},-1)^{\intercal}$.

The different iterations of GRA can then be described as follows. At iteration
$\displaystyle k$, we assume that the sub-matrix $\displaystyle
W^{k}_{[1:k,1:k]}$ is positive definite, but $\displaystyle
W^{k}_{[1:k+1,1:k+1]}$ can have a non-trivial kernel. The idea is first to
identify the kernel of $\displaystyle W^{k}_{[1:k+1,1:k+1]}$ by solving (10)
and then to compute a new control $\displaystyle\textbf{u}_{k+1}$ such that
the new updated matrix $\displaystyle W^{k+1}=W^{k}+W(\textbf{u}_{k+1})$ has a
positive definite upper-left block $\displaystyle W^{k+1}_{[1:k+1,1:k+1]}$.
The convergence of the algorithm follows from this iterative process.

The following two technical lemmas describe the optimizations used in the two
steps of the algorithm. In particular, Lemma 1 shows that the fitting step
identifies the kernel of the matrix $\displaystyle W^{k}_{[1:k+1,1:k+1]}$.

###### Lemma 1

Assume that $\displaystyle W^{k}_{[1:k,1:k]}$ is positive definite and
$\displaystyle W^{k}_{[1:k+1,1:k+1]}$ has a non-trivial kernel. Then the
vector $\displaystyle{\bf v}=(\beta_{k}^{\intercal},-1)^{\intercal}$, where
$\displaystyle\beta_{k}$ is the solution to (10), is in the kernel of
$\displaystyle W^{k}_{[1:k+1,1:k+1]}$.

The second lemma is the basis of the discriminatory-step algorithm and shows
that this step corrects the rank deficiency of $\displaystyle
W^{k}_{[1:k+1,1:k+1]}$.

###### Lemma 2

Let $\displaystyle W^{k}_{[1:k,1:k]}$ be a positive definite matrix and
$\displaystyle\beta_{k}$ a solution of the fitting-step problem (10). Any
solution $\displaystyle\textbf{u}_{k+1}$ of (11) satisfies:

$\displaystyle\langle{\bf v}|W_{[1:k+1,1:k+1]}|{\bf v}\rangle>0$

for $\displaystyle k=0,1,\cdots,K-1$, where $\displaystyle{\bf
v}=(\beta_{k}^{\intercal},-1)^{\intercal}$.

The mathematical proofs of these results and a detailed numerical analysis of
the GRA is beyond the scope of this work and will be presented elsewhere.

##  Appendix B The optimized greedy algorithm

We discuss in this paragraph the optimized version of GRA. It can been shown
numerically that the behavior and the efficiency of GRA is strongly affected
by the choice of the elements $\displaystyle\varphi_{k}$ and their ordering.
GRA is essentially a sweep over the set
$\displaystyle(\varphi_{k})_{k=1}^{K}$. However, a wrong choice of the
elements $\displaystyle\varphi_{k}$ and their ordering can lead to the
stagnation of the algorithm and to the computation of many useless control
functions. Note that the stagnation of the algorithm can be measured in terms
of rank corrections, i.e. if for consecutive iterations the rank of
$\displaystyle W$ does not increases. These reasons are at the origin of an
optimized algorithm [57]. OGRA takes as input a set $\displaystyle\Phi$,
possibly larger than $\displaystyle(\varphi_{k})_{k=1}^{K}$ with linearly
dependent elements, and returns as output not only a set of
$\displaystyle\widetilde{K}$ control functions, but also a set of linearly
independent functions
$\displaystyle(\widehat{\varphi}_{k})_{k=1}^{\widehat{K}}$. The integers
$\displaystyle\widetilde{K}$ and $\displaystyle\widehat{K}$ are not
necessarily equal and may be smaller than $\displaystyle K$ (in contrast to
GRA). The extension of the OGR method of [57] to the distribution
reconstruction problem is detailed below, where we use the map
$\displaystyle\textbf{h}_{\mathcal{S}}$ defined as

$\displaystyle\textbf{h}_{\mathcal{S}}(\beta,\textbf{u})=\sum_{\ell=1}^{K}\sum_{j=1}^{\textrm{card}[\mathcal{S}]}\beta_{j}\varphi_{j}(\ell)\textbf{Y}(\textbf{u},\alpha_{\ell}),$

where $\displaystyle\mathcal{S}=(\varphi_{1},\cdots,\varphi_{k})$. Note that
for the fitting-step problem, we do not have any constraint for the choice of
coefficients $\displaystyle\beta$. This is due to the fact that, during the
algorithm, we are not trying to reconstruct a distribution but to make the
respective sub-matrix positive definite.

Optimized Greedy Reconstruction Algorithm (OGRA): Given a set of
$\displaystyle K_{+}\geq K$ linearly independent functions
$\displaystyle(\varphi_{1},\ldots,\varphi_{K_{+}})$ and a tolerance
$\displaystyle\textrm{tol}>0$.
Solve the initialization problem

$\max_{n\in\\{1,\cdots,K_{+}\\}}\max_{\textbf{u}\in\mathcal{U}}\|\textbf{h}_{\varphi_{n}}^{(1)}(1,\textbf{u})\|^{2},$
(15)

which gives the control $\displaystyle\textbf{u}_{1}$, and the control
$\displaystyle\ell_{1}$. Set $\displaystyle k=1$ and
$\displaystyle\mathcal{S}=\\{\varphi_{\ell_{1}}\\}$,
$\displaystyle\widetilde{K}=K_{+}$, and update
$\displaystyle\Phi=\Phi\setminus\\{\varphi_{\ell_{1}}\\}$. The algorithm is
stopped if
$\displaystyle\|\textbf{h}_{\mathcal{S}}^{(1)}(1,\textbf{u})\|^{2}<\textrm{tol}$.
While $\displaystyle k\leq K-1$ do

  1. 1.

Remove elements from $\displaystyle\Phi$ that are linearly dependent on the
ones in $\displaystyle\mathcal{S}$. Shift the indices of the remaining
elements in $\displaystyle\Phi$. Update
$\displaystyle\textrm{card}[\Phi]\rightarrow\widetilde{K}$.

  2. 2.

for $\displaystyle\ell=1,\cdots,\widetilde{K}$ do
Fitting step: Find $\displaystyle(\beta^{\ell}_{j})_{j=1,\dots,k}$ that solve
the problem

$\min_{\beta\in\mathbb{R}^{k}}\sum_{m=1}^{k}\|\textbf{h}_{\varphi_{\ell}}(1,\textbf{u}_{m})-\textbf{h}_{\mathcal{S}}(\beta,\textbf{u}_{m})\|^{2},$
(16)

end for

  3. 3.

Discriminatory step: Find $\displaystyle\textbf{u}_{k+1}$ and
$\displaystyle\ell_{k+1}$ that solve the problem

$\max_{\ell\in\\{1,\cdots,K_{+}\\}}\max_{\textbf{u}\in\mathcal{U}}\|\textbf{h}_{\varphi_{\ell}}(1,\textbf{u})-\textbf{h}_{\mathcal{S}}(\beta^{\ell},\textbf{u})\|^{2}.$
(17)

If
$\displaystyle\|\textbf{h}_{\varphi_{\ell}}(1,\textbf{u}^{k+1})-\textbf{h}_{\mathcal{S}}(\beta^{\ell_{k+1}},\textbf{u}_{k+1})\|^{2}<\textrm{tol}$
then stop and return $\displaystyle\mathcal{S}$ and the computed
$\displaystyle(\textbf{u})_{m=1}^{k}$.

  4. 4.

Orthogonalize the function $\displaystyle\varphi_{\ell_{k+1}}$ with respect to
$\displaystyle\mathcal{S}$ and update
$\displaystyle\mathcal{S}\cup\\{\varphi_{\ell_{k+1}}\\}\rightarrow\mathcal{S}$,
$\displaystyle\Phi\setminus\\{\varphi_{\ell_{k+1}}\\}\rightarrow\Phi$ and
$\displaystyle k+1\rightarrow k$.

End while

## References

  * [1] C. Brif, R. Chakrabarti, and R. Rabitz, New J. Phys. 12, 075008 (2010).
  * [2] S. J. Glaser, U. Boscain, T. Calarco, C. P. Koch, W. Köckenberger, R. Kosloff, I. Kuprov, B. Luy, S. Schirmer, T. Schulte-Herbrüggen, D. Sugny, and F. K. Wilhelm, Eur. Phys. J. D 69, 279 (2015).
  * [3] C. P. Koch, M. Lemeshko, and D. Sugny, Rev. Mod. Phys. 91, 035005 (2019).
  * [4] C. P. Koch, J. Phys. Condens. Matter 28, 213001 (2016).
  * [5] A. Borzì, G. Ciaramella, and M. Sprengel,  _Formulation and Numerical Solution of Quantum Control Problems_ (SIAM, Philadelphia, PA, 2017).
  * [6] Y. Maday, J. Salomon, and G. Turinici, Numer. Math. 103, 323 (2006).
  * [7] M. Sprengel, G. Ciaramella, and A. Borzì, Journal of Dynamical and Control Systems 24, 657 (2018).
  * [8] U. Boscain, M. Sigalotti, and D. Sugny, to be published in PRX Quantum (2021).
  * [9] A. Acín, I. Bloch, H. Buhrman, T. Calarco, C. Eichler, J. Eisert, D. Esteve, N. Gisin, S. J. Glaser, F. Jelezko, S. Kuhr, M. Lewenstein, M. F. Riedel, P. O. Schmidt, R. Thew, A. Wallraff, I. Walmsley, and F. K. Wilhelm, New J. Phys. 20, 080201 (2018).
  * [10] D. D’Alessandro,  _Introduction to Quantum Control and Dynamics_ (Chapman & Hall/CRC, Boca Raton, 2007).
  * [11] A. E. Bryson and Y.-C. Ho,  _Applied Optimal Control: Optimization, Estimation, and Control_ (Hemisphere, Washington, DC, 1975).
  * [12] J. M. Geremia and H. Rabitz, The Journal of Chemical Physics 118, 5369 (2003).
  * [13] A. Feng, H. Rabitz, G. Turinici, and C. Le Bris, J. Phys. Chem. A 110, 7755 (2006).
  * [14] J. M. Geremia and H. Rabitz, Phys. Rev. Lett. 89, 263902 (2002).
  * [15] D. Ma, V. Gulani, and N. Seiberlich, Nature 495, 187 (2013).
  * [16] Q. Ansel, M. Tesch, S. J. Glaser, and D. Sugny, Phys. Rev. A 96, 053419 (2017).
  * [17] E. Y. Pierre, D. Ma, Y. Chen, C. Badve, and M. A. Griswold, Magnetic Resonance in Medicine 75, 2481 (2016).
  * [18] G. Ciaramella, A. Borzì, G. Dirr, and D. Wachsmuth, SIAM Journal on Scientific Computing 37, A319 (2015).
  * [19] S. G. Schirmer, A. Kolli, and D. K. L. Oi, Phys. Rev. A 69, 050306 (2004).
  * [20] S. G. Schirmer and D. K. L. Oi, Phys. Rev. A 80, 022333 (2009).
  * [21] S. G. Schirmer and F. C. Langbein, Phys. Rev. A 91, 022125 (2015).
  * [22] H. Yuan and C. F. Fung, npj Quantum Inf. 3, 14 (2017).
  * [23] J. Liu and H. Yuan, Phys. Rev. A 96, 012117 (2017).
  * [24] N. Wittler, F. Roy, K. Pack, M. Werninghaus, A. S. Roy, D. J. Egger, S. Filipp, F. K. Wilhelm, and S. Machnes, Phys. Rev. Applied 15, 034080 (2021).
  * [25] J. Zhang and M. Sarovar, Phys. Rev. Lett. 113, 080401 (2014).
  * [26] A. Sone and P. Cappellaro, Phys. Rev. A 95, 022335 (2017).
  * [27] J. Kiukas, K. Yuasa, and D. Burgarth, Phys. Rev. A 95, 052132 (2017).
  * [28] D. Burgarth and A. Ajoy, Phys. Rev. Lett. 119, 030402 (2017).
  * [29] S.-T. Wang, D.-L. Deng, and L.-M. Duan, New Journal of Physics 17, 093017 (2015).
  * [30] S. Xue, R. Wu, S. Ma, D. Li, and M. Jiang, Phys. Rev. A 103, 022604 (2021).
  * [31] S. Xue, L. Tan, R. Wu, M. Jiang, and I. R. Petersen, Phys. Rev. A 102, 042227 (2020).
  * [32] J. Zhang and M. Sarovar, Phys. Rev. A 91, 052121 (2015).
  * [33] R. S. Bennink and P. Lougovski, New Journal of Physics 21, 083013 (2019).
  * [34] L. Ljung, Annual Reviews in Control 34, 1 (2010).
  * [35] Le Bris, Claude, Mirrahimi, Mazyar, Rabitz, Herschel, and Turinici, Gabriel, ESAIM: COCV 13, 378 (2007).
  * [36] O. F. Alis, H. Rabiz, M. Q. Rosenthal, C. Phan, and M. Pence, J. Math. Chem. 35, 65 (2004).
  * [37] M. Ndong, J. Salomon, and D. Sugny, J. Phys. A: Math. Theor. 47, 265302 (2014).
  * [38] C. R. Rojas, J. S. Welsh, G. C. Goodwin, and A. Feuer, Automatica 43, 993 (2007).
  * [39] L. Baudouin and A. Mercado, Appl. Anal. 87, 1145 (2008).
  * [40] S. Bonnabel, M. Mirrahimi, and P. Rouchon, Automatica 45, 1144 (2009).
  * [41] Y. Fu and G. Turinici, Appl. Anal. 23, 1129 (2017).
  * [42] Y. Wang, D. Dong, B. Qi, J. Zhang, I. R. Petersen, and H. Yonezawa, IEEE Trans. Autom. Control 63, 1388 (2018).
  * [43] J. Li and N. Khaneja, Phys. Rev. A 73, 030302 (2006).
  * [44] J. S. Li and N. Khaneja, IEEE Trans. Auto. Control 54, 528 (2009).
  * [45] K. Kobzar, T. E. Skinner, N. Khaneja, S. J. Glaser, and B. Luy, J. Magn. Reson. 194, 58 (2008).
  * [46] M. Lapert, Y. Zhang, M. A. Janich, S. J. Glaser, and D. Sugny, Sci. Rep. 2, 589 (2012).
  * [47] G. Turinici, Phys. Rev. A 100, 053403 (2019).
  * [48] L. Van Damme, Q. Ansel, S. J. Glaser, and D. Sugny, Phys. Rev. A 95, 063403 (2017).
  * [49] A. Ruschhaupt, X. Chen, D. Alonso, and J. G. Muga, New Journal of Physics 14, 093040 (2012).
  * [50] L. Van-Damme, D. Schraft, G. T. Genov, D. Sugny, T. Halfmann, and S. Guérin, Phys. Rev. A 96, 022309 (2017).
  * [51] D. Daems, A. Ruschhaupt, D. Sugny, and S. Guérin, Phys. Rev. Lett. 111, 050404 (2013).
  * [52] G. T. Genov, D. Schraft, T. Halfmann, and N. V. Vitanov, Phys. Rev. Lett. 113, 043001 (2014).
  * [53] D. Buterakos, S. Das Sarma, and E. Barnes, PRX Quantum 2, 010341 (2021).
  * [54] J. Zeng and E. Barnes, Phys. Rev. A 98, 012301 (2018).
  * [55] J. Zeng, C. H. Yang, A. S. Dzurak, and E. Barnes, Phys. Rev. A 99, 052321 (2019).
  * [56] Y. Maday and J. Salomon, in  _Proceedings of the 48th IEEE Conference on Decision and Control, 2009, Held jointly whit the 28th Chinese Control Conference (CDC/CCC 2009)_ , IEEE Conference on Decision and Control (2009) pp. 375–379.
  * [57] S. Buchwald, G. Ciaramella, and J. Samonon, to appear in SIAM J. Control Optim. (2021).
  * [58] S. Conolly, D. Nishimura, and A. Macovski, IEEE Trans. Med. Imaging 5, 106 (1986).
  * [59] M. Lapert, Y. Zhang, M. Braun, S. J. Glaser, and D. Sugny, Phys. Rev. Lett. 104, 083001 (2010).
  * [60] M. Levitt,  _Spin Dynamics: Basics of Nuclear Magnetic Resonance_ (Wiley, Chichester, UK, 2008).
  * [61] T. E. Skinner, T. O. Reiss, B. Luy, N. Khaneja, and S. J. Glaser, J. Magn. Reson. 172, 17 (2005).
  * [62] Numerical results are available upon request to the corresponding author


# Quantum description of reality is empirically incomplete

Anubhav Chaturvedi<EMAIL_ADDRESS>Institute of
Theoretical Physics and Astrophysics, National Quantum Information Center,
Faculty of Mathematics, Physics and Informatics, University of Gdansk, 80-952
Gdańsk, Poland  International Centre for Theory of Quantum Technologies
(ICTQT), University of Gdansk, 80-308 Gdańsk, Poland  Marcin Pawłowski
Institute of Theoretical Physics and Astrophysics, National Quantum
Information Center, Faculty of Mathematics, Physics and Informatics,
University of Gdansk, 80-952 Gdańsk, Poland  International Centre for Theory
of Quantum Technologies (ICTQT), University of Gdansk, 80-308 Gdańsk, Poland
Debashis Saha  S. N. Bose National Centre for Basic Sciences, Block JD, Sector
III, Salt Lake, Kolkata 700098, India

###### Abstract

Empirical falsifiability of the predictions of physical theories is the
cornerstone of the scientific method. Physical theories attribute empirically
falsifiable operational properties to _sets_ of physical preparations. A
theory is said to be _empirically complete_ if such properties allow for a
_not_ fine-tuned realist explanation, as properties of underlying probability
distributions over states of reality. Such theories satisfy a family of
_equalities_ among fundamental operational properties, characterized
exclusively by the number of preparations. Quantum preparations deviate from
these equalities, and the maximal quantum deviation increases with the number
of preparations. These deviations not only signify the _incompleteness_ of the
operational quantum formalism, but they simultaneously imply quantum over
classical _advantage_ in suitably constrained one-way communication tasks,
highlighting the delicate interplay between the two.

††preprint: APS/123-QED

Quantum theory is all set to fuel the key technological advances of the 21st
century. However, even after almost a century since its conception, there is
no consensus about the precise sense in which the structure of reality that
quantum theory posits conflicts with classical worldviews. Addressing such
questions necessitates formal notions of classicality, and the realist
(ontological) framework provides a vital ground for such notions [1, 2, 3, 4].

Typically, these notions of classicality ascribe certain operational phenomena
a _not_ fine-tuned realist basis by requiring these phenomena to hold intact
at the level of potentially inaccessible underlying reality [5]. These
phenomena double as empirically falsifiable operational prerequisites for
tests of the notions of classicality. For instance, Bell’s local causality [1,
6] ascribes to non-signalling correlations a parameter-independent realist
explanation 111Bell’s local-causality entails two distinct realist
assumptions. While _parameter independence_ forms the not fined tuned realist
basis for the operational no-signaling condition, _outcome independence_ is an
auxiliary assumption of purely causal nature., and the no-signaling condition
ensured by sufficient spatial separation forms the operational prerequisite of
Bell tests. Generalised noncontextuality attributes identical realist
counterparts to operationally equivalent physical entities, and the
operational indistinguishability of these entities forms the prerequisite of
experimental tests of noncontextuality [3]. The recently introduced notion of
classicality, bounded ontological distinctness (BOD) explains the
distinguishability of operational physical entities by the distinctness of
their realist counterparts, and $p$-distinguishability of the physical
entities forms a prerequisite of the experimental tests of BOD [4].

On the other hand, the realist notions of classicality yield empirically
falsifiable operational consequences, typically in the form of statistical
inequalities. The quantum violation of these inequalities not only highlights
the necessity of realist fine-tuning, discarding a large class of realist
explanations, but also powers quantum advantage in a plethora of
computational, communication and information processing tasks [7, 8, 4, 9].
Therefore, empirically falsifiable phenomena feature as the operational
prerequisites, as well as the operational consequences of the realist notions
of classicality. Moreover, to distill the spectrum of non-classical
operational predictions of quantum theory, all empirically falsifiable
operational phenomena warrant a _not_ fine-tuned realist basis. Stemming from
this unifying perception, in this Letter, we introduce a realist notion of
classicality, termed _empirical completeness_.

Operational theories attribute _empirically falsifiable operational
properties_ to sets of preparations. Here, we consider an operationally
relevant class of such properties, namely, the maximum success metrics of one-
way communication tasks. If these properties allow for a _not_ fine-tuned
realist explanation, i.e., if there exists corresponding sets of probability
distributions over states of reality (referred to as epistemic states) which
exactly and exclusively explain the operational properties, the theory or the
fragment thereof is said to be empirically complete.

In contrast to other notions of classicality, we demonstrate that empirically
complete theories satisfy a family of _equalities_ among elemental operational
properties of set of preparations. These equalities hold irrespective of the
particulars of the preparations and are characterised solely by their number.
The corresponding properties of sets of quantum preparations deviate from
these equalities, signifying the incompleteness of the theory. Moreover, as
classical communication is empirically complete, these deviations imply
quantum over classical advantage in suitably constrained one-way communication
tasks. The fact that classical preparations adhere to statistical _equalities_
implies that quantum advantage can be obtained with very inefficient
detectors, making the experimental demonstrations much more feasible [10].
Finally, employing state of the art semidefinite programming (SDP) techniques
we characterise the quantum deviations for the case of three preparations, and
provide evidence of increasing maximal quantum deviation with the number of
preparations.

Let us begin by revisiting the requisite preliminaries. Operational physical
theories such as quantum theory serve a two fold purpose, (i) they _prescribe_
mathematical objects to experimental procedures, and (ii) _predict_ the
consequent observations. Here we consider simple experiments entailing a
preparation $P\in\mathcal{P}_{\mathcal{O}}$ of a physical system, followed by
a $K$-outcome measurement $M\in\mathcal{M}_{\mathcal{O}}$, where
$\mathcal{P}_{\mathcal{O}}$ and $\mathcal{M}_{\mathcal{O}}$ are sets of
preparations and measurements featuring in an operational theory. The theory
yields predictions of the form $p(k|P,M)$ specifying the probability of
observing an outcome $k\in[K]$ (where $[a]\equiv\\{1,\ldots,a\\}$) as a result
of performing a measurement $M$ on a preparation $P$.

In quantum theory preparations are described by density operators
$\rho_{P}\in\\{\rho\in\mathcal{B}_{+}(\mathcal{H})|\Tr\rho=1\\}$, and
measurements are defined by POVMs
$\\{M_{k}\in\mathcal{B}_{+}(\mathcal{H})\\}_{k}$ such that
$\sum_{k}M_{k}=\mathbb{I}$, where $\mathcal{B}_{+}(\mathcal{H})$ is the set of
bounded positive semidefinite operators on some Hilbert space $\mathcal{H}$
which is assumed to be finite dimensional, and $\mathbb{I}$ is the identity
operator on $\mathcal{H}$. The predictions of quantum theory are specified by
the Born rule $p(k|P,M)=\Tr(\rho_{P}M_{k})$. Yet another relevant class of
operational models is that of $d$-levelled classical communication. In such
theories an operational preparation $P$ corresponds to an _encoding scheme_
entailing a probability distribution $e_{P}$ over a $d$-levelled classical
message $\omega\in[d]$, and a $K$-outcome measurement $M$ is described by a
_decoding scheme_ specifying the probability $d_{M}(k|\omega)$ of producing an
outcome $k\in[K]$ upon receiving the message $\omega$. The predictions of such
models are given by $p(k|P,M)=\sum_{\omega}e_{P}(\omega)d_{M}(k|\omega)$.

Owing to the ubiquitous distinction between preparation and measurement
procedures, operational theories implicitly assign empirically falsifiable
properties to sets of physical preparations. In this Letter, we consider an
operationally relevant subclass of such properties defined as follows,

######  Definition 1 (Operational properties).

Associated to a given set of $n$ preparations
$\bm{P}\equiv\\{P_{x}\in\mathcal{P}_{\mathcal{O}}\\}^{n}_{x=1}$, an
empirically falsifiable operational property has the generic form,

$S^{(\mathcal{O})}_{n}(\bm{P})=\max_{M\in\mathcal{M}_{\mathcal{O}}}\\{\sum_{x,k}c^{x}_{k}p(k|P_{x},M)\\},$
(1)

where $\\{c^{x}_{k}\in\mathbb{R}\\}_{x,k}$ are real coefficients and the
maximization is over the set of all measurements $\mathcal{M}_{\mathcal{O}}$
allowed in the operational theory.

The operational properties (1) are closely linked to one-way communication
tasks, wherein the sender encodes the classical input $x\in[n]$ onto a
physical preparation $P_{x}\in\bm{P}$, and receiver decodes the message by
performing some measurement $M\in\mathcal{M}_{\mathcal{O}}$ to produce an
outcome $k$. We gauge the performance of these resources using a success
metric $S^{(\mathcal{O})}_{n}(\bm{P},M)=\sum_{x,k}c^{x}_{k}p(k|P_{x},M)$. For
a given set of preparations $\bm{P}$, each operational property (1)
corresponds to the maximal attainable success metric, i.e.,
$S^{(\mathcal{O})}_{n}(\bm{P})=\max_{M\in\mathcal{M}_{\mathcal{O}}}\\{S^{(\mathcal{O})}_{n}(\bm{P},M)\\}$.
The maximization relieves $S^{(\mathcal{O})}_{n}(\bm{P})$ of its dependence on
the measurements featuring in an operational theory deeming it to be an
exclusive operational property of the given set of preparations $\bm{P}$.

For a given set of preparations $\bm{P}$, if one can experimentally
demonstrate higher success than $S^{(\mathcal{O})}_{n}(\bm{P})$ in the
associated communication task, then the operational theory and its
prescriptions are falsified. For a given set of quantum preparations
$\bm{\rho}\equiv\\{\rho_{x}\\}^{n}_{x=1}$, evaluating an operational property
$S^{(\mathcal{Q})}_{n}(\bm{\rho})=\max_{M\in\mathcal{M}_{\mathcal{Q}}}\\{\sum_{x,k}c^{x}_{k}\tr(\rho_{x}M_{k})\\}$,
invariably constitutes a semidefinite program, which owing its to strong
duality yields the precisely the maximum success metric along with the optimal
POVM $\\{M_{k}\\}$ [11]. Whereas evaluating an operational propriety
$S^{(\mathcal{C})}_{n}(\bm{e})$ for a given set of encoding schemes
$\bm{e}\equiv\\{e_{x}\\}^{n}_{x=1}$ constitutes a linear program.

The realist framework seeks to explain the predictions of an operational
theory, whilst assuming the existence of observer independent attributes
associated with physical systems. The complete specification of each attribute
is referred to as the state of reality $\lambda$ of a physical system. A
realist model for a prepare and measure fragment of an operational theory
consists of the following three elements: (i) a measurable space, $\Lambda$,
known as the _realist state space_ , (ii) a probability measure $\mu_{P}$ on
$\Lambda$ describing the _epistemic state_ of the system for each preparation
$P\in\mathcal{P}_{\mathcal{O}}$, and (iii) for every state of reality
$\lambda\in\Lambda$ and measurement $M\in\mathcal{M}_{\mathcal{O}}$, a
probability distribution $\xi_{M}(\cdot|\lambda)$ over the possible outcomes
of $M$, referred to as a _response scheme_. A realist model explains the
operational predictions if
$p(k|P,M)=\int_{\Lambda}d\lambda\mu_{P}(\lambda)\xi_{M}(k|\lambda)$.

A realist explanation of an operational property (1) has the form
$S^{(\mathcal{O})}_{n}(\bm{P})=\max_{\\{M\in\mathcal{M}_{\mathcal{O}}\\}}\\{\sum_{x,k}c^{x}_{k}\int_{\Lambda}d\lambda\mu_{x}(\lambda)\xi_{M}(k|\lambda)\\}$,
where the maximization is over the possibly fine-tuned set of operationally
accessible response schemes. To distill the necessity of such fine-tuning, we
define _not_ fine-tuned exclusive properties of sets of epistemic states as,

######  Definition 2 (_Not_ fine-tuned realist properties).

Associated to a given set of $n$ epistemic states
$\bm{\mu}\equiv\\{\mu_{x}\\}^{n}_{x=1}$ a _not_ fine-tuned realist property
has the generic form,

$S^{(\Lambda)}_{n}(\bm{\mu})=\max_{\\{\xi(k|\lambda)\\}}\\{\sum_{x,k}c^{x}_{k}\int_{\Lambda}d\lambda\mu_{x}(\lambda)\xi(k|\lambda)\\},$
(2)

where $\\{c^{x}_{k}\in\mathbb{R}\\}_{x,k}$ are real coefficients and the
maximization is over the set of all valid response schemes which satisfy
_positivity_ ($\forall\lambda\in\Lambda,k\in[K]:\ \xi(k|\lambda)\geqslant 0$)
and _completeness_ ($\forall\lambda\in\Lambda:\sum_{k}\xi(k|\lambda)=1$).

As the set of response schemes constrained by only by positivity and
completeness forms a convex polytope with deterministic response functions as
extremal points, we can solve the maximization in (2) by selecting the
response functions that for each $\lambda$, yield the outcome $k$ which
maximises the function $\sum_{x}c^{x}_{k}\mu_{x}(\lambda)$ such that,

$S^{(\Lambda)}_{n}(\bm{\mu})=\int_{\Lambda}d\lambda\max_{k}\\{\sum_{x}c^{x}_{k}\mu_{x}(\lambda)\\}.$
(3)

This expression further substantiates the fact that the maximization over
response schemes relieves $S^{(\Lambda)}_{n}(\bm{\mu})$ from its dependence on
response schemes, deeming it to be an exclusive property of the set of
epistemic states $\bm{\mu}$.

Now we have all the necessary ingredients to formally introduce _empirical
completeness_ as a characterising feature of operation theories,

######  Definition 3 (Empirically complete theories).

An operational theory or a fragment thereof is said to be empirically complete
($\mathcal{EC}$) if for all sets of preparations
$\bm{P}\equiv\\{P_{x}\in\mathcal{P}_{\mathcal{O}}\\}^{n}_{x=1}$, and all
associated empirically falsifiable operational properties
$S^{(\mathcal{O})}_{n}(\bm{P})$ (1), there exists underlying sets of epistemic
states $\bm{u}\equiv\\{\mu_{x}\\}^{n}_{x=1}$ with _not_ fine-tuned realist
properties $S^{(\Lambda)}_{n}(\bm{\mu})$ (2) such that,

$S^{(\Lambda)}_{n}(\bm{\mu})=S^{(\mathcal{O})}_{n}(\bm{P}).$ (4)

As an example, $d$-levelled classical communication models are empirically
complete, owing to the fact that the $d$-levelled classical message $\omega$
forms an operationally accessible sufficient statistic for the underlying
state of reality $\lambda$. This implies that for all empirically falsifiable
operational properties $S^{\mathcal{C}}_{n}(\bm{e})$ associated with a set of
encoding schemes $\bm{e}\equiv\\{e_{x}\\}^{n}_{x=1}$, we can obtain not fine
tuned realist explanations by taking the message itself to be the state of
reality, i.e., $\lambda=\omega$, and the encoding schemes to be the epistemic
states, i.e., $\bm{e}\equiv\bm{\mu}$, such that
$S^{(\Lambda)}_{n}(\bm{\mu})=S^{(\mathcal{C})}_{n}(\bm{e})$.

Defined in this way, empirical completeness underlies other well-known notions
of classicality. As the distinguishability of a set of preparations forms an
empirically falsifiable operational property, empirical completeness directly
implies (symmetric) maximal $\psi$-epistemicity [12], bounded ontological
distinctness of preparations [4], and preparation noncontextuality [3].
Equipped with certain quantum theory dependent assumptions, empirical
completeness can also be shown to imply [13, 14, 4] generalized
noncotextuality [3], Kochen-Specker noncontextuality [2] and Bell local-
causality [1]. Consequently, the quantum violation of the operational
consequences of all other well-known notions of classicality which typically
constitute statistical inequalities, imply the operational incompleteness of
quantum theory.

To bring forth the characteristic empirically falsifiable operational
consequences of empirical completeness we consider the following elemental
operational properties associated with a set of $n\geqslant 2$ operational
preparations $\bm{P}\equiv\\{P_{x}\\}^{n}_{x=1}$ and its two member subsets,
respectively,

_$(i.)$ Average set distinguishability ($\bar{\cancel{D}}_{n}$)_ is the
average maximum success probability of correctly guessing which non-trivial
$m$-member subset a given preparation $P_{x}$ belongs to,

$\bar{\cancel{D}}^{(\mathcal{O})}_{n}(\bm{P})=\frac{1}{n-1}\sum^{n-1}_{m=1}\cancel{D}^{(\mathcal{O})}_{n,m}(\bm{P}),$
(5)

where
$\cancel{D}^{(\mathcal{O})}_{n,m}(\bm{P})=\frac{1}{n}\max_{M}\sum_{i_{1}<\ldots<i_{m}}\sum_{x\in\\{i_{1},\ldots,i_{m}\\}}\\{p(k=\\{i_{1},\ldots,i_{m}\\}|P_{x},M)\\}$
forms an operational property of entire set of preparations
$\bm{P}\equiv\\{P_{x}\\}^{n}_{x=1}$, where the first summation is over all
distinct $m$-member subsets $\\{i_{1},\ldots,i_{m}\\}\subset[n]$. We note here
that $\cancel{D}^{(\mathcal{O})}_{n,1}(\bm{P})$ corresponds to the
distinguishability of the set of preparations $\bm{P}$, whereas
$\cancel{D}^{(\mathcal{O})}_{n,n-1}(\bm{P})$ is equivalent to their anti-
distinguishability [15, 13].

_$(ii.)$ Average pair-wise distinguishability ($\bar{D}^{(\mathcal{O})}_{n}$)_
is the average of maximum success probability of perfectly distinguishing
distinct pairs of preparations $\\{P_{i},P_{j}\in\\{P_{x}\\}\\}$,

$\bar{D}^{(\mathcal{O})}_{n}(\bm{P})=\frac{1}{\binom{n}{2}}\sum_{i<j}\cancel{D}^{(\mathcal{O})}_{2,1}(\\{P_{i},P_{j}\\}),$
(6)

where
$\cancel{D}^{(\mathcal{O})}_{2,1}(\\{P_{i},P_{j}\\})=\frac{1}{2}\max_{M\in\mathcal{M}_{\mathcal{O}}}\\{\sum_{x\in\\{i,j\\}}p(k=x|P_{x},M)\\}$
forms an operational property of the pair $\\{P_{i},P_{j}\\}$, commonly
referred to as their distinguishability. Now we are prepared to present the
operational consequences of _empirical completeness_ , namely, a family of
empirically falsifiable statistical _equalities_ ,

######  Theorem 1.

For any _empirically complete_ theory $(\mathcal{EC})$, for any given set of
$n$ preparations $\bm{P}\equiv\\{P_{x}\\}^{n}_{x=1}$, the average set
distinguishability is exactly equal to average pair-wise distinguishability,
i.e.,

$\bar{D}^{(\mathcal{EC})}_{n}(\bm{P})=\bar{\cancel{D}}^{(\mathcal{EC})}_{n}(\bm{P}).$
(7)

The proof has been deferred to the supplementary material for brevity. These
equalities imply that for any empirically complete theory, and any set of $n$
preparations $\bm{P}$, if there exists measurements that achieve an average
pair-wise distinguishability $\bar{D}^{(\mathcal{EC})}_{n}(\bm{P})=p$, there
must exist measurements that attain an average set distinguishability with
exactly the same efficiency, i.e.,
$\bar{\cancel{D}}^{(\mathcal{EC})}_{n}(\bm{P})=p$, and vice versa. On the
contrary, if either average pair-wise distinguishability or average set
distinguishability falls short of the other, i.e., there is non-zero deviation
from the equalities
$\Delta^{(\mathcal{O})}_{n}(\bm{P})=\bar{D}^{(\mathcal{O})}_{n}(\bm{P})-\bar{\cancel{D}}^{(\mathcal{O})}_{n}(\bm{P})\neq
0$, then this signifies the lack of operational measurements which would
saturate the equalities. The deviation $\Delta^{(\mathcal{O})}_{n}(\bm{P})$
forms a measure of incompleteness of the operational theory with respect to
the set of preparations $\bm{P}$. Moreover, it is easy to see that any realist
model underlying such an operational theory must fine-tune the set of
operationally accessible response schemes to explain the operational
deviation.

We now demonstrate that such instances of incompleteness of an operational
theory fuel advantage over classical communication in one-way communication
tasks. To demonstrate this, we cast set distinguishability as an one way
communication task, wherein the sender (Alice) encodes her classical input
$x\in[n]$ onto an operational preparation
$P_{x}\in\bm{P}\equiv\\{P_{x}\in\mathcal{P}_{\mathcal{O}}\\}^{n}_{x=1}$, and
transmits it to the receiver (Bob). Bob upon receiving the transmission
performs a $\binom{n}{y}$-outcome measurement
$M^{y}\in\bm{M}\equiv\\{M^{y}\in\mathcal{M}_{\mathcal{O}}\\}$ based on her
input $y\in[n-1]$. They aim to maximize the success metric
$\bar{\cancel{D}}^{(\mathcal{O})}_{n}(\bm{P},\bm{M})$. Unlike communication
complexity problems, here the amount of communication remains unconstrained,
instead Alice’s preparations are constrained such that their average pair-wise
distinguishability $\bar{D}^{(\mathcal{O})}_{n}(\bm{P})$ is at most $p$. Now
as classical communication is empirically complete, (7) implies that the
success metric is bounded by $p$, i.e.,
$\bar{\cancel{D}}^{(\mathcal{C})}_{n}=\max_{\\{\bm{e},\bm{d}\\}}\bar{\cancel{D}}^{(\mathcal{C})}_{n}(\bm{e},\bm{d})\leqslant
p$. Similarly, average pair-wise distinguishability can also be cast as a
promised communication task with the success metric
$\bar{D}^{(\mathcal{O})}_{n}(\bm{P},\bm{M})$ and the channel constraint
$\bar{\cancel{D}}^{(\mathcal{O})}_{n}(\bm{P})\leqslant p$, such that for such
tasks $\bar{D}^{(\mathcal{C})}_{n}\leqslant p$. If a set of operational
preparations $\bm{P}$, have $\Delta^{(\mathcal{O})}_{n}(\bm{P})<0$ they fuel
advantage in the former task, and when $\Delta^{(\mathcal{O})}_{n}(\bm{P})>0$
an advantage is obtained in the latter.

Consider the set of three qubit preparations
$\bm{\rho}_{\triangle}\equiv\\{\rho_{x}=\frac{\mathbb{I}+\bm{n_{x}}\cdot\bm{\sigma}}{2}\\}^{3}_{x=1}$
which form an equilateral triangle on the Bloch sphere, specifically with
Bloch vectors
$\bf{n_{1}}=[1,0,0]^{T},\bf{n_{2}}=[\cos{\frac{2\pi}{3}},\sin{\frac{2\pi}{3}},0]^{T}$
and $\bf{n_{3}}=[\cos{\frac{4\pi}{3}},\sin{\frac{4\pi}{3}},0]^{T}$. While
their average set distinguishability turns out to be
$\bar{\cancel{D}}^{(\mathcal{Q})}_{3}(\bm{\rho}_{\triangle})=\frac{5}{6}\approx
0.833$, their average pairwise distinguishability is
$\bar{D}^{(\mathcal{Q})}_{3}(\bm{\rho}_{\triangle})=\frac{1}{2}(1+\frac{\sqrt{3}}{2})\approx
0.933$, such that the deviation from equality is
$\Delta^{(\mathcal{Q})}_{3}(\bm{\rho}_{\triangle})=\frac{3\sqrt{3}-4}{12}\approx
0.0997$. The set of four qubit preparations
$\bm{\rho}_{\raisebox{-1.05487pt}{\raisebox{3.01389pt}{\rotatebox{200.0}{\scalebox{1.0}[0.3]{$\bigtriangledown$}}}\raisebox{-0.60275pt}{\rotatebox{80.0}{\scalebox{1.0}[0.3]{$\bigtriangledown$}}}\raisebox{4.76193pt}{\rotatebox{320.0}{\scalebox{1.0}[0.3]{$\bigtriangledown$}}}
}}$ which form a regular tetrahedron on the Bloch sphere attains a higher
deviation from the four preparation equality, i.e.,
$\Delta^{\mathcal{Q}}_{4}(\bm{\rho}_{\raisebox{-1.05487pt}{\raisebox{3.01389pt}{\rotatebox{200.0}{\scalebox{1.0}[0.3]{$\bigtriangledown$}}}\raisebox{-0.60275pt}{\rotatebox{80.0}{\scalebox{1.0}[0.3]{$\bigtriangledown$}}}\raisebox{4.76193pt}{\rotatebox{320.0}{\scalebox{1.0}[0.3]{$\bigtriangledown$}}}
}})\approx 0.1453$.

While these form lower bounds specific to the set of states $\bm{\rho}$,
finding out the maximal quantum deviation from the equalities (7)
$|\Delta^{(\mathcal{Q})}_{n}|=\max_{\bm{\rho}}\\{|\Delta^{(\mathcal{Q})}_{n}(\bm{\rho})|\\}$
boils down to finding the maximal quantum values of a $n$ preparation
preparation success metric, i.e.,
$S^{(\mathcal{Q})}_{n}=\max_{\bm{\rho}}\\{S^{(\mathcal{Q})}_{n}(\bm{\rho})\\}$
given an upper bound on another generic operational property, i.e.,
$T^{(\mathcal{Q})}_{n}(\bm{\rho})\leqslant p$, where
$S^{(\mathcal{Q})}_{n}(\bm{\rho})$, and $T^{(\mathcal{Q})}_{n}(\bm{\rho})$ are
arbitrary operational properties of $n$ density operators of the form (1).
Such optimization problems are very arduous to solve as the dimension of the
quantum systems remains unconstrained. To this end, in a forthcoming article
[16], we device a hierarchy of SDP relaxations for such problems 222Similar to
the Navascués–Pironio–Acín hierarchy for nonlocal quantum correlations [17],
and [18, 19] for generalized contextuality scenarios, [20] for informationally
restricted correlations., to obtain efficient dimension independent tightening
upper bounds $S^{(\mathcal{Q}_{\mathcal{L}})}_{n}$, where
$\mathcal{L}\in\mathbb{N}_{+}$ is the level of the relaxation such that
$\forall\mathcal{L}\in\mathbb{N}_{+}:\
S^{(\mathcal{Q}_{\mathcal{L}})}\geqslant
S^{(\mathcal{Q}_{\mathcal{L}+1})_{n}}\geqslant S^{(\mathcal{Q})}_{n}$. We also
device a see-saw SDP algorithm to obtain dimension dependent efficient lower
bounds $S^{(\mathcal{Q}_{LB})}_{n}\leqslant S^{(\mathcal{Q})}_{n}$. Whenever
the upper bounds from the relaxations coincide with the dimension dependent
lower bounds, we obtain the maximal quantum value $S^{(\mathcal{Q})}_{n}$
(upto machine precision).

We employed the hierarchy of SDP relaxations for the optimization problem
$\bar{D}^{(\mathcal{Q})}_{n}=\max_{\bm{\rho}}\bar{D}^{(\mathcal{Q})}_{n}(\bm{\rho})$
such that $\bar{\cancel{D}}^{(\mathcal{Q})}_{n}(\bm{\rho})\leqslant p$.
Remarkably, the second level of the hierarchy yields a maximum deviation from
the three preparation equality saturated by $\bm{\rho}_{\triangle}$, i.e.,
$\Delta^{(\mathcal{Q}_{2})}_{3}=\Delta^{(\mathcal{Q})}_{3}(\bm{\rho}_{\triangle})=\frac{3\sqrt{3}-4}{12}$
when
$p=\bar{\cancel{D}}^{(\mathcal{Q})}_{3}(\bm{\rho}_{\triangle})=\frac{5}{6}$.
Quantum preparations can also violate the inequality the other way. We
employed the SDP hierarchy and the see-saw technique to obtain upper and lower
bounds for the solution of the inverse optimization problem,
$\bar{\cancel{D}}^{(\mathcal{Q})}_{n}=\max_{\bm{\rho}}\bar{\cancel{D}}^{(\mathcal{Q})}_{n}(\bm{\rho})$
such that $\bar{D}^{(\mathcal{Q})}_{n}(\bm{\rho})\leqslant p$. We find that
the third level of the hierarchy yields the maximum deviation
$-\Delta^{(\mathcal{Q}_{3})}_{3}=D^{(\mathcal{Q})}_{3}(\bm{\cancel{\rho}})\approx
0.0277$ when $p=\bar{D}^{(\mathcal{Q})}_{3}(\bm{\cancel{\rho}})\approx 0.8214$
saturated by a triplet of qutrits $\bm{\cancel{\rho}}$. This implies that the
maximal quantum deviation from the three preparation equality is
$|\Delta^{(\mathcal{Q})}_{3}|=\Delta^{(\mathcal{Q})}_{3}(\bm{\rho}_{\triangle})=\frac{3\sqrt{3}-4}{12}$.
Using these SDP techniques, we characterized the quantum deviations from the
three preparation equality (see FIG. 1). Moreover, the second level of the
hierarchy yields a higher maximum deviation from the four preparation equality
$|\Delta^{(\mathcal{Q})}_{4}|=\Delta^{(\mathcal{Q})}_{4}(\bm{\rho}_{\raisebox{-1.05487pt}{\raisebox{3.01389pt}{\rotatebox{200.0}{\scalebox{1.0}[0.3]{$\bigtriangledown$}}}\raisebox{-0.60275pt}{\rotatebox{80.0}{\scalebox{1.0}[0.3]{$\bigtriangledown$}}}\raisebox{4.76193pt}{\rotatebox{320.0}{\scalebox{1.0}[0.3]{$\bigtriangledown$}}}
}})\approx 0.1453$. In fact, we present evidence that the maximal quantum
deviation from the equalities increases with the number of preparations in
FIG. 2.

In summary, we introduced a realist notion of classicality, termed empirical
completeness, which requires empirically falsifiable operational properties
associated with sets of preparations, to have a not fine-tuned realist
explanation, as properties of stochastic distributions over states of reality.
However, unlike the other notions, the distinguishing experimentally
falsifiable operational consequence of empirical completeness constitute a
family of equalities among elemental operational properties of sets of
preparations. Quantum preparations violate these equalities both ways,
implying the theory’s incompleteness. In general, operational violation of
empirical incompleteness highlights the operational properties which require
realist fine-tuning, and fuels advantage over classical communication in
communication tasks wherein the sender’s preparations are constrained such
that they attain bounded success in other communication tasks.

Figure 1:  _Characterising quantum deviation from the three preparation
equality:_ The average pair-wise distiguishability
$\bar{D}^{(\mathcal{Q})}_{3}(\cdot)$ is plotted against the average set
distinguishability $\bar{\cancel{D}}^{(\mathcal{Q})}_{3}(\cdot)$ for randomly
sampled (from an Haar uniform distribution) triplets of density operators
$\bm{\rho}$ (green crosses), pure states $\bm{\psi}$ (orange solid circles).
The solid black line represents the three preparation equality
$\bar{D}^{(\mathcal{EC})}_{3}(\cdot)=\bar{\cancel{D}}^{(\mathcal{EC})}_{3}(\cdot)$
satisfied by all empirically complete ($\mathcal{EC}$) theories (7). All
random samples deviated from the equality with over $99\%$ exhibiting a
deviation $|\Delta^{(\mathcal{Q})}_{3}(\cdot)|\geqslant 0.00005$. The red
solid line represents the dimension independent upper bounds
$\bar{D}^{(\mathcal{Q}_{3})}_{3}$ such that
$\bar{\cancel{D}}^{(\mathcal{Q})}_{3}(\bm{\rho})=p$ and
$\bar{\cancel{D}}^{(\mathcal{Q}_{3})}_{3}$ when
$\bar{D}^{(\mathcal{Q})}_{3}(\bm{\rho})=p$ from the third level of the
hierarchy of SDP relaxations. The blue circles represent lower bounds
$\bar{D}^{(\mathcal{Q}_{LB})}_{3}$ and
$\bar{\cancel{D}}^{(\mathcal{Q}_{LB})}_{3}$ from the see-saw technique. The
triplet of qubit states $\bm{\rho}_{\triangle}$ which form an equilateral
triangle on the Bloch sphere attain the maximal quantum deviation from the
equality
$|\Delta^{(\mathcal{Q})}_{3}|=\bar{D}^{(\mathcal{Q})}_{3}(\bm{\rho}_{\triangle})-\bar{\cancel{D}}^{(\mathcal{Q})}_{3}(\bm{\rho}_{\triangle})=\frac{3\sqrt{3}-4}{12}\approx
0.0997$. Whereas, a triplet of qutrit states $\bm{\cancel{\rho}}$ attain the
maximal quantum difference
$\bar{\cancel{D}}^{(\mathcal{Q})}_{3}(\bm{\cancel{\rho}})-\bar{D}^{(\mathcal{Q})}_{3}(\bm{\cancel{\rho}})\approx
0.0277$.  Figure 2:  _Increasing quantum deviation:_ This graphic provides
evidence for increasing maximum quantum deviation
$|\Delta^{\mathcal{Q}}_{n}|=\Delta^{\mathcal{Q}}_{n}$ from the equalities (7)
with the number of preparations $n$. The larger blue bars represent upper
bounds $\Delta^{\mathcal{Q}_{1}}_{n}$ obtained from the first level of the
hierarchy of SDP relaxations, while the smaller orange bars represent lower
bounds $\Delta^{\mathcal{Q}_{LB}}_{n}$ obtained via sea-saw SDP technique. For
$n=\\{3,4,5\\}$ the yellow circles represent upper bounds from the second
level $\Delta^{\mathcal{Q}_{2}}_{n}$ of the hierarchy of SDP relaxations which
saturate the lower bounds for maximal quantum deviation
$|\Delta^{\mathcal{Q}_{2}}_{3}|=\Delta^{\mathcal{Q}}_{3}(\bm{\rho}_{\triangle})$,
$|\Delta^{\mathcal{Q}_{2}}_{4}|=\Delta^{\mathcal{Q}}_{4}(\bm{\rho}_{\raisebox{-1.05487pt}{\raisebox{3.01389pt}{\rotatebox{200.0}{\scalebox{1.0}[0.3]{$\bigtriangledown$}}}\raisebox{-0.60275pt}{\rotatebox{80.0}{\scalebox{1.0}[0.3]{$\bigtriangledown$}}}\raisebox{4.76193pt}{\rotatebox{320.0}{\scalebox{1.0}[0.3]{$\bigtriangledown$}}}
}})$, $|\Delta^{\mathcal{Q}_{2}}_{5}|=\Delta^{\mathcal{Q}_{LB}}_{5}$, and
hence constitute the maximal quantum deviations $|\Delta^{\mathcal{Q}}_{n}|$.

###### Acknowledgements.

AC and MP acknowledge financial support by the Foundation for Polish Science
(IRAP project, ICTQT, contract No. 2018/MAB/5, co-financed by EU within Smart
Growth Operational Programme). AC acknowledges financial support by NCN grant
SHENG 2018/30/Q/ST2/00625. DS acknowledges financial support by National Post-
doctoral Fellowship, India. MP financial support by QuantERA grant funded by
NCBiR QUANTERA/2/2020. AC acknowledges Racecar by Periphery for creative
support.

## References

  * Bell [1964] J. S. Bell, On the einstein podolsky rosen paradox, Physics Physique Fizika 1, 195 (1964).
  * KOCHEN and SPECKER [1967] S. KOCHEN and E. P. SPECKER, The problem of hidden variables in quantum mechanics, Journal of Mathematics and Mechanics 17, 59 (1967).
  * Spekkens [2005] R. W. Spekkens, Contextuality for preparations, transformations, and unsharp measurements, Phys. Rev. A 71, 052108 (2005).
  * Chaturvedi and Saha [2020] A. Chaturvedi and D. Saha, Quantum prescriptions are more ontologically distinct than they are operationally distinguishable, Quantum 4, 345 (2020).
  * Catani and Leifer [2020] L. Catani and M. Leifer, A mathematical framework for operational fine tunings, preprint arXiv:2003.10050 (2020).
  * BELL [1966] J. S. BELL, On the problem of hidden variables in quantum mechanics, Rev. Mod. Phys. 38, 447 (1966).
  * Brunner  _et al._ [2014] N. Brunner, D. Cavalcanti, S. Pironio, V. Scarani, and S. Wehner, Bell nonlocality, Rev. Mod. Phys. 86, 419 (2014).
  * Spekkens  _et al._ [2009] R. W. Spekkens, D. H. Buzacott, A. J. Keehn, B. Toner, and G. J. Pryde, Preparation contextuality powers parity-oblivious multiplexing, Phys. Rev. Lett. 102, 010401 (2009).
  * Saha and Chaturvedi [2019] D. Saha and A. Chaturvedi, Preparation contextuality as an essential feature underlying quantum communication advantage, Phys. Rev. A 100, 022108 (2019).
  * Hameedi  _et al._ [2020] A. Hameedi, B. Marques, P. Mironowicz, D. Saha, M. Pawłowski, and M. Bourennane, Experimental test of nonclassicality with arbitrarily low detection efficiency, Phys. Rev. A 102, 032621 (2020).
  * Boyd and Vandenberghe [2004] S. Boyd and L. Vandenberghe,  _Convex optimization_ (Cambridge university press, 2004).
  * Barrett  _et al._ [2014] J. Barrett, E. G. Cavalcanti, R. Lal, and O. J. E. Maroney, No $\psi$-epistemic model can fully explain the indistinguishability of quantum states, Phys. Rev. Lett. 112, 250403 (2014).
  * Leifer [2014] M. Leifer, Is the quantum state real? an extended review of $\psi$-ontology theorems, Quanta 3, 67 (2014).
  * Leifer and Maroney [2013] M. S. Leifer and O. J. E. Maroney, Maximally epistemic interpretations of the quantum state and contextuality, Phys. Rev. Lett. 110, 120401 (2013).
  * Caves  _et al._ [2002] C. M. Caves, C. A. Fuchs, and R. Schack, Conditions for compatibility of quantum-state assignments, Phys. Rev. A 66, 062111 (2002).
  * [16] E. Panwar, A. Chaturvedi, and et al., Characterising classical and quantum correlations constrained by bounded success metrics of one-way communication tasks, forthcoming.
  * Navascués  _et al._ [2008] M. Navascués, S. Pironio, and A. Acín, A convergent hierarchy of semidefinite programs characterizing the set of quantum correlations, New J. Phys. 10, 073013 (2008).
  * Chaturvedi  _et al._ [2021] A. Chaturvedi, M. Farkas, and V. J. Wright, Characterising and bounding the set of quantum behaviours in contextuality scenarios, Quantum 5, 484 (2021).
  * Tavakoli  _et al._ [2021] A. Tavakoli, E. Z. Cruzeiro, R. Uola, and A. A. Abbott, Bounding and simulating contextual correlations in quantum theory, PRX Quantum 2, 020334 (2021).
  * Tavakoli  _et al._ [2020] A. Tavakoli, E. Z. Cruzeiro, E. Woodhead, and S. Pironio, Informationally restricted correlations: a general framework for classical and quantum systems, arXiv preprint arXiv:2007.16145 (2020).
  * Ross [1976] S. M. Ross,  _A first course in probability_ , Tech. Rep. (1976).

## Supplemental material

##  Appendix A Proof of Theorem 1

Theorem 1 states that for any _empirically complete_ theory, for any given set
of $n$ preparations $\bm{P}\equiv\\{P_{x}\\}^{n}_{x=1}$, the average set
distinguishability $\bar{D}^{(\mathcal{EC})}_{n}(\bm{P})$ is exactly equal to
average pair-wise distinguishability
$\bar{\cancel{D}}^{(\mathcal{EC})}_{n}(\bm{P})$, i.e.
$\bar{D}^{(\mathcal{EC})}_{n}(\bm{P})=\bar{\cancel{D}}^{(\mathcal{EC})}_{n}(\bm{P})$.
To prove this thesis, first, for a set of epistemic states
$\bm{\mu}\equiv\\{\mu_{x}\\}^{n}_{x=1}$, we obtain expressions for the not
fine-tuned realist counterparts of average set distinguishability
$\bar{\cancel{D}}^{(\Lambda)}_{n}(\bm{\mu})$ and average pair-wise
distinguishability $\bar{D}^{(\Lambda)}_{n}(\bm{\mu})$,

$\displaystyle\bar{D}^{(\Lambda)}_{n}(\bm{\mu})$
$\displaystyle=\frac{1}{\binom{n}{2}}\sum_{i<j}\cancel{D}^{(\Lambda)}_{2,1}(\\{\mu_{i},\mu_{j}\\})$
(8)
$\displaystyle=\frac{1}{n(n-1)}\sum_{i<j}\int_{\Lambda}d\lambda\max\\{\mu_{i}(\lambda),\mu_{j}(\lambda)\\},$
$\displaystyle\bar{\cancel{D}}^{(\Lambda}_{n}(\bm{\mu})$
$\displaystyle=\frac{1}{n-1}\sum^{n-1}_{m=1}\cancel{D}^{(\Lambda)}_{n,m}(\bm{\mu})$
$\displaystyle=\frac{1}{n(n-1)}\sum^{n-1}_{m=1}\int_{\Lambda}d\lambda\max_{i_{1}<\ldots<i_{m}\in[n]}\sum_{x\in\\{i_{1},\ldots,i_{m}\\}}\\{\mu_{x}(\lambda)\\},$

where the second and fourth equality follows from
$\cancel{D}^{(\Lambda)}_{n,m}=\frac{1}{n}\int_{\Lambda}d\lambda\max_{i_{1}<\ldots<i_{m}\in[n]}\sum_{x\in\\{i_{1},\ldots,i_{m}\\}}\\{\mu_{x}(\lambda)\\}$
which in-turn follows from (3). Now we introduce following elementary identity
[21] which forms the key ingredient of our proof technique,

######  Lemma 1.

For any set of $n$ real numbers $\\{u_{x}\in\mathbb{R}\\}^{n}_{x=1}$, the
following identity holds,

$\sum_{i<j}\max\\{u_{i},u_{j}\\}=\sum^{n-1}_{m=1}\max_{i_{1}<\ldots<i_{m}\in[n]}\sum_{x\in\\{i_{1},\ldots,i_{m}\\}}\\{u_{x}\\}.$
(9)

###### Proof.

To prove the above thesis, we consider an ordered list $\\{a_{x}\\}^{n}_{x=1}$
such that $a_{1}\geqslant a_{2}\ldots\geqslant a_{n-1}\geqslant a_{n}$
associated to a given set of real numbers $\\{u_{x}\\}^{n}_{x=1}$.
Consequently, we re-express $\sum_{i<j}\max\\{u_{i},u_{j}\\}$ in terms of
members of the ordered list $\\{a_{x}\\}^{n}_{x=1}$ as,

$\sum_{i<j}\max\\{u_{i},u_{j}\\}=\sum^{n-1}_{x=1}(n-x)a_{x}\\\ $ (10)

Notice that the right hand side of the above equation can be regrouped into
sums of $m\leqslant n-1$ largest numbers such that,

$\sum_{i<j}\max\\{u_{i},u_{j}\\}=(n-1)a_{1}+(n-2)a_{2}+\ldots+2a_{n-2}+a_{n-1}=\sum^{n-1}_{m=1}\sum^{m}_{x=1}a_{x},$
(11)

Finally, observe that for any $m\leqslant n$,

$\max_{i_{1}<\ldots<i_{m}\in[n]}\sum_{x\in\\{i_{1},\ldots,i_{m}\\}}\\{u_{x}\\}=\sum^{m}_{x=1}a_{x},$
(12)

which when plugged back into (11) yields the desired thesis (9). ∎

As any set of epistemic states $\bm{\mu}\equiv\\{\mu_{x}\\}^{n}_{x=1}$ yields
a set of positive numbers $\\{\mu_{x}{\lambda}\\}^{n}_{x=1}$ for each ontic
state $\lambda\in\Lambda$, Lemma 1 implies,

$\sum_{i<j}\max\\{\mu_{i}(\lambda),\mu_{j}(\lambda)\\}=\sum^{n-1}_{m=1}\max_{i_{1}<\ldots<i_{m}\in[n]}\sum_{x\in\\{i_{1},\ldots,i_{m}\\}}\\{\mu_{x}(\lambda)\\},\
\forall\lambda\in\Lambda.$ (13)

Integrating both sides of (13) over the ontic-state space, and multiplying by
$\frac{1}{n(n-1)}$ yields,

$\bar{D}^{(\Lambda)}_{n}(\bm{\mu})=\frac{1}{n(n-1)}\int_{\Lambda}d\lambda\sum_{i<j}\max\\{\mu_{i}(\lambda),\mu_{j}(\lambda)\\}=\int_{\Lambda}d\lambda\sum^{n-1}_{m=1}\max_{i_{1}<\ldots<i_{m}\in[n]}\sum_{x\in\\{i_{1},\ldots,i_{m}\\}}\\{\mu_{x}(\lambda)\\}=\bar{D}^{(\Lambda)}_{n}(\bm{\mu}),$
(14)

where we have used the relations (8). Finally, as for any _empirically
complete_ theory, for any set of preparations
$\bm{P}\equiv\\{P_{x}\\}^{n}_{x=1}$ there exists an underlying set of
epistemic states $\bm{\mu}\equiv\\{\mu_{x}\\}^{n}_{x=1}$, such all that
operational properties have a not fine-tuned realist explanation, i.e.,
$S^{(\mathcal{EC})}_{n}=S^{(\Lambda)}_{n}$. In particular, empirical
completeness allows us to port the relation (14) to the operational level such
that
$\bar{D}^{(\mathcal{EC})}_{n}(\bm{P})=\bar{\cancel{D}}^{(\mathcal{EC})}_{n}(\bm{P})$,
which concludes the proof.


\prod_{\rho \in \Sigma(1)}\lambda_\rho^{{b_\rho}} \cdot \prod_{\rho \in \Sigma(1)}z_\rho^{a_\rho}=\prod_{\rho \in \Sigma(1)}(z_\rho \lambda_\rho)^{a_\rho}~\text{for all $\rho \in \Sigma(1)$},
\]

which is equivalent to the equations $\prod_{\rho \in \Sigma(1)}\lambda_\rho^{b_\rho-a_\rho}=1$ for all $\rho \in \Sigma(1)$. By exactness of (<ref>), this is equivalent to the existence of $m \in {\rm M}$ such that $a_\rho-{b_\rho}={|E_\rho|}\langle m,u_\rho \rangle$ for all $\rho \in \Sigma(1)$.

§ COMBINATORIAL RULES FOR THE LINEAR SYSTEM OF EQUATIONS VBW

In this appendix, we collect some combinatorial rules that facilitate the computation of the small polygons and equations in $\mathbb V_{\bw \w}$.

§.§ Equivalent description of the small polygons

Consider the lines
L_ρ:= { m ∈M_: ⟨m,u_ρ⟩=-b_ρ}
that form the boundary of the small Newton polygon $N_{\bw \w}$.  We give an alternate description of these lines.  Recall that $\widetilde \Gamma$ be the biperiodic graph on the plane given by the lift of $\Gamma$ to the universal cover of $\T$. The zig-zag paths in $\widetilde \Gamma$ for a given $\rho$ divide the plane into an infinite collection of strips $\mathscr S_\rho(d)$  parameterized by $d \in \frac 1 {|E_\rho|} \mathbb Z$ such that

\[
\mathscr S_\rho(d) \cap V(\widetilde \Gamma)=\{{\rm v} \in V(\widetilde \Gamma): {[D_\rho]}\dd({\rm v})=d\},
\]

where for a divisor $D$, $[D_\rho] D$ denotes the coefficient of $D_\rho$ in $D$. We assign to each strip $\mathscr S_\rho(d)$  a line $L_\rho(d)$ in ${\rm M}_\R$ parallel to $E_\rho$, using the following rule illustrated on  Figure <ref>:

* The line associated to a strip  ${\mathscr  S}_\rho(d)$  contains the side $E_\rho$ if and only if either

i)  The strip ${\mathscr S}_\rho(d)$ is on the right (when facing in the direction of the path) of a zig-zag path $\alpha_1 \in Z_\rho$, and   $\alpha_1$ contains $\bw$.

ii) The strip ${\mathscr  S}_\rho(d)$ contains   $\bw$, and  $\bw$ is not in a zig-zag path in $Z_\rho$, or

* Moving to the strip to the left shifts the line $1/|E_\rho|$ steps to the left.

We call the strip to the left of the one whose line contains $E_\rho$, and all strips obtained by  its translations by $H_1(\T,\Z)$, exceptional strips.

If we associate lines to strips as above, the boundary of the small Newton polygon $N_{\bw \w}$ is given by the lines $\{L_\rho(d_\rho)\}$, where $d_\rho \in \frac 1 {|E_\rho|}\Z$ is determined by the condition $\w \in \mathscr S_\rho(d_\rho)$, that is, it is the index of the strip containing $\w$ in the direction $\rho$.

In order for the line $L_\rho$ in (<ref>) to contain $E_\rho$, we must have $b_\rho=0,$ where $b_\rho$ is the coefficient of $D_\rho$ in (<ref>). We have to consider two cases.

* There is a zig-zag path $\alpha \in Z_\rho$ such that $\bw $ is contained in $\alpha$. We need $[D_\rho](\dd(\bw))=\frac 1 {|E_{\rho}|}+ [D_\rho](\dd(\w))$, which means $\w$ is contained in the strip $\mathscr S$ to the right of the one containing $\bw$, with $\alpha$ separating the two strips.

* No zig-zag path in $Z_\rho$ contains $\bw$. In this case, we need the coefficients  $[D_\rho](\dd(\bw))=[D_\rho](D(\w))$, which means $\w$ is in the strip $\mathscr S$ containing $\bw$.

If $\w_2$ is a white vertex in the strip to the left of the strip containing a white $\w_1$, then   $[D_\rho](\dd(\w_2))=[D_\rho](\dd(\w_1))+\frac{1}{|E_\rho|}$. So if we define $b_\rho(\w_1)$ and $b_\rho(\w_2)$ as in (<ref>) with $\w=\w_1$ and $\w=\w_2$ respectively, then $b_\rho(\w_2)=b_\rho(\w_1)+\frac 1 {|E_\rho|}$. Note that the line (<ref>) which bounds $N_{\bw \w_2}$ is given by

\[
L_\rho(\w_2):=\{m \in {\rm M}_\R : \langle m,u_\rho \rangle = b_\rho(\w_2)\}.
\]

The similar line  which bounds  $N_{\bw \w_1}$ is

\[
L_\rho(\w_1):=\{m \in {\rm M}_\R : \langle m,u_\rho \rangle = b_\rho(\w_1)\},
\]

so the line  $L_\rho(\w_2)$ is obtained from the line $L_\rho(\w_1)$ by shifting $1/{|E_\rho|}$ steps to the left.

§.§ The equations in Vbw

We describe the equations of type 2 in Section <ref>.

[gray,line width=0.01, fill opacity=0.5, fill=gray] (-1,0) rectangle (4,1);
[gray,line width=0.01, fill opacity=0.5, fill=gray] (-1,5) rectangle (4,6);
[fill=black] (0,0) circle (4pt);
[fill=black] (0,5) circle (4pt);
[] (no) at (-1.5,0)$\alpha_1$;
[] (no) at (-1.5,1)$\alpha_k$;
[] (no) at (-1.5,2)$\alpha_{k-1}$;
[] (no) at (-1.5,4)$\alpha_{2}$;
[] (no) at (-1.5,5)$\alpha_{1}$;
[] (no) at (-1.5,6)$\alpha_k$;
[->](-1,0) – (4,0);
[->](-1,1) – (4,1);
[->](-1,2) – (4,2);
[->](-1,4) – (4,4);
[->](-1,5) – (4,5);
[->](-1,6) – (4,6);
[blue] (no) at (3,0.5)empty;
[blue] (no) at (3,1.5)$\alpha_k$;
[blue] (no) at (3,4.5)$\alpha_2,\dots,\alpha_k$;
[blue] (no) at (3,2.5)$\alpha_{k-1},\alpha_k$;
[blue] (no) at (3,5.5)empty;
[blue] (no) at (3,6.5)$\alpha_k$;
[](no) at (1.5,3)$\vdots$;

[fill=black] (0,0.5) circle (4pt);
[fill=black] (0,5.5) circle (4pt);
[gray,line width=0.01, fill opacity=0.5, fill=gray] (-1,1) rectangle (4,2);
[gray,line width=0.01, fill opacity=0.5, fill=gray] (-1,6) rectangle (4,7);
[] (no) at (-1.5,0)$\alpha_1$;
[] (no) at (-1.5,1)$\alpha_k$;
[] (no) at (-1.5,2)$\alpha_{k-1}$;
[] (no) at (-1.5,4)$\alpha_{2}$;
[] (no) at (-1.5,5)$\alpha_{1}$;
[] (no) at (-1.5,6)$\alpha_k$;
[->](-1,0) – (4,0);
[->](-1,1) – (4,1);
[->](-1,2) – (4,2);
[->](-1,4) – (4,4);
[->](-1,5) – (4,5);
[->](-1,6) – (4,6);
[blue] (no) at (3,0.5)$\alpha_1,\dots,\alpha_k$;
[blue] (no) at (3,1.5)$\alpha_k$;
[blue] (no) at (3,2.5)$\alpha_{k-1},\alpha_k$;
[blue] (no) at (3,4.5)$\alpha_2,\dots,\alpha_k$;
[blue] (no) at (3,5.5)$\alpha_1,\dots,\alpha_k$;
[blue] (no) at (3,6.5)$\alpha_k$;
[](no) at (1.5,3)$\vdots$;

The lifts of zig-zag paths $\alpha_1,\dots,\alpha_k$ in $Z_\rho$ divide the plane into strips. The side $L_{\rho}$ of the  small polygon $N_{\bw \w}$ and the columns of the matrix  $\mathbb V_{\bw \w}$ are determined by the strips containing $\bw$ and $\w$. The black vertex ${\bw}$ is the black dot. On the left panel, $\bw$ is on a zig-zag path, and on the right, it is between two zig-zag paths. Written inside each strip in blue is the subset of $Z_\rho$ that gives rise to equations in $\mathbb V_{\bw \w}$ if $\w$ is contained in that strip.  Exceptional strips are shaded.

Let $\rho \in \Sigma(1)$ be a ray and let $Z_\rho=\{\alpha_1,\dots,\alpha_k\}$, where $\alpha_1,\dots,\alpha_k$ are labeled in cyclic order. Their lifts  to the universal cover of the torus divides it
into strips, see Figure <ref>. We denote  by ${\mathscr S}_i$ the strip immediately to the right of $\alpha_i$.

The set of extra linear equations is described as follows:

* One of these zig-zag paths contains  ${\bw}$. We can assume it is $\alpha_1$. Then the subset of $Z_\rho$ that contributes an equation to $\mathbb V_{\bw \w}$ is:

\begin{align}
\text{empty} &\text{ if }\w \in {\mathscr S}_k; \nonumber\\
\alpha_{i+1}, \dots, \alpha_k  &\text{ if }\w \in {\mathscr S}_i \text{ for some }i \neq k.\label{colcase1}
\end{align}

* The vertex ${\bw}$ is not in any of  zig-zag paths  in $Z_\rho$. Then the   subset of $Z_\rho$ is

\begin{align}
\alpha_1,\dots, \alpha_k &\text{ if }\w \in {\mathscr S}_k; \nonumber \\
\alpha_{i+1},\dots \alpha_k  &\text{ if }\w \in {\mathscr S}_i, \text{ for some }i \neq k.\label{colcase2}
\end{align}

\[
\divebw_{\bw \w}=D_N-\dd(\w)+\dd(\bw)-\sum_{\rho \in \Sigma(1)}\sum_{\alpha \in Z_\rho : \bw \in \alpha } \frac{1}{|E_{\rho}|}D_{\rho}.
\]
into (<ref>), we first observe that (<ref>) does not change if we replace $\w$ by its any translate  on the universal cover because ${\bf d}(\w)-{\bf d}({ \bw})$ changes by the same amount as $\restr{[{E_{{\bw}\w}}]}{\spectralcurve}$ but with the opposite sign. Therefore we may assume that among all its possible translates in the universal cover, the strip $\mathscr S_i$ is the one that is immediately to the right of ${\bw}$. Then we have

\begin{align*}
\restr{({\bf d}(\w)-{\bf d}({\bw}))}{\spectralcurve \cap D_\rho}&=-\nu(\alpha_1)- \cdots -\nu{(\alpha_i)}\end{align*}

and the coefficient of $D_\rho$ in ${(\dd(\w)-\dd({\bw}))}$ is $ -\frac{i}{k} $.

Now we distinguish two cases:

* Suppose $\bw$ is contained in $\alpha_1$, so that the coefficient of $D_\rho$ in $ \sum_{\rho \in \Sigma(1)}\sum_{\alpha \in Z_\rho : \bw \in \alpha } \frac{1}{|E_{\rho}|}D_{\rho}$ is $\frac 1 k . $ Then we have

\[
\restr{\floor{-\dd(\w)+\dd(\bw)-\sum_{\rho \in \Sigma(1)}\sum_{\alpha \in Z_\rho : \bw \in \alpha } \frac{1}{|E_{\rho}|}D_{\rho}}}{ \spectralcurve \cap D_\rho }=0,
\]

so that (<ref>) is
$\nu({\alpha_{i+1}})+\cdots+\nu({\alpha_k})$, which proves  (<ref>).

* Suppose $\bw$ is not contained in any of the $\alpha_i$, so that the coefficient of $D_\rho$ in  $\sum_{\rho \in \Sigma(1)}\sum_{\alpha \in Z_\rho : \bw \in \alpha } \frac{1}{|E_{\rho}|}D_{\rho}$ is $0.$ We have

\[
\restr{\floor{-\dd(\w)+\dd(\bw)-\sum_{\rho \in \Sigma(1)}\sum_{\alpha \in Z_\rho : \bw \in \alpha } \frac{1}{|E_{\rho}|}D_{\rho}}}{ \spectralcurve \cap D_\rho }=\begin{cases}0 &\text{ if }i \neq k,\\
\sum_{j=1}^k \nu({\alpha_j}) &\text{ if }i=k.
\end{cases}
\]

This gives

\begin{align*}
\nu({\alpha_{i+1}})+\cdots +\nu(\alpha_k) &\text{ if }i \neq k,\\
\sum_{j=1}^k \nu({\alpha_j}) &\text{ if }i=k.\end{cases}\end{align*}

We obtain (<ref>).

author=Affolter, Niklas,
author=George, Terrence,
author=Ramassamy, Sanjay,
title=Cross-ratio dynamics and the dimer cluster integrable system,

author=Boutillier, Cédric,
author=Cimasoni, David,
author=de Tilière, Béatrice,
title=Elliptic dimers on minimal graphs and genus 1 harnack curves,

author=Boutillier, Cédric,
author=Cimasoni, David,
author=de Tilière, Béatrice,
title=Minimal bipartite dimers and higher genus harnack curves,

author=Beauville, Arnaud,
title=Theta functions, old and new,
booktitle=Open problems and surveys of contemporary mathematics,
series=Surv. Mod. Math.,
publisher=Int. Press, Somerville, MA,

author=Borisov, Lev,
author=Hua, Zheng,
title=On the conjecture of King for smooth toric Deligne-Mumford
journal=Adv. Math.,

author=Cohn, Henry,
author=Kenyon, Richard,
author=Propp, James,
title=A variational principle for domino tilings,
journal=J. Amer. Math. Soc.,

author=Cox, David A.,
author=Little, John B.,
author=Schenck, Henry K.,
title=Toric varieties,
series=Graduate Studies in Mathematics,
publisher=American Mathematical Society, Providence, RI,

author=Elkies, Noam,
author=Kuperberg, Greg,
author=Larsen, Michael,
author=Propp, James,
title=Alternating-sign matrices and domino tilings. I,
journal=J. Algebraic Combin.,

author=Elkies, Noam,
author=Kuperberg, Greg,
author=Larsen, Michael,
author=Propp, James,
title=Alternating-sign matrices and domino tilings. II,
journal=J. Algebraic Combin.,

author=Fock, Vladimir V.,
author=Marshakov, Andrey,
title=Loop groups, clusters, dimers and integrable systems,
booktitle=Geometry and quantization of moduli spaces,
series=Adv. Courses Math. CRM Barcelona,
publisher=Birkhäuser/Springer, Cham,

author=Fock, V. V.,
title=Inverse spectral problem for GK integrable systems,

author=Fulton, William,
title=Introduction to toric varieties,
series=Annals of Mathematics Studies,
publisher=Princeton University Press, Princeton, NJ,
note=The William H. Roever Lectures in Geometry,

author=Goncharov, Alexander B.,
author=Kenyon, Richard,
title=Dimers and cluster integrable systems,
journal=Ann. Sci. Éc. Norm. Supér. (4),

author=Hanany, Amihay,
author=Kennaway, Kristian D.,
title=Dimer models and toric diagrams,

author=Kasteleyn, P. W.,
title=The statistics of dimers on a lattice: I. The number of dimer
arrangements on a quadratic lattice,

author=Kasteleyn, P. W.,
title=Dimer statistics and phase transitions,
journal=J. Mathematical Phys.,

author=Kenyon, Richard,
title=Conformal invariance of domino tiling,
journal=Ann. Probab.,

author=Kenyon, R.,
title=The Laplacian and Dirac operators on critical planar graphs,
journal=Invent. Math.,

author=Kenyon, Richard,
title=Local statistics of lattice dimers,
journal=Ann. Inst. H. Poincaré Probab. Statist.,

author=Kenyon, Richard,
author=Okounkov, Andrei,
title=Planar dimers and Harnack curves,
journal=Duke Math. J.,

author=Kenyon, Richard,
author=Okounkov, Andrei,
author=Sheffield, Scott,
title=Dimers and amoebae,
journal=Ann. of Math. (2),

author=Temperley, H. N. V.,
author=Fisher, Michael E.,
title=Dimer problem in statistical mechanics—an exact result,
journal=Philos. Mag. (8),

author=Treumann, David,
author=Williams, Harold,
author=Zaslow, Eric,
title=Kasteleyn operators from mirror symmetry,
journal=Selecta Math. (N.S.),
pages=Paper No. 60, 36,