Properties of atomic pairs produced in the collision of Bose-Einstein condensates

TL;DR

This paper investigates atomic pair correlations resulting from Bose-Einstein condensate collisions, deriving analytical and numerical tools to understand their quantum properties and potential for nonclassicality demonstrations.

Contribution

It provides a detailed Bogoliubov-based analysis of atomic correlations in condensate collisions, including new analytical expressions and a semiclassical model.

Findings

Derived explicit formulas for one and two-particle correlation functions.

Compared quantum and semiclassical models, highlighting their differences.

Results are applicable to recent helium atom experiments and future nonclassicality tests.

Abstract

Collisions of Bose-Einstein condensates can be used as a mean to generate correlated pairs of atoms. The scattered massive particles, in analogy to photon pairs in quantum optics, might be used in the violation of Bell's inequalities, demonstration of Einstein-Podolsky-Rosen correlations, or sub-shot noise atomic interferometry. Usually, a theoretical description of the collision relies either on stochastic numerical methods or on analytical treatments involving various approximations. Here, we investigate elastic scattering of atoms from colliding elongated Bose-Einstein condensates within Bogoliubov method, carefully controlling performed approximations at every stage of the analysis. We derive expressions for the one and two particle correlation functions. The obtained formulas, which relate the correlation functions with condensate wavefunction, are convenient for numerical calculations. We employ variational approach for condensate wavefunctions to obtain analytical expressions for the correlation functions, which properties we analyze in details. We also present a useful semiclassical model of the process, and compare its results with the quantum one. The results are relevant for recent experiments with excited helium atoms, as well as for planned experiments aimed at investigating the nonclassicality of the system.

Equations

$$i\hslash {\partial }_t\psi (\mathbf{r},t)=\left(-\frac{{\hslash }^2}{2m}△+g|\psi (\mathbf{r},t){|}^2\right)\psi (\mathbf{r},t),$$

$$i\hslash {\partial }_t\hat{\delta }(\mathbf{r},t)=H_0(\mathbf{r},t)\hat{\delta }(\mathbf{r},t)+B(\mathbf{r},t){\hat{\delta }}^{†}(\mathbf{r},t)$$

$$H_0(\mathbf{r},t)$$

$$B(\mathbf{r},t)$$

$$\hat{\delta }(\mathbf{r},0)|0⟩=0.$$

$$M\left({\mathbf{r}}_1,{\mathbf{r}}_2,T\right)=⟨\hat{\delta }({\mathbf{r}}_1,T)\hat{\delta }({\mathbf{r}}_2,T)⟩,$$

$$G^{(1)}\left({\mathbf{r}}_1,{\mathbf{r}}_2,T\right)=⟨{\hat{\delta }}^{†}({\mathbf{r}}_1,T)\hat{\delta }({\mathbf{r}}_2,T)⟩.$$

$$G^{(2)}\left({\mathbf{r}}_1,{\mathbf{r}}_2,T\right)=⟨{\hat{\delta }}^{†}({\mathbf{r}}_1,T){\hat{\delta }}^{†}({\mathbf{r}}_2,T)\hat{\delta }({\mathbf{r}}_2,T)\hat{\delta }({\mathbf{r}}_1,T)⟩$$

$$=G^{(1)}\left({\mathbf{r}}_1,{\mathbf{r}}_1,T\right)G^{(1)}\left({\mathbf{r}}_2,{\mathbf{r}}_2,T\right)+{\left|G^{(1)}\left({\mathbf{r}}_1,{\mathbf{r}}_2,T\right)\right|}^2+$$

$$+{\left|M\left({\mathbf{r}}_1,{\mathbf{r}}_2,T\right)\right|}^2.$$

$$M({\mathbf{r}}_1,{\mathbf{r}}_2,T)$$

$$G^{(1)}\left({\mathbf{r}}_1,{\mathbf{r}}_2,T\right)=\int \text{mbox}d\mathbf{r}{M}^{\ast }\left({\mathbf{r}}_1,\mathbf{r},T\right)M\left(\mathbf{r},{\mathbf{r}}_2,T\right).$$

$$\mu \psi (\mathbf{r})=\left(-\frac{{\hslash }^2}{2m}△+V(\mathbf{r})+g|\psi (\mathbf{r}){|}^2\right)\psi (\mathbf{r}),$$

$$V({\bf r})=\frac{1}{2}m\big{(}\omega_{r}^{2}(x^{2}+y^{2})+\omega_{z}^{2}z^{2}\big{)},$$

$$\psi({\bf r},0)=\frac{C}{2}\big{(}\psi({\bf r})e^{iQz}\!+\!\psi({\bf r})e^{-iQz}\big{)}\!=\!C\psi({\bf r})\cos(Qz),$$

$$\frac{2{\sigma }_r}{Q{a}_{hor}^2}\ll 1,$$

$$\psi (\mathbf{r},0)=\sqrt{2}\psi (\mathbf{r})\cos (Qz).$$

$$\psi (\mathbf{r},t)=\left({\psi }_{+Q}(\mathbf{r},t)e^{iQz}+{\psi }_{-Q}(\mathbf{r},t)e^{-iQz}\right)e^{-i\frac{\hslash Q^2}{2m}t}.$$

$$i\hslash {\partial }_t{\psi }_{\pm Q}(\mathbf{r},t)=\left(\mp i\frac{{\hslash }^2}{m}Q{\partial }_z-\frac{{\hslash }^2}{2m}△\right){\psi }_{\pm Q}(\mathbf{r},t)$$

$$+g\left(|{\psi }_{\pm Q}(\mathbf{r},t){|}^2+2|{\psi }_{\mp Q}(\mathbf{r},t){|}^2\right){\psi }_{\pm Q}(\mathbf{r},t).$$

$${\psi }_{\pm Q}(\mathbf{r},0)=\sqrt{\frac{N}{2{\pi }^{3/2}{\sigma }_z{\sigma }_r^2}}\exp \left(-\frac{x^2+y^2}{2{\sigma }_r^2}-\frac{z^2}{2{\sigma }_z^2}\right).$$

$${\psi }_{\pm Q}(\mathbf{r},t)$$

$${\sigma }_r^2(t)$$

$$a_r(t)$$

$$a_z(t)$$

$$\varphi (t)$$

$$M({\mathbf{k}}_1,{\mathbf{k}}_2)$$

$$B(\mathbf{r},t)=g{\psi }^2(\mathbf{r},t)$$

$$=g\left({\psi }_{+Q}^2{e}^{i2Qz}+2{\psi }_{+Q}{\psi }_{-Q}+{\psi }_{-Q}^2{e}^{-i2Qz}\right)e^{-i\frac{\hslash Q^2}{m}t}.$$

$$B(\mathbf{r},t)=2g{\psi }_{+Q}(\mathbf{r},t){\psi }_{-Q}(\mathbf{r},t)\exp \left(-i\frac{\hslash Q^2}{m}t\right).$$

$$36{\left(\frac{gn}{\frac{{\hslash }^2{Q}^2}{m}}\right)}^{2}Q{\sigma }_r\ll 1,$$

$$M({\mathbf{k}}_1,{\mathbf{k}}_2)=\frac{1}{\hslash (2\pi {)}^3}{\int }_0^{\infty }\text{mbox}dt\int \text{mbox}d\mathbf{r}$$

$$\exp \left(-i({\mathbf{k}}_1+{\mathbf{k}}_2)\mathbf{r}+i\frac{\hslash (k_1^2+k_2^2)}{2m}t\right)B(\mathbf{K},\mathbf{r},t),$$

$$B(\mathbf{K},\mathbf{r},t)$$

Full text

I Introduction

The production of correlated pairs of particles is an important requirement for probing the foundations of quantum mechanics. For example, in quantum optics correlated pairs of photons were used to demonstrate violation of Bell’s inequalities bellAspect or the Hong-Ou-Mandel effect oum . In atomic physics, pairs of atoms were produced in the process of the four-wave mixing in optical lattice tunableParis , collision-induced deexcitation of the cloud viennaTB , and collision of Bose-Einstein condensates keterle1 ; paryz0 ; paryz1 . These many-body systems can be used to demonstrate sub-Poissonian relative atom number statistics paryz2 , the violation of the Cauchy-Schwarz inequality for matter waves paryz3 ; wasakCSI , atomic Hong-Ou-Mandel effect atomicHOMtheo ; atomicHOM , or ghost-imaging ghost . The nonclassicality of the states of the system can be potentially employed in the violation of Bell’s inequality for atoms atomBell or atomic interferometry njpWasak .

In order to have a high degree of control of the fragile states of the pairs of atoms, it is necessary to posess an accurate model describing the system as well as the processes underlying generation of the correlated pairs. In this paper we analyze elastic scattering of atoms from a pair of colliding Bose-Einstein condensates (BECs). Such condensate collisions were investigated in theoretical papers bach ; zin1 ; zin2 ; zin3 ; zin4 ; deuar1 ; deuar2 ; gardiner1 ; gardiner2 ; paryz10 ; karen0 ; karen1 ; wasakdeuar . As a result of binary collisions between the particles that constitute the counter-propagating clouds, atomic pair scatter out from the condensates with opposite velocities. In the spontaneous regime, where bosonic enhancement does not influence single collision event, the direction of velocity of outgoing particles is random. As a result of superposition principle, the quantum state of single atomic pair is entangled in momentum directions. Such systems, with non-classical correlations between massive particles entangled in external degrees of freedom, are ideal for fundamental experimental tests of quantum mechanics such as Einstein-Podolsky-Rosen gedanken experiment Wieden .

In order to prove the presence of non-classical correlations between scattered atoms, or to utilize the entanglement in quantum-enhanced atomic interferometry, properties of the state of the system need to be well known. In particular, an important quantities are the widths of the correlation volume, which describes the extent of the correlation between the two particles. To find those widths one has to analyze properties of the two particle correlation function. This function was recently measured experimentally paryz0 ; paryz3 , and also was the subject of theoretical studies based on perturbation theory bach ; zin2 ; zin3 ; karen1 or stochastic calculations deuar1 ; deuar2 ; gardiner1 ; gardiner2 ; karen0 .

In the case of perturbation theory, a number of uncontrolled approximations were introduced in order to obtain analytical results. The main approximations were to neglect the interaction between atoms within BECs during the collision process, and omission of the influence of the mother-condensates on the scattered atoms. Therefore, it is important to know how the interactions between atoms change the state of the system, and how they modify the non-classicality of the correlations between the particles. The results are of importance for the planned experiments that might test EPR-correlations, non-locality, or utilize such systems in quantum-enhanced metrology.

Here, we derive, within the perturbation theory, approximate analytical expressions for the two particle correlation function taking into account the contribution from the interaction between atoms. We show that there are two types of correlations: “back to back” – related to the fact that particle are scattered in pairs, and “local” – which describe bosonic bunching effect. The obtained formulas connect correlations functions with colliding condensate wavefunctions being a solution of time-dependent Gross-Pitaevskii (GP) equation. The derived expressions have the form convenient for numerical analysis and as such can be useful for interpretation of the experimental data. Additionally, we derive approximate analytical form of the colliding condensates wavefunctions using time-dependent variational approximation. Using these wavefunctions we obtain analytical expressions for the two-body correlation functions. We analyze properties of the correlation functions among which the most important are the correlation widths for which we present explicit expressions. Finally, we present the properties of the correlation volumes.

The paper is constructed as follows. In Section II we present the system and the methods to describe it based on the Bogoliubov theory. In this Section we derive general expressions for a two particle correlation function which structure can be divided into the two types of correlations. Also, we present the results of variational approximation for static and colliding Bose-Einstein condensates. In Section III we analyze the “back to back” correlations, and present approximate expressions for the pair correlation function. Then using the approximate condensate wave functions given by variational approach, we derive explicit formulas for the pair correlation function and analyze its properties in important cases. We unify the obtained results in explicit expressions, and present the properties of the correlation volume. Additionally, we introduce a semiclassical model and compare its results with the results of the quantum one. In Section IV we analyze the “local” pair correlations, which are directly related to single particle correlation function $G^{(1)}$. We introduce a classically motivated function describing properties of the source of atoms and relate it with the single particle correlation function. Further, we present explicit formula for this function obtained from a quantum model of the process. We calculate $G^{(1)}$ using the approximate condensate wave functions given by variational approach and analyze its properties in certain cases. Also in this section, we unify the obtained results in explicit expressions, and present the properties of the correlation volume. The details of the derivations of most of the formulas as well as the conditions presented in the main body of the paper can be found in the Appendices.

II Theoretical model

In the limit of low energies the collision of bosonic atoms can be described by a single parameter called the scattering length zderzenia . In our considerations we assume that the interaction potential is effectively described by this parameter. Additionally, we assume that the atomic gas is in the dilute limit, i.e., $na^{3}\ll 1$, where $n$ is a maximal density of the colliding clouds and $a$ is the scattering length. Furthermore, we restrict our considerations to the so-called collisionless regime, in which the probability of a secondary collision of the scattered atom with the atoms from the condensate is much smaller than unity. This condition requires that the size $\sigma$ of the cloud in the direction of the atomic velocity is much smaller than the mean free path of the scattered particles $\ell_{\mathrm{mfp}}=1/8\pi a^{2}n$, that is $\sigma/\ell_{\mathrm{mfp}}\ll 1$. Another assumption we use is that the total number of scattered atoms is much smaller than the number of atoms in the moving clouds.

We remark one fact that determines our interests in specific observables. After the collision the system consists of the two condensates and the halo of scattered atoms. Some of the scattered atoms are located on the condensates, because their mean velocity points in the direction of the velocities of the condensates, and both of these vectors are of the same magnitude. However, the analysis of those scattered atoms is difficult experimentally, because the density of the halo is small compared to the density of the condensates. For this reason, in all the considerations below, we restrict our study to the atoms that are scattered away from the collision direction of two condensates (we take it to be $z$ axis). Specifically, we assume that the part ${\bf k}_{r}=k_{x}{\bf e}_{x}+k_{y}{\bf e}_{y}$ of the wavevector ${\bf k}$ of the scattered atom perpendicular to the long axis of the condensate satisfies the following conditions: $|{\bf k}_{r}|/k>1/2$.

II.1 Bogoliubov method

The system that obeys all of the above restrictions is well described by the Bogoliubov method przeglad . In this method moving condensates are described by single particle wave-function $\psi({\bf r},t)$ which satisfy Gross Pitaevskii equation

[TABLE]

where $g=\frac{4\pi\hbar^{2}a}{m}$ parametrizes the interaction strength between atoms. The properties of scattered atoms are described by the field operator $\hat{\delta}({\bf r},t)$ which undergoes time evolution given by

[TABLE]

with:

[TABLE]

We assume that the initial state of noncondensed particles is vacuum dopiska , i.e.,

[TABLE]

In the calculation of mean values of the product of $\hat{\delta}$ operators, taken for various positions, the Wick’s theorem can be applied. This is a consequence of the lineary of the equation of motion, Eq. (2), and the fact that the initial state is vacuum. Therefore, all the correlation functions of arbitrary order decompose into products of anomalous density

[TABLE]

and single particle correlation function

[TABLE]

For example, the two particle correlation function is of the following form:

[TABLE]

In this equation, the first term, given by $G^{(1)}\left({\bf r}_{1},{\bf r}_{1},T\right)G^{(1)}\left({\bf r}_{2},{\bf r}_{2},T\right)$, is a product of single particle densities, and so it represents uncorrelated particles. The presence of second and third terms are responsible for nontrivial corelations between particles.

In the next sections, we show that the terms $\left|M\left({\bf r}_{1},{\bf r}_{2},T\right)\right|^{2}$ and $\left|G^{(1)}\left({\bf r}_{1},{\bf r}_{2},T\right)\right|^{2}$ represent correlations of particles with opposite and collinear velocities, respectively. The appearance of correlation of particles with opposite velocities, called “back to back” or “cross” correlation, is a consequence of the fact that atoms are scattered in pairs of opposite momenta. The correlation of particles with collinear velocities, called “local” correlation, is related to the bosonic bunching effect.

In the considered system a bosonic enhancement effect can take place. This effect was both predicted theoretically keterle1 ; zin1 ; zin2 and observed experimentally keterle1 . In this paper, we restrict to the regime where the effects of the bosonic enhancement are negligible. Then, the Heisenberg equation of motion, Eq. (2), can be approximately solved with help of perturbation theory. In this case, the formula for the anomalous density reads

[TABLE]

where $K({\bf r}_{1},t_{1};{\bf r}_{2},t_{2})$ is a single body propagator of Hamiltonian given in Eq. (3). Furthermore, on the grounds of perturbative approach, the following relation between one body correlation function and anomalous density can be established:

[TABLE]

The details of the derivation of the above formulas are presented in Appendix A. The conditions for the validity of the first order perturbation calculus are derived and discussed in Appendix G.

II.2 Properties of the condensates

Let us now further specify the properties of the considered system. We initially deal with a single condensate described by a wavefunction $\psi({\bf r})$ satisfying stationary GP equation

[TABLE]

where the trapping potential

[TABLE]

and the normalization condition $\int\mbox{d}{\bf r}|\psi({\bf r})|^{2}=N$. Below, we focus on elongated cigar shaped condensate, for which $\omega_{z}\ll\omega_{r}$. The two counterpropagating condensates are created from the stationary one by applying Bragg pulse and switching off the trapping potential keterle1 . After the pulse, the wavefunction takes the following form

[TABLE]

where $C$ is the normalization coefficient. This wavefunction represents two wavepackets, $\psi({\bf r})^{\pm iQz}$, each propagatin with mean velocity $\pm\hbar Q/m$ alogn $z$-axis. The collision takes place along longitudinal axis, the $z$-axis, of the condensate.

The decomposition of the wavefunction into two counterpropagating wavepackes is permissible also for later times. To this end, we assume that the width of two counterpropagating condensates in the momentum space during the collision are much smaller then their mean momenta equal to $\pm\hbar Q$. We additionally assume that the mean field potential $gn$ is much smaller than the kinetic energy $\frac{\hbar^{2}Q^{2}}{2m}$. In Appendix B.2 we show that these two conditions can be replaced by a single one which reads

[TABLE]

where $\sigma_{r}$ denotes the radial width of the initial condensate, and $a_{hor}=\sqrt{\frac{\hbar}{m\omega_{r}}}$ is the radial harmonic oscillator length. The above condition combined with the fact that the system is elongated, implies $Q\gg\frac{1}{\sigma_{z}}$, where $\sigma_{z}$ denotes the longitudinal size of the initial condensate. As a consequence, $C\simeq\sqrt{2}$ which makes Eq. (13) to take the following simple form:

