Multi-Qubit Bose–Einstein Condensate Trap for Atomic Boson Sampling

Abstract

We propose a multi-qubit Bose–Einstein-condensate (BEC) trap as a platform for studies of quantum statistical phenomena in many-body interacting systems. In particular, it could facilitate testing atomic boson sampling of the excited-state occupations and its quantum advantage over classical computing in a full, controllable and clear way. Contrary to a linear interferometer enabling Gaussian boson sampling of non-interacting non-equilibrium photons, the BEC trap platform pertains to an interacting equilibrium many-body system of atoms. We discuss a basic model and the main features of such a multi-qubit BEC trap.

1. Introduction to Quantum Statistical Physics of Atomic Boson Sampling in a BEC Trap

1.1. The Essence of the Problem

Recently, a stationary stochastic process of many-body fluctuations of the excited-atom occupations in a trapped Bose–Einstein-condensed gas has been suggested for quantum simulation of the ♯P-hard problem of boson sampling [1]. Such an atomic boson sampling, based on the Bose–Einstein-condensate (BEC) platform, is an alternative to a well-known photonic boson sampling based on the linear interferometer platform [2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25]. It has the potential to demonstrate quantum advantage [26,27,28,29,30] of the many-body interacting systems over classical computers. For a full and clear demonstration of a ♯P-hardness of computing atom-excitation sampling, a condensate should be nonuniformly spread over an entire BEC trap and provide, via an interparticle interaction, multimode Bogoliubov coupling between a large number of excited atom states. Moreover, all of the above parameters of the many-body system should be controllable in a wide range to ensure sufficient variability of the observed joint occupation statistics of the excited states or coarse-grained groups of excited states. So, there is an open problem of designing BEC traps most suitable for experimental studies of various phenomena associated with atomic boson sampling.

The present paper is devoted to this problem: We discuss a basic model of a potential design of the multi-qubit BEC trap that could provide the required conditions and be particularly suitable for atomic-boson-sampling experiments. It is inspired by an analogy with multi-qubit or multi-qudit systems [31,32] and could look like a system of a finite number, Q, of single-qubit or -qudit cells shown in Figure 1 in a two-dimensional (2D) case.

Remarkably, a direct measurement of fluctuations in a total occupation of the noncondensate in cold dilute gases has already been achieved [33,34]. Splitting the noncondensate into some parts associated with the groups of excited states and measuring atom-number fluctuations in the occupations of those parts is the next important step in the many-body statistical physics toward testing quantum advantage. It is beyond the bulk of previous studies of the BEC phenomena, which is devoted to the mean properties of the condensate and quasiparticles, and it is also beyond the previous studies of fluctuations in the total occupation of the condensate (see, e.g., [34,35,36,37,38,39,40,41,42]). The atom-number fluctuations are especially important for the applications related to quantum information science and matter-wave interferometers [43], including Ramsey [44,45] and Mach–Zehnder [46] on-chip interferometers. In the literature, there are also other interesting discussions of the atom-number fluctuations associated with a subvolume of a BEC trap [47,48], BEC collapse [49], and squeezed states [44].

1.2. What Is the Atomic Boson Sampling?

In statistical theory, sampling is a selection of events (subsets) from within a sample space of all possible outcomes (or results or sample points) to mimic the characteristics of the probability distribution in a probability space (a probabilistic model). Atomic boson sampling means sampling from the excited-state occupations of identical Bose atoms subject to interparticle scattering (interaction) in a trapped Bose–Einstein condensed gas within a statistical ensemble of a given experimental setup. The atoms forming the condensate are not counted. One can consider the integer occupations $n_k=0,1,2,\dots$ and their joint probability distribution $\rho \left(\left\{n_k\right\}\right)$ for the individual orthonormal excited states $\{{\varphi }_k\left(\mathbf{r}\right)|k=1,2,\dots \}$, orthogonal to the condensate wave function, or for groups of such excited states. The simultaneous measurement of their occupations has to be completed by multiple detectors via projecting atoms onto preselected subsets (groups) of the excited states. The latter subsets determine the sampling probability distribution in question.

Condensed-matter statistical physics of a mesoscopic system of N atoms confined in a trap is highly nontrivial due to an interaction between massive atoms taking place on a background of the Bose–Einstein condensate formed by the same interacting atoms via spontaneous symmetry breaking at a critical temperature $T_c$. Quantum many-body fluctuations in this system remain ♯P-hard for computing [1] even in equilibrium and even within the grand-canonical-ensemble [50] and Bogoliubov–Popov approximations [51,52,53]. For simplicity’s sake, we adopt these approximations in the present paper and assume that the temperature is well below the critical region of the BEC phase transition, $T\ll T_c$.

The computational ♯P-hardness of atomic boson sampling is a real property of the interacting BEC gas, not just a feature of the Bogoliubov–Popov approximation. It follows from the exact non-perturbative theory of critical fluctuations in BEC, which is based on the non-polynomial diagram technique [54,55,56] and also leads to the representation of the joint probability distribution of the occupations of the bare excited atomic states in terms of the hafnian of a matrix associated with a correlation matrix which is similar to the one given by the stated approximation. In other words, the ♯P-hardness is a robust property of atomic boson sampling in the sense that it manifests itself (survives) even in the Bogoliubov–Popov, i.e., the first order with respect to the interaction parameter, mean-field approximation.

Each act of atomic boson sampling occurs by means of a measurement of the occupations of excited states or groups of them. One can employ, say, a simultaneous optical multi-detector imaging. Then, the system of interacting atoms returns back to its equilibrium state and becomes ready for the next act of sampling/measurement. In a sense, the system of atoms resets itself. This is true both if the atoms were reloaded into the trap after being released from the trap in the case of a trap-destruction measurement or if the atoms were not removed from the trap. Thus, the atomic boson sampler is not a quantum simulator of some input signal or some artificial, controlled process. The system of atoms in the BEC trap equipped with the atom-number detectors is just a quantum generator of random strings of atom-excitation occupations based on the natural process of persistent equilibrium fluctuations.

Surprisingly, atomic boson sampling has the potential to demonstrate a quantum advantage that is similar to the one of Gaussian boson sampling of non-interacting photons in a linear interferometer.

1.3. Comparison with Photonic Boson Sampling in a Linear Interferometer

Physics of photonic boson sampling in a linear interferometer, widely studied in the past decade [2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25], looks significantly simpler since photons are non-interacting, massless and enter the interferometer in a given quantum (e.g., Fock or squeezed) state from the external sophisticated on-demand light sources. Importantly, atomic boson sampling from an equilibrium BEC trap does not require external sources of atoms in a prescribed quantum state because the interacting atoms generate squeezing by themselves. Thus, in order to demonstrate the quantum advantage in a linear interferometer, one has to prepare the system of photons in a prescribed nonequilibrium state, whereas a usual thermal state of atoms in the BEC trap is enough in the case of atomic boson sampling. The reason for switching from the original proposal [3] to the Gaussian boson sampling protocol was an absence of the on-demand single-photon sources and availability of the reliable on-demand sources of input photons in the Gaussian/squeezed states based on the process of a parametric down-conversion [4,5,23].

The computational ♯P-hardness of the joint excited-state occupation fluctuations in the BEC trap exists by itself and does not require an adjustment or fine tuning of any input, processing, interaction or coupling parameters. On the contrary, Gaussian boson sampling in a linear interferometer, like other usually discussed nonequilibrium quantum processors and simulators, requires a lossless propagation through a system of beam couplers, splitters, and phase shifters as well as external sources of photons in squeezed states. Moreover, the main limiting factor for a demonstration of quantum advantage in photonic sampling is an exponential growth of photon losses on the input–output propagation with an increasing number of optical channels and coupling elements. Such a limiting factor is absent in the case of atomic boson sampling.

1.4. Why Atomic Boson Sampling Is ♯P-Hard for Computing?—A Brief Theory

Despite a number of outstanding differences stated above, both atomic and photonic boson samplings belong to the same ♯P-hard complexity class. The fact is that in both cases, the sampling (i.e., joint boson-occupation) probability distribution is determined by a hafnian (or, sometimes, a permanent) of matrices built from an appropriate covariance matrix G of the boson creation and annihilation operators. The above fact follows from the analytical calculation of the characteristic function (that is, Fourier transform) of the joint occupation probability distribution by means of the Wigner transform technique and application of the Hafnian Master Theorem that gives an explicit Taylor expansion of this characteristic function via the aforementioned hafnian [1,57]. It is well known that for a general $n\times n$ matrix, computing the hafnian or permanent is ♯P-hard [58,59,60], that is, requires an exponentially large number of operations $\sim n^3{2}^{n/2}$ or $\sim n^2{2}^n$, respectively, for the best general-purpose algorithms [60,61,62,63]. Contrary to a determinant that can be computed in polynomial time $\sim n^3$ via Gaussian elimination, the hafnian and permanent matrix functions are intimately related to the analysis of the ♯P-hard problems [26,64].

Yet, the computational ♯P-hardness (which is the basis of quantum advantage) of atomic boson sampling pertains only if the matrices under the hafnian and the covariance matrix G from which they are built are variable/controllable in a wide range by varying the trapping potential, interparticle interaction (via Feshbach resonance [65]), temperature, number of trapped atoms and other parameters of the BEC trap [66]. In other words, those matrices should not be restricted to a narrow subset of matrices for which the hafnian could be computed in polynomial time by some efficient approximation or algorithm, such as a fully polynomial random approximation scheme (FPRAS) [59], a recurrence method [67], etc.

In order to understand how to satisfy this requirement in the design of the BEC trap, we need to know how the covariance matrix G depends on the trap parameters. Fortunately, we have an analytical formula for the G via the matrix R of the Bogoliubov transformation,

$$G=RD{R}^{†}+\frac{1}{2}\left(R{R}^{†}-\mathbb{1}\right);G\equiv 〈:\left(\begin{array}{c}{\widehat{\mathbf{a}}}^{†}\\ \widehat{\mathbf{a}}\end{array}\right){\left(\begin{array}{c}{\widehat{\mathbf{a}}}^{†}\\ \widehat{\mathbf{a}}\end{array}\right)}^T:〉,D=\left[\begin{array}{cc}{D}_1& 0\\ 0& D_1\end{array}\right],D_1=\underset{j}{⨁}\frac{1}{e^{E_j/T}-1}.$$

(1)

Here, the boldface operator vectors ${\widehat{\mathbf{a}}}^{†}={({\widehat{a}}_1^{†},{\widehat{a}}_2^{†},\dots )}^T$ and $\widehat{\mathbf{a}}={({\widehat{a}}_1,{\widehat{a}}_2,\dots )}^T$ are the column-vectors of the creation and annihilation operators for bare excited atom states, respectively. The superscripts ${}^T$ and ${}^{†}$ denote transpose and conjugate transpose, respectively. The angles stand for the statistical average. The colons denote the normal ordering of operators, meaning that all creation operators stand to the left from the annihilation ones. The $2\times 2$ block-diagonal matrix D include the Bose–Einstein thermal occupation numbers for the quasiparticles of eigenenergies $\left\{E_j\right\}$. The matrix R performs the Bogoliubov transformation

$$\left(\begin{array}{c}{\widehat{\mathbf{a}}}^{†}\\ \widehat{\mathbf{a}}\end{array}\right)=R\left(\begin{array}{c}{\widehat{\mathbf{b}}}^{†}\\ \widehat{\mathbf{b}}\end{array}\right),R=\left[\begin{array}{cc}{U}^{*}& 0\\ 0& U\end{array}\right]\left[\begin{array}{cc}cosh{r}^{*}& {(e^{i\theta }sinhr)}^{*}\\ e^{i\theta }sinhr& coshr\end{array}\right],$$

(2)

