#  Application of the gradient method to Hartree-Fock-Bogoliubov theory

L.M. Robledo  Departamento de Fisica Teorica, Univeridad Autonoma de Madrid,
E-28049 Madrid, Spain  G.F. Bertsch  Institute for Nuclear Theory and Dept. of
Physics, University of Washington, Seattle, Washington

###### Abstract

A computer code is presented for solving the equations of Hartree-Fock-
Bogoliubov (HFB) theory by the gradient method, motivated by the need for
efficient and robust codes to calculate the configurations required by
extensions of HFB such as the generator coordinate method. The code is
organized with a separation between the parts that are specific to the details
of the Hamiltonian and the parts that are generic to the gradient method. This
permits total flexibility in choosing the symmetries to be imposed on the HFB
solutions. The code solves for both even and odd particle number ground
states, the choice determined by the input data stream. Application is made to
the nuclei in the $sd$-shell using the USDB shell-model Hamiltonian.

##  I Introduction

An important goal of nuclear structure theory is to develop the computational
tools for a systematic description of nuclei across the chart of the nuclides.
There is hardly any alternative to self-consistent mean-field (SCMF) for the
starting point of a global theory, but the SCMF has to be extended by the
generator coordinate method (GCM) or other means to calculate spectroscopic
observables. There is a need for computational tools to carry out the SCMF
efficiently in the presence of the multiple constraints to be used for the
GCM. Besides particle number, quantities that may be constrained include
moments of the density, angular momentum, and in the Hartree-Fock-Bogoliubov
(HFB) theory, characteristics of the anomalous densities.

The gradient method described by Ring and Schuck (RS , Section 7.3.3) is very
suitable for this purpose: it is robust and easily deals with multiple
constraints. However, the actual computational aspects of the method as
applied to HFB have not been well documented in the literature. This is in
contrast to methods based on diagonalizing the HFB matrix eigenvalue equation.
Here there are several codes available in the literature, eg. po97 ; bo05 ;
be05 ; do05 ; st05 . Other, less used, methods to solve the HFB equation with
multiple constraints can be found in the literature; for example the method
described in Ref. eg95  is close in spirit to the one presented here. We note
also that the computational issues for using the gradient method in nuclear
Hartree-Fock theory have been discussed in detail in Ref. re82 . That paper
also contains references to related techniques such as the imaginary time step
method.

Here we will describe an implementation of the gradient algorithm for HFB
following the iterative method used by Robledo and collaborators wa02 . The
code presented here, hfb_shell, is available as supplementary material to this
article (see Appendix). The code has separated out the parts that are basic to
the gradient method and the parts that are specific to the details of the
Hamiltonian. As an example, the code here contains a module for application to
the $sd$-shell with a shell-model Hamiltonian containing one-body and two-body
terms. There is a long-term motivation for this application as well. The
$sd$-shell could be a good testing ground for the extensions of SCMF such as
the GCM and approximations derived from GCM. Since one has a Hamiltonian for
the $sd$-shell that describes the structure very well, one could test the
approximations to introduce correlations, such as projection, the random-phase
approximation, etc and compare them with the exact results from the Shell
Model. Preliminary results along this line are discussed in ro08 ; ma11 . As a
first step in this program, one needs a robust SCMF code that treats shell-
model Hamiltonians. Extensions to other shell model configuration spaces are
straightforward and only limited by the availability of computational
resources.

The code described here is more general than earlier published codes in that
it can treat even or odd systems equally well. The formalism for the extension
to odd systems and to a statistical density matrix will be presented elsewhere
Robledo-Bertsch . We also mention that the present code (with a different
Hamiltonian module) has already been applied to investigate neutron-proton
pairing in heavy nucleige11 .

##  II Summary of the gradient method

The fundamental numerical problem to be addressed is the minimization of a
one- plus two-body Hamiltonian under the set of Bogoliubov transformations in
a finite-dimensional Fock space. We remind the reader of the most essential
equations, using the notation of Ring and Schuck RS . The basic variables are
the $U$ and $V$ matrices defining the Bogoliubov transformation. The main
physical variables are the one-body matrices for the density $\rho$ and the
anomalous density $\kappa$, given by

$\rho=V^{*}V^{t};\,\,\,\,\,\,\kappa=V^{*}U^{t}.$ (1)

The Hamiltonian may be defined in the Fock-space representation as

$\hat{H}=\sum_{12}\varepsilon_{12}c^{\dagger}_{1}c_{2}+{1\over
4}\sum_{1234}v_{1234}c^{\dagger}_{1}c^{\dagger}_{2}c_{4}c_{3}.$ (2)

The expectation value of the Hamiltonian under a Bogoliubov transformation of
the vacuum is given by

$H^{00}\equiv\langle\hat{H}\rangle={\rm Tr}(\varepsilon\rho+\hbox{${1\over
2}$}\Gamma\rho-\hbox{\rm${1\over 2}$}\Delta\kappa^{*}).$ (3)

in terms of the fields for the ordinary potential $\Gamma$ and the pairing
potential $\Delta$. These are defined as

$\Gamma_{12}=\sum_{34}v_{1423}\rho_{34};\,\,\,\,\Delta_{12}=\hbox{${1\over
2}$}\sum_{34}v_{1234}\kappa_{34}.$ (4)

The gradient method makes extensive use of the quasiparticle representation
for operators related to the ordinary and anomalous densities. For a single-
particle operator $\hat{F}=\sum_{ij}F_{ij}c^{\dagger}_{i}c_{j}$ we write

$\sum_{ij}F_{ij}c^{\dagger}_{i}c_{j}\equiv
c^{\dagger}Fc=F^{00}+\beta^{\dagger}F^{11}\beta^{\dagger}+\hbox{${1\over
2}$}\left(\beta F^{02}\beta+\beta^{\dagger}F^{20}\beta^{\dagger}\right).$ (5)

where $\beta,\beta^{\dagger}$ are quasiparticle annihilation and creation
operators. The gradients will be constructed from the skew-symmetric matrix
$F^{20}$, which for a normal one-body operator is given by

$F^{20}=U^{\dagger}FV^{*}-V^{\dagger}F^{t}U^{*}.$ (6)

The corresponding representation for an operator $\hat{G}$ of the anomalous
density is

$\hbox{${1\over
2}$}(c^{\dagger}Gc^{\dagger}-cG^{*}c)=G^{00}+\beta^{\dagger}G^{11}\beta+\hbox{${1\over
2}$}(\beta^{\dagger}G^{20}\beta^{\dagger}+\beta G^{02}\beta)$ (7)

The skew-symmetric matrix $G^{20}$ is given by

$G^{20}=U^{\dagger}GU^{*}-V^{\dagger}G^{*}V^{*}.$ (8)

Two operators that are particularly useful to characterize the HFB states are
the axial quadrupole operator $Q_{Q}$ and the number fluctuation operator
$\Delta N^{2}$. We define $Q_{Q}$ as

$Q_{Q}=2z^{2}-x^{2}-y^{2};$ (9)

its expectation value distinguishes spherical and deformed minima. The number
fluctuation is an indicator of the strength of pairing condensates and is zero
in the absence of a condensate. It depends on the two-body operator
$\hat{N}^{2}$, but like the Hamiltonian can be expressed in terms of one-body
densities. We define it as

$\Delta
N^{2}\equiv\langle\hat{N}^{2}\rangle-\langle\hat{N}\rangle^{2}=\frac{1}{2}{\rm
Tr}\left(N^{20}N^{02}\right)=2{\rm Tr}\left(\rho(1-\rho)\right)=-2{\rm
Tr}\left(\kappa^{*}\kappa\right).$ (10)

The full expansion of the Hamiltonian in the quasiparticle basis is given in
Eqs. (E.20-E.25) of RS . Here we will mainly need $H^{20}$, given by

$H^{20}=h^{20}+\Delta^{20}=U^{\dagger}hV^{*}-V^{\dagger}h^{t}U^{*}-V^{\dagger}\Delta^{*}V^{*}+U^{\dagger}\Delta
U^{*}.$ (11)

where $h=\epsilon+\Gamma$. Starting from any HFB configuration $U,V$ one can
construct a new configuration $U^{\prime},V^{\prime}$ by the generalized
Thouless transformation. The transformation is defined by a skew-symmetric
matrix $Z$ having the same dimensions as $U,V$. One often assumes that the
transformation preserves one or more symmetries such as parity or axial
rotational symmetry. Then the $U,V$ matrices are block diagonal and $Z$ has
the same block structure. Otherwise the elements of $Z$ are arbitrary and can
be real or complex. The transformation is given by

$U^{\prime}=(U+V^{*}Z^{*})(1-ZZ^{*})^{-1/2}=U+V^{*}Z^{*}+{\cal O}(Z^{2})$ (12)
$V^{\prime}=(V+U^{*}Z^{*})(1-ZZ^{*})^{-1/2}=V+U^{*}Z^{*}+{\cal O}(Z^{2}).$

The last factor, $(1-ZZ^{*})^{-1/2}$, ensures that the transformed set
$U^{\prime},V^{\prime}$ satisfies the required unitarity conditions for the
Bogoliubov transformation. We now ask how the expectation value of some
bilinear operator $\hat{Q}$ changes when the Thouless transformation is
applied. The result is very simple, to linear order in $Z$:

$Q_{new}^{00}=Q^{00}-\frac{1}{2}({\rm Tr}(Q^{20}Z^{*})+\textrm{h.c.})+{\cal
O}(Z^{2}).$ (13)

The same formula applies to the Hamiltonian as well,

$H_{new}^{00}=H^{00}-\frac{1}{2}({\rm Tr}(H^{20}Z^{*})+\textrm{h.c.})+{\cal
O}(Z^{2}).$ (14)