[TABLE]

The assumption the widths of the momentum distribution of both colliding condensates during the whole collision are much smaller than $Q$, leads to two distinct momentum distributions centered around $\pm Q{\bf e}_{z}$. Therefore, it is natural to decompose the condensate wavefunction $\psi({\bf r},t)$ into two wavefunctions, denoted by $\psi_{\pm Q}({\bf r},t)$, in the following way:

[TABLE]

These functions, $\psi_{\pm Q}$, describe the two counterprogating parts of the condensate. Note that, as implied by Eq. (15), initial conditions are $\psi_{\pm Q}({\bf r},0)=\frac{1}{\sqrt{2}}\psi({\bf r})$ for these components.

The stated assumption leads to great simplification of the solution of the GP equation, Eq. (1), on the basis of the slowly varying envelope approximation slowlyvarying . Within this approximation, the GP equation for $\psi({\bf r},t)$ decouples into set of two equations for $\psi_{\pm Q}({\bf r},t)$

[TABLE]

Let us remark, that Eq. (17) is much better for numerical implementation of our problem than the initial GP equation, Eq. (1). The reason is that in Eq. (17) the highly oscillatory behavior in position, due to $e^{\pm iQz}$, as well as in time, due to $e^{-iQ^{2}\hbar t/2m}$, is removed. Consequently, the window in the momentum representation, required for numerical simulation, needs to take the momentum width of the $\psi_{\pm Q}$ alone, and this is much smaller than the momentum window required for the solution of Eq. (1), the latter being of the order of $2Q$.

In order to describe the properties of the wavefunctions $\psi_{\pm Q}$, we use variational method, to solve stationary and time dependent GP equations, Eqs. (11) and (17). The details of the solution are described in Appendix B. Here, we just state the obtained results.

We assume that initially the wavefeunction can be approximated by a gaussian ansatz of the form:

[TABLE]

After time $t$, the wavefunctions evolve into:

[TABLE]

where

[TABLE]

where $v_{0}=\hbar Q/m$, and $\beta=\frac{\sigma_{r}^{2}}{a_{hor}^{2}}\geqslant 1$. Notice, that the final form is also a gaussian function, but with centres moving in opposite directions with velocities $\pm v_{0}\mathbf{e}_{z}$, and with time-dependent widths.

The components of the wavefunctions, given by $\psi_{\pm Q}$ in Eq. (19), can be investigated further in order to determine important timescales in the problem. First, notice that there is a characteristic time $\tau_{c}=\frac{\sigma_{z}}{v_{0}}$ during which the wavepackets cross each other. This can be defined as the time of the collision. The second characteristic time $\tau_{ex}=\omega_{r}^{-1}$ is equal to the time needed for $\sigma_{r}(t)$ to change its width by a factor of $\sqrt{2}$. This two times describe the density properties of the system. Next, characteristic times correspond to changes in the phase of the wavepacket. It is natural to define following three timescales: $\tau_{r}=\frac{1}{\beta\omega_{r}}$, $\tau_{z}=\frac{1}{\omega_{r}}\tan\left(\beta-\frac{1}{\beta}\right)^{-1}$, and $\tau_{\phi}=\frac{1}{\omega_{r}}\tan\left(\frac{7\beta}{4}-\frac{3}{4\beta}\right)^{-1}$. Note, that $\tau_{ex}\geqslant\tau_{r}$, $\tau_{z}\geqslant\tau_{r}$, and $\tau_{\phi}\geqslant\frac{4}{7}\tau_{r}$. Therefore, we see the three times, $\tau_{ex}$, $\tau_{z}$, and $\tau_{r}$ are all larger than $\tau_{r}/2$.

In this section we described the colliding condensates. In the next one, we characterize the properties of the scattered atoms.

III Back to back correlations

III.1 General considerations

Let us now show that the back to back correlations are given by the term $\left|M\left({\bf r}_{1},{\bf r}_{2},T\right)\right|^{2}$. To this end, we analyze the structure of the anomalous density $M\left({\bf r}_{1},{\bf r}_{2},T\right)$. However, notice first that in the experimental situation the atoms are measured at time $T$ which is usually much larger than the time of the scattering process. It is thus permissible to investigate the limit $T\rightarrow\infty$. We introduce new variables: ${\bf r}_{1,2}=\frac{\hbar{\bf k}_{1,2}}{m}T$, and define

[TABLE]

In the above we have introduced additional phase $\exp[-i\hbar(k_{1}^{2}+k_{2}^{2})T/2m]$ to get finite limit and additional factor $(\hbar T/m)^{3}$ to satisfy normalization condition $\int\mbox{d}{\bf k}_{1}\mbox{d}{\bf k}_{2}\,|M({\bf k}_{1},{\bf k}_{2})|^{2}=\int\mbox{d}{\bf r}_{1}\mbox{d}{\bf r}_{2}\,|M({\bf r}_{1},{\bf r}_{2},T)|^{2}$. Let us now continue taking $K$ as free propagator. Upon inserting Eq. (9) into Eq. (21) and using the explicit form of the free propagator we arrtive at

[TABLE]

Substituting now Eq. (16) into Eq. (4) we obtain

[TABLE]

We shall now make use of the assumptions stated in Section II.2. We assumed there that the widths in the momentum space of each of wavepacket $\psi_{\pm Q}$ at any time are much smaller than $Q$ and that the kinetic energy $\frac{\hbar^{2}Q^{2}}{2m}$ is much larger than than the mean field interaction energy $gn$. This assumptions have two consequences. First, the dominant temporal phase is given by factor $\exp(-i\hbar Q^{2}t/2m)$. Second, the dominant spatial phase factor in $\psi_{\pm Q}e^{\pm iQz}$ is $\exp\left(\pm iQz\right)$. Inserting Eq. (23) into Eq. (III.1) and using the above stated assumptions we notice that the temporal integral is vanishingly small unless $k_{1}^{2}+k_{2}^{2}\simeq 2Q^{2}$. As the width of $\psi_{\pm Q}$ in momentum space is much smaller than $Q$ the spatial integral gives nonzero values for the term $\psi_{+Q}\psi_{-Q}$ if $|{\bf k}_{1}+{\bf k}_{2}|\ll Q$, and for the terms $\psi_{\pm Q}^{2}$ if $|{\bf k}_{1}+{\bf k}_{2}\mp 2Q{\bf e}_{z}|\ll Q$. Consequently, the term responsible for the scattering of atoms into the observed halo is due to the term $\psi_{+Q}\psi_{-Q}$. As we are interested only in the properties of atoms appearing in the halo we neglect the contribution from the other two terms, i.e.,

[TABLE]

The above analysis also shows that the anomalous density with $B$ given by Eq. (24) leads to ${\bf k}_{1}$ and ${\bf k}_{2}$ being practically antiparallel with length equal approximately to $Q$. For other choices of ${\bf k}_{1}$ and ${\bf k}_{2}$ the value of anomalous density is vanishing. Thus, as long as it is permissible to exploit free propagator for $K$, we have shown that the back to back correlation is represented by the term $G^{(2)}_{bb}({\bf k}_{1},{\bf k}_{2})\equiv\left|M\left({\bf k}_{1},{\bf k}_{2}\right)\right|^{2}$ in the two particle correlation function in Eq. (8).

The situation with the true propagator, the one which takes into account the interaction of scattered particles with the atoms from the condensates, is much more complicated as there is no analytical formula for $K$. However, under assumptions stated in Section II together with the additional condition

[TABLE]

a semiclassical approximation can be used, which results in

[TABLE]

where

[TABLE]

and ${\bf K}=\frac{{\bf k}_{1}-{\bf k}_{2}}{2}$, ${\bf e}_{{\bf K}}=\frac{{\bf K}}{K}$. The details of the derivation are given in Appendix C.2. There we also show that $|M|^{2}$ represents the back to back correlation, the same result which we obtained for the free propagator case. The expression in Eq. (26) has exactly the same form as the free propagator formula, Eq. (III.1), the only difference being the change from $B({\bf K},{\bf r},t)$ to $B({\bf r},t)$.

The semiclassical approximation leads to the same results as the free propagator approximation if $|\Phi|\ll 1$, for which $B({\bf K},{\bf r},t)\simeq B({\bf r},t)$. In Appendix C.2 we show that the contribution from this phase can be neglected if

[TABLE]

From this condition, it is straightforward to obtain the one given in Eq. (14).

The set of equations: (11), (12), (17), (24), (26), (27), and (28) allow for calculation of the anomalous density $M$. This function includes the mean field effects of the interaction between condensates and scattered atoms, if the condition given in Eq. (25) is satisfied. For a given system configuration the calculations are to be done numerically. However, below we exploit the approximate solution of the colliding condensates wavefunctions, given by Eq. (19), to calculate and analyze the properties of the back to back correlation function. But before going on, we pause to introduce a semiclassical model of the back to back correlation function, which will serve as a probe of quantum characteristics of the collision process.

III.2 Semiclassical model

From a theoretical point of view, it is interesting to compare the results of a quantum model to a classical one. To this end, we consider a semiclassical model of colliding clouds. We describe the two wavepackets by a single-particle phase space densities, denoted by $W_{\pm Q}({\bf r},{\bf k},t)$, and apply all the assumptions stated previously that defined our system, i.e.,

• •

dilute gas limit, so that only two body collisions are of importance,

• •

neglect secondary collision between the scattered atoms and the atoms in the condensates,

• •

neglect depletion of the condensate due to the scattering.

As we compare this model with the quantum one we put ${\bf k}$ vector instead of velocity in the phase space density definition, these are simply related by ${\bf k}=m\mathbf{v}/\hbar$. For this semiclassical model the formula for the back to back part of the second order correlation function takes the following form:

[TABLE]

where the total cross-section $\sigma_{tot}=8\pi a^{2}$, whereas the wavevectors ${\bf K}=\frac{{\bf k}_{1}-{\bf k}_{2}}{2}$ and $\Delta\mathbf{K}={\bf k}_{1}+{\bf k}_{2}$. This formula resembles the production term in the collision integral of the Boltzmann equation, in which apropriate substitutions for the physical quantities are made. Now, we take $W_{\pm Q}({\bf r},{\bf k},t)$ to be the Wigner distribution:

[TABLE]

where we introduced

[TABLE]

The back to back correlation function is expressed, by Eq. (30) and Eq. (31), in terms of the condensate wavefunctions only.

III.3 Examples of $G^{(2)}_{bb}$ function

Let us now calculate the anomalous density for some cases using the approximate analytical form of $\psi_{\pm Q}$ given by the variational method presented in Appendix B. However, in the considerations below in this Section we assume that the condition in Eq. (29) is satisfied, and so $B({\bf k},{\bf r},t)\simeq B({\bf r},t)$. Using Eqs. (19), (20), (24), and (26), performing the spatial integral, and introducing the dimensionless time $\tau=\omega_{r}t$, we arrive at the following formula

[TABLE]

where $\Delta{\bf K}={\bf k}_{1}+{\bf k}_{2}$, $A=\frac{gN}{\hbar(2\pi)^{3}\omega_{r}}$, $\alpha=\frac{\tau_{ex}}{\tau_{c}}=\frac{Qa_{hor}^{2}}{\sigma_{z}}$, $\Delta K_{r}^{2}=\Delta K_{x}^{2}+\Delta K_{y}^{2}$; the dimensionless frequency is

[TABLE]

and the time dependent functions are

[TABLE]

Now, we center $K$ at $Q$, introducing $\delta K=K-Q$, and rewrite $\omega$ in terms of $\delta K$ as $\omega=\left(2Q\delta K\left(1+\frac{\delta K}{2Q}\right)+\frac{\Delta K^{2}}{4}\right)a_{hor}^{2}$. According to general consideration above $\delta K\ll Q$, which leads to approximate form of $\omega$:

[TABLE]

In Eq. (33) the term $\exp\left(-\frac{\Delta K_{z}^{2}\sigma_{z}^{2}}{4c_{z}(\tau)}\right)$ implies that $|\Delta K_{z}\sigma_{z}|$ is maximally of the order of unity. This gives $\frac{1}{4}\Delta K_{z}^{2}a_{hor}^{2}$ present in the above formula to be maximally equal to $\frac{a_{hor}^{2}}{\sigma_{z}^{2}}\ll 1$ and can be neglected. As a result we end up with

[TABLE]

Now the anomalous density given by (33) and (34) is a function of dimensionless variables $\Delta K_{r}\sigma_{r}$, $\Delta K_{z}\sigma_{z}$ and $2Q\delta Ka_{hor}^{2}$ that depends on two dimensionless parameters $\alpha$ and $\beta$. In fact a complete analysis should investigate the anomalous density for all values of $\alpha$ and $\beta$. However, this is in practice impossible. We therefore choose few examples for different values of the parameters.

The measurement of correlation function was already performed for metastable helium atoms in the Palaiseau group paryz0 and is planned to be performed in the Vienna group Wieden . In Appendix E we calculate $\alpha$ and $\beta$ present in the above formulas using the parameters of these experiments. We obtain $\alpha\simeq 0.22$, $\beta\simeq 3.3$ for the Palaiseau, and $\alpha\simeq 0.2$, $\beta\simeq 11$ for the Vienna setup. In the experiment, the number of condensate atoms can be reduced by the use of radio frequency “knife” which results in the increase of $\alpha$ and decrease of $\beta$. Thus, in all the calculations presented below we take $\alpha=0.2$ as the minimal value and $\beta=10$ as the maximal value.

III.3.1 Fast collision

The “fast collision” case realize when the velocity of the wavepackets is large enough so that the change of the wavefunctions $\psi_{\pm Q}$, apart from movement along the collisional axis, is negligible during the whole collision. This means that $\tau_{c}$ is much smaller than all the other characteristic timescales. As $\tau_{r}/2$ is the smallest of the characteristic times (apart from $\tau_{c}$) this condition can be restated as

[TABLE]

In this case the anomalous density reads

[TABLE]

The details of the derivation can be found in Appendix D. Notice first that the width in $\Delta K_{r,z}$ of the $|M({\bf K},\Delta{\bf K})|^{2}$ function is $\sqrt{2}$ larger than the momentum density width of $|\psi_{\pm Q}({\bf k},t)|^{2}\propto\exp\left(-(k_{x}^{2}+k_{y}^{2})\sigma_{r}^{2}-k_{z}^{2}\sigma_{z}^{2}\right)$. Also, the width in $K$, denoted by $\Delta_{K}$ and equal approxmately to $\sigma_{z}^{-1}$, is much smaller than the width in $\Delta K_{r}$, equal approximately to $\sigma_{r}^{-1}$. The analogous semiclassical expression for the $G^{(2)}_{bb}({\bf K},\Delta{\bf K})$ function (calculated also in Appendix D) reads:

[TABLE]

Comparing the above formula with $|M({\bf K},\Delta{\bf K})|^{2}$ given by Eq. (36), we clearly see that the $\Delta{\bf K}$ dependence is the same in both cases. However, the semiclassical and quantum formulas differ in $K$ dependence. In Fig. 1 we plot both the semiclassical dependence $2\exp(-2x^{2})$, where $x=\delta K\sigma_{z}$, and the quantum dependence $|\exp(-x^{2})(1+\mbox{erf}(ix))|^{2}$. We observe that the quantum dependence is wider with respect to the semiclassical one, and has a long tail which is absent in the semiclassical case. What is worth noticing both functions integrated over $\delta K$ give the same result, i.e., $\int\mbox{d}x\,2\exp(-2x^{2})=\int\mbox{d}x\,|\exp(-x^{2})(1+\mbox{erf}(ix))|^{2}$.

III.3.2 Strong radial confinement

Let us now consider another case when the mean-field energy $gn$ is much smaller than the kinetic energy along radial direction, $gn\ll\frac{\hbar^{2}}{2m\sigma_{r}^{2}}$. According to Eq. (99), in such “strong radial confinement” case $\sigma_{r}\simeq a_{hor}$ or, alternatively, $\beta\simeq 1$. Then the anomalous density given by (33) and (34) takes the form:

[TABLE]

where

[TABLE]

For this collision configuration, the semiclassical result reads

[TABLE]

where

[TABLE]

We observe that the dependence on $\Delta{\bf K}$ and $K$ decouples in the quantum as well as in the semiclassical model. Additionally, the $\Delta{\bf K}$ dependence is the same in both models. The widths in $\Delta K_{r}$ and $\Delta K_{z}$ are equal approximately to $a_{hor}^{-1}$ and $\sigma_{z}^{-1}$, respectively. As in the fast collision case, the widths in $\Delta{\bf K}$ of the $|M({\bf K},\Delta{\bf K})|^{2}$ function are $\sqrt{2}$ larger than the momentum density widths, for which $|\psi_{\pm Q}({\bf k},t)|^{2}\propto\exp\left(-(k_{x}^{2}+k_{y}^{2})a_{hor}^{2}-k_{z}^{2}\sigma_{z}^{2}\right)$. The dependence on $\delta K$ of $G^{(2)}_{bb}({\bf K},\Delta{\bf K})=|M({\bf K},\Delta{\bf K})|^{2}$ in both models is determined by $|f_{q}(x,\alpha)|^{2}$ and $f_{cl}(x,\alpha)$ functions (where $x=\delta KQa_{hor}^{2}$), respectively. We note here that, as in the fast collision case, $f_{q}$ and $f_{cl}$ satisfy normalization condition: $\int\!\mbox{d}x\,|f_{q}(x,\alpha)|^{2}=\int\mbox{d}x\,f_{cl}(x,\alpha)$.

The two function $f_{q}$ and $f_{cl}$ depend on parameter $\alpha$, the ratio of the expansion time to the collision time. The case $\alpha\gg 1$ describes the fast collision analyzed above, and, therefore, we focus only on $\alpha<1$. As mentioned previously, the experimentally important minimal value of $\alpha$ is $0.2$. In Fig. 2 (panels a and b), we plot the functions $|f_{q}(x,\alpha)|^{2}$ and $f_{cl}(x,\alpha)$ for two values $\alpha=1/2,1/5$.

