Dynamical formation and interaction of bright solitary waves and solitons in the collapse of Bose-Einstein condensates with attractive interactions

We model the dynamics of formation of multiple, long-lived, bright solitary waves in the collapse of Bose-Einstein condensates with attractive interactions as studied in the experiment of Cornish et al. [Phys. Rev. Lett. 96 (2006) 170401]. Using both mean-field and quantum field simulation techniques, we find that while a number of separated wave packets form as observed in the experiment, they do not have a repulsive \pi phase difference that has been previously inferred. We observe that the inclusion of quantum fluctuations causes soliton dynamics to be predominantly repulsive in one dimensional simulations independent of their initial relative phase. However, indicative three-dimensional simulations do not support this conclusion and in fact show that quantum noise has a negative impact on bright solitary wave lifetimes. Finally, we show that condensate oscillations, after the collapse, may serve to deduce three-body recombination rates, and that the remnant atom number may still exceed the critical number for collapse for as long as three seconds independent of the relative phases of the bright solitary waves.


Introduction
At low temperatures it is well known that the wave function of an interacting Bose-Einstein condensate can be described by solutions of the nonlinear Gross-Pitaevskii equation (GPE) for both repulsive and attractive atomic interactions [1]. In a one dimensional system a Bose-Einstein condensate (BEC) with intrinsic attractive interactions can form localised states for any number of particles. These are known as bright matter-wave solitons, and demonstrate particle-like behaviour in collision events [2]. A bright matter-wave soliton has been observed in a quasi-1D geometry at ENS through the tuning of atomic interactions from repulsive to attractive in a 7 Li condensate with a Feshbach resonance [3]. A bright soliton was also observed by Eiermann et al. by tuning the effective mass of a repulsive 87 Rb condensate to be negative with a 1D optical lattice [4].
Corresponding soliton solutions for condensates with attractive interactions are not possible in free space in three dimensions. However for a harmonically trapped BEC a condensate wave function can be found for particle numbers less than a critical value N cr [5]. Condensates with a particle number exceeding this value are unstable against collapse as demonstrated in the Bosenova experiment at JILA [6] and in the observation of the formation of quasi-1D soliton trains at Rice University [7]. Bright solitary wave (BSW) solutions to the 3D GPE with moderately tight harmonic trapping in two dimensions can be found for a condensate with attractive interactions as demonstrated by Parker et al. [8,9]. These are analogous to the soliton solutions of the 1D GPE, but are only self-trapped in one dimension. The formation and collisions of apparent multiple BSWs in condensate collapse have been observed recently by Cornish et al. [10].
In the Rice experiment, repulsive interactions between neighboring solitons have been inferred directly from the soliton trajectories [7]. In the JILA experiment, repulsive interactions between the BSWs were concluded indirectly as a consequence of the condensate atom numbers remaining in excess of N cr [10]. The long-time evolution of soliton trains [11,12,13] can qualitatively be accounted for by assuming that neighbouring solitons have a repulsive phase difference ∆ϕ. In mean-field theory this implies π/2 < ∆ϕ < 3π/2 [14]. However, the development mechanism of a repulsive phase-relation between the soliton-like structures in the experiments has not been explained.
Here we attempt to understand the physics of how 3D bright solitary waves form in condensate collapse, by using mean-field and quantum dynamical simulation methods incorporating the effects of three-body loss. We focus on the results of the experiment by Cornish et al. [10], who caused an initially stable 85 Rb condensate to collapse by switching the inter-atomic interaction from repulsive to attractive using a Feshbach resonance [15] whilst keeping the trapping potential on. At this point the number of atoms in the condensate exceeds the critical atom number for stability, N cr , which causes the macroscopic wave function to collapse [6,16,17,18], in the process losing atoms due to three-body recombination. As was found in earlier experiments, the onset of loss is sudden, and this allows the precise definition of the collapse time, t collapse [6], which is largely independent of the three-body recombination rate [6,16,17]. Following this primary collapse, it might be expected that repeated collapse processes and subsequent loss occur until the atom number decreases below the critical value, N cr [19]. In contrast, experiments [6,10] have reported remnant atom numbers, N remn , several times larger than N cr even after long evolution times. In Ref. [10] this is attributed to the creation of multiple bright solitons with mutual repulsive phase differences preventing any further collapse of the condensate. Parker et al. have studied the collisions of 3D bright solitary waves for the conditions of the experiment of Cornish et al. [10] in mean-field theory as a function of relative phase and velocity [8,20] . They have found that while the survival of the BSWs requires a repulsive relative phase at low velocities, for sufficiently large relative velocity the BSWs survive single collision events independent of the relative phase. They also found that a ∆ϕ in the repulsive range can inhibit collapse induced destruction of the soliton pair for many oscillation periods involving more than 40 collisions [8,20].
In this article we models the initial collapse and structure formation of the experiment of Cornish et al. [10] by solving the mean-field Gross-Pitaevskii equation in three dimensions with realistic parameters. We find that the dynamical creation of solitons as in Ref. [10] does not favour repulsive ∆ϕ. For the Rice experiment [7] the absence of fixed repulsive phase-relations in Gross-Pitaevskii (GP) simulations of soliton creation was pointed out in Ref. [21]. For the JILA experiment [10], the fact that the initial GP wave function has even symmetry prevents repulsive phase relations in the final state for an even number of BSWs as mean-field theory preserves this symmetry. In this situation the central two members of the train must therefore have ∆ϕ = 0. In order for this statement to be false, there must be some form of symmetry-breaking mechanism.
To go beyond mean-field theory, we account for quantum field corrections to the GPE by modelling the collapse using the truncated Wigner approximation (TWA) [22,23,24] which incorporates the effects of quantum noise. However, we find that for conditions close to the experiment [10] the relative phase between the central two dynamically formed BSWs is still near zero. An example of BSWs with attractive phase relations found within truncated Wigner quantum field theory is shown in figure 1.
In order to better understand the effects of quantum noise, we simulate a 1D analogue of the JILA experiment [10], and find that quantum fluctuations have a significant effect on the resulting soliton interactions. We perform controlled collisions between solitons of varying relative phase in 1D, and find that quantum noise tends to render soliton interactions more repulsive than the corresponding mean-field solution, irrespective of their relative phase. Three-dimensional (3D) simulations indicate that N remn can remain above N cr for long times in the presence of quantum and thermal fluctuations even if the interactions of the corresponding mean-field BSWs are attractive. However, we find the effects of quantum noise tend to destroy the BSW structures for all initial relative phases, in contradiction to the results of mean field simulations [8,20] and experiment [10].

Three dimensional BEC collapse
We first model collapsing condensates using the time-dependent Gross-Pitaevskii with phenomenological two-and three-body loss terms [19,16,25]. Here the condensate wave function Φ(r, z) is written in cylindrical co-ordinates. We assume 85 Rb atoms with mass m = 1.411 × 10 −25 kg, interaction strengh g = 4πh 2 a s (t)/m, where the scattering length a s (t) varies in time. We define a s (t = 0) = a ini ≥ 0 and a s (t > 0) = a collapse < 0. Further, ω r /2π = 17.3 Hz, ω z /2π = 6.8 Hz [10]. The threebody loss rate, K 3 , has not been precisely determined for the conditions of Ref. [10]. For this reason, and for numerical necessity, we vary K 3 throughout this work. We take K 2 = 1.87 × 10 −20 m 3 s −1 [26]. Our initial condition is the solution to the time-independent Gross-Pitaevskii equation with a s = a ini found via imaginary time evolution. To reduce the computational demands of the problem we use a spatial grid size that does not fully accommodate burst atoms ejected during collapse [6]. As these are not the focus of our studies, we employ absorbing potentials near all numerical grid edges §.
In figure 2 we compare our numerical solutions of (1) with the experimental data from Ref. [10] for the case a collapse = −11.4a 0 and a ini = 9a 0 , where a 0 is the Bohr radius. For these initial conditions the experiment observes the formation of two BSWlike structures. We solve (1) in cylindrical co-ordinates using a Fourier transform based adaptive step-size Runge-Kutta algorithm. The collapse time [6] for this scenario is t collapse = 12 ms. Near this time the peak-density rises dramatically with a maximal value that scales as U/K 3 [27]. Directly after the collapse burst atoms are ejected from the condensate . A so-called burst focus [6] takes place at a later time, when all ejected atoms that did not hit the absorber have classically completed half of an oscillation § The absorbers have the generic form , where ∆x is the distance to the nearest grid edge, d the absorber width, θ the Heavyside step function and γ is the absorber strength.
In our simulations most of these are eliminated with the absorbing potential. in the radial potential and have returned to the trap symmetry axis. This happens at t bf = t collapse + π/ω r ≈ 41 ms. We see the onset of the burst focus around 35 ms. The burst focus develops into two axially separate clouds, in agreement with experimental observations [28]. In our simulations these become precursors of a BSW pair. It takes until t = 72 ms before these clouds develop two distinct BSW-like shapes. The distance by which the atoms are axially displaced by the collapse, and hence the width found in figure 2, increases for lower K 3 due to the larger interaction energy that is liberated during the collapse. Beyond t ∼ 50 ms the absorbing potentials become relevant for the dynamics.
It is generally expected that the atom bursts have a strong uncondensed component [29,30], which in principle necessitates a theoretical treatment beyond the mean-field for t > t collapse . The numerical results in figure 2 obtained using GP theory are hence approximate. Nonetheless, they describe the initial evolution of the experimentally observed full-width half-maximum (FWHM) rather well. Furthermore, the amplitude of the condensate oscillations depends strongly on the three-body loss rate K 3 and the agreement between theory and experiment is better for smaller K 3 . We suggest that this may allow for a better experimental measurement of K 3 than has been possible so far [26]. The reliability of this approach would depend on the influence of the uncondensed burst atoms. Most importantly, we confirm that the BSWs are created in-phase as required by symmetry.
We now consider the effect of quantum fluctuations on the collapse. It has been suggested that quantum fluctuations can imprint a staggered phase pattern onto the condensate cloud during the development of instability [12]. To investigate this we employ the truncated Wigner approximation [18,22,23,24,31] of the full quantum evolution of the condensate. In this approach quantum field effects enter through the rigorous addition of noise to the initial condensate wave function Here Φ(x, t = 0) is the condensate initial state as used in mean field simulations. α(x) is termed stochastic field, ϕ n (x) are the basis elements employed numerically and η n is complex Gaussian random noise with correlations η n η m = 0, η * n η m = δ nm . In essence, this noise mimics the effect of vacuum fluctuations. The stochastic field α(x) then evolves according to (1). For a more accurate treatment of the three-body losses we optionally include dynamical noise terms in our simulations [32]. The theoretical and numerical methods employed here are described in Ref. [18]. The full truncated Wigner equation of motion for the stochastic field is [32] dα with dynamical noise dξ(x, t) that has correlations dξ( Here, the cylindrical symmetry is broken by the noise in the initial state (2). For a consistent treatment of the noise the propagation algorithm uses the harmonic oscillator basis [31,18]. For the problem to remain computationally tractable in this basis we increase K 3 compared to the values in figure 2 (see e.g. figure 1), thus sacrificing direct correspondence to the experiment. We also apply the quantum treatment to a low energy subset of the full mode space in order to justify the TWA. Within these constraints our simulations show no sign of a strong preference for repulsive ∆ϕ within a large range of a collapse . Therefore, the results of our 3D simulations question the assertion that BSWs originating from a condensate collapse as in [10] are created with a repulsive phase relation between neighbours.

One dimensional collapse with quantum noise
While the investigation of quantum effects in three dimensions for realistic experimental parameters remains computationally intractable ¶, we can potentially gain insight into the effects of quantum noise using a one dimensional model. To do so we integrate out ¶ The strong localisation of the condensate at the moment of collapse necessitates a fine numerical grid and the ejected burst atoms require a large spatial grid domain. Additionally, the experimental conclusions stem from BSW behaviour at long time. Altogether these requirements mean that quantum simulations for the precise experimental parameters would be extremely demanding. We start from an initial total number of N ini = 8000 atoms in the trap ground state (a ini = 0), which is then exposed to an attractive nonlinearity with a collapse = −40a 0 . A collapse occurs around t = 5 ms. Other parameters are ω z /2π = 6.8 Hz, a r = 1 µm (ω r /2π = 100 Hz) and K 3 = 2 × 10 −40 m 6 s −1 . Bright (dark) regions correspond to high (low) atomic density. We truncate the density scale at 10% of its peak value to improve visualisation. the transverse dimension r of (1) for K 2 = 0, using the approximation that it remains in the ground harmonic oscillator state [33]. The 1D equation becomes where g 1D = g/(2πa 2 r ) andK 3,1D = K 3 /(3π 2 a 4 r ) with a r = h/(mω r ). Figure 3 shows typical simulations of a one dimensional collapse + . Figure 3(a) is the mean-field GPE solution. Close inspection of the high density regions reveals solitonlike structures for most times. Their number, however, does not remain conserved at all times but varies after subsequent collapses. We find that no fixed repulsive phaserelations form and that the overall system exhibits breathing oscillations with repeated increase of density accompanied by three-body loss whenever solitons closely approach one another.
Applying the dimensionality reduction used above to (3) yields the one dimensional stochastic differential equation for the truncated Wigner framework Figure 3(c) shows a trajectory of a 1D collapse simulation using the TWA with initial and dynamical noise, implemented only on the central half of Fourier space to avoid aliasing and keep the noise density much lower than the condensate density [34]. A qualitative difference in the dynamics in the presence of vacuum fluctuations [figure 3(c)] can be noted compared to the mean-field simulations: The overall dynamics of the solitons is more repulsive leading to a decrease in the overall atom loss due to lower densities at the moments of solitons closest approach. This behaviour is typical for all realisations of noise. In the following section we investigate the phenomenon of noise induced repulsion more closely.

Quantum soliton collisions
While we do not find any evidence that quantum noise facilitates the creation of solitons (in 1D) or BSWs (in 3D) with repulsive phase relations during the BEC collapse, our 1D model suggests that the noise strongly affects the collisional dynamics of the observed solitons. Strong effects of incoherence on the interaction properties of solitons are for example known from nonlinear optics [35,36].
To understand this in more detail, we focus on the collision stage and investigate the effect of quantum noise on the interactions of analytically prepared bright solitons. + We refer to a one dimensional collapse as the situation where a narrow density spike has developed and then exploded due to the onset of strong losses. Unlike in 3D, a 1D setting does not develop a mathematical density singularity. For this we make Eq. (4) dimensionless with energy, time and length scales given byhω r , ω −1 r and a r respectively, and setK 3,1D = 0. The initial state of these 1D simulations is a superposition of a pair of separated bright solitons, individually multiplied by a phase factor e iknx such that the solitons propagate towards one other with equal speed. Explicitly with A = 20, γ = |g|A, p 1,2 = ±10, k 1,2 = ∓1 and ϕ 1 = 0. Within mean-field theory the character of the ensuing collision is determined by the independent parameter of relative phase ∆ϕ = ϕ 2 . The density evolution for ∆ϕ = 0 (∆ϕ = π) is shown in figure 4(c) [ figure 4(d)]. The most significant difference between collisions with opposite relative phase occurs at the moment of closest approach of the solitons, which is labelled by t * . The solitons merge into a single peak for ∆ϕ = 0, while for ∆ϕ = π they preserve a double peak structure. Figure 4(a,b) demonstrates a Gaussian fit to the overall density profile at the time of collision, giving a quantitative criterion to distinguish the attractive and repulsive situations. In the attractive case the width of the Gaussian, σ fit , is given by the single narrow solitonic peak [ figure 4(a)], while in the repulsive case it is set by the separation between the solitons, which is much larger [ figure 4(b)]. A second distinguishing feature is the peak density at the moment of closest approach, which is also the peak density during the collision n peak = max x,t |Φ(x, t)| 2 . The dependence of both these indicators on the relative phase ∆ϕ is shown in figure 4(e,f). In the mean-field treatment they behave as expected, i.e. there is a clear increase in σ fit around ∆φ = π. However, the situation becomes completely different if the effects of vacuum fluctuations on these collisions are accounted for.
We find that in the presence of noise the relative phase between the solitons in the mean-field component of the condensate wave function loses its significance. When averaged over 200 realisations of the atom-field with different realisations of noise, the overall peak of the density, n peak , and the width of the Gaussian fit, σ fit , are no longer distinct for in-phase and out-of-phase solitons * . Instead, both variables are independent of the initial relative phase. For all relatives phases the width of the Gaussian fit corresponds closely to that of the repulsive mean field simulations.

Three dimensional bright solitary wave collisions
Finally, we investigate whether observations regarding soliton interactions in the one dimensional simulations carry over to three dimensional collisions of bright solitary waves. As stated earlier, full simulations of the entire collapse experiment are numerically intractable. It is possible, however, to study binary collisions of three dimensional BSWs using the Gross-Pitaevskii equation [8,20]. We extend the work of Ref. [8,20] by performing the same simulations incorporating quantum noise in the truncated Wigner approximation, as well as the effects of three-body loss. It has previously been shown in Ref. [20] that collisions of bright solitary waves in attractive BECs depend on the collision velocity. In particular, at low velocities a phase difference ∆ϕ = π can increase the possible number of repeated collisions before a significant degradation of the soliton shape. In this section we show that quantum and thermal noise have the opposite effect, reducing the number of collisions with preserved soliton shape.
Our initial binary soliton state is chosen to most closely fit our findings of section 2. Thus, we employ K 3 = 2 × 10 −40 m 6 s −1 , an initial total atom number N ini = 3500 and a soliton separation of about d = 32 µm for a collapse = −11.4a 0 . We find the initial soliton state using a conjugate gradient method, while neglecting the axial trap. For * In detail, we first obtain σ fit and n peak from each trajectory and then take the average. Although, this does not have a direct interpretation within the TWA it provides an indication of how individual experimental realisations might behave and serves as a measure of the effect of quantum fluctuations on the soliton collisions.  Fig. 4 of Ref. [8] for corresponding data without three-body loss). (a) Total number of atoms remaining in the trap: meanfield theory (solid), TWA with initial noise (dashed) and TWA with dynamical noise (dot-dashed). Results for different relative phases are shown: ∆ϕ = 0 (black, lower curves), ∆ϕ = π (red, upper curves). Panels (b-g) display atomic density along the zaxis with upper (lower) panels showing the evolution at early (late) times. We present mean-field simulations (b-e) and TWA with initial noise (f,g). The BSWs are in-phase ∆ϕ = 0 (b,c) or out of phase ∆ϕ = π (d-g). For BSWs within the single trajectory TWA we enforce the relative phase before quantum fluctuations are added. Darker areas correspond to higher densities. Note that for better visibility a nonlinear grey scale that exaggerates low density regions is used.
the dynamical simulation to replicas of it are symmetrically placed around the origin on the z-axis with a relative phase ∆ϕ. The axial trap is then switched on and the solitary waves are allowed to propagate towards the origin to collide [8,20].
We determine the remnant atom number of such a scenario after the evolution time of three seconds (as in the experiment of Cornish et al. [10]). Three theoretical approaches were employed for the dynamical evolution: mean-field theory, truncated Wigner formalism with initial noise only and TWA with the addition of dynamical noise [18]. All simulations use the harmonic oscillator basis. The results are summarized in figure 5. It can be seen that for both initial phases the atom number remaining after 3 seconds exceeds the critical number, which for this scenario is N crit ≈ 2800.♯ Furthermore, in the presence of three-body loss [figure 5(c,e)] the relative phase ∆ϕ = π allows oscillations of the BSW pair in the trap for longer times than ∆ϕ = 0, as was shown without the inclusion of three-body loss in Ref. [20].
The inclusion of vacuum fluctuations via the truncated Wigner simulation does not significantly alter the remnant atom numbers around three seconds. It does, however, strongly affect the spatial BSW evolution. As shown in figure 5(g), even for ∆ϕ = π the two clearly separated clouds no longer survive at long evolution times, in contrast to mean-field studies [figure 5(e)]. This also disagrees with the experimental results, which did observe BSW-like structures after three seconds and more than 40 collisions ♯ The critical number is N crit = kā/|a s |, with k = 0.55,ā = h/mω,ω = (ω 2 r ω z ) 1/3 [10].
between the two wavepackets [10]. Therefore, from our 3D simulations results we cannot conclude that the quantum noise renders BSW collisions more repulsive -in fact quantum noise appears to destabilise the case of ∆ϕ = π. We note, however, that this qualitative difference between 3D and 1D results could stem from the reduced amount of noise that we are forced to apply in 3D due to numerical reasons. † † The apparent agreement between mean field theory and the experimental results in [20] is compelling; however there is no obvious mechanism for the development of the π phase difference, and we can offer no explanation for why a calculation that we expect to be more accurate than GP simulation gives results that so obviously differ from observations in the lab.

Conclusions
We have studied the short and long time dynamics of the collapse of attractive Bose-Einstein condensates, applying 3D and effective 1D models that include three-body losses and quantum noise. We find that due to the kinetic energy liberated, the condensate exhibits violent oscillations in the trap after the collapse, with an amplitude that strongly depends on three-body loss rates, K 3 . We propose that the sensitivity of these oscillations after the collapse may improve measurements of K 3 . After the collapse highly excited BSW-like structures form from the remnant condensate. These also exhibit violent dynamics and disappear, and reform multiple times throughout the evolution. The relative phase between dynamically formed BSWs does not usually take repulsive values (π/2 < ∆ϕ < 3π/2) as has been previously postulated [10]. However, we show in one dimension that the addition of quantum fluctuations can render interactions repulsive, even for solitons with a phase relationship for which standard Gross-Pitaevskii theory predicts them to exhibit attractive interactions. Our three dimensional studies of interactions with quantum noise do not show a similar effect at our level of approximation; instead they suggest that quantum noise results in a shorter BSW lifetime. Our results raise questions about the interpretation of the JILA collapse experiment [10] in terms of the creation of BSWs with mutually repulsive relative phases. In light of this we suggest that the system warrants further theoretical and experimental study. In particular, it will be interesting to compare the results of the present studies with more direct treatments of quantum effects in soliton or BSW collisions. Experimentally, one could aim to measure the soliton phase directly [37] rather than inferring it from collisional behaviour. † † Care must be taken when interpreting these simulations as quantum field evolution in the sense of the TWA. Due to the application of noise to a restricted set of modes and, in particular, the long evolution times, the TWA is not strictly justified. Nonetheless, the noise present in the simulation can closely mimic the effect of hot uncondensed atoms that are known to be present after the condensate collapse. In this sense we would expect the simulations with noise to offer a more realistic description than a pure GP treatment.