Next Article in Journal
Dynamics of a Quantum Common-Pool Resource Game with Homogeneous Players’ Expectations
Next Article in Special Issue
Daily Streamflow of Argentine Rivers Analysis Using Information Theory Quantifiers
Previous Article in Journal
Performance Optimization and Exergy Analysis of Thermoelectric Heat Recovery System for Gas Turbine Power Plants
Previous Article in Special Issue
Quasi-Hyperbolically Symmetric γ-Metric
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

by
Vitaly V. Kocharovsky
*,
Vladimir V. Kocharovsky
,
William D. Shannon
and
Sergey V. Tarasov
Department of Physics and Astronomy, Institute for Quantum Science and Engineering, Texas A&M University, College Station, TX 77843, USA
*
Author to whom correspondence should be addressed.
Entropy 2023, 25(12), 1584; https://doi.org/10.3390/e25121584
Submission received: 6 October 2023 / Accepted: 21 November 2023 / Published: 25 November 2023
(This article belongs to the Special Issue Selected Featured Papers from Entropy Editorial Board Members)

Abstract

:
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.

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 V = L 3 with periodic boundary conditions. Excited atoms are described by a field operator
ψ ^ ex ( r ) = l ϕ l ( r ) b ^ l ,
where b ^ l denotes an operator annihilating an atom in the bare-atom 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-Popov approximation [45,46] 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:
H ^ = 1 2 b ^ b ^ T H b ^ b ^ , H = Δ ˜ ϵ + Δ ϵ + Δ * Δ ˜ * , ϵ = ϕ l | ϵ ^ | ϕ l , ϵ ^ = 2 2 / ( 2 m ) , Δ = g N 0 ϕ l * | ϕ 0 | 2 ϕ l d 3 r , Δ ˜ = g N 0 ϕ l * ϕ l * ϕ 0 2 d 3 r .
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 vector b ^ on top of the column vector b ^ , and multiplied from the left by a 2-block row vector ( b ^ , b ^ ) = b ^ b ^ T . The bold-faced operator b ^ = { b ^ l } T or b ^ = { b ^ l } T denotes a column vector of annihilation or creation operators in the { ϕ l } basis; ϵ ^ = 2 2 / ( 2 m ) 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 g = 4 π 2 a s / m via the s-wave scattering length a s .
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 { ϕ l } the field operator and Hamiltonian acquire more specific and transparent form. In a usual basis of plane traveling waves
ϕ k = e i kr V , k = 2 π j L , j = { j x , j y , j z } N 3 ,
with the wave vectors k enumerated by the integer 3-vector j , we have
ψ ^ ex ( r ) = k 0 e i k r V a ^ k ,
H ^ = 1 2 a ^ a ^ T Δ ˜ ϵ + Δ ϵ + Δ * Δ ˜ * a ^ a ^ , ϵ k = 2 k 2 2 m , ϵ k k = ϵ k δ k k , Δ k , k = g N 0 δ k , k , Δ ˜ k , k = g N 0 δ k , k .
Hereinafter, we denote the creation and annihilation operators of bare atoms in the traveling-plane-wave states e i kr V by symbols a ^ = { a ^ k | k 0 } T and a ^ = { a ^ k | k 0 } T , respectively, as opposed to the above general-case operators b ^ and b ^ .
Equilibrium quantum many-body statistics of atoms in the BEC trap at a temperature T is determined by the statistical operator ρ ^ = e H ^ / T / tr e H ^ / 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 atom states { ϕ l } 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,
ρ { n l } = π π π π e i l τ l n l Θ { τ l } l d τ l 2 π , n l = 0 , 1 , ,
of the characteristic function of the atom-number operators:
Θ { τ l } = e i l τ l n ^ l tr e i l τ l n ^ l ρ ^ , n ^ l = b ^ l b ^ l .
The symbol tr { } 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:
Θ = 1 det ( 1 + G Z G ) ; Z = diag { z l } 0 0 diag { z l } .
Here, the symbol 0 stands for the matrix block with zero entries.
The characteristic function is determined by a covariance matrix
G = : b ^ b ^ ( b ^ , b ^ ) : = b ^ l b ^ l b ^ l b ^ l b ^ l b ^ l b ^ l b ^ l
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 ^ l b ^ l , and anomalous, b ^ l b ^ l , b ^ l b ^ l , correlators. Each variable τ l of the characteristic function Θ ( { τ l } ) appears in the matrix of variables Z in Equation (8) twice—via the entry z l = e i τ l in each of the two identical diagonal blocks diag { z l } ; the symbol 1 denotes the identity matrix. Note, that in the literature on Gaussian states [47], the covariance matrix is often defined with a half anti-commutator, ( b ^ l b ^ l + b ^ l b ^ l ) / 2 , replacing the normal product, b ^ l b ^ l , of the creation/annihilation operators, that adds a half identity matrix 1 / 2 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 τ l 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 { z l } within such a group.
Consider the quasiparticle creation and annihilation operators, which constitute the column vectors b ˜ ^ = { b ˜ ^ l } T and b ˜ ^ = { b ˜ ^ l } T and diagonalize the Hamiltonian in Equation (2),
H ^ = l E l b ˜ ^ l b ˜ ^ l .
The quasiparticle eigen energy is denoted as E l 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:
ψ ^ ex ( r ) = l u l ( r ) b ˜ ^ l + v l * ( r ) b ˜ ^ l .
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 ˜ ^ l b ˜ ^ l b ˜ ^ l b ˜ ^ 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 quasiparticle wave function
V | u l | 2 | v l | 2 d 3 r = 1 .
Suppose one chooses plane traveling waves e i kr V | k 0 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,
H ^ = k 0 E k a ˜ ^ k a ˜ ^ k , E k = ϵ k 2 + 2 g N 0 ϵ k ,
for the quasiparticle creation, annihilation operators a ˜ ^ = { a ˜ ^ k | k 0 } T , a ˜ ^ = { a ˜ ^ k | k 0 } T related to the corresponding bare-atom operators a ^ = { a ^ k | k 0 } T , a ^ = { a ^ k | k 0 } T as follows
a ^ k = u ¯ k a ˜ ^ k + v ¯ k * a ˜ ^ k , a ˜ ^ k = u ¯ k a ^ k v ¯ k * a ^ k .
Here, we introduced the amplitudes u ¯ k and v ¯ k * ,
u ¯ k = ϵ k + E k 2 ϵ k E k 1 2 ξ k + 1 ξ k , v ¯ k * = ϵ k E k 2 ϵ k E k 1 2 ξ k 1 ξ k ; ξ k = E k ϵ k ,
of the functions constituting the two-component wave function of those traveling-plane-wave quasiparticles
u k ( r ) = u ¯ k e i kr V , v k * ( r ) = v ¯ k * e i kr V .
In terms of such quasiparticles, the excited-atom field operator in Equation (1) acquires the following explicit form
ψ ^ ex ( r ) = k 0 u k ( r ) a ˜ ^ k + v k * ( r ) a ˜ ^ k .
Generally, the quasiparticles are completely independent of each other. Thus, the correlations between the quasiparticle creation/annihilation operators b ˜ ^ , b ˜ ^ , analogous to the bare-atom correlations in Equation (9), are given by the ( 2 × 2 ) -block diagonal matrix of the thermal occupations of quasiparticles
D = diag ( e E l / T 1 ) 1 0 0 diag ( e E l / T 1 ) 1 .
Employing the canonical Bose commutation relations b ^ l b ^ l b ^ l b ^ 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,
G = R D + 1 2 R 1 2 ,
via the Bogoliubov transformation which relates the uncorrelated quasiparticles to the squeezed and interfering bare-atom excitations. It is described by the ( 2 × 2 ) -block symplectic matrix [9] R ˜ and its inverse matrix R = R ˜ 1 as follows
b ˜ ^ b ˜ ^ = R ˜ b ^ b ^ , R ˜ = A B B T A T ; b ^ b ^ = R b ˜ ^ b ˜ ^ , R = A B * B A * , R ˜ = R 1 .
Below we assume that some finite number M of orthogonal excited states (i.e., atom wave functions) { ϕ l | l = 1 , , M } 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 ( 2 M × 2 M ) -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 Δ l l , Δ ˜ l l , and ϵ l l 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 k and k . Then the analysis gets easier. All bold-faced vectors (such as b ^ ) 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 ( 2 M × 2 M ) -matrices containing corresponding ( M × M ) -blocks such as Δ , Δ ˜ , ϵ , A , B , etc. For the sake of simplicity, we will denote any such matrix of a finite, reduced dimension 2 M × 2 M or M × M by the same symbol that stands for its infinite-dimensional counterpart.

3. The Hafnian Master Theorem and Sampling Probabilities