As we see from Fig. 2, in strong radial confinement case the differences between the semiclassical and quantum model are similar to that in the fast collision. As before the quantum function is wider than the semiclassical one. However in this case, there is a shift of the position of the maximum to positive values of $\delta K$, which is, however, much smaller than the width in $\delta K$. The curves $|f_{q}(x,0.2)|^{2}$ and $|f_{q}(x,0.5)|^{2}$, as can be seen in Fig. 2 (panel c), approach each other at a certain point. The universal curve to which all $|f(q,\alpha)|^{2}$ converge is $|f_{q}(x,0)|^{2}$. Both curves are almost the same for $x>1.7\alpha$, but for smaller $x$ they start to differ. For small value of $x$, $|f_{q}(x,0)|^{2}\propto\log^{2}(2x)$, and thus tends to infinity as $x\to 0$.

The maximum of $|f(q,\alpha)|^{2}$ grows with its position tending to zero as $\alpha$ gets smaller. Thus, we cannot choose halfwidth as the characteristic width of $|f(q,\alpha)|^{2}$. Instead, we define it as the value of $x_{0}$ for which the normalized integral under the curve is $1/2$, i.e,

[TABLE]

This definition is motivated by the fact that the detectors, on which particles fall, have finite sizes. The measurement of two particle correlation function is always accompanied by integration of that quantity over the size of the detector. As $x\propto\delta K$, which is directly related to the position of the detectors, the measurement results in integration over $x$. Using the above formula, we find $x_{0}=0.26$ for $\alpha=0$ (universal curve), $x_{0}=0.32$ for $\alpha=0.2$ (the smallest realistic value considered in the paper), $x_{0}=0.60$ for $\alpha=0.6$, and $x_{0}=0.86$ for $\alpha=1$. We see that these values are of the same order.

As a consequence, the characteristic width $\Delta_{K}$ for $\alpha<1$ is approximately equal to $1/2Qa_{hor}^{2}$ while the position of the maximum, denoted as $\delta K_{max}$, is always smaller than $1/2Qa_{hor}^{2}$. For $\alpha\rightarrow\infty$, we obtain the fast collision case where the width of $\delta K$ is given by $1/\sigma_{z}$, and this turns out to be true for $\alpha>1$. According to the condition given by Eq. (14) and the fact that $\sigma_{z}\gg\sigma_{r}$ both of this widths are much smaller than $1/\sigma_{r}$. Therefore, for all considered values of $\alpha$, the width in $\delta K$ is much smaller than the width in $\Delta K_{r}$ equals approximately to $1/\sigma_{r}$.

III.3.3 Largest mean-field energy impact

According to Eq. (99) mean field energy divided by characteristic kinetic energy connected with the radial confinement $\frac{\hbar^{2}}{m\sigma_{r}^{2}}$ is proportional to $\beta^{2}-1$. Thus the larger the mean field energy (we mean $gn/(\hbar^{2}/m\sigma_{r}^{2})$) the larger is $\beta$. On the other hand the effective time of integration is given by $1/\alpha$. The mean field energy impact shall be the largest for largest possible value of $\beta$ and smallest possible value of $\alpha$. Therefore we call such case ”largest mean-field energy impact”.

As discussed above, we take specific values of $\alpha=0.2$ and $\beta=10$ as extremal values that are experimentally feasible. The anomalous density given by Eq. (33) is a function of $\omega$, $\Delta K_{r}\sigma_{r}$ and $\Delta K_{z}\sigma_{z}$, and can be written as $|M(\omega,\Delta K_{r}\sigma_{r},\Delta K_{z}\sigma_{z})|^{2}$, that explicitly depends on three independent parameters. In Fig. 3 we plot its cuts $|M(\omega,\Delta K_{r}\sigma_{r},0)|^{2}$ and $|M(\omega,0,\Delta K_{z}\sigma_{z})|^{2}$. We observe that the maximum is for $\omega=25$ and $\Delta K_{r}=\Delta K_{z}=0$. The characteristic width in $\Delta K_{z}$ and $\Delta K_{r}$ is approximately $2\sigma_{z}^{-1}$ and $2\sigma_{r}^{-1}$, respectively. The characteristic halfwidth in $\omega$ is roughly equal to $6$. With this values, the term $\Delta K_{r}^{2}a_{hor}^{2}/4$ in Eq. (34) for $\omega$ can be estimated to be maximally

[TABLE]

As $\beta=10$ this term is much smaller than unity and, as the width in $\omega$ is $6$, can be neglected. This results in $\omega\simeq 2\delta KQa_{hor}^{2}$, and thus in the halfwidth in variable $\delta K$ equal to ${3}/{Qa_{hor}^{2}}$. According to the condition Eq. (14), it is much smaller than $1/\sigma_{r}$. Note, that the shift of the maximum in $\omega$ is larger than the width in $\omega$.

III.3.4 The limiting case of $\alpha\to 0$ and $\beta\to\infty$

Here, we analyze the limiting case for which $\beta\to\infty$ and $\alpha\to 0$. The most important contribution to the integral in Eq. (33) comes from the times $\tau\ll 1$. The anomalous density, given by Eq. (33), takes then the following form

[TABLE]

in which we changed the variable from $\beta\tau$ to $\tau$.

The anomalous density $M$ in this situation is a function of $\tilde{\omega}=\left(\frac{\omega}{\beta}-\frac{7}{2}\right)$ and $\Delta^{2}=\Delta K_{r}^{2}\sigma_{r}^{2}+\Delta K_{z}^{2}\sigma_{z}^{2}$, so the number of parameters were significantly reduced as compared to the general formula given in Eq. (33).

In Fig. 4 we plot $|M(\tilde{\omega},\Delta)|^{2}$, from which we observe that the width in $\Delta$ at the maximum value of $\tilde{\omega}\simeq 0$ is equal to $2$. The shape of the function in $\tilde{\omega}$ is highly asymmetric. For $\tilde{\omega}>0$ the function rapidly decays, while for $\tilde{\omega}<0$ the decay is much slower, with a width of the tail approximately equal to unity.

Remembering that

[TABLE]

we notice that the maximum of $\omega$ is shifted to $7\beta/2$, and the width approximately equals $\beta$. In the same way as in largest mean-field energy impact case, we estimate in Eq. (34) the term ${\Delta K_{r}^{2}a_{hor}^{2}}/{4}\simeq{1}/{\beta}$. As we see, it can be neglected yielding

[TABLE]

Substituting these results into Eq. (39), we obtain the position of the maximum, $\delta K_{max}$, and the width $\Delta_{K}$:

[TABLE]

Due to condition (29), both $\delta K_{max}$ and $\Delta_{K}$ are much smaller than the width in $\Delta K_{r}\simeq 2\sigma_{r}^{-1}$.

Finally, now we comment on the results ofthe previous example for which $\beta=10$. The shift of the maximum was equal to $25$, which is smaller than $7\beta/2=35$. The width in $\omega$ was found to be equal to $6$, which is quite close to $10$ as predicted above. The small discrepancies are caused by the fact that $\beta$ is still not large enough.

III.4 General correlation properties

Above we have shown few examples of the anomalous density. We saw that in the fast collision and strong radial confinement cases the $K$ and $\Delta{\bf K}$ dependence decoupled. In the tw other cases as seen in the figures the $K$ and $\Delta{\bf K}$ dependence almost decouple. the general observation (although not strict) is that for constant $\beta$ this dependence tends to decouple better while enlarging $\alpha$. It can be see even in the fast collision case when $\alpha\gg\beta$ and the dependence is fully decoupled. The additional observation is that the width in $\Delta{\bf K}$ gets smaller with enlargement of $\alpha$. This change is the largest for $\beta\gg 1$ when in the case the case of small $\alpha$ the width in $\Delta K_{r,z}$ is roughly $2/\sigma_{r,z}$ and goes down to $1/\sigma_{r,z}$ for $\alpha\gg\beta$ which is the fast collision case.

III.4.1 Properties of averaged function $\int\!\!\mbox{d}{\bf K}\ G^{(2)}_{bb}({\bf K},\Delta{\bf K})$

The fact that in the two of the analyzed examples the integrals over $\delta K$ gave the same result for both quantum and semiclassical model suggests that it may be a general property. This is indeed the case and in Appendix D we prove that $\int\mbox{d}{\bf K}\ G^{(2)}_{bb}({\bf K},\Delta{\bf K})$ is the same for both classical and quantum models as long as $B({\bf k},{\bf r},t)\simeq B({\bf r},t)$.

Now, let us investigate the properties of the averaged back to back correlations. In Appendix D we show that under condition stated in Section II it takes the following form

[TABLE]

where $\Omega_{\bf K}$ denotes the solid angle coordinates of ${\bf K}$ vector. Note, that the width in $\Delta{\bf K}$ is given by the width of the $B({\bf K},{\bf r},t)$ averaged over time and the directions of ${\bf K}$. Let us now turn to the regime for which $B({\bf k},{\bf r},t)\simeq B({\bf r},t)$. Emploing Eq. (24) we find that

[TABLE]

From this formula, we see that the width of $\int\mbox{d}{\bf K}\,|M({\bf K},\Delta{\bf K})|^{2}$ in $\Delta{\bf K}$ is directly related to the momentum width of $\psi_{+Q}\psi_{-Q}$ averaged over time. This relation can be evaluated exactly in gaussian ansatz the case, for which we have Eq. (19). In Appendix D we prove that

[TABLE]

where $\psi({\bf k},t)$ is the Fourier transform of the wavefunction, $\int\mbox{d}{\bf r}\,e^{-i{\bf k}{\bf r}}\psi({\bf r},t)$. The width of averaged $|M|^{2}$ is increased by a factor of $\sqrt{2}$ with respect to $\psi_{\pm Q}$ momentum density. The same effect is visible in the two of the above calculated examples: the fast collision and the strong axial confinement cases.

Let us now analyze the shape of the averaged function. After the integration over ${\bf K}$ it depends on $\Delta K_{r}$ and $\Delta K_{z}$. To obtain dependence on a single parameter only, we perform additional integration over variable $\Delta K_{z}$. In Appendix D we show that the final functions takes the following form

[TABLE]

The evaluation of this integral requires addition (over time domain) of gaussian functions $\exp(-{\Delta K_{r}^{2}\sigma_{r}^{2}}/{2w^{2}(\tau)})$ with a time depended width $w(\tau)=\surd[{{(1+\tau^{2})}/{(1+\tau^{2}/\beta^{2})}}]$ and weight $\exp(-2{\alpha^{2}\tau^{2}}/{\beta^{2}})/(1+\tau^{2})$. Below, we investigate this function in more details.

Note first, that in the case $\beta\simeq 1$, the width $w(\tau)\simeq 1$ and the averaging of the correlation function yields $\exp(-{\Delta K_{r}^{2}\sigma_{r}^{2}}/{2})$. This result can be seen in the calculation of the anomalous density in the fast collision and strong radial confinement cases. The obtained gaussian function is the narrowest one of all the possibilities. In all the other cases $\beta>1$ and the width $w(t)$ grows in time. The widest possible functions in effective variable $\frac{\Delta K_{r}^{2}\sigma_{r}^{2}}{2}$ is attained in the limit $\beta\rightarrow\infty$, and reads

[TABLE]

where $F(x)=(-i{\sqrt{\pi}}/{2})\exp(-x^{2})\mbox{erf}(ix)$ is the Dawson integral.

In Fig. 5 we plot the function given in Eq. (44), together with the one in Eq. (43) for the narrowest function case, $\beta=1$. We notice that the width of the function given by (44) is about $3/2$ larger than that of gaussian function. The large difference between these functions results from the tail of the Dawson integral. In Appendix B.2 we show that the axial momentum width of $\psi_{\pm Q}$ grows in time from $\frac{1}{\sigma_{r}}$ to $\frac{1}{\sigma_{r}}\frac{\sigma_{r}^{2}}{a_{hor}^{2}}=\frac{\beta}{\sigma_{r}}$. When $\beta\gg 1$, the momentum width grows substantially. As seen in expression from Eq. (42) the shape of $\int\!\mbox{d}{\bf K}\,|M({\bf K},\Delta{\bf K})|^{2}$ is directly related with the momentum density of the wavefunctions $\psi_{\pm Q}$ integrated over time. The long tail seen in Fig. 5 results from the addition of the momentum densities with the width growing in time.

III.4.2 Properties of $G^{(2)}_{bb}(K,\Delta{\bf K})$ in variable $K$

As we have seen above in all of the analyzed examples, the width in $\delta K$ is much smaller than the width in $\Delta K_{r}\simeq\sigma_{r}^{-1}$. Here, we give a simple explanation of this fact by analyzing Eqs. (III.1) and (24).

The temporal dependence of the integrand in these equations is given by

[TABLE]

where $\phi(t)$ is the spatially independent phase of the $\psi_{\pm Q}$ functions approximately given by $\phi(t)\approx\frac{\mu}{\hbar}t$, where $\mu$ is the chemical potential. As analyzed above $|\delta K|\ll Q$ so the term $\delta K^{2}$ can be neglected. Additionally, we take $\Delta{\bf K}=0$. The the shift of the maximum $\delta K_{max}$ is equal to $\frac{\hbar}{m}2Q\delta K_{max}-2\frac{\mu}{\hbar}=0$, and the width $\Delta_{K}$ can be estimated by setting $\frac{\hbar}{m}2Q\Delta_{K}\tau_{d}=1$, where $\tau_{d}$ is some characteristic time of the process. Consequently, we arrive at

[TABLE]

The widest possible $\delta K$ is reached for the shortest $\tau_{d}$. The characteristic times described in the Section II provide the shortest timescale represented by

[TABLE]

Substituting the above in Eq. (45), we arrive at

[TABLE]

where we made use of Eq. (100).

Let us see, whether these formulas agree with the examples presented above remembering that they are only rough estimation of true values. In the fast collision case the approximate formula (47) predicts nonzero value of $\delta K_{max}$ while Eq. (36) gives $\delta K_{max}=0$. However due to condition given by Eq. (35), $\delta K_{max}\ll{1}/{\sigma_{z}}$, which is the $\delta K$ width. This means that within approximations undertaken in the paper both formulas give the same value. The width in $\delta K$ predicted by (48) is $\Delta_{K}={1}/{2\sigma_{z}}$ while the value given by (36) is about twice larger.

In the strong radial confinement case the above formulas give $\delta K_{max}={1}/{Qa_{hor}^{2}}$ and $\Delta_{K}={1}/{Qa_{hor}^{2}}$. While comparing them with true values we notice that the scaling $\frac{1}{Qa_{hor}^{2}}$ is correct. The difference is the prefactors which in fact are smaller than predicted by Eq. (47) and (48). The same situation happens in the third example where the width predicted by Eq. (48) is about twice larger then the true value while the true shift is about 30 percent smaller than given by Eq. (47). In the last example, both the shift and width are correctly predicted by expressions in Eq. (47) and (48). To conclude, the above formulas are in fair agreement with all analyzed examples.

We remark, that the simple derivation of the shift $\delta K_{max}$ can be also interpreted as energy conservation law during the collision of two particles. Specifically, we have two particles with incoming kinetic energy ${\hbar^{2}Q^{2}}/{2m}$, they collide and leave the condensate. During the process, each particle gain additional energy equal to the chemical potential. Finally, their kinetic energies (due to the assumption $\Delta{\bf K}=0$ they are the same for both particles) are ${\hbar^{2}(Q+\delta K_{max})^{2}}/{2m}$. The energy conservation law requires

[TABLE]

which, after omission of the term $\delta K_{max}^{2}$, leads to Eq. (47) as expected.

The width in $\Delta{\bf K}$ is given by the spatial Fourier transform, the width in axial direction is $\sigma_{r}$. Therefore, the minimal width in $\Delta K_{r}$ is ${1}/{\sigma_{r}}$. According to Eq. (29), both $\delta K_{max}$ and $\Delta_{K}$, given by Eqs. (47) and (48), are much smaller than the width in $\Delta K_{r}$.

The above analysis was based on quantum considerations. It is instructive to present a simple classical reasoning, though. We denote by $\sigma_{kr,kz}$ the widths of the momentum densities of $\psi_{\pm Q}$. The momenta of the two atoms, before the collision, may be written as

[TABLE]

where the coefficients $|c|\leqslant 1$. We define ${\bf K}^{\prime}=\frac{{\bf k}_{1}^{\prime}-{\bf k}_{2}^{\prime}}{2}$ and $\Delta{\bf K}^{\prime}={\bf k}_{1}^{\prime}+{\bf k}_{2}^{\prime}$. The energy and momentum conservation laws require

[TABLE]

Substituting ${\bf k}_{1}^{\prime}$ and ${\bf k}_{2}^{\prime}$ into the above equations, we obtain

[TABLE]

where the new coefficients are $C_{j}=c_{1,j}+c_{2,j}$ and $\Delta c_{j}=c_{1,j}-c_{2,j}$, with $j=x,y,z$. The second of the equations can be rewritten as

[TABLE]

where $\delta K=K-Q$. Neglecting $\delta K^{2}$ term we obtain

[TABLE]

The constraints $|\Delta c_{j}|\leqslant 2$ restrict the values of $\delta K$, yielding the width in $\delta K$ to be approximately equal to

[TABLE]

As $Q\gg\sigma_{kr}$ and $\sigma_{kr}\gg\sigma_{kz}$, the with is much smaller than the width in $\Delta K_{r}\simeq 2\sigma_{kr}$. Thus, as in the quantum model the width in $\delta K$ is much smaller than the width in $\Delta K_{r}$.

III.4.3 Correlation volume

Here we show that the fact $\Delta_{K}$ is much smaller than the width in $\Delta K_{r}$ has crucial consequences for the particle correlation properties, i.e., for the correlation volume. To this end, we analyze measurement of two atoms, one at ${\bf k}$ and another at ${\bf k}^{\prime}$. We define the correlation volume, as the volume for which ${\bf k}^{\prime}$, with ${\bf k}$ being constant, the particles are still significantly correlated. In order to calculate this quantity, we first find ${\bf k}_{0}^{\prime}$ for which $G^{(2)}_{bb}({\bf K}_{0},\Delta{\bf K}_{0})$, where ${\bf K}_{0}=\frac{{\bf k}-{\bf k}_{0}^{\prime}}{2}$, $\Delta{\bf K}={\bf k}+{\bf k}_{0}^{\prime}$, takes maximum value. By changing the variables to $\delta{\bf k}^{\prime}={\bf k}^{\prime}-{\bf k}_{0}^{\prime}$, from which we have

[TABLE]

we arrive at

[TABLE]