from the representation of the excited-particle field operator in terms of the quasiparticle creation and annihilation operators ${\widehat{\mathbf{b}}}^{†}={({\widehat{b}}_1^{†},{\widehat{b}}_2^{†},\dots )}^T$ and $\widehat{\mathbf{b}}={({\widehat{b}}_1,{\widehat{b}}_2,\dots )}^T$ to the representation in terms of the bare-particle creation and annihilation operators, that is,

$${\widehat{\psi }}_{\mathrm{ex}}\left(\mathbf{r}\right)=\sum _{k\ne 0}{\varphi }_k\left(\mathbf{r}\right){\widehat{a}}_k=\sum _j\left(u_j\left(\mathbf{r}\right){\widehat{b}}_j+v_j^{*}\left(\mathbf{r}\right){\widehat{b}}_j^{†}\right).$$

(3)

Within the mean-field Bogoliubov–Popov approximation [51,52,53], adopted in the present paper, a Bose–Einstein-condensed gas is described via the Hamiltonian given by a quadratic form of the bare-particle creation and annihilation operators

$$\begin{array}{cc}\hfill & \widehat{H}=\sum _{k,k^{\prime }}\left({\widehat{a}}_k^{†}{\widehat{a}}_{k^{\prime }}+\frac{1}{2}\right)\int {\varphi }_k^{*}\left(-\frac{{ħ}^2\Delta }{2m}+U\left(\mathbf{r}\right)-\mu +2g\left(N_0{|{\varphi }_0\left(\mathbf{r}\right)|}^2+n_{\mathrm{ex}}\left(\mathbf{r}\right)\right)\right){\varphi }_{k^{\prime }}{d}^3\mathbf{r}\hfill \\ \hfill & +\frac{g{N}_0}{2}\sum _{k,k^{\prime }}{\widehat{a}}_k^{†}{\widehat{a}}_{k^{\prime }}^{†}\int {\varphi }_k^{*}{\varphi }_0^2{\varphi }_{k^{\prime }}^{*}{d}^3\mathbf{r}+\frac{g{N}_0}{2}\sum _{k,k^{\prime }}{\widehat{a}}_k{\widehat{a}}_{k^{\prime }}\int {\varphi }_k{\left({\varphi }_0^{*}\right)}^2{\varphi }_{k^{\prime }}{d}^3\mathbf{r}.\hfill \end{array}$$

(4)

Here, $\Delta$ is the three-dimensional Laplace operator, m is a particle mass, $U\left(\mathbf{r}\right)$ is an external potential, $g=4\pi {ħ}^{2}a/m$ is an interaction constant, a is a s-scattering length, $\mu$ is a chemical potential, $N_0$ is the mean number of particles in the condensate, $n_{\mathrm{ex}}\left(\mathbf{r}\right)$ is the mean density of the excited particle fraction (the noncondensate), and ${\varphi }_0\left(\mathbf{r}\right)$ is a condensate wave function normalized to unity, ${\int }_V|{\varphi }_0^2|d^3\mathbf{r}=1$. The superscript * means complex conjugation. The Hamilton operator (4) can be equivalently rewritten in the matrix form via the $(2\times 2)$-block Hamiltonian matrix H as follows

$$\widehat{H}={\left(\begin{array}{c}{\widehat{\mathbf{a}}}^{†}\\ \widehat{\mathbf{a}}\end{array}\right)}^{T}H\left(\begin{array}{c}{\widehat{\mathbf{a}}}^{†}\\ \widehat{\mathbf{a}}\end{array}\right),H=\left[\begin{array}{cc}\tilde{K}& K\\ K^{*}& {\tilde{K}}^{*}\end{array}\right].$$

(5)

The Bogoliubov transformation (2) diagonalizes the blocks K responsible for the co-rotating contributions ${\widehat{a}}_k^{†}{\widehat{a}}_{k^{\prime }}$ or ${\widehat{a}}_k{\widehat{a}}_{k^{\prime }}^{†}$ to the Hamiltonian and nullifies the blocks $\tilde{K}$ responsible for the counter-rotating contributions ${\widehat{a}}_k^{†}{\widehat{a}}_{k^{\prime }}^{†}$ or ${\widehat{a}}_k{\widehat{a}}_{k^{\prime }}$ to the Hamiltonian,

$$R^{T}HR=\left[\begin{array}{cc}0& E\\ E& 0\end{array}\right],E=\mathrm{diag}\{E_j|j=1,2,\dots \}.$$

(6)

The point is that the quasiparticles are the eigenstates of the Bogoliubov Hamiltonian with the eigenenergies $\{E_j|j=1,2,\dots \}$.

The multimode squeezing matrix [1,68,69,70,71,72,73] $r=\left(r_{k,k^{\prime }}\right)$, which is a positive semi-definite Hermitian matrix, as well as the unitary matrices U and $e^{i\theta }$, are calculated in [1]. Any additional unitary transformation V to another complete orthonormal set of excited states $\{{\varphi }_k^{\prime }\left(\mathbf{r}\right)|k=1,2,\dots \}$ in the single-particle Hilbert space,

$${\varphi }_k=\sum _{k^{\prime }\ne 0}{V}_{k^{\prime },k}{\varphi }_{k^{\prime }}^{\prime },$$

(7)

results in the Bogoliubov transformation $R^{\prime }$ which differs from the R in Equation (2) just by replacement of the unitary U with the composite unitary transformation $U^{\prime }=VU$.

Clearly, the Bogoliubov transformation (2) is a superposition of the squeezing and unitary transformations. In essence, their matrices r and U determine the complexity and variability of the covariance matrix G, since the classical thermal occupations of quasiparticles entering the matrix D are easy to compute. As a result, the ultimate reason for the ♯P-hardness (quantum advantage) of atomic boson sampling is an interplay between the squeezing (found in [74]) due to the interparticle interactions and interference due to the unitary mixing of bare-particle excited states. If either the squeezing or interference vanishes, then the joint occupation probabilities can be computed classically in polynomial time.

1.5. The Content of the Paper

Based on the main aspects of a truly hard for computing and innovative problem of atomic boson sampling in a BEC trap explained in the Introduction (Section 1), we formulate, in Section 2, general requirements for the BEC trap design facilitating testing the quantum advantage of atomic boson sampling in a full, controllable and clear way. A basic model of a multi-qubit BEC trap is devised in Section 3. In Section 4 and Section 5, we present the analytical and numerical results for the single-particle energy spectrum and eigenstates in the one- and two-dimensional multi-qubit traps, respectively. The corresponding solutions to the Gross–Pitaevskii equation for the condensate wave function are presented in Section 6. In Section 7, they are employed in the analysis of the Bogoliubov transformation and couplings for estimates of the multimode squeezing parameters. We discuss also the multimode dimensionality, squeezing and interference, which determine an asymptotic parameter for computational ♯P-hardness of atomic boson sampling and require analysis of the Bogoliubov–De Gennes equations for the quasiparticle spectrum and eigenfunctions. The experimental aspects of atomic boson sampling are discussed in Section 8. Section 9 contains concluding remarks.

2. General Requirements for the BEC Trap Design Facilitating Atomic Boson Sampling

Although any general-case BEC trap can be employed for studying manifestations of its potential quantum advantage over classical computing of boson sampling [1], in order to test this advantage in a controllable and unambiguous way, one should better use a specially designed trapping potential (see an example in Figure 1).

The challenge of the BEC trap design is twofold. On the one hand, it is desirable to have a trap with a finite (or even mesoscopic) number, M, of lower split-off excited states or groups of states, which are predominantly populated and strongly coupled to each other by means of Bogoliubov coupling. The atomic boson sampling could refer to the so-called marginal statistics—the quantum statistics of the occupations of these excited states irrespective to the occupations of all other states. It is especially informative if all of the higher excited states are separated from such a split-off miniband or sub-miniband of the selected M lower excited states or some groups of them by an energy gap $\Delta E$ wider than the temperature T and are not significantly coupled to the lower energy states. Then, these higher states are relatively poor populated, do not contribute to ♯P-hard complexity and can be skipped or accounted for as a kind of perturbation.

On the other hand, it is required to provide a way to simultaneously measure, that is to sample, the occupations of those M excited particle states or groups of them, say, by means of multi-detector optical imaging based on the light transmission through or scattering from the atomic cloud. Each detector should measure an appropriate occupation by projecting upon a prescribed state or group of states. Moreover, this subset of states or groups of states should be variable and controllable by means of tuning the detectors.

The geometry of a 2D multi-qubit BEC trap, such as the one shown in Figure 1, looks especially convenient for such boson sampling experiments since it allows one to implement multi-detector imaging by means of the laser light passing through the trap perpendicular to its plane. A controllable reconfiguration of the system of detectors aimed at varying the states or groups of states prescribed for occupation sampling also looks easier in the 2D geometry.

Suppose we design a multi-qubit BEC trap with a confining potential $U\left(\mathbf{r}\right)$ supporting a finite number Q of single-qubit cells, which form a 1D, 2D or 3D lattice and have two split-off lower energy levels each. Those two levels appear when a twofold-degenerate ground level is split by a certain perturbation. Then, such a lattice of single-qubit cells should be placed on top of a slightly varying in space background potential with high walls at the trap borders, and the inter-cell potential walls should be adjusted to be relatively high but narrow enough to allow for a quantum tunneling of atoms between cells. All this is necessary for establishment of the common to all cells nonuniform condensate and significant interaction between atoms from different single-qubit cells. The last two conditions are required for the existence of significant Bogoliubov coupling between a large number of excited states without which the multimode squeezing as well as interference via dressed quasiparticles are not well pronounced in the Bogoliubov transformation matrix R in Equation (2) and, hence, the computational ♯P-hardness disappears.

In fact, building a confining potential in the form of a single-qubit cell and duplicating it into a lattice is a relatively straightforward enterprise since such potentials are reminiscent of a double-well potential and an optical lattice potential, in which the BEC as well as the Bogoliubov excitations had been studied a lot [35,43,44,75,76,77,78,79,80,81]. The size of the multi-qubit trap depends on its dimensionality. In the 2D case of Figure 1, an overall dimension of the BEC trap is about $\sqrt{Q}\mathsf{\mu }$m, since each single-qubit cell has a scale of a de Broglie wavelength ∼1 $\mathsf{\mu }$m.

The starting point of our analysis is the limiting case of infinitely high inter-cell barriers and identical single-qubit cells, each with two single-particle eigenfunctions corresponding to the first and second energy levels, $e_1$ and $e_2$. These eigenfunctions form a natural basis for constructing the single-particle excited states of the actual trap. There are $2^Q$ combinations of these single-qubit states which are the eigenfunctions of the whole multi-qubit trap. Their $Q+1$ different energy levels $\{{\epsilon }_q=(Q-q)e_1+q{e}_2;q=0,1,\dots ,Q\}$ constitute a lower energy miniband. Degeneracies of levels are given by binomial coefficients $g_q=\left(\genfrac{}{}{0pt}{}{Q}{q}\right)$. Their sum coincides with the number of the single-qubit-state combinations: ${\sum }_{q=0}^Q\left(\genfrac{}{}{0pt}{}{Q}{q}\right)=2^Q$.

Those limiting-case eigenfunctions of the multi-qubit trap emerge adiabatically from the wavefunctions of an empty flat box trap when the inter-cell potential barriers are gradually introduced. This process can be easily understood within a 1D model of a flat background potential and almost equally spaced delta-function potential barriers. Its analysis shows that the limiting-case eigenfunctions of the finite lattice of independent qubit cells correspond to some superpositions of the $2Q$ lower-energy eigenfunctions of the whole trap with finite barriers. Hence, the $2Q$ lower-energy eigenfunctions ${\psi }_n\equiv {\psi }_{s,p}$ of order $n=p+sQ,p=1,\dots ,Q,s=0,1$, of the actual trap with finite barriers can be considered as a system of the “generating” eigenfunctions constituting two bands ($s=0,1$) and enumerated by the intra-band index $p=1,\dots ,Q$ and the band index $s=0,1,\dots$.