The Hafnian master theorem recently found in [9,48] 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 [49,50], it means that computing the hafnian and using it as an oracle is enough for polynomial-time 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 Equation (8),
Θ { z l } = { n l } haf C ˜ { n l } det ( 1 + G ) l z l n l n l ! .
It is expressed via the hafnian of an extended covariance-related ( 2 n × 2 n ) -matrix C ˜ { 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
C = P G ( 1 + G ) 1 , P = 0 1 1 0 .
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 l-th and ( M + l ) -th columns, respectively. If n l = 0 , then the corresponding rows and columns should be erased. The matrix P just permutes the diagonal and off-diagonal blocks of the ( 2 × 2 ) -block matrix G ( 1 + G ) 1 in a way that is appropriate for the hafnian.
The concept of the matrix hafnian had been introduced in the quantum field theory [51,52]. It expresses the Wick’s, or Isserlis’, theorem [53,54] stating that the mean value, denoted below by angles, of the product of an even number, 2 M , of the centered Gaussian random variables y 1 , y 2 , , y 2 M 1 , y 2 M is equal to the hafnian of their covariance 2 M × 2 M matrix G ¯ ,
haf G ¯ = y 1 y 2 y 2 M 1 y 2 M .
Employing the identity
ρ { n l } = l n l n l ! z l n l Θ | { z l = 0 } ,
following from the definition of the characteristic function, we obtain an explicit analytical formula for the joint probabilities of the atom numbers in Equation (6) via the hafnian as follows
ρ { n l } = haf C ˜ ( { n l } ) det ( 1 + G ) l n l ! .
A similar formula for the probabilities of output occupation numbers appears also in the theory of photonic boson sampling of Gaussian states in a linear interferometer [11,12].

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 C ˜ { n l } 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]
b ˜ ^ b ˜ ^ = R ˜ W R ˜ r R ˜ V b ^ b ^ ; R ˜ W = W * 0 0 W , R ˜ r = cosh Λ r sinh Λ r sinh Λ r cosh Λ r , R ˜ V = V * 0 0 V .
It gives the blocks of the Bogoliubov-transformation matrix in Equation (20) in the form of a singular value decomposition,
A = V T cosh Λ r W T , B = V sinh Λ r W T ,
and involves two unitary matrices, W and V, as well as the diagonal matrix of M single-mode squeezing parameters
Λ r = diag { r l | l = 1 , , M } , r l 0 .
The Bloch–Messiah reduction (26) can be written in the form
cosh r * ( sinh r * ) W * W ( sinh r ) W W T cosh r b ˜ ^ b ˜ ^ = ( W V ) * 0 0 W V b ^ b ^
that states that the effective evolution (resembling the Hamiltonian one, e i H ^ t ) of the quasiparticle and bare-atom creation/annihilation operators under the action of the unitary multimode squeeze, S ^ , and rotation, Φ ^ , operators, respectively, yields the same creation/annihilation operators:
S ^ b ˜ ^ S ^ = Φ ^ b ^ Φ ^ , S ^ b ˜ ^ S ^ = Φ ^ b ^ Φ ^ .
Such an interpretation, typical for quantum optics [55,56,57,58,59,60], employs the multimode squeeze and rotation operators,
S ^ = exp 1 2 b ˜ ^ T S b ˜ ^ b ˜ ^ T S b ˜ ^ ,
Φ ^ = exp i b ^ T Φ b ^ ,
determined by a symmetric matrix S = r e i θ = ( W Λ r W ) ( W W T ) = W Λ r W T (built of the Hermitian matrices r and θ ; the unitary matrix e i θ = W W T is symmetric) and a Hermitian matrix Φ = Φ generating the unitary e i Φ = W V , respectively. The Hermitian factor, r, of the multimode squeeze matrix, S, is a positive semi-definite Hermitian M × M matrix whose diagonal representation Λ r is determined by the unitary W as follows
r = W Λ r W , r w l = r l w l , W l l = ( w l ) l ,
so that cosh r = W cosh Λ r W , sinh r = W sinh Λ r W .
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 { w l | l = 1 , , M } and values { r l 0 | l = 1 , , M } 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, { r l } 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 ( 2 × 2 ) -block transformations R ˜ W R ˜ r R ˜ V 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,
ρ ^ = l e E l b ˜ ^ l b ˜ ^ l / T / tr e E l b ˜ ^ l b ˜ ^ l / T ,
expressed in terms of the eigen-energy quasiparticle creation and annihilation operators.
The three 2 × 2 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, R ˜ V , corresponds to the unitary rotation V * ,
φ l ( r ) = l V l l * ϕ l ( r ) , β ^ l = l V l l b ^ l , ψ ^ ex ( r ) = l φ l ( r ) β ^ l ,
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 eigen-squeeze modes { φ l } associated with the eigenvectors { w l } 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, { b ^ l } , to the annihilation operators of the eigen-squeeze modes, { β ^ l } , is performed by the unitary V and does not involve creation operators.
The second part of the Bogoliubov transformation, R ˜ r , is associated with the presence in the Hamiltonian the terms beyond the resonant-wave approximation (non-RWA terms), b ^ l b ^ l and b ^ l b ^ l , 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 φ l ( r ) constitutes the same spatial profile for both wave-function components u l ( r ) , v l * ( r ) of the l-th eigen-squeeze two-component excitation which acquires the corresponding nontrivial squeezing parameter r l 0 . Note, according to Equation (38), those two components have different amplitudes, + cosh r l and sinh r l .
The third part of the Bogoliubov transformation, R ˜ W , relates the creation, { b ^ l } , and annihilation, { b ^ l } , operators of the eigen-squeeze two-component excitations, formed on the second step, to the eigen-energy-quasiparticle creation, { b ˜ ^ l } , and annihilation, { b ˜ ^ l } , operators:
b ^ l = l W l l b ˜ ^ l , b ^ l = l W l l * b ˜ ^ l .
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 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 Equation (11) constitute the bases of the first and second components of the two-component functional space of the quasiparticle wave functions ( u l ( r ) , v l * ( r ) ,
u l ( r ) = + l W l l * ( cosh r l ) φ l ( r ) , v l * ( r ) = l W l l ( sinh r l ) φ l ( r ) .
Since the quasiparticle wave functions { u l } and { v l * } are fully predetermined, fixed by the coupling parameters (interactions) in the Hamiltonian (2) (of course, up to a possible degeneracy), the unitary matrix W = ( W ) 1 determines the unique (up to a possible degeneracy) eigen-squeeze modes { φ l ( r ) } which are the same for both components of the quasiparticle functional space as per equations inverse to Equations (37):
+ ( cosh r l ) φ l ( r ) = u l ( r ) = l W l l u l ( r ) , ( sinh r l ) φ l ( r ) = v l * ( r ) = l W l l * v l * ( r ) .
Each single-squeeze two-component excitation originating from the mode φ l ( r ) owns the single-mode squeezing parameter r l 0 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 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 eigen-squeeze mode φ l ( r ) just multiplied by 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 [9],
R Ω R T = Ω , Ω = 0 1 1 0 ,
of the Bogoliubov transformation, that is, canonical Bose commutation relations, 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 , 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 Δ ˜ l l , Δ l l , ϵ l l 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 R ˜ W is almost completely fixed. Yet, the first Bogoliubov transformation R ˜ V 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 { ϕ l | l = 1 , , M } 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 W , V 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 R ˜ r determined exclusively by the eigenvalues { r l } 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, { r l = 0 } , 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 P . 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 r implies annihilation of a superposition of quasiparticles with the amplitudes given by the first component of the quasiparticle wave functions { u l ( r ) | l = 1 , , M } and simultaneous creation of a superposition of quasiparticles with the amplitudes given by the second component of the quasiparticle wave functions { v l * ( r ) | l = 1 , , M } .
In the observational basis coinciding with the eigen-squeeze modes { φ l } , when R ˜ V = 1 , the covariance matrix (19) acquires a unique irreducible form G = G Q + G T 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 β ^ l , β ^ l of the eigen-squeeze modes:
G Q = R R 1 2 = sinh 2 Λ r sinh Λ r cosh Λ r sinh Λ r cosh Λ r sinh 2 Λ r ,
G T = R D R = cosh Λ r sinh Λ r sinh Λ r cosh Λ r q * 0 0 q cosh Λ r sinh Λ r sinh Λ r cosh Λ r .
Here, q = W diag { ( e E j / T 1 ) 1 | j = 1 , , M } W is the M × M 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 { ϕ l } matrices G and C, Equation (22), appears as the corresponding unitary transform of the irreducible one as per unitary rotation in Equation (35),
G = V T 0 0 V ( G Q + G T ) V * 0 0 V , C = V 0 0 V T P G Q + G T 1 + G Q + G T V * 0 0 V .
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 { ϕ l | l = 1 , , M } .

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 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 eigen-squeeze excited states would lead to the same atom-number sampling statistics.
Thus, the excited atoms are described by the field operator
ψ ^ ex ( r ) = 2 V k 0 sin ( k r ) s ^ k + cos ( kr ) c ^ k ,
written via operators s ^ k , c ^ k annihilating a bare atom in the corresponding sinusoidal states with a wave vector k = 2 π j L , j = { j x , j y , j z } N 3 . 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, Z + 3 = { j | j z > 0 ( j z = 0 & j y > 0 ) ( j z = j y = 0 & j x > 0 ) } , 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 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,
{ ϕ l } = 2 sin ( kr ) V , 2 cos ( kr ) V | k = 2 π j L , j Z + 3 ,
the Hamiltonian (2) splits into two independent parts involving separate operators related to sin ( kr ) or cos ( kr ) modes,
H ^ = k , k 0 [ ϵ kk + Δ kk s ^ k s ^ k + 1 2 Δ ˜ kk s ^ k s ^ k + s ^ k s ^ k + ϵ kk + Δ kk c ^ k c ^ k + 1 2 Δ ˜ kk c ^ k c ^ k + c ^ k c ^ k ] .
It is in contrast to the Hamiltonian in Equation (5) in which the modes e i kr and e i kr of the traveling-plane-wave basis e i kr V | k 0 were coupled to each other. The Hamiltonian in Equation (45), compared to the Hamiltonian in Equation (5), has the same bare-atom energies ϵ k and normal overlapping integrals Δ kk , but different, now diagonal anomalous overlapping integrals
ϵ k k = ϵ k δ k , k , Δ k k = Δ ˜ k k = g N 0 δ k , k .
In this case, the quasiparticle annihilation operators { c ˜ ^ k , s ˜ ^ k } are given by a similar Bogoliubov transformation via the corresponding bare-particle annihilation and creation operators:
s ^ k = u ¯ k s ˜ ^ k + v ¯ k s ˜ ^ k , c ^ k = u ¯ k c ˜ ^ k + v ¯ k c ˜ ^ k ; s ˜ ^ k = u ¯ k s ^ k v ¯ k s ^ k , c ˜ ^ k = u ¯ k c ^ k v ¯ k c ^ k .
The quasiparticle operators diagonalize the Hamiltonian (45),
H ^ = k 0 E k s ˜ ^ k s ˜ ^ k + c ˜ ^ k c ˜ ^ k , E k = ϵ k 2 + 2 g N 0 ϵ k ,
and turn the field operator (43), annihilating bare atoms, into a sum of both annihilation and creation quasiparticle operators:
ψ ^ ex ( r ) = k 0 sin ( kr ) ( u ¯ k s ˜ ^ k + v ¯ k s ˜ ^ k ) + cos ( kr ) ( u ¯ k c ˜ ^ k + v ¯ k c ˜ ^ k ) V / 2 .
In both exponential and sinusoidal, bases, the quasiparticle energies E k and factors u ¯ k , v ¯ k are the same, given in Equations (13) and (15). The energy E k is larger than the bare-atom energy ϵ k .
Thus, in the observational basis of the eigen-squeeze modes (44) 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 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 R ˜ W = R ˜ V = 1 , that is
s ˜ ^ k s ˜ ^ k = cosh r k sinh r k sinh r k cosh r k s ^ k s ^ k , c ˜ ^ k c ˜ ^ k = cosh r k sinh r k sinh r k cosh r k c ^ k c ^ k ; r k = 1 2 ln E k ϵ k .
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 given by Equations (40) and (41) as follows
G ( k ) = N ˜ k + 1 2 cosh 2 r k sinh 2 r k sinh 2 r k cosh 2 r k 1 2 = η α α η ; N ˜ k = 1 / e E k / T 1 , η s ^ k s ^ k = N ˜ k + 1 / 2 cosh 2 r k 1 / 2 , α s ^ k s ^ k = N ˜ k + 1 / 2 sinh 2 r k ; | α | η ( 1 + η ) .
Here, N ˜ 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 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 N ˜ k 0 . The maximum value of the anomalous correlator, max | α | = η ( 1 + η ) , is achieved when N ˜ k = 0 , that is, when the quasiparticles are in the vacuum state. The relations between the average quasiparticle occupation N ˜ k , the single-mode squeezing parameter r k 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, η = n , corresponds to a circle on the plane | α | , N ˜ k with the origin at the point ( 0 , 1 / 2 ) , i.e., ( N ˜ k + 1 / 2 ) 2 + | α | 2 = ( η + 1 / 2 ) 2 .
The eigenvalues of this covariance matrix depend on the single-mode squeezing parameter r k and the mean quasiparticle occupation N ˜ k ,
λ 1 , 2 = η ± | α | = N ˜ k e ± r k + e ± r k 1 / 2 .
The atom-number probability distribution in a sinusoidal mode is determined, as per Equation (22), by the covariance-related matrix C = P G ( k ) 1 + G ( k ) 1 which has the following explicit form:
C α η η α ; η = η 2 + η | α | 2 ( η + 1 ) 2 | α | 2 , α = α ( η + 1 ) 2 | α | 2 , det 1 + G ( k ) = ( η + 1 ) 2 | α | 2 = e 2 E k / T tanh 2 r k ( e E k / T 1 ) 2 cosh 2 r k .
It is convenient to denote the entries α , η of the matrix C = P G ( k ) 1 + 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 ) 1 + G ( k ) 1 look similar to the eigenvalues λ 1 , 2 of the covariance matrix G ( k ) , namely,
λ 1 , 2 = η ± | α | = η ± | α | 1 + η ± | α | = 1 ± e E k / T tanh r k e E k / T ± tanh r k .
Note, that the maximal value of the absolute value of the anomalous correlator, max | α | = η ( 1 + η ) , is achieved when η = 0 .
Now we calculate the atom-number probabilities. The result is immediately given by the Hafnian master theorem (25),
ρ n = haf C ˜ ( n ) n ! det ( 1 + G ( k ) ) , C ˜ ( n ) = α J n × n η J n × n η J n × n α J n × n ,
where the extended covariance-related matrix C ˜ ( 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 Equation (55) is reduced to the known permanent of the unity matrix J n × n , that is, haf C ˜ ( n ) = per η J n × n = n ! ( η ) n . In this case, the atom-number probability distribution is reduced to the simple exponential law
ρ n = 1 N ˜ k + 1 N ˜ k N ˜ k + 1 n , α = 0 ,
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
haf C ˜ ( n ) = ( 2 n 1 ) η haf C ˜ ( n 1 ) ( n 1 ) 2 ( η ) 2 | α | 2 haf C ˜ ( n 2 ) .
It starts from the plain values of the hafnian haf C ˜ ( 1 ) = η at n = 1 and haf C ˜ ( 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 [63] P n ,
n P n ( x ) = ( 2 n 1 ) x P n 1 ( x ) ( n 1 ) P n 1 ( x ) ,
we at once reduce the hafnian to the Legendre polynomials:
haf C ˜ ( n ) = n ! ( η ) 2 | α | 2 n / 2 P n η / ( η ) 2 | α | 2 .
Thus, the atom-number probabilities for sampling the occupation of the single-mode-squeezed standing plane wave are as follows
ρ n = ( η ) 2 | α | 2 n / 2 det ( 1 + G ( k ) ) P n η ( η ) 2 | α | 2 = 1 ( η + 1 ) 2 | α | 2 η 2 | α | 2 ( η + 1 ) 2 | α | 2 n / 2 P n η 2 + η | α | 2 ( η 2 | α | 2 ) ( η + 1 ) 2 | α | 2 .
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 Equation (55) can be calculated in a straightforward way via elementary combinatorics as follows
haf C ˜ ( n ) = ( η ) n k = 0 n / 2 ( n ! ) 2 4 k ( n 2 k ) ! ( k ! ) 2 α η 2 k .
Here, the symbol n / 2 stands for the largest integer less than or equal to n / 2 , and the essential argument is
| α | / η = e E k / T e E k / T sinh 2 r k .
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 η = 0 , or | α | = η ( 1 + η ) , the recursive relation becomes purely two-step,
haf C ˜ ( n ) = ( n 1 ) 2 | α | 2 haf C ˜ ( n 2 ) , η = 0 .
Equation (62) yields the result
ρ 2 m = 1 η + 1 ( 2 m 1 ) ! ! 2 ( 2 m ) ! η η + 1 m , ρ 2 m + 1 = 0 ,
which is equivalent (due to identities ( 2 m 1 ) ! ! 2 / ( 2 m ) ! = ( 2 m ) ! / 4 m / ( m ! ) 2 and max | α | = tanh r k ) to the one given by the reduction of the hafnian of a ( 2 × 2 ) -block-diagonal matrix to the square of the hafnian of just one of matrix’s blocks,
haf α J n × n 0 n × n 0 n × n α J n × n = haf α J n × n 2 = 0 , odd n n ! 2 n / 2 ( n / 2 ) ! ( α ) n 2 , even n ,
where the symbol 0 n 1 × n 2 denotes the n 1 × n 2 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],
ρ 2 m = 1 cosh r ( 2 m ) ! 4 m ( m ! ) 2 tanh 2 m r k , ρ 2 m + 1 = 0 .
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, N ˜ 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 Equation (59) for sampling from eigen-squeeze bare-atom excited states (44) are illustrated in Figure 3. In the range | α | α c r , where the absolute value of the anomalous correlator is less than the critical value α c r = η , an increasing absolute value 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 absolute value of the anomalous correlator exceeds the critical value, | α | > α c r = η , 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, max | α | = η ( 1 + η ) , which corresponds to the squeezed vacuum state, N ˜ k = 0 , all probabilities for even occupation numbers become zero.
The statistics of the eigen-squeeze-mode occupation can be deduced directly from its characteristic function Θ ( z ) . It 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,
ln Θ ( z ) = m = 1 κ ˜ m ( z 1 ) m m ! = m = 1 κ m ( ln z ) m m ! .
This characteristic function is given by the general result in Equation (8) as follows
Θ ( z ) = 1 det 1 + G ( k ) z G ( k ) = 1 det 1 + G ( k ) | 1 z tr ( P C ) + z 2 det ( P C ) .
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,
κ ˜ m = ( m 1 ) ! 2 tr G ( k ) m = ( m 1 ) ! 2 λ 1 m + λ 2 m ,
where λ 1 , 2 , Equation (52), are the eigenvalues of the covariance matrix (51). The value of the first generating cumulant κ ˜ m = 1 yields the mean occupation n s ^ k s ^ k = η whose thermal and quantum contributions are equal to ( e E k / T 1 ) 1 cosh 2 r k and sinh 2 r k , respectively. Adding the value of the second generating cumulant, we immediately obtain 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 the nontrivial anomalous correlator makes this characteristic function in Equation (67) very different from the one, Θ I B G ( z ) = [ 1 ( z 1 ) η ] 1 , associated with the atom-number sampling statistics in the limit of non-interacting, ideal Bogoliubov gas (IBG, α = 0 ) 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 { ϕ l } (which is the set of bare-atom excited states for atom-number measurements) to be the traveling plane waves e ± i kr , 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 a ^ k , a ^ k which differ from the annihilation operators s ^ k , c ^ k of the eigen-squeeze modes due to the unitary transformation
V = 1 2 + i i + 1 + 1 , s ^ k c ^ k = V a ^ + k a ^ k .
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 R ˜ r in Equation (50), the right-side unitary block R ˜ V 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, R ˜ W = 1 , 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 V * as per Equation (35),
2 sin ( k r ) cos ( k r ) = V * e + i k r e i k r = 1 2 i + i + 1 + 1 e + i k r e i k r .
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 ( a ^ k , a ^ k , a ^ k , a ^ k ) T . The related covariance matrix defined in Equation (9) is equal to
G exp ( k ) = η 0 0 α 0 η α 0 0 α η 0 α 0 0 η = N ˜ k + 1 2 × cosh 2 r k 0 0 sinh 2 r k 0 cosh 2 r k sinh 2 r k 0 0 sinh 2 r k cosh 2 r k 0 sinh 2 r k 0 0 cosh 2 r k 1 2 .
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 ( c ^ k , s ^ k , c ^ k , s ^ k ) T of the creation and annihilation operators of the eigen-squeeze modes in the form of the ( 4 × 4 ) -matrix
G sin ( k ) = η 0 α 0 0 η 0 α α 0 η 0 0 α 0 η = N ˜ k + 1 2 × cosh 2 r k 0 sinh 2 r k 0 0 cosh 2 r k 0 sinh 2 r k sinh 2 r k 0 cosh 2 r k 0 0 sinh 2 r k 0 cosh 2 r k 1 2 .
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,
C P G exp ( k ) 1 + G exp ( k ) 1 = 0 α η 0 α 0 0 η η 0 0 α 0 η α 0 ,
where the parameters η , α are the entries of the covariance-related matrix C in Equation (53) describing the eigen-squeeze mode 2 / V sin ( kr ) (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 e i kr / V and e i kr / V for a given wave vector k . The result is provided by the Hafnian master theorem (25),
ρ n 1 , n 2 = haf C ˜ ( n 1 , n 2 ) n 1 ! n 2 ! ( η + 1 ) 2 | α | 2 , C ˜ ( n 1 , n 2 ) = 0 n 1 × n 1 α J n 1 × n 2 η J n 1 × n 1 0 n 1 × n 2 α J n 2 × n 1 0 n 2 × n 2 0 n 2 × n 1 η J n 2 × n 2 η J n 1 × n 1 0 n 1 × n 2 0 n 1 × n 1 α J n 1 × n 2 0 n 2 × n 1 η J n 2 × n 2 α J n 2 × n 1 0 n 2 × n 2 .
Here, the hafnian matrix function is applied to the extended covariance-related matrix C ˜ ( 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 n 1 × n 2 is a zero matrix of the same size. The determinant in the denominator has been calculated explicitly as follows: det 1 + G exp ( k ) = ( η + 1 ) 2 | α | 2 2 .
The matrix C ˜ ( n 1 , n 2 ) consists of blocks of unequal dimensions. However, its total dimension 2 n 1 + 2 n 2 is even. Such a hafnian haf C ˜ ( n 1 , n 2 ) 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
haf 0 A A T 0 = per A
valid for any square matrix A. As a result, we calculate the required hafnian as follows
haf C ˜ ( 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 = k = 0 min ( n 1 , n 2 ) C n 1 k C n 2 k n 1 ! n 2 ! ( η ) n 1 + n 2 2 k | α | 2 k = k = 0 min ( n 1 , n 2 ) ( n 1 ! ) 2 ( n 2 ! ) 2 ( η ) n 1 + n 2 2 k | α | 2 k ( n 1 k ) ! ( n 2 k ) ! ( k ! ) 2 = n 1 ! n 2 ! ( η ) n 1 + n 2 2 F 1 n 1 , n 2 , 1 , ( α / η ) 2 .
Here, the permanent has been calculated explicitly by simple combinatorial means via binomial coefficients C n k = n ! k ! ( n k ) ! . Indeed, one just calculates the number of permutations of n 1 + n 2 elements that swap k elements between the block of the first n 1 elements and the block of the last n 2 elements, which corresponds to the summand ( η ) n 1 + n 2 2 k | α | 2 k . 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 [65]; however, the hypergeometrical representation is more convenient. The argument of the hypergeometric function in terms of the eigen-energies E k and eigen-squeezing parameters r k is given in Equation (61).
The hafnian in Equation (74) can be calculated also directly from its combinatorial definition [54]. 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 2 n 1 + 2 n 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 the 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 between the third and fourth groups. There are C n 1 k C n 2 k 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 of C ˜ . 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 . 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:
ρ n 1 , n 2 = n 1 ! n 2 ! ( η ) n 1 + n 2 ( η + 1 ) 2 | α | 2 k = 0 min ( n 1 , n 2 ) ( α / η ) 2 k ( n 1 k ) ! ( n 2 k ) ! ( k ! ) 2 = ( η ) n 1 + n 2 ( η + 1 ) 2 | α | 2 2 F 1 n 1 , n 2 , 1 , ( α / η ) 2 .
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, α = α = 0 , 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 C ˜ ( 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 = 1 η + 1 η η + 1 n 1 × 1 η + 1 η η + 1 n 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 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 n 1 and n 2 and growing the probability of the completely entangled occupation states n 1 = n 2 . This occurs in accord with a simultaneous increase of the pure quantum contribution G Q , Equation (40), associated with the quantum depletion of the condensate and growing up single-mode squeezing parameter r k , to the covariance matrix in Equation (42). The point is that the thermal quasiparticle occupation N ˜ k and, hence, the relative value of the complimentary thermal contribution G T , 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 k -modes have a close-to-maximum anomalous correlator | α | and small mean quasiparticle occupation N ˜ k . Achieving this in the BEC experiments requires a proper choice of the sampling k -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 max | α | = η ( 1 + η ) can be calculated directly from the hafnian formula in Equation (74) by means of setting η = 0 and employing the identity (75) as follows
ρ n 1 , n 2 1 n 1 ! n 2 ! [ ( η + 1 ) 2 | α | 2 ] per 0 n 1 × n 1 α J n 1 × n 2 α J n 2 × n 1 0 n 2 × n 2 = δ n 1 , n 2 ( α ) n 1 + n 2 ( η + 1 ) 2 | α | 2 .
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 ρ 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 (79) is reduced to a well-known result for the two-mode squeezed vacuum state [47,62]
ρ n 1 , n 2 = δ n 1 , n 2 cosh 2 r k tanh n 1 + n 2 r k .

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 sin ( kr ) , cos ( kr ) by means of an arbitrary unitary mixing via the unitary matrix
V = cos ξ e i β sin ξ e i β sin ξ cos ξ , s ^ k c ^ k = V b ^ k , 1 b ^ k , 2 .
It corresponds to the first part of the Bogoliubov transformation in Figure 1, from the annihilation operators b ^ k , 1 , b ^ k , 2 of these two excited atom states to the annihilation operators s ^ k , c ^ k 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, U = x y e i γ y * e i γ x * ,   | x | 2 + | y | 2 = 1 , γ R , 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 ξ [ 0 , π 4 ] , β [ 0 , π 2 ] . In view of Equation (35), wave functions of the selected measurement basis,
ϕ k , 1 ( r ) ϕ k , 2 ( r ) = V T 2 / V sin ( kr ) 2 / V cos ( kr ) ,
are the superpositions of standing and traveling plane waves.
The covariance matrix for the operators b ^ k , 1 , b ^ k , 2 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:
G = V T 0 0 V η 1 α 1 α 1 η 1 V * 0 0 V = η 1 α V T V α V V * η 1 .
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
C P G ( 1 + G ) 1 = α V V * η 1 η 1 α V T V .
Here, the amplitudes α and η (see Equation (53)) are the same as in Section 5 and Section 6.
Note, that if the 2 × 2 symmetric matrix V T V appeared in Equation (84) 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 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 V = V * and V T V = 1 .
The sampling probabilities are given by the Hafnian master theorem (25),
ρ n 1 , n 2 = haf C ˜ ( n 1 , n 2 ) n 1 ! n 2 ! ( η + 1 ) 2 | α | 2 , C ˜ ( n 1 , n 2 ) = α c 1 * J n 1 × n 1 α c 2 * J n 1 × n 2 η J n 1 × n 1 0 n 1 × n 2 α c 2 * J n 2 × n 1 α c 3 * J n 2 × n 2 0 n 2 × n 1 η J n 2 × n 2 η J n 1 × n 1 0 n 1 × n 2 α c 1 J n 1 × n 1 α c 2 J n 1 × n 2 0 n 2 × n 1 η J n 2 × n 2 α c 2 J n 2 × n 1 α c 3 J n 2 × n 2 .
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 V T V and V V * in Equation (84). Let us denote these entries as follows
V T V = c 1 c 2 c 2 c 3 , V V * = c 1 * c 2 * c 2 * c 3 * , c 1 = cos 2 ξ + e 2 i β sin 2 ξ , c 2 = 2 i cos ξ sin ξ sin β , c 3 = c 1 * .
We can calculate the hafnian directly from its combinatorial definition implementing its further reduction to a combinatorial problem of counting different partitions of 2 n 1 + 2 n 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 obtain the hafnian as the following sum
haf C ˜ ( n 1 , n 2 ) = ( n 1 ! n 2 ! ) 2 d 1 = 0 n 1 d 2 = 0 n 2 ( η ) d 1 + d 2 ( α ) n 1 + n 2 d 1 d 2 d 1 ! d 2 ! × j 1 , j 3 = 0 min ( n 1 d 1 , n 2 d 2 ) ( c 2 ) j 1 j 1 ! ( c 2 * ) j 3 j 3 ! c 1 2 n 1 d 1 j 1 2 ( n 1 d 1 j 1 2 ) ! c 3 2 n 2 d 2 j 1 2 ( n 2 d 2 j 1 2 ) ! c 1 * 2 n 1 d 1 j 3 2 ( n 1 d 1 j 3 2 ) ! c 3 * 2 n 2 d 2 j 3 2 ( n 2 d 2 j 3 2 ) ! .
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 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.
Figure 5 exemplifies how correlations between occupation numbers n 1 and n 2 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 ( n 1 , n 2 ) . Switching to purely traveling waves, which is the case considered in the previous section, ultimately leads to forming 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 + η ) , since then the pure quantum, condensate-depletion-based contribution G Q , Equation (40), to the covariance matrix (42) dominates the thermal one G T , Equation (41). Otherwise, for smaller | α | and larger η , the nontrivial correlation pattern is masked by thermal contributions due to large thermal population N ˜ k of quasiparticles and small single-mode squeezing parameter r k as per Figure 2.
All of the probability distribution patterns on the plane of stochastic variables n 1 and n 2 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, c 1 = c 3 * .
For a squeezed vacuum state, which implies zero mean quasiparticle occupation, the off-diagonal super-block of the extended covariance-related matrix C ˜ , Equation (85), is zero since η = 0 . 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,
haf C ˜ n 1 , n 2 = haf α c 1 J n 1 × n 1 α c 2 J n 1 × n 2 α c 2 J n 2 × n 1 α c 3 J n 2 × n 2 2 = ( n 1 ! n 2 ! ) 2 | α | n 1 + n 2 k = 0 μ / 2 c 2 μ 2 k ( c 1 c 3 / 4 ) k k ! ( μ 2 k ) ! k + | n 2 n 1 | / 2 ! 2 × c 1 2 2 max ( n 1 n 2 , 0 ) c 3 2 2 max ( n 2 n 1 , 0 ) , μ min ( n 1 , n 2 ) ,
for even n 1 + n 2 , and zero otherwise. Here, μ / 2 stands for the maximal integer less or equal μ / 2 , where μ is a minimal of two atom numbers n 1 and n 2 . 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 n 1 and n 2 . 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 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 η = 10 and the unitary mixing matrix V = 1 2 3 i i 3 . 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 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 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 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 G T , Equation (41), to the covariance matrix (42) due to increasing quasiparticle population N ˜ k , but also decrease and restructuring of the pure quantum contribution G Q , Equation (40), due to simultaneous decrease of the single-mode squeezing parameter r k 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 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 . 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,
ρ n n 1 = 0 n ρ n 1 , n n 1 ,
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, z j = z . (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 z 1 , 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
θ ( z , z ) = det 1 + η α α η det 1 z η α α η 1 .
Its generating cumulants are easy to calculate as follows
κ ˜ m = ( m 1 ) ! tr G m = ( m 1 ) ! ( η + | α | ) m + ( η | α | ) m = Γ ( m ) N ˜ k e r k + e r k 1 2 m + N ˜ k e r k + e r k 1 2 m ,
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 λ 1 , 2 of the renormalized covariance matrix G ( k ) 1 + G ( k ) 1 ,
θ ( z , z ) = 1 1 z λ 1 1 z λ 2 det 1 + η α α η = λ 1 1 z λ 1 λ 2 1 z λ 2 1 λ 1 λ 2 det 1 + η α α η .
Recalling Equation (54), λ 1 , 2 = η ± | α | , and taking into account the relation | α | det 1 + η α α η = | α | following from Equation (53), we finally obtain a simple formula for the probability of sampling n atoms total:
ρ n = ( λ 1 ) n + 1 ( λ 2 ) n + 1 2 | α | , λ 1 , 2 = 1 ± e E k / T tanh r k e E k / T ± tanh r k .
This result has been obtained for a pair of counter-propagating plane waves in [36,65]. Note, that the combinations of parameters Z ( ± A k ) or Y ± introduced in [36] or [65] 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, 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, 0 < λ 1 1 , the second eigenvalue λ 2 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 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 6), and the total-occupation probability ρ n = n 1 = 0 n ρ n 1 , n n 1 hits this diagonal only for even n. In the vacuum state, when the anomalous correlator achieves its maximal possible absolute value, | α | = η 2 + η , we have λ 1 = λ 2 = tanh r k , 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 R ˜ r , 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 k and k . 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, R ˜ W R ˜ r R ˜ V , 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 R ˜ r of the Bogoliubov transformation. The eigenvalues { r l } 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 R ˜ V 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 R ˜ W 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 R ˜ W , 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 C ˜ ( { n l } ) 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 e ± i kr , see Equation (35). Subsequent mixing of the first wave components of these pairs ( e + i kr waves) via a unitary V + and independent mixing of the second wave components of these pairs ( e i kr waves) via a different unitary V yield the bipartite protocol by involving a large-size, M 2 × M 2 matrix block V + ( tanh Λ r ) V 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 1 / M 1 / 2 (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, ∼ 1 / M 4 / 5 , 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.

Author Contributions

All coauthors, V.V.K. (Vitaly V. Kocharovsky), V.V.K. (Vladimir V. Kocharovsky), W.D.S. and S.V.T., contributed equally to deriving the results and writing the paper. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest

The authors have no conflict to disclose.

References

  1. Aaronson, S.; Arkhipov, A. The computational complexity of linear optics. In Proceedings of the Forty-Third Annual ACM Symposium on Theory of Computing (Association for Computing Machinery), New York, NY, USA, 6–8 June 2011; pp. 333–342. [Google Scholar] [CrossRef]
  2. Aaronson, S.; Arkhipov, A. The computational complexity of linear optics. Theory Comput. 2013, 9, 143–252. [Google Scholar] [CrossRef]
  3. Harrow, A.W.; Montanaro, A. Quantum computational supremacy. Nature 2017, 549, 203. [Google Scholar] [CrossRef] [PubMed]
  4. Brod, D.J.; Galvão, E.F.; Crespi, A.; Osellame, R.; Spagnolo, N.; Sciarrino, F. Photonic implementation of boson sampling: A review. Adv. Photonics 2019, 1, 034001. [Google Scholar] [CrossRef]
  5. Zhong, H.S.; Wang, H.; Deng, Y.H.; Chen, M.C.; Peng, L.C.; Luo, Y.H.; Qin, J.; Wu, D.; Ding, X.; Hu, Y.; et al. Quantum computational advantage using photons. Science 2020, 370, 1460–1463. [Google Scholar] [CrossRef]
  6. Boixo, S.; Isakov, S.V.; Smelyanskiy, V.N.; Babbush, R.; Ding, N.; Jiang, Z.; Bremner, M.J.; Martinis, J.M.; Neven, H. Characterizing quantum supremacy in near-term devices. Nat. Phys. 2018, 14, 595–600. [Google Scholar] [CrossRef]
  7. Arute, F.; Arya, K.; Babbush, R.; Bacon, D.; Bardin, J.C.; Barends, R.; Biswas, R.; Boixo, S.; Brandao, F.G.S.L.; Buell, D.A.; et al. Quantum supremacy using a programmable superconducting processor. Nature 2019, 574, 505–510. [Google Scholar] [CrossRef] [PubMed]
  8. Dalzell, A.M.; Harrow, A.W.; Koh, D.E.; Placa, R.L.L. How many qubits are needed for quantum computational supremacy? Quantum 2020, 4, 264. [Google Scholar] [CrossRef]
  9. Kocharovsky, V.V.; Kocharovsky, V.V.; Tarasov, S.V. Atomic boson sampling in a Bose–Einstein-condensed gas. Phys. Rev. A 2022, 106, 063312. [Google Scholar] [CrossRef]
  10. Kocharovsky, V.V.; Kocharovsky, V.V.; Shannon, W.D.; Tarasov, S.V. Multi-Qubit Bose–Einstein Condensate Trap for Atomic Boson Sampling. Entropy 2022, 24, 1771. [Google Scholar] [CrossRef]
  11. Hamilton, C.S.; Kruse, R.; Sansoni, L.; Barkhofen, S.; Silberhorn, C.; Jex, I. Gaussian Boson Sampling. Phys. Rev. Lett. 2017, 119, 170501. [Google Scholar] [CrossRef]
  12. Kruse, R.; Hamilton, C.S.; Sansoni, L.; Barkhofen, S.; Silberhorn, C.; Jex, I. Detailed study of Gaussian boson sampling. Phys. Rev. A 2019, 100, 032326. [Google Scholar] [CrossRef]
  13. Grier, D.; Brod, D.J.; Arrazola, J.M.; Alonso, M.B.d.; Quesada, N. The complexity of bipartite Gaussian boson sampling. Quantum 2022, 6, 863. [Google Scholar] [CrossRef]
  14. Lund, A.P.; Laing, A.; Rahimi-Keshari, S.; Rudolph, T.; O’Brien, J.L.; Ralph, T.C. Boson Sampling from a Gaussian State. Phys. Rev. Lett. 2014, 113, 100502. [Google Scholar] [CrossRef] [PubMed]
  15. Shi, J.; Byrnes, T. Effect of partial distinguishability on quantum supremacy in Gaussian Boson sampling. NPJ Quantum Inf. 2022, 8, 54. [Google Scholar] [CrossRef]
  16. Chin, S.; Huh, J. Generalized concurrence in boson sampling. Sci. Rep. 2018, 8, 6101. [Google Scholar] [CrossRef]
  17. Quesada, N.; Arrazola, J.M.; Killoran, N. Gaussian boson sampling using threshold detectors. Phys. Rev. A 2018, 98, 062322. [Google Scholar] [CrossRef]
  18. Yung, M.-H.; Gao, X.; Huh, J. Universal bound on sampling bosons in linear optics and its computational implications. Natl. Sci. Rev. 2019, 6, 719–729. [Google Scholar] [CrossRef]
  19. Kim, Y.; Hong, K.-H.; Kim, Y.-H.; Huh, J. Connection between BosonSampling with quantum and classical input states. Opt. Express 2020, 28, 6929–6936. [Google Scholar] [CrossRef]
  20. Villalonga, B.; Niu, M.Y.; Li, L.; Neven, H.; Platt, J.C.; Smelyanskiy, V.N.; Boixo, S. Efficient approximation of experimental Gaussian boson sampling. arXiv 2021, arXiv:2109.11525v1. [Google Scholar]
  21. Bentivegna, M.; Spagnolo, N.; Vitelli, C.; Brod, D.J.; Crespi, A.; Flamini, F.; Ramponi, R.; Mataloni, P.; Osellame, R.; Galvao, E.F.; et al. Bayesian approach to boson sampling validation. Int. J. Quantum. Inform. 2015, 12, 1560028. [Google Scholar] [CrossRef]
  22. Renema, J.J.; Menssen, A.; Clements, W.R.; Triginer, G.; Kolthammer, W.S.; Walmsley, I.A. Efficient Classical Algorithm for Boson Sampling with Partially Distinguishable Photons. Phys. Rev. Lett. 2018, 120, 220502. [Google Scholar] [CrossRef] [PubMed]
  23. Renema, J.J. Simulability of partially distinguishable superposition and Gaussian boson sampling. Phys. Rev. A 2020, 101, 063840. [Google Scholar] [CrossRef]
  24. Popova, A.; Rubtsov, A. Cracking the Quantum Advantage threshold for Gaussian Boson Sampling. arXiv 2021, arXiv:2106.01445. [Google Scholar]
  25. Qi, H.; Brod, D.J.; Quesada, N.; García-Patrón, R. Regimes of Classical Simulability for Noisy Gaussian Boson Sampling. Phys. Rev. Lett. 2020, 124, 100502. [Google Scholar] [CrossRef] [PubMed]
  26. Rahimi-Keshari, S.; Lund, A.P.; Ralph, T.C. What can quantum optics say about computational complexity theory? Phys. Rev. Lett. 2015, 114, 060501. [Google Scholar] [CrossRef]
  27. Lim, Y.; Oh, C. Approximating outcome probabilities of linear optical circuits. arXiv 2022, arXiv:2211.07184v1. [Google Scholar]
  28. Shchesnovich, V.S. Noise in boson sampling and the threshold of efficient classical simulatability. Phys. Rev. A 2019, 100, 012340. [Google Scholar] [CrossRef]
  29. Wang, H.; Qin, J.; Ding, X.; Chen, M.C.; Chen, S.; You, X.; He, T.-M.; Jiang, X.; You, L.; Wang, Z.; et al. Boson Sampling with 20 input photons and a 60-mode interferometer in a 1014-dimensional Hilbert space. Phys. Rev. Lett. 2019, 123, 250503. [Google Scholar] [CrossRef]
  30. Zhong, H.S.; Deng, Y.H.; Qin, J.; Wang, H.; Chen, M.C.; Peng, L.C.; Luo, Y.H.; Wu, D.; Gong, S.Q.; Su, H.; et al. Phase-Programmable Gaussian Boson Sampling Using Stimulated Squeezed Light. Phys. Rev. Lett. 2021, 127, 180502. [Google Scholar] [CrossRef]
  31. Madsen, L.S.; Laudenbach, F.; Askarani, M.F.; Rortais, F.; Vincent, T.; Bulmer, J.F.; Miatto, F.M.; Neuhaus, L.; Helt, L.G.; Collins, M.J.; et al. Quantum computational advantage with a programmable photonic processor. Nature 2022, 606, 75–81. [Google Scholar] [CrossRef]
  32. Bentivegna, M.; Spagnolo, N.; Vitelli, C.; Flamini, F.; Viggianiello, N.; Latmiral, L.; Mataloni, P.; Brod, D.J.; Galvao, E.F.; Crespi, A.; et al. Experimental scattershot b Experimental scattershot boson sampling. Sci. Adv. 2015, 1, e1400255. [Google Scholar] [CrossRef]
  33. Zhong, H.S.; Peng, L.C.; Li, Y.; Hu, Y.; Li, W.; Qin, J.; Wu, D.; Zhang, W.; Li, H.; Zhang, L.; et al. Experimental Gaussian Boson sampling. Sci. Bull. 2019, 64, 511–515. [Google Scholar] [CrossRef] [PubMed]
  34. Wang, H.; He, Y.; Li, Y.H.; Su, Z.E.; Li, B.; Huang, H.L.; Ding, X.; Chen, M.C.; Liu, C.; Qin, J.; et al. High-efficiency multiphoton boson sampling. Nat. Photonics 2017, 11, 361–365. [Google Scholar] [CrossRef]
  35. Loredo, J.C.; Broome, M.A.; Hilaire, P.; Gazzano, O.; Sagnes, I.; Lemaitre, A.; Almeida, M.P.; Senellart, P.; White, A.G. Boson Sampling with Single-Photon Fock States from a Bright Solid-State Source. Phys. Rev. Lett. 2017, 118, 130503. [Google Scholar] [CrossRef]
  36. Kocharovsky, V.V.; Kocharovsky, V.V.; Scully, M.O. Condensation of N bosons. III. Analytical results for all higher moments of condensate fluctuations in interacting and ideal dilute Bose gases via the canonical ensemble quasiparticle formulation. Phys. Rev. A 2000, 61, 053606. [Google Scholar] [CrossRef]
  37. Tarasov, S.V.; Kocharovsky, V.V.; Kocharovsky, V.V. Bose–Einstein condensate fluctuations versus an interparticle interaction. Phys. Rev. A 2020, 102, 043315. [Google Scholar] [CrossRef]
  38. Sinatra, A.; Castin, Y.; Li, Y. Particle number fluctuations in a cloven trapped Bose gas at finite temperature. Phys. Rev. A 2010, 81, 053623. [Google Scholar] [CrossRef]
  39. Klawunn, M.; Recati, A.; Pitaevskii, L.P.; Stringari, S. Local atom-number fluctuations in quantum gases at finite temperature. Phys. Rev. A 2011, 84, 033612. [Google Scholar] [CrossRef]
  40. Calzetta, E.A.; Hu, B.L. Bose–Einstein condensate collapse and dynamical squeezing of vacuum fluctuations. Phys. Rev. A 2003, 68, 043625. [Google Scholar] [CrossRef]
  41. Opanchuk, B.; Rosales-Zárate, L.; Teh, R.Y.; Dalton, B.J.; Sidorov, A.; Drummond, P.D.; Reid, M.D. Mesoscopic two-mode entangled and steerable states of 40 000 atoms in a Bose–Einstein-condensate interferometer. Phys. Rev. A 2019, 100, 060102(R). [Google Scholar] [CrossRef]
  42. Shin, Y.; Saba, M.; Pasquini, T.A.; Ketterle, W.; Pritchard, D.E.; Leanhardt, A.E. Atom Interferometry with Bose–Einstein Condensation in a Double-Well Potential. Phys. Rev. Lett. 2004, 92, 050405-1. [Google Scholar] [CrossRef] [PubMed]
  43. Egorov, M.; Anderson, R.P.; Ivannikov, V.; Opanchuk, B.; Drummond, P.; Hall, B.V.; Sidorov, A.I. Long-lived periodic revivals of coherence in an interacting Bose–Einstein condensate. Phys. Rev. A 2011, 84, 021605(R). [Google Scholar] [CrossRef]
  44. Berrada, T.; van Frank, S.; Bucker, R.; Schumm, T.; Schaff, J.-F.; Schmiedmayer, J. Integrated Mach–Zehnder interferometer for Bose–Einstein condensates. Nat. Commun. 2013, 4, 2077. [Google Scholar] [CrossRef] [PubMed]
  45. Shi, H.; Griffin, A. Finite-temperature excitations in a dilute Bose-condensed gas. Phys. Rep. 1998, 304, 1–87. [Google Scholar] [CrossRef]
  46. Zagrebnov, V.A.; Bru, J.B. The Bogoliubov model of weakly imperfect Bose gas. Phys. Rep. 2001, 350, 291–434. [Google Scholar] [CrossRef]
  47. Weedbrook, C.; Pirandola, S.; García-Patrón, R.; Cerf, N.J.; Ralph, T.C.; Shapiro, J.H.; Lloyd, S. Gaussian quantum information. Rev. Mod. Phys. 2012, 84, 621–669. [Google Scholar] [CrossRef]
  48. Kocharovsky, V.V.; Kocharovsky, V.V.; Tarasov, S.V. The Hafnian Master Theorem. Linear Algebra Appl. 2022, 651, 144–161. [Google Scholar] [CrossRef]
  49. Toda, S. PP is as hard as the polynomial-time hierarchy. SIAM J. Comput. 1991, 20, 865–877. [Google Scholar] [CrossRef]
  50. Basu, S. A complex analog of Toda’s theorem. Found. Comput. Math. 2012, 12, 327–362. [Google Scholar] [CrossRef]
  51. Caianiello, E.R. On quantum field theory-I: Explicit solution of Dyson’s equation in electrodynamics without use of Feynman graphs. Nuovo Cim. 1953, 10, 1634–1652. [Google Scholar] [CrossRef]
  52. Caianiello, E.R. Combinatorics and Renormalization in Quantum Field Theory. In Frontiers in Physics; W. A. Benjamin Inc.: London, UK, 1973. [Google Scholar]
  53. Wick, G.C. The Evaluation of the collision matrix. Phys. Rev. 1950, 80, 268–272. [Google Scholar] [CrossRef]
  54. Barvinok, A. Combinatorics and Complexity of Partition Functions, Algorithms and Combinatorics 30; Springer International Publishing AG: Cham, Switzerland, 2016. [Google Scholar]
  55. Braunstein, S.L. Squeezing as an irreducible resource. Phys. Rev. A 2005, 71, 055801. [Google Scholar] [CrossRef]
  56. Cariolaro, G.; Pierobon, G. Reexamination of Bloch–Messiah reduction. Phys. Rev. A 2016, 93, 062115. [Google Scholar] [CrossRef]
  57. Vogel, W.; Welsch, D.-G. Quantum Optics, 3rd ed.; Wiley-VCH Verlag GmbH: Berlin, Germany, 2006. [Google Scholar]
  58. Huh, J.; Yung, M.-H. Vibronic Boson Sampling: Generalized Gaussian Boson Sampling for Molecular Vibronic Spectra at Finite Temperature. Sci. Rep. 2017, 7, 7462. [Google Scholar] [CrossRef] [PubMed]
  59. Huh, J. Multimode Bogoliubov transformation and Husimi’s Q-function. J. Phys. Conf. Ser. 2020, 1612, 012015. [Google Scholar] [CrossRef]
  60. Ma, X.; Rhodes, W. Multimode squeeze operators and squeezed states. Phys. Rev. A 1990, 41, 4625–4631. [Google Scholar] [CrossRef]
  61. Stockmeyer, L. On approximation algorithms for ♯P. SIAM J. Comput. 1985, 14, 849–861. [Google Scholar] [CrossRef]
  62. Barnett, S.M.; Radmore, P. Methods in Theoretical Quantum Optics; Oxford University Press: Oxford, UK, 1996. [Google Scholar]
  63. Abramowitz, M.; Stegun, I.A. Handbook of Mathematical Functions; National Bureau of Standards (NBS): Washington, DC, USA, 1964.
  64. Gould, H.W. Combinatorial Identities; Morgantown Printing and Binding Co.: Morgantown, WV, USA, 1972. [Google Scholar]
  65. Englert, B.-G.; Fulling, S.A.; Pilloff, M.D. Statistics of dressed modes in a thermal state. Opt. Commun. 2002, 208, 139–144. [Google Scholar] [CrossRef]
  66. Schrade, G.; Akulin, V.M.; Man’ko, V.I.; Schleich, W.P. Photon statistics of a two-mode squeezed vacuum. Phys. Rev. A 1993, 48, 2398. [Google Scholar] [CrossRef]
  67. Kaufman, A.M.; Tichy, M.C.; Mintert, F.; Rey, A.M.; Regal, C.A. The Hong–Ou–Mandel effect with atoms. Adv. At. Mol. Opt. Phys. 2018, 67, 377–427. [Google Scholar] [CrossRef]
  68. Lopes, R.; Imanaliev, A.; Aspect, A.; Cheneau, M.; Boiron, D.; Westbrook, C.I. Atomic Hong–Ou–Mandel experiment. Nature 2015, 520, 66–68. [Google Scholar] [CrossRef]
  69. Valiant, L.G. The complexity of computing the permanent. Theor. Comput. Sci. 1979, 8, 189–201. [Google Scholar] [CrossRef]
  70. Jerrum, M.; Sinclair, A.; Vigoda, E. A polynomial-time approximation algorithm for the permanent of a matrix with nonnegative entries. J. ACM 2004, 51, 671–697. [Google Scholar] [CrossRef]
  71. Bjorklund, A.; Gupt, B.; Quesada, N. A faster hafnian formula for complex matrices and its benchmarking on a supercomputer. ACM J. Exp. Algorithmics 2019, 24, 1–17. [Google Scholar] [CrossRef]
  72. Kocharovsky, V.V.; Kocharovsky, V.V.; Tarasov, S.V. Unification of the nature’s complexities via a matrix permanent—Critical phenomena, fractals, quantum computing, ♯P-complexity. Entropy 2020, 22, 322. [Google Scholar] [CrossRef] [PubMed]
  73. Kristensen, M.A.; Christensen, M.B.; Gajdacz, M.; Iglicki, M.; Pawłowski, K.; Klempt, C.; Sherson, J.F.; Rzążewski, K.; Hilliard, A.J.; Arlt, J.J. Observation of atom number fluctuations in a Bose–Einstein condensate. Phys. Rev. Lett. 2019, 122, 163601. [Google Scholar] [CrossRef]
  74. Jiang, T. How many entries of a typical orthogonal matrix can be approximated by independent normals? Ann. Probab. 2006, 34, 1497–1529. [Google Scholar] [CrossRef]
  75. Jiang, T. The entries of circular orthogonal ensembles. J. Math. Phys. 2009, 50, 063302. [Google Scholar] [CrossRef]
  76. Tenart, A.; Hercé, G.; Bureik, J.-P.; Dareau, A.; Clément, D. Observation of pairs of atoms at opposite momenta in an equilibrium interacting Bose gas. Nat. Phys. 2021, 17, 1364–1368. [Google Scholar] [CrossRef]
  77. Hercé, G.; Bureik, J.-P.; Ténart, A.; Aspect, A.; Dareau, A.; Clément, D. Full counting statistics of interacting lattice gases after an expansion: The role of condensate depletion in many-body coherence. Phys. Rev. Res. 2023, 5, L012037. [Google Scholar] [CrossRef]
  78. Armijo, J.; Jacqmin, T.; Kheruntsyan, K.V.; Bouchoule, I. Probing three-body correlations in a quantum gas using the measurement of the third moment of density fluctuations. Phys. Rev. Lett. 2010, 105, 230402. [Google Scholar] [CrossRef]
  79. Jacqmin, T.; Armijo, J.; Berrada, T.; Kheruntsyan, K.V.; Bouchoule, I. Sub-Poissonian fluctuations in a 1D Bose gas: From the quantum quasicondensate to the strongly interacting regime. Phys. Rev. Lett. 2010, 105, 230405. [Google Scholar] [CrossRef]
  80. Esteve, J.; Trebbia, J.-B.; Schumm, T.; Aspect, A.; Westbrook, C.I.; Bouchoule, I. Observations of density fluctuations in an elongated Bose gas: Ideal gas and quasicondensate regimes. Phys. Rev. Lett. 2006, 96, 130403. [Google Scholar] [CrossRef] [PubMed]
  81. Chuu, C.-S.; Schreck, F.; Meyrath, T.P.; Hanssen, J.L.; Price, G.N.; Raizen, M.G. Direct observation of sub-Poissonian number statistics in a degenerate Bose gas. Phys. Rev. Lett. 2005, 95, 260403. [Google Scholar] [CrossRef] [PubMed]
  82. Dotsenko, I.; Alt, W.; Khudaverdyan, M.; Kuhr, S.; Meschede, D.; Miroshnychenko, Y.; Schrader, D.; Rauschenbeutel, A. Submicrometer Position Control of Single Trapped Neutral Atoms. Phys. Rev. Lett. 2005, 95, 033002. [Google Scholar] [CrossRef] [PubMed]
  83. Schlosser, N.; Reymond, G.; Grangier, P. Collisional Blockade in Microscopic Optical Dipole Traps. Phys. Rev. Lett. 2002, 89, 023005. [Google Scholar] [CrossRef] [PubMed]
  84. Pons, M.; del Campo, A.; Muga, J.G.; Raizen, M.G. Preparation of atomic Fock states by trap reduction. Phys. Rev. A 2009, 79, 033629. [Google Scholar] [CrossRef]
Figure 1. Schematic diagram of the Bloch–Messiah reduction in Equation (26): Three irreducible steps of the Bogoliubov transformation (26) of the creation/annihilation operators and wave functions from the observational bare-atom basis to the quasiparticle basis.
Figure 1. Schematic diagram of the Bloch–Messiah reduction in Equation (26): Three irreducible steps of the Bogoliubov transformation (26) of the creation/annihilation operators and wave functions from the observational bare-atom basis to the quasiparticle basis.
Entropy 25 01584 g001
Figure 2. An average quasiparticle occupation, N ˜ k , (red solid curves) and the single-mode squeezing parameter, r k , (blue dashed curves) vs the absolute value of the anomalous correlator | s ^ k s ^ k | | α | for different fixed values of the normal correlator s ^ k s ^ 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 ) .
Figure 2. An average quasiparticle occupation, N ˜ k , (red solid curves) and the single-mode squeezing parameter, r k , (blue dashed curves) vs the absolute value of the anomalous correlator | s ^ k s ^ k | | α | for different fixed values of the normal correlator s ^ k s ^ 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 ) .
Entropy 25 01584 g002
Figure 3. Probability ρ n for sampling n atoms in a single eigen-squeeze mode (44): Orange points correspond to the moderate values of the anomalous correlator | α | = 9.5 (left panel) or | α | = 10.4 (central panel); blue and green points on the left and central panels show probabilities ρ n in the case of zero ( α = 0 ) and maximal ( | α | = η ( 1 + η ) 10.5 ) squeezing, respectively. The right panel exemplifies a dependence of the probabilities ρ 0 , , ρ 5 on the anomalous correlator | α | . The normal correlator for all panels has the same value η = 10 .
Figure 3. Probability ρ n for sampling n atoms in a single eigen-squeeze mode (44): Orange points correspond to the moderate values of the anomalous correlator | α | = 9.5 (left panel) or | α | = 10.4 (central panel); blue and green points on the left and central panels show probabilities ρ n in the case of zero ( α = 0 ) and maximal ( | α | = η ( 1 + η ) 10.5 ) squeezing, respectively. The right panel exemplifies a dependence of the probabilities ρ 0 , , ρ 5 on the anomalous correlator | α | . The normal correlator for all panels has the same value η = 10 .
Entropy 25 01584 g003
Figure 4. 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 η (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.7 , (c) | α | = 10.4 , and (d) max | α | = η ( 1 + η ) (which corresponds to the two-mode squeezed vacuum state). Increasing | α | leads to growing up correlations between the random numbers n 1 and n 2 . Note, that the scale of the values ρ n 1 , n 2 changes significantly from (ad).
Figure 4. 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 η (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.7 , (c) | α | = 10.4 , and (d) max | α | = η ( 1 + η ) (which corresponds to the two-mode squeezed vacuum state). Increasing | α | leads to growing up correlations between the random numbers n 1 and n 2 . Note, that the scale of the values ρ n 1 , n 2 changes significantly from (ad).
Entropy 25 01584 g004
Figure 5. The joint probability distribution ρ n 1 , n 2 for atomic boson sampling from two excited atom states formed by a general-case unitary mixing of two eigen-squeeze modes (44) with the same wave vector k . The normal and anomalous correlators in Equation (51) are η = 10 and | α | = 10.4 , respectively, which amounts to the mean quasiparticle occupation N ˜ k 0.95 and the single-mode squeezing parameter r k 1.34 . The panels show the evolution of the joint probability distribution by adjusting the matrix V of the unitary mixing in Equation (26): (a) V = 1 0 0 1 (the sampling states coincide with the eigen-squeeze modes and the probability distribution factorizes into the product of two single mode distributions given in Equation (59)), (b) V = 1 2 2 3 + 1 ( 1 3 ) i ( 1 3 ) i 3 + 1 , (c) V = 1 2 3 i i 3 , (d) V = 1 2 + 1 i i + 1 (the sampling states coincide with the two counter-propagating plane waves and the probability distribution is given by the hypergeometric function in Equation (77)). The scaling on all panels is the same because the most probable outcome n 1 = n 2 = 0 represents a Fock state invariant under unitary transformations.
Figure 5. The joint probability distribution ρ n 1 , n 2 for atomic boson sampling from two excited atom states formed by a general-case unitary mixing of two eigen-squeeze modes (44) with the same wave vector k . The normal and anomalous correlators in Equation (51) are η = 10 and | α | = 10.4 , respectively, which amounts to the mean quasiparticle occupation N ˜ k 0.95 and the single-mode squeezing parameter r k 1.34 . The panels show the evolution of the joint probability distribution by adjusting the matrix V of the unitary mixing in Equation (26): (a) V = 1 0 0 1 (the sampling states coincide with the eigen-squeeze modes and the probability distribution factorizes into the product of two single mode distributions given in Equation (59)), (b) V = 1 2 2 3 + 1 ( 1 3 ) i ( 1 3 ) i 3 + 1 , (c) V = 1 2 3 i i 3 , (d) V = 1 2 + 1 i i + 1 (the sampling states coincide with the two counter-propagating plane waves and the probability distribution is given by the hypergeometric function in Equation (77)). The scaling on all panels is the same because the most probable outcome n 1 = n 2 = 0 represents a Fock state invariant under unitary transformations.
Entropy 25 01584 g005
Figure 6. The joint probability distribution ρ n 1 , n 2 for atomic boson sampling from two excited atom states formed by a unitary mixing of two eigen-squeeze modes (44) with the same wave vector k . The unitary matrix is V = 1 2 3 i i 3 , which coincides with the matrix corresponding to the panel (c) in Figure 5. The system of excited atoms is in a squeezed vacuum state with zero quasiparticle occupations, N ˜ k = 0 . The single-mode squeezing parameter, normal and anomalous correlators in Equation (51) are r k 1.87 , η = 10 and | α | = max | α | = 110 , respectively.
Figure 6. The joint probability distribution ρ n 1 , n 2 for atomic boson sampling from two excited atom states formed by a unitary mixing of two eigen-squeeze modes (44) with the same wave vector k . The unitary matrix is V = 1 2 3 i i 3 , which coincides with the matrix corresponding to the panel (c) in Figure 5. The system of excited atoms is in a squeezed vacuum state with zero quasiparticle occupations, N ˜ k = 0 . The single-mode squeezing parameter, normal and anomalous correlators in Equation (51) are r k 1.87 , η = 10 and | α | = max | α | = 110 , respectively.
Entropy 25 01584 g006
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Kocharovsky, V.V.; Kocharovsky, V.V.; Shannon, W.D.; Tarasov, S.V. Towards the Simplest Model of Quantum Supremacy: Atomic Boson Sampling in a Box Trap. Entropy 2023, 25, 1584. https://doi.org/10.3390/e25121584

AMA Style

Kocharovsky VV, Kocharovsky VV, Shannon WD, Tarasov SV. Towards the Simplest Model of Quantum Supremacy: Atomic Boson Sampling in a Box Trap. Entropy. 2023; 25(12):1584. https://doi.org/10.3390/e25121584

Chicago/Turabian Style

Kocharovsky, Vitaly V., Vladimir V. Kocharovsky, William D. Shannon, and Sergey V. Tarasov. 2023. "Towards the Simplest Model of Quantum Supremacy: Atomic Boson Sampling in a Box Trap" Entropy 25, no. 12: 1584. https://doi.org/10.3390/e25121584

APA Style

Kocharovsky, V. V., Kocharovsky, V. V., Shannon, W. D., & Tarasov, S. V. (2023). Towards the Simplest Model of Quantum Supremacy: Atomic Boson Sampling in a Box Trap. Entropy, 25(12), 1584. https://doi.org/10.3390/e25121584

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop