1. Introduction: The Simplest Quantum Many-Body Model Showing ♯P-Hard Complexity
The analysis of various quantum many-body systems capable of simulating ♯P-hard computational problems is one of the main topics of modern research in quantum physics and computing (see, for example, [
1,
2,
3,
4,
5,
6,
7,
8] and references therein). Recently, an atomic boson sampling of excited atom occupations in an equilibrium gas with a Bose–Einstein condensate (BEC) has been suggested as a process that could be ♯P-hard for classical computing [
9]. An example of a multi-qubit BEC trap [
10] shows that it could serve as a rich and, at the same time, convenient platform for studies of various phenomena associated with atomic boson sampling and quantum supremacy.
The present paper aims at the simplest possible model of the BEC trap that would allow one to greatly simplify the general theory outlined in [
9,
10] and explicitly disclose the mechanism behind the ♯P-hard computational complexity of quantum many-body systems. As a result, many general formulas and the entire theory, based on the method of characteristic function and Hafnian master theorem, acquire an explicit, transparent form clearly revealing the origin of the ♯P-hardness of computing the joint probability distribution of the excited-atom occupations. This ♯P-hardness is the ultimate reason for a potential quantum supremacy of the atomic boson sampling over classical-computing-based implementations of the generation of such random numbers. Remarkably, as is shown in
Section 8, the ♯P-hard complexity of quantum systems is equivalent to an intuitively obvious complexity of computing the multivariate integral for Fourier series coefficients of a sign-indefinite strongly-oscillating function.
Due to the interatomic interaction, nonzero mass of atoms, equilibrium, and an absence of external sources of bosons, the proposed atomic boson sampling has substantially different physics compared to photonic boson sampling in a linear interferometer, widely studied in the last decade, both theoretically [
2,
3,
4,
11,
12,
13,
14,
15,
16,
17,
18,
19,
20,
21,
22,
23,
24,
25,
26,
27,
28] and experimentally [
5,
29,
30,
31,
32,
33,
34,
35]. The interatomic interaction is especially important. It greatly complicates and changes the behavior of the quantum system of many bosons leading, in particular, to the fundamental phenomenon of two-mode squeezing of excited atom states predicted in [
36] and strongly pronounced in the statistics of the total noncondensate occupation [
37]. Nevertheless, the simplicity of the BEC-in-the-box model allows us to disclose analytically the quantum-statistical physics of sampling the interacting Bose atoms. Our results show that the joint probability distribution of atom numbers in atomic boson sampling is expressed via the matrix hafnian in a way similar to that of the Gaussian boson sampling which utilizes external sources of squeezed photons. Thus, the same ♯P-hardness of computing the hafnian is involved. It allows one to transfer a large part of the previously developed analysis of the computational complexity and quantum supremacy from photonic to atomic boson sampling.
In addition to fundamentals of quantum simulations, physics of occupation fluctuations in the excited atom states of the BEC gas is very important for various applications and other BEC-gas setups, for instance, for studying trap cells [
38,
39], BEC collapse [
40], squeezed states [
41], matter-wave interferometers [
42], including Ramsey [
41,
43] and Mach-Zehnder on-chip [
44] ones, etc.
The content of the paper is as follows. In
Section 2 we introduce the main quantities and notations as well as summarize some known facts relevant to the many-body BEC system of atoms in a box trap. In
Section 3 and
Section 4 we present the general solution for sampling probabilities obtained by means of the Hafnian master theorem and reveal two ingredients of quantum supremacy, squeezing and interference, via the Bloch–Messiah reduction of the Bogoliubov transformation. In
Section 5,
Section 6 and
Section 7 we illustrate the increasing apparent complexity of sampling probability patterns by changing the observational basis of excited atom states from the eigen-squeeze modes to more and more involved unitary mixtures of them.
Section 8 contains a discussion of the ♯P-hardness of computing the joint atom-number probabilities in a general case of arbitrary sampling states. The focus is on the origin of ♯P-hard complexity of atomic boson sampling due to interference and squeezing of the sampled atom states via their interplay with the eigen-squeeze modes and eigen-energy quasiparticles. Concluding remarks related to a possible experimental demonstration of atomic boson sampling and various manifestations of quantum supremacy constitute
Section 9.
2. Quantum Statistical Physics of Joint Fluctuations of the Excited Atom Numbers
Consider a box trap of volume
with periodic boundary conditions. Excited atoms are described by a field operator
where
denotes an operator annihilating an atom in the bare-atom excited state
of the single-particle Hilbert space
. The sum
denotes a summation over a basis of excited atom states in a box trap,
, enumerated by an integer
.
A condensate wave function in the box trap corresponds to the zero integer
and is uniform,
. We abide the Bogoliubov-Popov approximation [
45,
46] and replace the condensate annihilation operator by a c-number,
, where
is a mean number of atoms in the condensate.
The Bogoliubov Hamiltonian in the basis
is given, up to an insignificant additive c-valued constant, by a quadratic form in the creation and annihilation operators:
It is written via a (
)-block matrix
H which is built of matrices
, and
as its blocks, applied to a 2-block column vector
, consisting of the column vector
on top of the column vector
, and multiplied from the left by a 2-block row vector
. The bold-faced operator
or
denotes a column vector of annihilation or creation operators in the
basis;
is the single-particle energy operator for a bare atom of mass
m in the box trap. The nabla symbol ∇ stands for the three-dimensional (3D) vector differential operator. The interatomic interaction is determined by a constant
via the s-wave scattering length
.
The general notations for the field operator and Hamiltonian in Equations (
1) and (
2) are convenient for an abstract general analysis presented below in
Section 3,
Section 4 and
Section 8. For any particular choice of the basis of excited atom states
the field operator and Hamiltonian acquire more specific and transparent form. In a usual basis of plane traveling waves
with the wave vectors
enumerated by the integer 3-vector
, we have
Hereinafter, we denote the creation and annihilation operators of bare atoms in the traveling-plane-wave states
by symbols
and
, respectively, as opposed to the above general-case operators
and
.
Equilibrium quantum many-body statistics of atoms in the BEC trap at a temperature
T is determined by the statistical operator
. The atomic boson sampling implies sampling in accord with the joint probability distribution
of the occupation numbers
of the excited atom states
in Equation (
1) or some subset or groups of them preselected for measurenment by the appropriate detectors projecting atoms onto those states. This probability distribution is given by Fourier series coefficients,
of the characteristic function of the atom-number operators:
The symbol stands for a trace of an operator or matrix.
The characteristic function (
7) had been found analytically [
9] via a determinant function which can be easily computed in polynomial time for any finite-size matrix:
Here, the symbol
stands for the matrix block with zero entries.
The characteristic function is determined by a covariance matrix
whose entries are given by the quantum-statistical average of the normally ordered tensor product of the 2-block column vector
and 2-block row vector
of the atom creation and annihilation operators, that is, all possible self- and inter-mode normal,
, and anomalous,
,
, correlators. Each variable
of the characteristic function
appears in the matrix of variables
Z in Equation (
8) twice—via the entry
in each of the two identical diagonal blocks
; the symbol
denotes the identity matrix. Note, that in the literature on Gaussian states [
47], the covariance matrix is often defined with a half anti-commutator,
, replacing the normal product,
, of the creation/annihilation operators, that adds a half identity matrix
to our covariance matrix
G.
The characteristic function in Equation (
8) has the same form for an arbitrarily restricted, marginal subset of the excited states or coarse-grained groups of them, if they are considered irrespective of the other excited states. In order to average over all irrelevant excited-state occupations one just needs to nullify all irrelevant variables
and keep only those rows and columns in the covariance matrix
G that correspond to the chosen marginal subset of excited states. In order to combine some excited states into a coarse-grained group, it is necessary to set equal all of the variables
within such a group.
Consider the quasiparticle creation and annihilation operators, which constitute the column vectors
and
and diagonalize the Hamiltonian in Equation (
2),
The quasiparticle eigen energy is denoted as
and is given in Equation (
13) or (
48). The field operator of excited atoms in Equation (
1) also can be represented via the quasiparticle operators as a sum of both annihilation and creation quasiparticle operators:
The functions
and
constitute the two-component wave function of the
l-th quasiparticle. Canonical Bose commutation relations for the quasiparticle operators,
, where
is the Kronecker delta, imply the following normalization of the two components,
and
, of each quasiparticle wave function
Suppose one chooses plane traveling waves
as the bare-atom excited-state basis. Then the Hamiltonian in Equation (
5) acquires a diagonal form with the canonical Bogoliubov spectrum of eigen energies,
for the quasiparticle creation, annihilation operators
,
related to the corresponding bare-atom operators
,
as follows
Here, we introduced the amplitudes
and
,
of the functions constituting the two-component wave function of those traveling-plane-wave quasiparticles
In terms of such quasiparticles, the excited-atom field operator in Equation (
1) acquires the following explicit form
Generally, the quasiparticles are completely independent of each other. Thus, the correlations between the quasiparticle creation/annihilation operators
, analogous to the bare-atom correlations in Equation (
9), are given by the
-block diagonal matrix of the thermal occupations of quasiparticles
Employing the canonical Bose commutation relations
for bare atoms, it is easy to see that the covariance matrix (
9) of bare-atom creation/annihilation operators can be obtained from the matrix
D in a compact form,
via the Bogoliubov transformation which relates the uncorrelated quasiparticles to the squeezed and interfering bare-atom excitations. It is described by the
-block symplectic matrix [
9]
and its inverse matrix
as follows
Below we assume that some finite number
M of orthogonal excited states (i.e., atom wave functions)
are preselected for sampling, that is, for a multi-detector measurement of the atom numbers, and constitute a basis of a finite-dimensional subspace of the single-particle Hilbert space
. Their normal and anomalous correlations are given by the corresponding
-submatrix of the covariance matrix in Equation (
9). Suppose the preselected wave functions are coupled via the Bogoliubov Hamiltonian only between themselves. That is, their off-diagonal couplings
, and
in Equation (
2) with any wave functions outside the preselected subspace are zero or negligible. For example, in the uniform box trap the s-wave scattering in the Bogoliubov Hamiltonian couples just plane waves with opposite wave vectors
and
. Then the analysis gets easier. All bold-faced vectors (such as
) are reduced to
M-dimensional vectors. All
-block matrices (including the Hamiltonian,
H, Bogoliubov,
R, quasiparticle occupation,
D, covariance,
G, and variable,
Z, matrices) are reduced to the
-matrices containing corresponding
-blocks such as
, etc. For the sake of simplicity, we will denote any such matrix of a finite, reduced dimension
or
by the same symbol that stands for its infinite-dimensional counterpart.
4. Two Ingredients of Quantum Supremacy and ♯P-Hardness: Squeezing and Interference
The result in Equation (
25) explicitly shows that the ♯P-hard complexity and, hence, potential quantum supremacy of atomic boson sampling from the ♯P-hardness of computing the hafnian of the extended covariance-related matrix
that depends only on the nontrivial Bogoliubov-transformation matrix
R and mean quasiparticle occupations constituting the diagonal matrix
D in Equation (
18). Thus, the mystery of quantum supremacy is encoded in the structure of the Bogoliubov-transformation matrix in Equation (
20) and can be revealed via its unique, irreducible Bloch–Messiah representation [
55,
56,
57,
58,
59]
It gives the blocks of the Bogoliubov-transformation matrix in Equation (
20) in the form of a singular value decomposition,
and involves two unitary matrices,
W and
V, as well as the diagonal matrix of
M single-mode squeezing parameters
The Bloch–Messiah reduction (
26) can be written in the form
that states that the effective evolution (resembling the Hamiltonian one,
) of the quasiparticle and bare-atom creation/annihilation operators under the action of the unitary multimode squeeze,
, and rotation,
, operators, respectively, yields the same creation/annihilation operators:
Such an interpretation, typical for quantum optics [
55,
56,
57,
58,
59,
60], employs the multimode squeeze and rotation operators,
determined by a symmetric matrix
(built of the Hermitian matrices
r and
; the unitary matrix
is symmetric) and a Hermitian matrix
generating the unitary
, respectively. The Hermitian factor,
r, of the multimode squeeze matrix,
S, is a positive semi-definite Hermitian
matrix whose diagonal representation
is determined by the unitary
W as follows
so that
,
.
The unitaries
W and
V are chosen to satisfy the so-called rotation condition emphasized in [
56] in view of a possible nonuniqueness of the singular value decomposition, particularly in the presence of degenerate singular values. The singular vectors
and values
are the eigenvectors (comprising the unitary
W as columns) and the eigenvalues of the Hermitian factor
r of the squeeze matrix
S, respectively. The squeezing parameters, or eigenvalues,
constitute the unique, irreducible resource of the many-body interacting system and do not depend on the choice of bases or unitaries.
Equation (
26) describes the overall Bogoliubov transformation as a sequence of three
-block transformations
from the creation/annihilation operators in the observational basis of the excited atom states preselected for sampling measurement to the creation/annihilation operators of completely independent quasiparticles which correspond to the diagonalized Hamiltonian in Equation (
10) and stay in totally separable, disentangled equilibrium states. Their combined state is described by the density matrix, or statistical operator,
expressed in terms of the eigen-energy quasiparticle creation and annihilation operators.
The three
blocks of the Bogoliubov reduction in Equation (
26) could be thought of as the matrices transforming the annihilation and creation operators or the corresponding wave functions in Equation (
11). This is explained below by means of the schematic diagram in
Figure 1 and Equations (
35)–(
38).
The first part of the Bogoliubov transformation,
, corresponds to the unitary rotation
,
of the basis of the bare-atom excited states
, prescribed for measurement of atom numbers by means of multi-detector imaging in the process of sampling, into the basis of the eigen-squeeze modes
associated with the eigenvectors
of the Hermitian factor of the squeeze matrix
r in Equation (
33) and explicitly expressed below in Equation (
38) via the quasiparticle wave functions. The corresponding transformation from annihilation operators of the observable states,
, to the annihilation operators of the eigen-squeeze modes,
, is performed by the unitary
V and does not involve creation operators.
The second part of the Bogoliubov transformation,
, is associated with the presence in the Hamiltonian the terms beyond the resonant-wave approximation (non-RWA terms),
and
, which create or annihilate a pair of excited atoms. It converts the state of the system into the squeezed state in which each single-particle eigen-squeeze mode
constitutes the same spatial profile for both wave-function components
of the
l-th eigen-squeeze two-component excitation which acquires the corresponding nontrivial squeezing parameter
. Note, according to Equation (
38), those two components have different amplitudes,
and
.
The third part of the Bogoliubov transformation,
, relates the creation,
, and annihilation,
, operators of the eigen-squeeze two-component excitations, formed on the second step, to the eigen-energy-quasiparticle creation,
, and annihilation,
, operators:
According to the form of the quasiparticle field operator in Equation (
11), it means two simultaneous unitary rotations (under the action of one and the same unitary matrix
) of the single-squeeze modes
into the wave functions
and
, which according to Equation (
11) constitute the bases of the first and second components of the two-component functional space of the quasiparticle wave functions
,
Since the quasiparticle wave functions
and
are fully predetermined, fixed by the coupling parameters (interactions) in the Hamiltonian (
2) (of course, up to a possible degeneracy), the unitary matrix
determines the unique (up to a possible degeneracy) eigen-squeeze modes
which are the same for both components of the quasiparticle functional space as per equations inverse to Equations (
37):
Each single-squeeze two-component excitation originating from the mode
owns the single-mode squeezing parameter
and is not subject to inter-mode squeezing with other eigen-squeeze two-component excitations.
The existence of such a unique unitary
W, simultaneously converting the basis wave functions
and
of both components of the two-component functional space of quasiparticles into the basis wave functions
and
which are equal to the same eigen-squeeze mode
just multiplied by different constant factors
and
, respectively, is a nontrivial and important property of the Bogoliubov transformation. It is a consequence of the symplectic property [
9],
of the Bogoliubov transformation, that is, canonical Bose commutation relations, and highlights the fact that the two components of the quasiparticle eigen function
are not independent, but, on the contrary, are deeply inter-correlated.
Both the first,
, and the third,
, factors of the Bogoliubov transformation (
26) introduce unitary interference, that is entanglement, into the quantum many-body state (statistical operator
) of
M excited atom modes chosen for atomic boson sampling. They are controlled by two different means. The trapping potential and parameters of the Hamiltonian in Equation (
2), that is, couplings
and atom interactions, control both unitaries
W and
V since they determine the composition of the eigen-squeeze modes both with respect to the quasiparticle wave functions and the bare-atom wave functions. The unitary
V of the first factor is additionally controlled via a choice of different excited atom states by means of reconfiguring the multiple detectors for atomic sampling. In the particular case of a box trap with a uniform condensate, discussed in the present paper, the third Bogoliubov transformation
is almost completely fixed. Yet, the first Bogoliubov transformation
provides access to practically arbitrary unitary matrices
V, controllable in a wide range of their parameters. This should be enough for ensuring ♯P-hardness of computing the hafnian in Equation (
25) on average. Such complexity on average [
3,
12] is important for the possibility of demonstrating the quantum supremacy of atomic boson sampling.
Thus, the Bloch–Messiah reduction in Equation (
26) unambiguously specifies two preferred bases: (i) the basis of the quasiparticle operators in Equation (
11) ensuring the diagonal form of the Hamiltonian (
10) and (ii) the basis of the eigen-squeeze single-particle excited states in Equation (
38) diagonalizing the Hermitian factor of the multimode squeeze matrix, Equation (
33). Moreover, the Bogoliubov transformation in the Bloch–Messiah representation explicitly relates both above-mentioned bases to the observational basis of excited bare-atom states
which can be arbitrarily selected by a reconfiguration of atom detectors.
The interference between wave functions of those bases, appearing in the statistical operator (state) of the atomic many-body (multimode) system due to the unitaries and leading to the entanglement of different Bose bare-atom modes as opposed to the separability of quasiparticle states, constitutes the first of the two ingredients of the computational ♯P-hardness and potential quantum supremacy of the atomic boson sampling in equilibrium.
The second ingredient is the squeezing of the equilibrium state of the excited-atom modes originating from the reduced, canonical Bogoliubov transformation
determined exclusively by the eigenvalues
of the Hermitian factor of the multimode squeeze matrix
r. If it was given by the identity matrix with all single-squeezing parameters equal zero,
, then the hafnian in Equation (
25) would be reduced to the permanent of a positive matrix, which according to [
2,
14,
26] could be approximated by the Stockmeyer’s approximating algorithm [
61] in a computational-complexity class simpler than
. A physical mechanism of squeezing can be seen from Equation (
11) for the atomic field operator. Any quantum or thermal fluctuation associated with disappearance of an atom at a point
implies annihilation of a superposition of quasiparticles with the amplitudes given by the first component of the quasiparticle wave functions
and simultaneous creation of a superposition of quasiparticles with the amplitudes given by the second component of the quasiparticle wave functions
.
In the observational basis coinciding with the eigen-squeeze modes
, when
, the covariance matrix (
19) acquires a unique irreducible form
of a sum of the pure quantum (independent on temperature and associated with a quantum depletion of the condensate due to interatomic interaction) and complimentary thermal correlations between creation/annihilation operators
of the eigen-squeeze modes:
Here,
is the
matrix defined by thermal occupations of quasiparticles. These irreducible contributions are determined exclusively by the two intrinsic entities of the BEC gas, the eigen-squeeze modes and the eigen-energy quasiparticles, and are not subjected to the arbitrariness of choosing any observational basis. In any observational basis
matrices
G and
C, Equation (
22), appears as the corresponding unitary transform of the irreducible one as per unitary rotation in Equation (
35),
Below we illustrate and discuss the interplay and manifestations of two aforestated ingredients of quantum supremacy in atomic boson sampling for different choices of the observational basis of excited bare-atom states .
5. Separability of Atomic Boson Sampling in the Basis of Single-Mode-Squeezed Standing Plane Waves: Atom-Number Statistics via the Hafnian and Legendre Polynomials
A minimal complexity of joint atom-number probabilities in the atomic boson sampling is achieved if one chooses the eigen-squeeze modes in Equation (
38) as the bare-atom excited states for atom-number measurements by the multi-detector imaging system. For the considered model of the uniform condensate in a box trap with the periodic boundary conditions, the wave functions of the eigen-energy quasiparticles could be chosen in such a way that they coincide with the eigen-squeeze modes. Let us choose them to be
and
spatial modes with a wave vector
. In fact, any orthogonal (but not any unitary) transformation of a pair of energy-wise degenerate eigen-squeeze excited states would lead to the same atom-number sampling statistics.
Thus, the excited atoms are described by the field operator
written via operators
annihilating a bare atom in the corresponding sinusoidal states with a wave vector
. Hereinafter, an apostrophe in the symbol of the sum
means summation over nonzero wave vectors modulo multiplication by
, that is, the integer vector
is running over a half of a three-dimensional (3D) integer space,
, with its origin excluded. Such an enumeration is convenient for the subsequent analysis of the particular cases of standing plane waves, corresponding to the opposite wave vectors
and
, because it allows one to explicitly take into account the fact of a two-fold degeneracy of their energies in the box trap. In other words, for each pair of terms corresponding to the opposite nonzero wave vectors
and
only one belongs to the sum. (Note, that the sum in Equations (
1) and (
4) does not include the apostrophe.)
In such a sinusoidal, standing-wave basis of excited states,
the Hamiltonian (
2) splits into two independent parts involving separate operators related to
or
modes,
It is in contrast to the Hamiltonian in Equation (
5) in which the modes
and
of the traveling-plane-wave basis
were coupled to each other. The Hamiltonian in Equation (
45), compared to the Hamiltonian in Equation (
5), has the same bare-atom energies
and normal overlapping integrals
, but different, now diagonal anomalous overlapping integrals
In this case, the quasiparticle annihilation operators
are given by a similar Bogoliubov transformation via the corresponding bare-particle annihilation and creation operators:
The quasiparticle operators diagonalize the Hamiltonian (
45),
and turn the field operator (
43), annihilating bare atoms, into a sum of both annihilation and creation quasiparticle operators:
In both exponential and sinusoidal, bases, the quasiparticle energies
and factors
are the same, given in Equations (
13) and (
15). The energy
is larger than the bare-atom energy
.
Thus, in the observational basis of the eigen-squeeze modes (
44) the sampling joint probability distribution factorizes into the product of probabilities for single
or
modes. Each of such single-mode-squeezed probabilities can be calculated analytically by means of the general solution in Equation (
25). We just need to appreciate the fact that in this case the general Bloch–Messiah representation of the Bogoliubov transformation in Equation (
26) is reduced to the simple independent blocks with trivial unitary parts
, that is
This is the Bogoliubov single-mode squeezing in the simplest, pure form with a single-mode squeezing parameter
.
The corresponding single-mode covariance matrix is given by Equations (
40) and (
41) as follows
Here,
and
stand for the average number of quasiparticles, normal correlator and anomalous correlator of the single sinusoidal mode, respectively. The anomalous correlator is negative for the chosen phases of the basis wave functions and its absolute value is bounded from above by the inequality in Equation (
51) which is equivalent to the obvious inequality
. The maximum value of the anomalous correlator,
, is achieved when
, that is, when the quasiparticles are in the vacuum state. The relations between the average quasiparticle occupation
, the single-mode squeezing parameter
and the normal and anomalous correlators
and
are illustrated in
Figure 2. Any fixed value of the normal correlator which is equal to the mean number of atoms in the eigen-squeeze mode,
, corresponds to a circle on the plane
with the origin at the point
, i.e.,
.
The eigenvalues of this covariance matrix depend on the single-mode squeezing parameter
and the mean quasiparticle occupation
,
The atom-number probability distribution in a sinusoidal mode is determined, as per Equation (
22), by the covariance-related matrix
which has the following explicit form:
It is convenient to denote the entries
of the matrix
by adding apostrophe to the symbols
denoting the entries of the covariance matrix
. The eigenvalues
of the renormalized covariance matrix
look similar to the eigenvalues
of the covariance matrix
, namely,
Note, that the maximal value of the absolute value of the anomalous correlator,
, is achieved when
.
Now we calculate the atom-number probabilities. The result is immediately given by the Hafnian master theorem (
25),
where the extended covariance-related matrix
involves the
matrix
with all entries equal unity.
In the absence of anomalous correlations, when
, the eigen-squeeze modes are in a non-squeezed thermal state and the hafnian in Equation (
55) is reduced to the known permanent of the unity matrix
, that is,
. In this case, the atom-number probability distribution is reduced to the simple exponential law
typical for a counting statistics of a mode in the thermal state [
47,
57,
62].
The presence of a nonzero anomalous correlator
introduces squeezing and makes this distribution nontrivial. Fortunately, the hafnian in Equation (
55) in the case of an arbitrary
is easy to calculate via the known recursive relation for hafnians [
54]. Performing two recursive steps and excluding the hafnian of an auxiliary matrix, one finds the following second-order recursive formula
It starts from the plain values of the hafnian
at
and
at
(as per convention in the Hafnian master theorem). Comparing it with the well-known recursive relation for the Legendre polynomials [
63]
,
we at once reduce the hafnian to the Legendre polynomials:
Thus, the atom-number probabilities for sampling the occupation of the single-mode-squeezed standing plane wave are as follows
The second equality is due to substitution
. Note, that all Legendre polynomials
have the same, independent on the atom number
n, argument which is determined only by the normal and anomalous correlators. Instead, the atom number
n appears in the order of the Legendre polynomials.
Alternatively, the hafnian in Equation (
55) can be calculated in a straightforward way via elementary combinatorics as follows
Here, the symbol
stands for the largest integer less than or equal to
, and the essential argument is
According to the identity (3.137) in [
64], the sum in Equation (
60) is proportional to the Legendre polynomial of the order
n. Thus, the result stated above in Equation (
59) is rederived.
The recursive relation in Equation (
57) suggests the suppression of odd-occupation probabilities at small values of
. In the extreme case of zero mean occupation of quasiparticles, that is when
, or
, the recursive relation becomes purely two-step,
Equation (
62) yields the result
which is equivalent (due to identities
and
) to the one given by the reduction of the hafnian of a
-block-diagonal matrix to the square of the hafnian of just one of matrix’s blocks,
where the symbol
denotes the
matrix with zero entries.
It leads precisely to the formula for the well-known occupation probabilities for a mode in the squeezed vacuum state [
47,
62],
Thus, all odd-occupation probabilities are zero, and only even occupation numbers show up in the atomic boson sampling. This is not surprising since the zero average quasiparticle occupation,
, means that the mode is in the vacuum Fock state
with respect to the quasiparticle operators.
Typical occupation probability patterns given by Equation (
59) for sampling from eigen-squeeze bare-atom excited states (
44) are illustrated in
Figure 3. In the range
, where the absolute value of the anomalous correlator is less than the critical value
, an increasing absolute value of the anomalous correlator
leads to growing probabilities of occupations which are less than some number
and suppressing probabilities of occupations which are larger than that number
. Such a behavior drastically switches to a completely different behavior when the absolute value of the anomalous correlator exceeds the critical value,
, at which the argument of the Legendre polynomial jumps from the pure real to pure imaginary values through the infinity. Namely, then the probabilities of even occupation numbers start to grow, while the probabilities of odd occupation numbers start to drop down. At the maximal possible value of the anomalous correlator,
, which corresponds to the squeezed vacuum state,
, all probabilities for even occupation numbers become zero.
The statistics of the eigen-squeeze-mode occupation can be deduced directly from its characteristic function
. It gives explicit information on the moments as well as ordinary,
, and generating,
cumulants of the sampling probabilities via the Taylor series of its logarithm,
This characteristic function is given by the general result in Equation (
8) as follows
In Equation (
67) we provide two expressions for it. The first one is convenient for the cumulant analysis, while the second one—for the straightforward calculation of probabilities. In particular, we see that the second expression (for its equivalent form, see [
65]) is proportional to the well-known generating function of the Legendre polynomials [
63] that gives another proof of Equation (
59). The first expression yields the exact result for the generating cumulants,
where
, Equation (
52), are the eigenvalues of the covariance matrix (
51). The value of the first generating cumulant
yields the mean occupation
whose thermal and quantum contributions are equal to
and
, respectively. Adding the value of the second generating cumulant, we immediately obtain the standard deviation,
. The higher moments of the atom-number probability distribution can be found in a similar way.
It is worth noting that the presence of the nontrivial anomalous correlator makes this characteristic function in Equation (
67) very different from the one,
, associated with the atom-number sampling statistics in the limit of non-interacting, ideal Bogoliubov gas (IBG,
) when all excited atom modes remain in a non-squeezed state.
6. Two-Mode Squeezing of the Atomic Boson Sampling in the Basis of Traveling Plane Waves: Reduction of the Hafnian to the Permanent and Hypergeometric Function
Let us see now how the sampling statistics are getting more complex due to adding the effect of interference on top of the effect of pure squeezing considered in the previous section. The minimal complication occurs when we choose the observational basis for sampling
(which is the set of bare-atom excited states for atom-number measurements) to be the traveling plane waves
, Equation (
3). This is the simplest possible unitary mixing of the eigen-squeeze modes (
44) employed as the observational basis in the previous section. So, the excited atoms are described now by the field operator (
4) via the annihilation operators
which differ from the annihilation operators
of the eigen-squeeze modes due to the unitary transformation
As a result, the Bogoliubov transformation to the eigen-energy quasiparticles in its Bloch–Messiah reduction form (
26) acquires, in addition to the central pure squeezing block
in Equation (
50), the right-side unitary block
with the unitary matrix
V which is given in Equation (
69) above and is not equal to the identity matrix anymore. The left-side unitary block in the Bloch–Messiah reduction (
26) remains trivial,
, since it is the unitary transformation from the eigen-squeeze two-component excitations to the eigen-energy quasiparticles (see
Figure 1 in
Section 4) which in the case of the uniform BEC in a box trap coincide with the eigen-squeeze two-component excitations. The unitary transformation of the basis excited states corresponding to the operator transformation in Equation (
69) is performed by the unitary
as per Equation (
35),
Obviously, the introduced interference proceeds independently within each
-block of the atom excited states with the wave vectors
and
. Hence, it suffices to consider its effect on the sampling statistics just for one of such blocks—the block corresponding to the following 4-vector of creation and annihilation operators
. The related covariance matrix defined in Equation (
9) is equal to
It immediately follows from the general formula in Equation (
19) after plugging in the Bloch–Messiah reduction of the Bogoliubov transformation stated above. Equivalently, it can be easily obtained by means of the unitary transformation
V from the covariance matrix in Equations (
40)–(
42) or (
51) written for the 4-vector
of the creation and annihilation operators of the eigen-squeeze modes in the form of the
-matrix
The normal and anomalous correlators in the exponential-function basis remain the same as they were in the sinusoidal-function basis and are equal to
and
, respectively (see Equation (
51)).
It is easy to find the covariance-related matrix in Equation (
22) via Equation (
71) in the explicit form,
where the parameters
are the entries of the covariance-related matrix
C in Equation (
53) describing the eigen-squeeze mode
(see Equation (
44)).
Now we can calculate the joint probability distribution of atom numbers in the two excited states which are the traveling plane waves
and
for a given wave vector
. The result is provided by the Hafnian master theorem (
25),
Here, the hafnian matrix function is applied to the extended covariance-related matrix
which is a block matrix. The block
is a
matrix with all entries equal unity,
is a zero matrix of the same size. The determinant in the denominator has been calculated explicitly as follows:
.
The matrix
consists of blocks of unequal dimensions. However, its total dimension
is even. Such a hafnian
is easy to compute. The simplest way to do so is via consecutive swapping the 2-nd and 4-th block-rows and then the 2-nd and 4-th block-columns. This operation keeps the hafnian invariant and allows us to convert the matrix into a block-counter-diagonal form for which the hafnian is reduced to the permanent as per the well-known identity
valid for any square matrix
A. As a result, we calculate the required hafnian as follows
Here, the permanent has been calculated explicitly by simple combinatorial means via binomial coefficients
. Indeed, one just calculates the number of permutations of
elements that swap
k elements between the block of the first
elements and the block of the last
elements, which corresponds to the summand
. The obtained sum is an ordinary (Gaussian) hypergeometric function
with integer parameters which is a polynomial. It could be also expressed in terms of the Jacobi polynomial [
65]; however, the hypergeometrical representation is more convenient. The argument of the hypergeometric function in terms of the eigen-energies
and eigen-squeezing parameters
is given in Equation (
61).
The hafnian in Equation (
74) can be calculated also directly from its combinatorial definition [
54]. Each product of
entries contributing to the sum in the hafnian’s definition is encoded by a division of the matrix dimension
into
unordered pairs. Elements from the first group of
elements could be paired either to the elements from the second group of
elements, which corresponds to picking an
entry from the non-diagonal block, or to the elements from the third group of
elements, which corresponds to picking an
entry from the non-diagonal block. Pairing to the same group or the fourth group means picking a zero element, which vanishes the summand. Thus, each permutation encoding nonzero summand swaps
k elements between the first and second blocks (of
and
elements, respectively), as well as swaps
k elements between the third and fourth blocks. The rest of the elements should be swapped between the first and third blocks, or between the second and fourth blocks.
Let
be a number of pairings between the first and second groups. It is also the number of pairings between the third and fourth groups. There are
alternative ways to choose a pairing between the first and second groups, and each variant corresponds to the multiplier
accumulated from a non-diagonal block of
. The same holds for pairing between the third and fourth groups. Both the first and the third groups have
unpaired elements in rest. There are
alternative ways to pair them correspondingly, each corresponds to the multiplier
. Similarly, both the second and the fourth groups have
unpaired elements in rest. There are
alternative ways to pair them correspondingly, each corresponds to the multiplier
. Thus reaching the result stated above in Equation (
76).
Combining Equations (
74) and (
76), we obtain the explicit formula for the joint probability distribution in the case of sampling from the atom excited states given by two counter-propagating traveling plane waves via the ordinary hypergeometric function:
An equivalent result may be derived via a straightforward differentiation of the characteristic function, but the corresponding calculations are more cumbersome [
65].
In the particular case of an ideal, non-interacting BEC gas, when the anomalous correlator equals zero,
, and the excited atom states are not squeezed, the calculation of the joint probability distribution in Equation (
74) becomes elementary since the hafnian of the corresponding, extremely degenerate matrix
is reduced to the product of two permanents:
Thus, the joint probability distribution for occupations of the two degenerate counter-propagating waves separates into the product of two independent exponential single-mode distributions similar to the thermal counting statistics, Equation (
56).
In the presence of nonzero anomalous correlators, the joint probability distribution could not be factorized anymore. This is due to the appearance of the intra-modal squeezing. As illustrated in
Figure 4, increasing the absolute value of the anomalous correlator results in enhanced correlations between the sampled atom numbers
and
and growing the probability of the completely entangled occupation states
. This occurs in accord with a simultaneous increase of the pure quantum contribution
, Equation (
40), associated with the quantum depletion of the condensate and growing up single-mode squeezing parameter
, to the covariance matrix in Equation (
42). The point is that the thermal quasiparticle occupation
and, hence, the relative value of the complimentary thermal contribution
, Equation (
41), tend to zero when the anomalous correlator approaches it maximum value as per
Figure 2. As a result, the evolution of the probability distribution pattern in
Figure 4 looks like the formation and appearance of the quantum entangled ridge out of the disappearing thermal covering layer, or background. Obviously, the full complexity of the sampling probability patterns becomes clearly visible when the relevant excited bare-atom
-modes have a close-to-maximum anomalous correlator
and small mean quasiparticle occupation
. Achieving this in the BEC experiments requires a proper choice of the sampling
-modes as well as increasing the interatomic interaction, that is, the quantum depletion of the condensate, and decreasing the temperature, respectively.
The leading term in the asymptotics of the joint probability distribution in Equation (
77) when the anomalous correlator approaches its extremum
can be calculated directly from the hafnian formula in Equation (
74) by means of setting
and employing the identity (
75) as follows
It coincides with the contribution of the highest-order monomial in the hypergeometric-function polynomial in Equation (
77) and yields purely diagonal joint probability distribution in which the probability
exponentially decreases with increasing atom numbers
.
A well-known case of the squeezed vacuum state corresponds to the point at the very expremum,
, when the asymptotics (
79) is reduced to a well-known result for the two-mode squeezed vacuum state [
47,
62]
7. Nontrivial Effect of Interference on the Atomic Boson Sampling in the Basis of Any Two-Mode-Squeezed Unitary-Mixed Degenerate Standing Plane Waves
Let us look at further complications in the sampling probability patterns due to the nontrivial effect of interference in a more general set of the observational excited atom states. Consider atomic boson sampling from two excited atom states formed from the two eigen-squeeze modes
,
by means of an arbitrary unitary mixing via the unitary matrix
It corresponds to the first part of the Bogoliubov transformation in
Figure 1, from the annihilation operators
of these two excited atom states to the annihilation operators
of the two eigen-squeeze modes. Contrary to the particular case of two counter-propagating plane waves in Equation (
69), now the unitary
V involves two arbitrary real-valued angles
and
. The matrix in Equation (
81) has only two free parameters opposite to an arbitrary unitary matrix,
,
, involving four parameters. Here, the number of parameters is halved since the gauge phase factors of the two excited wave functions chosen to constitute the observational basis do not play any physical role. In fact, the whole set of observational bases which correspond to physically different joint statistical distributions of atom numbers is parameterized by matrix
V in the form of Equation (
81) with
,
. In view of Equation (
35), wave functions of the selected measurement basis,
are the superpositions of standing and traveling plane waves.
The covariance matrix for the operators
could be easily obtained from the covariance matrix (
72) formed by the operators of the eigen-squeeze sine and cosine modes via the following transformation generated by the unitary
V:
Here, the diagonal block containing normal commutators is proportional to the identity matrix, while the unitary-induced interference affects the block with anomalous correlators. The covariance-related matrix, which determines the matrix in the Hafnian master theorem, takes the following form
Here, the amplitudes
and
(see Equation (
53)) are the same as in
Section 5 and
Section 6.
Note, that if the
symmetric matrix
appeared in Equation (
84) is proportional to the identity matrix and, hence, the same is true for the complex conjugated matrix
, computing probabilities is reduced to the simplest case of two independent eigen-squeeze modes solved in
Section 5. This case corresponds to selecting two orthogonal standing waves as a measurement basis. In particular, such simplification happens if the matrix
V is chosen to be orthogonal, so that
and
.
The sampling probabilities are given by the Hafnian master theorem (
25),
In the present case of a general-type unitary mixing the hafnian in Equation (
85) strongly depends on the variable entries of the nontrivial symmetric matrices
and
in Equation (
84). Let us denote these entries as follows
We can calculate the hafnian directly from its combinatorial definition implementing its further reduction to a combinatorial problem of counting different partitions of
elements belonging to one of the four different groups (corresponding to different entries on the main diagonal) into pairs, either by linking an element from one group with an element from another group or by pairing elements within the same group. Finally, we obtain the hafnian as the following sum
It runs only over such subsets of indices that yield integer numbers under factorials in the denominator of Equation (
87); otherwise, the summand should not be included in the sum. In other words, all numbers
and
and
and
should be even.
Figure 5 exemplifies how correlations between occupation numbers
and
arise while one varies the unitary
V making each observational-basis bare-atom excited state more inter-correlated combination of the eigen-squeeze modes. The occupation statistics of the two eigen-squeeze modes are separable (see panel (a)). While one switches the observational basis from the eigen-squeeze standing waves to a basis consisting of partially-traveling waves, there appear nontrivial regions of the enlarged probabilities for some pairs of atom numbers
. Switching to purely traveling waves, which is the case considered in the previous section, ultimately leads to forming a ridge extending along the diagonal
direction (see panel (d)).
In general, the probability pattern has a nontrivial structure determined by how the unitary mixing distributes the anomalous correlator, brought in by squeezing, over the entries
, and
. An example is a two-crest structure whose divergence angle depends on the unitary angles
(see panels (b), (c)). It is strongly pronounced for the anomalous correlator values
close to its extremum,
, since then the pure quantum, condensate-depletion-based contribution
, Equation (
40), to the covariance matrix (
42) dominates the thermal one
, Equation (
41). Otherwise, for smaller
and larger
, the nontrivial correlation pattern is masked by thermal contributions due to large thermal population
of quasiparticles and small single-mode squeezing parameter
as per
Figure 2.
All of the probability distribution patterns on the plane of stochastic variables
and
shown in
Figure 5 are symmetric because we consider here the unitary mixture of two degenerate eigen-squeeze modes. They have the same entries
and
in their covariance-related matrices (see Equation (
53)), which implies that the “diagonal” anomalous correlators have the same absolute value for any orthogonal states of the observational basis, in particular,
.
For a squeezed vacuum state, which implies zero mean quasiparticle occupation, the off-diagonal super-block of the extended covariance-related matrix
, Equation (
85), is zero since
. Thus, the hafnian determining the probabilities is simplified to the product of the hafnians of two diagonal blocks which are complex conjugated to each other,
for even
, and zero otherwise. Here,
stands for the maximal integer less or equal
, where
is a minimal of two atom numbers
and
. Such a case of unitary mixed squeezed vacuum in optics had been described in [
66] via integrals of Hermite polynomials, and the probability distribution had been finally reduced to the associated Legendre polynomial which parameters are determined by
and
. Note, that in the squeezed vacuum state the probability of sampling the atom numbers
and
of different parity is always zero. This is an immediate consequence of the fact that the hafnian of the matrix of odd size
is identically equal to zero.
The squeezed vacuum state corresponds to the most contrasting correlation pattern of the joint probability distribution. We illustrate this thesis by
Figure 6 plotted for the case of the vacuum state with the normal correlator
and the unitary mixing matrix
. The wave functions of the selected measurement basis represent some mixtures of standing and traveling plane waves. On the plane
of sampling outcomes one can clearly see the dedicated directions along which the probabilities are relatively large. This global landscape is also essentially decorated by a check-mate pattern representing the strong effect of probability cancellation for odd-parity total occupation stated above. The probabilities of adjacent outcomes typically dramatically differ from each other even if they are all close to the directions of large probabilities. There are even gaps of unlikely outcomes surrounded on all sides by outcomes with higher probabilities (for example, around the point
).
The presence of thermal excitations, i.e., the presence of a nonzero number of quasiparticles in the system, makes the correlation pattern not that sharp, as is seen from comparison of
Figure 5 and
Figure 6, panel (c). Both of them refer to exactly the same matrix
V and observational basis, while the anomalous correlator
in
Figure 6 is larger. The dedicated directions of more probable outcomes in
Figure 5 are still the same as in
Figure 6 and can be clearly seen, but they turn into gently sloping hills, and the abrupt checkerboard pattern has disappeared.
Note, however, that thermal excitations do not simply raise the background and the pattern of vacuum-state statistics does not just draw into them. One should take into account not just relative increase of the thermal contribution
, Equation (
41), to the covariance matrix (
42) due to increasing quasiparticle population
, but also decrease and restructuring of the pure quantum contribution
, Equation (
40), due to simultaneous decrease of the single-mode squeezing parameter
as per
Figure 2, and mostly important a nontrivial mixing of the thermal and quantum contributions in the covariance matrix (
42) due to the unitary rotation
V. The redistribution of probability which makes the landscape more smooth happens first among adjacent (or close enough) outcomes. This is also seen from comparison of (d) panels in
Figure 4 and
Figure 5, representing the joint occupation statistics for two counter-propagating plane waves for the states with no and few quasiparticles, respectively. While in
Figure 4d, corresponding to the vacuum state, there is a well-pronounced and sharp diagonal ridge at
, the thermal excitations make it wider and enlarge, first of all, the probabilities of the adjacent states,
. In general such behavior could be described by the formula for hafnians in Equation (
87) which essentially allows one to employ a perturbation approach with respect to small value of the entries involving the parameter
(which are exactly zero for the vacuum state).
Importantly, while controlling the observational basis via the unitary
V may lead to nontrivial joint statistics of the atom numbers with widely ranging correlation patterns as described above, it does not affect the total noncondensate occupation statistics at all. This fact is not obvious from the formula for the joint probabilities since the direct sum representing the total-occupation probability,
may look nontrivial in virtue of Equation (
85) involving complicated hafnians.
The distribution of the total number of atoms in any selected subset of excited states is easier to described via the characteristic function in Equation (
8) with all arguments within this subset set to be equal,
. (Setting some groups of arguments equal to each other amounts to calculating a coarse-grained statistics.) Then, the corresponding submatrix of variables
Z in Equation (
8) turns into the scaled identity matrix
, which commutes with any matrix
V. Thus, the unitary matrix, switching the observational basis and transforming the covariance matrix
G, does not change the determinant in the expression for the characteristic function.
In particular, for the considered case of mixing two eigen-squeeze modes the characteristic function of the total-occupation statistics is
Its generating cumulants are easy to calculate as follows
where
denotes a gamma-function. They are exactly the same, up to an obvious common factor of 2, as the generating cumulants of a total occupation in independent eigen-squeeze modes with sinusoidal wave functions. Knowing these cumulants, one may restore in a simple way the central moments of the distribution as per comments after Equation (
68) and detailed discussion in [
36].
The probabilities of sampling
n atoms total in the considered excited states are easy to calculate expressing the determinant in the characteristic function in terms of the eigenvalues
of the renormalized covariance matrix
,
Recalling Equation (
54),
, and taking into account the relation
following from Equation (
53), we finally obtain a simple formula for the probability of sampling
n atoms total:
This result has been obtained for a pair of counter-propagating plane waves in [
36,
65]. Note, that the combinations of parameters
or
introduced in [
36] or [
65] via some algebraic manipulations are, in fact, nothing else but inverse eigenvalues of the renormalized covariance matrix
.
In fact, the total noncondensate occupation statistics do not inherit the sophisticated, related to the ♯P-hardness behavior of the joint probability distribution. However, it inherits the nontrivial property associated with the parity effect discussed above. While the first eigenvalue in Equation (
93) is positive for any values of the normal and anomalous correlators,
, the second eigenvalue
becomes negative when the absolute value of the anomalous correlator gets larger than the normal correlator,
. Then the probabilities of sampling an even total number of excited atoms are enhanced. For a particular choice of the observational basis consisting of standing plane waves (see
Section 5), this property may be interpreted as a simple consequence of the fact that both independent occupation numbers
and
have strongly suppressed probabilities to be odd. For the observational basis consisting of traveling plane waves the joint statistics has a strong correlation at
(see
Section 6), and the total-occupation probability
hits this diagonal only for even
n. In the vacuum state, when the anomalous correlator achieves its maximal possible absolute value,
, we have
, and the probabilities of odd total occupations completely vanish.
8. The Nature of the ♯P-Hard Complexity: Squeezing and Interference of Atom Sampling States via Their Interplay with Eigen-Squeeze Modes and Eigen-Energy Quasiparticles
Let us briefly overview the aforestated analysis of the ♯P-hard problem of atomic boson sampling from an interacting BEC gas trapped in a box. The analysis is conducted by means of a new approach based on the recently found Hafnian master theorem [
9,
48]. We intentionally chose a textbook quantum many-body model aiming to explain how to use the hafnian approach as a regular method for dealing with various ♯P-hard problems. We infer that an equilibrium BEC gas in a box with periodic boundary conditions is one of the simplest models that has the potential to demonstrate quantum supremacy over classical computing. It allows us to greatly simplify and clarify the general formulas and reveal an explicit analytical description of the mechanism leading to the ♯P-hardness of computing quantum properties of many-body systems.
The general theory is formulated in
Section 2,
Section 3 and
Section 4. The simple explicit formulas and their numerical illustrations are given in
Section 5,
Section 6 and
Section 7 which show the increasing complexity of the joint probability distribution of the occupations of excited atom states with the growing complexity of the unitary mixing used to form those excited atom states out of the eigen-squeeze modes (
44) (compare patterns in
Figure 3,
Figure 4,
Figure 5 and
Figure 6). Of course, while this unitary mixing occurs separately in different low-dimensional blocks of eigen-squeeze modes the joint probability distribution remains easily computable by analytical, recursive or numerical means. Only with increasing dimensionality of those mixing blocks of eigen-squeeze modes computing sampling probabilities requires exponential time and becomes ♯P-hard.
We find that there are two ingredients of the ♯P-hardness of the atomic boson sampling—the squeezing and the interference. Moreover, there are two corresponding unique entities existing in the BEC gas of interacting atoms, the eigen-squeeze modes and the eigen-energy quasiparticles, which are directly responsible for the above-mentioned squeezing and interference (see the schematic diagram in
Figure 1,
Section 4). The eigen-squeeze modes are the eigenvectors of the Hermitian factor of the multimode squeeze matrix. They are the elementary, intrinsic carriers of its eigenvalues—the single-mode squeezing parameters, Equation (
33). The quasiparticles are the eigenvectors of the Hamiltonian, Equation (
10), and are described by two-component wave functions, Equation (
11). They are the elementary collective excitations carrying quanta of collective energy—the eigenvalues of the Hamiltonian.
As such each eigen-energy quasiparticle lives completely independent of other quasiparticles, but it is a superposition of many two-component eigen-squeezed quasiparticles formed by the irreducible, eigen-squeezing Bogoliubov transformation
, Equation (
26), from the one-component eigen-squeeze modes, each of which is itself a superposition of many excited bare-atom wave functions (see the schematic diagram in
Figure 1). As a result, the excited bare-atom wave functions are persistently interfering with each other, and the joint probability distribution of their atom numbers turns out to be ♯P-hard for computing.
Only in some very special cases, the above-mentioned probability distribution can be computed faster than in exponential time. An example is the case when each bare-atom excited state chosen for atom number measurements by a multi-detector imaging system coincides with an eigen-squeeze mode which, at the same time, gives the same spatial profile for both components of the quasiparticle wave function. Then the joint probability distribution turns into a separable product of the occupation probabilities of the single eigen-squeeze modes. Each of those probabilities describes nontrivial squeezed-state statistics but is computable via Legendre polynomials as per Equation (
59).
Another more involved example is the case when the bare-atom excited states are chosen to be traveling plane waves. Then atomic boson sampling splits into independent samplings within separable blocks of two counter-propagating plane waves with the wave vectors
and
. The sampling probability distribution for each of such blocks shows statistics of two-mode squeezing, which is more complex, but again is easily computable via the ordinary hypergeometric function as per Equation (
77). Note, that the presence of squeezing makes the sampling statistics of Bose-atom numbers, even in the two-mode case, very different from and much more involved than the joint occupation statistics of purely interfering non-squeezed bosons, including a well-known two-mode Hong-Ou-Mandel statistics of just interfering bosons [
67,
68].
One more, extreme example is the case when the squeezing is absent, say, due to the absence of interatomic interactions, like in an ideal Bose gas, or due to the absence of the condensate, like in a classical gas above the critical temperature. In this case all single-squeezing parameters in Equation (
28) vanish. Hence, the hafnian in Equation (
25) determining the joint sampling probabilities reduces to the permanent of a positive matrix since the extended covariance-related matrix has vanishing anomalous-correlator blocks and one can employ Stockmeyer’s approximating algorithm [
2,
14,
26,
61] for such a permanent.
A pivotal key for understanding the origin of the ♯P-hardness of atomic boson sampling is provided by the irreducible Bloch–Messiah reduction of the Bogoliubov transformation into the product of three blocks,
, as per Equation (
26) and
Figure 1. Via its direct relation to the hafnian of the extended covariance-related matrix constituting the general result for the joint probability distribution in Equation (
25), the Bloch–Messiah reduction explicitly reveals the mechanism of squeezing and two mechanisms of interference responsible for the ♯P-hard complexity. The squeezing is attributed to the single-mode squeezing block
of the Bogoliubov transformation. The eigenvalues
of the Hermitian factor of the multimode squeeze matrix, that is, the squeezing parameters of the eigen-squeeze modes, form an irreducible resource of the system of many interacting atoms in the BEC trap. The existence of squeezing in the BEC gas is known since [
36]. The interference, controlled by the block
and its unitary
V as per Equation (
35), between the observational bare-atom excited states and the eigen-squeeze modes constitutes the first mechanism of interference. The interference between the eigen-squeeze two-component excitations and the eigen-energy quasiparticles, controlled by the block
and its unitary
W as per Equation (
38), constitutes the second mechanism of interference (see the schematic diagram in
Figure 1).
In other words, the interacting BEC gas in a trap has two naturally built-in, intrinsic interferometers associated with the two interference mechanisms disclosed above.
On this basis, we conclude that even if just one of the interference mechanisms is available for controlling parameters of sampling in a wide range, then the ♯P-hard complexity for the average case still exists. This is the case for the presented model of atomic boson sampling in a box trap with a uniform condensate for which the parameters of the quasiparticles and eigen-squeeze modes, including the unitary
W defining their interference block
, are almost fixed by a given trapping potential and couplings in Equation (
2) and cannot be varied in a wide range. So, a wide variability of the BEC parameters and Bogoliubov couplings provided by the multi-qubit BEC trap [
10] is useful, but not necessary for demonstrating computational ♯P-hardness and potential quantum supremacy of atomic boson sampling.
Mathematically, this ♯P-hardness is the property of the hafnian (or permanent) [
69,
70,
71] of the extended covariance-related matrix in Equation (
25) which, according to the Hafnian master theorem [
9,
48], determines the atom-number sampling probabilities. These probabilities are calculated as the Fourier series coefficients (
6) of the easy-to-compute characteristic function in Equation (
8). Thus, the ♯P-hardness is due to an intuitively obvious complexity of computing the multivariate Fourier integral in Equation (
6) for a sign-indefinite strongly-oscillating function (its analog is a lacunary or fractal function with an exponentially wide spectrum) [
9,
72].
The ♯P-hardness of computing the hafnian in Equation (
25) follows [
2,
12] from two known facts: (a) the Haar randomness of the unitary matrices yields the Gaussian randomness of the extended covariance-related matrix
and (b) computing the hafnian of a random Gaussian matrix is a ♯P-complete problem. The ♯P-complete problems constitute the top-level complexity class within the class of ♯P-hard problems. In virtue of Toda’s theorem [
49,
50], the solution of any ♯P-complete problem is reducible in polynomial time to the solution of any other problem in this class. Therefore, the multivariate Fourier integration can be viewed as the universal origin, or source, of the computational ♯P-hardness and potential quantum supremacy of the many-body quantum systems. The point is that the quantum many-body systems process the multivariate Fourier-series transform naturally, as a routine part of their life, that is in linear or, at least, polynomial time, while the classical simulators or computers can do this only in exponential time.
In fact, the multivariate Fourier transform described above reveals a certain duality of quantum and classical computational complexity. Both the characteristic function and its Fourier transform, that is its Fourier-series coefficients which constitute the joint probability distribution of sampling probabilities, contain full information about the atomic sampling statistics. However, the former, as any matrix-determinant function, can be easily calculated by classical computers while the latter cannot. There is no contradiction hidden in this statement since the calculation of the multivariate Fourier integral for the sampling probabilities requires, in a general average case, computing the characteristic function under the integral in the exponentially large number of points. First, it means that there is nothing mysterious in ♯P-hardness since it is just an ordinary property of multivariate integration. Second, it means that deriving the sampling probabilities from experimentally accumulated probabilities of atom-number samples, each of which is easily given by the BEC-gas quantum simulator almost in no time, requires an exponentially large number of samples (experimental runs) and, hence, an exponential time. Thus, only some special problems, such as a generation of strings of random numbers obeying the hafnian-based probability distribution (
25), are easy for the quantum many-body BEC-gas system but exponentially hard for classical computing.
9. Towards Experiments on Manifestations of ♯P-Hard Complexity and qUantum Supremacy of Atomic Boson Sampling
We emphasize that contrary to a widely discussed Gaussian boson sampling of noninteracting photons in a linear interferometer, the proposed atomic boson sampling does not require sophisticated synchronized external sources of bosons in squeezed states. The squeezing and interference of atom excited states, both of which are necessary for the computational ♯P-hardness of boson sampling, are self-generated even in an equilibrium BEC gas. Hence, the major limitation factor for achieving quantum supremacy via boson sampling in a deep linear interferometer, which is an exponential loss of photons due to scattering and absorption on coupling elements (beam splitters, phase shifters, etc.) during propagation through the interferometer, is not an issue for the atomic boson sampling.
Conceptually, the experiments on sampling are simple. In fact, the experiments on the statistics of the total occupation of excited atom states, that is, the total noncondensate occupation, has been successfully performed [
73]. In order to pioneer atomic boson sampling one just needs to split the noncondensate into fractions and to measure many times the atom occupation numbers in a preselected subset of excited wave functions by means of some multi-detector imaging system. Then, (a) to reconfigure detectors and/or trapping potential and other parameters of the BEC-gas in a trap and (b) to measure sampling statistics for such a unitary-transformed subset of excited wave functions and new system’s parameters, and so on. There is no need neither for any controllable non-equilibrium unitary-evolution processes typical for most quantum-computing experiments nor in the suppression of various concomitant processes of relaxation and decoherence. The quantum system of interacting atoms in a BEC trap just simulates its own equilibrium life which consists of persistent quantum-statistical fluctuations. It is described by the statistical operator that intrinsically involves properties that are ♯P-hard for computing.
The ♯P-hardness of computing the hafnian-based sampling statistics is a fundamental reason for, but of course not equivalent to or sufficient for quantum supremacy of the BEC gas in a trap over classical simulators with respect to the generation of the strings of random numbers obeying the predetermined joint probability distribution in Equation (
25). Whether a classical computer/simulator can or cannot provide the generation of such random numbers in polynomial time is an open question. Obviously, the quantum system of many interacting atoms in a BEC trap is in a privileged position because its equilibrium state naturally provides the required sampling statistics. As is always the case in the discussion of quantum supremacy, the very choice of the problem for simulation is intentionally unfair with regard to classical computing. Surely, the relation between classical and quantum computing is asymmetric.
The aforestated analysis of atomic boson sampling for the interacting atoms in the BEC trap, in particular, the result for sampling probabilities in Equation (
25), reveals its close similarity to the Gaussian boson sampling of noninteracting photons in a linear interferometer. This fact allows one to transfer an existing extensive analysis of the prospects and requirements for demonstrating quantum supremacy of photonic boson sampling [
3,
4,
5,
11,
12,
13,
14,
15,
16,
17,
18,
19,
20,
21,
22,
23,
24,
25,
26,
29,
30,
31,
32,
33] to the case of atomic boson sampling. For instance, two previously suggested schemes of photonic Gaussian boson sampling—the scheme that smuggles a random Gaussian matrix as a submatrix of the covariance matrix [
11] and another scheme that smuggles an arbitrary symmetric matrix [
12]—could be, in principle, imitated within the atomic BEC platform by a proper choice of the system parameters and unitary
V. Demonstration of the recently suggested bipartite protocol of photonic Gaussian boson sampling [
13] is also possible since the required balanced two-mode squeezing is naturally generated in the atomic BEC box trap for the pairs of degenerate counter-propagating plane waves
, see Equation (
35). Subsequent mixing of the first wave components of these pairs (
waves) via a unitary
and independent mixing of the second wave components of these pairs (
waves) via a different unitary
yield the bipartite protocol by involving a large-size,
matrix block
which could be set arbitrary in view of the theorem on the singular value decomposition of a matrix. However, such an analysis goes beyond the scope of the present paper.
We just note a recent analysis [
2,
11,
12,
13,
27,
74,
75] suggesting a possibility of achieving quantum supremacy in the regimes with the mean number of bosons per each of
M sampled modes (channels) on the order of
(dilute sampling regime) or even ∼1 (high-collision regime), rather than only in a deeply unitary-hiding regime with a very small mean occupation per mode, ∼
, as had been assumed previously [
1,
3,
6]. Two well-known recent experiments on photonic boson sampling [
30,
31] had been also implemented at the order-of-unity mean occupation per mode, featuring up to 113 or 219 photon detection events out of a 144 or 216-mode photonic circuit, respectively. Hopefully, the atomic boson sampling and other experiments based on such a BEC platform for studying computational ♯P-hardness in quantum many-body systems will become available soon. Especially promising in this regard could be experiments similar to the measurement of the full counting statistics of excited, noncondensed atoms in the momentum space after their release from a trap and subsequent free-fall expansion [
76,
77]. Remarkably, in these experiments, the detectors with a large quantum efficiency and single-atom resolution have been demonstrated.
As is discussed above and illustrated in
Figure 4,
Figure 5 and
Figure 6, a clear demonstration of a full complexity of the sampling probability patterns in the experiments with atomic BEC-gas requires a proper choice of the observational excited bare-atom states, a strong enough interatomic interaction (i.e., quantum depletion of the condensate) and low enough temperature so that the squeezing would be strongly pronounced and quantum statistics would not be hidden under thermal fluctuations.
Discussion of various closely related techniques for imaging of the local atom-number fluctuations and achieving nearly single-atom resolution in BEC-gas experiments can be found in [
38,
39,
67,
78,
79,
80,
81,
82,
83,
84].
An ultimate demonstration of quantum supremacy for the average case [
2,
3,
6] can be achieved only if one obtains access to a wide-range unitary mixing, interference of the eigen-squeeze modes via either the observational basis states (the unitary
V in Equation (
26)) or the quasiparticle states (the unitary
W in Equation (
33)). The former can be conducted solely by reconfiguring atom-number detectors for projecting onto the unitary mixed excited bare-atom states. Reconfiguration of the trapping potential and other parameters controlling the condensate profile and interatomic interactions provides control on the interference via both unitaries
W and
V. In other words, one needs to obtain control on either the first or second intrinsic interferometers built by nature in the interacting BEC gas which correspond to the first or second interference mechanisms revealed in
Section 4 (see
Figure 1) and VIII. However, such a full control is not necessary for pioneering experiments on atomic boson sampling. Obtaining nontrivial patterns (see
Figure 4) of the joint probability distribution of atom numbers for two counter-propagating traveling plane waves, for example, based on the full counting statistics accumulated in the available experiments [
76,
77], would be already a proof-of-principle demonstration of atomic boson sampling.
Compared to the experiments on fluctuations of the total noncondensate occupation [
36,
37,
73], the experiments on atomic boson sampling do not imply counting all noncondensed atoms. In this respect, the latter experiments are even simpler than the former ones since an exact separation of a relatively small fraction of noncondensed atoms from much more occupied condensate, that is, drawing a precise borderline between the two, is the main challenge for the former experiments. For implementing atomic boson sampling, it suffices to measure joint occupation statistics just for some excited atom states (or coarse-grained groups of them) all of which could be far from the condensate wave function in the momentum space or, more generally, in the functional space and, therefore, easily distinguishable from the condensate.