Squeezing in Bose-Einstein condensates with large numbers of atoms

We examine the feasibility of creating and measuring large relative number squeezing in multicomponent trapped Bose-Einstein condensates. In the absence of multimode effects, this squeezing can be arbitrarily large for arbitrarily large condensates, but a range of processes limit the measurable squeezing in realistic trap configurations. We examine these processes, and suggest methods to mitigate them. We conclude that high levels of squeezing with large numbers of atoms is feasible, but can realistically only be achieved in particular trap geometries. We also introduce a method of maximising the measurable squeezing by using a $\pi$-pulse during the process to improve spatial mode-matching.

Ultracold atomic interferometers are also reaching state of the art sensitivities in precision measurements [30,31]. While all current atom interferometry is limited by technical noise sources, it appears that BEC-based coherent sources offer the best chance of overcoming some of these issues, for the same reasons that optical lasers are typically preferred for optical interferometry. In particular, it appears that Bosecondensed sources allow for higher fidelity mirrors and beamsplitters [32], are less sensitive to the effects of distortions in the optical wavefronts of the beamsplitters and mirrors and the Coriolis effect [30,33], and are robust to the adverse effects of strong outcoupling when feedback-cooled [34][35][36].
If these technical noise sources can be reduced and/or controlled, atom interferometers will be limited by atomic shot noise, which leads to a sensitivity that scales as the square root of the atomic flux. From that point, sensitivity will only be improved by increasing the atomic flux and brightness [37][38][39], and/or using nonclassical quantum states such as squeezed states to go below the shot noise limit [30].
Producing squeezed states for metrological purposes is therefore only relevant when it can be achieved in a context of large numbers of atoms. This paper examines the generation of atomic squeezing using the inherent Kerr-like nonlinearities of ultracold atoms, and demonstrates a trap geometry in which metrologically relevant levels of squeezing are present in trap geometries that are compatible with large atomic number.
The creation, detection and application of squeezing requires precise modematching, and is therefore very sensitive to uncontrolled coupling between spatial modes of the atomic field. Atomic squeezing experiments have typically attempted to operate in the single-mode regime, where the trap is tight enough such that all higher modes are frozen out [40]. However, tight traps limit the total number of atoms, as the 3-body recombination rate places an effective maximum atomic density in the gas, and well before that limit, the interatomic interactions break the single-mode operation. Both of these effects limit the number of atoms that can be squeezed effectively, so we will examine a broader range of trap geometries, and calculate the dynamics in the presence of multiple spatial modes. Many similar previous simulations have been performed in one or two spatial dimensions, by integrating out some of the remaining dimensions [5][6][7]12]. For simulations of semiclassical fields this can be a reasonable approximation, but it can ignore spontaneous scattering processes in full quantum field theory calculations, so we take care to include full 3D simulations.
In Section 2 we examine the ideal behaviour of a nonlinear bosonic two-state system when that system possesses only a single spatial mode, extracting full analytic solutions for the number squeezing in each of the two states, as well as the number difference squeezing between the two states. These solutions of the ideal case provide a check on the multimode cases we consider later in the paper, as well as providing a bound on the best possible squeezing for the system. In Section 3 we show how multimode behaviour becomes important in physical systems, and introduce our model for simulating higher dimensional trap geometries. Section 4 analyses one-dimensional models, highlighting some of the issues regarding optimising squeezing parameters and mode-matching in squeezing experiments. Section 5 shows that the application of a π-pulse during the process can significantly improve mode-matching issues. Section 6 investigates measurable squeezing in 2D and 3D, and uses a Bogoliubov analysis to explain the scaling with system size and dimensionality. Section 7 describes the conclusion that large volume, flat-bottomed, three-dimensional traps, which are obviously compatible with large atom number, can produce extremely high degrees of squeezing.

