Bogoliubov theory for atom scattering into separate regions

Our theoretical description of the pair-production process starts with a many-body Hamiltonian with contact two-body interactions

$$\begin{equation} \skew2\hat{H}=\int\!\! \mathrm{d}{\mathbf r}\,\hat\Psi^\dagger({\mathbf r})\left(-\frac{\hbar^2\nabla^2}{2m}+V({\mathbf r})\right)\hat\Psi({\mathbf r})+\frac g2\int\!\! \mathrm{d}{\mathbf r}\,\hat\Psi^\dagger({\mathbf r})\hat\Psi^\dagger({\mathbf r})\hat\Psi({\mathbf r})\hat\Psi({\mathbf r}). \end{equation} \tag{ 1 }$$

Here V (r) is an external trapping potential and $g=\frac {4\pi \,\hbar ^2a}m$ is the strength of the two-body interactions, a is the scattering length, m is the atomic mass and the field operator $\hat \Psi ({\mathbf r})$ satisfies the bosonic commutation relations. To derive the Bogoliubov equation, we first find the c-number (mean field) wave function of the BEC using the Gross–Pitaevskii equation (GPE)

$$\begin{equation} \mathrm{i}\,\hbar\partial_t\psi({\mathbf r})=\left(-\frac{\hbar^2\nabla^2}{2m}+V({\mathbf r})+g|\psi({\mathbf r})|^2\right)\psi({\mathbf r}). \end{equation} \tag{ 2 }$$

We then write the field operator as a sum of the c-number part and the Bogoliubov correction, $\hat \Psi ({\mathbf r})=\psi ({\mathbf r})+\hat \delta ({\mathbf r})$ and insert this expression into (1). By keeping only the terms up to quadratic in $\hat \delta$ we obtain the Bogoliubov Hamiltonian

$$\begin{equation} \fl \skew2\hat{H}_{\rm bog}=\int\!\! \mathrm{d}{\mathbf r}\,\hat\delta^\dagger({\mathbf r})\left(-\frac{\hbar^2\nabla^2}{2m}+V({\mathbf r})+2g|\psi({\mathbf r})|^2\right)\hat\delta({\mathbf r}) +g\int\!\! \mathrm{d}{\mathbf r}\left(\hat\delta^\dagger({\mathbf r})\hat\delta^\dagger({\mathbf r})\,\psi^2({\mathbf r})+{\rm h.c.}\right). \end{equation} \tag{ 3 }$$

The resulting Bogoliubov equation of motion is linear

$$\begin{equation} \mathrm{i}\,\hbar \partial_t \hat \delta({\mathbf r},t)=\left(-\frac{\hbar^2 \nabla^2}{2m} + V({\mathbf r}) + 2\,g |\psi({\mathbf r})|^2 \right) \hat \delta({\mathbf r},t)+g \psi^2({\mathbf r}) \hat \delta^\dagger({\mathbf r},t). \end{equation} \tag{ 4 }$$

Usually, a numerical solution of this equation is found in a following way. The field operator is expanded in a basis of wave-functions $\tilde \varphi _i({\mathbf r})$ which match the geometry of the scattering problem

$$\begin{equation} \hat\delta({\mathbf r},t)=\sum_i\tilde\varphi_i({\mathbf r})\hat a_i(t). \end{equation} \tag{ 5 }$$

This expression is inserted into equation (4), the resulting equation is multiplied by $\tilde \varphi ^*_j({\mathbf r})$ and the outcome is integrated by sides over the whole space. In effect, what we obtain is an equation of motion, which, through the matrices $\hat {\mathcal {A}}$ and $\hat {\mathcal {B}}$, couples the evolution of the jth operator $\hat a_j(t)$, with (in general) all others operators

$$\begin{equation} \mathrm{i}\,\partial_t\hat a_j(t)=\sum_k\mathcal{A}_{jk}\hat a_k(t) + \sum_k\mathcal{B}_{jk}\hat a_k^{\dagger}(t). \end{equation} \tag{ 6 }$$

This equation is linear—a consequence of the linearity of the Bogoliubov equation (4)—so the general solution of (6) reads

$$\begin{equation} \hat a_i(t)=\sum_j\mathcal{C}_{ij}(t)\hat a_j(0)+\sum_j\mathcal{S}_{ij}(t)\hat a_j^\dagger(0), \end{equation} \tag{ 7 }$$

where the matrices $\hat {\mathcal {C}}$ and $\hat {\mathcal {S}}$ satisfy $\hat {\mathcal {C}}\hat {\mathcal {C}}^\dagger -\hat {\mathcal {S}}\hat {\mathcal {S}}^\dagger =\hat {\mathbbm{1}}$ and $\hat {\mathcal {C}}\hat {\mathcal {S}}^{\mathrm { T}}-\hat {\mathcal {S}}\hat {\mathcal {C}}^{\mathrm { T}}=0$. Later, we will apply this method to solve the Bogoliubov dynamics of the twin-beam production. However, we will show in the following, that in cases where the detailed form of the Hamiltonian (3) drives the scattering of atomic pairs into opposite regions (as indeed happens in twin-beam experiments), the basic properties of the system can be deduced analytically if an appropriate set of mode functions φ_i(r) is chosen.

Let us denote the two separate regions into which the particles are scattered by L (left) and R (right). Particles populate L and R in a process of elastic scattering, so the regions are usually separated in momentum space. From this point of view, it is convenient to switch to the space of wave-vectors k and decompose the field operator as follows:

$$\begin{equation} \hat\delta({\mathbf k},t)=\sum_i\varphi_{\rm R}^{(i)}\!({\mathbf k},t)\,\hat a^{(i)}_{\rm R}(t)+\sum_i\varphi_{\rm L}^{(i)}\!({\mathbf k},t)\,\hat a^{(i)}_{\rm L}(t). \end{equation} \tag{ 8 }$$

The operators $\hat a^{(i)}_{\mathrm { R/L}}(t)$ annihilate a particle in a mode characterized by the time-dependent wave function φ⁽ⁱ⁾_R/L(k,t), which is localized in the right/left region in momentum space. We underline, that this kind of separation is also present in position space after expansion of the cloud, or after application of a Stern–Gerlach pulse in internal-state experiments, respectively. Moreover, the vector k might denote the quasi-momentum, if the scattering takes place in an optical lattice.