From these formulas it is apparent that the derivative of the expectation
value with respect to the variables $z_{ij}^{*}$ in $Z^{*}$ is111 The
derivative is taken with respect to the variables in the skew-symmetric
$Z^{*}$, ie. $z_{ji}^{*}=-z_{ij}^{*}$ and $z_{ij}$, $z_{ij}^{*}$ are treated
as independent variables.

${\partial\over\partial z_{ij}^{*}}Q^{00}=Q^{20}_{ij}.$ (15)

With a formula for the gradient of the quantity to be minimized, we have many
numerical tools at our disposal to carry out the minimization.

It is quite straightforward to introduce constraining fields in the
minimization process. As seen in Eq. (13) the transformation $Z$ will not
change the expectation value of $\hat{Q}$ to linear order provided ${\rm
Tr}{(Q^{20}Z^{*})}+\textrm{h.c.}=0$. Thus, one can change the configuration
without affecting the constraint (to linear order) by projecting $Z$ to
$Z_{c}$ as $Z_{c}=Z-\lambda Q^{20}$ with $\lambda=\frac{1}{2}({\rm
Tr}(Q^{20}Z^{*})+\textrm{h.c.})/{\rm Tr}(Q^{20}Q^{20\,*})$. With multiple
constraints, the projection has the form

$Z_{c}=Z-\sum_{\alpha}\lambda_{\alpha}Q^{20}_{\alpha}.$ (16)

The parameters $\lambda_{\alpha}$ are determined by solving the system of
linear equations,

$\sum_{\alpha}M_{\alpha\beta}\lambda_{\alpha}=\frac{1}{2}({\rm
Tr}(Q^{20}_{\beta}Z^{*})+\textrm{h.c.})$ (17)

where $M_{\alpha\beta}={\rm Tr}(Q^{20}_{\alpha}Q^{20\,*}_{\beta})$. Since we
want to minimize the energy, an obvious choice for the unprojected $Z$ is the
gradient of the Hamiltonian $H^{20}$. In this case the constraining parameters
$\lambda_{\alpha}$ are identical to the Lagrange multipliers in the usual HFB
equations. We will use the notation $H_{c}$ for the constrained Hamiltonian

$H_{c}=H-\sum_{\alpha}\lambda_{\alpha}Q_{\alpha}.$ (18)

###  II.1 Numerical aspects of the minimization

The most obvious way to apply the gradient method is to take the direction for
the change from Eq. (16,17), and take the length of the step as an adjustable
numerical parameter. We will call this the fixed gradient (FG) method. It is
implemented in the program as

$Z_{\eta}=\eta H_{c}^{20}.$ (19)

Typically the starting $U,V$ configuration will not satisfy the constraints,
and the $Z$ transformations must also bring the expectation values of the
operators to their target values $q_{\alpha}$. The error vector $\delta
q_{\alpha}$ to be reduced to zero is given by

$\delta q_{\alpha}=Q^{00}_{\alpha}-q_{\alpha}.$ (20)

We apply Eq. (13) to first order to obtain the desired transformation
$Z_{\delta q}$,

$Z_{\delta q}=-\sum_{\alpha\beta}M^{-1}_{\alpha\beta}\delta
q_{\alpha}Q^{20}_{\beta}.$ (21)

With these elements in hand, a new configuration is computed using the
transformation

$Z=Z_{c}+Z_{\delta q}.$ (22)

This process is continued until some criterion for convergence is achieved. We
shall measure the convergence by the norm of the gradient $|H_{c}^{20}|$. This
is calculated as

$|H_{c}^{20}|=\left({\rm Tr}[H_{c}^{20}(H_{c}^{20})^{\dagger}]\right)^{1/2}.$
(23)

An example using this method as given is shown in Fig. 1.

Figure 1:  Number of iterations required for convergence using Eq. (19) and
fixed $\eta$. At the point $\eta=0.12$ MeV-1 and beyond, the iteration process
is unstable. The converged solutions and their energies are the same for all
values of $\eta$ shown in the plot. All values producing converged solutions
The system is 24Mg with three constraints, $N$, $Z$, and $<Q_{Q}>=10$
$\hbar/m\omega_{0}$. The convergence criterion is $|H^{20}_{c}|<1.0\times
10^{-2}$ MeV. See Section VII.2 for further details.

The parameter $\eta$ is fixed to some value and the iterations are carried out
until convergence or some upper limit is reached. The required number of
iterations varies roughly inversely with $\eta$, up to some point where the
process is unable to find a minimum in a reasonable number of iterations.

There are a number of ways to speed up the iteration process. If the
constraints are satisfied, the parameter $\eta$ can be increased considerably.
Fig. 2 shows the change in $H^{00}_{c}$ from one iteration cycle as a function
of $\eta$ using $Z_{c}$ to update.

Figure 2:  Single-step energy change as a function of $\eta$ in Eq. (19). The
configuration that was updated is the 10th iteration step of the system in
Fig. 1.

For small values of $\eta$, the change in constrained energy is given by the
Taylor expansion Eq. (14), $\Delta H^{00}_{c}\approx-\eta
Tr\left(H_{c}^{00\,*}H_{c}^{00}\right)$. This function is shown as the
straight line in the Figure. The actual change is shown by the black circles.
One sees that $\eta$ could be doubled or tripled from the maximum value
permitted in Fig. 1. However, the constraints and other aspects of the new
$U,V$ become degraded so that such steps are not permissible for many
iterations re82 . Still, one can take advantage of the possible improvement by
choosing $\eta$ at each iteration taking account of the relevant information
from the previous iteration. This can be extracted from the ratio

$r=\frac{\Delta H^{00}_{c}}{\eta Tr\left(H_{c}^{00\,*}H_{c}^{00}\right)}$ (24)

which is close to one for too-small $\eta$ values and close to $\frac{1}{2}$
at the value corresponding to the steepest-descent minimum. We call such
methods variable gradient. We note that updates with $Z_{\delta q}$ alone are
relatively quick because there is no need to evaluation matrix elements of the
Hamiltonian. These considerations are implemented in the code of Ref. wa02  by
interspersing cycles of iteration by $Z_{\delta q}$ alone among the cycles
with updates by Eq. (22).

Another way to improve the efficiency of the iteration process is to divide
the elements of $H^{20}_{c}$ by preconditioning factors $p_{ij}$,

$(Z_{c})_{ij}=\eta{(H_{c}^{20})_{ij}\over p_{ij}}.$ (25)

The choice of the preconditioner is motivated by Newton’s method to find zeros
of a function (here $H_{c}^{20}$) based on knowledge of its derivative. This
could be accessible from the second-order term in Eq. (14), but unfortunately
it cannot be easily computed as it involves the HFB stability matrix. However
a reasonable approximation to it can be obtained from $H_{c}^{11}$, the one-
quasiparticle Hamiltonian that, when in diagonal form, is the dominant
component of the diagonal of the stability matrix. One first transforms $U,V$
to a basis that diagonalizes $H_{c}^{11}$. Call the eigenvalues of the matrix
$E_{i}$ and the transformation to diagonalize it $C$. The $U,V$ are
transformed to $U^{\prime},V^{\prime}$ in the diagonal quasiparticle basis by

$U^{\prime}=UC;\,\,\,\,\,V^{\prime}=V^{\prime}C$ (26)

In the new basis the preconditioner is given by

$p_{ij}=\max(E_{i}+E_{j},E_{min})$ (27)

where $E_{min}$ is a numerical parameter of the order of 1-2 MeV. The main
effect of the preconditioner is to damp away those components of the gradient
with high curvatures (i.e. second derivatives) which correspond to two-
quasiparticle excitations with large excitation energies. This is very
important for Hamiltonians that have a large range of single-particle
energies, such as the ones derived from commonly used nuclear energy density
functionals such as Skyrme and Gogny.

In Table I we show the number of iterations required to reach convergence for
a case calculated in Table II, to be described below.

Method | $\eta$ | $\eta_{min}$ | $\eta_{max}$ | $I_{conv}$
---|---|---|---|---
fixed gradient | 0.10 MeV-1 |  |  | 140
variable gradient |  | 0.08 MeV-1 | 0.3 MeV-1 | 65
fixed pr. | 0.7 |  |  | 72
variable pr. |  | 0.7 | 2.0 | 34
Table 1:  Number of iterations to convergence $I_{conv}$ with various
treatments of the update. Eq. (19) with fixed and variable gradients is used
for the top two lines and the preconditioned gradients Eq. (25) are used for
the lower two lines. The system is 21Ne as calculated in the top first entry
in Table II.

We see that there is a gain of more than a factor of 3 between the naive
steepest descent and the preconditioned gradient with a variable $\eta$.
Similar ideas have been used in a HF context in re82 ; um85  with similar
speedups.

##  III Odd-A nuclei

As discussed by Ring and SchuckRS , each $U,V$ set can be characterized by its
number parity, either even or odd. This means that when the wave function is
constructed and states of definite particle number are projected out, the
nonzero components will have either all even or all odd particle number.
Another important fact is that the generalized Thouless transformation does
not change the number parity of the Bogoliubov transformation. Thus, if we
start from a $U,V$ set of odd number parity, the final converged configuration
will only have components of odd nucleon number.

In fact, in the matrix-diagonalization method of solving the HFB equations,
the higher energy of the odd-A configurations requires some modification to
the Hamiltonian or to the iteration process. A common solution is to add
additional constraining fields so the that odd-A system has lower energyba73 ;
be09 . Typically the external field to be added breaks time reversal symmetry
in some way. But then one can no longer assert that a true minimum has been
found, because the extra constraints can affect the configuration. The
gradient method does not have this shortcoming. If the space of odd-number
parity Bogoliubov transformations is adequately sampled, it will find the
global minimum of the odd-A configurations. Moreover, with the gradient method
one does not need to modify the computer code to treat odd-$A$ systems. Only
the initial $U,V$ set is different for the two cases.