where ${\bf e}_{{\bf K}_{0}}=\frac{{\bf K}_{0}}{K_{0}}$. As found above, the widths in $\Delta{\bf K}$ are approximately equal to ${1}/{\sigma_{r}}$ and ${1}/{\sigma_{z}}$ in the axial and longitudinal direction, respectively. Thus, from Eq. (49) we find the same for $\delta{\bf k}^{\prime}$. On the other hand, from Eq. (50) we find that the width in ${\bf e}_{{\bf K}_{0}}\delta{\bf k}^{\prime}$ is approximately equal to $2\Delta_{K}$. As ${\bf e}_{{\bf K}_{0}}\simeq{\bf e}_{\bf k}$, we find that the width in $\delta k^{\prime}$ in a given direction ${\bf e}_{\delta{\bf k}^{\prime}}=(\cos\tilde{\phi}\sin\tilde{\theta},\sin\tilde{\phi}\sin\tilde{\theta},\cos\tilde{\theta})$ can be estimated as

[TABLE]

where ${\bf e}_{\bf k}=(\sin\theta,0,\cos\theta)$.

Now, let us briefly analyze Eq. (51). As ${1}/{\sigma_{r}}\gg\Delta_{K}$ the largest width $\frac{1}{\sigma_{r}}$ is possible only for a small region around $\cos\tilde{\theta}=\cos\tilde{\phi}=0$. In the remaining area, we have a competition between second and third term of the above formula. The third term is smaller than the second as long as $2\Delta_{K}\sigma_{z}\frac{|\cos\tilde{\theta}|}{{\bf e}_{\bf k}{\bf e}_{\delta{\bf k}^{\prime}}}<1$. Thus, the value $2\Delta_{K}\sigma_{z}$, which is not smaller than unity, cf. Eq. (48), defines the range of angles where the above inequality holds. But, independently of this value, for $\cos\tilde{\theta}=0$ and ${\bf e}_{\bf k}{\bf e}_{\delta{\bf k}^{\prime}}\neq 0$ the inequality is satisfied. For example, for ${\bf e}_{\bf k}={\bf e}_{\delta{\bf k}^{\prime}}={\bf e}_{x}$ ) the width in $\delta k^{\prime}=2\Delta_{K}$. In Section IV.2 we show that in the radial direction the minimal width of the single particle density equals approximately $1/\sigma_{r}$. Therefore, the correlation width along $x$-axis is much smaller than the density width.

IV Local correlations and single particle density

IV.1 General considerations

In this section we show that the term $\left|G^{(1)}\left({\bf r}_{1},{\bf r}_{2},T\right)\right|^{2}$ in the two particle correlation function, Eq. (8), represents collinear correlations of the particles with aligned velocities.

To this end, we investigate the limit $T\rightarrow\infty$ and, along the lines of the study of the back to back correlations, we define:

[TABLE]

Now, we insert Eq. (10) and (21) in Eq. (52) arriving at

[TABLE]

This formula is the central subject of the analysis in this section.

According to previous considerations, the anomalous density $M({\bf k}_{1},{\bf k}_{2})$ is nonzero only in the region where ${\bf k}_{1}$ and ${\bf k}_{2}$ are practically antiparralel with length approximately equal to $Q$. The Eq. (53) implies that ${\bf k}_{1}$ and ${\bf k}_{2}$ are practically parallel, with the length approximately equal to $Q$. This shows that the correlations represented by $G^{(2)}_{loc}({\bf k}_{1},{\bf k}_{2})=\left|G^{(1)}({\bf k}_{1},{\bf k}_{2})\right|^{2}$ are the “local” ones, i.e., nonzero only if ${\bf k}_{1}\approx{\bf k}_{2}$. As a consequence, the width of the scattered atoms halo, described by the single particle density $G^{(1)}({\bf k}_{1},{\bf k}_{1})$, is much smaller than its radius being close to $Q$.

Below we derive an approximate formula for the first order correlation function $G^{(1)}$ that is well suited for numerical treatment. First, however, it is convenient to introduce into the quantum problem the quantity that in the semiclassical limit describes the source of atoms. This object, denoted by $f\left({\bf r},{\bf k},t\right)$, characterizes the distribution of atomic momenta $\hbar{\bf k}$ at every point in space ${\bf r}$ and time $t$ emitted by the source. The function $f$ attains its semiclassical meaning through the relation with the single particle Wigner function $W({\bf r},{\bf k},T)$ of atoms emitted by the source:

[TABLE]

This formula is classical in a sense that it assumes the particles to travel with a velocity $\hbar{\bf k}/m$. In Appendix F.1 we show that $G^{(1)}$ can be expressed in terms of the source $f$ function, by the following formula:

[TABLE]

This equation is useful, because it is the source $f$ that is well suited for various approximations. Below, taking this equation as a starting point, we exploit specific properties of the system to simplify the formula for $G^{(1)}$.

First, we shall approximate $k\simeq Q$ in Eq. (55), and find

[TABLE]

where we renamed the variables $G^{(1)}\left({\bf k},\Delta{\bf k}\right)=G^{(1)}\left({\bf k}+\frac{\Delta{\bf k}}{2},{\bf k}-\frac{\Delta{\bf k}}{2}\right)$, which should not lead to confusion. In Appendix F.2 we show that, as long as Eq. (29) is satisfied, the source function $f$ takes the approximate form

[TABLE]

where $\delta k=k-Q$, and $\sigma_{tot}=8\pi a^{2}$ is the total cross section for collisions of two atoms.

If the source function is a semiclassical quantity, one would expect to find it from classical considerations. In Appendix F.3 we show that under the condition

[TABLE]

this turns out to be true, i.e., the expression in Eq. (57) can be derived from semiclassical model presented in Section III. Furthermore, in Appendix F.2 we show that in a such case the formula for the source function can be simplified to

[TABLE]

where we took ${\bf k}=k(\sin\theta,0,\cos\theta)$ without the lost of generality. Now, we insert Eq. (59) into Eq. (56), use the identity

[TABLE]

where $\delta{\bf k}=(\delta k,\delta k_{y},\delta k_{z})$, introduce ${\bf r}_{1}={\bf r}+\delta{\bf r}$, ${\bf r}_{2}={\bf r}-\delta{\bf r}$, and, finally, we obtain

[TABLE]

This expression be be further simplified, when we substitute $\tilde{\Psi}({\bf k},t)=\int\mbox{d}{\bf r}\,\exp\left(-i{\bf k}{\bf r}\right)\psi_{+Q}\psi_{-Q}({\bf r},t)$. As a result, we obtain

[TABLE]

From this equation, we obtain particularly simple form of the single particle density, $\rho({\bf k})=G^{(1)}\left({\bf k},0\right)$, which reads

[TABLE]

The above formulae are well suited for numerical computation using FFT routines.

Finally, let us notice an important property of single particle density given by (56) and (59). It has the functional form $\rho(\delta k)=\frac{1}{\sin\theta}g\left(\frac{\delta k}{\sin\theta}\right)$, where $g$ is a function. It follows that upon integration over radial variable $k$ we obtain the spherical angle density,

[TABLE]

that is an $\theta$ angle independent value. This is caused by the fact that atoms scatter only in the s-wave, which has an angle independent differential cross section.

IV.2 Single particle density and number of scattered atoms

The formulae derived above can be used to calculate the total number of scattered atoms, and to investigate the properties of the density of the scattered atoms.

First, let us concentrate on total number of scattered atoms, $N_{sc}=\int\mbox{d}{\bf k}\,\rho({\bf k})$. The source $f$ function depends only on $\delta k$, and is nonvanishing only if $|\delta k|\ll Q$, so we approximate $\int\!\mbox{d}{\bf k}\simeq 4\pi Q^{2}\int_{-\infty}^{\infty}\mbox{d}\delta k$. From Eqs. (56) and (57) we obtain that

[TABLE]

This result shows the production rate of scattered atoms is directly proportional to the product of the densities of the counterpropagating clouds. In the case of a gaussian ansatz, see Eq. (19), the number of scattered atoms is

[TABLE]

Therefore, the number of atoms scattered is an integral of a time dependent rate of production of atoms, which vanishes for times larger than min$(2\tau_{ex},\tau_{c})$. This timescale can thus be interpreted as an effective production time, during which particles are scattered from the colliding clouds.

Let us now analyze the density $\rho({\bf k})$ of scattered atoms in the case of a gaussian ansatz. Inserting formula from Eq. (19) into Eq. (60), and performing the gaussian integral we obtain

[TABLE]

where $\delta k=k-Q$, $\tau=\omega_{r}t$ and, due to cylindrical symmetry, we took specific ${\bf k}=k(\sin\theta,0,\cos\theta)$. As seen from the formula, the integral above is a function of $\frac{\delta k\sigma_{r}}{|\sin\theta|}$ and parameters $\alpha$ and $\beta$. In the fast collision $\alpha\rightarrow\infty$ or strong axial confinement $\beta=1$ cases, the density takes the form $\rho({\bf k})\propto\frac{1}{\sin\theta}\exp\left(-2\frac{\delta k^{2}\sigma_{r}^{2}}{\sin^{2}\theta}\right)$. In all the other cases the above formula shows that we sum gaussian functions $d(\tau)\exp\left(-2\frac{x^{2}}{w^{2}(\tau)}\right)$ with the widths, $w^{2}(\tau)={(1+\beta^{2}\tau^{2})}/{(1+\tau^{2})}$, growing in time, and weights, $d(\tau)=[(1+\tau^{2})(1+\beta^{2}\tau^{2})]^{1/2}\exp\left(-2\alpha^{2}\tau^{2}\right)$, decreasing in time.

Therefore, the density in the fast collision and strong radial confinement cases takes the possible narrowest shape and the the widest density is for largest possible value of $\beta$ and smallest possible value of $\alpha$. As discussed in Section III this is the case for which $\alpha=0.2$ and $\beta=10$. The density $\rho({\bf k})$ normalized to its maximal value is a function of $x=\frac{\delta k\sigma_{r}}{\sin\theta}$, and is plotted in Fig. 6 for these values of $\alpha$ and $\beta$. In the plot, we show additionally the narrowest possible case. As we see for $\beta=10$ the density distribution has a long tail coming from the latest $\tau$, where $w(\tau)$ is the largest. However, the halfwidth changed only roughly twice with respect to the $\beta=1$ case.

IV.3 Examples of $G^{(2)}_{loc}$ function

Let us now analyze the properties of the source function $f$ and $G^{(1)}$ in the case of gaussian ansatz given by Eq. (19). Inserting the gaussian ansatz into Eq. (59) and performing the integral one arrives at

[TABLE]

where $C_{f}=\frac{\sigma_{tot}\hbar^{2}N^{2}}{2^{3/2}\pi^{9/2}m^{2}v_{0}\sigma_{z}^{2}\sigma_{r}^{3}}$ and, without lost of generality, we took ${\bf k}=k(\sin\theta,0,\cos\theta)$. Note that $\frac{\delta k\sigma_{r}(t)}{\sin\theta}$ dependence is given by the gaussian with a shift of the maximum equal to $\beta\omega_{r}t{x}/{\sigma_{r}(t)}$.

With the source function $f$ at hand, we can now calculate single particle correlation function $G^{(1)}$. To this end, we insert $f$ from Eq. (63) into Eq. (56), and, after performing the spatial integral, we finally obtain

[TABLE]

where $D_{f}=\frac{(Naa_{hor})^{2}}{\pi^{2}Q\sigma_{z}\sigma_{r}}$, and

[TABLE]

with $\Delta{\bf k}{\bf e}_{\bf k}=\Delta k_{x}\sin\theta+\Delta k_{z}\cos\theta$.

From the plot of the density shown in Fig. 6 in the case of $\beta=10$ and $\alpha=0.2$, we see its halwidth is reached for $\left|\frac{\delta k\sigma_{r}}{\sin\theta}\right|\simeq 2$. As we are not interested in the structure of the tails of the $G^{(1)}$, and, therefore, we shall restrict our considerations to the region where $\left|\frac{\delta k\sigma_{r}}{\sin\theta}\right|\leqslant 1$, which is a bulk region of high density for all values of $\beta$ and $\alpha$.

Let us now analyze Eq. (64) together with Eq. (65) in the bulk density region. Note first that the condition in Eq. (29) together with $\left|\frac{\delta k\sigma_{r}}{\sin\theta}\right|\leqslant 1$ and $\sin\theta>\frac{1}{2}$ implies that $|Qa_{hor}^{2}\sin\theta|\gg\left|\beta\frac{\delta k\sigma_{r}^{2}}{\sin\theta}\right|$. Consequently, as long as $|{\bf e}_{\Delta{\bf k}}{\bf e}_{\bf k}|\gg\frac{\sigma_{r}^{3}}{Qa_{hor}^{4}}$ (which according to Eq. (29) is much smaller than unity) we can neglect the term $\beta\frac{\Delta k_{x}\delta k\sigma_{r}^{2}}{|\sin\theta|}\frac{1+\tau^{2}}{1+\beta^{2}\tau^{2}}$ in $\omega$. Other values of $|{\bf e}_{\Delta{\bf k}}{\bf e}_{\bf k}|$ are much smaller than unity. Therefore, if we set ${\bf e}_{\bf k}=(\sin\theta,0,\cos\theta)$, we can approximate ${\bf e}_{\Delta{\bf k}}\simeq(\cos\theta\cos\tilde{\varphi},\sin\tilde{\varphi},-\sin\theta\cos\tilde{\varphi})$. We note that in such a case the width in $\Delta k$, given by the exponent $\exp\left(-{\Delta k_{z}^{2}\sigma_{z}^{2}}/{8}\right)$ term, is equal to ${1}/{(\sigma_{z}\sin\theta\cos\tilde{\varphi})}$. Therefore, we have the following equality

[TABLE]

The condition in Eq. (58) together with $\left|\frac{\delta k\sigma_{r}}{\sin\theta}\right|\leqslant 1$ results in the following inequality: $\left|\beta\frac{\Delta k_{x}\delta k\sigma_{r}^{2}}{\sin\theta}\frac{1+\tau^{2}}{1+\beta^{2}\tau^{2}}\right|\ll 1$. The largest characteristic time $\tau$ is equal to unity. Thus, the term $\beta\frac{\Delta k_{x}\delta k\sigma_{r}^{2}}{\sin\theta}\frac{1+\tau^{2}}{1+\beta^{2}\tau^{2}}\tau$ can be neglected. We have thus shown that for any direction of the vectors ${\bf e}_{\Delta{\bf k}}$ and ${\bf e}_{{\bf k}}$ we can approximate

[TABLE]

In what follows, we show that ${\Delta k_{x}^{2}\sigma_{r}^{2}}/{8}\ll 1$, which further simplifies the single particle correlation function. We take ${\bf e}_{\Delta{\bf k}}=(\cos\tilde{\theta}\cos\tilde{\varphi},\sin\tilde{\varphi},-\sin\tilde{\theta}\cos\tilde{\varphi})$. The width in $\Delta k$ given by the the exponent in $\exp(-{\Delta k_{z}^{2}\sigma_{z}^{2}}/{8})$ term is equal to ${1}/{(\sigma_{z}\sin\tilde{\theta}\cos\tilde{\varphi})}$, which yields

[TABLE]

This is much smaller than unity for $\tilde{\theta}\gg\frac{\sigma_{r}}{\sigma_{z}}$, which, for $\sigma_{r}\ll\sigma_{z}$, is satisfied for almost all $\tilde{\theta}$ except an excluded region around $\tilde{\theta}=0$. There we note that $|\Delta k_{z}|=|\sin\tilde{\theta}\cos\tilde{\varphi}|\ll|\cos\tilde{\theta}\cos\tilde{\varphi}|=|\Delta k_{x}|$, and so we approximate $\omega\simeq Qa_{hor}^{2}\Delta k_{x}\sin\theta$. By inspecting expression in Eq. (64), we notice that the minimal characteristic time $\tau$ is approximately equal to $\mbox{min}\!\left({1}/{2\alpha},{1}/{\beta}\right)$. Consequently, the width in $\omega$ is equal to $\mbox{max}(2\alpha,\beta)$. As a result, the maximal width in $\Delta k$ can be estimated as $\frac{1}{\sin\theta}\left(\frac{\sigma_{r}^{2}}{Qa_{hor}^{4}},\frac{2}{\sigma_{z}}\right)$. According to the condition in Eq. (29), this is much smaller than ${1}/{\sigma_{r}}$, which results in $\Delta k_{x}^{2}\sigma_{r}^{2}\ll 1$. Finally, the single particle correlation function takes the following form

[TABLE]

The formula for the single particle correlation function in Eq. (68) can be further simplified. To this end, we note that from Eq. (64) the term $(1+\tau^{2})$ in $\exp\left(-{\Delta k_{y}^{2}\sigma_{r}^{2}(1+\tau^{2})}/{8}\right)$ may effectively vary between $1$ and $5$, always giving the width in $\Delta k$ being of the order of ${1}/{\sigma_{r}}$. Then, we may approximate $\Delta k_{y}^{2}\sigma_{r}^{2}(1+\tau^{2})\simeq\Delta k_{y}^{2}\sigma_{r}^{2}$. Furthermore, to restore the cylindrical symmetry of the correlation function, we multiply the right hand side in Eq. (68) by the factor $\exp\left(-{\Delta k_{x}^{2}\sigma_{r}^{2}}/{8}\right)\simeq 1$, Therefore, the correlation function can be written as a product of two factors:

[TABLE]

where the two functions can be conveniently defined according to:

[TABLE]

In the paragraph preceding Eq. (68), we estimated the width in $\Delta{\bf k}{\bf e}_{\bf k}$ to be maximally given by

[TABLE]

Below, we analyze few examples to see how to apply the obtained formulas, and investigate the width of the correlation function.

IV.3.1 Fast collision and strong radial confinement

In the case of fast collision and strong axial confinement case, the function $G_{1}$ decouples into a product of two factors:

[TABLE]

where the density is given by the expression in Eq. (IV.2), which in this case is:

[TABLE]

In the fast collision case, we have

[TABLE]

In Fig. 7 we plot $|G_{3}|$ as an universal function of variable $\Delta{\bf k}\cdot{\bf e}_{\bf k}\sigma_{z}$ . From the figure, we observe that the halfwidth in $\Delta{\bf k}\cdot{\bf e}_{\bf k}$ is equal to ${4.2}/{\sigma_{z}}$.

In the strong radial confinement configuration, the function $G_{3}$ takes the following form:

[TABLE]