Model and single-mode solutions
Our goal is to describe the number squeezing possibilities of a BEC comprising two relevant internal states denoted by |a and |b . The Hamiltonian for a pair of coupled atomic internal modes is the spatial integral of the Hamiltonian density: whereψ j (r) is the annihilation field operator for a particle at position r and in internal state |j , m is the mass of the atoms, V j (r) is the trapping potential for atoms in internal state |j , κ describes the coupling between the two states, and the U ij describe the various inter-and intra-state nonlinearities and are given by where a ij are the inter-and intra-species s-wave scattering lengths. We begin by considering a simplified version of the problem, assuming that only one of these spatial modes is relevant for each of the condensate's components. This is done by writing the multimode atomic field operators in terms of a set of basis functionsψ a (r, t) = n u n (r)â n (t) andψ b (r, t) = n u n (r)b n (t), whereâ n andb n annihilate a particle with normalised spatial mode u n (r) and internal state |a or |b respectively. With this expansion, the nonlinear term in the multimode Hamiltonian for the self-interaction of the |a internal state becomeŝ When the system remains predominantly in spatial mode u 0 (r), we can approximate this sum as the single termĤ = χ aaâ †â †ââ /2 where We proceed similarly for the other terms. In the parameter regimes we are interested in, the detuning between the two spatial modes will be irrelevant, so setting our zero of energy appropriately, we are left with only a standard Kerr-type nonlinearity arising from atom-atom interactions, as well as a linear coupling between the two modes. The Hamiltonian governing this reduced system is given bŷ where the χ ij describe the various inter-and intra-mode nonlinearities, and Ω couples the two modes and allows for population transfer between them. When considering squeezing, the Hamiltonian (5) has typically been handled within the framework of spin squeezing, which has a long history [9,19,25,28,41,42]. This approach involves exploiting the equivalence between the algebra of two harmonic oscillators and that of angular momentum, and writinĝ Using these variables Eq. (5) becomeŝ whereN =â †â +b †b is the total number operator, ∆ω(N ) = (N − 1) (χ aa − χ bb ) /2 and χ + = χ aa + χ bb − 2χ ab . For initial states which are total number Fock states, the last two terms are physically meaningless phase shifts and the first term can be removed by an appropriate choice of detuning. In this case Eq. (11) reduces to the Josephson Hamiltonian with zero detuning [42] The system described byĤ J is then solved using angular momentum algebra. The reduction of Eq. (5) to the Josephson Hamiltonian has led to the common misconception that no squeezing can be generated in this system if χ + = 0 as the effects of the ∆ω term are usually dismissed as an irrelevant detuning. This is not true in general for initial states other than a total number Fock state as the ∆ω(N )Ĵ z term leads to a number-dependent detuning. We demonstrate in this paper that for the case of an initial coherent state this term leads to squeezing. A spin-echo pulse applied half-way through the experiment which exchanges the population of modes |a and |b (which is the protocol used in most spin squeezing experiments) also cancels the effects of this term, and in this case squeezing is only generated if χ + = 0. In this work we directly solve Eq. (5) demonstrating that squeezing is generated (in the absence of a spin-echo pulse) even if χ + = 0.
The system is prepared with all the population initially in |a , and vacuum in |b . We choose the initial state |a to be a coherent state with average number â †â = N . As the system is insensitive to the initial phase of the coherent state, this is equivalent to a mixture of coherent states of uncertain phase, or a mixture of Poisson-distributed number states. Such a state is consistent with BEC coherence experiments [43].
At time t = 0, the coupling Ω is applied until time t = t 1 , resulting in a portion of the population being transferred into mode |b . We assume this transfer process is fast relative to the nonlinear energy time scale, ensuring that, as the initial states were coherent states, after the population transfer both |a and |b are still described by coherent states with mean number â †â = n a and b †b = n b respectively. The coupling Ω is then switched off until time t 2 , while the atoms interact solely through the nonlinear terms. After this hold time τ hold = t 2 − t 1 , the coupling Ω is switched back on until time t 3 , with a phase shift φ compared to its first application. During this last stage, the two modes exchange population and the quadrature variances are converted into number variance, in the same way that homodyne measurements are used in quantum optics to convert quadrature squeezing into number squeezing, which can be directly measured. In this interpretation, φ is the relative phase of the strong local oscillator, which allows specific phase angles of the quadrature squeezing to be examined. This experiment can also be interpreted as a Ramsey interferometer with a final beam splitter phase of φ. The entire sequence of pulses and the resulting populations are shown schematically in Figure 1. We solve this system analytically, and derive expressions for absolute and relative number squeezing. We present the solutions here; the full derivation can be found in Appendix A. If we define then the expectation values for number, number variance, and number difference variance after applying the final coupling field are given by = n a s 2 + n b c 2 + 2n a n b c 2 s 2 In the case where all the nonlinearities are equal, such that χ aa = χ ab = χ bb , the number variances all reduce to those of a coherent state, i.e. Var[N a ] = N a , and no number squeezing is possible. It should also be noted that, although there are three independent nonlinearities χ ij in the Hamiltonian (5), the expressions for the number variances only depend on differences. This means that number squeezing is parameterised by only two independent quantities, for example χ aa −χ ab and χ aa −χ bb .
In Figure 2 we plot a solution for the normalised number variance in mode |a as a function of φ, choosing a typical set of nonlinearities. This example illustrates the basic features of the system: during the final coupling period, the normalised number variance undergoes cyclic variation with a period τ = π/Ω, and has a series of minima. It is necessary to choose both the correct phase φ as well as the correct length of time for the coupling to be applied in order to obtain optimum squeezing. As mentioned previously, squeezing is possible in this system even when the nonlinearity in the reduced Josephson Hamiltonian χ + = χ aa + χ bb − 2χ ab = 0. Figure 3 demonstrates that significant squeezing is possible in this system if a spin-echo pulse is not used.
This two-mode analytic model allows arbitrarily good number squeezing, provided there are no limits to the number of particles in the system or to the hold time. The two-mode approximation can be valid in tight traps with very low particle number, and this result correctly predicts the squeezing demonstrated by Gross et al. [10]. Unfortunately, for larger numbers of atoms in realistic traps there are more than two modes present, and these multimode effects limit the level of squeezing that can be achieved.