We note the $H^{11}_{c}$ has negative quasiparticle eigenenergies in the odd
number-parity space, assuming that the true minimum of the HFB functional is
an even number-parity configuration.

##  IV Other special cases

The variational minimum might not be directly reachable by the generalized
Thouless transformation, but it always is a limit of a succession of
transformations. This is the case if the condensate vanishes at the minimum
while the starting configuration has a finite condensate. This does not cause
any practical difficulties except for reducing the rate of convergence. Still,
in such cases it is more direct to start with a $U,V$ configuration of the
pure Hartree-Fock form. It is not possible to use the gradient method in the
other direction, to go to a minimum having a finite condensate from a starting
$U,V$ of Hartree-Fock form, as explained below.

##  V Imposed symmetries

The $U,V$ matrices have a dimension of the size of the Fock space of nucleon
orbitals and in principle can be dense matrices. However, one often imposes
symmetries on the wave function by assuming that the $U,V$ have a block
structure with all elements zero outside the blocks. For example, most codes
assume separate blocks for neutrons and protons. This is well-justified when
there is a significant difference in neutron and proton numbers but in general
it is better to allow them to mix. Other quantum numbers that are commonly
imposed on the orbital wave functions are parity and axial symmetry. There are
only a few exceptional nuclei that have HFB ground states breaking these
symmetries. For the parity, there are the Ra nuclei and Th nuclei. Concerning
axial symmetry, a global study of even-even nuclei with the Gogny functional
de10  found only three cases of nonaxial HFB minima among 1712 nuclei.

The number of orthogonal minima that can be easily calculated in the gradient
method depends on the assumed block structure. In the even number-parity space
there is just one global minimum. But in the odd number-parity space the
number parity of each block is conserved in the iteration process, so there
will be one state for each block. For example, states of different $K$-quantum
number may be calculated by imposing a block structure that imposes axial
symmetry. Thus for odd-A nuclei, the quasiparticle can be in any of the
$K$-blocks, giving a spectrum of states with $K$ specified by the block.

A more subtle form of possible imposed symmetries is those contained in the
starting $U,V$ configuration. The energy $H^{00}$ is essentially a quadratic
function of symmetry-breaking densities because the products of densities in
the functional must respect the symmetries of the Hamiltonian. If these
components are zero in the initial configuration, the energy is stationary at
that point and there is no gradient to generate nonzero field values. The
typical cases are quadrupole deformation in the ordinary density and any form
of anomalous densities. Fortunately, it is very easy to avoid unwanted
symmetries in the starting $U,V$ as discussed below.

##  VI The code hfb_shell

The code hfb_shell presented in this paper is described in more detail in the
Appendix. The main point we want emphasize about the code is that it is
organized in modules that separate out the functions that are independent of
the Hamiltonian from those that are specific to it. Also, the block structure
is specified only by the code input, and can easily be changed. The examples
we show are for the $sd$-shell using the USDB Hamiltonian br06 . Since that
Hamiltonian is specified by the fitted numerical values of the 3 single-
particle energies and the 63 $JT$-coupled two-particle interaction energies,
it does not have any symmetries beyond those demanded by the physics. In
particular, the HFB fields obtained with it should provide a realistic
description of aspects such as the time-odd fields, that are difficult to
assess with the commonly used energy functionals such as those in the Skyrme
family.

###  VI.1 Application to the $sd$-shell

The $sd$ shell-model space has a dimension of 24 and the principal matrices
$U,V,Z,...$ have the same dimension. In the application presented here, we
assume axial symmetry which splits the matrices in blocks of dimension 12, 8
and 4 for $m$-quantum numbers $\pm\frac{1}{2}$, $\pm\frac{3}{2}$, and
$\pm\frac{5}{2}$ respectively. Neutron and proton orbitals are in the same
blocks, so the basis is sufficiently general to exhibit neutron-proton
pairing, if that is energetically favorable. We also assume that the matrices
are real.

We often start with a $U,V$ configuration of canonical form, namely $U$
diagonal, $U_{ij}=u\delta_{ij}$. The nonzero entries of the $V$ are all equal
to $\pm v=\pm\sqrt{1-u^{2}}$, and are in positions corresponding to pairing in
the neutron-neutron channel and the proton-proton channel. We arbitrarily take
$u=0.8$ and $v=0.6$ for the starting configuration $U_{0},V_{0}$. This may be
modified in a number of ways before it is used as a starting configuration in
the gradient minimization. When calculating a nucleus for which $N$ or $Z$ is
zero or 12, it is more efficient to use $U,V$ matrices that have those
orbitals empty or completed filled in the starting configuration. This is
carried out by changing $u,v$ to zero or one for the appropriate orbitals. The
particle number of that species is then fixed and is not constrained in the
gradient search.

For odd-number parity configurations, the $U,V$ is changed in the usual way by
interchanging a column in the $U$ matrix with the corresponding column in $V$.
The space that will be searched in the gradient method then depends on the
block where the interchange was made. In principle it does not depend on which
column of the block was changed. However, there is some subtlety is making use
of this independence which will be discussed below.

We may also apply a random $Z$ transformation to the starting configurations.
Since all the entries in the upper triangle of the $Z$ matrix are independent,
we can populate them with random numbers. This seems to be a good way to break
unwanted symmetries in the starting configuration that would be preserved by
the gradient update. We denote by $U_{r},V_{r}$ the configuration generated
from $U_{0},V_{0}$ by a randomly generated $Z$.

In principle one could also start from the $U,V$ configuration of the vacuum:
$U=1,V=0$. We have tried this and found, as might be expected, that the
proportion of false minima is larger than is obtained with $U_{0},V_{0}$.

##  VII Three examples

In this section we will describe the HFB calculations for three nuclei, 32Mg,
24Mg, and 21Ne. The first one is typical of a spherical nucleus that exhibits
identical-particle pairing. The second is a well-deformed nucleus. The third
illustrates the method for an odd-A system.

For calculating matrix elements of the quadrupole operator $Q_{Q}$, we will
treat the single-particle wave functions as harmonic oscillator functions of
frequency $\omega_{0}$, and report the quadrupole moments in units of
$\hbar/m\omega_{0}$.

###  VII.1 32Mg

The nucleus 32Mg ($(N,Z)=(12,4)$ in the $sd$-shell) behaves as expected of a
semimagic nucleus in HFB. Please note that we do not include in our
configuration space the $f_{7/2}$ intruder shell required to explain the
deformation properties of this nucleus mo95 ; ro02 . We calculate the HFB
ground state in two ways, illustrating the role of the starting configuration.
The first is to use a randomized $U_{r},V_{r}$ configuration and constraining
the particle numbers to the above values. Another way is to start with a
prolate configuration similar to $U_{0},V_{0}$ for the protons and with all
the neutron orbitals filled. In that case, only the proton number is
constrained. Both iteration sets converge to the same minimum, a spherical
configuration having a strong proton pairing condensate. The output
characteristics are $E_{HFB}=-135.641$ MeV, $Q_{Q}^{00}=0.00$ and $\Delta
Z^{2}=2.93$. The zero value for $Q_{Q}^{00}$ shows that the configuration is
spherical, and the nonzero value for $\Delta Z^{2}$ shows that protons are in
a condensate. Next we calculate the condensation energy, defined as the
difference between $E_{HFB}$ and the Hartree-Fock minimum $E_{HF}$. The
easiest way to find the HF minimum is to repeat the calculation with an
additional constraint that forces the condensate to zero. This is done by
adding a $G$-type operator that is sensitive to the presence of a condensate.
Carrying this out, we find a minimum at $E_{HF}=-134.460$ MeV and
$Q_{Q}^{00}=5.08$. The extracted correlation energy is $E_{HF}-E_{HFB}=1.18$
MeV, which is much smaller than what one would obtain with schematic
Hamiltonians fitted to pairing gap. It is also interesting to extract the
quasiparticle energies, since they provide the BCS measure of the odd-even
mass differences. These are obtained by diagonalizing $H_{c}^{11}$. The
results for the HFB ground state range from 1.5 to 9 MeV, with the lowest
giving the BCS estimate of the pairing gap.

###  VII.2 24Mg

The next nucleus we consider, 24Mg with $N=4$ and $Z=4$, is strongly deformed
in the HFB ground state. We find that the converged minimum has a quadrupole
moment $\langle Q_{Q}\rangle=12.8$, close to the maximum allowed in the space.
More surprisingly, the pairing condensate vanishes at the HFB convergence. We
now make a set of constrained calculations to display the energy as a function
of quadrupole moment. The starting configuration is generated by applying a
random transformation to $U_{0},V_{0}$. The gradient code carries out the
iterations with the constraints $N=4$, $Z=4$, and the chosen value of $Q$. The
convergence of the constraints to their target values is very rapid, using the
update in Eq. (21). This is illustrated in Fig. 3, showing the

Figure 3:  Error in constrained quantities as a function of iteration number
for the $\eta=0.1$ run of the 24Mg iterations in Fig. 1. Quantities
constrained are: $N$, open circles; $Z$, filled squares; and $Q_{Q}$, filled
circles.

deviation from the target values as a function of iteration number in one of
the cases ($Q=10$). On the other hand, the convergence to the minimum of the
HFB energy can be slow, using a fixed-$\eta$ update with Eq. (19). The
calculations were carried out setting the convergence criterion
$|H_{c}^{20}|<0.01$ MeV. Fig. 4 shows the number of iterations required to
reach convergence for the various deformations.

Figure 4:  Number of iterations required to convergence for the calculated
configurations on the deformation energy curve Fig. 5.