Although one could consider all $2^Q$ energy levels associated with Q qubits, it is more convenient to operate with a smaller number of the single-particle energy levels $M+1$ which constitute some miniband. In Section 4 and Section 5, we show that it is possible to choose the trap parameters in such a way that $M+1=2Q$ levels will form a lower miniband separated from the higher energy levels by an energy gap $\Delta E$ wider than the temperature T. Below, we mainly discuss the multi-qubit trap properties associated with such a miniband.

Analysis of the case when each cell has a larger number, $d>2$, of the lower split-off energy levels is very similar. Then, a finite lattice of such cells forms a multi-qudit trap.

In any case, the eigenfunctions ${\psi }_{s,p}$ and energy levels of the actual trap with a finite trapping potential can be easily controlled and varied in a wide range by means of controlling the background and barrier potentials as well as the dimensions of the single-qubit cells. For instance, the relative occupations of cells, i.e., the relative wavefunction amplitudes in different cells, within an eigenfunction of a given order can be individually controlled by tuning the cell background potentials. The intra-cell qubit properties, including the energy splittings $\delta E_j,j=1,\dots ,Q$, can be addressed by adjusting the intra-cell barriers.

The ground-state properties also can be controlled in this way. Implementing also control of the interparticle interactions via the Feshbach resonance [65], one can adjust the condensate wave function as needed. Below, we consider a favorable for the atomic boson sampling regime of a common condensate which is macroscopically occupied and inhomogeneously spread over the entire trap at a low temperature $T\ll T_c$ and a relatively large number of trapped atoms $N\gg Q$. At certain conditions, a particular number M (for example, $M=2Q-1$) of lower-miniband excited states can be considered as being decoupled from the continuum of excited states of the total infinite-size Hilbert space and constituting a finite-size Hilbert subspace. The situation could become especially clean and favorable for atomic boson sampling experiments if, in addition, the Bogoliubov couplings are adjusted to be spread over the whole lower miniband but not above the energy gap.

Apparently, the multi-qubit trap is capable of providing a whole series of other BEC regimes [78], starting from a strongly correlated regime and a regime of anomalous fluctuations in the critical region at $T\approx T_c$ to the regimes of fragmented condensates of the individual single-qubit cells and a quasi-condensate. However, their discussion goes beyond the scope of the present paper.

3. A Basic Model of a Multi-Qubit BEC Trap

The first step in designing the multi-qubit trap is to find its single-particle energy spectrum $\{{\epsilon }_n|n=0,1,\dots \}$ and eigenstates $\left\{{\psi }_n\right\}$, given by the linear Schr$\ddot{\mathrm{o}}$dinger equation

$$\left(-\frac{{ħ}^2}{2m}\Delta +U\left(\mathbf{r}\right)\right){\psi }_n\left(\mathbf{r}\right)={\epsilon }_n{\psi }_n\left(\mathbf{r}\right),$$

(8)

and adjust the trap parameters in order to fulfill the requirements on the split-off lower-energy miniband formulated above. The energies may be counted from the energy ${\epsilon }_0$ of a nondegenerate ground state $n=0$. An integer n orders all eigenstates in increasing energies ${\epsilon }_0<{\epsilon }_1\le {\epsilon }_2\le \cdots$ . Solutions to the single-particle Schr$\ddot{\mathrm{o}}$dinger Equation (8) provide a valuable starting point for the design of the multi-qubit trap—the zero-order approximation for the energies and wave functions of the excited states ($n=1,2,\cdots$) as well as the wave function of the ground state ($n=0$).

The second step is to find how the repulsive interparticle interaction modifies the ground state, that is, to find the condensate wave function ${\varphi }_0\left(\mathbf{r}\right)$ which obeys the Gross–Pitaevskii equation (the nonlinear Schr$\ddot{\mathrm{o}}$dinger equation) [35,52,53]

$$\left(-\frac{{ħ}^2\Delta }{2m}+U\left(\mathbf{r}\right)+g{N}_0{|{\varphi }_0\left(\mathbf{r}\right)|}^2+2g{n}_{\mathrm{ex}}\left(\mathbf{r}\right)-\mu \right){\varphi }_0=0,g=\frac{4\pi {ħ}^{2}a}{m}.$$

(9)

The goal is to verify the presence of a non-fragmented condensate, which is common for the entire trap and spreading over all single-qubit cells. A non-uniformity of the condensate should be controllable by adjusting the trap parameters. Accurate knowledge of the condensate wave function is necessary for calculating the Bogoliubov couplings

$${\Delta }_{k{k}^{\prime }}=g{N}_0\int {\varphi }_k^{*}\left(\mathbf{r}\right){|{\varphi }_0\left(\mathbf{r}\right)|}^2{\varphi }_{k^{\prime }}\left(\mathbf{r}\right)d^3\mathbf{r},{\tilde{\Delta }}_{k{k}^{\prime }}=g{N}_0\int {\varphi }_k^{*}\left(\mathbf{r}\right){\varphi }_0{\left(\mathbf{r}\right)}^2{\varphi }_{k^{\prime }}^{*}\left(\mathbf{r}\right)d^3\mathbf{r}$$

(10)

between the preselected bare-particle excited states $\{{\varphi }_k|k=1,2,\dots \}$ and making sure that they are well pronounced for a large enough number of these states as per requirements stated in Section 2. If most of atoms are in the condensate, $N_0\approx N$, then a characteristic length of a condensate inhomogeneity is equal to a so-called healing length

$$\xi =\frac{ħ}{\sqrt{2mgN/V}}=\frac{1}{\sqrt{8\pi aN/V}}.$$

(11)

The Gross–Pitaevskii equation, as a mean field approximation, is valid if an average distance d between atoms is small compared to the healing length,

$$d\ll \xi .$$

(12)

The next step involves solving the Bogoliubov–De Gennes equations for the quasiparticle spectrum and eigenfunctions as well as calculation of the Bogoliubov transformation matrix, squeezing and other parameters describing the joint probability distribution of the excited atom occupations and atomic boson sampling. We just briefly comment on this step in Section 7 and Section 8, since the analysis of quasiparticles goes beyond the scope of this article.

In the present paper, we limit ourselves to the first two steps and calculation of Bogoliubov couplings responsible for interparticle interactions in the Bogoliubov Hamiltonian.

For the sake of clarity and simplicity, we consider only a simple basic model of the multi-qubit BEC trap illustrated in Figure 2: namely, a one-dimensional (1D) or two-dimensional (2D) array of a finite number Q of the single-qubit cells. In the case of a 1D chain of the single-qubit cells, each q-th single-qubit cell includes two flat background potentials $U_{2q-1},U_{2q}$ and a delta-function potential ${\beta }_q\delta (x-x_q)$ located near its center, while the cells are separated by the delta-function potential walls $\{{\alpha }_q\delta (x-X_q)\ge 0|q=1,2,\dots ,Q-1\}$ and ordered along the x axis so that $0=X_0\le x_1\le X_1\le x_2\le X_2\le \dots \le x_Q\le X_Q=QL$. The corresponding 1D trapping potential is modeled as follows

$$\begin{array}{cc}\hfill & U\left(x\right)=\sum _{q=1}^Q\{U_{2q-1}[\theta (x-X_{q-1})-\theta (x-x_q)]+U_{2q}[\theta (x-x_q)-\theta (x-X_q)]\hfill \\ \hfill & +{\beta }_q\delta (x-x_q)+{\alpha }_q\delta (x-X_q)\}\mathrm{if}x\in (0,QL);\hfill \\ \hfill & U\left(x\right)=\infty \mathrm{if}x\le 0\mathrm{or}x\ge QL.\hfill \end{array}$$

(13)

The amplitudes of the background potentials $\{U_j|j=1,\dots ,2Q\}$ and all delta-function potentials $\left\{{\alpha }_q\right\},\left\{{\beta }_q\right\}$ as well as their locations $\left\{x_q\right\},\left\{X_q\right\}$ could be different for different single-qubit cells and constitute a set of controllable parameters of the multi-qubit BEC trap; $\delta \left(x\right)$ is the Dirac delta function, and $\theta \left(x\right)$ is the unit step function: $\theta \left(x\right)=0$ if $x<0$, $\theta \left(x\right)=1$ if $x\ge 0$.

In the case of a 2D square $Q_1\times Q_1$ array of $Q={\left(Q_1\right)}^2$ single-qubit cells, we adopt a model potential $U(x,y)=U\left(x\right)+U^{\prime }\left(y\right)$ given by a sum of two 1D potentials along the axes x and y, each of which being similar to the 1D potential in Equation (13):

$$\begin{array}{cc}\hfill & U\left(x\right)=\sum _{q=1}^{Q_1}\{U_{2q-1}[\theta (x-X_{q-1})-\theta (x-x_q)]+U_{2q}[\theta (x-x_q)-\theta (x-X_q)]\hfill \\ \hfill & +{\beta }_q\delta (x-x_q)+{\alpha }_q\delta (x-X_q)\}\mathrm{if}x\in (0,Q_{1}L);\hfill \\ \hfill & U^{\prime }\left(y\right)=\sum _{q=1}^{Q_1}\{U_{2q-1}^{\prime }[\theta (y-Y_{q-1})-\theta (y-y_q)]+U_{2q}^{\prime }[\theta (y-y_q)-\theta (y-Y_q)]\hfill \\ \hfill & +{\beta }_q^{\prime }\delta (y-y_q)+{\alpha }_q^{\prime }\delta (y-Y_q)\}\mathrm{if}y\in (0,Q_1{L}^{\prime }).\hfill \end{array}$$

(14)

Again, we set the potential to be infinitely high beyond the outer borders of the entire multi-qubit trap: $U\left(x\right)=\infty \mathrm{if}x\le 0\mathrm{or}x\ge Q_{1}L,U^{\prime }\left(y\right)=\infty \mathrm{if}y\le 0\mathrm{or}y\ge Q_1{L}^{\prime }$. The amplitudes of the background potentials $\{U_j|j=1,\dots ,2Q\},\{U_j^{\prime }|j=1,\dots ,2Q\}$ and all delta-function potentials $\left\{{\alpha }_q\right\},\left\{{\beta }_q\right\},\left\{{\alpha }_q^{\prime }\right\},\left\{{\beta }_q^{\prime }\right\}$ as well as their locations $\left\{x_q\right\},\left\{X_q\right\},\left\{y_q\right\},\left\{Y_q\right\}$ constitute a set of controllable parameters of the 2D multi-qubit BEC trap.

Modeling the confining potential by piecewise flat and delta-function potentials is a well-justified textbook approach pertinent to the analysis of the effects of tunneling, reflection and trapping of particles by potential barriers and walls on the wave functions and energy spectrum in quantum mechanics (see, e.g., [82,83,84,85] and references therein). It is consistent with the well-known facts that (a) the Rayleigh–Ritz characterization of the eigen energies involves only weighted averages of the potential and (b) the multiple-scale perturbation theory yields the correct leading-order asymptotics within the piecewise-flat-potentials approximation [86]. The main quantities in question for the analysis in the present paper are the condensate wave function and Bogoliubov couplings, which determine the ultimate result for the covariance matrix and statistics of atomic boson sampling. Their representativeness and robustness with respect to the adopted modeling by the piecewise flat and delta-function potentials are predetermined by the nature of the Bogoliubov couplings (10) as the overlapping integrals which do not depend significantly on a jump in the value of the first or second derivative of the wave function originated from the presence of the delta- or step-function, respectively, in the external potential. Furthermore, the actual potential in the Gross–Pitaevskii and Bogoliubov–de Gennes Equations (9) and (34) is always curved by the interparticle–interaction contribution $g{N}_0{\varphi }_0^2\left(\mathbf{r}\right)$ proportional to the continuous condensate occupation $|{\varphi }_0^2\left(\mathbf{r}\right)|$. Obviously, in an experimental setting, a non-flat background potential will lead to qualitatively the same results.

4. One-Dimensional Multi-Qubit Trap: Single-Particle Eigen Functions and Energies

