Towards the Simplest Model of Quantum Supremacy: Atomic Boson Sampling in a Box Trap

We describe boson sampling of interacting atoms from the noncondensed fraction of Bose–Einstein-condensed (BEC) gas confined in a box trap as a new platform for studying computational ♯P-hardness and quantum supremacy of many-body systems. We calculate the characteristic function and statistics of atom numbers via the newly found Hafnian master theorem. Using Bloch–Messiah reduction, we find that interatomic interactions give rise to two equally important entities—eigen-squeeze modes and eigen-energy quasiparticles—whose interplay with sampling atom states determines the behavior of the BEC gas. We infer that two necessary ingredients of ♯P-hardness, squeezing and interference, are self-generated in the gas and, contrary to Gaussian boson sampling in linear interferometers, external sources of squeezed bosons are not required.

We describe boson sampling of interacting atoms from the noncondensed fraction of Bose-Einstein-condensed (BEC) gas confined in a box trap as a new platform for studying computational ♯P-hardness of quantum many-body systems. In this case the theory becomes really simple and transparent. We calculate analytically the characteristic function and statistics of the excited-state atom occupations in the Bogoliubov approximation by means of the newly found hafnian master theorem and show their analogy to those of the Gaussian boson sampling of noninteracting photons in a linear interferometer. Importantly, due to interatomic interactions, the squeezing and interference of atom excited states, both of which are necessary for the computational ♯P-hardness of boson sampling, are self-generated in the gas and do not require neither sophisticated external sources of bosons in squeezed states nor controlled couplers, beam splitters and phase shifters needed for boson sampling in optical interferometers. On this basis, we discuss how to get manifestations of quantum supremacy of boson sampling over classical computing in the experiments with BEC gas.

I. INTRODUCTION: THE SIMPLEST QUANTUM MANY-BODY MODEL SHOWING ♯P-HARD COMPLEXITY
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] 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 8 . A multi-qubit BEC trap that could serve as a rich and, at the same time, convenient platform for experimental studies of various phenomena associated with atomic boson sampling has been discussed in 9 .
The purpose of the present paper is somewhat opposite to that of 9 . Namely, we aim at the simplest possible model of the BEC trap that would allow us to greatly simplify the general theory outlined in 8,9 and explicitly disclose the mechanisms behind the ♯P-hard computational complexity of quantum many-body systems. In 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 ♯Phardness of computing the joint probability distribution of the excited-atom occupations.
The content of the paper is as follows. In section II 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 sections III and IV 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, by means of the Bloch-Messiah reduction of the Bogoliubov transformation. In sections V through VII we illustrate increasing complexity of sampling probability patterns with changing the observational basis of excited atom states from the eigensqueeze modes to more and more involved unitary mixtures of them. Section VIII contains discussion of the ♯P-hardness of computing the joint atom-number probabilities in a gen-eral case of arbitrary sampling states. The focus is on the origin and manifestations of ♯P-hard complexity of atomic boson sampling due to interference and squeezing of sampling atom states via their interplay with the eigen-squeeze modes and eigen-energy quasiparticles.