They range from $~{}40$ to $~{}250$. In a number of cases, the iterations seem
to be approaching convergence, but the system is actually in a long valley,
and eventually a lower minimum is found. It may also happen that the gradient
method finds a local minimum that is not the global one. Perhaps 10% of the
runs end at a false minimum. This can often be recognized when carrying
constrained calculations for a range of constraint values, as it gives rise to
discontinuities in the energy curves. The only systematic way we have to deal
with the false minima is to run the searches with different randomly generated
starting configurations, and select the case that gives the lowest energy. The
resulting deformation plot combining two runs is shown in Fig. 5.

Figure 5:  HFB energies as a function of deformation, using the $Q_{Q}$
quadrupole constraint. The nucleus is 24Mg, $N=Z=4$ in the $sd$-shell.

The global minimum is at a large prolate deformation as mentioned earlier.
There is also a secondary minimum at a large oblate deformation. For all
deformations, the ordinary neutron-neutron and proton-proton pairing
condensates are small or vanish.

###  VII.3 21Ne

The next nucleus we discuss, 21Ne with $(N,Z)_{sd}=(3,2)$, illustrates how the
gradient method makes use of the conserved number parity to find the minimum
of odd-A systems. We start with the $U_{0},V_{0}$ configuration, and convert
it to an odd-number parity configuration by exchanging two columns in the
$m=\pm\frac{1}{2}$ block. There are 6 possible columns with $m=+\frac{1}{2}$
that can be exchanged. The results for the converged energies are shown in the
top row of Table 2. All of the neutron exchanges give the same final energy,
$-40.837$ MeV. However, the energy is different for proton exchanges. The
reason is that the starting configurations do not mix neutrons and protons,
and for reasons discussed earlier the corresponding gradients are zero. This
unwanted symmetry can be broken by making a random transformation of the
initial configuration. The results are shown in the second row. Now all the
energies are equal, showing that the minimum can be accessed from any column
exchange. Interestingly, the energy is lower than in the previous set of
minimizations. This shows that there is a significant neutron-proton mixing in
the condensate for 21Ne.

$U,V$ | $d^{n}_{5/2,1/2}$ | $d^{n}_{3/2,1/2}$ | $s^{n}_{1/2,1/2}$ | $d^{p}_{5/2,1/2}$ | $d^{p}_{3/2,1/2}$ | $s^{p}_{1/2,1/2}$
---|---|---|---|---|---|---
$U_{0},V_{0}$ | -40.837 | -40.837 | -40.837 | -40.215 | -40.176 | -40.176
$U_{r},V_{r}$ | -41.715 | -41.715 | -41.715 | -41.715 | -41.715 | -41.715
Table 2: HFB energies of 21Ne, with different starting configurations. For the
top row, the starting configuration is $U_{0},V_{0}$ with the indicated column
in the $m=\pm\frac{1}{2}$ block interchanged. The second row starts from a
randomized configuration $U_{r},V_{r}$ as discussed in Sect. VI.1.

## Acknowledgments

The authors thank A. Gezerlis and P. Ring for discussions, T. Lesinski and J.
Dobaczewski for comments on the manuscript, and M. Forbes for comments on the
code. This work (GFB) was supported in part by the U.S. Department of Energy
under Grant DE-FG02-00ER41132, and by the National Science Foundation under
Grant PHY-0835543. The work of LMR was supported by MICINN (Spain) under
grants Nos. FPA2009-08958, and FIS2009-07277, as well as by Consolider-Ingenio
2010 Programs CPAN CSD2007-00042 and MULTIDARK CSD2009-00064.

## References

  * (1) P. Ring and P. Schuck, The nuclear many-body problem, (Springer, 1980).
  * (2) P.-G. Reinhard and R.Y. Cusson, Nucl. Phys. A378 418 (1982).
  * (3) K.T.R. Davies, H. Flocard, S. Krieger, and M.S. Weiss, Nucl. Phys. A342 111 (1980).
  * (4) P. Bonche, H. Flocard, and P.-H. Heenen, Comput. Phys. Commun. 171 49 (2005).
  * (5) J. Dobaczewski, and P. Olbratowski, Comput. Phys. Commun. 167 214 (2005).
  * (6) K. Bennaceur and J. Dobaczewski, Comp. Phys. Commun. 168 96 (2005)
  * (7) W. Pöschl, D. Vretenar, A. Rummel, and P. Ring, Comput. Phys. Commun. 101 75 (1997).
  * (8) M. Stoitsov, et al., Comput. Phys. Commun. 167 43 (2005).
  * (9) J.L. Egido, J. Lessing, V. Martin, and L.M. Robledo Nucl. Phys. A594 70 (1995)
  * (10) M. Warda, J.L. Egido, L.M. Robledo, and K. Pomorski, Phys. Rev. C 66 014310 (2002).
  * (11) I. Maqbool, J.A. Seikh, P.A. Ganai, and P.Ring, J. Phys. G: Nucl. Part. Phys. 38 045101 (2011).
  * (12) R. Rodríguez-Guzmán, Y. Alhassid, and G.F. Bertsch, Phys. Rev C77, 064308 (2008)
  * (13) L.M. Robledo and G.F. Bertsch, in preparation.
  * (14) A. Gezerlis, G.F. Bertsch, and L. Luo, arXiv:1103.5793 (2011).
  * (15) A.S. Umar, et al. Phys. Rev. C32 172 (1985).
  * (16) G.F. Bertsch, J. Dobaczewski, W. Nazarewicz, and J. Pei, Phys. Rev. A 79 043662 (2009).
  * (17) B. Banerjee, P. Ring, and H.J. Mang, Nucl. Phys. A 215 266 (1973).
  * (18) J.-P. Delaroche, et al., Phys. Rev. C 81 014303 (2010).
  * (19) B.A. Brown and W.A. Richter, Phys. Rev. C 74 034315 (2006).
  * (20) T. Motobayashi, et al., Phys. Lett. B346 9 (1995).
  * (21) R.Rodríguez-Guzmán, J.L. Egido, and L.M. Robledo Nucl. Phys. A709 201 (2002).
  * (22) K. J. Millman and M. Aivazis, Comp. Sci. Eng. 13 9 (2011)

## Appendix: explanation of the code

The code hfb_shell that accompanies this article implements the gradient
method discussed in the text222The code may be downloaded from
http://www.phys.washington.edu/users/bertsch/hfb-shell.21.tar until it has
been published in a journal repository.. The code is written in Python and
requires the Python numerical library numpy to run (see py11  and accompanying
papers for a description of Python in a scientific environment). The main
program is the file hfb.py. It first carries out the initialization using
information from the primary input data file that in turn contains links to
other needed data files. There are three of these, one for the Hamiltonian
parameters, one for the correspondence between orbitals and rows of the $U,V$
matrices include the assumed block structure, and one for the input $U,V$
configuration. The input data format is explained in the readme.txt of the
code distribution.

Following initialization, program enters the iteration loop, calling the
various functions used to carry out the iteration. The loop terminates when
either a maximum number of iterations itmax is reached or the convergence
parameter $|H^{20}_{c}|$ go below a set value converge.

The function calls that are specific to the $sd$-shell application are
collected in the module sd_specific.py. The tasks carried out by these
functions include:

  * •

initialization of matrix sizes and block structures

  * •

setting up the matrices representing single-particle operators in the shell-
model basis.

  * •

calculation of the fields $\Gamma,\Delta$ from the densities $\rho,\kappa$.
This function makes use of a table of interaction matrix elements $v_{ijkl}$
that are read in from a file. The present distribution of the code only
provides the Hamiltonian data for the USDB interaction br06 .

The functions that are generic to the gradient method are collected in the
module hfb_utilities.py. Many of these functions are defined by equations in
the text; the correspondence is given in Table III.

Function call | Equation in text
---|---
rho_kappa | (1)
F20 | (6)
G20 | (8)
H20 | (11)
H00 | (3)
Ztransform | (12)
Table 3:  Python functions in hfb_utilities.py corresponding to equations in
the text.

The output of hfb.py reports the expectation values of the Hamiltonian and the
single-particle operators $N,Z$ and $Q_{Q}$ at each iteration step, together
with the convergence parameter $|H^{20}_{c}|$. After the final iteration, the
values are reported for the expectation values of constraining parameters
$\lambda_{\alpha}$ and the number fluctuations $\Delta N^{2},\Delta Z^{2}$.
The final $U,V$ configuration is written to the file uv.out. Thus additional
iterations can be performed simply by specifying uv.out as the new input file.

In addition, there is a set of functions collected in the module hfb_tools.py.
These are useful for making input $U,V$ configurations and for analyzing the
output $U,V$ configuration, but are not needed to run hfb.py. For example, a
randomizing transformation can be applied to a $U,V$ configuration by the
function randomize. Another useful function is canonical, used to extract the
eigenvalues of the $\rho$ operator needed for the canonical representation.


# Quantifying Unknown Quantum Entanglement via a Hybrid Quantum-Classical
Machine Learning Framework

Xiaodie Lin1 Zhenyu Chen1 Zhaohui Wei2,3, Email<EMAIL_ADDRESS>1Institute for Interdisciplinary Information Sciences, Tsinghua University,
Beijing 100084, China
2Yau Mathematical Sciences Center, Tsinghua University, Beijing 100084, China
3Yanqi Lake Beijing Institute of Mathematical Sciences and Applications,
101407, China

###### Abstract

Quantifying unknown quantum entanglement experimentally is a difficult task,
but also becomes more and more necessary because of the fast development of
quantum engineering. Machine learning provides practical solutions to this
fundamental problem, where one has to train a proper machine learning model to
predict entanglement measures of unknown quantum states based on
experimentally measurable data, say moments or correlation data produced by
local measurements. In this paper, we compare the performance of these two
different machine learning approaches systematically. Particularly, we first
show that the approach based on moments enjoys a remarkable advantage over
that based on correlation data, though the cost of measuring moments is much
higher. Next, since correlation data is much easier to obtain experimentally,
we try to better its performance by proposing a hybrid quantum-classical
machine learning framework for this problem, where the key is to train optimal
local measurements to generate more informative correlation data. Our
numerical simulations show that the new framework brings us comparable
performance with the approach based on moments to quantify unknown
entanglement. Our work implies that it is already practical to fulfill such
tasks on near-term quantum devices.