Consider a 1D trap with the model potential (13). The basic model adopted above allows one to solve the 1D Schr$\ddot{\mathrm{o}}$dinger Equation (8),

$$\left(-\frac{{ħ}^2}{2m}\frac{d^2}{d^{2}x}+U\left(x\right)\right){\psi }_n\left(x\right)={\epsilon }_n{\psi }_n\left(x\right),$$

(15)

analytically and easily find the single-particle energy spectrum and wave eigenfunctions. In this section, we demonstrate the single-particle properties of the 1D multi-qubit traps in a series of generic examples.

4.1. Asymmetric 1D Single-Qubit Trap: Explicit Solution for a Double-Well Trap

The solution to Equation (15) for the eigen functions and energies of an asymmetric 1D single-qubit trap described by the model (13) of a double-well trap with the intra-cell delta-function potential of a magnitude $\beta$ located at a position $x_1=\eta L,\eta \in (0,1)$, is elementary:

$${\psi }_n\left(x\right)=Asin\left(k_{n}x\right)\mathrm{if}0\le x\le \eta L,{\psi }_n\left(x\right)=\frac{Asin\left(\eta k_{n}L\right)}{sin\left[(1-\eta )k_{n}L\right]}sin\left[k_n(L-x)\right]\mathrm{if}\eta L\le x\le L,$$

(16)

$$sin\left(\eta k_{n}L\right)sin\left[(1-\eta )k_{n}L\right]=-\frac{{ħ}^2{k}_n}{2m\beta }sin\left(k_{n}L\right),{\epsilon }_n=\frac{{ħ}^2{k}_n^2}{2m}.$$

(17)

Here, A is an appropriate normalization constant. The dependence of the first six eigen wave numbers $k_n$ on the asymmetry parameter $\eta \in (0,1)$ is illustrated in Figure 3. Note a very narrow energy splitting ${\epsilon }_2-{\epsilon }_1\ll {\epsilon }_2$ between the two lower excited states $n=1,2$ and a very wide energy gap ${\epsilon }_3-{\epsilon }_2\approx 3{\epsilon }_2$ separating them from the next two excited states $n=3,4$ in the case of the central, symmetric location of the intra-cell delta-function potential, $\eta =1/2$. With an increasing asymmetry, the energy ladder experiences a significant restructuring. For example, if the asymmetry is $\eta \approx 1/3$ or $2/3$, then already, the three lower (the lowest first, ${\epsilon }_1$, and two very close second and third, ${\epsilon }_2\approx {\epsilon }_3$) energy levels are separated from the higher energy levels by an energy gap ${\epsilon }_4-{\epsilon }_3\approx 1.3{\epsilon }_3$.

4.2. Symmetric 1D Two-Qubit Trap: Even versus Odd Eigenfunctions and Their Eigenenergies

Consider the symmetric two-qubit trap (13) with the central locations of the intra-cell delta-function potentials of equal magnitude ${\beta }_1={\beta }_2\equiv \beta$ at $x_1=L/2$, $x_2=3L/2$ and the inter-cell delta-function potential of the magnitude ${\alpha }_1\equiv \alpha$ at $X_1=L$ in the absence of the background potential, $U_1=U_2=0$. The odd wave functions, which have the odd spatial symmetry relative to the center of the trap, equal zero at the trap center and are not affected by the inter-cell potential wall. Obviously, solutions for them are reduced to the single-qubit trap solution (16) in each of two single-qubit cells. For example, the odd wave function in the left single-qubit cell is

$${\psi }_n\left(x\right)=Asin\left(k_{n}x\right)\mathrm{if}0\le x\le L/2,{\psi }_n\left(x\right)=Asin\left[k_n(L-x)\right]\mathrm{if}L/2\le x\le L.$$

(18)

Hence, the dimensionless energy spectrum, ${\overline{\epsilon }}_n=(2m{L}^2/{ħ}^2){\epsilon }_n$, of the odd eigenfunctions is given by Equation (17), that is

$$k_{n}Lcos(k_{n}L/2)+\overline{\beta }sin(k_{n}L/2)=0,{\overline{\epsilon }}_n={\left(k_{n}L\right)}^2;\overline{\alpha }=\frac{\alpha mL}{{ħ}^2},\overline{\beta }=\frac{\beta mL}{{ħ}^2}.$$

(19)

A dimensionless parameter $\overline{\beta }$ describes the effect of the intra-cell delta-function potential.

The solution to Equation (15) for the even eigenfunctions is more involved:

$$\begin{array}{cc}\hfill & {\psi }_n\left(x\right)=A_{1}sin\left(k_{n}x\right)\mathrm{if}0\le x\le L/2,\hfill \\ \hfill & {\psi }_n\left(x\right)=A_{2}sin[k_n(L-x)+\phi ]\mathrm{if}L/2\le x\le L,\hfill \\ \hfill & {\psi }_n\left(x\right)=A_{2}sin[k_n(x-L)+\phi ]\mathrm{if}L\le x\le 3L/2,\hfill \\ \hfill & {\psi }_n\left(x\right)=A_{1}sin\left[k_n(2L-x)\right]\mathrm{if}3L/2\le x\le 2L.\hfill \end{array}$$

(20)

It includes two normalization constants $A_1,A_2$, the phase shift $\phi$ (such that $tan\phi =k_{n}L/\overline{\alpha }$) and depends on the inter-cell delta-function potential via the dimensionless parameter $\overline{\alpha }=\alpha mL/{ħ}^2$. The energy spectrum ${\overline{\epsilon }}_n={\left(k_{n}L\right)}^2$ of the even eigenfunctions is determined by the eigen wave number $k_n$ that can be found from the explicit transcendental equation:

$$k_{n}L\left[k_{n}Lcos\left(k_{n}L\right)+\overline{\beta }sin\left(k_{n}L\right)\right]+2\overline{\alpha }sin\left(\frac{k_{n}L}{2}\right)\left[k_{n}Lcos\left(\frac{k_{n}L}{2}\right)+\overline{\beta }sin\left(\frac{k_{n}L}{2}\right)\right]=0.$$

(21)

Figure 4 shows clearly the full structure of energy spectrum of a two-qubit trap. Firstly, one can see the unperturbed energy level spread for an empty rectangular well when the inter-cell and intra-cell potentials are zero. This behavior can be modulated in two ways, by increasing either of the two potentials. There is almost a complete symmetry between how these two potentials effect the energy level structure, with the only difference coming in the even-numbered energy levels. While every fourth level is totally unperturbed by both potentials, the other even energy levels only see the intra-cell potentials, as these functions are zero at the center of the well. This leads to an asymmetry in the structure, which affects the orange-colored energy levels in the figure.

In addition to this asymmetry, the overall structure of the two-qubit-trap energy spectrum is largely determined by the formation of minibands. If either of the two potentials are individually raised to be large, four sub-minibands are formed, while the raising of both potentials leads to the formation of two minibands, with a large energy gap between the first four and second four energy levels. This can be quantitatively measured by taking the ratio of the energy separation between the first and second miniband with the energy width of the first miniband (see Figure 5). This will demonstrate how easy it will be to set the temperature such that the lower miniband is fully populated while the higher has little to no occupation.This figure of merit must be balanced against the necessity that atoms are still able to easily move between cells, requiring that the potential barriers not be too high.

4.3. Four-Qubit Chain of Identical Symmetric Single-Qubit Cells: Hierarchy of Even/Odd Solutions

The analysis presented in Section 4.2 can be easily generalized to the case of the four-qubit trap (13) with similarly symmetric parameters ${\alpha }_q=\alpha ,{\beta }_q=\beta ,x_q=(q-1/2)L,$$X_q=qL,U_{2q-1}=U_{2q}=0$$\forall q=1,2,3,4$ and total length $4L$. Again, the solutions for the wave functions with the odd spatial symmetry relative to the center of the trap and their energy spectrum are reduced to the solutions for the half trap, that is, for the two-qubit trap and, hence, are given (say, for the left half of the four-qubit trap) by Equations (18)–(21). The only novel element of the analysis is the solution for the even eigenfunctions. It has the following form in the left half of the four-qubit trap

$$\begin{array}{cc}\hfill & {\psi }_n\left(x\right)=A_{1}sin\left(k_{n}x\right)\mathrm{if}0\le x\le L/2,\hfill \\ \hfill & {\psi }_n\left(x\right)=A_{2}sin[k_n(L-x)+{\phi }_2]\mathrm{if}L/2\le x\le L,\hfill \\ \hfill & {\psi }_n\left(x\right)=A_{3}sin[k_n(x-L)+{\phi }_3]\mathrm{if}L\le x\le 3L/2,\hfill \\ \hfill & {\psi }_n\left(x\right)=A_{4}sin[k_n(2L-x)+{\phi }_4]\mathrm{if}3L/2\le x\le 2L,\hfill \end{array}$$

(22)

with the same form of equations being found reflected across the center of the trap at $x=2L$. Now, it includes four normalization constants $A_1,A_2,A_3,A_4$, three phase shifts ${\phi }_2,{\phi }_3,{\phi }_4$ and the eigen wave number $k_n$. The latter four quantities can be found from the following four equations expressing a discontinuity of the wave-function derivative across each delta-function potential barrier:

$$\begin{array}{cc}\hfill & k_{n}L[cot\left({\phi }_4\right)+cot\left({\phi }_4\right)]=2\overline{\alpha },\hfill \\ \hfill & k_{n}L[cot({\phi }_3+k_{n}L/2)+cot({\phi }_4+k_{n}L/2)]=-2\overline{\beta },\hfill \\ \hfill & k_{n}L[cot\left({\phi }_2\right)+cot\left({\phi }_3\right)]=2\overline{\alpha },\hfill \\ \hfill & k_{n}L[cot(k_{n}L/2)+cot(k_{n}L/2+{\phi }_2)]=-2\overline{\beta }.\hfill \end{array}$$

(23)

Excluding the phase shifts, we arrive to the explicit transcendental equation,

$$\begin{array}{l}{(k_{n}L/2)}^3[\overline{\beta }sin\left(2k_{n}L\right)+k_{n}Lcos\left(2k_{n}L\right)]\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{(k_{n}L/2)}^2(2\overline{\alpha }+\overline{\beta })sin\left(k_{n}L\right)[\overline{\beta }sin\left(k_{n}L\right)+k_{n}Lcos\left(k_{n}L\right)]\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\overline{\alpha }}^2{sin}^2(k_{n}L/2){[\overline{\beta }sin(k_{n}L/2)+k_{n}Lsin(k_{n}L/2)]}^2\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+k_{n}L\overline{\alpha }\overline{\beta }{sin}^2(k_{n}L/2)[\overline{\beta }sin\left(k_{n}L\right)+k_{n}Lcos\left(k_{n}L\right)]=0,\end{array}$$

(24)

for the eigen wave number $k_n$ which determines the dimensionless energy spectrum ${\overline{\epsilon }}_n={\left(k_{n}L\right)}^2$ of the even eigenfunctions.

The entire energy spectrum is illustrated in Figure 6 by dependence of the first sixteen lower energy levels on the inter-cell and intra-cell delta-function potentials. As expected, it is similar to the analogous dependence for the two-qubit trap shown in Figure 4. Again, on the far left, the unperturbed energy levels can be observed and compared to the energy levels on the far right, which show clearly the expected miniband behavior, with a larger gap between the strongly grouped first eight and second eight energy levels.

Let us look again at the ratio of the energy gap between the first and second minibands and the energy width of the first miniband in Figure 5. Although this ratio is smaller for similar delta-function potentials of the two-qubit trap, the doubling of the available energy levels is a strong advantage. The degradation is not significant, as to once again achieve a ratio of about $3.5$, we only need to go from a dimensionless delta-function potential magnitude of 10 to about 12.

Apparently, the even/odd hierarchy of solutions revealed above is suggestive for an extension to any 1D chain of $Q=2^p,p=1,2,3,\dots$, identical symmetric single-qubit cells.