If $\alpha\gg 1$, the function reduces to the one in the fast collision case. In the other limiting case, if $\alpha\ll 1$ the term $2\alpha^{2}\tau^{2}$ can be neglected resulting in a simple integral. In Fig. 7 we plot that function, i.e. $|G_{3}|$, normalized to unity at the maximum given by Eq. (74) with the term $2\alpha^{2}\tau^{2}$ dropped. From the figure, we observe that the halfwidth in $\Delta{\bf k}\cdot{\bf e}_{\bf k}$ is equal to ${1.2}/{Qa_{hor}^{2}}$.

IV.3.2 Largest mean-field energy impact considered

As in the back to back correlation we consider the largest mean-field energy impact case $\alpha=0.2$, $\beta=10$. It is additionally motived by the fact that in such case the density takes the widest possible form as found in Subsection IV.2. In Fig. 8 we plot $\left|G_{1}\left(\frac{\delta k}{\sin\theta},\Delta{\bf k}\,{\bf e}_{\bf k}\right)\right|$ for $\left|\frac{\delta k\sigma_{r}}{\sin\theta}\right|<1$. From the figure we obseve that the halfwidth in $Qa_{hor}^{2}\Delta{\bf k}\,{\bf e}_{\bf k}$ depends monotonically on the value of $\frac{\delta k\sigma_{r}}{\sin\theta}$ taking values equal to $4.25$ and 2.91 for zero and unity, respectively. Note, that it is about $3$ times larger than for $\beta=1$. The rough estimation gives the width to be $\beta=10$ times larger. Thus, we see that, although overestimated, it is still a correct upper estimate.

IV.4 General correlation properties

IV.4.1 The $\Delta{\bf k}$ and $\Delta{\bf k}{\bf e}_{\bf k}$ widths

The gaussian ansatz provided us with Eq. (69), which expressed the single particle correlation function as a product of two factors, $G_{1}$ and $G_{2}$. The function $G_{2}$ is the Fourier transform of the initial condensate density squared. Basing on the analysis of the presented examples, we can improve the formula in Eq. (72) for the halwidth in $\Delta{\bf k}{\bf e}_{\bf k}$ to the following form:

[TABLE]

where the term $4/\sigma_{z}$ is an approximation of $4.2/\sigma_{z}$ obtained in the fast collision case.

It is interesting to note that the same result can be obtained by the inspection of Eq. (56), which we rewrite below for convenience,

[TABLE]

From this expression, we see the $\Delta{\bf k}$ dependence enters through $\Delta{\bf k}\cdot{\bf e}_{\bf k}v_{0}t$ and $\Delta{\bf k}\cdot{\bf r}$. Intuitively, the source size is of the order of the condensate size. We can therefore estimate the width in $\Delta{\bf k}$ given by the spatial integral as equal to the inverse size of the condensate width in respective direction. On the other hand, the width in $\Delta{\bf k}\,{\bf e}_{\bf k}$ can be estimated from expression in Eq. (76); it is equal to $\frac{1}{v_{0}\tau_{d}}$, where $\tau_{d}$ is the characteristic time given by (46) which results in the width of $\Delta{\bf k}\,{\bf e}_{\bf k}$ given by:

[TABLE]

The appearance of $\tau_{d}$ in this derivation is not surprising. As discussed in Section II.2, the time $\tau_{d}$ is the minimal characteristic time present in the evolution of functions $\psi_{\pm Q}$, which, through Eq. (59), are directly related to the source function $f$. Note that the above formula is the same estimate as given by Eq. (72) up to factor $2$ (which can appear in such rough estimates).

It is interesting to notice, that the width in $\Delta{\bf k}{\bf e}_{\bf k}$ given by Eq. (75) is close to the width in $K$ given by Eq. (48). In fact, this can be understood from the analysis of Eq. (53), which can be rewritten as

[TABLE]

where we used the symmetry $M({\bf k}^{\prime},{\bf k}_{2})=M({\bf k}_{2},{\bf k}^{\prime})$. We introduce now a convenient representation of variables, i.e., we rewrite the above in (${\bf k}$, $\Delta{\bf k}$, ${\bf K}$, $\Delta{\bf K}$): $M({\bf K},\Delta{\bf K})=M({\bf k}_{1},{\bf k}_{2})$, $G^{(1)}\left({\bf k},\Delta{\bf k}\right)=G^{(1)}\left({\bf k}_{1},{\bf k}_{2}\right)$. In addition to this, we change variable ${\bf k}^{\prime}$ into $\delta{\bf k}$ according to: ${\bf k}^{\prime}=-{\bf k}-\delta{\bf k}$. We then obtain the following expression:

[TABLE]

The two anomalous densities are taken for the vectors which magnitudes are given by:

[TABLE]

The difference of the two vectors is then $K_{1}-K_{2}\simeq\Delta{\bf k}\,{\bf e}_{\bf k}/2$. If $|K_{1}-K_{2}|$ is larger than $2\Delta_{K}$, the above integral is vanishingly small. Therefore, the width in $\Delta{\bf k}\,{\bf e}_{\bf k}$, denoted as $\Delta_{k}$, is approximately $4\Delta_{K}$. This relation is confirmed by equations (48) and (75) when we take into account that $\Delta_{K}$ and $\Delta_{k}$ given by these equations are rough estimates.

IV.4.2 Correlation volume

In analogy to the definition from the previous Section, we define here the correlation volume as the volume in $\Delta{\bf k}$ for which two particles are still significantly correlated.

The correlation volume is in fact related to the width in $\Delta k$ of $|G^{(1)}|$ function, calculated for arbitrary direction $(\sin\tilde{\theta}\cos\tilde{\phi},\sin\tilde{\theta}\sin\tilde{\phi},\cos\tilde{\theta})$. Combining all the results at hand regarding the local correlations, we can estimate this width as

[TABLE]

where ${\bf e}_{\Delta{\bf k}}\cdot{\bf e}_{\bf k}=\sin\tilde{\theta}\cos\tilde{\phi}\sin\theta+\cos\tilde{\theta}\cos\theta$ and $\Delta_{k}$ is given by Eq. (75). The factor $2.5$ comes directly from the halfwidth of the $G_{2}$ function given by (71).

Let us briefly analyze this formula. Due to the fact that ${1}/{\sigma_{r}}\gg\Delta_{k}$, the largest possible width, known to be ${2.5}/{\sigma_{r}}$, is possible only for small region around $\cos\tilde{\theta}=0$ and $\cos\tilde{\phi}=0$. In the remaining area we have a competition between second and third term of the above formula. However, as $\Delta_{k}\sigma_{z}\geqslant 4$, which is implied by Eq. (75), the second term of the above formula is less than the third one if ${|\cos\tilde{\theta}|}/{{\bf e}_{\Delta{\bf k}}\cdot{\bf e}_{\bf k}}\geqslant{2.5}/{4}$. For $|\cos\tilde{\theta}|\geqslant{2.5}/{4}$, the condition is satisfied, which gives us the region where $\delta_{\Delta k}={2.5}/{\sigma_{z}|\cos\tilde{\theta}|}$. On the other hand, for $\cos\tilde{\theta}\simeq 0$, excluding the region for which ${\bf e}_{\Delta{\bf k}}\cdot{\bf e}_{\bf k}\simeq 0$, we obtain $\delta_{\Delta k}={\Delta_{k}}/{{\bf e}_{\Delta{\bf k}}\cdot{\bf e}_{\bf k}}$.

Note that $\Delta_{k}$ is similar to $\Delta_{K}$ and furthermore the formulas for $\delta_{\delta k^{\prime}}$ and $\delta_{\Delta k}$ given by Eqs. (51) and (77) are similar as well. Thus the back to back and local correlation volumes take similar values.

Let us finally remark, that in the case ${\bf e}_{\Delta{\bf k}}\cdot{\bf e}_{\bf k}=0$ the local correlation width is given by the spatial Fourier transform of the source, which is given by the inverse of the condensate widths in respective directions. This is in fact the famous Handbury Brown and Twiss effect HBT0 . In the original work this effect was associated with the measurement of the star diameter HBT . Here, the star is replaced by a condensate – a source of particles.

V Summary

In this paper we analyzed the elastic scattering of atoms from elongated Bose-Einstein condensates colliding in the direction of the long axis which we choose to be $z$ axis. Our theory is valid in the collisionless regime, where multiple scattering processes are negligible. Additionally we focused our considerations on the spontaneous regime, in which the bosonic enhancement effect is negligible and the use of perturbation theory is justified. We showed that the two particle correlation functions decomposes into “back to back” and “local” parts which describe different aspects of the system. The “back to back” part $G^{(2)}_{bb}$ describes the correlation of two particles with almost opposite velocities. It reflects the fact that particles, due to binary collisions from counter propagating condensates, are scattered in pairs with almost opposite velocities. The “local” part $G^{(2)}_{loc}$ characterizes correlation of particles with velocities being almost the same, and is related to the bosonic bunching effect. Within perturbation theory, we derived approximate expressions for the one and two particle correlations functions of scattered atoms connecting these functions with the wavefunctions of the colliding condensates.

The formulas that we obtained are convenient for numerical computation when time dependent Gross-Pitaevskii equation is solved numerically. Furthermore, we introduced time dependent variational approach and obtained approximate form of the condensate wavefunctions. Having those, we calculated and analyzed back to back and local parts of the two particle correlation function in different regimes of parameters of the system. We found that the correlation function depends on two dimensionless parameters of the system: $\alpha=\frac{Qa_{hor}^{2}}{\sigma_{z}}$ and $\beta=\frac{\sigma_{r}^{2}}{a_{hor}^{2}}$ where $a_{hor}=\sqrt{\frac{\hbar}{m\omega_{r}^{2}}}$ is the radial harmonic oscillator length, Q is the mean wavevector of the colliding clouds and $\sigma_{r,z}$ denotes the radial and longitudinal sizes of the initial condensate.

In the chosen regime of parameters, the back to back part, $G^{(2)}_{bb}$, is a function of $\Delta{\bf K}={({\bf k}_{1}+{\bf k}_{2})}/{2}$ and $K=|{\bf K}|\equiv{|{\bf k}_{1}-{\bf k}_{2}|}/{2}$, where ${\bf k}_{1}$ and ${\bf k}_{2}$ are the wavevectors of scattered atoms. Analyzing few examples we have found that the $\Delta{\bf K}$ and $K$ dependence of the $G^{(2)}_{bb}(K,\Delta{\bf K})$ almost decouples. For $\beta=1$ and $\alpha\gg\beta$ cases the width in $\Delta K_{x,y}$ and $\Delta K_{z}$ are approximately equal to ${1}/{\sigma_{r}}$ and ${1}/{\sigma_{z}}$. This widths broaden maximally by factor of two with growing value of $\beta$ and decreasing value of $\alpha$. We found that $G^{(2)}_{bb}$ has the maximum at $K=Q+\delta K_{max}$ where $\delta K_{max}\approx\frac{m\mu}{\hbar^{2}Q}$ where $\mu$ is the chemical potential of the initial condensate. The width in $K$ of $G^{(2)}_{bb}$ can be estimated as $\frac{m}{2\hbar Q\tau_{d}}$ where $\tau_{d}$ is the characteristic time on which the wavefunctions of the colliding clouds changes substantially.

The local correlation, $G^{(2)}_{loc}({\bf k}_{1},{\bf k}_{2})$, of the two particle correlation function is directly related to a single particle correlation function, $G^{(2)}_{loc}=|G^{(1)}|^{2}$. In this study the convenient variables are ${\bf k}={({\bf k}_{1}+{\bf k}_{2})}/{2}$ and $\Delta{\bf k}={\bf k}_{1}-{\bf k}_{2}$. Exploiting the condensate wavefunctions given by the variational ansatz we have calculated $G^{(1)}({\bf k},\Delta{\bf k})$. Analyzing the problem on the single particle level, we found that the density is the narrowest in the $\beta=1$ and $\alpha\gg\beta$ cases. In other cases, the bulk region of the density stays practically the same while the tails of the density distribution start to grow together with the increase of the $\beta$ and decrease of $\alpha$. We have analyzed the dependence of $G^{(1)}({\bf k},\Delta{\bf k})$ in variable $\Delta{\bf k}$ in the bulk region of the density. We found that the ${\bf k}$ and $\Delta{\bf k}$ dependence decouples. We found that the dependence in $\Delta{\bf k}$ is given by two contributions. First comes from spatial Fourier transform of the source of particles, and is in fact the famous Hanbury Brown and Twiss effect. As a result the width in $\Delta{\bf k}$ is given by the inverse size of the initial condensate density in the respective direction. The second arises as a dependence in $\Delta{\bf k}\cdot{\bf e}_{\bf k}$ variable. We estimated the width in $\Delta{\bf k}\cdot{\bf e}_{\bf k}$ of $|G^{(1)}|$ to be roughly equal to $\frac{m}{\hbar Q\tau_{d}}$ and showed its direct relation to the width in $K$ of the $G^{(2)}_{bb}$.

Having all these results we found the back to back and local correlation volume properties and showed that they are similar in both cases.

Finally, we presented semiclassical models of both types of correlations. In the case of the back to back correlations we employed Wigner functions as a mean to describe the phase space densities of the colliding clouds. We showed that $G^{(2)}_{bb}$, taken from the semiclassical model, averaged over ${\bf K}$ yields the same results as the quantum model. Furthermore, we showed that large differences between quantum and semiclassical model appear in ${\bf K}$ dependence. In the case of the local correlations, we showed that, under specific requirements, the semicalssical formula for the source function is the same as the quantum one.

Acknowledgements.

We acknowledge the discussions with Denis Boiron, Jan Chwedeńczuk, Piotr Deuar, Karen Kheruntsyan Marek Trippenbach and Chris Westbrook. Paweł Ziń acknowledges the support of the Mobility Plus programme. Tomasz Wasak acknowledges the support of the Polish Ministry of Science and Higher Education programme “Iuventus Plus” (2015-2017) Project No. IP2014 050073, and the National Science Center Grant No. 2014/14/M/ST2/00015.

Appendix A Bogoliubov method: perturbative approach

Let us introduce the propagator defined by the equation:

[TABLE]

with the boundary condition:

[TABLE]

We further introduce the operator $\hat{\delta}^{\prime}({\bf r},t)$ by:

[TABLE]

The time $t=0$ is chosen in the propagator since at that time the evolution of the system starts. Substituting Eq. (80) into Eq. (2) we obtain:

[TABLE]

We multiply both sides of the above equation by $K({\bf r}^{\prime\prime},0;{\bf r},t)$ and integrate over ${\bf r}$ obtaining

[TABLE]

With help of the property

[TABLE]

Eq. (81) simplifies to

[TABLE]

We expand the field operator $\hat{\delta}^{\prime}({\bf r},t)$ in the following series:

[TABLE]

By substituting the above expansion into Eq. (82), we obtain an infinite hierarchy of equations. The first two of them reads:

[TABLE]

Eq. (84) can be solved straightforwardly: $\hat{\delta}^{\prime(0)}({\bf r},t)=\hat{\delta}^{\prime(0)}({\bf r},0)$. Substituting this into Eq. (85) and integrating over $t$ we obtain

[TABLE]

Subsituting from Eq. (86) into Eq. (80) we obtain the formal expression for the field operator $\hat{\delta}$, valid to the first order of expansion:

[TABLE]

where

[TABLE]

Now, let us find quantum averages of the field operator on the vaccum state from Eq. (5). From Eqs. (80), (79), (83) and (85) we obtain that

[TABLE]

Using the above, together with the definition of the vacuum state, see Eq. (5), we obtain

[TABLE]

which implies

[TABLE]

Combining these formulas, with the bosonic commutation relation, $[\hat{\delta}({\bf r},0),\hat{\delta}^{\dagger}({\bf r}^{\prime},0)]=\delta({\bf r}-{\bf r}^{\prime})$, we arrive at

[TABLE]

Now, we substitute expression from Eq. (87) into Eq. (6), which define anomalous density, and obtain

[TABLE]

The right hand side of this equation consists of four terms. As we see from Eq. (87) $\hat{\delta}^{(0)}({\bf r},T)$ is a function of $\hat{\delta}^{\prime(0)}({\bf r}^{\prime},0)$ while $\hat{\delta}^{(0)}({\bf r},T)$ a function of $\hat{\delta}^{\prime(0)\dagger}({\bf r}^{\prime},0)$. This fact, together with Eqs. (88) and (89), implies that

[TABLE]

As a consequence, we obtain the following form for the anomalous dnsity:

[TABLE]

Substituting now the forms of $\hat{\delta}^{(0)}$ and $\hat{\delta}^{(1)}$, given by Eq. (87), we obtain

[TABLE]

Using now the property of the operator $\delta^{\prime(0)}$ at the initial time, see Eq. (90), we can perform one integral over ${\bf r}^{\prime\prime\prime}$, which yields

[TABLE]

In the final step we invoke the properties of the propagator:

[TABLE]

which results in the following for the anomalous density:

[TABLE]

It can be interpreted in the following way. The probability amplitude of observing a pair of particles at ${\bf r}_{1}$ and ${\bf r}_{2}$ is the superposition of amplitudes of emitting a pair of particles, $B({\bf r},t)$, at position ${\bf r}$ and time $t$, which further propagate to final positions ${\bf r}_{1}$ and ${\bf r}_{2}$ and time $T$.

Let us calculate the one body correlation function given by Eq .(7). Substituting the form of $\hat{\delta}$ given by Eq. (87) into Eq. (7) we obtain:

[TABLE]

Now, using the same reasoning which led us to the anomalous density, we obtain that the only nonzero term in the above equation is $\langle 0|\hat{\delta}^{\dagger(1)}({\bf r}_{1},T)\hat{\delta}^{(1)}({\bf r}_{2},T)|0\rangle$. We substitute the form of $\hat{\delta}^{(0)}$ and $\hat{\delta}^{(1)}$ given by Eq. (87) to obtain:

[TABLE]

With the help of the Eqs. (90) and (91) we obtain the final formula for the single particle correlation function:

[TABLE]

Appendix B Approximate solution of the GP equation

In this appendix we find approximate initial state of the Bose-Einstein condensate and its subsequent time evolution.

B.1 Initial state

To find the initial state of cloud we rely on the variational method. To this end, we assume a gaussian profile:

[TABLE]

with the norm equal to $N$, and search for parameters $\sigma_{r}$ and $\sigma_{z}$ that minimize the Hamiltonian of the Gross-Pitaevskii Eq. (11):

[TABLE]