††preprint: APS/123-QED

##  I Introduction

Quantum entanglement is a crucial resource for many quantum schemes and
quantum protocols, such as quantum superdense coding Bennett and Wiesner
(1992), quantum teleportation Bennett  _et al._ (1993) and quantum
cryptography Ekert (1991). As a result, numerous measures have been raised to
quantify the amount of entanglement contained in quantum states Horodecki  _et
al._ (2009); Amico  _et al._ (2008). For bipartite pure states, entanglement
measure is uniquely defined by the von Neumann entropy of subsystems. However,
the landscape is far more complex for bipartite mixed states, where many
important questions on entanglement quantifications have not been answered. An
even more complicated case is the quantifications of multipartite quantum
entanglement, for which many different measures have been proposed but most of
them are very hard to calculate. It can be said that quantum entanglement has
not been understood well theoretically.

Meanwhile, in recent years many subareas of quantum engineering are under fast
development, and because of this, detecting and even quantifying _unknown
entanglement_ involved in quantum experiments is becoming more and more
realistic and necessary, though this is also a difficult task to fulfill. For
example, in a future factory that produces maximally entangled states by
entanglement distillations, an efficient technique that can roughly quantify
unknown entanglement will be helpful to select high quality raw materials. We
stress that unknown entanglement means that we do not have any prior knowledge
of target quantum states, and therefore we have to characterize them and
particularly quantify their entanglement by performing quantum measurements.

Existing quantification methods for unknown entanglement can be mainly divided
into four branches. First, quantum tomography is the most popular approach
adopted by quantum experimentalists when the size of target quantum states is
small Chuang and Nielsen (1997); Poyatos  _et al._ (1997), where
reconstructing quantum states allows us to look into the underlying
entanglement. However, quantum tomography requires exponential cost, which is
unbearable in high dimensional quantum systems. Second, recently a new
technique of quantifying entanglement, say estimating the Rényi entropies, has
been proposed, which performs random measurements on quantum states and then
analyzes the outcome statistics Van Enk and Beenakker (2012); Brydges  _et
al._ (2019). Nevertheless, the cost of this technique is very high, where a
large number of measurement settings are necessary. Third, device-independent
protocols have been proposed to lower bound various entanglement measures
Moroder  _et al._ (2013); Wei and Lin (2021); Lin and Wei (2021). These
device-independent methods quantify entanglement exclusively from the observed
measurement statistics on subsystems, thus independent of any assumptions on
the interested quantum systems. In the realm of noisy intermediate scale
quantum (NISQ) era, device-independent protocols are attractive due to their
efficiency and reliability. However, such device-independent protocols usually
have limited applications in practical as they can provide nontrivial results
only when the observed quantum nonlocality is very strong. Fourth, in addition
to these analytical methods, machine learning has also been utilized to
provide practical solutions to quantify entanglement experimentally Gray  _et
al._ (2018); Lin  _et al._ (2021); Roik  _et al._ (2022). In such methods,
different experimentally accessible data on quantum states is collected and
fed into machine learning models such that the mappings from experimental data
to target entanglement measures are learned, by which one can predict the
entanglement measures of unknown quantum states. In the current paper, we will
focus on the fourth kind of methods.

In fact, according to the sorts of involved experimental data, there exist
already two different approaches to apply machine learning onto quantifying
unknown quantum entanglement experimentally Gray  _et al._ (2018); Lin  _et
al._ (2021). Specifically, in Ref.Gray  _et al._ (2018) the target
entanglement measure is negativity, and for this the moment data of partially
transposed states is fed into machine learning models as data features, which
is usually very costly to obtain Zhou  _et al._ (2020); Elben  _et al._
(2020). In Ref.Lin  _et al._ (2021), correlation data serves as data features
to quantify unknown entanglement, which is relatively convenient to prepare as
one only needs to measure subsystems of target quantum states with a small set
of local measurements chosen beforehand.

In this paper, we first show that several entanglement measures, like the
relative entropy of entanglement, can also be quantified accurately using
machine learning models based on the original moments of quantum state $\rho$
defined as $\mu_{m}(\rho)={\rm Tr}\left(\rho^{m}\right)$. Then we compare this
new approach with the one in Ref.Lin  _et al._ (2021), and show that it can
beat the latter easily in performance, though usually moment data are much
harder to extract than correlation data. For example, a large number of
measurement devices are already needed even if only $\mu_{2}(\rho)$ is
measured Vermersch  _et al._ (2018); Brydges  _et al._ (2019) .

Meanwhile, since correlation data comes from a small set of local measurements
and is much easier to collect, it will be nice if we can somehow improve the
performance of the machine learning approach based on correlation data. Note
that in the above comparisons, correlation data is generated by a certain
fixed set of local measurements for all training and test quantum states.
Therefore, a possible way to improve the performance is to choose better local
measurements for correlation data generations. For this, we propose a hybrid
quantum-classical machine learning framework to quantify unknown entanglement
based on correlation data, where optimal local measurements are trained to
generate correlation data. Our numerical simulations show that the new
framework allows the correlation method to achieve a comparable performance
with the machine learning approach based on moments in quantifying unknown
entanglement task. Since the cost to produce correlation data is reasonable,
our results show that it is already practical to quantify unknown entanglement
on near-term quantum devices.

##  II Settings and Entanglement Measures

Consider a bipartite state $\rho$ shared by separated two parties, Alice and
Bob. Alice (respectively, Bob) has a set of measurement devices labeled by $X$
(respectively, $Y$) to measure her (his) subsystem, and the possible
measurement outcomes are labeled by $A$ (respectively, $B$). After repeating
the measurement many times, Alice and Bob calculate the joint conditional
probabilities $p(ab|xy)$, which indicates the probability of obtaining
outcomes $(a,b)\in A\times B$ upon selecting measurement settings $(x,y)\in
X\times Y$. Suppose $\\{M_{x}^{a}\\}$ is the operators for the quantum
measurement performed by Alice’s measurement device $x\in X$, where $a\in A$,
and analogously for $\\{N_{y}^{b}\\}$, then it holds that

$\displaystyle p(ab|xy)={\rm Tr}\left[\left(M_{x}^{a}\otimes
N_{y}^{b}\right)\rho\right].$ (1)

A _correlation_ $p=[p(ab|xy)]$ is a vector contains all the joint conditional
probabilities of form $p(ab|xy)$.

We now turn to _moments_ of quantum states, which are defined as Ekert  _et
al._ (2002); Horodecki (2003)

$\displaystyle\mu_{m}(\rho)={\rm Tr}\left(\rho^{m}\right).$ (2)

Clearly, $\mu_{1}(\rho)={\rm Tr}(\rho)=1$ and $\mu_{2}(\rho)$ is the _purity_
of $\rho$. Experimentally, $\mu_{m}(\rho)$ can be measured directly by
performing joint measurements on $m$ copies of the same state $\rho$ Ekert
_et al._ (2002), while this operation is very hard with current quantum
technologies, especially when $m$ is large. To overcome this difficulty,
techniques that can estimate $\mu_{m}(\rho)$ on single-copy states have also
been developed Van Enk and Beenakker (2012).

A complete supervised machine learning system contains three ingredients: data
features, data labels and a machine learning model Mohri  _et al._ (2018);
Shalev-Shwartz and Ben-David (2014). The training data features with correct
labels are fed into the machine learning model, where the parameters contained
in the learning model are trained properly such that the resulting system can
predict labels of unknown test data precisely. In our problem, correlation
data or moments serve as data features of quantum states, data labels are the
values of target entanglement measures, and the machine learning model is
supposed to learn an unknown nonlinear relationship between data features and
data labels. According to the no free lunch theorem Wolpert and Macready
(1997), a small target error rate implies that we need a large set of
representative training data. However, it is well-known that for bipartite
mixed quantum states and multipartite quantum states, most entanglement
measures are extremely hard to calculate, which means that it is hard for us
to prepare correct labels for training data sets. According to the computation
hardness and the importance of entanglement measures, in this paper we choose
the coherent information and the relative entropy of entanglement as our
target measures to quantify unknown quantum entanglement.

Coherent information is a fundamental quantity that measures the capability of
transition of quantum information Schumacher and Nielsen (1996); Lloyd (1997).
For an arbitrary bipartite quantum state
$\rho\in\mathcal{H}_{A}\otimes\mathcal{H}_{B}$, its coherent information is
defined as

$\displaystyle I_{C}(\rho)=S(\rho_{A})-S(\rho),$ (3)

where $S(\rho)$ is the von Neumann entropy of $\rho$ and $\rho_{A}={\rm
Tr}_{B}(\rho)$ is the subsystem in $\mathcal{H}_{A}$. A crucial property of
the coherent information is that for any bipartite $\rho$, it holds that
Cornelio  _et al._ (2011)

$\displaystyle E_{F}(\rho)\geq E_{D}(\rho)\geq I_{C}(\rho),$ (4)

where $E_{F}(\rho)$ and $E_{D}(\rho)$ are two most important measures of
entanglement, the entanglement of formation and the entanglement of
distillation respectively. Therefore, a good estimation of $I_{C}(\rho)$
implies that we obtain a very nontrivial piece of information on the amount of
entanglement for $\rho$. Furthermore, it is not hard to see that the coherent
information is very easy to calculate, which means that if we pick it as the
data label, we can generate a large amount of training or test data at low
costs.

Another quantity we will utilize is the relative entropy of entanglement
Vedral (2002), defined as