Multimode number squeezing analysis
The zero-dimensional, two-mode model described in the previous section shows that arbitrarily good squeezing can be achieved using the intrinsic nonlinearities in a BEC, but it ignores the existence of multiple spatial modes. There are three problems caused by the existence of multiple modes, and these will limit the achievable squeezing.
The first problem is present even if the spatial modes are uncoupled. The two-mode model shows that the parameters required to obtain best squeezing are dependent on the strength of the nonlinearity. In the multimode case, the effective nonlinearity of a mode is a function of the mode shape, as is shown by Eq. (4). The optimum hold time τ hold depends on the product of the effective mode nonlinearity and the number of particles in that mode, as does the recombination phase φ. Different modes in a multimode environment will typically have varying effective nonlinearityparticle number products and therefore squeeze at different rates, with each mode having a different optimum τ hold and optimum recombination phase φ. This means that the best squeezing will be achieved by choosing values of τ hold and φ that, averaged across all the spatial modes present, result in the overall lowest total number variance. Clearly, however, this averaged number variance will be still be higher than the specific mode that exhibits the best squeezing.
The second problem is due to the number-dependent dynamics of the spatial modes, which is observable even in a semiclassical simulation. Due to the nonlinear term in the Hamiltonian, the local phase evolution of the atomic field is density dependent. Consequently the mode shape of the atomic field changes depending on the nonlinearity-density product, which means that unless U aa |ψ a (r)| 2 + U ab |ψ b (r)| 2 = U bb |ψ b (r)| 2 + U ab |ψ a (r)| 2 for all r, the mode shapes will not overlap perfectly during the final coupling pulse. This degraded mode-matching will result in lower efficiency when converting the quadrature squeezing to number squeezing.
We might speculate that the effects of these first two problems would be lessened by trap geometries that lead to near-constant density profiles for the BEC. This should enable higher squeezing. We cannot make strong conclusions without also considering the effects of coupling between spatial modes, however, which leads to the third potential problem: coupling between modes is inevitable in the presence of a nonlinearity, meaning that the modes cannot be analysed independently. Coupling between modes will disturb mode-matching, as well as mixing fields with different phase evolution, which can rapidly destroy squeezing.
An analytic solution to the multimode, higher-dimensional problem is intractable, so we must turn to numerical solutions. We use stochastic methods based on phase space representations to simulate the dynamics of the multimode quantum fields. These methods alleviate the problem of directly simulating elements of the Hilbert space, the size of which increases exponentially with the number of spatial modes [44,45]. Stochastic methods achieve this by finding a sufficiently well-behaved quasiprobability representation for the density matrix, which can then be simulated as the average behaviour of a number of low-dimensional samples. These sample objects are the same size as the equivalent classical field. In our case, the only tractable phase space representation is the functional Wigner representation [45], which is well behaved when the quantum field is approximately Gaussian, as it is for coherent and squeezed states. Stochastic methods based on the functional Wigner representation have been used to model the behaviour of BECs and atom lasers [6,46,47], as the Wigner representation allows the mapping of the system to a Fokker-Planck equation with positive semi-definite diffusion matrix (unlike the P -representation [44]), and does not suffer exponential path-weighting problems that would constrain it to very short simulation intervals (as would the positive-P representation [44]).
The equation of motion for the density matrix is defined by the Hamiltonian in Eq. (1). This defines the evolution of the functional Wigner distribution of the system, which has a one-to-one correspondence with the density matrix. The equation of motion for the functional Wigner distribution contains derivatives of third order which we assume to be negligible. This uncontrolled approximation is called the Truncated Wigner Approximation (TWA), and while it has been used widely on ultracold gases [6,7,45,47,48], care must be taken that the simulation remains valid. For our simulations we checked the results of our simulations against the analytic twomode solution, as well as checking for typical indications of TWA breakdown such as the appearance of negative densities in lightly populated modes. We were also able to use the Bogoliubov analysis described in Section 6 as an independent estimate of the validity of the TWA. We found that the TWA remained valid except in some 3D simulations with high nonlinearities, where there were a large number of modes, and the number of particles per mode could drop to ten or less in the densest regions where the squeezing was generated. Over very long time scales, this low mode occupation began to result in TWA breakdown. Fortunately, the best squeezing was typically found in regions where the TWA was valid.
Under the TWA, the equation of motion for the functional Wigner distribution is a Fokker-Planck equation, and can therefore be sampled by a set of stochastic equations, as described in [44]. These stochastic fields φ(x) are related to field operator expectation values by where : ⋆ : sym denotes symmetric ordering of the operators, and E denotes a stochastic average over the variables φ(x). While we see below that the evolution of these fields is deterministic, the initial state still requires a random element, so multiple realisations are required. When approximating these equations on a discrete grid, the magnitude of initial noise, the relationship between stochastic averages and the expectation values, and the equations of motion all become grid-dependent. This is not a signature of a fundamental problem, as the Hamiltonian in Eq. (1) assumes a contact potential between the atoms, which is not correct below a length scale that can probe the details of the true potentials. Furthermore, it is not a problem in practice, as the calculation gives grid-independent predictions for physical observables well before that regime is reached.
These stochastic equations take the form where , dV is the volume element of the grid on which the simulation is carried out, and the terms proportional to 1/dV correspond to vacuum corrections. It is a peculiarity of the Wigner representation that, up to the vacuum correction terms, Eqs. (28) look identical to those of the coupled semi-classical nonlinear Schrödinger equation. The quantum statistical nature of the Wigner equations enters due to the fact that Eqs. (28) are run many times and stochastically averaged, with each run using different set of random initial conditions describing (in this case) the noise on a coherent state.
We integrate Eqs. (28) numerically in one, two and three spatial dimensions using the numerical package XMDS2 [49]. For the one-and two-dimensional simulations, the dimensional reduction is achieved by estimating the mode shape in three dimensions, and using dimensionally reduced values for U ij that match the zero-dimensional reduction given by Eq. (4). To make relevant comparisons between equivalent situations with different numbers of spatial dimensions, this mode shape is typically chosen to match the chemical potentials of the initial states.