Formally, the only difference between the formulation (5) and (8) is the splitting of the field operator into the R and L modes. However, for a linear equation of motion such as (4), there exists a unique basis of mode functions for which the evolution equations of the mode pairs decouple from each other:

$$\begin{equation} \hat a^{(i)}_{\rm R}(t)=c_i(t)\,\hat a^{(i)}_{\rm R}(0)+s_i(t)\,\hat a^{(i)\dagger}_{\rm L}(0),\end{equation} \tag{ 9 }$$

$$\begin{equation} \hat a^{(i)}_{\rm L}(t)=c_i(t)\,\hat a^{(i)}_{\rm L}(0)+s_i(t)\,\hat a^{(i)\dagger}_{\rm R}(0), \end{equation} \tag{ 10 }$$

where |c_i(t)|² − |s_i(t)|² = 1. This form of the Bogoliubov equation has a simple physical interpretation: atoms scatter pair-wise into opposite regions, and the total field operator (8) is a sum of independent mode pairs, which are squeezed in their relative population fluctuations, as will be explained in detail below.

Although the diagonal form (9) and (10) is much clearer than (7), it is not obvious at the moment how this particular basis (8) can be found. This is done in two steps, applying the procedure of the Bloch–Messiah reduction [33–35]. First, using equations (8)–(10), we evaluate the one-body density matrix (first-order correlation function) and obtain

$$\begin{eqnarray} G^{(1)}({\mathbf k}_1,{\mathbf k}_2;t)&\equiv&\langle\hat\delta^\dagger({\mathbf k}_1,t)\hat\delta({\mathbf k}_2,t)\rangle\\ &=&\sum_in_i\left(\varphi_{\rm R}^{(i)*}\!({\mathbf k}_1,t)\varphi_{\rm R}^{(i)}\!({\mathbf k}_2,t)+\varphi_{\rm L}^{(i)*}\!({\mathbf k}_1,t)\varphi_{\rm L}^{(i)}\!({\mathbf k}_2,t)\right),\nonumber \end{eqnarray} \tag{ 11 }$$

where n_i = |s_i(t)|². Note that n_i is the average occupation of the ith eigen-mode, and that a pair of modes φ⁽ⁱ⁾_R(k,t) and φ⁽ⁱ⁾_L(k,t) is degenerate (has the same eigen-value n_i) due to the assumed symmetry between the left and the right region. Since we are using the Heisenberg picture, the average value in equation (11) and all equations that follow are calculated in the initial vacuum state of scattered atoms.

In any practical approach, if the basis (8) is not known a priori, this step first requires a numerical evaluation of the density matrix (11) in any convenient basis (5), and subsequent diagonalization. Once this is done, then according to equation (11), the basis functions φ⁽ⁱ⁾_R/L(k,t) are the momentary eigen-functions of the one-body density matrix (natural orbitals). However, a second step is necessary to fully determine the functions φ⁽ⁱ⁾_R/L(k,t), because the density matrix—contrary to the field operator (8)—is insensitive to global phases of the mode functions. To retrieve this additional information, we calculate the anomalous density

$$\begin{eqnarray} &&\fl M({\mathbf k}_1,{\mathbf k}_2;t)\equiv\langle\hat\delta({\mathbf k}_1,t)\hat\delta({\mathbf k}_2,t)\rangle=\sum_i\sqrt{n_i(n_i+1)}\left(\varphi_{\rm R}^{(i)}\!({\mathbf k}_1,t)\varphi_{\rm L}^{(i)}\!({\mathbf k}_2,t)+\varphi_{\rm R}^{(i)}\!({\mathbf k}_2,t)\varphi_{\rm L}^{(i)}\!({\mathbf k}_1,t)\right)\nonumber\\ \end{eqnarray} \tag{ 12 }$$

multiply it by sides with the eigen-functions of the density matrix and integrate over space. As a result, we retrieve the information about the phases and obtain the full form of the mode functions φ⁽ⁱ⁾_R/L(k,t) of the diagonal basis.

To summarize, we have outlined the structure of the solution of the Bogoliubov equation for cases where atoms are scattered into two opposite regions. We will now show that the extra step, which is the transition from the 'numerical approach' (5) to the diagonal basis (8), allows to easily determine the basic properties of the system of scattered atoms, like its density or higher correlation functions.

The simplest observable characterizing the pair-production process is the density

$$\begin{equation} \rho({\mathbf k};t)\equiv G^{(1)}({\mathbf k},{\mathbf k};t)=\sum_in_i\left(|\varphi_{\rm R}^{(i)}\!({\mathbf k},t)|^2+|\varphi_{\rm L}^{(i)}\!({\mathbf k},t)|^2\right) \end{equation} \tag{ 13 }$$

which is, consistently with our derivation, localized in the two opposite regions. By integrating the above function over space, we obtain the information about the expected number of scattered atoms as a function of time

$$\begin{equation} \langle\hat N\rangle=\int\!{\mathrm d}{\mathbf k}\,\rho({\mathbf k};t)=2\sum_in_i. \end{equation} \tag{ 14 }$$

Additional information about the system is carried by the correlations between the scattered particles. The probability of simultaneous detection of two atoms at momenta k₁ and k₂ can be obtained from the normalized second-order correlation function

$$\begin{equation} g^{(2)}({\mathbf k}_1,{\mathbf k}_2;t)=\frac{\langle\hat\delta^\dagger({\mathbf k}_1,t)\hat\delta^\dagger({\mathbf k}_2,t)\hat\delta({\mathbf k}_2,t)\hat\delta({\mathbf k}_1,t)\rangle}{\rho({\mathbf k}_1;t)\rho({\mathbf k}_2;t)}. \end{equation} \tag{ 15 }$$

According to the Wick's theorem, this function can be written in terms of the one-body density matrix (11) and the anomalous density (12) as follows:

$$\begin{equation} g^{(2)}({\mathbf k}_1,{\mathbf k}_2;t)=1+\frac{|G^{(1)}({\mathbf k}_1,{\mathbf k}_2;t)|^2+|M({\mathbf k}_1,{\mathbf k}_2;t)|^2}{\rho({\mathbf k}_1;t)\rho({\mathbf k}_2;t)}. \end{equation} \tag{ 16 }$$

The transition from equations (15) to (16) might seem an unnecessary complication, however we will argue that it allows for a simple and intuitive interpretation of the second-order correlation function. According to equation (11), the density matrix is non-vanishing only when k₁ and k₂ are both either in the right or left region, so |G⁽¹⁾|² governs the Hanbury–Brown and Twiss (HBT) type of local correlations. On the other hand, as can be seen from equation (12), the anomalous density is non-zero only when k₁ and k₂ are in the opposite regions, so it describes the cross-correlations between the two members of the scattered pair. Clearly, this simple interpretation of the second order correlation function as a sum of local- and opposite-momentum correlations would have been much more difficult if we had not applied the diagonalization procedure and the Wick's theorem.

Another property characterizing the scattering process are the fluctuations of the population imbalance between the two regions. If these fluctuations are suppressed below the properly defined shot-noise level, the system is number squeezed, which proves that atoms scatter in pairs rather then independently to the left and to the right region. A quantitative description of the number squeezing involves the left and right atom number operators defined as the integrals of the density operators over the corresponding volumes, i.e.

$$\begin{equation} \hat N_{\rm R/L}=\int\nolimits_{\rm R/L}\!\!{\mathrm d}{\mathbf k}\,\hat\delta^\dagger({\mathbf k},t)\hat\delta({\mathbf k},t). \end{equation} \tag{ 17 }$$

The population imbalance operator is then simply defined as $\hat n=\hat N_{\mathrm { R}}-\hat N_{\mathrm { L}}$ and using equation (8) we obtain

$$\begin{equation} \hat n=\sum_i\left(\hat a_{\rm R}^{(i)\dagger}(0)\,\hat a_{\rm R}^{(i)}(0)-\hat a_{\rm L}^{(i)\dagger}(0)\,\hat a_{\rm L}^{(i)}(0)\right). \end{equation} \tag{ 18 }$$

The number squeezing factor is defined as

$$\begin{equation} \xi^2=\frac{\Delta^2\hat n}{\langle\hat N\rangle}, \end{equation} \tag{ 19 }$$

where $\Delta ^2\hat n=\langle \hat n^2\rangle -\langle \hat n\rangle ^2$ is the variance of the population imbalance operator. If the fluctuations between the two regions are suppressed below the shot-noise level defined as ξ² = 1, the system is called 'number-squeezed'. In our case, since $\hat n$ does not depend on time and the initial state is a vacuum, we obtain that ξ² ≡ 0. Therefore, the two-region Bogoliubov system is perfectly number-squeezed, as anticipated in the previous section.

The ideal number-squeezing is a result of clear separation of the two scattering regions. In such case, it is natural to define the local atom-number operators (17) and the population imbalance operator (18). It is important to note that not all systems, where particles are scattered in pairs are perfectly number squeezed. For instance, when two BECs collide, they produce a halo of atoms due to two-body elastic scattering into the initially unoccupied modes [17, 25]. In this system however, there is no simple way to define two separate regions. One can instead measure the number of atoms in two bins lying on the opposite sides of the halo. Moderate number-squeezing of the atom number difference between these bins has been observed experimentally [19], but it is impossible to reach the limit ξ² = 0 [36]. In contrast, the twin matter wave configurations [22, 24, 28–30], are ideal sources of correlated atomic pairs occupying two well-defined areas.

Apart from the number squeezing, the twin-region system can be characterized by another expression, which is called the Cauchy–Schwarz inequality. It relates the strength of the local and opposite correlations to witness the pair-scattering process. Following [20], we define averaged second-order correlations as

$$\begin{equation} \mathcal{G}^{(2)}_{\mu \nu} \equiv \int\nolimits_\mu \mathrm{d}{\mathbf k}_1 \int\nolimits_\nu \mathrm{d}{\mathbf k}_2 \langle\hat\delta^\dagger({\mathbf k}_1,t)\hat\delta^\dagger({\mathbf k}_2,t) \hat\delta({\mathbf k}_2,t)\hat\delta({\mathbf k}_1,t)\rangle, \end{equation} \tag{ 20 }$$

where μ,ν∈{R,L}. In the symmetric case, the Cauchy–Schwarz inequality $\mathcal G^{(2)}_{\mathrm {RL}} \leqslant ( \mathcal G^{(2)}_{\mathrm {RR}} \mathcal G^{(2)}_{\mathrm {LL}})^{1/2}$ can now be re-written as

$$\begin{equation} \int\nolimits_{\rm R}\!\!\mathrm{d}{\mathbf k}_1\!\!\int\nolimits_{\rm L}\!\!\mathrm{d}{\mathbf k}_2 |M({\mathbf k}_1,{\mathbf k}_2;t)|^2 \leqslant \int\nolimits_{\rm R}\!\!\mathrm{d}{\mathbf k}_1\!\!\int\nolimits_{\rm R}\!\!\mathrm{d}{\mathbf k}_2|G^{(1)}({\mathbf k}_1,{\mathbf k}_2;t)|^2. \end{equation} \tag{ 21 }$$

Using expressions (11) and (12) we obtain

$$\begin{equation} \mathcal{G}^{(2)}_{\mathrm{RL}}=\left(\sum_in_i\right)^2+\sum_in_i(n_i+1),\end{equation} \tag{ 22 }$$

$$\begin{equation} \mathcal{G}^{(2)}_{\mathrm{RR}}=\left(\sum_in_i\right)^2+\sum_in_i^2, \end{equation} \tag{ 23 }$$

thus the Cauchy–Schwarz inequality reads

$$\begin{equation} \sum_in_i(n_i+1)\leqslant\sum_in_i^2 \end{equation} \tag{ 24 }$$

which is true only for all n_i = 0. As soon as particles start to scatter into the two regions, the Cauchy–Schwarz inequality is clearly violated. To quantify the degree of violation, a coefficient $\mathcal {C}$ was introduced in [20], which reads

$$\begin{equation} \mathcal{C}=\frac{\mathcal{G}^{(2)}_{\mathrm{RL}}}{\mathcal{G}^{(2)}_{\mathrm{RR}}}. \end{equation} \tag{ 25 }$$

When $\mathcal {C}\leqslant 1$ the system is in the 'classical' regime, while $\mathcal {C}>1$ signify correlations which are stronger than allowed by the classical physics. In our case this coefficient reads

$$\begin{equation} \mathcal{C}=1+\frac{\sum_in_i}{\left(\sum_in_i\right)^2+\sum_in_i^2}. \end{equation} \tag{ 26 }$$

Clearly, always $\mathcal {C}>1$, because it is a sum of unity and a non-negative part. For high mode populations n_i, the second term, which is inversely proportional to the number of scattered particles tends to zero, restoring the classical limit. Nevertheless, as demonstrated with photons in [37], the confidence by which the Cauchy–Schwarz inequality can be violated in the presence of classical noise still increases with more strongly populated modes.

It has been recently demonstrated that the violation of the Cauchy–Schwarz inequality is a proof of particle entanglement in all systems of identical bosons [21]. The relation holds for every case, when either the number of particles is fixed or fluctuates from shot to shot, as happens in systems described by the Bogoliubov theory, for as long as coherences between different number states are absent. However in the high-gain regime, it is the Fisher information, which is the quantity more sensitive to particle entanglement then the Cauchy–Schwarz criterion, as we show in the following section. The Fisher information quantifies the potential for sub-shot-noise interferometry, and increases with rising mode population, in spite of the decreasing 'granularity' of the matter wave [35] that leads to all second-order correlation functions approaching equal values.

We now show that atoms occupying the two regions are entangled, and could be used as an input of a quantum interferometer operating below the shot-noise level. We first recall how the precision of the phase estimation is related to the entanglement of input states using as an example the standard two-mode Mach–Zehnder interferometer (MZI). Then, we extend these concepts to the case, where the interferometer operates between two regions, each having a multi-mode structure determined by the Bogoliubov equations.

When speaking about two-mode interferometers, it is convenient to introduce a set of three operators

$$\begin{equation} \skew4\hat{J}_x=\textstyle\frac{1}{2}\left(\hat a^\dagger_{\rm R}\hat a^{\phantom{\dagger}}_{\rm L}+\hat a^\dagger_{\rm L}\hat a^{\phantom{\dagger}}_{\rm R}\right),\end{equation} \tag{ 27 }$$

$$\begin{equation} \skew4\hat{J}_y=\frac{1}{2i}\left(\hat a^\dagger_{\rm R}\hat a^{\phantom{\dagger}}_{\rm L}-\hat a^\dagger_{\rm L}\hat a^{\phantom{\dagger}}_{\rm R}\right),\end{equation} \tag{ 28 }$$

$$\begin{equation} \skew4\hat{J}_z=\textstyle\frac{1}{2}\left(\hat a^\dagger_{\rm R}\hat a^{\phantom{\dagger}}_{\rm R}-\hat a^\dagger_{\rm L}\hat a^{\phantom{\dagger}}_{\rm L}\right) \end{equation} \tag{ 29 }$$

which obey the same commutation relations as the angular momentum operators. The MZI, which is an interferometric device, where the imprint of the phase θ onto the input state is preceded and followed by a pair of symmetric beam-splitters, can be represented by a unitary evolution operator $\hat U(\theta )=\mathrm {e}^{-\mathrm {i}\theta \skew4\hat{J}_y}$. If the phase is estimated in a series of ν ≫ 1 measurements performed on the output state, the precision of the phase estimation is limited by the Cramer–Rao lower bound (CRLB) [38, 39]

$$\begin{equation} \Delta\theta\geqslant\frac1{\sqrt\nu}\frac1{\sqrt F_{\rm Q}}. \end{equation} \tag{ 30 }$$

Here, F_Q is the quantum Fisher information (QFI), which is related to the unitary transformation $\hat U(\theta )$. For pure states transformed by the MZI it is equal to $F_{\mathrm { Q}}=4\,\Delta ^2\!\skew4\hat{J}_y$, where the variance is calculated in the input state of the interferometer [40]. The CRLB states, that if θ is determined using any possible type of measurement and estimator, then the precision Δθ is bounded as in equation (30).

Apart from providing a lower bound for the error of the phase estimation, the F_Q is an entanglement measure. Namely, when the input state has an average number of $\langle \hat N\rangle$ particles, then if $F_{\mathrm { Q}}>\langle \hat N\rangle$, the state is particle-entangled [41].

We now show, that a natural extension of the two-mode picture allows to employ the concept of the QFI as an entanglement measure also in our multi-mode system of interest. To this end, we introduce the following analogue of the two-mode angular momentum operators (27)–(29):

$$\begin{equation} \skew4\hat{J}_x=\frac{1}{2}\int\nolimits_{\rm R}\!\!{\mathrm d}{\mathbf k}\,\left(\hat\delta^\dagger({\mathbf k})\hat\delta(-{\mathbf k})+\hat\delta^\dagger(-{\mathbf k})\hat\delta({\mathbf k})\right),\end{equation} \tag{ 31 }$$

$$\begin{equation} \skew4\hat{J}_y=\frac{1}{2\,{i}}\int\nolimits_{\rm R}\!\!{\mathrm d}{\mathbf k}\,\left(\hat\delta^\dagger({\mathbf k})\hat\delta(-{\mathbf k})-\hat\delta^\dagger(-{\mathbf k})\hat\delta({\mathbf k})\right),\end{equation} \tag{ 32 }$$

$$\begin{equation} \skew4\hat{J}_z=\frac{1}{2}\int\nolimits_{\rm R}\!\!{\mathrm d}{\mathbf k}\,\left(\hat\delta^\dagger({\mathbf k})\hat\delta({\mathbf k})-\hat\delta^\dagger(-{\mathbf k})\hat\delta(-{\mathbf k})\right), \end{equation} \tag{ 33 }$$

where we dropped the explicit time-dependence of the $\hat \delta ({\mathbf k},t)$ to simplify the notation. Also, for simplicity, we choose the well-separated regions R and L to be localized symmetrically on the opposite sites of k = 0. The construction of these operators, which satisfy the same commutation relations as (27)–(29), is based on the analogy between the two-mode systems and the twin-beam configuration. In the former case, the operators connect the right and left modes, while in the latter the left and right sub-spaces. Such a definition (31)–(33) is meaningful only in situations, where the system consists of two well-separated sub-systems.

Using the decomposition of the field operator into the set of independent modes, equations (8)–(10), the above integrals yield, that each angular momentum operator is a sum of operators acting on each mode independently, that is

$$\begin{equation} \skew4\hat{J}_x=\sum_i\frac{1}{2}\left(\hat a^{(i)\dagger}_{\rm R}\hat a^{(i)}_{\rm L}+\hat a^{(i)\dagger}_{\rm L}\hat a^{(i)}_{\rm R}\right)\equiv\sum_i\skew4\hat{J}_x^{(i)},\end{equation} \tag{ 34 }$$

$$\begin{equation} \skew4\hat{J}_y=\sum_i\frac{1}{2{i}}\left(\hat a^{(i)\dagger}_{\rm R}\hat a^{(i)}_{\rm L}-\hat a^{(i)\dagger}_{\rm L}\hat a^{(i)}_{\rm R}\right)\equiv\sum_i\skew4\hat{J}_y^{(i)},\end{equation} \tag{ 35 }$$

$$\begin{equation} \skew4\hat{J}_z=\sum_i\frac{1}{2}\left(\hat a^{(i)\dagger}_{\rm R}\hat a^{(i)}_{\rm R}-\hat a^{(i)\dagger}_{\rm L}\hat a^{(i)}_{\rm L}\right)\equiv\sum_i\skew4\hat{J}_z^{(i)}. \end{equation} \tag{ 36 }$$

These expressions show again that it is natural to describe the two-region system using the diagonal basis (8). In this language, the angular momentum operators are simply a sum of operators acting on each pair of modes independently, which vastly simplifies the further analysis.

To establish a direct relation between the two-mode and two-region case, we now assume that the system is transformed in the multi-mode analogue of the MZI. As outlined above, to demonstrate the presence of useful entanglement between atoms in the left and in the right, it is necessary to calculate the QFI. Using equation (34) we obtain that

$$\begin{equation} F_{\rm Q}=4\,\Delta^2\!\skew4\hat{J}_y=4\left\langle\left(\sum_i\skew4\hat{J}_y^{(i)}\right)^2\right\rangle= 4\sum_i\left\langle\left(\skew4\hat{J}_y^{(i)}\right)^2\right\rangle+4\sum_{i\neq j}\left\langle\skew4\hat{J}_y^{(i)}\skew4\hat{J}_y^{(j)}\right\rangle. \end{equation} \tag{ 37 }$$

Since the construction of the basis (8) explicitly assumes that each mode is independent from all others, the second term in the last equality is $4\sum _{i\neq j}\langle \skew4\hat{J}_y^{(i)}\rangle \langle \skew4\hat{J}_y^{(j)}\rangle =0$, because the symmetry between the R and L regions implies that $\langle \skew4\hat{J}_y^{(i)}\rangle =0$ for all i. Therefore we obtain that the QFI is equal to

$$\begin{equation} F_{\rm Q}=4\sum_i\left\langle\left(\skew4\hat{J}_y^{(i)}\right)^2\right\rangle=4\sum_in_i^2+2\left\langle\hat N\right\rangle, \end{equation} \tag{ 38 }$$

where the last equality comes directly from the substitution of (9) and (10) into the definition of the $\skew4\hat{J}_y^{(i)}$ operator. Also, we used $\langle \hat N\rangle =2\sum _in_i$, according to equation (14). Clearly $F_{\mathrm { Q}}>\langle \hat N\rangle$, so the system is entangled. Moreover, one can refer the QFI to the ultimate bound for the precision of the parameter estimation, which is the Heisenberg limit. For a system with fluctuating number of particles, this upper bound is equal to $\langle \hat N^2\rangle$. Using (9) and (10) again, we obtain that [41]

$$\begin{equation} \left\langle\hat N^2\right\rangle=\left\langle\left(\sum_i\left(\hat a^{(i)\dagger}_{\rm R}\hat a^{(i)}_{\rm R}+\hat a^{(i)\dagger}_{\rm L}\hat a^{(i)}_{\rm L}\right)\right)^2\right\rangle= 8\sum_in_i^2+\left\langle\hat N\right\rangle. \end{equation} \tag{ 39 }$$

For a large number of scattered particles, when $\langle \hat N\rangle \ll \sum _in_i^2$, we obtain $F_{\mathrm { Q}}\simeq \frac 12\langle \hat N^2\rangle$. The value of the QFI, which is only one-half smaller than the Heisenberg limit is a clear indication of very strong entanglement present in the system in the high-gain regime. At intermediate times, $F_{\mathrm { Q}}<\frac 12\langle \hat N^2\rangle$ due to mode competition, which has a negative impact on the entanglement as witnessed by the QFI [42]. To picture this, consider a 'frustrated case', where all atoms scatter uniformly into M pairs of modes, so that all n_i ≡ n are equal. In this case, the number of scattered atoms is simply $\langle \hat N\rangle =2nM$, and the QFI is $F_{\mathrm { Q}}=4n^2M+2\langle \hat N\rangle$. The QFI normalized to the SNL is

$$\begin{equation} \frac{F_{\rm Q}}{\langle\hat N\rangle}=2+\frac{\langle\hat N\rangle}{M}. \end{equation} \tag{ 40 }$$

When, on average, there is less than a particle per a set of modes, i.e. $\frac {\left \langle \hat N\right \rangle }{M}\ll 1$, the QFI surpasses the SNL only by a factor of 2, a natural reminiscence of atoms being scattered in pairs. Equation (40) is a simple yet intuitive estimation of the lower bound of useful entanglement in terms of the number of scattered atoms and occupied modes.

Although usually it is impossible to directly measure the QFI, one can employ some other estimation scheme, which gives Δθ close to the CRLB from equation (30). For instance, in a two-region setup as in [30], the entanglement was detected through the phase sensitivity in a two-step protocol. First, the particles from the two modes were mixed using a radio-frequency pulse, which can be represented by the unitary operation

$$\begin{equation} \hat U_{\rm mix}(\theta)=\mathrm{e}^{-\mathrm{i}\theta\skew4\hat{J}_x}, \end{equation} \tag{ 41 }$$

where the $\skew4\hat{J}_x$ operator is defined in equation (31). Then, the population imbalance between the two regions was measured for different values of θ and the phase sensitivity was estimated using the error propagation formula basing on the second moment of the population imbalance operator. Introducing the operator $\hat A=\skew4\hat{J}_z^2$ the relation reads