4.4. Multi-Qubit Chain of Q Identical Single-Qubit Cells: Asymptotics of Zeroes and Miniband of 2Q Energy Levels

Consider a 1D chain of Q identical symmetric single-qubit cells, each with a zero flat potential and of length L, separated by delta-function potential walls of the same amplitude ${\alpha }_q=\alpha$. The system is placed inside an infinitely high box potential well of length $QL$. Assume that each single-qubit cell contains a delta-function potential of the equal amplitude, ${\beta }_q=\beta$, placed at the center of the cell and perturbing its energy levels.

As was explained in Section 2, the eigenfunctions of the multi-qubit trap {${\psi }_{s,p}\left(x\right)|p=1,\cdots ,Q;s=0,1,\cdots$} can be considered as arising adiabatically with increasing delta-function potentials $\alpha$ and $\beta$ from the sinusoidal wave eigenfunctions of a box trap with a zero flat potential and of length $QL$,

$${\psi }_n^{\left(0\right)}\left(x\right)=\sqrt{\frac{2}{QL}}sin\left(\frac{n\pi x}{QL}\right),{\epsilon }_n^{\left(0\right)}=\frac{{(ħ\pi n)}^2}{2m{\left(QL\right)}^2},n=p+sQ;p=1,\cdots ,Q;s=0,1,\cdots .$$

(25)

The index n is equal to the number of half-wavelength variations between the borders of the entire box trap and orders the wave functions in accord with the linearly growing numbers of zeroes, $n-1$, and quadratically growing energies, ${\epsilon }_n^{\left(0\right)}$. The index $s=0,1,\cdots$ enumerates the bands. Each band consists of Q eigenfunctions enumerated by the intra-band index $p=1,\cdots ,Q$.

We find a general asymptotic rule: When ${\alpha }_q\to \infty$, all Q eigenfunctions ${\psi }_{s,p}$ within a given s-band have exactly the same (equal to the band order s) number of zeroes inside each single-qubit cell. The only exception constitutes such single-qubit cells and such eigenfunctions ${\psi }_{s,p}$ for which there is just one zero of the corresponding sinusoidal eigenfunction ${\psi }_{p+sQ}^{\left(0\right)}$ located exactly at the center of a single-qubit cell. This is an exceptional, degenerate case of a node frustration when neither of the two delta-function walls of the single-qubit cell are able to shift the location of this zero toward (underneath) its (wall’s) location with increasing delta-potential ${\alpha }_q\to \infty$. The amplitude of this eigenfunction tends to zero everywhere inside such an exceptional single-qubit cell.

The remarkable asymptotic behavior stated above is a consequence of the fact that the eigenfunctions ${\psi }_{s,p}\left(x\right)$ tend to zero at the positions of the inter-cell delta-potential walls with increasing magnitude of the delta-function potential: ${\psi }_{s,p}(x=jL)\to 0\mathrm{at}{\alpha }_q\to \infty$ for $j=1,\dots ,Q-1$. This occurs via two mechanisms. A delta-potential wall either (i) gradually digs a deep dip forcing the eigenfunctions to approach zero at the wall location, or (ii) gradually shifts the closest-to-the-wall zero of the sinusoidal wave eigenfunction ${\psi }_{p+sQ}^{\left(0\right)}$ to (underneath) the wall location. Accordingly, the eigenfunctions do not or do change their sign across the delta-potential wall. It is illustrated in Figure 7(left), where both mechanisms of asymptotics formation are clearly represented. In particular, the inter-cell delta-potential walls at the dimensionless positions $x/L=2$ and $x/L=6$ implement the first mechanism on the eigenfunctions ${\psi }_{s,p=1}$ (blue) and ${\psi }_{s,p=2}$ (yellow), the second mechanism on the eigenfunction ${\psi }_{s,p=3}$ (green), and do not affect the eigenfunction ${\psi }_{s,p=4}$ (red) whose sinusoidal counterpart ${\psi }_{4+sQ}^{\left(0\right)}$ is already equal to zero at the wall locations. Another situation when the above two mechanisms clearly manifest themselves is discussed in the next Section 4.5 in regard to Figure 11.

The sinusoidal wave functions ${\psi }_n^{\left(0\right)}$ of the higher band orders $s\ge 2$, i.e., $n>2Q$, have three or more zeroes at least within one, say j-th, single-qubit cell, that is, within the interval $x\in [qL,(q+1\left)L\right]$. So, at ${\alpha }_q\to \infty \forall q$, they turn into the eigenfunctions ${\psi }_{s,p},s\ge 2$, which has the s (two or more) zeroes inside each single-qubit cell and cannot be associated with the $2^Q$ eigenfunctions of the multi-qubit lowest miniband that have no more than one zero inside a single-qubit cell. Hence, only the first two bands of the eigenfunctions, ${\psi }_{0,p}$ and ${\psi }_{1,p}$, are relevant to the wave function superpositions that asymptotically yield the $2^Q$ eigenfunctions of the multi-qubit trap (i.e., combinations of the single-qubit states with the energies within the miniband).

Thus, we focus below on the analysis of the miniband of the first $2Q$ energy levels corresponding to the first two bands, $s=0,1$, of the eigenfunctions ${\psi }_{s,p}$. In principle, we can build a new system of Q qubits assigning an arbitrary pair of eigenfunctions ${\psi }_{0,p}$ and ${\psi }_{1,p^{\prime }}$ to be the lower and upper energy states of any new qubit. The most natural system of Q qubits will be formed by the pairs of eigenfunctions $\{{\psi }_{0,p},{\psi }_{1,p}|p=1,\dots ,Q\}$ with equal indices $p^{\prime }=p$. Then, we again can consider their $2^Q$ multi-qubit combinations of the eigenfunctions of the miniband of the lowest $2Q$ energy levels in the multi-qubit trap with arbitrary finite (not necessary infinite) inter-cell potential walls. These qubits are not identical anymore, even if their lengths $L_q,q=1,\dots ,Q$ are the same.

4.5. Multi-Qubit Chain of Significantly Different Single-Qubit Cells: Control of Occupations, Energies

If the single-qubit cells in the multi-qubit chain are not identical, then the excited-state wave functions become less symmetric. However, by controlling the background flat potentials $\{U_j|j=1,\dots ,2Q\}$ in Equation (13), one can make the ground state more uniform. A typical example of spatial profiles of the ground-state wave function and three lower-energy excited-state eigenfunctions is shown in Figure 8. In this figure and throughout the present paper as in Equation (19), a bar above a symbol of the potential or other energy quantity denotes its dimensionless value in terms of the energy unit ${ħ}^2/\left(2m{L}^2\right)$ where L is the length of a typical single-qubit cell, that is,

$${\overline{U}}_j=(2m{L}^2/{ħ}^2)U_j,{\overline{\epsilon }}_n=(2m{L}^2/{ħ}^2){\epsilon }_n.$$

(26)

Relatively large variations in the lengths, $L_q$, and background flat potentials, $U_{2q-1},U_{2q}$, of different single-qubit cells allow one to control and vary the trap eigenfunctions ${\psi }_{s,p}$ and energy levels in a wide range. Adjusting separately the flat potentials $U_j$ in different single-qubit cells allows one to control individually the relative amplitudes of eigenfunctions in different cells, that is, in particular, the relative occupations of different single-qubit cells as is illustrated in the right column of Figure 7.

The general structure and asymptotic behavior of the trap eigenfunctions described above for the chain of identical single-qubit cells is robust (remains qualitatively the same) with respect to small variations of the trap parameters. However, large variations change the picture. In particular, the single-qubit cells of significantly different lengths could acquire different numbers of eigenfunction zeroes per a single-qubit cell even within the same band of eigenfunctions as is illustrated in Figure 9a.

At last, tuning the intra-cell delta-function potentials ${\beta }_q,q=1,\dots ,Q$, provides one more tool for controlling and varying the profile and energy spectrum of the multi-qubit trap eigenfunctions ${\psi }_{s,p}\left(x\right)$. As is illustrated in Figure 9b, it affects the eigenfunctions in the first, $s=0$, band stronger than the eigenfunctions in the second, $s=1$, band. Thus, it is an efficient tool for controlling the intra-qubit properties, in particular, the qubit energy splittings $\delta E_q,q=1,\dots ,Q$.

One other possibility when constructing these multi-qubit traps is to place the inter-cell and intra-cell delta potentials in such a way as to break the symmetry between each qubit cell. This can act as yet another knob by which to control the diversity of the system along with the heights of the delta potentials and the modulation of the background potential. Delta-function potentials can be moved individually or following some group pattern. The position of the potential can be represented by a number from 0 to 1, essentially what percentage of the cell is traveled starting from the center of the cell before dropping down the delta-function potential. Suppose the intra-cell barriers are shifted leftwards on the left side of the four-qubit trap and an equal distance rightward on the right side. The effect of this shifting on the energy levels of a four-qubit trap can be seen in Figure 10. It is not symmetric and can be understood by observing where the delta-function potentials are modulating the unperturbed wave functions for each energy level. For example, take the ninth energy level seen in cyan in Figure 10. This energy level is maximized around $0.18$ and minimized around $0.36$. When the positions of the delta-function potentials are overlaid on the unperturbed ninth eigenfunction, it becomes clear why this is the case. At the maximum effect, the delta-function potentials are positioned near the maxima of the unperturbed function, while at the minimum effect, the delta potentials are placed instead near the zeros of the unperturbed function, as seen in Figure 11. The latter figure illustrates also the two mechanisms of a wave-function pertubation stated in Section 4.4. In the former case, the first mechanism, namely, digging a dip in a function profile underneath the delta-function potential, takes place. In the latter case, the second mechanism, namely, dragging the nearest node underneath the delta-function potential, takes place. The respective effects on the energy of the eigenfunctions are very different as is shown in Figure 11. Note that the digging effect of the central inter-cell delta-function potential is the same in both cases.

5. Two-Dimensional Multi-Qubit Trap: Single-Particle Eigen Functions and Energies

Here, we describe a simple 2D model of the multi-qubit trap formed by a potential, $U(x,y)=U_x\left(x\right)+U_y\left(y\right)$, which is the sum of the two 1D potentials considered above. In this case, the solution to the 2D single-particle Schr$\ddot{\mathrm{o}}$dinger Equation (8) is reduced, via a factorization, to the solutions to the 1D Schr$\ddot{\mathrm{o}}$dinger equation described above.

As a generic example, let us consider the 2D symmetric four-qubit trap, in which both the potential in the x and y plane are exactly the same. Both will share the same inter-cell and intra-cell delta-function potential strengths. As the energy levels are known for the 1D case (see, for example, Figure 6), constructing the energy levels for the 2D case is a simple task, requiring only for the individual energies to be added in every possible combination. A visualization of the energy levels created by combining the first nine energy levels in each dimension is shown in Figure 12 for the symmetric 2D $(4\times 4)$-qubit trap. This fully covers the first miniband of each dimension plus the first energy level in the second miniband.

This structure is interesting. Firstly, let us understand why there is a noticeable asymmetry between the effects of the inter-cell and intra-cell potentials. We see that for the two different axes, the first set of energy levels that are bunched together contains either nine levels for intra-cell potentials or 16 levels for inter-cell potentials. This comes from the fact that for the 1D case of a large intra-cell potential with little inter-cell potential, there is a splitting of the energy levels in the first miniband into two parts, one containing the first three with the other containing the remaining five. These two sub-minibands have a significant enough energy gap between them so that when the 2D plot is created, we see the formation of three sub-minibands, one containing combinations consisting of only energy levels in the first half of the sub-miniband, one containing the crossover terms, and the final piece containing the combinations consisting of only energy levels in the second half of the sub-miniband. Because there are nine possible combinations for the first half, we see nine energy levels in the first band of the plot. However, for the case of the large inter-cell potential, we instead see a 4–4 split of the energy levels in the miniband rather than the 3–5 split in the intra-cell case. Thus, 16 possible combinations of the energy levels in the first half of the miniband.