where $V({\bf r})$ is given by Eq. (12). Since the trapping frequency in radial direction is much larger than in the axial,$\omega_{r}\gg\omega_{z}$, the ground state wavefunction is very elongated. This allows us to approximate the Hamiltonian by neglecting the kinetic energy along $z$-direction:

[TABLE]

Now we insert the ansatz from Eq. (93) into the Hamiltonian in Eq. (95, which now becomes a function of two parameters $\sigma_{r}$ and $\sigma_{z}$. Taking the first derivatives, equating to zero, and solving, yields

[TABLE]

where $a$ is the scattering length, $a_{hor}=\sqrt{\frac{\hbar}{m\omega_{r}}}$, and $a_{hoz}=\sqrt{\frac{\hbar}{m\omega_{z}}}$. Note, that according to Eq. (96)

[TABLE]

Finally, Let us calculate a quantity that often appears in this paper, namely the meanfield energy $gn$, $n$ being the maximal density of the cloud. According to Eq. (93) and Eq. (96) this quantity is

[TABLE]

We also need chemical potential $\mu$, which by definition is the derivative of the ground state energy $E_{0}$ with respect to the number of atoms, $\mu=\frac{\partial E_{0}}{\partial N}$. Inserting the variational ansatz from Eq. (93) into Eq. (95), and using the formulas obtained above, we arrive at:

[TABLE]

B.2 Time evolution

Due to the symmetry of the trapping potential $V({\bf r})$, see Eq. (12), we have $\psi(x,y,z)=\psi(x,y,-z)$. This fact supplemented by Eq. (16) implies that $\psi_{+Q}({\bf r},0)=\psi_{-Q}({\bf r},0)$. Using now Eq. (17), we find the following property of the initial state:

[TABLE]

The term $\mp i\frac{\hbar^{2}}{m}Q\partial_{z}$ in Eq. (17) is responsible for the movement of the wave-packet with a constant velocity $\pm v_{0}=\pm\hbar Q/m$. It is thus convenient to define a new function by

[TABLE]

Exploiting Eqs. (101) and (102), we find that Eq. (17) can be rewritten as

[TABLE]

This equation describes the expansion of the wave-packet $\tilde{\psi}$ in the presence of the mean-field potential $g\left(|\tilde{\psi}({\bf r},t)|^{2}+2|\tilde{\psi}(x,y,-z-2v_{0}t,t)|^{2}\right)$. The term $g|\tilde{\psi}({\bf r},t)|^{2}$ has the same symmetry as $\tilde{\psi}$ and, thus, it increases the rate of the expansion. However, the second term in this potential is caused by the cloud $-Q$ which is moving with velocity $-2v_{0}$ with respect to the $+Q$ wavepacket, which results in the asymmetry in the expansion of the $\tilde{\psi}$ function. The goal of the present analysis is to obtain an approximate form of $\tilde{\psi}$. We therefore approximate the term $|\tilde{\psi}(x,y,-z-2v_{0}t,t)|^{2}$ by $|\tilde{\psi}(x,y,z,t)|^{2}$ restoring the symmetry. This approach leads to an “upper” limit on the expansion rate. Conversely, a “lower” limit can be obtained by completely neglecting the term $|\tilde{\psi}(x,y,-z-2v_{0}t,t)|^{2}$. In both cases, Eq. (103) takes the simple form

[TABLE]

where $\tilde{g}=3g$ or $\tilde{g}=g$ in the “upper” or “lower” limiting case, respectivelly.

We now use the variational method to solve the above equation, assuming that the norm $\tilde{N}=\int\!\!\mbox{d}{\bf r}|\psi({\bf r},t)|^{2}$; we dropped tilde from the wavefunction for the sake of notational clarity. The Eq. (103) can be formally derived from minimizing the action

[TABLE]

with the Lagrangian

[TABLE]

To obtain approximate evolution we assume a time-dependent variational ansatz:

[TABLE]

where $a_{r,z}(t)={1}/{2\sigma_{r,z}^{2}(t)}-ib_{r,z}(t)$, and $\sigma_{r,z}(t)$, $b_{r,z}(t)$ and $\phi(t)$ are real time dependent variational parameters. The norm of the profile is equal to $\tilde{N}$. We insert the above profile into the Lagrangian given by Eq. (105), and integrate over space variables to obtain

[TABLE]

The Euler-Lagrange equations of motion,

[TABLE]

lead to

[TABLE]

These can be transformed into more useful form:

[TABLE]

where we substituted $\tilde{g}=\frac{4\pi\hbar^{2}\tilde{a}}{m}$. The phase $\phi(t)$ is obtained from integrating GP equation:

[TABLE]

The primary object of our study is a strongly elongated system for which initially $\sigma_{z}(0)\gg\sigma_{r}(0)$. We can expect that in the course of time $\sigma_{r}(t)$ changes substantially, whereas the change of $\sigma_{z}(t)/\sigma_{z}(0)$ is small. We therefore approximate $\sigma_{z}(t)$ by $\sigma_{z}(0)$ in Eqs. (107) and (108), which supports analytical solutions of the following form:

[TABLE]

where

[TABLE]

and

[TABLE]

The phase is derived from Eq. (109) using Eqs. (107) and (108), and takes the form:

[TABLE]

After neglecting the kinetic energy along $z$-direction in the above equation, and taking $\sigma_{z}(t)\simeq\sigma_{z}(0)$, we obtain

[TABLE]

or its alternative form

[TABLE]

To derive the latter, we used Eqs. (114) and (115).

As can be seen from Eq. (110), the characteristic time on which $\sigma_{r}(t)$ changes its width is equal to $1/\tilde{\omega}$. To be consistent with the assumption stated in the derivation of the above formulas, the change of $\sigma_{z}(t)$ during that time has to be much smaller than $\sigma_{z}(0)$. Invoking Eq. (112), this condition takes the form:

[TABLE]

Let us now concentrate on the upper limit and take $\tilde{N}=\frac{N}{2}$, $\tilde{a}=3a$, with the initial condition $\sigma_{r,z}(0)=\sigma_{r,z}$ given by Eqs. (96) and (97). Then, the condition from Eq. (117) reads

[TABLE]

Now, from Eqs. (96), (114) and (115) we obtain

[TABLE]

Due to the following inequality,

[TABLE]

we obtain that

[TABLE]

where we have used Eq. (115). The high anisotropy, for which $\sigma_{r}\ll\sigma_{z}$, combined with Eq. (98) results in $a_{hor}\ll\sigma_{z}$. As a consequence, the first factor in condition from Eq. (118) is small, $\frac{\tilde{a}_{hor}^{4}}{2\sigma_{z}^{4}}\ll 1$. From Eqs. (96) and (97), we obtain

[TABLE]

Thus, we have shown that the condition in Eq. (117) is satisfied, making the above derivation self-consistent.

Now, we simplify the expression for $b_{z}(t)$ given in Eq. (113). To this end, notice that the maximal phase that appears in the variational ansatz is roughly equal to $b_{z}(t)\sigma_{z}^{2}$. Then, the first term in expression (113) takes the form $\frac{\tilde{a}_{hor}^{2}}{4\sigma_{z}^{2}}\tilde{\omega}t$. As $\omega t$ is maximally of the order of few, and due to elongation of the system, $\tilde{a}_{hor}\ll\sigma_{z}$, this term can be neglected. As a result, we obtain

[TABLE]

with the alternative form

[TABLE]

where Eqs. (114) and (115) were used.

Let us notice that in the case of “lower” limit, for which $\tilde{g}=g$, we obtain different values of $\tilde{a}_{hor}$ and $\tilde{\omega}$, that satisfy following chain of inequalities:

[TABLE]

We may conclude that there is no substantial qualitative difference between these two cases. Therefore, in the rest of the paper we take $\tilde{g}=2g$ as a interpolation between the lower and upper limits. In such a case, we have $\tilde{a}_{hor}=a_{hor}$ and $\tilde{\omega}=\omega_{r}$. Finally, using Eqs. (101), (102), (106) together with Eqs. (110), (111), (116) and (121) we obtain approximate forms for $\psi_{\pm Q}$, which are given by:

[TABLE]

where $\beta=\frac{\sigma_{r}^{2}}{a_{hor}^{2}}$, which, according to Eq. (98), cannot be less than unity, $\beta\geqslant 1$.

The momentum density corresponding to the wavefunction is given by

[TABLE]

where $k_{r}^{2}=k_{x}^{2}+k_{y}^{2}$, and

[TABLE]

We note following useful relations: initially $\sigma_{kr,kz}(0)={1}/{\sigma_{r,z}}$, the final value of the width in radial is $\sigma_{kr}(\infty)={\beta}/{\sigma_{r}}={\sigma_{r}}/{a_{hor}^{2}}$, and the axial width is $\sigma_{kz}(\infty)=({1}/{\sigma_{z}})\sqrt{1+\left(\beta-{1}/{\beta}\right)^{2}{\pi^{2}}/{4}}$.

In Section II we stated two assumptions. Now, with help of the derived formulas, they can be written as

[TABLE]

Using Eq. (99), the first of the above conditions can be rewritten as

[TABLE]

Thus, if

[TABLE]

both of the conditions present in Eq. (125) are satisfied.

Appendix C Inclusion of the mean field propagator

C.1 Construction of the mean-field propagator

Here we describe how to construct approximate formula for the propagator $K$ of the Hamiltonian given in Eq. (3). The mean field potential present $\hat{H}_{0}$ equals $2g|\psi({\bf r},t)|^{2}$. Using the decomposition from Eq. (16) we arrive at

[TABLE]

The mean field potential decomposes into two parts: a slowly varying envelope part $V_{en}({\bf r},t)=2g(|\psi_{Q}({\bf r},t)|^{2}+|\psi_{-Q}({\bf r},t)|^{2})$, and an oscillating part $V_{osc}({\bf r},t)$ with the fringes oscillating as $\cos(2Qz)$.

We begin with neglecting the oscillating part of the potential $V_{osc}$ in the mean field potential, the scattered particles are then influenced only by $V_{en}$. We comment on this approximation below. To remind, we restricted our analysis to the collision of highly elongated cigar shaped condensates along the longitudinal $z$-direction. Also, we are interested only in the atoms that scattered away from the condensates with velocities distant from the $z$-axis. The time needed for these atoms to leave the cloud is approximately equal to ${\sigma_{r}}/{v_{0}}$. During this interval each of the condensates moves by the distance equal to $v_{0}({\sigma_{r}}/{v_{0}})=\sigma_{r}$, which is much smaller than the longitudinal size $\sigma_{z}$ of the cloud. On the other hand, during this time the radial width of the condensates, given by Eq. (20), increases by a factor

[TABLE]

which follows directly from the condition in Eq. (14). Thus, the condensate densities, and similarly $V_{en}$, do not change appreciably during the time that takes the scattered atoms to escape the clouds. In such a case, we are legitmate to approximate the time dependent scattering one-body problem

[TABLE]

by

[TABLE]

where $\hbar\omega_{\bf k}=\frac{\hbar^{2}k^{2}}{2m}$ and $\varphi({\bf r},t)$ is a plane wave far away for the potential $\lim_{r\rightarrow\infty}\varphi({\bf r},t)=\frac{1}{(2\pi)^{3/2}}e^{i{\bf k}{\bf r}}$. The solution of the above equation is not unique. However, in scattering theory we identify two solutions denoted as $\varphi_{\bf k}^{(+)}$ and $\varphi_{\bf k}^{(-)}$ with the boundary condition

[TABLE]

The function $\varphi_{\bf k}^{(+)}$ describes a physical situation where an incident particle comes from $r=-\infty$ and scatters on the potential, resulting in the scattered waves that are directed outside the potential. The $\varphi_{\bf k}^{(-)}$ describes the time-reversed situation where the scattered waves are directed toward the potential. For a given $t$ both of these sets of functions form a complete and orthogonal set zderzenia . Thus, from both of these sets we can construct a propagator

[TABLE]

This is our approximation to the true propagator of the Hamiltonian $\hat{H}_{0}$ with neglected $V_{osc}$. As we will see below, ${\bf r}_{2}$ and $t_{2}$ is the position and time the scattered particle is produced in the condensate whereas ${\bf r}_{1}$ and $t_{1}$ is the time and position of the measurement. As the detection is far away from the condensate then $\varphi({\bf r}_{1},t)\simeq\frac{1}{(2\pi)^{3/2}}e^{i{\bf k}{\bf r}}$, and it does not depend on $t$. The time $t$ has to be taken as some mean time between the time of the birth of the scattered particle and the time the particle leaves the cloud. We have shown above that on the time the particle leaves the cloud $V_{en}$ practically does not change. Thus, we can take $t=t_{2}$, which results in

[TABLE]

We are now ready to derive an approximate analytical formula for $\varphi_{\bf k}^{(+)}$. In Section IV we showed that under presented approximations the width of the halo of scattered atoms is much smaller than its radius being close to $Q$. Thus, we search for $\varphi_{\bf k}^{(+)}$ only for $k$ close to $Q$. The characteristic length on which the potential $V_{en}$ changes, which is of the order of $\sigma_{r}$, is much larger than the mean wavelength of the scattered atom $\frac{2\pi}{Q}$ (this fact follows from Eq. (14) and the fact that $\sigma_{r}>a_{hor}$ as stated in Eq. (98)). In such a situation the use of semiclassical approximation is justified. In subsection C.3 of this Appendix we show that under condition given in Eq. (145) the wavefunction $\varphi_{\bf k}^{(+)}$ takes the approximate form

[TABLE]

where ${\bf e}_{\bf k}=\frac{{\bf k}}{k}$. The above formula is derived in zderzenia and is known as the “eikonal approximation”. The expression for $\varphi^{(+)}_{{\bf k}}({\bf r})$ is a correct approximation in the part of space before, understood as ${\bf r}$ that satisfy ${\bf k}\cdot{\bf r}<0$, as well as on the potential. After the potential the scattering part appears which is clearly not present in the above formula. This means that the form of $\varphi^{(+)}_{{\bf k}}({\bf r})$ is no longer given by Eq. (131). The problem is that in our calculations we need to know the form of $\varphi^{(+)}_{{\bf k}}({\bf r})$ on and also after the potential. The way to overcome this issue is to realize that in the construction of the propagator $K$ we can use the states $\varphi^{(-)}$ instead of $\varphi^{(+)}$. These states satisfy

[TABLE]

so $\varphi_{\bf k}^{(+)}$ before the potential is $\varphi_{\bf k}^{(-)}$ after the potential. As a result we take

[TABLE]

together with

[TABLE]

which is defined on and after the potential $V_{en}$. In deriving the above, we used the fact that ${\bf e}_{-{\bf k}}=-{\bf e}_{\bf k}$.

Finally, let us comment on the omission of $V_{osc}$ part of the total mean-field potential. This potential has fringes represented by the term $\cos 2Qz$. Thus, in terms of the time dependent perturbation theory the potential couples incoming plane wave $e^{i{\bf k}{\bf r}}$ to the plane wave $\exp\left(i{\bf k}{\bf r}\pm i2Qz\right)$. The matrix coupling elements of $V_{osc}\propto gn$, and so the probability amplitude of such coupled wavefunction should be proportional to

[TABLE]

which is much less that unity. Therefore, we neglect $V_{osc}$ in our considerations.

C.2 Anomalous density

Here we analyze the formula for the anomalous density. Inserting Eq. (133) into Eq. (9), omitting the superscript $(-)$, we obtain

[TABLE]

Taking this as a starting point, we derive a formula for $M({\bf k}_{1},{\bf k}_{2})$ defined in Eq. (21). For large ${\bf r}_{1}$ and ${\bf r}_{2}$ the wavefunctions $\varphi_{{\bf k}_{1}^{\prime}}({\bf r}_{1},t)$ and $\varphi_{{\bf k}_{2}^{\prime}}({\bf r}_{2},t)$ are plane waves (close to the potential the scattered part of the wave-functions may still be present):

[TABLE]

where we have used ${\bf r}_{1,2}=\frac{\hbar{\bf k}_{1,2}}{m}T$. Introducing ${\bf k}_{1}^{\prime}={\bf k}_{1}+\delta{\bf k}_{1}$ and ${\bf k}_{2}^{\prime}={\bf k}_{2}+\delta{\bf k}_{2}$ we obtain

[TABLE]

Using the above together with Eqs. (135) and (21) we obtain

[TABLE]

Here the term $\exp\left(-i\frac{\hbar}{2m}(\delta k_{1}^{2}+\delta k_{2}^{2})T\right)$ for $T\rightarrow\infty$ serves as an effective Dirac delta function:

[TABLE]

As a result, we obtain

[TABLE]

Now, from Eq. (134) we obtain

[TABLE]

The phase in square brackets is an integral over two straight lines meeting at point ${\bf r}$. In Section III we showed that when the free propagator is used the wavevctors ${\bf k}_{1}$ and ${\bf k}_{2}$ for which $M({\bf k}_{1},{\bf k}_{2})$ has non-vanishing value are practically anti-parallel, with the length approximately equal to $Q$. Let us show that the same applies here as well. To this end, let us analyze the temporal and spatial dependence of the integrand in Eq. (136). We have two terms with such a dependence: $B({\bf r},t)$ and the phase $\frac{m}{\hbar^{2}Q}\int^{0}_{-\infty}\!\!\mbox{d}s\,\left(V_{en}({\bf r}-s{\bf e}_{{\bf k}_{1}},t)+V_{en}({\bf r}-s{\bf e}_{{\bf k}_{2}},t)\right)$. As found in the Section III, the dominant temporal phase present in $B({\bf r},t)$ is given by $\exp\left(-{\hbar^{2}Q^{2}t}/{m}\right)$. On the other hand, the integral $\int_{-\infty}^{0}\mbox{d}s\,\left(V_{en}({\bf r}-s{\bf e}_{{\bf k}_{1}},t)+V_{en}({\bf r}-s{\bf e}_{{\bf k}_{2}},t)\right)$ can be estimated as $gn\sigma_{r}$. The temporal change of this integral can be estimated as $gn\sigma_{r}\frac{t}{t_{c}}$, where $t_{c}$ is the characteristic time of the change of the potential $V_{en}$. As we have shown in subsection C.1 of this Appendix, the change of the potential $V_{en}$ takes place during time which is much larger than $\frac{\sigma_{r}}{v_{0}}$. Thus, we have $t_{c}=\frac{1}{\epsilon}\frac{\sigma_{r}}{v_{0}}$ with $\epsilon\ll 1$. As a result, the temporal dependence of $\frac{m}{\hbar^{2}Q}\int^{0}_{-\infty}\mbox{d}s\,\left(V_{en}({\bf r}-s{\bf e}_{{\bf k}_{1}},t)+V_{en}({\bf r}-s{\bf e}_{{\bf k}_{2}},t)\right)$ can be approximated by

[TABLE]

As $\frac{\hbar^{2}Q^{2}}{2m}\gg gn$, which is the assumption stated in Section II, the above is much smaller than the temporal phase of $B$ equal to $-i{\hbar^{2}Q^{2}t}/{m}$. When we integrate this two terms with $\exp\left(i{\hbar(k_{1}^{2}+k_{2}^{2})t}/{2m}\right)$, present in Eq. (136), we obtain $k_{1}^{2}+k_{2}^{2}\simeq 2Q^{2}$.

Let us now analyze the spatial dependence of the integral $\int_{-\infty}^{0}\mbox{d}s\,\left(V_{en}({\bf r}-s{\bf e}_{{\bf k}_{1}},t)+V_{en}({\bf r}-s{\bf e}_{{\bf k}_{2}},t)\right)$. We approximate it as a drop from the maximal value $gn$ to zero on a distance equal to $\sigma_{r}$. For simplicity, we take $x$ in the direction of the drop which results in $\int_{-\infty}^{0}\mbox{d}s\,\left(V_{en}({\bf r}-s{\bf e}_{{\bf k}_{1}},t)+V_{en}({\bf r}-s{\bf e}_{{\bf k}_{2}},t)\right)\approx gn\sigma_{r}\frac{x}{\sigma_{r}}$. Thus, the phase factor of the analyzed term can be approximated by $\exp\left(-i{mgnx}/{\hbar^{2}Q}\right)$. As $gn\ll\frac{\hbar^{2}Q^{2}}{2m}$, the above term results in a much smaller shift of the wavevector than $Q$. On the other hand, according to the assumption stated in Section II the momentum width of the function $\psi_{\pm Q}$ is much smaller than $Q$. This two facts together with the definition of $B({\bf r},t)$, given by Eq. (24), imply that $|{\bf k}_{1}+{\bf k}_{2}|\ll Q$. Thus, we have shown that indeed the mean-field term coming from the propagator does not change the fact that ${\bf k}_{1}$ and ${\bf k}_{2}$ have lengths close to $Q$ and are approximatelly anti-parallel.

We can now simplify the expression for $M({\bf k}_{1},{\bf k}_{2})$, given by Eqs. (136) and (137). As ${\bf k}_{1}$ and ${\bf k}_{2}$ are almost anti-parallel, we approximate

[TABLE]

where ${\bf K}=\frac{{\bf k}_{1}-{\bf k}_{2}}{2}$. Consequently, from Eqs. (136) and (137) we finally obtain

[TABLE]

where

[TABLE]

From the results of subsection C.3 of this Appendix, the maximal value of $|\Phi|$ can be estimated as $\frac{m2gn}{\hbar^{2}Q}4\sigma_{r}$. Therefore, if

[TABLE]

the phase can be neglected and then we have $B({\bf K},{\bf r},t)\simeq B({\bf r},t)$. With help of Eq. (99) this condition can be rewritten in the following form

[TABLE]

If the condition

[TABLE]

is true, the one in Eq. (140) is satisfied.

C.3 Approximate solution of the scattering problem

Let us now analyze the scattering problem given by Eq. (126). It can be rewritten as

[TABLE]

To simplify the notation, in this subsection omit superscript $(+)$, subscript ${\bf k}$ and time $t$ in $\varphi_{\bf k}^{(+)}({\bf r},t)$ and $V_{en}({\bf r},t)$. To find the approximate solution of this Schrödinger equation we substitute $\varphi({\bf r})=e^{i{\bf k}{\bf r}+i\phi({\bf r})}/{(2\pi)^{3}}$. Then, the Schrödinger equation takes the following form

[TABLE]

We solve this equation by a perturbation series $\phi=\phi^{(0)}+\phi^{(1)}+...$. We have the set of equations

[TABLE]

The solution of the above is

[TABLE]

where $W_{0}=\frac{2m}{\hbar^{2}}V_{en}$ and $W_{1}=(\nabla\phi^{(0)})^{2}-i\triangle\phi^{(0)}$. The maximal value of $|\phi^{(j)}|$ can be estimated as

[TABLE]

where $|W_{j}|_{max}$ is the maximal value of the function $|W_{j}|$. In the above, we have taken the path through the center of the condensate with the smallest possible value of $\sin\theta=\frac{1}{2}$ and the effective width equal to $\frac{2\sigma_{r}}{|\sin\theta|}=4\sigma_{r}$. In this way we arrive at

[TABLE]

where we took $|V_{en}|_{max}=2gn$. As the width of $\phi^{(0)}_{1}$ is of the order of $\sigma_{r}$ in the axial direction, we estimate $\nabla$ and $\triangle$ operators acting on $\phi^{(0)}$ as $\frac{1}{\sigma_{r}}$ and $\frac{1}{\sigma_{r}^{2}}$. Then, using Eqs. (142) and (143), we obtain

[TABLE]

According to the assumption stated in Section II, $\frac{\hbar^{2}Q^{2}}{2m}\gg gn$, so the second term on the righthand side of the above equation can be neglected. Thus, from Eq. (141), approximating $\frac{1}{2k}\simeq\frac{1}{2Q}$, we obtain

[TABLE]

as long as

[TABLE]

is satisfied.

Appendix D Derivation of expressions used in Section III

In this Appendix we calculate the expressions presented in Section III, in the order they appear in the text.

First, we investigate the semiclassical model. The Wigner function in the case of gaussian ansatz, see Eq. (19), takes the form

[TABLE]

With this Wigner function we calculate the function $G^{(2)}_{bb}$, given by Eq. (30),

[TABLE]

where we introduced dimensionless time $\tau=\omega_{r}t$.

Now, we derive the anomalous density, given by Eq. (33), in the fast collision case. As $\tau_{c}$ is much smaller than all the other characteristic times we can approximate $c_{z}(\tau)\simeq 1$, and $1+\tau^{2}\simeq 1$, and $1-i\beta\tau\simeq 1$. Additionally, $\phi(\tau)\ll 1$ and we can neglect it. Thus, performing the temporal integral we arrive at

[TABLE]

where

[TABLE]

The term $\exp\left(-\frac{\Delta K_{r}^{2}\sigma_{r}^{2}}{4}\right)$ gives the width $\Delta K_{r}$ to be approximately equal to $1/\sigma_{r}$. Then, we have $\frac{\Delta K_{r}^{2}\sigma_{z}}{4Q}\approx\frac{\sigma_{z}}{4Q\sigma_{r}^{2}}$. According to condition in Eq. (35), and the facts that $\sigma_{r}\geqslant a_{hor}$, see Eq. 98, the above term is much smaller than unity and can be neglected. As a result we obtain

[TABLE]

Now, we derive the semiclassical expression in the fast collision case. In the formula from Eq. (147) all the temporal dependence, apart from $-2\alpha^{2}\tau^{2}$, can be neglected. We therefore obtain

[TABLE]

where we approximated $|K^{\prime}|\simeq Q$ and performed the temporal integral. Taking $\mbox{d}{\bf K}^{\prime}={K^{\prime}}^{2}\mbox{d}K^{\prime}\mbox{d}\cos\theta\mbox{d}\phi$ and integrating over $K^{\prime}$ and $\phi$ we arrive at

[TABLE]

where $z=\cos\theta$. We now introduce $y=1-z$ and $\delta K=K-Q$. In these new variables, the integrand takes the form $\exp\left(-2K^{2}\sigma_{r}^{2}(2y-y^{2})+2(\delta K(1-y)-Qy)^{2}\sigma_{z}^{2}\right)$. As shown in Section III the value of $K$ is close to $Q$, so the term $-4K^{2}\sigma_{r}^{2}y$ gives the characteristic width in $y$ being approximately equal to $\frac{1}{4Q^{2}\sigma_{r}^{2}}$. Due to the condition from Eq. (14), and the fact that $\sigma_{r}\geqslant a_{hor}$, given by Eq. (98), this width is much smaller than unity, and thus the term $2K^{2}\sigma_{r}^{2}y^{2}$ can be neglected. The last term reads $-2(\delta K(1-y)\sigma_{z}-Qy\sigma_{z})^{2}$. We have $Qy\sigma_{z}\approx\frac{\sigma_{z}}{4Q\sigma_{r}^{2}}\leq\frac{\sigma_{z}\sigma_{r}^{2}}{Qa_{hor}^{4}}\ll 1$, where we used $\sigma_{r}\geq a_{hor}$ and the fast collision condition given in Eq. (35). Additionally, as $y\ll 1$ we can approximate $\delta K(1-y)\sigma_{z}\simeq\delta K\sigma_{z}$. Therefore, we obtain

[TABLE]

As a result of all the approximations, we finally obtain

[TABLE]

The semiclassical expression from Eq. (147) in the strong confinement case takes the form

[TABLE]

where we put $K^{\prime}\simeq Q$. We now perform the temporal integral and, as before, partly perform the integral over ${\bf K}^{\prime}$ obtaining:

[TABLE]

We again introduce $y=1-z$ and $\delta K=K-Q$, which can be used to rewrite the term $\exp\left(-2K^{2}a_{hor}^{2}(1-z^{2})-2(Kz-Q)^{2}\sigma_{z}^{2}\right)=\exp\left(-2K^{2}a_{hor}^{2}(2y-y^{2})-2(\delta K(1-y)-Qy)^{2}\sigma_{z}^{2}\right)$. As before, the term $-4K^{2}a_{hor}^{2}y$ implies that the width of $y$ is maximally ${1}/{4Q^{2}a_{hor}^{2}}\ll 1$. Consequently, we neglect the term $2K^{2}a_{hor}^{2}y^{2}$ and approximate $-2(\delta K(1-y)-Qy)^{2}\sigma_{z}^{2}\simeq-2(\delta K-Qy)^{2}\sigma_{z}^{2}$. As $K\simeq Q$, we also approximate the term $-4K^{2}a_{hor}^{2}y\simeq-4Q^{2}a_{hor}^{2}y$. Furthermore, we approximate the term present in the denominator $K^{2}a_{hor}^{2}(1-z^{2})\simeq 2Q^{2}a_{hor}^{2}y$ and extend the limits of integration from $2$ to $\infty$. As a result we obtain

[TABLE]

where we changed the variables to $z=2Q^{2}a_{hor}^{2}y$.

Below, we derive expression for the back to back part of the two particle correlation function $|M({\bf K},\Delta{\bf K})|^{2}$ averaged over ${\bf K}$. From Eq. (26) we have

[TABLE]

where $\tilde{B}({\bf K},{\bf r},t)=B({\bf K},{\bf r},t)\exp\left(i\frac{\hbar Q^{2}}{m}t\right)$, $\Delta{\bf K}={\bf k}_{1}+{\bf k}_{2}$, and we used ${(k_{1}^{2}+k_{2}^{2})}/{2}=K^{2}+{\Delta K^{2}}/{4}$. As shown in Section III the value of $K$ is close to $Q$, so we substitute $K=Q+\delta K$. Then, the Eq. (149) takes the following form

[TABLE]

where we decomposed the integral over ${\bf K}$; note that $B({\bf K},{\bf r},t)$ does not depend on $K$ but only on the direction ${\bf e}_{\bf K}$. Due to the fact that $\delta K\ll Q$, we approximate $\int_{0}^{\infty}K^{2}\mbox{d}K\simeq Q^{2}\int_{-\delta Q}^{\delta Q}\mbox{d}\delta K$, where $\delta Q$ is much smaller than $Q$ but still bigger than the width in $\delta K$. Then, the last integral in Eq. (150) takes the form

[TABLE]

where $x=\delta K/Q$ and $x_{0}=\delta Q/Q$. If we take $x_{0}=0.1$ the above integral gives us a peaked function in $t_{1}-t_{2}$ with a width of $(t_{1}-t_{2})={5m}/{\hbar Q^{2}}$. For larger $(t_{1}-t_{2})$ the function oscillates rapidly. The condition in Eq. (14) together with the fact that $Q\sigma_{z}\gg 1$ imply that $\frac{5m}{\hbar Q^{2}}\ll\tau_{c},\tau_{r}$, which shows that during the time $(t_{1}-t_{2})={5m}/{\hbar Q^{2}}$ the change of the wavefunctions $\tilde{\psi}_{\pm Q}({\bf r},t)$ are negligible. Thus, the integration over $\delta K$ results in the Dirac delta function $\int_{-\Delta Q}^{\Delta Q}\mbox{d}\delta K\,\exp\left(i\frac{\hbar}{m}(2Q\delta K+\delta K^{2})(t_{1}-t_{2})\right)\simeq\frac{\pi m}{\hbar Q}\delta(t_{1}-t_{2})$. As a result, we obtain

[TABLE]

The classical expression in Eq. (30) integrated over ${\bf K}$ takes the form

[TABLE]

Approximating now $K^{\prime}\simeq Q$, substituting Wigner functions given by Eq. (31) and performing the integral over ${\bf K}^{\prime},$ which results in the Dirac delta function, one obtains exactly the same formula as in Eq. (151) with $B({\bf K},{\bf r},t)=B({\bf r},t)$ given by Eq. (24).

Now, we calculate the expression in Eq. (151) subject to the condition presented in Eq. (29), when $B({\bf K},{\bf r},t)\simeq B({\bf r},t)$, in the case of gaussian ansatz, see Eqs. (19) and (20). Using Eq. (24) and performing gaussian integrals we obtain

[TABLE]

This expression can be compared to the momentum density of $\psi_{\pm Q}$, given by Eq. (19), which takes the form

[TABLE]

where $\psi({\bf k},t)=\int\mbox{d}{\bf r}\,e^{-i{\bf k}{\bf r}}\psi({\bf r},t)$, and we used Eq. (20). By comparing the two expressions above we find that

[TABLE]

The expression in Eq. (153) integrated over $\Delta K_{z}$ takes the form:

[TABLE]

where $C_{b}=\frac{4Q(Na)^{2}}{(2\pi)^{3/2}\sigma_{z}^{2}\sigma_{r}^{2}}$ and $\tau=\beta\omega_{r}t$. In deriving this last expression we used Eq. (20).

Appendix E Metastable helium experiment parameters

In the metastable helium experiment the excited states of atoms are used for which the scattering length is $a=7.51$ nm. The parameters of the Palaiseau experiment are paryz0

[TABLE]

while the parameters of Vienna experiment are Wieden

[TABLE]

With these parameters we solve numerically Eq. (97) to obtain $\sigma_{z}$, and substitute it to Eq. (96) to obtain $\sigma_{r}$. With these results, we calculate $\alpha$ and $\beta$, which in the case of Palaiseau experiment take the values $\alpha\simeq 0.22$, $\beta\simeq 3.3$, while for Vienna experiment $\alpha\simeq 0.2$, $\beta\simeq 11$.

Appendix F Formulas needed in Section IV

F.1 Derivation of Eq. (56)

From the definition of the Wigner function we have

[TABLE]

Inserting this expression into Eq. (52), we obtain

[TABLE]

where we introduced ${\bf k}=\frac{{\bf k}_{1}+{\bf k}_{2}}{2}$ and $\Delta{\bf k}={\bf k}_{1}-{\bf k}_{2}$. Introducing new variable ${\bf r}^{\prime}=\frac{\hbar({\bf k}-{\bf k}^{\prime})T}{m}$ and using Eq. (54), we obtain

[TABLE]

Introducing now ${\bf r}={\bf r}^{\prime}+\frac{\hbar}{m}{\bf k}t$, we finally obtain

[TABLE]

F.2 Derivation of the source function $f({\bf r},{\bf k},t)$

In this Appendix we derive Eq. (57), valid as long as the condition in Eq. (29) is satisfied. The condition yields $B({\bf k},{\bf r},t)\simeq B({\bf r},t)$, and the anomalous density is then given by Eq. (III.1). Using Eqs. (III.1) and (24) together with Eqs. (53) and (56), we find that

[TABLE]

where

[TABLE]

and $K_{f}$ denotes the free propagator.

The formulas given by Eqs. (154) and (155) are the basis for understanding the physics of the processes. We first show, that only for $|\Delta t|<\Delta t_{0}=2C{m\sigma_{r}}/{\hbar Q}$ the integrand in these equations is nonzero. Here, $C$ is of the order of unity; it means that $C$ can be equal to few but not to few tens. This estimation stems from considering the integration over $\Delta x$. As the system has cylindrical symmetry, we take ${\bf k}=(k_{x},0,k_{z})$ without the lost of generality. The part of the integrand with $\Delta x$ dependence consist of free propagator $K_{f}\left(\Delta x,\Delta t\right)\propto\exp\left(i\frac{m\Delta t^{2}}{2\hbar\Delta t}\right)$ multiplied by the source function $B$. We introduce the phase of $B$ as $B=|B|e^{i\phi}$ and write $K_{f}\left(\Delta x,\Delta t\right)B\propto|B|\exp\left(i\frac{m\Delta x^{2}}{2\hbar\Delta t}+i\phi\right)$. The function $B$ has the width in the radial direction approximately equal to $\sigma_{r}$, since it is proportional $\psi_{+Q}\psi_{-Q}$. This means that the function $|B|$ is nonzero only for

[TABLE]

where $C$ is of the order of a unity. Let us now note that in Eq. (154) we deal with the product of two $B_{p}$ functions. The analysis of the integral of $B_{p}^{*}$ results in analogous condition

[TABLE]

where $\Delta{\bf r}^{\prime}$ denotes the integration variable. The above conditions, Eqs. (156) and (157), can be rewritten as

[TABLE]

As shown in Appendix C.2, the phase $\phi$ is maximally equal to $\epsilon Qx$ where $\epsilon\ll 1$. In such case the function $\exp\left(i\frac{m\Delta x^{2}}{2\hbar\Delta t}+i\phi\right)$ is an oscillating function with decreasing period of oscillation when moving from the central point $\Delta x_{0}=\epsilon\frac{\hbar Q\Delta t}{m}$. The characteristic width of $\Delta x$ equals $2\sqrt{\frac{\hbar|\Delta t|}{m}}$. As a result, the second of the above inequality takes the form

[TABLE]

According to the assumption stated in Section II, we restrict our analysis to the situation for which $k_{x}>\frac{Q}{2}$. Thus, $\left|\frac{Q}{k_{x}}\frac{\epsilon-\epsilon^{\prime}}{2}\right|<|\epsilon-\epsilon^{\prime}|$ and can be neglected. Let us notice that for $\Delta t$ equal to $C\frac{m\sigma_{r}}{\hbar k_{x}}$, we have $\frac{\hbar k_{x}\Delta t}{m}=C\sigma_{r}$, and

[TABLE]

because $Q\sigma_{r}\gg 1$. The above inequality implies that $\Delta t<C\frac{m\sigma_{r}}{\hbar k_{x}}$. Using $k_{x}>\frac{Q}{2}$, we obtain

[TABLE]

We have shown that the timescale of $\Delta t$ is ${m\sigma_{r}}/{\hbar Q}$, a time that is needed for the scattered particle to leave the cloud.

Basing on the results presented above, we now continue with approximating the expression in Eq. (155). First, we notice that due to condition in Eq. (29) and the fact that $\sigma_{z}\gg\sigma_{r}$, we have $\Delta t_{0}\ll\tau_{c},\tau_{r}$. This means that the changes of the wavefunctions $\psi_{\pm Q}$ can be neglected during the time $\Delta t_{0}$. Thus, according to the Eq. (24), we have

[TABLE]

where $\tilde{\bf r}={\bf r}+{\hbar{\bf k}\Delta t}/{m}$. We found above that the effective width in $\Delta x$ given by the propagator $K_{f}(\Delta{\bf r},\Delta t)$ is $\Delta_{0}\ll\sigma_{r}$. As the propagator is a function of $\Delta{\bf r}$ thus the width in $\Delta y$ and $\Delta z$ is also $\Delta_{0}$. On such a distance the change of the function $|\psi_{+Q}\psi_{-Q}|$ can be neglected. As a result, we obtain

[TABLE]

Finally, we analyze the phase of $\psi_{+Q}\psi_{-Q}(\tilde{\bf r}+\Delta{\bf r},t)$. Using gaussian ansatz given by Eq. (19), together with Eq. (20), we find that the $x$ dependent phase of the integrand in Eq. (155) equals

[TABLE]

It can be rewritten as

[TABLE]

where $\epsilon={2\hbar\Delta t\beta\omega_{r}t}/{m\sigma_{r}^{2}(t)}$, and its maximal value can be estimated to be $\frac{2\hbar\Delta t_{0}}{ma_{hor}^{2}}=4C\frac{\sigma_{r}}{Qa_{hor}^{2}}$. According to the condition (14), this is much smaller than unity, and we obtain $1+\epsilon\simeq 1$. Thus the phase is equal to

[TABLE]

The same reasoning can be repeated in the case of $y$ and $z$ coordinates. Thus, the Eq. (159), together with Eq. (160), gives

[TABLE]

Using the above together with Eq. (158) make Eq. (155) to take the form

[TABLE]

In order to perform the above approximation we assumed the quadratic form of the phase of the wavefunctions $\psi_{\pm Q}$ with position independent coefficients as given by the gaussian ansatz. However we note that it is enough that the phase would be well approximated by quadratic function on a distance $\Delta_{0}\ll\sigma_{r}$. As this is generally true thus the above reasoning apllies for true (not only variational) solution of the GP equation. Using the above equation the expression in Eq. (154) for the source function simplifies to

[TABLE]

Furthermore, introducing new variable $\delta k=k-Q$ we obtain $k^{2}-Q^{2}=2Q\delta k\left(1+\frac{\delta k}{2Q}\right)$. As $|\delta k|\ll Q$ we can approximate $k^{2}-Q^{2}\simeq 2Q\delta k$ and, additionally, $\frac{\hbar{\bf k}\Delta t}{m}\simeq\frac{\hbar Q\Delta t}{m}{\bf e}_{\bf k}$, where ${\bf e}_{\bf k}=\frac{{\bf k}}{k}$. The timescale of $t$ is equal or larger to one of the characteristic times $\tau_{c}$ or $\tau_{r}$. As we noticed before, $\Delta t_{0}\ll\tau_{c},\tau_{r}$, thus, most of the collision takes place at times $t$ much larger than $\Delta t_{0}$. Therefore, we extend the limit of integration in the above integral to infinity,$\int_{-t}^{t}\simeq\int_{-\infty}^{\infty}$. With these approximations, we arrive at

[TABLE]

where we changed the variables $\delta r={\hbar Q\Delta t}/{m}$. Sill, this formula can be further simplified basing on appropriate approximations. Specifically, from Eq. (20) we obtain that the maximal characteristic change of the function $\psi_{\pm Q}$ in the $z$ direction on the distance $\sigma_{r}$ is approximately equal to

[TABLE]

As $\beta=\sigma_{r}^{2}/a_{hor}^{2}$, this quantity is much smaller than unity, as long as the condition given by Eq. (58) is satisfied. Then, we can neglect the dependence on ${\bf e}_{z}$ in term $\delta r{\bf e}_{\bf k}$, arriving at

[TABLE]

where ${\bf e}_{{\bf k},r}=\frac{k_{x}{\bf e}_{x}+k_{y}{\bf e}_{y}}{k}$.

F.3 Source function: semiclassical model

We show here that the expression from Eq. (57) can be obtained from a semiclassical model introduced in Section III. The formula for a source function in such a model takes the form:

[TABLE]

where the function

[TABLE]

with ${\bf K}=\frac{{\bf k}-{\bf k}^{\prime}}{2}$ and $\Delta{\bf K}={\bf k}+{\bf k}^{\prime}$.

We first analyze the function $h$ given by Eq. (165) in the case of a gaussian ansatz, see Eq. (19). Using the formula for the Wigner function from Eq. (31) we arrive at

[TABLE]

where:

[TABLE]

and

[TABLE]

We now analyze the above functions discussing the values of $K_{max}$ and $\Delta{\bf K}_{max}$ for which the function $h(K_{max},\Delta{\bf K}_{max},{\bf r},t)$ has maximum, together with the widths in $K$ and $\Delta K_{r,z}$ around the maximum. From the form of $h_{3}$ one can deduce that the width in $K$ is approximately equal to $\Delta_{K}=\frac{1}{\sigma_{z}}+\frac{1}{2Q\sigma_{r}^{2}}$ with the maximum located between $Q-\delta_{K}$, where $\delta_{K}<\frac{\beta}{\sigma_{z}}=\frac{\sigma_{r}^{2}}{a_{hor}^{2}\sigma_{z}}$, and $Q-\delta_{K}+\Delta_{K}$. Additionally, it can be seen from the $h_{2}$ that $|\Delta K_{max;x,y}|<\frac{2\beta}{\sigma_{r}}=\frac{2\sigma_{r}}{a_{hor}^{2}}$ and $|\Delta K_{max;z}|<\frac{2\beta}{\sigma_{z}}$, with the width in $\Delta K_{r,z}$ around the maximum equal to $\frac{1}{\sigma_{r}}$ and $\frac{1}{\sigma_{z}}$, respectively. The width in $K$ is much smaller than the width in $\Delta K_{x,y}$. We make use of this fact below in performing an approximation.

To this end, we first find ${\bf K}_{0}$ and $\Delta{\bf K}_{0}$ satisfying ${\bf k}={\bf K}_{0}+\frac{\Delta{\bf K}_{0}}{2}$, for which the function $h(K_{0},\Delta{\bf K}_{0},{\bf r},t)$ has maximum. Note, that without the restriction ${\bf k}={\bf K}_{0}+\frac{\Delta{\bf K}_{0}}{2}$, we would have $|{\bf K}_{0}|=K_{max}$ and $\Delta{\bf K}_{0}=\Delta{\bf K}_{max}$, but here it does not need to be the case. Therefore, we introduce $\delta{\bf k}^{\prime}={\bf k}^{\prime}+{\bf K}_{0}-\frac{\Delta{\bf K}_{0}}{2}$. We obtain

[TABLE]

where ${\bf e}_{0}=\frac{{\bf K}_{0}}{K_{0}}$. Next, we change the variables of integration in Eq. (164) from ${\bf k}^{\prime}$ to $\Delta{\bf K}$, and further on to $\delta{\bf k}^{\prime}$. Thus, the integral over $k_{x}^{\prime}$ takes the following form

[TABLE]

In the considerations presented above, we found that the width in $K$ is much smaller than the width in $\Delta K_{x}$. The integration over $\delta k^{\prime}_{x}$ can be effectively changed into integration over $K$ ($\int\mbox{d}\delta k^{\prime}_{x}\rightarrow\frac{2}{e_{0,x}}\int\mbox{d}K$) keeping $\Delta K_{x}$ constant, and equal to $\Delta K_{0,x}$. This is true, if $e_{0,x}={\bf e}_{0}\cdot{\bf e}_{x}$ is of the order of unity. We take ${\bf k}=(k_{x},0,k_{z})=k(\sin\theta,0,\cos\theta)$ without the lost of generality. As we shall see below, for such a choice $e_{0,x}>\frac{1}{2}$. Thus, we arrive at

[TABLE]

where we also approximated $1/K\simeq 1/Q$ and $e_{0,x}\simeq e_{k,x}$. Note, that the dependence on ${\bf k}$ is hidden in $\Delta K_{0,x}$. In what follows, we explicitly find this dependence.

To this end, we start from the equality

[TABLE]

Taking ${\bf e}_{0}=(\sin\theta_{0}\cos\phi_{0},\sin\theta_{0}\sin\phi_{0},\cos\theta_{0})$ each of the components of ${\bf k}$ takes the form:

[TABLE]

Using the fact that $|\Delta K_{0,z}|\ll Q$, $k\simeq K_{0}\simeq Q$, and introducing $\theta_{0}=\theta+\delta\theta$, the Eq. (173) takes the approximate form

[TABLE]

where $\delta k=k-Q$. Consistently with the approximations, that were undertaken in order to obtain this formula, we have $|\delta\theta|\ll 1$. We then have $|\Delta K_{0,z}|\simeq|\Delta K_{max;z}|<\frac{2\beta}{\sigma_{z}}$, which together with $|\sin\theta|>\frac{1}{2}$ results in $\left|\frac{\Delta K_{0,z}}{2Q\sin\theta}\right|\leq\frac{2\sigma_{r}^{2}}{Qa_{hor}^{2}\sigma_{z}}$. Estimating now $|K_{0}-Q|\leq\delta_{K}+\Delta_{K}$, we obtain $\left|\frac{K_{0}-Q}{Q}\cot\theta\right|\leq\frac{\sigma_{r}^{2}}{Qa_{hor}^{2}\sigma_{z}}+\frac{1}{2Q^{2}\sigma_{r}^{2}}+\frac{1}{Q\sigma_{z}}$, where we used $|\cot\theta|<1$. Thus, we obtain

[TABLE]

As $|\Delta K_{0,y}|\ll K_{0}\simeq Q$, Eq. (172) can be solved approximately, leading to

[TABLE]

Thus, we have $|\delta\theta|\ll 1$ and $|\phi|\ll 1$, which results in $e_{0,x}=e_{k,x}={\bf e}_{\bf k}\cdot{\bf e}_{x}$. According to the assumption stated in Section II, we have $e_{0,x}>\frac{1}{2}$, and so the same applies to ${\bf e}_{0,x}$. Using $|\delta\theta|\ll 1$ and $|\phi|\ll 1$, Eq. (171) can be approximately rewritten as

[TABLE]

With the help of Eqs. (174) and (176), we rewrite the above equation as

[TABLE]

Proceeding in the same way as in the case of $r_{1}$, additionally invoking Eq. (175), we obtain that

[TABLE]

If the righthand side of the above is smaller than the width od $\frac{\Delta K_{0,x}}{2}$, equal to $\frac{1}{2\sigma_{r}}$, which gives

[TABLE]

then we obtain

[TABLE]

Let us now note that, as $\sigma_{r}\geqslant a_{hor}$, $\sigma_{z}\gg\sigma_{r}$, and from the condition given in Eq. (29) we obtain

[TABLE]

Using this results, the condition in Eq. (177) can be brought to the following form

[TABLE]

Inserting Eq. (178) into Eq. (170), using Eqs. (165) and (31), and performing the integrals, we arrive at the semiclassical form of the source function,

[TABLE]

In deriving this results, we additionally used Eq. (32). If we change the variable $\delta re_{k,x}=-\frac{\Delta x}{2}$ in this equation, we arrive at the quantum formula given in Eq. (163).

Appendix G Validity of the perturbative approach

In this Appendix we inspect the validity of the perturbative approach. The assumption presented in Eq. (29) implies that we can effectivelly neglect the mean-field potential $2g|\psi({\bf r},t)|^{2}$ in the Hamiltonian $H_{0}$, given by Eq. (3). As a result, Eq. (2) takes the form

[TABLE]

Substituting here

[TABLE]

we obtain

[TABLE]

where

[TABLE]

With the definitions $M({\bf k}_{1},{\bf k}_{2},t)=\langle\hat{\delta}({\bf k}_{1},t)\hat{\delta}({\bf k}_{2},t)\rangle$ and $G^{(1)}({\bf k}_{1},{\bf k}_{2},t)=\langle\hat{\delta}^{\dagger}({\bf k}_{1},t)\hat{\delta}({\bf k}_{2},t)\rangle$, from Eq. (179) we obtain that

[TABLE]

It can be proved that the anomalous density $M({\bf k}_{1},{\bf k}_{2})=iM({\bf k}_{1},{\bf k}_{2},\infty)$. Therefore, after integrating the equation above we obtain

[TABLE]

Now we apply the first order perturbation, which amounts to neglecting the second line of this equation. As a result, we obtain

[TABLE]

which is exactly formula given in Eq. (III.1). The perturbation theory, expressed in Eqs. (III.1) and (53), used in the paper are valid as long as $\frac{1}{\hbar}\int_{0}^{\infty}\mbox{d}t\int\mbox{d}{\bf k}^{\prime}(A({\bf k}_{1},{\bf k}^{\prime},t)G^{(1)}({\bf k}^{\prime},{\bf k}_{2},t)+A({\bf k}_{2},{\bf k}^{\prime},t)G^{(1)}({\bf k}^{\prime},{\bf k}_{1},t))$ is much smaller than $\frac{1}{\hbar}\int_{0}^{\infty}\mbox{d}t\,A({\bf k}_{1},{\bf k}_{2},t)$.

We now estimate the value of

[TABLE]

where we used ${\bf k}^{\prime}={\bf k}_{2}+\delta{\bf k}^{\prime}$. For simplicity of the calculation, we take ${\bf k}_{1}=-{\bf k}_{2}={\bf k}=Q{\bf e}_{\bf k}$. Using that, together with Eqs. (24) and (180), we obtain

[TABLE]

In rough approximation $\psi_{+Q}\psi_{-Q}({\bf r},t)$ decomposes into ${\bf r}$ and $t$ dependent parts. Additionally as $|\delta{\bf k}^{\prime}|\ll Q$ we have in most directions of ${\bf k}$, $\left|\frac{\hbar}{2m}(\delta{\bf k}^{\prime})^{2}\right|\ll|v_{0}{\bf e}_{\bf k}\delta{\bf k}^{\prime}|$ which makes us to neglect the term $\exp\left(i\frac{\hbar}{2m}(\delta{\bf k}^{\prime})^{2}t\right)$. As a result the approximate form of $A$ reads

[TABLE]

Inserting the above into Eq. (181) we obtain

[TABLE]

As $\Delta{\bf K}=\delta{\bf k}^{\prime}$, the function $A_{0}(\delta{\bf k}^{\prime})$ has the same width as the anomalous density, equal to ${1}/{\sigma_{r,z}}$ in the respective directions. The fact that the width in $K=\left|{\bf k}-\frac{\delta{\bf k}^{\prime}}{2}\right|\simeq Q-\frac{1}{2}{\bf e}_{\bf k}\delta{\bf k}^{\prime}$ of the anomalous density is equal to $\Delta_{K}$, implies that the width in $t$ of the function $A_{1}(t)$ is equal to ${1}/{2v_{0}\Delta_{K}}$. We normalized $A_{0}(0)=1$ and $\int_{0}^{\infty}\mbox{d}t\,A_{1}(t)=1$, which results in that $M_{0}({\bf k})$ is equal to the anomalous density $M({\bf k},-{\bf k})$ in the first order perturbation theory.

We now focus our attention on $G^{(1)}({\bf k}_{1},{\bf k}_{2},t)$. Using derivation analogous to the ones presented in Appendix F, one can show that

[TABLE]

In a rough approximation

[TABLE]

where $\rho$ is the density and $G_{2}$ is given by Eq. (71). The width in ${\bf e}_{-{\bf k}}\delta{\bf k}^{\prime}$ of $G_{1}$ is equal to $\Delta_{k}$, which implies that the width in $t^{\prime}$ of the function $f_{0}(t^{\prime})$ is ${1}/{v_{0}\Delta_{k}}$. The fact that $G^{(1)}(-{\bf k},-{\bf k},\infty)=\rho(-{\bf k})$ implies that $\int_{0}^{\infty}\mbox{d}t^{\prime}\,f_{0}(t^{\prime})=1$. Inserting now Eqs. (183) and (184) into Eq. (182) we obtain

[TABLE]

Due to the fact that $\Delta_{K}\approx\Delta_{k}$, the widths of both functions $A_{1}$ and $f_{0}$ are similar. Therefore, the temporal integrals give the width in ${\bf e}_{\bf k}\delta{\bf k}^{\prime}$ approximately equal to $\Delta_{k}$. Furthermore, the widths in $\delta{\bf k}^{\prime}$ of the functions $A_{0}$ and $G_{2}$ is similar. As a result, we obtain

[TABLE]

where $V_{c}$ is the correlation volume. Note that back to back and local correlation volumes are similar thus in rough approximation we take them to be the same and denote by $V_{c}$.

We conclude that the considered term is approximately equal to

[TABLE]

where the first order perturbation term is equal to $M_{0}({\bf k})$. The term $\frac{1}{\hbar}\int_{0}^{\infty}\mbox{d}t\int\mbox{d}{\bf k}^{\prime}A({\bf k}_{2},{\bf k}^{\prime},t)G^{(1)}({\bf k}^{\prime},{\bf k}_{1},t)$ shall take similar value as written above. Thus the condition for the validity of the perturbation approach is that $\rho(-{\bf k})V_{c}\ll 1$, i.e., the mean number of particles scattered into the correlation volume has to be much smaller than unity.