$\displaystyle E_{R}(\rho)=\min_{\sigma\in{\rm
SEP}}S(\rho\|\sigma)=\min_{\sigma\in{\rm SEP}}{\rm
Tr}(\rho\log\rho-\rho\log\sigma),$ (5)

where ${\rm SEP}$ denotes the set of all separable states. Relative entropy of
entanglement has a good geometric interpretation, as $E_{R}(\rho)$ measures
certain “distance” between $\rho$ and the set of separable states. Based on a
machine learning model called _active learning_ , we can numerically compute
the relative entropy of entanglement for some particular quantum states with
high accuracy Hou  _et al._ (2020). At the same time, it turns out that the
relative entropy of entanglement satisfies that Vedral and Plenio (1998)

$\displaystyle E_{F}(\rho)\geq E_{R}(\rho)\geq E_{D}(\rho).$ (6)

Therefore, similar to $I_{C}(\rho)$, estimating $E_{R}(\rho)$ of an unknown
quantum state $\rho$ also helps us to obtain nontrivial information of
$E_{F}(\rho)$ and $E_{D}(\rho)$.

##  III Comparisons between machine learning models based on orrelation data
and moments

In this section, we first present machine learning models that take
correlation data and moments as data features respectively, then we compare
their performance in predicting both the coherent information and the relative
entropy of entanglement systematically. For convenience, we call the former
the correlation method, and the latter the moment method.

###  III.1 Quantifying the Coherent Information

We begin with quantifying the coherent information of arbitrary 3-dimensional
bipartite quantum states $\rho\in\mathcal{H}^{3}\otimes\mathcal{H}^{3}$. To
generate more representative training data, we try to evenly sample quantum
states according to the distribution of their coherent information. However,
due to the low efficiency of sampling quantum states with coherent information
less than -1.5, in fact we only sample quantum states with coherent
information ranging from -1.5 to $\log_{2}3$. Specifically, we divide this
range into 31 intervals of size 0.1, and in each of these intervals we sample
1,291 quantum states randomly. Eventually, totally 40,021 states are randomly
sampled to compose the set of quantum states for training.

After sampling the training quantum states, we generate the corresponding
training data features for the correlation method, which is achieved by
performing the local measurements that maximize the violation of the Collins-
Gisin-Linden-Masser-Popescu (CGLMP) inequality on the training quantum states
Collins  _et al._ (2002); Zohren and Gill (2008), and then record the outcome
statistics. More concretely, Alice’s measurement $A_{k}$ can be characterized
by the eigenvectors

$\displaystyle|r\rangle_{A_{k}}=\frac{1}{\sqrt{d}}\sum_{q=0}^{d-1}{\rm
exp}\left(\frac{2\pi i}{d}q(r-\alpha_{k})\right)|q\rangle_{A},$ (7)

and Bob’s measurement $B_{l}$ can be characterized by the eigenvectors

$\displaystyle|r\rangle_{B_{l}}=\frac{1}{\sqrt{d}}\sum_{q=0}^{d-1}{\rm
exp}\left(-\frac{2\pi i}{d}q(r-\beta_{l})\right)|q\rangle_{B},$ (8)

where $0\leq r\leq d-1$, $1\leq k,l\leq N$,
$\alpha_{k}=(k-1/2)/N,\beta_{l}=l/N$, and $N=|X|=|Y|$ Barrett  _et al._
(2006). In the current task, we let $d=3$ and $N=2$. That is, both Alice and
Bob have $2$ different measurement devices.

Meanwhile, to apply the moment method, we generate two different sets of
training data, which contain different orders of moments. Specifically, one
set has $\\{\mu_{2}(\rho_{A}),\mu_{2}(\rho)\\}$ as data features, and the
other has
$\\{\mu_{2}(\rho_{A}),\mu_{2}(\rho),\mu_{3}(\rho_{A}),\mu_{3}(\rho)\\}$ as
data features.

To test the performance of the above three cases (one for the correlation
method and two for the moment method), we evenly sample around 2,000 quantum
states, and then produce the corresponding test correlation data or moments of
these states.

During the training stage, each training data set is fed into a 4-hidden-layer
fully connected neural network (FNN) with 400, 200, 100 and 50 neurons in each
layer, respectively. After training, we test its performance with the
corresponding test set sampled above. The results of all the three cases are
shown in Fig.1. As we can see, the moment method that utilizes only the second
order of moment already beats the correlation method in this task, whose
corresponding mean squared errors (MSEs) are 0.0582 and 0.0035, respectively.
In some sense this is not surprising since the number of measurement devices
required by the moment method is much larger than that of the correlation
method, and thus the moment data is more informative.

Table 1: The structure and configuration details of the convolutional neural network (see Refs.LeCun  _et al._ (1998); Albawi  _et al._ (2017) for introductions to convolutional neural networks). For each ‘$\cdot/\cdot$’, the former parameter represents for $d=5$ and the latter parameter represents for $d=8,10$. Layers | Type | Neurons | Filters | Kernel size | Strides | Pool size
---|---|---|---|---|---|---
0-1 | Convolution2D | (None, 9/18, 9/18, 32) | 32 |  $2\times 2$/$3\times 3$ | $1\times 1$ | -
1-2 | Max-pooling2D | (None, 8/16, 8/16, 32) | - | - | $1\times 1$ |  $2\times 2$/$3\times 3$
2-3 | Convolution2D | (None, 7/14, 7/14, 64) | 64 |  $2\times 2$/$3\times 3$ | $1\times 1$ | -
3-4 | Max-pooling2D | (None, 6/12, 6/12, 64) | - | - | $1\times 1$ |  $2\times 2$/$3\times 3$
4-5 | Convolution2D | (None, 5/10, 5/10, 64) | 64 |  $2\times 2$/$3\times 3$ | $1\times 1$ | -
5-6 | Fully-connected | (None, 64) | - | - | - | -
6-7 | Fully-connected | (None, 32) | - | - | - | -
7-8 | Fully-connected | (None, 1) | - | - | - | -

In addition, if we strengthen the moment method by also factoring in the third
order of moment, apparent further improvements can be observed, where the MSE
decreases to 0.0004, implying an excellent performance. As a result, compared
with the correlation method, the moment method enjoys a remarkable advantage
in coherent information quantification tasks.

Figure 1: The neural network predictions for $I_{C}(\rho)$. The blue stars
represent the prediction values of the correlation method and the
corresponding MSE is 0.0582. The green dots and orange diamonds represent the
prediction values of the moment method with $m=2$ and $m=2,3$, whose MSE are
$0.0035$ and $0.0004$, respectively.

We now increase the dimension of target quantum states
$\rho\in\mathcal{H}^{d}\otimes\mathcal{H}^{d}$ to $d=5,8,10$, and still aim to
quantify the coherent information. In this case, the quantum state sampling
procedure is similar to the dimension 3 case. Firstly, around 40,000 quantum
states are evenly sampled for training according to the distribution of
coherent information, and then around 2,000 quantum states are sampled for
test similarly. The generations of data features for these quantum states are
exactly the same as the previous task.

The training models for the moment method keep unchanged, i.e., a 4-hidden-
layer FNN with 400, 200, 100 and 50 neurons in each hidden layer is adopted.
While for the correlation method, since convolutional neural networks (CNNs)
behave better than FNNs when dimension increases, in the current task we
utilize CNNs to predict the coherent information instead of FNNs. Table 1
describes the structure and configuration details of our CNN model.

Applying the well-traind models on the test data sets we choose, the
experimental results for each case are listed in Table 2.

Table 2: The MSEs of predicting the coherent information for different methods and dimensions. Dimension | $d=3$ | $d=5$ | $d=8$ | $d=10$
---|---|---|---|---
Correlation, $N=2$ | 0.0582 | 0.0362 | 0.0212 | 0.0146
Moments, $m=2$ | 0.0035 | 0.0034 | 0.0025 | 0.0026
Moments, $m=2,3$ | 0.0004 | 0.0010 | 0.0011 | 0.0010

Similar to the dimension 3 case, there still exists an obvious gap between the
performance of the correlation method and that of the moment method.
Furthermore, it is interesting to see that the improvements of introducing the
third order of moment ($m=2,3$) over only using the second order ($m=2$)
decreases as the quantum dimension goes up. Therefore, considering the
experimental difficulty of measuring moments, it is a good choice to set $m=2$
when predicting coherent information for high dimensional quantum states.

###  III.2 Quantifying the Relative Entropy of Entanglement

Recall that the concept of the relative entropy of entanglement is defined as
an optimization problem over the set of separable states. Despite being
convex, the set of separable states is still very hard to fully characterize
Lu  _et al._ (2018); Doherty  _et al._ (2002). In Ref.Hou  _et al._ (2020), a
technique to upper bound the relative entropy of entanglement was proposed
based on active learning. Even though this method only provides upper bounds,
Ref.Hou  _et al._ (2020) demonstrated that these upper bounds are quite tight
for many quantum states, such as Werner states, isotropic states and random
bipartite quantum states with low dimensions. However, it should be stressed
that the method introduced in Ref.Hou  _et al._ (2020) requires full
descriptions of quantum states, i.e., the density matrix of the state, while
our mission is quantifying unknown entanglement based on experimentally
measuring quantities. In our task, we only utilize the active learning method
to provide labels for our data sets, due to its high accuracy.

As usual, to generate representative training data, we would like to sample
quantum states with evenly distributed relative entropy of entanglement.
However, when the target values of relative entropy of entanglement are high,
the sampling efficiency is very low. In addition, to provide correct labels
for the sampling quantum states, we have to calculate their relative entropy
of entanglement using the active learning method introduced above, which is
very costly. Both of these two facts make it challenging for us to generate
proper training and test data.

Here to prepare training data we focus on quantum states with the form

$\displaystyle\rho=(1-\epsilon)\rho_{0}+\epsilon|\psi_{+}\rangle\langle\psi_{+}|,$
(9)