To understand why there is a difference in splitting of the energy levels depending on if the inter-cell or intra-cell potentials dominate the trap, we need only to look at the fourth energy level and how it is affected by each of the two different types of traps. The easiest to understand is the inter-cell potential dominant traps. In this case, the fourth eigenfunction is almost totally unperturbed by the delta-function potentials, as its natural nodes are already placed at the locations of the inter-cell potential barriers, while the lower three energy levels are pulled up toward the fourth. Likewise, the fifth, sixth, and seventh energy levels are bought up to the totally unperturbed eighth, leading to the 4–4 split structure we observe. However, in the case of the intra-cell potentials dominating the trap, the natural nodes of the fourth energy level wavefunction are placed directly between two of the intra-cell potential walls. This, as described in Section 4.4, leads to a situation where new nodes must be created to accommodate the large potentials. These new nodes drastically increase the average derivative, bringing the energy level much further above the third energy level below it, whose eigenfunction is able to shift its nodes to fall under the existing large potentials without much effect to its energy level. This effect is illustrated in Figure 13.

The last notable aspect of Figure 12 is the crossing of energy levels that can be seen near this origin. This behavior arises from the fact that while initially, the energy levels are approximately evenly spread, once the potentials start ramping up, there are some significant gaps created in the miniband structure. Thus, for low potentials, the energy level created by combining the first and fifth energy levels may be lower than that created by the fourth energy level combined with itself. However, once the fifth energy level is drastically raised by the introduction of the potentials, the combination of the first and fifth levels will increase its total energy above that of the double fourth, at least for the inter-cell dominant case where the fourth energy level is mostly unperturbed. This sort of behavior is only seen near the origin where the energy levels change drastically with the introduction of the delta-function potentials, as once the overall structure begins to form, there is not a significant enough change to create more crossings.

One can also plot the occupation probability distribution for a specific single-particle state in a 2D multi-qubit trap. The spatial profile of the occupation distribution for the single-particle ground state is illustrated in Figure 14. This figure also demonstrates an important property of the multi-qubit trap: namely, that the parameters for the background potential and inter-/intra-cell walls can be tuned to achieve a desirable, in particular, relatively uniform distribution of occupation probability over the entire trap.

6. Controlling the Condensate in the Multi-Qubit Trap: The Gross–Pitaevskii Equation

In real interacting gases, the result for the macroscopic condensate wave function given by the Gross–Pitaevskii Equation (9) significantly deviates from the single-particle ground state of the linear Schr$\ddot{\mathrm{o}}$dinger Equation (8). A difference between the Schr$\ddot{\mathrm{o}}$dinger and Gross–Pitaevskii equations originates due to a collective effect of interparticle interactions described by the nonlinear self-interaction term $g{N}_0{|{\varphi }_0|}^2{\varphi }_0$ in Equation (9). The main features of the condensate are correctly described already in an approximation neglecting the interaction with the noncondensed fraction (the term $2g{n}_{\mathrm{ex}}{\varphi }_0$ in Equation (9)) and assuming $N_0\approx N$. For simplicity’s sake, we adopt the above approximation and limit discussion to 1D and 2D models of the multi-qubit BEC trap.

A 1D model implies a situation when atoms are tightly confined in the transverse to the x axis directions, for example, in a cylinder of length L, ($y,z$)-cross-section area $l_{\perp }^2$ with a small transverse size $l_{\perp }\ll L$, and volume $V=L{l}_{\perp }^2$. Then, in view of the normalization condition ${\int }_V|{\varphi }_0^2|d^3\mathbf{r}=1$ and averaging ${\varphi }_0^2(x,y,z)$ over the small $y,z$ cross-section, the 3D condensate wave function ${\varphi }_0\left(\mathbf{r}\right)$ can be efficiently replaced by a 1D function ${\varphi }_0\left(x\right)/l_{\perp }$ that corresponds to a rescaled interaction parameter $g_1=g/l_{\perp }^2$. Note that the mean-field condition (12) for validity of the 1D model of the Gross–Pitaevskii Equation (9) remains the same as in the usual 3D case,

$$8\pi a\ll {\left(\frac{V}{N}\right)}^{1/3},\frac{\xi }{d}=\frac{1}{\sqrt{8\pi a{(N/V)}^{1/3}}},g=\frac{4\pi {ħ}^{2}a}{m}(3\mathrm{D}\mathrm{mean}-\mathrm{field}\mathrm{regime}),$$

(27)

only if the average distance between atoms is $d={(L{l}_{\perp }^2/N)}^{1/3}$, i.e., the volume density of atoms in the trap is large enough: $N/\left(L{l}_{\perp }^2\right)>1/l_{\perp }^3$. Otherwise, the average distance between atoms becomes equal to $d=L/N$ and the mean-field validity condition (12) imposes a requirement on the scattering length

$$8\pi a\ll N{l}_{\perp }^2/L,\frac{\xi }{d}=\sqrt{\frac{N{l}_{\perp }^2}{8\pi aL}},g_1=\frac{g}{l_{\perp }^2}\mathrm{if}N{l}_{\perp }<L(1\mathrm{D}\mathrm{mean}-\mathrm{field}\mathrm{regime}),$$

(28)

which is getting more stringent with a decreasing number of atoms. It is worth noting that when $aN\gg L$, the system locally retains the original 3D character despite its 1D geometrical appearance, $L\gg l_{\perp }$. Only in the opposite case, when $aN\ll L$, the system approaches the ground state in the transverse directions and enters the so-called 1D mean-field regime (see [35] and references therein). In the low-density limit $8\pi a\gg N{l}_{\perp }^2/L$, which corresponds to the strong-coupling 1D limit $g_1\equiv g/l_{\perp }^2\to \infty$ and is opposite to (28), the mean-field approach fails and the system becomes the so-called Tonks–Girardeau gas of impenetrable bosons.

Similarly, a 2D model implies a situation when atoms are tightly confined just in one, axial direction, say, along the z-axis within a small linear dimension $l_z$, while the cross-section area of the trap of volume $V=L{L}^{\prime }{l}_z$ is relatively large, $L{L}^{\prime }\gg l_z^2$. In this case, the 3D condensate wave function ${\varphi }_0\left(\mathbf{r}\right)$ can be efficiently replaced by a 2D function ${\varphi }_0(x,y)/\sqrt{l_z}$ that corresponds to a rescaled interaction parameter $g_2=g/l_z$. The mean-field condition (12) for validity of the 2D model of the Gross–Pitaevskii Equation (9) retains the usual 3D form (27) only if the average distance between atoms is $d={(L{L}^{\prime }{l}_z/N)}^{1/3}$, i.e., the volume density of atoms in the trap is large enough: $N/\left(L{L}^{\prime }{l}_z\right)\gg 1/l_z^3$. Otherwise, that is when $N{l}_z^2<L{L}^{\prime }$, the average distance between atoms becomes equal to $d=\sqrt{L{L}^{\prime }/N}$ and the mean-field validity condition (12) is reduced to the requirement on the scattering length

$$8\pi a\ll l_z,\frac{\xi }{d}=\sqrt{\frac{l_z}{8\pi a}},g_2=\frac{g}{l_z}\mathrm{if}N{l}_z^2<L{L}^{\prime }(2\mathrm{D}\mathrm{mean}-\mathrm{field}\mathrm{regime}),$$

(29)

which is independent on the number of trapped atoms. Again, at extremely low densities, when $|ln(N{l}_z^2/L{L}^{\prime })|>l_z/a$, the mean-field approach fails, the interaction constant $g_2=g/l_z$ should be replaced by the density-dependent parameter ${\tilde{g}}_2=4\pi {ħ}^2/[m|ln(N{l}_z^2/L{L}^{\prime })|]$, and the system enters the regime analogous to the Tonks–Girardeau 1D regime [35].

6.1. Single-Qubit Trap: 1D Analytical and 2D Numerical Solutions to the Gross–Piraevskii Equation

The solution to the 1D nonlinear Schr$\ddot{\mathrm{o}}$dinger, that is Gross–Piraevskii, Equation (9) in the stated simple model can be found similar to the solution to the 1D linear Schr$\ddot{\mathrm{o}}$dinger Equation (15), described above, if one employs the elliptic Jacobi function $\mathrm{sn}(x|p)$ instead of the exponential function $exp\left(x\right)$. For simplicity’s sake, we adopt the Bogoliubov approximation at very low temperature $T\to 0$ assuming that practically all atoms are condensed, $N_0\approx N$, and the effect of the noncondensed atoms on the condensate is negligible.

Then, the condensate wave function in a box trap with zero potential, $U\left(x\right)=0$, and Dirichlet (zero) boundary conditions is given by the elliptic Jacobi function,

$${\varphi }_0\left(x\right)=\sqrt{\frac{pK\left(p\right)}{K\left(p\right)-E\left(p\right)}}\frac{\mathrm{sn}\left(\left.2K\left(p\right)\frac{x}{L}\right|p\right)}{\sqrt{L}},\frac{L}{\xi }=\sqrt{8K\left(p\right)\left(K\left(p\right)-E\left(p\right)\right)}.$$

(30)

It varies from the half-period sine to an almost constant function (quickly decreasing to zero just in the narrow boundary regions) with the interaction g increasing from zero to the larger values. The characteristic scale of the condensate is determined by the healing length (11). The solution includes complete elliptic integrals of the first and second kinds:

$$K\left(p\right)={\int }_0^{\pi /2}{(1-p{sin}^2\theta )}^{-1/2}d\theta ,E\left(p\right)={\int }_0^{\pi /2}{(1-p{sin}^2\theta )}^{1/2}d\theta ,$$

(31)

According to Equation (30), the range of the parameter p is from 0 to 1. The chemical potential is determined by the normalization condition $\int |{\varphi }_0{\left(x\right)|}^{2}dx=1$ as follows

$$\mu =4(1+p)K^2\left(p\right){ħ}^2/\left(2m{L}^2\right).$$

(32)

The analytical solution in Equation (30) fully describes the effect of the interparticle interaction on the condensate profile in each single-qubit cell if the inter-cell potential walls are infinitely high, the background potentials are the same in both halves of each cell and the intra-cell delta-function potentials are absent. It is illustrated in Figure 15 for the single-qubit cell. When the healing length is much longer than the cell’s length, $\xi \gg L$, that is, the interaction is very weak and the gas is almost ideal, the parameter p is very close to zero. As a result, the condensate profile in each cell is very close to the ground-state solution to the single-particle Schr$\ddot{\mathrm{o}}$dinger Equation (15), that is, a half of the sine function, ${\varphi }_0\left(x\right)=\sqrt{2/L}sin(\pi x/L)$, and $\mu \approx {\pi }^2{ħ}^2/\left(2m{L}^2\right)$. In the opposite case of a very short healing length, $\xi \ll L$, the parameter p approaches 1 and the strong interparticle interaction makes the condensate profile more flat and spread over the entire cell, except for narrow boundary layers of thickness $\xi$ near the walls. Obviously, a similar situation takes place in each half of the single-qubit cells if the intra-cell potential walls are also infinitely high.

The effect of the repulsive interparticle interaction on the condensate profile in the single-qubit box cell in 2D is shown in Figure 16. In the case of an ideal gas, the atoms condense into the ground state ${\psi }_0(x,y)=\frac{2}{L}sin\frac{\pi x}{L}\times sin\frac{\pi y}{L}$ of the single-particle Schr$\ddot{\mathrm{o}}$dinger Equation (15), as shown in Figure 16a. In the case of an interacting gas, the condensate profile ${\varphi }_0(x,y)$ is given by the numerical solution to the 2D Gross–Pitaevskii Equation (9), as shown in Figure 16b. Comparison of the two plots clearly shows that the particle repulsion flattens the peak of the ground-state wave function and forces the condensate to spread over the entire single-qubit cell. Just the boundary layers of a healing-length thickness remain unoccupied by the condensate.