Effects of trap geometries on squeezing in 1D
A one-dimensional simulation is sufficient to examine the hypothesis that the squeezing will be higher for a BEC where the spatial mode has a more constant density profile. We compare squeezing for fields starting in the ground states of a harmonic trap with negligible nonlinearity (a Gaussian), a harmonic trap with strong nonlinearity (Thomas-Fermi), and a constant potential (constant density). Using Eq. (4) we find  Figure 4: Effects of mode shape on squeezing. The plot shows the normalised number variance in mode |a as a function of the time the final coupling pulse is applied. We compare the analytic two-mode solution (thick red), the constant density mode (solid blue), the Thomas-Fermi mode (dashed) and the Gaussian mode (dash-dotted). The nonlinearities of the various 1D simulations have been adjusted for their mode shape so that they are equivalent to that of the 0D, two-mode mode nonlinearity. Parameters used were n a = n b = 1 × 10 5 , χ aa = 0.03s −1 , χ ab = χ bb = 0, τ hold = 3 × 10 −4 s, φ = 4.0. Maximum squeezing is achieved when the atomic density profile is the most uniform.
the following equivalences where ω x , ω y , ω z are the angular frequencies of the harmonic trap, and N is the number of particles in the mode of interest. Choosing parameters such that the squeezing occurs fast enough that the mode shape is reasonably stable over the simulation gives the best squeezing, and allows us to make a fair comparison between the squeezing achievable by different mode shapes without being concerned that different modes will change shapes at different rates, ruining mode-matching. The results of this simulation are shown in Figure 4, which clearly confirms the hypothesis that maximum squeezing is obtained when the mode shape is closest to constant density.
One way to avoid any requirement to engineer and maintain specific atomic modes is to post-select a spatial mode. Rather than considering the variance in the total number of atoms in each internal state, it may be possible to consider only the statistics in a spatially filtered area of the atomic cloud. For example, suppose the atoms were in a Gaussian atomic distribution, but we consider a mode that only includes the atoms within p standard deviations of the centre. For low p, this mode shape is considerably closer to the constant density case than a standard Gaussian. The equivalence between the single mode case and this mode is given by Such spatial filtering can easily produce higher squeezing, but there is the obvious disadvantage that it significantly reduces the effective number of atoms. This is a serious issue for most applications of squeezed sources, such as interferometry, where the signal to noise scales with the square root of the flux. We will focus on methods for producing squeezing with large atomic number.