where $\epsilon\in[0,1]$, $\rho_{0}$ is a random quantum state in
$\mathcal{H}^{d}\otimes\mathcal{H}^{d}$ and
$|\psi_{+}\rangle=\frac{1}{\sqrt{d}}\sum_{i}|ii\rangle$ is the maximally
entangled state. For each case of the dimension $d=2,3,4$, 3,000 states are
sampled by randomly selecting $\epsilon$ and $\rho_{0}$. Together with 1,000
random separable states, totally 4,000 states are sampled to serve as the
training quantum states.

We next generate training data based on the sampled quantum states. For the
correlation method, we prepare the training data by measuring the quantum
states via the local measurements given in Eqs.(7) and (8), where for each
party the number of measurement devices $N$ is fixed as 2. For the moment
method, the moments $\\{\mu_{m}(\rho_{A}),\mu_{m}(\rho_{B}),\mu_{m}(\rho)\\}$
are chosen as data features. Again, two sets of training data with $m=2$ and
$m=2,3$ are generated, in order to compare the power of moments of different
orders.

As mentioned, the labels of the training data are provided by the active
learning method. Meanwhile, the machine learning models for all the three
cases are the same, which is a 4-hidden-layer FNN with 400, 200, 100 and 50
neurons in each layer.

After training, to estimate the performance of these models, we first run them
on isotropic states, whose parameterized form is given by

$\displaystyle\rho_{d_{\epsilon}}^{me}=\frac{(1-\epsilon)}{d^{2}}I_{d^{2}}+\epsilon|\psi_{+}\rangle\langle\psi_{+}|,$
(10)

where $\epsilon\in[0,1]$. Given an isotropic state, its relative entropy of
entanglement can be analytically calculated Wei (2008). We make the
comparisons between the results given by all the three cases and the exact
results for $d=2,3,4$. The results are depicted in Fig.2, and the
corresponding MSEs are shown in the first part of Table 3.

Figure 2: The neural network predictions for $E_{R}(\rho_{d_{\epsilon}}^{me})$. The blue stars, grown dots and orange diamonds represent the prediction values of the correlation method, the moment method with $m=2$ and the moment method with $m=2,3$, respectively. The exact values are represented by the green triangles. All the corresponding MSEs are listed in the first part of Table 3. Table 3: The MSEs of predicting the relative entropy of entanglement for different methods and dimensions. Dimension | $d=2$ | $d=3$ | $d=4$
---|---|---|---
Isotropic |  Correlation, $N=2$ | $1.71\times 10^{-4}$ | $4.14\times 10^{-4}$ | $6.45\times 10^{-4}$
Moments, $m=2$ | $3.57\times 10^{-5}$ | $2.57\times 10^{-3}$ | $9.08\times 10^{-4}$
Moments, $m=2,3$ | $4.88\times 10^{-5}$ | $6.92\times 10^{-5}$ | $1.37\times 10^{-4}$
General |  Correlation, $N=2$ | $2.99\times 10^{-3}$ | $6.90\times 10^{-3}$ | $7.59\times 10^{-3}$
Moments, $m=2$ | $2.26\times 10^{-4}$ | $9.26\times 10^{-4}$ | $1.45\times 10^{-3}$
Moments, $m=2,3$ | $2.01\times 10^{-4}$ | $5.75\times 10^{-4}$ | $9.12\times 10^{-4}$

As illustrated, all the three cases (one for the correlation method and two
for the moment method) predict the relative entropy of entanglement of
isotropic states correctly with high precision. Particularly, the results
given by the active learning method, whose MESs are in the order of
$10^{-7}\sim 10^{-5}$, match the analytical results given by Ref.Wei (2008)
accurately.

Since the active learning method provides us a reliable way to calculate
relative entropy of entanglement, we now apply it to provide labels for random
quantum states, which allows us to test our models on more general quantum
states, not just isotropic ones.

Specifically, for each dimension, we sample 300 quantum states admit the form
Eq.(9), and combined with 100 random separable states, these 400 quantum
states serve as test quantum states. The generations of data features for
different methods are the same as before. Then we apply our trained FNNs to
predict the relative entropy of entanglement of these quantum states, the
corresponding MSEs are listed in the second part of Table 3.

It can be seen that in this task the overall behaviors of the MSEs of the
correlation method and the moment method are similar to those in the coherent
information prediction tasks. Actually in this task the performance of the
correlation method is even better, but a stable advantage of the moment method
can still be observed.

##  IV Hybrid Quantum-Classical Framework assisted Correlation Method

In Subsection.III, we fixed local measurements as the ones that achieve the
maximal violation of the CGLMP inequality, and fixed the number of measurement
devices for each party to be $N=2$. Due to these two constraints, the power of
the correlation method may be underestimated. Hence, a natural question is:
Can we improve the performance of the correlation method by relaxing the
constraints, say, enlarging the number of measurement devices $N$ or changing
the measurement devices?

###  IV.1 More Measurement Devices

Intuitively, enlarging the set of available measurement devices can probably
provide more information about target quantum states, and therefore may
improve the performance of the correlation method. To check whether this idea
works, we demonstrate it by setting $N=3,4$ in the coherent information
prediction tasks.

For each quantum dimension $d=3,5,8,10$, the sampled training and test quantum
states, and the mathematical structures of machine learning models keep the
same as before. The only difference is that now the training and test data
sets are composed of measurement outcome statistics involving $N=3,4$
measurement devices, rather than 2. It turns out that very limited
improvements are achieved by this change in predicting coherent information.
Table 4 lists the corresponding MSEs.

Table 4: The MSEs of predicting the coherent information with fixed measurement devices. Dimension | $d=3$ | $d=5$ | $d=8$ | $d=10$
---|---|---|---|---
CGLMP, $N=3$ | 0.0584 | 0.0356 | 0.0180 | 0.0170
CGLMP, $N=4$ | 0.0540 | 0.0364 | 0.0183 | 0.0138

The limited improvements given by increasing the number of measurement devices
are unexpected, because more measurement devices should have revealed more
information. A possible reason is that though the number of the measurement
devices we have utilized were increased, they are still of the form in Eqs.(7)
and (8). Therefore, to improve the performance of the correlation method
further, we need to find out whether this form is optimal or not.

###  IV.2 Learnable Measurement Devices

Figure 3: The hybrid quantum-classical framework for the correlation method.
(1) The input quantum state $\rho$, for which the full information is
required. (2) The generation of correlation, where the local measurements
$M_{x}^{a}$ and $N_{y}^{b}$ contain trainable variables. (3) The classical
training part, which can be, for example, a CNN.

Looking back at all our previous machine learning models that have been
discussed, we will see that all of them are classical models that deal with
pure classical information, where all data features of involved quantum states
are about correlation data or moments, which are essentially classical.
However, we now need to choose better local measurements to generate more
informative correlation data, which means that we have new quantum structures
to learn. For this purpose, below we introduce a hybrid quantum-classical
framework for our machine learning tasks.

In fact, a series of hybrid quantum-classical algorithms have been proposed
Peruzzo  _et al._ (2014); Mitarai  _et al._ (2018); McClean  _et al._ (2016);
Mari  _et al._ (2020); Liu  _et al._ (2021), where the concept of
parameterized quantum circuits (PQCs) is widely used to optimize a target
function by iteratively tuning parameters contained in underlying quantum
circuits. While in our hybrid quantum-classical framework, the tuning target
is quantum measurements. Since a general quantum measurement can be realized
by first performing a unitary operation and then measuring in the
computational basis, what we are aiming at is essentially learning $2N$ such
unitary operators.

Our hybrid quantum-classical framework works as follows. First, since we want
to find the best local measurements for our tasks, the entries of observables
$\\{M_{x}\\}$ and $\\{N_{y}\\}$ are now regarded as trainable variables and
will be updated repeatedly, which is the quantum part of our hybrid quantum-
classical machine learning framework. Second, after choosing the training
quantum states, we need to generate their training data features by measuring
the parameterized observables $\\{M_{x}\\}$ and $\\{N_{y}\\}$, and then
collecting the outcome statistics. Along with the correct labels, these data
features will be fed into the classical learning part, which is the same as
the previous models we have discussed. Fig.3 illustrates the whole framework.
It is worth mentioning that according to our previous experience, different
permutations of the measurement devices can result in quite different
prediction performance, especially when the training model is a CNN.
Therefore, introducing trainable measurements will bring us the optimized
permutation automatically.

We apply the hybrid quantum-classical machine learning framework to the
coherent information prediction tasks for the cases of $N=2,3,4,5$. The
sampled training and test quantum states, and the mathematical structure of
the classical machine learning part remain the same as Section.III.1, where
the dimension $d=3,5,8,10$. The experimental results are listed in Fig.4 and
Table 5.

Figure 4: Comparisons between the correlation method using the CGLMP measurements, learnable measurements and the moment method. The MSEs of the correlation method using learnable measurements decrease as $N$ increases. Particularly, when $N=3$ an apparent improvement is achieved. When $N=4$, the hybrid model even has the same order of MSE as the moment method with $m=2,3$ in high dimensions. The values of each MSE are listed in Table 5. Table 5: The MSEs of predicting the coherent information with learnable measurement devices. Dimension | $d=3$ | $d=5$ | $d=8$ | $d=10$
---|---|---|---|---
Learnable, $N=2$ | 0.0398 | 0.0345 | 0.0205 | 0.0143
Learnable, $N=3$ | 0.0149 | 0.0206 | 0.0121 | 0.0091
Learnable, $N=4$ | 0.0088 | 0.0184 | 0.0096 | 0.0081
Learnable, $N=5$ | 0.0085 | 0.0126 | 0.0080 | 0.0069