These effects can be approximately described analytically by a product of the exact analytical solution (30) to the Gross–Pitaevskii equation in the 1D box trap along the x-axis and the similar solution along the y-axis,

$${\varphi }_0(x,y)\approx \frac{pK\left(p\right)}{\left[K\right(p)-E(p\left)\right]L}\mathrm{sn}\left(\left.2K\left(p\right)\frac{x}{L}\right|p\right)\mathrm{sn}\left(\left.2K\left(p\right)\frac{y}{L}\right|p\right).$$

(33)

Such a 1D-factorization approximation is shown in Figure 16c. It takes into account the interparticle interaction only partially via its separate manifestations along the x and y dimensions. Comparing Figure 16b and Figure 16c, we conclude that the above approximation slightly underestimates the effect of 2D nonlinear diffusion of the condensate due to the self-interaction $g{N}_0{|{\varphi }_0(x,y)|}^2{\varphi }_0(x,y)$, which is opposite to the phenomenon of self-focusing of an intensive laser light beam in a nonlinear medium. Nevertheless, the 1D-factorization approximation represents the effect of the interparticle interaction on the condensate profile in a box trap in a qualitatively correct fashion.

6.2. Condensate Wave Function vs. Single-Particle Ground State in a Multi-Qubit Trap

For a nontrivial multi-qubit BEC trap, due to the presence of the trapping potential $U\left(x\right)\ne 0$, as shown in Equation (13), the Gross–Pitaevskii Equation (9), that is, the nonlinear Schr$\ddot{\mathrm{o}}$dinger equation, needs to be solved numerically, for instance, by the method of an imaginary-time evolution (see, e.g., [87]). It is illustrated in Figure 17 and Figure 18 for the case of a four-qubit 1D and 2D trap, respectively.

As a result of the interparticle repulsion, the condensate tends to spread more uniformly over all single-qubit cells. This tendency works against condensate fragmentation [78] and in favor of the formation of a common condensate occupying the entire BEC trap. Moreover, with increasing interaction, the regions of low condensate occupation near the inter- and intra-cell potential walls begin shrinking as well. Both of the above effects significantly increase the number of Bogoliubov-coupled excited states and magnitude of their Bogoliubov couplings in Equation (10) that favors manifestation of the ♯P-hardness of the atomic boson sampling as is explained in Section 2.

It is worth noting that such a considerable expansion of the condensate shown in Figure 17 and Figure 18 is provided by means of the interparticle interaction alone, without employment of the background potential, which also allows one to control the condensate profile in a similar direction via restructuring the ground-state wave function as is shown in Figure 8 and Figure 14.

Moreover, if the background potential makes the trap asymmetrical, the increasing repulsive interaction tends to restore the trap’s symmetry by converting an asymmetrical single-particle ground state into a more symmetrical condensate wave function. Such evolution of the condensate in the asymmetrical trap is illustrated in Figure 17b and should be compared against the condensate evolution in the symmetrical trap shown in Figure 17a. Clearly, a strong interparticle interaction makes the condensate profiles in both traps almost indistinguishable, while the profiles of the ground state in these traps in the absence of interaction are very different.

7. Controlling Multimode Squeezing of Bogoliubov Transform via Bogoliubov Couplings

In equilibrium, the statistics of the many-body system of atoms in the BEC trap is determined by the independent fluctuations of quasiparticles which form the eigenstates of the Bogoliubov Hamiltonian with the eigenenergies $\left\{E_j\right\}$ and have the Bose–Einstein occupation number statistics with an average occupation number ${\overline{n}}_j={\left(e^{(E_j-\mu )/T}-1\right)}^{-1}$. The two-component quasiparticle wave function $\{u_j,v_j\}$ determine the excited-particle field operator (3) and obeys the Bogoliubov–de Gennes equations:

$$\begin{array}{cc}\hfill & \widehat{L}{u}_j+g{N}_0{\varphi }_0^2\left(\mathbf{r}\right)v_j=+E_j{u}_j,\hfill \\ \hfill & \widehat{L}{v}_j+g{N}_0{\left({\varphi }_0^{*}\right)}^2\left(\mathbf{r}\right)u_j=-E_j{v}_j,\hfill \end{array}$$

(34)

where

$$\widehat{L}\equiv -\frac{{ħ}^2\Delta }{2m}+U\left(\mathbf{r}\right)+2g\left(N_0{|{\varphi }_0\left(\mathbf{r}\right)|}^2+n_{\mathrm{ex}}\left(\mathbf{r}\right)\right)-\mu ;n_{\mathrm{ex}}=\sum _j\left[|v_j{\left(\mathbf{r}\right)|}^2+\frac{|u_j{\left(\mathbf{r}\right)|}^2+{|v_j\left(\mathbf{r}\right)|}^2}{exp(E_j/T)-1}\right].$$

In essence, the Bogoliubov–de Gennes equations express the fact of diagonalization of the Hamiltonian (4) by the Bogoliubov transformation to quasiparticle creation/annihilation operators, as stated in the matrix form in Equation (6), in the form of differential equations for the coefficients $\{u_j,v_j\}$ in the expansion of the field operator (3) via the quasiparticle operators. The wave functions are normalized to unity: ${\int }_V{|{\varphi }_0|}^2{d}^3\mathbf{r}=1$, ${\int }_V\left(|u_j{|}^2-{|v_j|}^2\right)d^3\mathbf{r}=1$; $j=1,2,\dots$. For simplicity’s sake, hereinafter, we assume that all wave functions ${\varphi }_0,u_j,v_j$ are real valued. Below, we again neglect by temperature-dependent, Popov’s corrections, that is, skip the contribution due to the noncondensate density $n_{\mathrm{ex}}$ and assume $N_0\approx N$.

The matrix R of the Bogoliubov transformation (2) can be found from the equation

$$\left(\begin{array}{c}{\widehat{\mathbf{a}}}^{†}\\ \widehat{\mathbf{a}}\end{array}\right)=R\left(\begin{array}{c}{\widehat{\mathbf{b}}}^{†}\\ \widehat{\mathbf{b}}\end{array}\right)=\left[\begin{array}{cc}A& B^{*}\\ B& A^{*}\end{array}\right]\left(\begin{array}{c}{\widehat{\mathbf{b}}}^{†}\\ \widehat{\mathbf{b}}\end{array}\right)$$

(35)

that relates creation and annihilation operators ${\widehat{a}}^{†},\widehat{a}$ describing bare particles to the operators ${\widehat{b}}^{†},\widehat{b}$ describing quasiparticles as per Equation (3). Since the Bogoliubov transformation (2) leaves the canonical Bose commutation relations invariant, it obeys the symplectic property

$$R\mathrm{\Omega }{R}^T=\mathrm{\Omega },\mathrm{\Omega }=\left[\begin{array}{cc}0& \mathbb{1}\\ -\mathbb{1}& 0\end{array}\right].$$

(36)

Another, equivalent to (35), representation of the Bogoliubov transformation (2) can be written in terms of the wave functions, rather than the operators, determining the particle field operator in Equation (3):

$$\left(\begin{array}{c}\mathbf{\varphi }\\ 0\end{array}\right)=R\left(\begin{array}{c}\mathbf{u}\\ -{\mathbf{v}}^{*}\end{array}\right)\equiv \left[\begin{array}{cc}A& B^{*}\\ B& A^{*}\end{array}\right]\left(\begin{array}{c}\mathbf{u}\\ -{\mathbf{v}}^{*}\end{array}\right),$$

(37)

$$R^T\left(\begin{array}{c}0\\ \mathbf{\varphi }\end{array}\right)\equiv \left[\begin{array}{cc}{A}^T& B^T\\ B^{†}& A^{†}\end{array}\right]\left(\begin{array}{c}0\\ \mathbf{\varphi }\end{array}\right)=\left(\begin{array}{c}{\mathbf{v}}^{*}\\ \mathbf{u}\end{array}\right).$$

(38)

The column-vectors $\mathbf{\varphi }={({\varphi }_1,{\varphi }_2,\dots )}^T$ and $\mathbf{u}={(u_1,u_2,\dots )}^T$, $\mathbf{v}={(v_1,v_2,\dots )}^T$ are composed of the excited states $\left\{{\varphi }_k\right\}$ and quasiparticle wave functions $\left\{u_j\right\}$, $\left\{v_j\right\}$, respectively.

Projecting Equation (38) onto a set of the orthonormal excited states $\{{\varphi }_k|k=1,2,\dots \}$ which are also orthogonal to the condensate wave function ${\varphi }_0$, we obtain the explicit formulae for the entries of the Bogoliubov block matrices $A=\left(A_{kj}\right)$ and $B=\left(B_{kj}\right)$,

$$A_{kj}=\int u_j^{*}\left(\mathbf{r}\right){\varphi }_k\left(\mathbf{r}\right)d^3\mathbf{r},B_{kj}=\int v_j^{*}\left(\mathbf{r}\right){\varphi }_k^{*}\left(\mathbf{r}\right)d^3\mathbf{r},$$

(39)

as overlapping integrals between those bare-particle wave functions and the quasiparticle wave functions given by the solution to the Bogoliubov–De Gennes equations (34).

The Bogoliubov matrix R can be expressed explicitly also via the Bogoliubov couplings in Equation (10) by means of a pure algebraic diagonalization of the Bogoliubov Hamiltonian in the sense of the matrix Equation (6). Indeed, in any basis $\{{\varphi }_k|k=1,2,\dots \}$ of excited states, orthogonal to the condensate wave function and constituting the excitation field operator ${\widehat{\psi }}_{\mathrm{ex}}={\sum }_{k\ne 0}{\varphi }_k{\widehat{a}}_k$ as in Equation (3), the blocks of the Hamiltonian matrix in Equation (5) are explicitly given by the Bogoliubov couplings in Equation (10) as follows

$$K=\left({\epsilon }_{k{k}^{\prime }}-\mu {\delta }_{k,k^{\prime }}+2{\Delta }_{k{k}^{\prime }}\right),\tilde{K}=\frac{1}{2}\left({\tilde{\Delta }}_{k{k}^{\prime }}\right),H=\left[\begin{array}{cc}\tilde{K}& K\\ K^{*}& {\tilde{K}}^{*}\end{array}\right].$$

(40)

Here, ${\epsilon }_{k{k}^{\prime }}=\langle {\varphi }_k|\widehat{\epsilon }|{\varphi }_{k^{\prime }}\rangle$ is the matrix of the single-particle energy operator $\widehat{\epsilon }=-{ħ}^2\Delta /\left(2m\right)+U\left(\mathbf{r}\right)$ which constitutes the single-particle Schr$\ddot{\mathrm{o}}$dinger Equation (8). In particular, the basis $\{{\varphi }_k|k=1,2,\dots \}$ can be constructed out of the excited-state eigenfunctions $\{{\psi }_n|n=1,2,\dots \}$ of the Schr$\ddot{\mathrm{o}}$dinger Equation (8) by means of the standard Gram–Schmidt orthonormalization starting from making these functions orthogonal to the condensate wave function ${\varphi }_0$. Then, by means of the symplectic property (36), Equation (6) determining the Bogoliubov transformation can be rewritten as the following equation

$$\mathrm{\Omega }HR=R\left[\begin{array}{cc}E& 0\\ 0& -E\end{array}\right],E=\mathrm{diag}\{E_j|j=1,2,\dots \}.$$

(41)

It states that the j-th column of the Bogoliubov matrix, ${\mathbf{R}}_{\mathbf{j}}={\{A_{1j},A_{2j},\dots ,B_{1j},B_{2j},\dots \}}^T$, is the eigenvector of the matrix $\mathrm{\Omega }H$ corresponding to the quasiparticle eigenenergy $E_j$, that is