$$\begin{equation} \Delta^2\theta=\frac1{m}\frac{\langle(\Delta\hat A)^2\rangle}{(\frac{\partial\langle \hat A\rangle}{\partial\theta})^2}. \end{equation} \tag{ 42 }$$

The drop of Δ²θ below the shot-noise level in the experiment of [30] signaled entanglement between the scattered particles. A similar method could be employed for the twin-beam experiment [29], where the particle mixing operation should be implemented using for instance an optical Bragg pulse rather than a radio-frequency field.

The problem of finding the optimal measurement, which saturates the CRLB (30), is one of the central issues of quantum estimation theory. For the twin-beam setup and in absence of decohering influence of the environment, the estimation scheme presented above is optimal for small mixing angles. However, as soon as the noise starts to affect the system, one must seek for another phase determination protocols. Nevertheless, usually it is very difficult to find a truely optimal scheme, which would be easy to implement in the laboratory.

We begin by describing the procedure employed in [29] to produce correlated atom pairs. First, an almost pure BEC of N₀ ≈ 800 ⁸⁷Rb atoms with scattering length equal to a = 5.3 nm was created at temperature T ≈ 25 nK. The cloud was trapped in a slightly unharmonic potential, which for our purposes can well be approximated by

$$\begin{equation} V({\mathbf r}) \simeq \textstyle\frac{1}{2}m\omega_x^2x^2+\textstyle\frac{1}{2}m\omega_y^2y^2+\textstyle\frac{1}{2}m\omega_z^2z^2, \end{equation} \tag{ 43 }$$

where atomic mass is equal to m = 1.44 × 10⁻²⁵ kg, and the frequency ω_x = 2π × 16.3 Hz is much smaller than ω_y = 2π × 1.83 kHz and ω_z = 2π × 2.50 kHz, so the BEC is strongly elongated along the x-axis.

After the BEC was created, the trapping potential was shaken in a controlled way, so the atoms were transferred to the first excited state along the y-direction. In order to achieve the maximal transfer efficiency, the shaking was optimized using quantum optimal control theory [43]. Afterwards, binary collisions transfered particle pairs to the ground state of the potential, and the excess energy 2ℏω_y was converted into back-to-back movement of the two atoms along x. Momentum conservation ensured, that their wave vectors had equal lengths $k_0\approx \sqrt {2m\omega _y/\hbar }$ and point in opposite directions. Small corrections to the value of k₀ may arise from an effective mean-field potential, as will be discussed below.

Neglecting thermal phase fluctuations along the elongated direction x, which is valid at very low temperatures only [44], the condensate wave function acting as a source for the pair-production can be found by solving the stationary GPE

$$\begin{equation} \mu\,\psi({\mathbf r})=\left(-\frac{\hbar^2\nabla^2}{2m}+V({\mathbf r})+g|\psi({\mathbf r})|^2\right)\psi({\mathbf r}), \end{equation} \tag{ 44 }$$

where μ is the chemical potential. This function can be evaluated numerically, by referring to the description of the experiment from the previous section, and noting that after the shaking of the trap, the BEC is in the first excited state n_y = 1 along the y-axis and in the ground state n_z = 0 along z. However, this can be approximated by an analytical expression, as argued below.

First note, that since the characteristic energies ℏω_y and ℏω_z are large, and the number of atoms in the BEC is small, the nonlinear term can be safely neglected in evaluation of the eigenstates along y and z. As a result, assuming that the total wave-function ψ(r) separates in three directions (which has been confirmed numerically), we obtain

$$\begin{equation} \psi({\mathbf r}) = \phi(x)\times \psi^{\rm (ho)}_{n_y=1}(y)\times \psi^{\rm (ho)}_{n_z=0}(z), \end{equation} \tag{ 45 }$$

where the functions ψ^(ho)_ny = 1(y) and ψ^(ho)_nz = 0(z) are the eigen-states of the 1D harmonic potential in y and z correspondingly. The function ϕ(x) is found by inserting the above expression into equation (44) and integrating out the orthogonal directions. As a result, we obtain an effective equation

$$\begin{equation} (\mu-\epsilon_\perp)\, \phi(x) = \left( -\frac{\hbar^2}{2m}\frac{ \partial^2}{\partial x^2} + \tilde g |\phi(x)|^2 + \frac{1}{2}\,m \omega_x^2x^2\right)\phi(x), \end{equation} \tag{ 46 }$$

where zero-point energy equals $\epsilon _\perp = \frac 32\omega _y + \frac 12\omega _z$ and the nonlinearity reads

$$\begin{equation} \tilde g=g \left[\int\!\! \mathrm{d}y\, \left(\psi^{\rm (ho)}_{n_y=1}(y)\right)^4\right]\times \left[\int\!\! \mathrm{d}z\, \left(\psi^{\rm (ho)}_{n_z=0}(z)\right)^4\right]= \frac{3}{8\pi}\frac g{a_{\mathrm{ho},y}\, a_{\mathrm{ho},z}}. \end{equation} \tag{ 47 }$$

Here $a_{\mathrm {ho},i}=\sqrt {\frac {\hbar }{m\omega _i}}$ are the harmonic oscillator lengths for i = y, z. Since the trap is shallow in the x-direction, the solution of the stationary GPE (46) can be well approximated by the Thomas–Fermi (TF) formula [45]

$$\begin{equation} \phi(x) =\sqrt{ \frac{\tilde\mu}{\tilde g}}\sqrt{1 - \frac{x^2}{R_{\rm tf}^2}}, \end{equation} \tag{ 48 }$$

where the effective chemical potential is

$$\begin{equation} \tilde \mu =\mu-\epsilon_\perp =\left( \frac{3 \tilde g N_0 \sqrt{m \omega_x^2 }}{ 4\sqrt{2} } \right)^{2/3}=492\,{\rm Hz}\times h, \end{equation} \tag{ 49 }$$

leading to a TF radius of .

Within the approximation of neglecting thermal phase fluctuations, we have fully determined the wave-function of the BEC, which we insert into the Bogoliubov Hamiltonian (3). Next, we expand the field operator $\hat \delta ({\mathbf r})$ in an orthonormal basis. Along the y- and z-directions, it is natural to use the eigen-states of the harmonic oscillator as the basis functions, since it matches the geometry of the source BEC. Along the x-direction, we use a plane-wave basis, and get