Using a π-pulse to improve mode-matching
Modal mismatch is generated by the different atomic potentials seen by the two internal states. A straightforward method of minimising this effect is to apply a π-pulse to swap the populations of the internal states after half of the evolution. As can be seen by Eqs. (A.8) and (A.9), such a pulse also leaves the quantum statistics of each state unchanged. The idea of this scheme is to make the effective potential seen by the atoms in state |a during the second half of the hold time very similar to that seen by the atoms that were in state |a during the first half of the hold time, and will only work if the modes |a and |b have an equal occupation at the start of the hold time, or if χ + = χ aa + χ bb − 2χ ab = 0. In the latter case, this reduces to a spin-echo pulse, and as the squeezing effects of the ∆ω term in Eq. (11) are cancelled in this case, there will be no squeezing. However, in the limit of short hold times, for equal occupations of the two modes, and when χ + = 0, this technique reduces the difference in the mode shapes after the hold time. In practice this will never work perfectly, as ψ a (r) and ψ b (r) will typically have evolved by the time the π-pulse is applied, and so the spatial dynamics in the second half of the hold time will not be quite the same as those during the first half. Nonetheless, in many regimes, such a π-pulse can produce some improvement. More importantly, as this scheme accepts that the two modes will change shape, and attempts to make them change in a similar fashion, it allows the possibility of much longer hold times. That is, if we are no longer constrained to hold times short enough such that the mode shapes of the atoms in states |a and |b do not change significantly, much more squeezing can be extracted from the system. These effects are illustrated in Fig. 5, which shows the best squeezing obtainable from a specific system with and without the π-pulse, as well as the extra squeezing that can be obtained by allowing a longer hold time as well as the π-pulse.

Squeezing in 2D and 3D
The simulations in the previous two sections were one dimensional, and show how having spatially constant density and using a π-pulse will improve the resultant squeezing. While these conclusions remain valid for simulations in two or three dimensions, adding each additional dimension increases the number of spatial modes in a given energy range. For BECs of different dimension at the same energy per particle, this leads to a degradation of squeezing, even when the potential (and the resulting atomic density) is constant across the bulk of the trap. This is demonstrated in Fig. 6, where we show the squeezing generated in 1D, 2D and 3D BECs for two   different parameter regimes. In both cases we take the BEC to be trapped in a box potential, which in the Thomas-Fermi limit leads to a constant density profile across the BEC. In Figure 6a we consider the case where U aa = U bb = 0, U ab = 0, while in Figure 6b we consider the case U aa = 0, U ab = U bb = 0. In both cases it is clear that number squeezing is best in lower dimensions, and becomes degraded or nonexistent as we move to a three dimensional system; the effect is worse if more than one nonlinearity is non-zero.
To find the trap geometries that produce the best squeezing, it is important to understand the root cause of this loss of squeezing in higher dimensional traps. The dependence on dimensionality is not explained by the processes described in section 3, which focus on the effects of inhomogeneity within the trap. Fundamentally, this degradation of the squeezing is due to nonlinearity-induced coupling between different momentum modes.
We can understand the origin and scaling of this degradation of squeezing by considering the response of the condensate to small fluctuations about the mean field. We do this using Bogoliubov theory [50], which considers the first order quantummechanical fluctuations about the mean field. To simplify the analysis we consider the case in which U aa = 0, U bb = U ab = 0 and the condensate has a spatially-constant density in a box of side lengths L x , L y and L z . In this limit the occupation of the non-zero momentum modes (the non-condensed fraction) in the |a internal state at the end of the pulse sequence is (a derivation is given in Appendix B) where the Bogoliubov mode frequency is ω k = ω 0 k (ω 0 k + 2χ aa n a ), and the free particle frequency is ω 0 k = k 2 /2M . The non-zero momentum modes can be expected to have an impact on the squeezing of the system when their occupation is a nonnegligible fraction of the total number of atoms. The fractional occupation of these modes is where the sum over k = 0 is taken over the available non-zero momentum modes k = (2π/L x )n xx + (2π/L y )n yŷ + (2π/L z )n zẑ with n x , n y and n z arbitrary integers that are not all zero. Figure 7 demonstrates the agreement between truncated Wigner simulations and the Bogoliubov theory for the non-condensed fraction for parameters where the occupation of these modes approach 15% of the total occupation of the system. The disagreement for larger non-condensed fractions occurs because the Bogoliubov theory used in deriving (34) is not number-conserving as it neglects the effect of the non-zero momentum modes on the condensate.
To minimise the damaging effects of the non-condensed fraction it is better to operate in a regime where the total occupation of the non-zero momentum modes N k =0 is small, which is also fortuitously the regime in which the Bogoliubov theory is a good approximation, so we can investigate in more detail. In 1D Eq. (34) is We approximate this sum by the integral where ∆k = 2π/L and L is the length of the 1D box condensate. This approximation will be valid in the limit that ω kmin τ hold ≪ 1, where k min = 2π/L as for small k the allowed values k = n x 2π/L will closely sample the region ξ = ω 0 k (ω 0 k + 2χ aa n a )τ hold < π where sinc 2 (ξ) is significant. Changing integration variables to Φ = ω 0 k τ hold and defining Λ a = n a χ aa τ hold , the non-condensed fraction (36) can be written as where In 2D and 3D the equivalent expressions are where and A is the area of the two-dimensional system and V is the volume of the threedimensional system. Figure 7 demonstrates the agreement in 1D between the integral approximation to N k =0 /N in (37), the summation expression in (34) and a truncated Wigner simulation. We wish to examine the scaling of the non-condensed fraction N k =0 /N with the physical size of the condensate. Other system parameters such as the total number N and the occupation of the |a and |b states at t 1 (n a and n b , respectively) are kept constant and chosen to maximise the squeezing of the equivalent two-mode model presented in Section 2. In particular, this fixes the values of θ, φ, and λ ij = χ ij τ hold . As the size of the condensate is scaled, ij /V will change, but we will vary τ hold to keep λ ij constant. Therefore Λ a = n a χ aa τ hold will also be constant. With this in mind, we can write τ hold in terms of the system size and the nonlinearity U aa . The non-condensed fraction can then be shown to scale with the physical size of the condensate as: The non-condensed fraction N k =0 /N thus scales differently with trap size in each of 1D, 2D and 3D. For one-dimensional traps, the non-condensed fraction increases with system size, whereas for 2D traps it is constant. For three-dimensional traps, the non-condensed fraction decreases with volume. This suggests that the best system to use will be a three-dimensional box potential, and the system size can be increased until the condensate depletion becomes negligibly small. Figures 8-10 show the non-condensed fraction predicted by the Bogoliubov analysis compared to that predicted by a full numerical simulation in one, two and three dimensions, as a function of the trap size. Also shown in the 1D and 2D cases is the maximum amount of number squeezing obtainable, as predicted by the numerical simulation. The full squeezing scaling curve is not shown in 3D, as for the smaller trap sizes the TWA method breaks down, and results are unreliable. For the larger trap sizes with small non-condensate fractions, however, the TWA simulation remained valid, and we obtained relative number squeezing of 12.7 dB. The traps are assumed to have a constant potential, which minimises the effects of the multimode issues described in Section 3. We see that the squeezing does indeed degrade at large trap sizes in 1D, but remains constant in 2D. For small trap sizes, the squeezing in 1D asymptotically approaches that of the ideal 0D system.
We note that the non-smooth nature of the squeezing prediction in 2D is due to the stochastic nature of the simulations, and the necessity of averaging over fewer quantum trajectories due to the very high number grid points in this regime. The nonsmooth nature of the non-condensate mode fraction in 3D, on the other hand, is due to the number of modes being both discrete and low in this regime, and counting only the non-condensed fraction in modes with an energy less than than µ, the chemical potential.
The simulation shown in Fig. 10 demonstrates that provided particle number is kept fixed, the number of accessible states for three-dimensional traps decreases with volume in 3D. This is because the energy of the ground state of the BEC (the chemical potential, µ) is reducing faster than the density of states at low energies is increasing, which in turn results in fewer states that can become dynamically occupied. This scaling of the estimated Bogoliubov fraction in 3D suggests that strong squeezing will be found in very large, weakly trapped condensates, provided the trapping potentials are box-like. Such box-like traps have already been demonstrated and shown to be capable of producing BECs [51].
Although systematic analysis of squeezing across a wide range of large 3D traps and interaction strengths is impossible due to the breakdown of the TWA method, fortunately the method is accurate when the highest squeezing is produced. While inter-state scattering lengths of most atomic species and states are not well characterised, extremely high squeezing can be found in systems with large, threedimensional box potentials, and asymmetry in the three scattering lengths.
As a demonstration of the plausibility of such a scheme, we carried out a full 3D simulation using the geometry of the optical box potential described in Ref. [51], using the same potentials (a cylindrical "box" with length 63 µm and diameter 30 µm) and condensate number (N = 10 5 87 Rb atoms). We set one s-wave scattering length to 5 nm and the others to zero: a 00 = 5 nm, a 10 = a 11 = 0. The associated nonlinearities are given by Eq. (2). This system produced number squeezing of 15dB, and an extremely low non-condensed fraction of 4 × 10 −4 . Both of these numbers are in agreement with the Bogoliubov calculations.

Conclusion
The large Kerr-type nonlinearities present between neutral atoms in BEC present the possibility of significant squeezing, but they also cause potentially complicated multimode behaviour when large BECs are used. Furthermore, the nonlinearities between multiple internal atomic states are not yet known for many atomic species, and often have sensitive dependence on external magnetic fields due to Feshbach resonances [52,53]. This paper analyses the requirements for producing useful squeezing via these nonlinearities for systems with large numbers of atoms.
We have presented an analytic solution to the two-mode model of atomic fields with a strong Kerr nonlinearity, allowing for a rapid estimate of the atomic number difference squeezing that can be obtained from these systems as scattering length data is obtained. In general, it predicts arbitrarily good squeezing in the two-mode limit. We then examined trap geometries in which this limit can be approached. Using a constant potential in the trap (a 'box' trap) allows the BEC density to be spatially uniform throughout the process. This is shown to alleviate the difficulties with choosing correct hold times and phases for the coupling pulses, and also alleviates the majority of the mode-matching issues for the final recombination of the two modes. If such a box potential is not available, the application of a π-pulse part way through the hold period can also improve mode matching by imparting the same modal dynamics to both states.
Even in the best case, where a box potential is used, the existence of multiple possible momentum modes allows the possibility of cross-mode scattering via nonlinear processes, which can degrade or totally remove any number squeezing. This degradation is visible only in full 3D simulations, and not lower dimensional approximations, due to its dependence on the number of available modes. It is worth noting that despite their difficulty, in multimode systems with large atom number, full 3D simulations need to be carried out to accurately determine the presence of number or spin squeezing.
Excluding these modes by making them energetically unfavourable can reduce the system to a single mode problem, making our analytic solution valid, but requires extremely tight trapping potentials. Such a regime also strongly limits the number of atoms that can be trapped, partly due to the rate at which the required trapping frequencies grow with atom number, and partly due to inelastic scattering rates.
Using a Bogoliubov analysis, however, it was shown that the coupling between different modes could be minimised in a second regime, namely the weak trapping limit for 3D traps. This latter regime is in fact ideally compatible with large atomic numbers, and shows the prospect of extremely high squeezing. A 3D TWA simulation using an experimentally demonstrated trapping potential yielded 15 dB of squeezing with 10 5 particles. In practice, the only limit to the number of atoms in such traps would be the stability of the trap, as the BEC would become increasingly sensitive to fluctuations as the energy per particle is reduced as the volume increases.