II. QUANTUM STATISTICAL PHYSICS OF JOINT FLUCTUATIONS OF THE EXCITED ATOM NUMBERS
Consider a box trap of volume V = L 3 with periodic boundary conditions. Excited atoms are described by a field operator whereb l denotes an operator annihilating an atom in the bareatom excited state φ l (r) of the single-particle Hilbert space H . The sum ∑ l denotes a summation over a basis of excited atom states in a box trap, {φ l }, enumerated by an integer l = 0. A condensate wave function in the box trap corresponds to the zero integer l = 0 and is uniform, φ 0 = V −1/2 . We abide the Bogoliubov approximation and replace the condensate annihilation operator by a c-number,b 0 ≈ √ N 0 , where N 0 is a mean number of atoms in the condensate.
The Bogoliubov Hamiltonian in the basis {φ l } is given, up to an insignificant additive c-valued constant, by a quadratic form in the creation and annihilation operators 10,11 : It is written via a (2 × 2)-block matrix H which is built of matrices ε, ∆, and∆ as its blocks, applied to a 2-block column vector b † b , consisting of the column vectorb † on top of the column vectorb, and multiplied from the left by a 2-block row The bold-faced operatorb = {b l } T orb † = {b † l } T denotes a column vector of annihilation or creation operators in the {φ l } basis;ε =h 2 ∇ 2 /(2m) is the singleparticle energy operator for a bare particle of mass m in the box trap. The nabla symbol ∇ stands for the 3D-vector differential operator. The inter-atomic interaction is determined by a constant g = 4πh 2 a s /m via the s-wave scattering length a s .
The general notations for the field operator and Hamiltonian in Eqs. (1), (2) are convenient for an abstract general analysis presented below in sections III, IV, and VIII. For any particular choice of the basis of excited atom states {φ l } the field operator and Hamiltonian acquire more specific and transparent form. In a usual basis of plane traveling waves with the wave vectors k enumerated by an integer 3-vector j, we haveψ Hereinafter, we denote the creation and annihilation operators of bare atoms in the traveling-plane wave states e ikr √ V by symbolsâ † = {â † k | k = 0} andâ = {â k | k = 0}, respectively, as opposed to the above general-case operatorsb † andb.
Equilibrium quantum many-body statistics of atoms in the BEC trap at a temperature T is determined by the statistical operatorρ = e −Ĥ/T tr e −Ĥ/T . The atomic boson sampling implies sampling in accord with the joint probability distribution ρ({n l }) of the occupation numbers {n l } of the excited atomic states {φ l } in Eq. (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 tr{. . .} stands for a trace of an operator or matrix. The characteristic function (7) had been found analytically 8 via a determinant function which can be easily computed in a polynomial time, It 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 b † b and 2-block row vector (b,b † ) of the atom creation and annihilation operators, that is, all possible self-and inter-mode normal, b † lb l ′ , and anomalous, b † lb † l ′ , b lbl ′ , correlators. Each variable τ l of the characteristic function Θ({τ l }) appears in the matrix of variables Z in Eq. (8) twice -via the entry z l = e iτ l in each of the two identical diagonal blocks diag({z l }); the symbol ½ denotes the identity matrix. (Note that in the literature on Gaussian states 12 , the covariance matrix is often defined with a half anti-commutator, (b † lb l ′ +b l ′b † l )/2, replacing the normal product,b † lb l ′ , of the creation/annihilation operators, that adds a half identity matrix ½/2 to our covariance matrix G.) Consider the quasiparticle creation and annihilation operators, which constitute the column vectorsb † = {b † l } T and b = {b l } T and diagonalize the Hamiltonian in Eq. (2), The quasiparticle eigen energy is denoted as E l and will be given in Sect. V. The field operator of excited atoms in Eq. (1) can be represented also via the quasiparticle operators as a sum of both annihilation and creation quasiparticle operators: The functions u l (r) and v * l (r) constitute the two-component wave function of the l-th quasiparticle. Canonical Bose commutation relations for the quasiparticle operators,b lb † l ′ − b † lb l ′ = δ l,l ′ , where δ l,l ′ is the Kronecker delta, imply the following normalization of the two components, u l and v * l , of each basis' quasiparticle wave function Suppose one chooses plane traveling waves e ikr √ V | k = 0 as the bare-atom excited-state basis. Then the Hamiltonian in Eq. (5) acquires a diagonal form with the canonical Bogoliubov spectrum of eigen energies, for the quasiparticle creation and annihilation operatorsâ † = {â † k | k = 0},â = {â k | k = 0} related to the corresponding bare- Here we introduced the amplitudesū k andv k , 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 Eq. (1) acquires the following explicit form In any case, the quasiparticles are completely independent on each other. Thus, the correlations between the quasiparticle creation/annihilation operatorsb † ,b, analogous to the bareatom correlations in Eq. (9), are given by the (2 × 2)-block diagonal matrix of the thermal occupations of quasiparticles Employing the canonical Bose commutation relationŝ b lb † l ′ −b † lb l ′ = δ l,l ′ 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, by means of the Bogoliubov transformation which relates the uncorrelated quasiparticles with the squeezed and interfering bare-atom excitations. It is described by the (2 × 2)-block symplectic matrix 8R and its inverse matrix R =R −1 as follows Below we assume that some finite number M of orthogonal excited states (i.e., atom wave functions) {φ l } 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 H . Their normal and anomalous correlations are given by the corresponding (2M × 2M)-submatrix of the covariance matrix in Eq. (9). Suppose the preselected wave functions are coupled via the Bogoliubov Hamiltonian only between themselves, that is, their off-diagonal couplings ∆ ll ′ ,∆ ll ′ , and ε ll ′ in Eq. (2) with any wave functions outside the preselected subspace are zero. This is easy to achieve since in the uniform box trap the swave scattering in the Bogoliubov Hamiltonian couples just plane waves with opposite wave vectors k and −k. This simplifies the analysis. All bold-faced vectorsb and the like are reduced to M-dimensional vectors. All (2 × 2)-block matrices (including the Hamiltonian, H, Bogoliubov, R, quasiparticle occupation, D, covariance, G, and variable, Z, matrices) are reduced to the (2M × 2M)-matrices containing corresponding (M × M)-blocks such as ∆,∆, ε, A, B, etc. For the sake of simplicity, we'll denote any such matrix of a finite, reduced dimension 2M × 2M or M × M by the same symbol that stands for its infinite-dimensional counterpart.

III. THE HAFNIAN MASTER THEOREM AND SAMPLING PROBABILITIES
The hafnian master theorem recently found in 8,13 provides the most convenient and powerful regular method for the analysis of the atomic, gaussian boson sampling and other problems associated with the ♯P-hard computational complexity. The point is that it directly reduces a ♯P-hard-for-computing quantity in question to a rigorously defined, canonical mathematical function -a matrix hafnian (or its particular case -a matrix permanent) computation of which belongs to the hardest, ♯P-complete class of computational complexity. In accord with the famous Toda's theorem 14,15 , it means that computing the hafnian and using it as an oracle is enough for polynomialtime reduction of every other ♯P-complete or ♯P-hard problem to an easy, polynomial-time computational problem.
In our case the hafnian master theorem gives an explicit Fourier series (that could be viewed also as a Taylor expansion) of the characteristic function in Eq. (8), It is expressed via the hafnian of an extended covariancerelated (2n × 2n)-matrixC({n l }), which has a dimension determined by the total atom number n = ∑ l n l in the sample of occupations of the preselected excited states {φ l } and is built from the covariance-related matrix One has to replace the l-th and (M + l)-th rows with the n l copies of the l-th and (M + l)-th rows, respectively, and then the l-th and (M + l)-th columns with the n l copies of the lth and (M + l)-th columns, respectively. The matrix P just permutes the diagonal and off-diagonal blocks of the (2 × 2)block matrix G(½ + G) −1 in a way appropriate for the hafnian. Employing the identity following from the definition of the characteristic function, we get an explicit analytical formula for the joint probabilities of the atom numbers in Eq. (6) via the hafnian as follows

IV. TWO INGREDIENTS OF QUANTUM SUPREMACY: SQUEEZING AND INTERFERENCE
The result in Eq. (24) explicitly shows that the ♯P-hard complexity and, hence, quantum supremacy of atomic boson sampling come from the ♯P-hardness of computing the hafnian of the extended covariance-related matrixC({n l }) that depends only on the nontrivial Bogoliubov-transformation matrix R and mean quasiparticle occupations constituting the diagonal matrix D in Eq. (18). Thus, the mystery of quantum supremacy is encoded in the structure of the Bogoliubovtransformation matrix in Eq. (20) and can be revealed via its unique Bloch-Messiah representation 19,20 It involves two unitary matrices, W and V , as well as the diagonal matrix of M single-mode squeezing parameters The Bloch-Messiah reduction (25) can be derived starting from the polar decomposition of the blocks A = U cosh r and B = Ue iθ sinh r of the matrix R in Eq. (20) via the unitaries U and Ue iθ as well as the Hermitian amplitudes determined by the M × M multimode squeezing matrix r = (r ll ′ ) by means of the hyperbolic functions cosh r and sinh r, respectively. In quantum optics [19][20][21][22][23][24] this polar decomposition is expressed usually as an effective evolution resembling the Hamiltonian one, e −iĤt , but under the action of the unitary multimode squeeze,Ŝ, and rotation,Φ, operators: The multimode squeeze and rotation operators, are determined by a symmetric matrix S = e iθ r (built of the Hermitian matrices r, θ ) and a Hermitian matrix Φ = Φ † generating the unitary U = e iΦ , respectively. The relation S = e iθ r can be viewed as the definition of the multimode squeeze matrix r.
Both above-mentioned approaches yield the same polar decomposition of the Bogoliubov transformation: Its further reduction yields the Bloch-Messiah representation in Eq. (25) by means of the singular value decomposition and the relations cosh r = W cosh Λ r W † , sinh r = W sinh Λ r W † , which follow from the diagonal representation of the multimode squeezing matrix The unitaries W and V are chosen to satisfy the so-called rotation condition emphasized in 20  Eq. (25) describes the overall Bogoliubov transformation as a sequence of three (2 × 2)-block transformationsR WRrRV from the creation/annihilation operators in the observational basis of the excited atom states preselected for sampling measurement to the creation/annihilation operators of M completely independent quasiparticles which correspond to the diagonalized Hamiltonian in Eq. (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 quasiparticle creation and annihilation operators. The first Bogoliubov transformationR V corresponds to the unitary rotation V * of the basis of the bare-atom excited states {φ l | l = 1, ..., M}, prescribed for measurement of atom numbers by means of multi-detector imaging in the process of sampling, into the basis of the single-squeeze modes {ϕ l } associated with the eigenvectors {w l } of the multimode squeezing matrix r in Eq. (31) and explicitly expressed below in Eq. (35) via the quasiparticle wave functions.
The second Bogoliubov transformationR r corresponds to a modification of the Hamiltonian by introducing the terms beyond the resonant-wave approximation (non-RWA terms), b † lb † l ′ andb lbl ′ , which create or annihilate a pair of excited atoms, and, thus, converts the state of the system into the squeezed state in which each single-squeezed mode acquires the corresponding nontrivial squeezing parameter r l ≥ 0.
The third Bogoliubov transformationR W converts the creation, {b ′ l ′ † }, and annihilation, {b ′ l ′ }, operators of the singlesqueeze modes, obtained on the second step, into the quasiparticle creation, {b † l }, and annihilation, {b l }, operators: According to the form of the quasiparticle field operator in Eq. (11), it means two simultaneous unitary rotations (under the action of one and the same unitary matrix W † ) of the single-squeeze modes {ϕ l } into the wave functions {u l | l = 1, ..., M} and {v * l | l = 1, ..., M}, which according to Eq. (11) constitute the bases of the first and second components of the two-component functional space of the quasiparticle wave functions (u(r), v * (r)), Each single-squeeze mode ϕ l ′ (r) owns the single-mode squeezing parameter r l ′ ≥ 0 and is not subject to inter-mode squeezing with other eigen-squeeze modes.
The existence of such a unique unitary W , simultaneously converting the basis wave functions u l (r) and v * l (r) of both components of the two-component functional space of quasiparticles into the basis wave functions u ′ l ′ (r) = (cosh r l ′ )ϕ l ′ (r) and v ′ * l ′ (r) = (sinh r l ′ )ϕ l ′ (r) which are equal to the same eigensqueeze modes ϕ l ′ (r) just multiplied by the different constant factors cosh r l and sinh r l , respectively, is a nontrivial and important property of the Bogoliubov transformation. It is a consequence of the symplectic property 8 , of the Bogoliubov transformation and highlights the fact that the two components of the quasiparticle eigen function (u l (r), v * l (r)) are not independent, but, on the contrary, are deeply inter-correlated.
Both the first,R V , and the third,R W , Bogoliubov transformations introduce unitary interference, that is entanglement, into the quantum many-body state (statistical operatorρ) of M excited atom modes chosen for atomic boson sampling. However, they are controlled by different means. The unitary W of the third one is controlled by the trapping potential and parameters of the Hamiltonian in Eq. (2), that is, by cou-plings∆ ll ′ , ∆ ll ′ , ε ll ′ and atom interactions. The unitary V of the first one can be controlled simply by a reconfiguration of the multi-detector imaging system for sampling different excited atom states. In the particular case of a box trap with a uniform condensate, discussed in the present paper, the third Bogoliubov transformationR W is almost completely fixed. Yet, the first Bogoliubov transformationR V provides an 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 Eq. (24) on average. Such complexity on average 3,16 is crucially important for possibility of demonstrating quantum supremacy of atomic boson sampling.
Thus, the Bloch-Messiah reduction in Eq. (25) unambiguously specifies two preferred bases: (i) the basis of the quasiparticle operators in Eq. (11) ensuring the diagonal form of the Hamiltonian (10) and (ii) the basis of the eigen-squeeze single-particle excited states in Eq. (35) diagonalizing the multi-squeezing matrix, Eq. (31). Moreover, the Bogoliubov transformation in the Bloch-Messiah representation explicitly relates both above-mentioned bases with the observational basis of excited bare-atom states {φ l | l = 1, ..., M} which can be arbitrarily selected by a reconfiguration of atom detectors.
The interference between those bases, appearing in the statistical operator (state) of the atomic multimode system via the unitaries W,V and leading to the entanglement of different 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 transformationR r determined exclusively by the eigenvalues r l of the multimode squeezing matrix r. If it was given by the identity matrix with all singlesqueezing parameters equal zero, r l = 0 ∀l = 1, ..., M, then the hafnian in Eq. (24) would reduce to the permanent of a positive-definite matrix, which could be approximated by the Stockmeyer's approximating algorithm 17,18 in a class simpler than ♯P. A physical mechanism of squeezing can be seen from Eq. (11) for the quasiparticle field operator. Any quantum or thermal fluctuation associated with a quasiparticle excitation implies creation of a bare atom in an excited state given by some superposition of the basis wave functions {u l | l = 1, ..., M} and simultaneous annihilation of a bare atom in another excited state given by another superposition of the basis wave functions {v * l | l = 1, ..., M}. Below we illustrate and discuss the interplay and manifestations of these two ingredients of quantum supremacy in atomic boson sampling for different choices of the observational basis of excited bare-atom states {φ l | l = 1, ..., M}.

V. FACTORIZATION 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 computational complexity of joint atom-number probabilities in the atomic boson sampling is achieved if one chooses the eigen-squeeze modes in Eq. (35) as the bare-atom excited states for atom-number measurements by the multidetector 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's choose them to be sin(kr) and cos(kr) spatial modes with a wave vector k. In fact, any orthogonal (but not any unitary) transformation of a pair of energy-wise degenerate single-particle excited states would lead to the same atom-number sampling statistics.
So, the excited atoms are described by the field operator written via operatorsŝ k ,ĉ k annihilating a bare atom in the corresponding sinusoidal states with a wave vector k = 2πj Hereinafter, an apostrophe in the symbol of the sum ∑ k =0 ′ means summation over nonzero wave vectors modulo multiplication by −1, that is, the integer vector j ≡ { j x , j y , j z } is running over a half of a three-dimensional (3D) integer space, 3 , 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 k = 2πj/L and −k = −2πj/L, 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 k and −k only one of the terms belongs to the sum. (Note that the sum in Eqs. (1), (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 separately operators related to sin(kr) or cos(kr) modes, It is in contrast to the Hamiltonian in Eq. (5) in which the modes e ikr and e −ikr of the running-plane-wave basis e ikr √ V k = 0 were coupled to each other. The Hamiltonian in Eq. (39), compared to the Hamiltonian in Eq. (5), has the same bare-atom energies ε k and normal overlapping integrals ∆ kk ′ , but different, now symmetric anomalous overlapping integrals In this case, the quasiparticle annihilation operators {ĉ k ,ŝ k } are given by a similar Bogoliubov transformation via the corresponding bare-particle annihilation and creation operators: The quasiparticle operators diagonalize the Hamiltonian (39), and turn the field operator (37), annihilating bare atoms, into a sum of both annihilation and creation quasiparticle operators: (43) In both, exponential and sinusoidal, bases, the quasiparticle energies E k and factorsū k ,v k are the same, given in Eqs. (13), (15). The energy E k is larger than the bare-atom energy ε k .
Thus, in the observational basis of the eigen-squeeze modes (38) the sampling joint probability distribution factorizes into the product of probabilities for single sin(kr) or cos(kr) modes. Each of such single-mode-squeezed probabilities can be calculated analytically by means of the general solution in Eq. (24). We just need to appreciate the fact that in this case the general Bloch-Messiah representation of the Bogoliubov transformation in Eq. (25) is reduced to the simple independent blocks with trivial unitary partsR This is the Bogoliubov single-mode squeezing in the simplest, pure form with a single-mode squeezing parameter r k ≥ 0. The corresponding single-mode covariance matrix is HereÑ k , η and α stand for the average number of quasiparticles, normal correlator and anomalous correlator of the single sinusoidal mode, respectively. The anomalous correlator is negative and its absolute value is bounded from above by the inequality in Eq. (45) which is equivalent to the obvious inequalityÑ k ≥ 0. The maximum value of the anomalous correlator, max |α| = η(1 + η), is achieved whenÑ k = 0, that is, when the quasiparticles are in the squeezed vacuum state. The relations between the average quasiparticle occupatioñ N k , the single-mode squeezing parameter r k and the normal and anomalous correlators η and α are illustrated by Fig. 1. Any fixed value of the normal correlator which is equal to the mean number of atoms in the eigen-squeeze mode, η = n , corresponds to a circle on the plane Ñ k , |α| with the origin at the point (−1/2, 0): (Ñ k + 1/2) 2 + |α| 2 = (η + 1/2) 2 . The eigenvalues of this covariance matrix depend on the singlemode squeezing parameter r k and the mean quasiparticle oc-cupationÑ k , λ 1,2 = η ± |α| =Ñ k e ±r k + (e ±r k − 1)/2.
The atom-number probability in a sinusoidal mode is determined, as per Eq.
1. An average quasiparticle occupation,Ñ k , (red solid curves) and the single-mode squeezing parameter, r k , (blue dashed curves) vs the absolute value of the anomalous correlator | ŝ kŝk | = |α| for different fixed values of the normal correlator ŝ † kŝ k ≡ η. Red solid curves representing the average quasiparticle occupation are, in fact, the arcs of a circle with a radius η + 1/2 centered at (0, −1/2). C = PG (k) ½+G (k) −1 which has the following explicit form: It is convenient to denote the entries α ′ , η ′ of the matrix C = PG (k) ½ + G (k) −1 by adding apostrophe to the symbols α, η denoting the entries of the covariance matrix G (k) . The eigenvalues λ ′ 1,2 of the renormalized covariance matrix G (k) ½ + G (k) −1 look similar to the eigenvalues λ 1,2 of the covariance matrix G (k) , namely, Note that the maximal value of the anomalous correlator, max |α| = η(1 + η), is achieved when η ′ becomes zero, η ′ = 0. Now we calculate the atom-number probabilities. The result is immediately given by the hafnian master theorem (24), where the extended covariance-related matrixC(n) involves the n × n matrix J n×n with all entries equal unity.
In the absence of anomalous correlations, when α ′ = 0, the eigen-squeeze modes are in a non-squeezed thermal state and the hafnian in Eq. (49) is reduced to the known permanent of the unity-matrix J n×n , that is, hafC(n) = per η ′ J n×n = n!(η ′ ) n . In this case the atom-number probability distribution is reduced to the simple exponential law The presence of a nonzero anomalous correlator α ′ introduces squeezing and makes this distribution nontrivial. Fortunately, the hafnian in Eq. (49) in the case of an arbitrary α ′ is easy to calculate via the known recursive relation for hafnians 25 . Performing two recursive steps and excluding the hafnian of an auxiliary matrix, one finds the following secondorder recursive formula It starts from the plain values of the hafnian hafC(1) = η ′ at n = 1 and hafC(0) = 1 at n = 0 (as per convention in the hafnian master theorem). Comparing it with the well-known recursive relation for the Legendre polynomials 26 P n : we at once reduce the hafnian to the Legendre polynomials: (52) Thus, the atom-number probabilities for sampling the occupation of the single-mode-squeezed standing plane wave are The second equality is due to substitution (η ′ ) 2 − |α ′ | 2 = η 2 − |α| 2 (η + 1) 2 − |α| 2 . Note that all Legendre polynomials P n 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 Eq. (49) can be calculated in a straightforward way via elementary combinatorics as follows Here the symbol ⌊n/2⌋ stands for the largest integer less than or equal to n/2, and the essential argument is In the extreme case of zero mean occupation of quasiparticles, that is when η ′ = 0, or |α| = η(1 + η), the recursive relation becomes purely two-step, Eq. (56) yields the result It leads precisely to the formula for the well-known occupation probabilities for a mode in the squeezed vacuum state 28 , 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,Ñ k = 0, means that the mode is in the vacuum Fock state |0 0| with respect to the quasiparticle operators. Typical occupation probability patterns given by Eq. (53) for sampling from eigen-squeeze bare-atom excited states (38) are illustrated in Fig. 2. In the range |α| ≤ α cr , where the absolute value of the anomalous correlator is less than the critical value α cr = η, an increase of the anomalous correlator α leads to growing probabilities of occupations which are less than some number n 0 and suppressing probabilities of occupations which are larger than that number n 0 . Such a behavior drastically switches to a completely different behavior when the anomalous correlator exceeds the critical value, |α| > α cr = η, 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 probabilities of odd occupation numbers start to drop down. At the maximal possible value of the anomalous correlator, max |α| = η(1 + η), which corresponds to the squeezed vacuum state,Ñ k = 0, all probabilities for even occupation numbers become zero.
The characteristic function Θ(z) of this probability distribution gives explicit information on the moments as well as ordinary, κ m , and generating,κ m cumulants of the sampling probabilities via the Taylor series of its logarithm: It is given by the general result in Eq. (8) as follows We see that it is proportional to the well-known generating function of the Legendre polynomials 26,30 . In Eq. (61) we provide two expressions for it. The first one is convenient for the cumulant analyzis, while the second one -for the calculation of probabilitites via the hafnian master theorem. In particular, the exact result for the generating cumulants is where λ 1,2 , Eq. (46), are the eigenvalues of the covariance matrix (45). The value of the first generating cumulantκ m=1 yields the mean occupation n ≡ ŝ † kŝ k = η whose thermal and quantum contributions are equal to (e E k /T − 1) −1 cosh 2r k and sinh 2 r k , respectively. Adding the value of the second generating cumulant, we immediately get the standard deviation, σ 2 ≡κ m=1 +κ m=2 = η 2 + η + |α|. The higher moments of the atom-number probability distribution can be found in a similar way. It is worth noting that the presence of nontrivial anomalous correlator makes this characteristic function in Eq. (61) very different from the one, θ IBG (z) = (1 − (z − 1)η IBG ) −1 , characterising the atom-number sampling statistics in the limit of non-interacting, ideal Bogoliubov gas (IBG) when all excited atom modes remain in the non-squeezed state.

VI. 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's see now how the sampling statistics is 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 {φ l } (which is the set of bare-atom excited states for atom-number measurements) to be the traveling plane waves ∝ e ±ikr , Eq. (3). This is the simplest possible unitary mixing of the eigen-squeeze modes (38) 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â k ,â −k which differ from the annihilation operatorsŝ k ,ĉ k 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 (25) e +ikr e −ikr .
(64) Obviously, the introduced interference proceeds independently within each (2 × 2)-block of the atom excited states with the wave vectors k and −k. 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 (â † k ,â † −k ,â k ,â −k ) T . The related covariance matrix defined in Eq. (9) is equal to It immediately follows from the general formula in Eq. (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 (45) written for the 4-vector (ĉ † ,ŝ † ,ĉ,ŝ) T of the creation and annihilation operators of the eigen-squeeze modes in the form of the (4 × 4)-matrix The normal and anomalous correlators in the exponentialfunction basis remain the same as they were in the sinusoidalfunction basis and are equal to η and α, respectively (see Eq. (45) in the previous section).

(68)
Here the hafnian matrix function is applied to the extended covariance-related matrixC(n 1 , n 2 ) which is a block matrix. The block J n 1 ×n 2 is a n 1 × n 2 matrix with all entries equal unity, 0 J n 1 ×n 2 is a zero matrix of the same size. The determinant in the denominator has been calculated explicitly as follows: det ½ + G (k) exp = (η + 1) 2 − |α| 2 2 . The matrixC(n 1 , n 2 ) consists of blocks of unequal dimensions, however its total dimension 2n 1 + 2n 2 is even. Such a hafnian hafC(n 1 , n 2 ) is easy to compute. The simplest way to do so is via consecutive swapping the 2-nd and 4-th blockrows 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 a well-known identity valid for any matrix A. As a result, we calculate the required hafnian as follows hafC(n 1 , n 2 ) = per η ′ J n 1 ×n 1 α ′ J n 1 ×n 2 α ′ J n 2 ×n 1 η ′ J n 2 ×n 2 = min(n 1 ,n 2 ) Here the permanent has been calculated explicitly by simple combinatorial means and binomial coefficients C k n = n! k!(n−k)! . Indeed, one just calculates the number of permutations of n 1 + n 2 elements which swap k elements between the block of the first n 1 elements and the block of the last n 2 elements, which corresponds to one summand (η ′ ) n 1 +n 2 −2k |α ′ | 2k . The obtained sum is an ordinary (Gaussian) hypergeometric function 2 F 1 with integer parameters which is a polynomial. It could be also expressed in terms of the Jacobi polynomial, however, the hypergeometrical representation is more convenient. The argument of the hypergeometric function in terms of E k and r k is given in Eq. (55).
The hafnian in Eq. (68) can be calculated also directly from its combinatorial definition 25 . Each product of n 1 + n 2 entries contributing to the sum in the hafnian's definition is encoded by a division of the matrix dimension 2n 1 + 2n 2 into n 1 + n 2 unordered pairs. Elements from the first group of n 1 elements could be paired either to the elements from the second group of n 2 elements, which corresponds to picking an α ′ entry from the non-diagonal block, or to the elements from the third group of n 1 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 n 1 and n 2 elements, respectively), as well as swaps k elements between the third and fourth blocks. The rest of elements should be swapped between the first and third blocks, or between the second and fourth blocks.
Let k ≤ min(n 1 , n 2 ) be a number of pairings between the first and second groups. It is also the number of pairings be-tween the third and fourth groups. There are C k n 1 C k n 2 k! 2 alternative ways to choose a pairing between the first and second groups, and each variant corresponds to the multiplier α k accumulated from a non-diagonal block ofC. The same holds for pairing between the third and fourth groups. Both the first and the third groups have n 1 − k unpaired elements in rest. There are (n 1 − k)! alternative ways to pair them correspondingly, each corresponds to the multiplier (η ′ ) n 1 −k . Similarly, both the second and the fourth groups have n 2 − k unpaired elements in rest. There are (n 2 − k)! alternative ways to pair them correspondingly, each corresponds to the multiplier (η ′ ) n 2 −k . Here comes the result stated above in Eq. (70).
Combining Eqs. (68) and (70), we get 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: An equivalent result may be derived via a straightforward differentiation of the characteristic function, but the corresponding calculations are more cumbersome 30 .
In the particular case of an ideal, non-interacting BEC gas, when the anomalous correlator equals zero, α = α ′ = 0, and the excited atom states are not squeezed, the calculation of the joint probability distribution in Eq. (68) becomes elementary since the hafnian of the corresponding, extremely degenerate matrixC(n 1 , n 2 ) is reduced to the product of two permanents: ρ n 1 ,n 2 = per(η ′ J n 1 ×n 1 ) per(η ′ J n 2 ×n 2 ) n 1 ! n 2 ! (η + 1) 2 = 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 Eq. (50).
In the presence of nonzero anomalous correlators the joint probability distribution could not be factorized anymore. This is due to appearance of the intra-modal squeezing. As is illustrated in Fig. 3, increasing absolute value of the anomalous correlator results in enhanced correlations between the sampled atom numbers n 1 and n 2 and growing up probability of the completely entangled occupation states n 1 = n 2 .
The leading term in the asymptotics of the joint probability distribution in Eq. (71) when the anomalous correlator approaches its extremum max|α| = η(1 + η) can be calculated directly from the hafnian formula in Eq. (68) by means 3. Typical dependence of the joint probability distribution ρ n 1 ,n 2 for sampling n 1 and n 2 atoms in the two counter-propagating plane waves on the absolute value of the anomalous correlator, |α|. The normal correlator eta (that is, the average number of atoms per one mode) is set to be η = 10. The anomalous correlator values are (a) α = 0 (which corresponds to an ideal, non-interacting gas), (b) |α| = 9.5, (c) |α| = 10, and (d) max |α| = η(1 + η) (which corresponds to the two-mode squeezed vacuum state). Enlarging |α| leads to growing correlations between the stochastic variables n 1 and n 2 . Note that scaling for values of ρ n 1 ,n 2 changes significantly from (a) to (d).
It coincides with the contribution of the highest order monomial in the hypergeometric-function polynomial in Eq. (71) and yields purely diagonal joint probability distribution in which the probability ρ n 1 ,n 2 =n 1 exponentially decreases with increasing atom numbers n 1 = n 2 .
A well-known case of the squeezed vacuum state corresponds to the point at the very expremum, |α| = max|α|, η ′ = 0, when the asymptotics (73) is reduced to a well-known result for the two-mode squeezed vacuum state 28 ρ n 1 ,n 2 = δ n 1 ,n 2 cosh −2 r k tanh n 1 +n 2 r k . (74)

VII. NONTRIVIAL EFFECT OF INTERFERENCE ON THE ATOMIC BOSON SAMPLING IN THE BASIS OF ANY TWO-MODE-SQUEEZED UNITARY-MIXED DEGENERATE STANDING PLANE WAVES
Let's look at further complication of the sampling probability patterns due to 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 sin(kr), cos(kr) via an arbitrary unitary mixing described by the unitary matrix V and corresponding Bogoliubov transformation from the annihilation operatorsb k,1 ,b k,2 of these two excited atom states to the annihilation operatorsŝ k ,ĉ k of the two eigen-squeeze modes: Contrary to the particular case of two counter-propagating plane waves in Eq. (63), now the unitary V involves two arbitrary real-valued angles ξ and β . The matrix has only two free parameters instead of four parameters, which allow one to encode an arbitrary unitary matrix in the following manner: U = x y −e iγ y * e iγ x * , |x| 2 + |y| 2 = 1, γ ∈ R. Here the number of parameters is halved since the gauge phase factors of the two excited wave functions chosen to constitute the observational basis are not important. In fact, the whole set of observational bases which correspond to physically different joint statistical distributions of occupation numbers is parametrized by matrix V in the form of Eq. (75) with ξ ∈ [0, π/4] and β ∈ [0, π/2]. Thus, the wave-functions of the selected measurement basis, are the superpositions of standing and traveling plane waves. The covariance matrix for operatorsb k,1 ,b k,2 could be easily obtained from the covariance matrix (66) formed by the operators of the eigen-squeeze sine and cosine modes via the following transformation generated by the unitary matrix V : (77) 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 Eq. (47)) are the same as in sections V, VI.
Note that if the involved 2 × 2 symmetric matrix V T V is proportional to the identity matrix and, hence, the same is true for the complex conjugated matrix V † V * , computing probabilities is reduced to the simplest case of two independent eigensqueeze modes solved in section V. 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 V = V * and V T V = ½.
The sampling probabilities are given by the hafnian master theorem (24), ρ n 1 ,n 2 = hafC(n 1 , n 2 ) In the present case of a general-type unitary mixing the hafnian in Eq. (79) strongly depends on the variable entries of the nontrivial symmetric matrices V T V and V † V * in Eq. (78). Let's 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 2n 1 + 2n 2 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 get the hafnian as the following sum hafC(n 1 , n 2 ) = (n 1 !n 2 !) 2 It runs over only such subsets of indices which yield integer numbers under factorials in the denominator of Eq. (81); otherwise, the summand should not be included in the sum. In other words, all numbers n 1 − d 1 − j 1 and n 2 − d 2 − j 1 and n 1 − d 1 − j 3 and n 2 − d 2 − j 3 should be even. Fig. 4 exemplifies how correlations between occupation numbers n 1 and n 2 arise while one varies the unitary V making the observational-basis bare-atom excited states more entangled. Occupation statistics of two eigen-squeeze modes is separable (see panel (a)). While one switches the measuring 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 (n 1 , n 2 ). Switching to fully traveling waves, which is the case considered in the previous section, ultimately leads to a ridge extending along the diagonal n 1 = n 2 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 c 1 , c 2 , and c 3 . 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, max |α| = η(1 + η). Otherwise, for smaller |α| and larger η, the nontrivial correlation pattern appeared to be masked by thermal contributions.
All correlation patterns on the plane of stochastic variables n 1 and n 2 shown in Fig. 4 are symmetric because we been finally reduced to the associated Legendre polynomial with n 1 and n 2 determining its parameters. Note that in the squeezed vacuum state the probability of sampling the atom numbers n 1 and n 2 of different parity is always zero. This is an immediate consequence of the fact that the hafnian of the matrix of odd size n 1 + n 2 is identically equals zero.
The squeezed vacuum state corresponds to the most contrasting correlation pattern of the joint probability distribution. We illustrate this thesis by Fig. 5 in the case of the vacuum state of the system with the normal correlator η = 10 and the unitary mixing matrix V = 1 2 The wave functions of the selected measurement basis represent some mixtures of standing and traveling plane waves. On the plane (n 1 , n 2 ) of sampling outcomes one can clearly see the dedicated directions along which the probabilities are relatively large. This global landscape is also crucially affected 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 n 1 = n 2 = 4). 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 Fig. 5 and Fig. 4, panel (c). Both of them refer to exactly the same matrix V and observational basis, while the anomalous correlator |α| in Fig. 5 is larger. The dedicated directions of more probable outcomes are still the same and can be clearly seen, but they turn into gently sloping hills, and the abrupt check-mate pattern has disappeared. One should note, however, that thermal excitations don't simply raise the background and the pattern of vacuum-state statistics doesn't just draw into them. 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 Fig. 3 and Fig. 4, representing the joint occupation statistics for two counter-propagating plane waves for the states with no and few quasiparticles, respectively. While in Fig. 3(d), corresponding to the vacuum state, there is a wellpronounced and sharp diagonal ridge at n 1 = n 2 , the thermal excitations make it wider and enlarge, first of all, the probabilities of the adjacent states, |n 1 − n 2 | = 1, 2. Such behavior in general could be described by Eq. (81) for hafnians, since they essentially allow one to employ a perturbation approach with respect to small value of the entries η ′ (which are exactly zero for the vacuum state).
Importantly, while controlling the observational basis via the unitary matrix V may lead to nontrivial joint statistics of the occupation numbers with widely ranging correlation patterns as described above, it doesn't affect the total-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 Eq. (79) involving complicated hafnians.
The distribution of the total number of atoms in any selected set of excited states is easier to described via the characteristic function in Eq. (8) with all arguments set to be equal, z j = z. (Setting some groups of arguments equal to each other corresponds to calculating a coarse-grained statistics.) Then, the matrix of variables Z in Eq. (8) turns into the scaled identity matrix z½, which commutes with any matrix V . Thus, the unitary matrix, switching the observational basis and transforming the covariance matrix G, doesn't 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 totaloccupation statistics is (84) Its generating cumulants are easy to calculate as follows 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 on Eq. (62) and detailed discussion in 29 .
The probabilities of sampling n atoms total in the considered excited states is easy to calculate expressing the determinant in the characteristic function in terms of the eigenvalues λ ′ 1,2 of the renormalized covariance matrix Recalling Eq. (48), λ ′ 1,2 = η ′ ± |α ′ |, and taking into account the relation |α ′ | det ½ + η α α η = |α| following from Eq. (47), we finally get 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 29,30 . Note that the combinations of parameters Z(±A k ) or Y ± introduced in 29 or 30 via some algebraic manipulations are, in fact, nothing else but inverse eigenvalues of the renormalized covariance matrix G (k) 1 + G (k) −1 .
In fact, this total-occupation statistics does 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 Eq. (87) is positive for any values of the normal and anomalous correlators, 0 < λ ′ 1 ≤ 1, the second eigenvalue λ ′ 2 becomes negative when the anomalous correlator gets larger than the normal correlator, |α| ≥ η. Then the probabilities of sampling an even total number of excited atoms is getting enhanced. For a particular choice of the observational basis consisting of standing plane waves (see section V), this property may be interpreted as a simple consequence of the fact that both independent occupation numbers n 1 and n 2 have strongly suppressed probabilities to be odd. For the observational basis consisting of traveling plane waves the joint statistics has a strong correlation at n 1 = n 2 (see section VI), and the total-occupation probability ρ n = ∑ n n 1 =0 ρ n 1 ,n−n 1 hits this diagonal only for even n. In the vacuum state, when the anomalous correlator |α| = η 2 + η achieves its maximal possible absolute value, we have λ ′ 1 = −λ ′ 2 = tanh r k , and the probabilities of odd total occupations completely vanish.

VIII. 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
In conclusion, we briefly overview the presented analysis of the ♯P-hard problem of atomic boson sampling from an interacting BEC gas trapped in a box. The analysis is done by means of a new approach based on the recently found hafhian master theorem 8,13 . We intentionally choose the simplest 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 has a potential for demonstrating quantum supremacy over classical computing and is one of the simplest models of quantum supremacy. 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 sections II-IV, while the simple explicit formulas and their numerical illustrations are given in sections V-VII.
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 eigensqueeze modes and the eigen-energy quasiparticles, which are directly responsible for the above-mentioned squeezing and interference. The eigen-squeeze modes are the eigenvectors of the multimode squeezing matrix. They are the elementary, intrinsic carriers of the squeezing matrix's eigenvalues -the single-mode squeezing parameters, Eq. (31). The quasiparticles are the eigenvectors of the Hamiltonian, Eq. (10), and are described by two-component wave functions, Eq. (11). They are the elementary collective excitations carrying quanta of collective energy -the eigenvalues of the Hamiltonian.
As such each quasiparticle lives completely independent on other quasiparticles, but it is a superposition of many twocomponent squeezed modes formed by the irreducible, eigensqueezing Bogoliubov transformationR r , Eq. (25), from the one-component eigen-squeeze modes, each of which is itself a superposition of many excited bare-atom wave functions. As a result, the excited bare-atom wave functions are permanently interfere 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 a case when each bare-atom excited state chosen for atom number measurements by a multi-detector imaging system coincides with the eigen-squeeze mode which simultaneously constitutes both components of the quasiparticle wave function. Then the joint probability distribution factorizes into the product of the occupation probabilities of single eigen-squeeze modes. Each such probability amounts to a nontrivial squeezed-state statistics, but is computable via Legendre polynomials as per Eq. (53).
Another example is a case when the squeezing is absent, say, due to the absence of inter-atomic 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 Eq. (26) vanish. Hence, the hafnian in Eq. (24) determining the joint sampling probabilities reduces to the permanent of a positive-definite matrix since the extended covariance-related matrix has vanishing anomalous-correlator blocks and one can employ the Stockmeyer's approximating algorithm 17,18 for such a permanent.
A pivotal key for understanding the origin of the ♯Phardness of atomic boson sampling is provided by the Bloch-Messiah reduction of the Bogoliubov transformation into the product of three irreducible blocks,R WRrRV , as per Eq. (25). Via its direct relation to the hafnian of the extended covariance-related matrix constituting the general result for the joint probability distribution in Eq. (24), 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 blockR r of the Bogoliubov transformation. The interference between the observational bare-atom excited states and the eigen-squeeze modes, controlled by the blockR V and its unitary V , constitutes the first mechanism of interference. The interference between the eigen-squeeze modes and the eigen-energy quasiparticles, controlled by the blockR W and its unitary W , constitutes the second mechanism of interference.
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 on average still exists. This is the case for the presented model of atomic boson sampling in a box trap with a uniform condensate when parameters of the quasiparticles and eigen-squeeze modes, including the unitary W defining their interference blockR W , are almost fixed by given trapping potential and couplings in Eq. (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 9 is useful, but not necessary for demonstrating quantum supremacy.
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 neither sophisticated external sources of bosons in squeezed states nor a controlled optical interferometer. 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.
Conceptually, the experiments are simple. In fact, the experiments on the statistics of the total occupation of excited atom states, that is the total noncondensate occupation, had been done already 32 . 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 multidetector imaging system. Then to reconfigure detectors and to measure sampling statistics for another, unitary-transformed subset of excited wave functions, and so on. There is no need neither in any controllable non-equilibrium unitary-evolution processes typical for most quantum-computing experiments nor in 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. It is described by the statistical operator that intrinsically involves properties which are ♯P-hard for computing.
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.