$$\begin{equation} \hat \delta({\mathbf r},t) = \sum_{n_y,n_z}\int\!\frac{\mathrm{d}k}{2\pi} \mathrm{e}^{\mathrm{i} k x} \psi^{\rm (ho)}_{n_y}(y)\psi^{\rm (ho)}_{n_z}(z) \hat \delta(k,n_y,n_z,t). \end{equation} \tag{ 50 }$$

Since the atom pairs are emitted into the ground state along y only (which is ensured by the anisotropy and anharmonicity of the potential), the sum over the eigen-states can be safely truncated at n_y = 0 and n_z = 0. This reduces the dynamics of the pair-production to 1D problem along the x-axis, with the orthogonal directions frozen out, i.e.

$$\begin{equation} \hat \delta({\mathbf r},t) \simeq \psi^{\rm (ho)}_{n_y=0}(y)\psi^{\rm (ho)}_{n_z=0}(z)\sum_k \frac{\mathrm{e}^{\mathrm{i} k x}}{\sqrt L}\, \hat a_k(t), \end{equation} \tag{ 51 }$$

where L is the quantization volume. We insert this field operator into equation (3), evaluate the spatial integrals and upon the change of variables $\hat a_k(t)\mathrm {e}^{ \mathrm {i} \tilde \mu t}\rightarrow \hat a_k(t)$ obtain