Appendix A. Derivation of two mode squeezing formula
In this Appendix we derive the equations for the number variances. We use the model described in Section 2, where we consider a two mode system with N particles, initially prepared as a coherent state with all the population in mode |a at time t = 0. A coupling pulse is then applied, assumed to be on resonance, so that the system evolves under the HamiltonianĤ until time t = t 1 . From time t 1 until time t = t 2 the coupling is turned off, and the system evolves under the Hamiltonian Finally, from time t 2 until time t = t 3 , the coupling is once again applied with a phase offset of φ relative to the first time, and the system evolves viâ We will work in the Heisenberg picture, and derive expressions forâ(t 3 ) andb(t 3 ) in terms ofâ(t 2 ) andb(t 2 ), which in turn can be expressed in terms ofâ(t 1 ) andb(t 1 ). We will show that the state of the system at time t 1 consists of separate coherent states in modes |a and |b , enabling us to calculate expectation values at time t 3 in terms of the known state at t 1 .
In order to have more compact expressions we will use the notation Beginning with the HamiltonianĤ 1 , during the period (0, t 1 ) modeâ will evolve asâ we obtainâ Similarly, we findb From the form of (A.8) and (A.9) it is clear that as both modesâ andb began in coherent states, they will remain in coherent states, with only their amplitudes changing. Specifically, at time t 1 the system is in a product of coherent states, which we denote |α, β with α = √ n a and β = −i √ n b with n a + n b = N . Next we consider evolution underĤ 2 . Taking modeâ first, we note that as it commutes with theb †b †bb term, and sinceâ †â †ââ commutes withâ †âb †b we havê where we have defined λ ij = χ ij (t 2 − t 1 ). Using (A.7) to moveâ 1 through the exponential, with some algebra one can show that To handle the cross-nonlinearity term in (A.10) we use the identity [54] e xâ †â f (â,â † )e −xâ †â = f (âe −x ,â † e x ) (A.12) to obtainâ Similarly, from permutation symmetry, we havê Finally, to obtainâ 3 andb 3 we use (A.8) and (A.9) with the phase factor attached to theb operator to obtain We are now in a position to evaluate the number variance of the system throughout the period of the final coupling. The number variances of the fields are defined as Making use of (A.15) and using the shorthand notation The remaining terms are not constants of motion, so we proceed by making use of (A.13) and (A.14). For â † 2b 2 we havê To calculate the expectation value of this expression, the creation and annihilation operators must be normally ordered. Making use of the identity [54] e xâ †â = where ω 0 k = k 2 /2M is the frequency of a free particle, and ω k = ω 0 k (ω 0 k + 2χ aa n a ) are the frequencies of the Bogoliubov modes: δb k (t) = exp (−iω k t) δb k (0). (B.5) The expectation value of the number of atoms in momentum mode k after the hold time is then found to be n k (t 2 ) = δâ † k (t 2 )δâ k (t 2 ) (B.6) = 1 + 4u 2 k (2u 2 k − 1) sin 2 (ω k τ hold ) δâ † k (t 1 )δâ k (t 1 ) − 4u k v k sin(ω k τ hold ) Im 1 − i2u 2 k sin(ω k τ hold ) δâ −k (t 1 )δâ k (t 1 ) + 4u 2 k v 2 k sin 2 (ω k τ hold ). (B.7) If the non-zero momentum modes are in a vacuum state at t = t 1 then this reduces to n k (t 2 ) = 4u 2 k v 2 k sin 2 (ω k τ hold ) (B.8) = [n a χ aa τ hold sinc(ω k τ hold )] 2 . (B.9) After the final Rabi coupling pulse the occupation of the non-zero momentum mode k is n a (k, t 3 ) = [n a χ aa τ hold cos(θ) sinc(ω k τ hold )] 2 (B.10) which is the expression given earlier in Eq. (33).