Bogoliubov theory for atom scattering into separate regions

We review the Bogoliubov theory in the context of recent experiments, where atoms are scattered from a Bose-Einstein Condensate into two well-separated regions. We find the full dynamics of the pair-production process, calculate the first and second order correlation functions and show that the system is ideally number-squeezed. We calculate the Fisher information to show how the entanglement between the atoms from the two regions changes in time. We also provide a simple expression for the lower bound of the useful entanglement in the system in terms of the average number of scattered atoms and the number of modes they occupy. We then apply our theory to a recent"twin-beam"experiment [R. B\"ucker {\it et al.}, Nat. Phys. {\bf 7}, 608 (2011)]. The only numerical step of our semi-analytical description can be easily solved and does not require implementation of any stochastic methods.

We review the Bogoliubov theory in the context of recent experiments, where atoms are scattered from a Bose-Einstein Condensate into two well-separated regions. We find the full dynamics of the pair-production process, calculate the first and second order correlation functions and show that the system is ideally number-squeezed. We calculate the Fisher information to show how the entanglement between the atoms from the two regions changes in time. We also provide a simple expression for the lower bound of the useful entanglement in the system in terms of the average number of scattered atoms and the number of modes they occupy. We then apply our theory to a recent "twin-beam" experiment [R. Bücker et al., Nat. Phys. 7, 608 (2011)]. The only numerical step of our semi-analytical description can be easily solved and does not require implementation of any stochastic methods.

I. INTRODUCTION
In recent years, systems where strong correlations between particles are induced by pair-wise scattering, have attracted much attention. In the canonical example, which is the parametric down-conversion, photon pairs are generated during the propagation of a laser beam through a non-linear medium. The outcoming pairs of photons are entangled, and can serve as a probe of fundamental properties of quantum mechanics [1,2], such as the Einstein-Podolsky-Rosen paradox, or violation of the Bell inequalities [1][2][3]. On the other hand, entanglement can be exploited in practical applications, such as teleportation [4,5] or metrology beyond the Shot-Noise Limit (SNL) [6,7].
In this latter context, recent experiments with entangled states of atoms were a major breakthrough [8]. In [8][9][10], two-body interactions were utilized to prepare nonclassical squeezed states of atoms trapped in a doublewell potential, which implies presence of many-body entanglement [11]. A similar idea was exploited to generate squeezing in the internal [12][13][14] degrees of freedom. In [15,16], squeezing of a large spin of a collection of two-level atoms was achieved, using an intense laser field interacting with particles trapped in an optical cavity.
Simultaneously, a substantial experimental effort was put in order to generate entangled pairs of atoms scattered out of a Bose-Einstein Condensate (BEC). In [17][18][19][20], a collision of two BECs lead to weak scattering of correlated atomic pairs onto a three-dimensional sphere of initially unoccupied modes. Although moderate number-squeezing between the opposite regions of the halo, and the related violation of the Cauchy-Schwarz inequality were experimentally demonstrated [19,20], entanglement was never directly observed. Alternatively, pair-production schemes were developed, where only few modes are strongly populated in a stimulated process, making the system somewhat easier to handle. Stimulated four-wave-mixing processes have been implemented using different spin states of atoms [21][22][23] or Bragg scat-tering [24,25]. Also, dynamic instabilities in moving optical lattices, populating modes with opposite quasimomenta, have been used [26,27]. In [28], a BEC was transferred into the first excited state of a trapping potential and subsequent two-body collisions created a "twin-beam" system, where stronger-than-classical correlations could directly be observed.
Analogous schemes have been implemented in internal atomic states, building upon spin-changing collisions [29][30][31]. Furthermore, in [29] it was shown that particles scattered in this process into a pair of m F = ±1 Zeeman sub-levels are usefully entangled from the metrological point of view.
In this work we develop a theoretical model for the generic type of experiments, where particles scatter in pairs into two well-separated regions. If these regions are separated in the momentum space, they could also be set apart by a sufficiently long expansion of the cloud. On the other hand, in cases when the regions are defined as two Zeeman sub-levels, they can be separated in space using a Stern-Gerlach scheme. Our model applies to any such possible configuration. Therefore, the general conclusions of this study, concerning the correlations, number-squeezing and entanglement, are valid for recent experiments [21][22][23][27][28][29] in the regime, where the depletion of the source BEC is low. We derive the Bogoliubov equations governing the dynamics of pair formation, and applying the Bloch-Messiah reduction [32][33][34], we write the state in terms of pairs of independently squeezed modes. We calculate the density and the number of scattered atoms, and the two body correlation between them. We demonstrate the presence of ideal number-squeezing between the opposite regions, violation of the Cauchy-Schwarz inequality and, using the Fisher information criterion known from quantum metrology, we show that the atoms from the twin-beam system are entangled. We also provide a simple yet useful lower bound for the Fisher information in terms of the average number of scattered atoms, and the number of modes they occupy. Finally, we apply the above formalism to the twin-beam experiment of [28]. This paper is organized as follows. In Sec. II A we discuss the general properties of the solutions of the Bogoliubov equations. These observation allow to easily calculate the density and the second order coherence of the system in Sec. II B and the fluctuations of the population imbalance between the opposite regions in Sec. II C. In Sec. II D we take the first step towards the demonstration of particle entanglement present in the system, by showing that the second order correlation function violates the Cauchy-Schwarz inequality. In Sec. II E we demonstrate that the scattered atoms are usefully entangled from the metrological point of view. In Sec. III A we briefly describe the experimental setup of [28] and in III B derive the corresponding effective Bogoliubov equations. Finally, in Sec. III C we review the most relevant properties of the twin-beam system by showing the results of the numerical simulation. We conclude in Sec. IV.

II. GENERAL PROPERTIES OF THE SCATTERED PARTICLES
We first present the general properties of the solution of the Bogoliubov equation, in cases where particles scatter pair-wise into well-separated regions.

A. Bogoliubov equation for pair scattering
Our theoretical description of the pair-production process starts with a many-body Hamiltonian with contact two-body interactionŝ Here V (r) is an external trapping potential and g = 4π 2 a m is the strength of the two-body interactions, a is the scattering length, m is the atomic mass and the field operator Ψ(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) We then write the field operator as a sum of the c-number part and the Bogoliubov correction,Ψ(r) = ψ(r) +δ(r) and insert this expression into (1). By keeping only the terms up to quadratic inδ we obtain the Bogoliubov Hamiltonian The resulting Bogoliubov equation of motion is linear Usually, a numerical solution of this equation is found in a following way. The field operator is expanded in a basis of wave-functionsφ i (r) which match the geometry of the scattering problem This expression is inserted into Eq. (4), the resulting equation is multiplied byφ * j (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Â andB, couples the evolution of the j-th operatorâ j (t), with (in general) all others operators This equation is linear -a consequence of the linearity of the Bogoliubov equation (4) -so the general solution of (6) readsâ where the matricesĈ andŜ satisfyĈĈ † −ŜŜ † =1 and CŜ T −ŜĈ 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 followŝ The operatorsâ (i) R/L (t) annihilate a particle in a mode characterized by the time-dependent wave function ϕ (i) 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: where |c i (t)| 2 − |s i (t)| 2 = 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) 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 [32][33][34]. First, using Equations (8) and (9), we evaluate the one-body density matrix (first-order correlation function) and obtain where n i = |s i (t)| 2 . Note that n i is the average occupation of the i-th eigen-mode, and that a pair of modes ϕ is degenerate (has the same eigenvalue 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 Eq. (10) 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 (10) in any convenient basis (5), and subsequent diagonalization. Once this is done, then according to Eq. (10), the basis functions ϕ (i) 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 ϕ (i) 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 multiply it by sides with the eigen-functions of the density matrix and integrate over space. As a result, we retrive the information about the phases and obtain the full form of the mode functions ϕ 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.

B. Density and correlations
The simplest observable characterizing the pairproduction process is the density (12) 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 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 1 and k 2 can be obtained from the normalized secondorder correlation function According to the Wick's theorem, this function can be written in terms of the one-body density matrix (10) and the anomalous density (11) as follows (15) The transition from Eq. (14) to (15) 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 Eq. (10), the density matrix is non-vanishing only when k 1 and k 2 are both either in the right or left region, so |G (1) | 2 governs the Hanbury-Brown and Twiss (HBT) type of local correlations. On the other hand, as can be seen from Eq. (11), the anomalous density is non-zero only when k 1 and k 2 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.

C. Number squeezing
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 the 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.
The population imbalance operator is then simply defined asn =N R −N L and using Eq. (8) we obtain The number squeezing factor is defined as where ∆ 2n = n 2 − n 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 ξ 2 = 1, the system is called "number-squeezed". In our case, sincen does not depend on time and the initial state is a vacuum, we obtain that ξ 2 ≡ 0. Therefore, the two-region Bogoliubov system is perfectly numbersqueezed, 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 (16) and the population imbalance operator (17). It is important to note that not all systems, where particles are scattered in pairs are perfectly number squeezed. For instance, when two Bose-Einstein condensates collide, they produce a halo of atoms due to two-body elastic scattering into the initially unoccupied modes [35,36]. 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 [37], but it is impossible to reach the limit ξ 2 = 0 [38]. In contrast, the twin matter wave configurations [21,23,[27][28][29], are ideal sources of correlated atomic pairs occupying two well-defined areas.

D. Violation of the Cauchy-Schwarz inequality
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 In the symmetric case, the Cauchy-Schwarz inequality G Using expressions (10) and (11) we obtain thus the Cauchy-Schwarz inequality reads 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 C was introduced in [20], which reads When C 1 the system is in the "classical" regime, while C > 1 signify correlations which are stronger than allowed by the classical physics. In our case this coefficient reads Clearly, always 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 [39], the confidence by which the Cauchy-Schwarz inequality can be violated in the presence of classical noise still increases with more strongly populated modes. However, it is the Fisher information, which is the quantity highly sensitive to particle entanglement in the highgain regime, as we show in the following section. This measure quantifies the potential for sub-shot-noise interferometry, and increases with rising mode population, in spite of the decreasing "granularity" of the matter wave [34] that leads to all second-order correlation functions approaching equal values.

E. Entanglement and interferometry
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 operatorŝ 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Û (θ) = e −iθĴ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) [40,41], Here, F Q is the Quantum Fisher Information (QFI), which is related to the unitary transformationÛ (θ). For pure states transformed by the MZI it is equal to F Q = 4 ∆ 2Ĵ y , where the variance is calculated in the input state of the interferometer [42]. The CRLB states, that if θ is determined using any possible type of measurement and estimator, then the precision ∆θ is bounded as in Eq. (26).
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 N particles, then if F Q > N , the state is particleentangled [6,7,43,44].
We now show, that a natural extension of the twomode 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 analog of the two-mode angular momentum operators (25), where we dropped the explicit time-dependence of thê δ(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 (25), 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 (27) is meaningful only in situations, where the system consists of two wellseparated sub-systems.
Using the decomposition of the field operator into the set of independent modes, Equations (8) and (9), the above integrals yield, that each angular momentum operator is a sum of operators acting on each mode independently, that iŝ 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 analog of the Mach-Zehnder interferometer. As outlined above, to demonstrate the presence of useful entanglement between the atoms in the left and in the right, it is necessary to calculate the QFI. Using Eq. (28) we obtain that 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 i =j Ĵ (i) y Ĵ (j) y = 0, because the symmetry between the R and L regions implies that Ĵ (i) y = 0 for all i. Therefore we obtain that the QFI is equal to where the last equality comes directly from the substitution of (9) into the definition of theĴ (i) y operator. Also, we used N = 2 i n i , according to Eq. (13). Clearly F Q > N , 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 N 2 . Using (9) again, we obtain, that [44] For a large number of scattered particles, when N ≪ i n 2 i , we obtain F Q ≃ 1 2 N 2 . 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 Q < 1 2 N 2 due to mode competition, which has a negative impact on the entanglement as witnessed by the QFI [45]. 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 N = 2nM , and the QFI is F Q = 4n 2 M + 2 N . The QFI normalized to the SNL is When, on average, there is less than a particle per a set of modes, i.e. N M ≪ 1, the QFI surpasses the SNL only by a factor of 2, a natural reminiscence of atoms being scattered in pairs. Equation (34) is a simple yet intuitive estimation of the lower bound of useful entanglement in terms of the number of scattered atoms and occupied modes.

III. APPLICATION: TWIN-BEAM SYSTEM
We now apply the above formalism to the twin-beam system of [28]. First, we describe the physical mechanism which leads to the creation of the two correlated beams. As shown below, some basic information about the dynamics of the pair production allow to construct a simple one-dimensional Bogoliubov model, which can be easily solved numerically.

A. Scheme of the experiment
The experimental sequence applied in [28] to produce correlated atom pairs was following. First, an almost pure Bose-Einstein Condensate (BEC) of N 0 ≈ 800 87 Rb atoms with scattering length equal to a = 5.3 nm was created at temperature T ≈ 25 nK. The cloud was trapped in an approximately harmonic potential where atomic mass is equal to m = 1.44 × 10 −25 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 [46]. 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 ≈ 2mω y / and point in opposite directions. Small corrections to the value of k 0 may arise from an effective mean-field potential, as will be discussed below.

B. Theoretical description
Neglecting thermal phase fluctuations along the elongated direction x, which is valid at very low temperatures only [47], the condensate wave function acting as a source for the pair-production can be found by solving the stationary GPE 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 non-linear 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 where the functions ψ (ho) ny=1 (y) and ψ (ho) nz =0 (z) are the eigenstates of the one-dimensional harmonic potential in y and z correspondingly. The function φ(x) is found by inserting the above expression into Eq. (36) and integrating out the orthogonal directions. As a result, we obtain an effective equation where zero-point energy equals ǫ ⊥ = 3 2 ω y + 1 2 ω z and the non-linearity reads Here a ho,i = mω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 (38) can be well approximated by the Thomas-Fermi (TF) formula [48] φ where the effective chemical potential is leading to a TF radius of R tf = 2μ mω 2 x = 20.75 µm.
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δ(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 planewave basis, and get δ(r, t) = ny,nz dk 2π e ikx ψ (ho) ny (y)ψ (ho) nz (z)δ(k, n y , n z , t).
(42) 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 eigenstates can be safely truncated at n y = 0 and n z = 0. This reduces the dynamics of the pair-production to onedimensional problem along the x axis, with the orthogonal directions frozen out, i.e.
where L is the quantization volume. We insert this field operator into Eq. (3), evaluate the spatial integrals and upon the change of variablesâ k (t)e iμt →â k (t) obtain where f (q) = 2 3g L dx e −iqx φ 2 (x). We solve the resulting Bogoliubov equation numerically, using a fourth order Runge-Kutta method, and find the matricesĈ andŜ as defined in Eq. (7).
Using the above Hamiltonian, one can also analytically determine k 0 , 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 where k 0 is shifted with respect to the harmonic excitation energy due to the mean-field repulsion and reads This result is in good agreement with the experimentally measured position of the peak density, i.e. k 0,exp = 5.55 µm −1 .

C. Numerical results
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 (10). In Fig. 1 we plot the first four eigen-values of the density matrix, as a function of time. The inset shows the total number of scattered atoms N normalized to the occupation of the BEC, as a function of time. The Bogoliubov approximation is valid for as long as N ≪ N 0 , so we interrupt the simulation at t = 1.2 ms, when N ≃ 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 [49] must be used.
In Fig. 2 we plot the first four eigen-vectors of G (1) localized in the right half-space, i.e. |ϕ This can be seen even more clearly, by plotting the density ρ(k; t) at these two instants, as shown in Fig. 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 (2)  plot the normalized second-order correlation function as defined in Eq. (14), with one of the arguments set equal to the resonant wave-vector k 0 , i.e. g (2) (k 1 , k 2 ≡ k 0 ; t). At t = 0.1 ms, the cross-correlation, which is governed by the anomalous density, is very large, i.e. g (2) (−k 0 , k 0 ; 0.1 ms) ≃ 40. This is a characteristic property of the Bogoliubov system in the low-occupation regime [38], and indicates strong violation of the Cauchy-Schwarz inequality (20). Also, for this early time, the width of both g (2) peaks are much more narrow than the beam size. This is consistent with the results shown in Fig. 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 (2) 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 Eq. (19) reaches the limit of G  (2) (k1, k2 ≡ k0; t) for fixed k2 (solid lines, left yaxis), and density profiles ρ(k1; t) (dashed lines, right axis) in momentum space. The results are calculated at t = 0.1 ms (a) and t = 1.2 ms (b). At early times, many momentum modes are occupied and the width of g (2) is much smaller than the beam size. Later, two distinct peaks emerge, which are almost single-mode.
The typical tool to capture local second-order correlations in experiments are collinearly integrated functions of the type g (2) cl (δk; t) = where the integrals run over an appropriately chosen momentum region [17,20]. For symmetric, non-local correlations, back-to-back integration of the type is used. The corresponding normalized functions for our system g (2) cl (δ k ; t), g c l , τ → ∞ g ( 2 ) b b , τ → ∞ g ( 2 ) c l , τ = 4 6 m s g In the next step we take towards 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 [28], 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 Fig. 5, the finite expansion time affects the back-to-back peak at (k 0 , −k 0 ) much more strongly than the collinear HBT peak, leading to smearing of the measured g (2) bb (k + ; t) (dash-dotted line) over the entire size of the twin-beam packets. This observation is consistent with some previous results [50]. Intuitively, the broadening effect is related to the random position of scattering events along x within the size of the initial cloud, which is non-vanishing with respect to the expanded size of the twin-beam peaks. On the other hand, the local correlation function g (2) cl (δk; t) (dashed line) remains largely unaffected.
Note that although at every instant of the evolution, the field operatorδ(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 N ≃ 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 N ≃ 8 scattered atoms, all the general considerations from Sec. II apply.
Finally, in Fig. 6 we plot the QFI from Eq. (31) as a function of time and normalized to the Heisenberg Limit, i.e. F Q N 2 . 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 0 . This is done solely to illustrate that, once the population of one of the modes dominates, F Q → 1 2 N 2 , as argued in Sec. II E. Note that the dominance of a single mode pair at long times is also predicted by the number-conserving theory [49], 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 5 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.

IV. CONCLUDING REMARKS
We have developed a simple Bogoliubov model describing twin-atom beam experiments similar to Refs. [27,28]. Due to the elongated geometry of the trapping potential, the dynamics is one-dimensional. As a consequence, the final step of our method can be easily solved numerically without the need for any stochastic method. Furthermore, basic information about the scattered particles can directly be drawn from the general properties of the solution of the Bogoliubov equations. In this way, we can quantitatively characterize the mode structure and correlation functions of the scattered atoms. Also, quite generally, we can show that the population imbalance between the two beams is ideally squeezed and that the system is strongly entangled. These general observations can be applied to most recent experiments, where the atomic pairs scatter into two well-separated regions. Finally, using the notion of the Quantum Fisher Information, we have derived a simple lower bound for the useful entanglement of the system. This expression employs only the average number of scattered particles and the number of occupied modes.
Having understood the fundamental properties of fewmode twin beams, further steps can be made to take into account more specific issues of experimental implementations. A general feature of strongly elongated Bose gases at realistic temperatures, such as the source cloud in [28], is the quasi-condensation [47], where the coherence along the x-axis is limited due to thermal phase fluctuations. Although these fluctuations do not alter the general considerations of Sec. II, they affect the emission dynamics [49], and also might have influence on the spatial properties of both the density and the correlation functions [50]. In future, our method could be applied to the analysis of some more complicated schemes, basing on the twin-beam setup, such as the Rarity-Tapster-type experiments [51]. Also, according to our results, the two-region state could be used as an input to the Mach-Zehnder-like interferometer, similarly to [29].
Finally, note that the Bogoliubov approximation neglects the secondary collisions between the scattered particles and the atoms from the source cloud. When a scattered atom propagates through the BEC, the number of secondary collisions is bounded from above by N col = 16πa 2 n R tf , where n is the peak density of the BEC. By plugging in the experimental numbers, we obtain N col = 0.38, which well justifies the use of the Bogoliubov approximation.