$$\begin{eqnarray} &&\fl \skew2\hat{H}_{\rm bog}\simeq\sum_k\!\left( \frac{\hbar^2 k^2}{2m}-(\tilde\mu+\hbar\omega_y)\!\right)\hat a^\dagger_k(t)\hat a^{\phantom{\dagger}}_k(t)+\!\sum_{k,k'}\!\left(2f_{k-k'}\,\hat a^\dagger_k(t)\hat a^{\phantom{\dagger}}_{k'}(t)+f_{k+k'}\!\left(\hat a^\dagger_k(t)\hat a^{{\dagger}}_{k'}(t)+{\rm h.c.}\right)\!\right)\!,\nonumber\\ \end{eqnarray} \tag{ 52 }$$

where $f_q=\frac 23\frac {\tilde g}L\!\int \!\mathrm {d}x\, \mathrm {e}^{-\mathrm {i}qx} \phi ^2(x)$. We solve the resulting Bogoliubov equation numerically, by taking the exponent of the evolution matrix, and find the matrices $\hat {\mathcal {C}}$ and $\hat {\mathcal {S}}$ as defined in equation (7).

Using the above Hamiltonian, one can also analytically determine k₀, i.e. the position of the central peak. To this end, we employ a two-mode approximation by replacing the function f_q with a Dirac delta, and obtain the Bogoliubov equation

$$\begin{equation} \mathrm{i}\,\hbar \partial_t \hat \delta_k(t)=\frac{\hbar^2}{2m}\left(k^2-k_0^2\right)\hat \delta_k(t)+\frac23\,\tilde\mu\,\hat\delta^\dagger_{-k}(t), \end{equation} \tag{ 53 }$$

where k₀ is shifted with respect to the harmonic excitation energy due to the mean-field repulsion and reads

$$\begin{equation} k_0= \sqrt{\frac{2m}{\hbar^2}\left(\hbar\omega_y-\frac{\tilde\mu}3\right)}=5.35\,\mu{\rm m}^{-1}. \end{equation} \tag{ 54 }$$

This result is in good agreement with the experimentally measured position of the peak density, i.e. k_0,exp = 5.55(5) μm⁻¹, where the error arises mainly from the mean-field fluctuations.

In this section, we display the most important characteristics of the twin-beam system, starting from the solution of the eigen-problem of the density matrix (11). In figure 1 we plot the first four eigen-values of the density matrix, as a function of time. The inset shows the average number of scattered atoms $\langle \hat N\rangle$ normalized to the occupation of the BEC, as a function of time. The Bogoliubov approximation is valid for as long as $\langle \hat N\rangle \ll N_0$, so we interrupt the simulation at t = 1.2 ms, when $\langle \hat N\rangle \simeq 15\%N_0$. For longer times, when the depletion of the BEC cannot be neglected, a atom-number conserving method, such as the one introduced in [46] must be used.

Zoom In Zoom Out Reset image size

In figure 2 we plot the first four eigen-vectors of G⁽¹⁾ localized in the right half-space, i.e. |φ⁽ⁱ⁾_R(k)|² with i = 1,2,3,4, calculated at an early time t = 0.1 ms and at t = 1.2 ms. Due to the time–energy uncertainty relation, the eigen-modes localize around k = k₀ at later times.

Zoom In Zoom Out Reset image size

This can be seen even more clearly, by plotting the density ρ(k;t) at these two instants, as shown in figure 3 (dashed lines). At t = 0.1 ms, two broad beams start to form on top of the uniform density. Later, at t = 1.2 ms, strongly localized peaks clearly dominate over the flat background. On top of these curves, we plot the normalized second-order correlation function as defined in equation (15), with one of the arguments set equal to the resonant wave-vector k₀, i.e. g⁽²⁾(k₁,k₂ ≡ k₀;t). At t = 0.1 ms, the cross-correlation, which is governed by the anomalous density, is very large, i.e. g⁽²⁾(− k₀,k₀;0.1 ms) ≃ 40. This is a characteristic property of the Bogoliubov system in the low-occupation regime [36], and indicates strong violation of the Cauchy–Schwarz inequality (21). Also, for this early time, the width of both g⁽²⁾ peaks are much more narrow than the beam size. This is consistent with the results shown in figure 1, where at early times many eigen-mode pairs of the density matrix are almost equally occupied. At later times, when a single pair of modes start to become dominant, the width of the peak in g⁽²⁾ and the system size approach each other. While this corresponds to beams that are single-mode with respect to their local one-body properties, the local averaged correlation function as introduced in equation (20) reaches the limit of $\mathcal {G}^{(2)}_{\mu \mu } \simeq 2(\sum _i n_i)^2 \equiv \frac 12 \langle \hat N \rangle ^2$, exceeding the number fluctuations of a coherent state by a factor of two.

Zoom In Zoom Out Reset image size

The full correlation function G⁽²⁾(k₁,k₂) at t = 1.2 ms, shown in figure 4, reveals both the local and cross correlations. In typical experiments, however, such a function is difficult to measure, and one uses collinearly integrated functions of the type

$$\begin{eqnarray} g_{\mathrm{cl}}^{(2)}(\delta {\mathbf k};t) &=& \frac{\int G^{(2)}({\mathbf k},{\mathbf k}+\delta {\mathbf k};t) \mathrm{d}{\mathbf k}}{ \int \rho({\mathbf k};t) \rho({\mathbf k}+\delta{\mathbf k};t) \mathrm{d}{\mathbf k}}, \end{eqnarray} \tag{ 55 }$$

where the integrals run over an appropriately chosen momentum region [17, 20]. For symmetric, non-local correlations, back-to-back integration of the type

$$\begin{eqnarray} g_{\mathrm{bb}}^{(2)}({\mathbf k}_+;t) &=& \frac{\int G^{(2)}({\mathbf k},{\mathbf k}_+- {\mathbf k};t) \mathrm{d}{\mathbf k}}{ \int \rho({\mathbf k};t) \rho({\mathbf k}_+ - {\mathbf k};t) \mathrm{d}{\mathbf k}} \end{eqnarray} \tag{ 56 }$$

is used. The corresponding normalized functions for our system g⁽²⁾_cl(δ_k;t),g⁽²⁾_bb(δ_k;t) are shown in figure 5 at t = 1.2 ms. For both local (solid) and non-local (dotted) functions, correlations peaks, which do not span the entire populated range (grey area) and follow a Gaussian shape, are clearly present.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the next step we take toward future comparison with experiments, we present the results not in momentum space, but rather using real-space data calculated after some finite time τ of ballistic expansion. Only in the limit of τ → ∞ (far field), the real-space data is equivalent to the initial momentum space distribution (if the expanding clouds are sufficiently dilute, so that the mean-field repulsion can be safely neglected). In [29], the expansion time was τ = 46 ms, which was sufficient to resolve the twin-beam peaks. Nevertheless, the system was not fully in the far-field regime yet, which may have some impact on the correlation functions. As shown in figure 5, the finite expansion time affects the back-to-back peak at (k₀,−k₀) much more strongly than the collinear HBT peak, leading to smearing of the measured g⁽²⁾_bb(k₊;t) (dash-dotted line) over the entire size of the twin-beam packets. This observation is consistent with some previous numerical results [47].

The different behavior of these two types of correlations in the near field can be explained using following semi-classical arguments. First, consider a single body function G⁽¹⁾(k₁,k₂), which according to equation (16) determines the local second-order correlations. The time, when the system enters the far-field regime, corresponds to an instant, when the measured density in position space samples the momentum distribution. If a particle is propagating freely, starting from x₀, for a time τ and having velocity v, then it is registered at x₀ + vτ, which recovers the initial momentum if x₀ can be neglected when compared to the characteristic velocity spread σ_vτ. From the width of momentum distribution shown in figure 3, we obtain that σ_v ≈ 290 μm s⁻¹ while the position density profile gives x₀ ≈ 6.3 μm. Therefore, σ_vτ/x₀ ≈ 2.12, which suggests that G⁽¹⁾ after τ = 46 ms of ballistic expansion is approximately in the far-field regime.

On the other hand, the anomalous density M(k₁,k₂) is a two-body function, which governs the opposite correlations between a pair of atoms. If two atoms are scattered from the source at position $\tilde x$ with velocities v₁ and v₂, then after the ballistic expansions time τ their total coordinate is $x_1+x_2 = (v_1+v_2)\tau + 2\tilde x$. In the far-field regime, the position measurement samples momentum correlations meaning that the first term v₊ = v₁ + v₂ dominates over $2\tilde x/\tau$. The characteristic range of v₊ can be estimated from figure 5 giving v₊ ≈ 146 μm s⁻¹. The range of allowed positions of the collisions is given by the size of the mother cloud, thus $\tilde x\approx R_{\mathrm {tf}} \approx 21\,\mu {}\mathrm {m}$. Therefore, we obtain that $v_+ \tau / 2\tilde x\approx 0.16$, which shows that the anomalous density after τ = 46 ms—and in consequence the opposite momentum correlations—is not in the far-field regime. We conclude that the broadening of the opposite correlations results mainly from the spread in the positions of the scattering events.

Note that although at every instant of the evolution, the field operator $\hat \delta ({\mathbf k},t)$ can be written as a sum of independently squeezed modes, at very early times the division between the right and left modes is unjustified, because the two peaks are not yet fully separated. However, at t = 0.1 ms when the density distribution is broad, the number of scattered atoms is $\langle \hat N\rangle \simeq 2$. Therefore, the system at such early time is hardly accessible experimentally so the quantum state of much less then a single particle is not of interest. As soon as the two peaks are well-formed, at t ≈ 0.3 ms, with $\langle \hat N\rangle \simeq 8$ scattered atoms, all the general considerations from section 2 apply.

Finally, in figure 6 we plot the QFI from equation (38) as a function of time and normalized to the Heisenberg limit, i.e. $F_{\mathrm { Q}}\big /\langle \hat N^2\rangle$. Instead of interrupting the simulation at 1.2 ms, where the scattered fraction of atoms becomes non-negligible and particle number conservation is strongly violated, we extend the calculation up to 7 ms, when the number of scattered atoms significantly exceeds 15% of N₀. This is done solely to illustrate that, once the population of one of the modes dominates, $F_{\mathrm { Q}}\rightarrow \frac 12\langle \hat N^2\rangle$, as argued in section 2.5. Note that the dominance of a single mode pair at long times is also predicted by the number-conserving theory [46], justifying this proceeding. Indeed, in the inset, we show the number of pairs of right/left modes which have an occupation bigger or equal to 10% of the largest mode. This approximately tells, how many modes are significantly occupied in the system. At early times, there are over 100 pairs of modes. At 1.2 ms there are still five significantly occupied pairs, and only around 4.2 ms a single pair of modes starts to dominate. At the same time the QFI approaches its upper bound.

Zoom In Zoom Out Reset image size