$$\mathrm{\Omega }H{\mathbf{R}}_{\mathbf{j}}=E_j{\mathbf{R}}_{\mathbf{j}},\mathrm{\Omega }H=\left[\begin{array}{cc}{K}^{*}& {\tilde{K}}^{*}\\ -\tilde{K}& -K\end{array}\right].$$

(42)

(There is also the nonphysical eigenvector counterpart ${\mathbf{R}}_{\mathbf{j}}^{(-)}={\{B_{1j}^{*},B_{2j}^{*},\dots ,A_{1j}^{*},A_{2j}^{*},\dots \}}^T$ corresponding to the negative eigenenergy $-E_j<0$.)

After calculating the Bogoliubov transformation matrix R as per Equation (6), one can find the multimode squeezing parameters [1,68,69,70,71,72,73] from Equation (2). Quantum statistics of the many-body fluctuations in a BEC trap and, in particular, the computational complexity of the atomic boson sampling are determined by two fundamental sets of eigenvectors and eigenvalues associated with the diagonalization of (a) the squeezing matrix and (b) the Bogoliubov Hamiltonian, as is explained in Section 1 and Section 2. Both of those sets of eigenvectors and eigenvalues are determined by the Bogoliubov couplings (10) via the Bogoliubov transformation matrix R. Thus, the key problem is to calculate the Bogoliubov couplings and understand how many of them can be essentially nonzero and controllable in a wide range within the multi-qubit BEC trap suggested and described in this paper.

Knowing the condensate wave function from the solution to the Gross–Pitaevskii equation outlined in Section 6 and choosing the bare-particle excited states, for example, as the purely harmonic, sine functions or the solutions to the single-particle Schr$\ddot{\mathrm{o}}$dinger Equation (8) (see Section 4 and Section 5), made orthogonal to the condensate and each other via the standard Gram–Schmidt orthonormalization, it is straightforward to calculate the integrals constituting the Bogoliubov couplings (10) and analyze their set for the multi-qubit BEC trap. The related numerical results are illustrated for the cases of symmetrical and asymmetrical 1D four-qubit traps in the plots shown in Figure 19 and Figure 20, respectively.

First, comparing Figure 20 against Figure 19 makes it clear that the trap asymmetry greatly enlarges the number of Bogoliubov-coupled bare-particle excited states. Indeed, in the asymmetrical trap, the essentially nonzero couplings spread much further from the main diagonal of the Bogoliubov-coupling matrix ${\Delta }_{k,k^{\prime }}$ than in the symmetrical trap where only narrow lanes of entries around the main diagonal and anti diagonal are essentially nonzero. In addition, the degeneracy of zero coupling between the bare-particle excited states of exactly odd and even spatial parity in the symmetrical trap (Figure 19), that results in exactly zero values of all entries in each diagonal of an odd number parallel to the main diagonal, is essentially broken in the asymmetrical trap (Figure 20). It is restored only in the limit of a very strong interparticle interaction.

Second, the maximum spread of essentially nonzero couplings occurs at a moderate interparticle interaction $\frac{L}{\xi }\sim 1$ corresponding to the healing length $\xi$ being on order of the single-qubit cell length L. Much stronger interaction, $\frac{L}{\xi }\gg 1$, tends to localize nonzero couplings just onto the main and anti diagonals. Both these effects are seen in each row of plots in Figure 19 and Figure 20 where the interaction strength is increasing from left to right.

Third, changing the bare-particle excited states, chosen for the simultaneous atom-number detecting within the atomic boson sampling, from the set generated by the Gram–Schmidt orthonormalization out of the sine functions (the upper rows in Figure 19 and Figure 20) to the set generated out of the solutions to the single-particle Schr$\ddot{\mathrm{o}}$dinger equation (the lower rows in Figure 19 and Figure 20) greatly affects the structure of the Bogoliubov-coupling matrix both in the symmetrical and asymmetrical traps.

All of the above observations confirm that the inference the multi-qubit BEC trap provides is an excellent opportunity for controlling the Bogoliubov couplings and, hence, the multi-mode squeezing and interference of bare-atom excited modes in a very wide range. Obviously, the more chaotic, messy, dense and wide the distribution of the essentially nonzero elements over the Bogoliubov-coupling matrix (10), the more favorable the set of trap’s parameters and bare-atom excited states chosen for detection of atom numbers for testing manifestations of the computational ♯P-hardness of atomic boson sampling. Among patterns shown in Figure 19 and Figure 20, the one in the center of the lower row in Figure 20 is the most representative picture of such a complexity.

The asymptotic parameter of this complexity is determined by the Bogoliubov transformation via a multimode dimensionality of the subspace of the excited-states involved in the squeezing-matrix eigenvectors with essentially nonzero squeezing parameters (see Equation (2)) and the Hamiltonian-matrix eigenvectors with low enough eigenenergies corresponding to quasiparticles with essentially nonzero populations (see Equation (42)). In general, this asymptotic parameter increases as the number of groups of excited states chosen for occupation sampling via multi-detector imaging is growing. However, for a given experimental setup with a BEC trap of a finite size, there is a maximum number M of modes/channels started from which a further increase of the number of sampled/detected occupations would not essentially increase the complexity of boson sampling.

8. Toward Experiments on Atomic Boson Sampling in a BEC Trap

Suppose one has an appropriate BEC trap. (A possible model/example of such a trap is discussed in the previous sections.) Then, as is explained in Section 1.1 and Section 1.2, the excited atoms, by themselves, naturally fluctuate and stay in the squeezed states inside the trap even at thermal equilibrium due to interactions with each other. This allows one to eliminate any nonequilibrium processes or dynamics, such as a precise time-dependent control of system parameters and gates or any other type of processing usually associated with quantum computers or simulators, as well as the sophisticated external sources of squeezed or single bosons (required for photonic sampling) from the atomic sampling experiments. It remains just to split the noncondensate into fractions based on the groups of excited states and to measure the distribution of atom numbers over the chosen groups of excited states by means of appropriate detectors.

For instance, one can divide the volume of the trap into a system of spatial cells. Another possibility is to separate atoms in accord with their velocities, that is, to deal with the cells in the momentum space. Anyway, the measurement of atom numbers could be completed by means of a multi-detector imaging. In a BEC destruction scheme, one switches off the confining trap and allows the cloud of trapped atoms to expand freely. In this case, following a standard time-of-flight measuring technique, it is required to take a few successive images of the expanding cloud and properly interpret them in terms of kinetic equations for expansion. In this way, different spatial or momentum subsets of atoms could be separated from each other, and sampling of their occupation numbers could be obtained.

The imaging technique implies an illumination of the atomic cloud with a laser pulse and measuring its transmitted or scattered components by multiple detectors. The transmitted signal carries information on the absorption, dispersion and polarization transformation of light caused by an atomic cloud [33,75,88,89,90]. The signal due to scattering and fluorescence [91] could be controlled and structured by employing special external cavities and laser sources that support light modes which mimic the excited states preselected for sampling. The optical imaging for atomic boson sampling has much in common with the experiments on the local atom-number fluctuations in BEC traps [47,48,88,90,91,92,93,94].

The spatial or momentum cells/modes represent groups of excited states selected for detection/sampling. Of course, the excited states can be described/composed via an arbitrary basis in the single-particle Hilbert space. Accordingly, the analytical formulae for their joint occupation probability distribution $\rho \left(\left\{n_k\right\}\right)$ and characteristic function, derived in [1], have a universal form, i.e., are valid for any choice of such a basis. A transition from one basis to another one just adds an extra unitary transformation (7) to the Bogoliubov matrix R in Equation (2). Moreover, the universality of the general result for the characteristic function obtained in [1] extends to the so-called marginal or coarse-grained statistics of occupations of any groups of excited states, that is, to the occupation statistics evaluated irrespective to the occupations of all other excited states. The corresponding “incomplete” experiments on atomic boson sampling are the ones to be devised and implemented in reality. Obviously, the condensed atoms, which constitute the macroscopic condensate wave function ${\varphi }_0\left(\mathbf{r}\right)$ orthogonal to the excited states ${\varphi }_k\left(\mathbf{r}\right)$, should not be countered during the sampling procedure.

A computational complexity of atomic boson sampling depends on the number M of groups of excited atomic states which are resolved by the multi-detector imaging and are subject to interference due to mixing through the quasiparticles and to squeezing due to interparticle interaction. This number M plays a part of the number of channels in the optical interferometer. A mean occupation of the groups of excited states scales as $(N-N_0)/M$. Obviously, by increasing the total number N of atoms loaded in the trap and, therefore, the number $N-N_0$ of atoms in the noncondensate, one can make larger and, hence, easier for detection the number of atoms in each of M groups of preselected excited states. Note, however, that the asymptotic parameter responsible for the ♯P-hardness of atomic boson sampling is not proportional to any of the numbers $N,N-N_0$ or M.

Modern technology allows one to measure the number of atoms in a specified volume or subset of atoms with nearly single atom resolution [91,93,94,95]. Yet, achieving the single atom accuracy is not absolutely necessary. In particular, the ♯P-hardness of boson sampling exists even in the case of the threshold detecting scheme in which an outcome of the measurement is just zero or non-zero occupation in each preselected group of excited states [16,24,25].

Finally, an experimental setup should provide the means to reconfigure detectors for projecting upon a vastly varying set of groups of excited states, i.e., to accumulate statistics of joint occupations for numerous different subsets of groups of excited states. Only in this way can the quantum advantage be demonstrated at the most challenging level of the average case.

9. Concluding Remarks

We introduce the multi-qubit BEC trap for studying manifestations of the quantum many-body statistical phenomena which are ♯P-hard for computing. In particular, we describe the basic properties of the multi-qubit trap, including the single-particle excited states and their energy spectrum via the single-particle Schr$\ddot{\mathrm{o}}$dinger Equation (8), the condensate wave function versus the single-particle ground state via the Gross–Pitaevskii Equation (9), and the Bogoliubov couplings (10) between excited states responsible for the formation of quasiparticles and multimode squeezing via the Bogoliubov–de Gennes equations (34). It is completed within the 1D and 2D models, as shown in Equations (13) and (14).

We show that the multi-qubit BEC trap offers a convenient and thoroughgoing control of the many-body system parameters essential for the interplay between excited states’ interference and squeezing. This interplay can be revealed via an apropriate decomposition of the Bogoliubov-transformation matrix in Equation (2) and is responsible for the computational ♯P-hardness which is the basis for a potential quantum advantage of atomic boson sampling over classical computing [1].

It would be very interesting to study experimentally various phenomena associated with the atomic boson sampling. The BEC trap is a boson-sampling platform alternative to a photonic interferometer. Both systems provide the output multivariate statistics which shows computational ♯P-hardness associated with the hafnian of complex-valued matrices. The proposed multi-qubit trap design discussed in the present paper allows one to vary those matrices and, hence, the output statistics over a wide range. Thus, the latter, major requirement for testing quantum advantage is fulfilled by the multi-qubit BEC trap. The remarkable fact is that classical computing of the hafnian of even relatively low-dimensional matrices corresponding to the number of sampled modes/channels of the order of $M=8\times 8=64$ is already inaccessible to modern supercomputers. Especially promising are boson-sampling experiments with the multi-qubit BEC trap containing a finite number M of lower-miniband split-off excited states or groups of them (see Figure 1).

The case of a few single-qubit cells with a relatively small number of sampled occupations $M=2,3,4,\dots$ promises the discovery of new quantum effects similar and beyond a particle analog of the simple Hong-Ou-Mandel interference effect. It can be accomplished by means of the current magneto-optical trapping and detection technology. The value of such experiments for the comprehension of the fundamental aspects of the many-body quantum systems responsible for their computational ♯P-hardness is difficult to overestimate.

The conclusive experiments with an asymptotically large numbers of single-qubit cells, $Q\ll 1$, and sampled excited states or groups of them, $M\gg 1$, addressing the computational ♯P-hardness of quantum many-body processes are very challenging. Yet, they seem to be within reach and could hit convincing manifestations of quantum advantage.