Comparing the case $N=2$ with the old results in Table 2, it can be seen that
almost no improvement is achieved for each dimension except $d=3$, which
indicates that in generating statistics data the local measurements in Eqs.(7)
and (8) are almost optimal for these cases. However, once we enlarge the
number of possible measurement devices by only one, that is $N=3$, the
situation becomes totally different, where the corresponding MSEs are 0.0206,
0.0121 and 0.0091 for dimensions $5,8$ and $10$, respectively. Compared with
the original version of the correlation method, $43.09\%,42.92\%$ and
$37.67\%$ improvements are achieved for dimensions $5,8$ and $10$,
respectively.

Furthermore, if we increase the number of measurement devices further, the
improvements are even more remarkable. As illustrated in Fig.4, when $d=3$,
our hybrid model with $N=4$ has the same order of MSE as the moment method
with $m=2$, and when $d=8,10$, it is even comparable with that of the moment
method with $m=2,3$. In fact, it turns out that even if we restrict Alice and
Bob to perform the same set of local measurements, the obtained MSEs differ
very little in high dimensions, e.g., $d=5,8,10$. Therefore, a half number of
the variables in the quantum part can be reduced to achieve a similar
performance, hence facilitating the training efficiency.

Recall that though the moment method has very good performance, the
experimental cost of measuring moments is much higher than that of generating
correlation data. Therefore, our result clearly implies that our hybrid
quantum-classical framework of machine learning largely overcomes this
difficulty, and achieves similar performance. This also shows that it is
already practical to quantify unknown entanglement on current NISQ devices.

##  V Conclusion

Quantifying unknown entanglement experimentally is a very challenging task
today. However, considering the profound importance of quantum entanglement
and the fast developments of quantum industries, fulfilling this kind of tasks
will probably become a daily routine in the future. As a result, finding
realistic and economical solutions for this problem from the viewpoints of
engineering is extremely necessary, and a typical one of this kind is using
machine learning methods.

In this paper, we focus on two known machine learning approaches for
quantifying unknown entanglement experimentally, one using moments as data
features, and the other using correlation data. We systematically compare
their performance in predicting the coherent information and the relative
entropy of entanglement. According to our results, the moment method behaves
much better than the correlation method, which is consistent with the
intuition. In fact, as mentioned, each moment is a quantity related to
eigenvalues of the corresponding density matrix, and all the moments together
in principle can pin down the essential information of the density matrix.
Hence, the moment method naturally enjoys remarkable advantage in quantifying
entanglement. However, estimating moments is much harder than correlation
data.

This motivates us to improve the performance of the correlation method in
quantifying entanglement. For this, we propose a hybrid quantum-classical
framework of machine learning models to generate more informative correlation
data. Specifically, we studied two possible directions to achieve this, where
one directly enlarges the number of local measurements fixed beforehand to
produce more informative correlation data, and the other adds a new quantum
module for the machine learning model to search for better local measurements
utilized in correlation data generations, and at the same time increases their
number. It turns out that the former direction only has little improvement,
while the latter is much better, and it even achieves comparable performance
with the moment method. Therefore, it is fair to say the hybrid quantum-
classical machine learning framework we propose has very good performance in
quantifying unknown entanglement experimentally. And we hope that in the
future such protocols can be deployed on NISQ devices.

###### Acknowledgements.

This work was supported by the National Key R&D Program of China, Grants No.
2018YFA0306703, 2021YFE0113100, and the National Natural Science Foundation of
China, Grant No. 61832015.

## References

  * Bennett and Wiesner (1992) C. H. Bennett and S. J. Wiesner, Physical review letters 69, 2881 (1992).
  * Bennett  _et al._ (1993) C. H. Bennett, G. Brassard, C. Crépeau, R. Jozsa, A. Peres, and W. K. Wootters, Physical review letters 70, 1895 (1993).
  * Ekert (1991) A. K. Ekert, Phys. Rev. Lett. 67, 661 (1991).
  * Horodecki  _et al._ (2009) R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki, Reviews of modern physics 81, 865 (2009).
  * Amico  _et al._ (2008) L. Amico, R. Fazio, A. Osterloh, and V. Vedral, Reviews of modern physics 80, 517 (2008).
  * Chuang and Nielsen (1997) I. L. Chuang and M. A. Nielsen, Journal of Modern Optics 44, 2455 (1997).
  * Poyatos  _et al._ (1997) J. Poyatos, J. I. Cirac, and P. Zoller, Physical Review Letters 78, 390 (1997).
  * Van Enk and Beenakker (2012) S. Van Enk and C. Beenakker, Physical review letters 108, 110503 (2012).
  * Brydges  _et al._ (2019) T. Brydges, A. Elben, P. Jurcevic, B. Vermersch, C. Maier, B. P. Lanyon, P. Zoller, R. Blatt, and C. F. Roos, Science 364, 260 (2019).
  * Moroder  _et al._ (2013) T. Moroder, J.-D. Bancal, Y.-C. Liang, M. Hofmann, and O. Gühne, Physical review letters 111, 030501 (2013).
  * Wei and Lin (2021) Z. Wei and L. Lin, Physical Review A 103, 032215 (2021).
  * Lin and Wei (2021) L. Lin and Z. Wei, Physical Review A 104, 062433 (2021).
  * Gray  _et al._ (2018) J. Gray, L. Banchi, A. Bayat, and S. Bose, Physical review letters 121, 150503 (2018).
  * Lin  _et al._ (2021) X. Lin, Z. Chen, and Z. Wei, arXiv preprint arXiv:2104.12527 (2021).
  * Roik  _et al._ (2022) J. Roik, K. Bartkiewicz, A. Černoch, and K. Lemr, arXiv preprint arXiv:2203.01607 (2022).
  * Zhou  _et al._ (2020) Y. Zhou, P. Zeng, and Z. Liu, Physical Review Letters 125, 200502 (2020).
  * Elben  _et al._ (2020) A. Elben, R. Kueng, H.-Y. R. Huang, R. van Bijnen, C. Kokail, M. Dalmonte, P. Calabrese, B. Kraus, J. Preskill, P. Zoller, _et al._ , Physical Review Letters 125, 200501 (2020).
  * Vermersch  _et al._ (2018) B. Vermersch, A. Elben, M. Dalmonte, J. I. Cirac, and P. Zoller, Physical Review A 97, 023604 (2018).
  * Ekert  _et al._ (2002) A. K. Ekert, C. M. Alves, D. K. Oi, M. Horodecki, P. Horodecki, and L. C. Kwek, Physical review letters 88, 217901 (2002).
  * Horodecki (2003) P. Horodecki, Physical review letters 90, 167901 (2003).
  * Mohri  _et al._ (2018) M. Mohri, A. Rostamizadeh, and A. Talwalkar,  _Foundations of machine learning_ (MIT press, 2018).
  * Shalev-Shwartz and Ben-David (2014) S. Shalev-Shwartz and S. Ben-David,  _Understanding machine learning: From theory to algorithms_ (Cambridge university press, 2014).
  * Wolpert and Macready (1997) D. Wolpert and W. Macready, IEEE Transactions on Evolutionary Computation 1, 67 (1997).
  * Schumacher and Nielsen (1996) B. Schumacher and M. A. Nielsen, Physical Review A 54, 2629 (1996).
  * Lloyd (1997) S. Lloyd, Physical Review A 55, 1613 (1997).
  * Cornelio  _et al._ (2011) M. F. Cornelio, M. C. de Oliveira, and F. F. Fanchini, Physical review letters 107, 020502 (2011).
  * Vedral (2002) V. Vedral, Reviews of Modern Physics 74, 197 (2002).
  * Hou  _et al._ (2020) S.-Y. Hou, C. Cao, D. Zhou, and B. Zeng, Quantum Science and Technology 5, 045019 (2020).
  * Vedral and Plenio (1998) V. Vedral and M. B. Plenio, Physical Review A 57, 1619 (1998).
  * Collins  _et al._ (2002) D. Collins, N. Gisin, N. Linden, S. Massar, and S. Popescu, Physical review letters 88, 040404 (2002).
  * Zohren and Gill (2008) S. Zohren and R. D. Gill, Physical review letters 100, 120406 (2008).
  * Barrett  _et al._ (2006) J. Barrett, A. Kent, and S. Pironio, Physical Review Letters 97, 170409 (2006).
  * LeCun  _et al._ (1998) Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner, Proceedings of the IEEE 86, 2278 (1998).
  * Albawi  _et al._ (2017) S. Albawi, T. A. Mohammed, and S. Al-Zawi, in  _2017 International Conference on Engineering and Technology (ICET)_ (Ieee, 2017) pp. 1–6.
  * Lu  _et al._ (2018) S. Lu, S. Huang, K. Li, J. Li, J. Chen, D. Lu, Z. Ji, Y. Shen, D. Zhou, and B. Zeng, Physical Review A 98, 012315 (2018).
  * Doherty  _et al._ (2002) A. C. Doherty, P. A. Parrilo, and F. M. Spedalieri, Physical Review Letters 88, 187904 (2002).
  * Wei (2008) T.-C. Wei, Physical Review A 78, 012327 (2008).
  * Peruzzo  _et al._ (2014) A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’brien, Nature communications 5, 1 (2014).
  * Mitarai  _et al._ (2018) K. Mitarai, M. Negoro, M. Kitagawa, and K. Fujii, Physical Review A 98, 032309 (2018).
  * McClean  _et al._ (2016) J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, New Journal of Physics 18, 023023 (2016).
  * Mari  _et al._ (2020) A. Mari, T. R. Bromley, J. Izaac, M. Schuld, and N. Killoran, Quantum 4, 340 (2020).
  * Liu  _et al._ (2021) J. Liu, K. H. Lim, K. L. Wood, W. Huang, C. Guo, and H.-L. Huang, Science China Physics, Mechanics & Astronomy 64, 1 (2021).