Nonequilibrium quantum dynamics of atomic dark solitons

We study quantum dynamics of bosonic atoms that are excited to form a phase kink, or several kinks, by an imprinting potential in a one-dimensional trap. We calculate dissipation due to quantum and thermal fluctuations in soliton trajectories, collisions and the core structure. Single-shot runs show weak filling of a soliton core, typically deeper solitons in the case of stronger fluctuations and spreading/disappearing solitons due to collisions. We also analyze a soliton system in an optical lattice that shows especially strong fluctuation-induced phenomena.


Introduction
Improved technology to manipulate ultra-cold atomic gases by optical lattice potentials [1], atom chips [2], magnetic fields, etc., has provided tools to trap atoms in strongly confined geometries and to adjust the dimensionality of the system. Such many-particle systems only interact weakly with their environment and, over short time-scales, they may in many situations be approximately considered as isolated and studied by models based on unitary quantum evolution. Moreover, the many-body atomic states may be engineered and controlled at high accuracy and the detailed dynamics of atoms in response to external perturbations can be experimentally investigated.
In a tightly-confined one-dimensional (1D) tube-potential quantum fluctuations of bosonic atoms can become strong [3,4,5]. The quantum effects may be further enhanced, e.g., by applying an optical lattice potential along the axial direction of the trap [6,7,8]. Stochastic phase-space methods provide a useful description of dissipative bosonic atom dynamics in 1D traps due to quantum and thermal fluctuations of the atoms when the particle number is not too low [9,10,11,12]. The approximate quantum dynamics of atoms can be modelled in the truncated Wigner approximation (TWA) [13,14,15,16] by unitary quantum evolution where the quantum statistical correlations of the initial state may be accurately synthesized for a variety of quantum states in the Wigner representation. The method can incorporate a very large phase-space, with a large number of degrees of freedom, in which case the atom dynamics, e.g., does not need to be restricted to the lowest energy band of an optical lattice. Dissipative dynamics emerges from a microscopic treatment of the unitary quantum evolution without any explicit damping terms and there is no need for the frequently problematic separation of quantum dynamics to 'system' degrees of freedom and 'environment' [17]. In addition, it is possible to study excitations of the system far from thermal equilibrium, quantum expectation values of various experimental observables and outcomes of singleshot measurements. Previous unitary TWA simulations [10] were qualitatively able to produce the experimentally observed damping rate of a dipolar centre-of-mass motion of bosonic atoms in a very shallow, strongly confined 1D optical lattice [7]. It was found that, due to quantum fluctuation induced momentum uncertainty, a small atom population reaches a critical velocity, representing an onset of a dynamical instability in the corresponding classical system. The calculated damping rate was approximately proportional to the atom population in the dynamically unstable velocity region.
Here we study dissipative quantum dynamics and relaxation of a 1D bosonic atomic gas that has been excited to form a moving phase kink by an optical imprinting potential. Phase kinks have been experimentally imprinted in this way in 3D atomic Bose-Einstein condensates (BECs) by imaging the atom cloud, e.g., through an absorption plate [18,19] or by using optical holograms [20]. We recently showed how such a setup could be used to study dynamics of dark solitons in situations where quantum fluctuations are important [21]. The integrability of the soliton dynamics is broken by a trap and quantum and thermal fluctuations promote sound wave emission that may dissipate and eventually equilibrate with the soliton. Numerically tracking soliton coordinates in individual stochastic realizations provided us with a tool to calculate quantum mechanical expectation values and uncertainties of the soliton position. We found, e.g., a surprising result that the quantum expectation value of the speed of a soliton is reduced due to enhanced quantum fluctuations, as a result of the nonlinear dependence of the soliton speed upon its phase distribution. Single-shot runs in an optical lattice revealed effects of dynamical instabilities, such as jittering oscillatory motion, splitting and disappearing solitons. In this paper, we study a detailed structure of a quantum soliton, its quantum statistics and dissipative dynamics due to quantum and thermal fluctuations in fairly large atom-number systems. In particular, we consider the effects of dissipative dynamics on a soliton core structure, predicting a distinct difference between the filling behaviour of a soliton core in single-shot experimental runs and in the quantum expectation value of atom density found by averaging over many experimental runs. Single-shot realizations show only weak filling of a soliton core at later times and typically deeper solitons for the case of stronger quantum fluctuations. On the other hand, quantum expectation values for the atom density and pair-correlation functions are smeared out due fluctuation-induced drifting of solitons along the trapthis is particularly visible for very slow solitons in an optical lattice system we consider. We also simulate quantum effects of collisions of soliton pairs and relaxation in a system of several solitons that may behave chaotically in the corresponding classical case. As in [21] quantum and thermal fluctuations are synthesized within the truncated Wigner approximation in the quasi-condensate description. Dark solitons in nonlinear optical fibers were previously modelled using TWA [22]. In atomic BECs, finite temperature dissipation has now been theoretically shown to manifest itself in an increasing oscillation amplitude of solitons in a harmonic trap [23,24,25,21]. Quantum properties of dark solitons [17] have recently started attracting considerable interest, including studies, e.g., of soliton core structure [26,27,28,29,30,31,32], phase kinks in a uniform space close to the Mott transition [33], and soliton statistics [34].

Dark soliton experiments
Dark solitons have been prepared in many experiments, and their subsequent dynamics probed. The first dark solitons in ultra-cold atomic Bose gases were generated by imprinting a sharp phase jump across the gas by optical potentials [18,19], or by imprinting density defects using slow light pulses [35,36]. Subsequently, dark solitons have been created, e.g., by merging of two BECs [37] and by superfluid flow around a barrier [38]. In typical experiments the atoms were confined in 3D traps, or in 1D traps with only a weak radial trapping frequency, and the systems are accurately described by the classical Gross-Pitaevskii equation (GPE). In the phase imprinting method, a soliton is generated by applying a constant 'light-sheet potential', of value V φ , to a part of the atom cloud, for time τ [18,19,20]. The light-sheet potential is obtained by shining a far-detuned laser, e.g., through an absorption plate, so that the resulting dipole potential for the atoms exhibits a sharp edge. In the classical GPE limit, the potential imprints a phase jump generating a dark (or grey) soliton (phase kink) with velocity v and the density dip at the phase kink (soliton core) with the minimum density ρ s [39] v/c = cos(φ c /2), where c is the speed of sound and ρ b is the density of the atomic background. The soliton is stationary (dark) for φ c = π, with a zero density at the kink. Other phase jumps produce moving (grey) solitons, with non-vanishing densities at the soliton core, so that |v| → c for φ c → 0. Soliton imprinting is accompanied by a density (sound) wave moving in the opposite direction to the soliton with speed approximately equal to c [18], created as a by-product of the perturbation of density by the imprinting potential. More complex defect structures may generally be imprinted on superfluids using optical phase holograms to shape laser fields, so that, via coupling with matter waves, the light acts as a hologram to shape the BEC [40]. A spatial light modulator was experimentally used to generate a multi-step potential and to prepare multiple phase kinks [41]. The potential (where θ is the Heaviside step-function) imprints phases φ j at position x j . The potential V φ j may be negative, such that V φ j = −Vφ j imprints a phase of 2π −φ j , corresponding to a soliton travelling in the direction of the negative x-axis. Hence the initial positions and velocities of solitons may be accurately controlled. Other experiments have created multiple solitons by colliding two BECs [42,43].
Subsequent to their formation, single solitons in a harmonic trap oscillate at angular frequency ω/ √ 2, where ω is the frequency of the trap [44]. Although many soliton experiments were performed in harmonic traps, these had a 3D character [45], such that solitons could decay into a hybrid of vortex lines and rings ('snake' instability) [46]; consequently, lifetimes of solitons were limited to less than one oscillation period. However, some experiments suppressed the snake instability by increasing the radial trapping potential [20,43], permitting observations of oscillations at the predicted frequency [20]. Recently, experiments probed the interactions between dark solitons [20,41,43] or between dark and dark-bright solitons [20]. Collisions between two dark solitons are accompanied by position-shifts [47,48], which change the soliton oscillation frequency [43]. During experiments, solitons repeatedly collided up to seven times without losing their integrity [43].
Although some soliton experiments were sufficiently 1D to suppress snake instability, residual radial dynamics affected the oscillation frequencies by approximately 3% [43]. In more tightly confined 1D traps the radial density fluctuations may be completely suppressed. For example, in an atom transport experiment, a 2D optical lattice divided a BEC into an array of decoupled 1D tubes [7] with radial and axial trap frequencies in each tube ω r = 2π×38kHz and ω r = 2π×60Hz, respectively, and about 70 atoms in the central tube. A wide variety of 1D trapping schemes are also provided on atom chips using electromagnetic fields [2]. In a tightly confined 1D limit, bosonic atoms become more strongly interacting at low atom density [49], characterized by the effective interaction parameter γ int = mg/h 2 ρ, where ρ is the 1D atom density. For increasing values of γ int , enhanced phase fluctuations destroy the long-range phase coherence across the quasi-condensate, but density fluctuations remain suppressed-eventually the atoms reach the fermionized Tonks-Girardeau regime [6,4].

Truncated Wigner approximation
In our analysis, quantum and thermal fluctuations of atoms are approximately included within TWA [13,14,15,11,16]. Our TWA formalism and noise generation closely follows [11], except that here we fix the total atom number and use a quasi-condensate description for quantum statistical atom correlations [21]. We replace the quantum field operators (ψ,ψ † ) by the classical fields (ψ W , ψ * W ), governed by a 1D unitary evolution where the interaction strength g = 2hω r a, the s-wave scattering length a, the total number of atoms N tot , and V = V ext + V i is the external 1D potential including the time-dependent imprinting potential. Equation (5) formally coincides with the classical GPE, but here quantum and thermal fluctuations are sampled in the initial state by an ensemble of stochastic fields ψ W (x, 0) in the Wigner representation. Initially, before the phase imprinting, the atoms are assumed to be in thermal equilibrium in a harmonic trap, or in a combined harmonic trap and optical lattice. In a tight 1D trap we expect significant phase fluctuations that we found particularly important to phase kinks [21] and in order to incorporate these we introduce a quasi-condensate description for the field operator [50] ψ( The density δρ(x) and phaseφ(x) operators may be written in the Bogoliubov-type expansion [50], requiring δρ/ρ 0 and |δl∆φ| to be much less than one, where δl is the spacing on the numerical grid on which we calculate the operators and ∆φ is the gradient of the phase operator across one gridpoint. Thus (for j > 0) where −h 2 2m where µ is the chemical potential. Here ρ 0 =N 0 |ψ 0 (x)| 2 and ψ 0 (x) is the ground state wavefunction withN 0 particles. We fix the total atom number to N tot and allow the ground-state and excitedstate populations, N 0 and N nc , to fluctuate. In the Bogoliubov approximation, the expectation value of the excited state population (j > 0) In order to sample the stochastic initial state of TWA with the correct quantum statistics, the operators (α † j ,α j ) in (7) and (8) are replaced by complex variables (α * j , α j ), which may be generated from the Wigner distribution for harmonic oscillators with energy ǫ j at temperature T [11], resulting in a stochastic Wigner representation (ϕ W (x), δρ W (x)) of phase and density operators. The stochastic initial state for the time evolution (5) then reads where ρ 0,W will be specified shortly. Due to the Wigner representation, the stochastic variables correspond to symmetrically ordered expectation values of operators (α † j ,α j ) and, e.g., the expectation values of the mode occupations differ from those calculated in the Wigner representation by half an atom per mode: α † jαj = α * j α j W − 1/2, where the notation ... W indicates the expectation values are calculated in the Wigner representation. The excited state population for each run therefore reads fluctuating around the correct expectation value (11). We set the ground-state population for each stochastic realization as N 0 = N tot − N nc , so we can set for the initial state in each run ρ 0, In simulations we vary the ground-state depletionN nc /N tot . At temperature T = 0 we keep the nonlinearity fixed at N tot g = 100hωl, but adjust the ratio g/N tot where l = (h/mω) 1/2 . This is tantamount to varying the effective interaction strength γ in . We also study the effects of thermal depletion by varying T , by fixing both N tot and g and setting g/N tot to be sufficiently small so that quantum fluctuations are not dominant.
After evolution of the stochastic fields ψ W , we extract physical quantities from the ensemble by calculating their quantum mechanical expectation values. Due to the symmetric ordering of the expectation values in the simulation data, we need to transform these to the normally ordered expectation values, corresponding to typical physical measurements [9,11]. We find that the soliton position measurements considered in this paper are unchanged in the normal-ordering of phonon-mode occupation numbers, but that densities, pair-correlations and number fluctuations are quantitatively affected by the ordering. Hence the soliton position coordinates in individual experimental runs can be accurately extracted from the Wigner density |ψ W | 2 .

Damping and dissipation in oscillatory dynamics
We study the imprinting of single phase kinks and their subsequent dynamics in a harmonic trap in TWA. In our previous study [21] we showed how quantum expectation values and uncertainties could be calculated for the dynamical variables of the soliton. We ran ensembles of hundreds of realizations (e.g. 400), and numerically tracked the soliton position coordinates x(t) from |ψ W | 2 in each run. This allowed us to calculate, e.g., the quantum mechanical position expectation value x and uncertainty δx = x 2 − x 2 . Based on our findings, we made, e.g., the following observations that are relevant to our present study: (i) soliton trajectories are damped by both quantum and thermal fluctuations, increasing the amplitude of soliton oscillation; (ii) quantum mechanical soliton position uncertainties increased as a function of time for systems with large quantum depletionN nc /N tot -this increase was due to the initial uncertainty in soliton velocity and deviations in the soliton trajectories as they interacted with sound waves; (iii) one of our most dramatic findings was that the quantum expectation values for the speed of a soliton were lower in systems with large quantum depletion which is attributable to a broad distribution of phase jump values across the kink caused by quantum fluctuations. We mapped the phase kink of each trajectory in the ensemble to a corresponding speed using the classical formulae (3) and found that due to the negative curvature of the soliton speed | cos(φ/2)|, a symmetric phase distribution always leads to a lower quantum expectation value for speed | cos(φ/2)| than the speed given by the phase expectation value | cos( φ/2 )|. Consequently, a greater phase uncertainty leads to a lower expectation value for the soliton speed.
Here we evaluate in more detail the dissipation in soliton dynamics due to quantum and thermal fluctuations. We also study the structure of the soliton core and phase in individual stochastic realizations of TWA that represent possible outcomes of singleshot experiments. We then compare the quantum statistics of a soliton core in singleshot outcomes to quantum expectation values of the atom density which would be obtained in an experiment by averaging the density over many runs. We simulate the optical imprinting of phase kinks according to (5) in the harmonic trapping potential V exp = mω 2 x 2 /2, with the imprinting potential (4), for V φ = 4166.7hω, applied over a period of τ = 4.78 × 10 −4 /ω. Such a potential imprints a phase kink of φ c = 2 in the classical (GPE) limit. Due to the stochastic sampling of quantum and thermal noise in the initial state of each realization of the Wigner wavefunction, the actual phase jump values across the soliton core fluctuate from run to run, resulting in a variation in the shape of the solitons and in their initial velocities. After the imprinting the solitons oscillate in the harmonic trap.
Previously [21], we observed that soliton oscillations in a harmonic trap exhibit damping; since dark solitons have negative mass and kinetic energy, damping of a soliton increases its amplitude of oscillation. We track each soliton's trajectory by locating the local density minimum around the phase kink in ψ W for each realization. We fit individual trajectories with the curve x(t) = f (t) exp(−γt), where f (t) is a sinusoid and γ < 0. We separately study the damping due to the ground-state depletion resulting from quantum fluctuations (at T = 0) and from finite-temperature atoms (for parameters for which the corresponding T = 0 damping is weak). We only fit trajectories for parameters for which fluctuations are not too large, such that the trajectories are sufficiently close to sinusoidal. Figures 1 (a) and (b) show that there is a faster-thanlinear increase in damping with depletion in both cases.
The source of this damping behaviour can be deduced by consideration of individual trajectories in ensembles with large depletion (shown in figures 2 and 3). Strongly dissipative behaviour occurs when energy is exchanged between the solitons and sound waves. This is particularly clear for large quantum depletion shown in figures 2(a)-(c). We find that at finite temperature large thermal depletion has an associated sound wave background [figure 3(d)-(f)], but the soliton trajectories are less perturbed by the sound waves than in cases of large quantum depletion. We contend that the breaking of integrability by the harmonic trap causes the energy to disperse amongst the excitations. Solitons with relatively less negative kinetic energy are more prone to dissipation, and the damping rate is greatest in these trajectories.
The limit g/N tot → 0, for fixed gN tot , corresponds toN nc /N tot → 0 and we expect the Bogoliubov approximation to become increasingly accurate (if the ground state and the excited state populations are not solved self-consistently,N nc is approximately constant and N tot → ∞ [51]). Surprisingly, in the limit that g/N tot → 0,γ converges to a more negative value than that found in classical GPE simulations [ figure 1(b)]. The classical soliton exhibits damping due to interaction with sound waves that originate from the imprinting process. Such interactions may represent dynamically unstable processes that require very weak numerical noise to be triggered which may be absent in GPE. Alternatively, the soliton trajectory may represent a state with non-classical statistics that does not converge to the classical GPE value, similarly to the atom number-squeezed states. Our simulations seem to indicate that the first case is a more likely explanation, since, e.g., an uncorrelated noise at each spatial grid point, with the magnitude ∼ 10 −4 weaker than quantum vacuum noise, reproduces the g/N tot → 0 TWA limit.

Soliton core
Many studies [52,26,30,28,29] suggest a link between dissipative soliton dynamics and filling of a soliton core. To investigate this, we calculate the ratio ρ s /ρ b of the minimum density in the soliton core to the background density around the soliton. For the classical soliton solution these variables are related by (3). We stress the distinction between the quantum expectation value of the soliton depth, which we infer from the expectation value ρ s /ρ b found by locating the soliton and evaluating its depth in each stochastic realization, and the value of the density notch in the quantum expectation value ρ(x) , obtained by ensemble averaging over densities in many realizations. We first calculate the quantum expectation value ρ s /ρ b at time t = 6.7 × 10 −3 /ω -very soon after the phase imprinting when the density dip is still forming and all solitons in different stochastic realizations are still approximately centered at x = 0. We approximate ρ b by ρ(0.39l) in each trajectory. Figure 4(a) shows how the filling of the emerging notch increases due to quantum fluctuations at T = 0, as density fluctuations obscure the soliton. However, by time t = 6 × 10 −2 /ω after the formation of the density dip, the filling behaviour reverses and ρ s /ρ b becomes smaller in systems with large depletion. This increase in soliton depth due to quantum fluctuations is consistent with the aforementioned slowing down of the soliton velocities found in our previous study [21] in which case the increased phase uncertainty reduced the quantum expectation value of soliton speed | cos(φ/2)| due to the nonlinear dependence of the speed upon the phase jump. The classical formula between the phase jump and soliton core density (3) also exhibits a negative curvature and we expect strong quantum fluctuations to reduce the quantum expectation value of soliton core density, reflected by cos 2 (φ/2) . This deepening of solitons is surprising considering recent studies [52,26,30] which suggest that a soliton core in individual runs may develop density peaks due to groundstate depletion, and we would naively expect shallower, faster solitons. Instead we find that the depth of a soliton in individual runs is dominated by large phase fluctuations which give rise to the opposite effect, so that, after the formation, a soliton core is deeper in systems with large depleted fraction. We now consider the behaviour of the soliton depth during the evolution in a harmonic trap. Figure 4(b) shows the quantum expectation value ρ s /ρ b as the soliton passes the point x = −0.4l for the first and second time during its oscillation, which occur at times t ≃ 4.8 and t ≃ 8. In this case we may determine ρ s in every trajectory by the density at position x at the time when the soliton passes, and ρ b by averaging the density at the point x between consecutive passes of the soliton. We find that for individual trajectories, the soliton depth decreases for each pass; i.e., there is a slight filling of the soliton core with time. This is accompanied by a gradual speeding of the solitons due to dissipation (γ < 0). Quantum fluctuations enhance filling of the soliton core. After formation, the soliton in systems with strong quantum fluctuations, however, is initially typically deeper than in the case of weak fluctuations, and figure 4 shows that they remain deeper despite experiencing greater filling.
Numerically tracking the location of the phase kink in each individual realization, we have established that the expectation value for the depth of the soliton core is larger for larger values ofN nc /N tot . Solitons may drift around due to quantum fluctuations and in different runs they are generally not located in the same spatial position at the same time -soliton position uncertainties grow as a function of time due to both initial velocity uncertainties and interactions of solitons with sound waves [21]. In an ensemble averaged density over many single realizations, corresponding to the quantum expectation value for the atom density ρ(x, t) , soliton core density profiles appear broader and shallower due to different soliton locations in individual runs, as displayed in figure 5. The soliton position has a particularly large uncertainty in the presence of strong fluctuations and we find that, despite the fact that the expectation value of the soliton core is deeper in individual realizations with stronger quantum fluctuations, the density notches become shallower in ρ(x, t) when the fluctuations are enhanced. The flattening of ρ(x, t) due to the wandering solitons is similar to the previous results obtained using a Bogoliubov analysis [28,29] around a stationary dark soliton state showing that ρ(x, t) is smeared out due to fluctuations in the soliton positions.
Matrix product state simulation of a dark soliton in a bosonic atomic gas in a lattice in the tight-binding approximation with unit filling was studied in references [31,32]. From the expectation value of the atom density it was shown that an instantaneously imprinted phase kink with a vanishing density at the soliton core at the center of the lattice decays due to quantum fluctuations, as the soliton core gets filled with atoms. This is consistent with the aforementioned Bogoliubov analysis [28,29]. Moreover, based on the time-evolution of the atomic pair-correlation function it was argued that the soliton core also in single-shot experimental realizations is filled. In the simulations the initial state of g 2 (0, k) = a † k a † 0 a 0 a k = 0, since the central site j = 0 was empty (the location of the density dip of the soliton) and the other sites had approximately one atom. (Here a k denotes the annihilation operator for the atoms at the site k.) At later times g 2 became a flatter (and non-vanishing) function of k. Whether the filling of a soliton core in single-shot experiments in such a system duly happens, however, is inconclusive. This is because the evolution of g 2 can in many cases equally well be explained by a soliton whose core does not get filled in individual single-shot runs, but which randomly drifts along the lattice with the standard deviation of the soliton position increasing as a function of time. As a simple, idealized example, consider a lattice system where the long-range correlations may be ignored so that g 2 (0, k) ∼ n 0 n k . Here n k represents the atom density at the site k that is obtained by an ensemble average over many single-shot runs. The soliton is initially assumed to be located at the central site, with n 0 = 0 and n j = 1, for j = 0, and g 2 (0, k) = 0 for all k. At some later time Increasing position uncertainty at low atom numbers create a filling effect when the density is averaged over many realizations.
we may have n k = (L − 1)/L and g 2 (0, k) ∼ (L − 1) 2 /L 2 , for all k, where L denotes the number of sites. This can represent an outcome where every single-shot run yields a constant atom density (L − 1)/L for all k, but equally well a situation where each individual run has a single soliton at a random location along the lattice, so that the probability of finding a soliton (with a vanishing density) at an arbitrary site k in each realization is 1/L and the probability of not finding a soliton at the same site (L − 1)/L.
Our numerical simulations of quantum dynamics of a soliton with the corresponding classical value of the imprinted phase jump φ c = π in a lattice demonstrate a similar phenomenon (as shown in Section 6): Individual stochastic realizations exhibit drifting soliton dynamics along the lattice due to quantum fluctuations. The ensemble average of the atom density and the pair-correlation function g 2 become flatter as a function of time when the solitons have more time to drift over longer distances in the lattice. At the same time, however, the soliton cores in single-shot runs show very little, if any, effects of filling.

Colliding solitons
Recent studies of quantum dynamics of solitons in an optical lattice in the tight-binding limit with one mode function per lattice site have shown that soliton density may spread after colliding due to quantum fluctuations [31,32]. We examine soliton collisions in the absence of a lattice, and consider the relationship of such inelastic collisions with the soliton structure in the presence of strong quantum fluctuations.
We simulate the generation of two solitons using an optical imprinting potential (4) with two steps V φ = 4166.7hω at x b = −x a = 0.49l which in a classical system would imprint phase kinks of φ c = 2 and φ c = 2π − 2 with opposite velocities ± cos(φ c /2). We study the effects of quantum fluctuations on the imprinting of the phase kinks and the subsequent collisions of two such kinks (at about t ≃ 0.2/ω). As before, we do this by varying the effective interaction strength, corresponding to different values of the ground-state depletion, by keeping the nonlinearity N tot g fixed but adjusting the ratio g/N tot . Figure 6 shows |ψ W (x)| 2 for individual realizations of soliton collisions, and, for comparison, trajectories with the same stochastic realization, but using an imprinting potential that imprints only one of the solitons. We observe splitting/spreading of the solitons as they emerge from collisions in the presence of large quantum fluctuations.
Since the corresponding trajectories of single solitons do not exhibit this behaviour, we conclude that the splitting is induced by quantum fluctuations during the soliton collisions. Despite the spreading of the region of low atom density around the soliton core in individual realizations, the minimum density in the soliton core does not increase in these trajectories. We infer that one soliton splits into many solitons of comparable depth and velocity. Due to the dipolar motion in a harmonic trap, two solitons can experience repeated collisions. Recent experiment [43] showed seven repeated collisions of two solitons generated by merging two BECs. The parameters were close to those in our system: ω ≃ 333Hz, ω r ≃ 5592Hz and N tot ≃ 1700 (gN tot ≃ 205hωl); or ω ≃ 364Hz, ω r ≃ 2563Hz and N tot ≃ 950 (gN tot ≃ 50.1hωl). The experiment was performed in an elongated 1D trap, but without sufficiently strong transverse confinement that would have suppressed radial density fluctuations, and quantum fluctuations are not expected to be important. In our 1D system we find little quantum effects at similar atom numbers (other than a reduction in the soliton velocities). However, after reducing the atom number to N tot = 100, we observe soliton trajectories becoming perturbed [figure 7(a)] or disappearing after several collisions that may be detectable in an experiment using a tighter transverse confinement than the one in [43].
Increasing the number of solitons from two introduces the possibility of far richer dynamics. Systems of three or more harmonically trapped bright solitons may exhibit chaotic dynamics; the dynamics are not integrable due to the interplay between the harmonic motion and the soliton interaction [53]. The similarity between the interactions of bright solitons and those of dark solitons suggests that multiple dark solitons in a harmonic potential may also display chaotic dynamics. Multiple soliton collisions may also appear, e.g., in self-focusing and revival dynamics of BECs [54]. As a signature of quantum fluctuations, we find that in such systems, solitons are more susceptible to disappearing. We use the imprinting potential (4) to imprint six phase kinks all with different initial positions and velocities. Figure 7(d)-(f) shows |ψ W (x)| 2 for individual realizations. The enhanced dissipative effects visible in 7(d)-(e) where solitons trajectories become rapidly perturbed or disappear, may be related to dynamical instability of the system associated with chaotic dynamics, or merely due to the increased number of collisions experienced by each soliton in a short period of time.

Soliton in a combined harmonic trap and optical lattice
The classical nonlinear dynamics of a phase kink in a combined harmonic trap and optical lattice exhibits dynamical instabilities and the soliton may dissipate energy via sound emission even without quantum and thermal fluctuations [55,56]. The interplay between quantum fluctuations and the dynamical instabilities of the corresponding classical system was previously studied by us [21]. We found that the fluctuations can enhance the effect of instabilities, in some situations resulting in a surprising reduction in the position uncertainty of the solitons. Single-shot soliton trajectories displayed splitting and disappearing solitons in an optical lattice, similar to those following collisions, studied in section 5. Here we will focus on solitons that retain their integrity during the evolution. We consider very slow solitons that would have zero initial velocity in the classical case. Such a system exhibits strong fluctuation-induced effects and we calculate quantum statistical properties of the soliton dynamics, such as atom number fluctuations in an individual sites, waiting time distributions for the soliton to exit the initial site, atom populations and pair-correlations.
In the simulations of imprinting and dynamics we include in (5) the optical lattice potential V exp = mω 2 x 2 /2+sE r sin 2 (πx/d), where E r =h 2 π 2 /2md 2 is the lattice photon recoil energy and d is the lattice period. Here we set d = πl/4 and s = 1. After the simulations, in order to translate Wigner distributed stochastic variables to normally ordered expectation values, following the approach in [9,11] we define the ground-state operators a j for the individual lattice sites as where ψ j is the Gaussian ground-state wavefunction (Wannier function) of the jth well and ψ W (x, t) is the numerically integrated full Wigner wavefunction. The population of the jth site is thus a † j a j = a * j a j W − 1/2, and the number fluctuations in the jth site The pair-correlation function between the kth and the 0th site (the central site in which the soliton is imprinted) reads = a * 0 a * k a k a 0 W − We apply an imprinting potential that would, in the classical (GPE) case, prepare a phase kink of φ c = π (a soliton with zero velocity). This state is unstable [55], as small oscillations in a lattice site become amplified by sound emission so the soliton can escape from the central site and begin to drift along the lattice [56]. Quantum fluctuations seed  Figure 8(a) shows that the quantum uncertainty in the time τ s that it takes for a soliton to exit the central site becomes comparable to the quantum expectation value τ s .
The population expectation value for the central lattice site n 0 initially shows a soliton oscillating at the frequency of the lattice site, indicating overlapping oscillations between this site and the adjacent sites. This is followed by a rapid increase in central site occupation n 0 after time t ≃ 4/ω when there is a high probability that the soliton has left the site. This behaviour is reflected in the expectation value of the atom density in the lowest energy band at x = 0 [figure 8(f)]. We also find an initial oscillatory behaviour followed by a rapid increase in the pair correlation function [ figure 8(e)]. Previous studies of solitons in lattices with small atomic occupations [31,32] cited such an increase as evidence of filling of a soliton core in individual realizations, but here in large atom number systems we find a similar effect caused by solitons drifting along the lattice, while a soliton core in none of the individual realizations is significantly filled.
As well as affecting the expectation values n 0 and g 2 , the drifting behaviour of solitons has dramatic effects on the atom number fluctuations. Figure 9 shows the number squeezing parameter in the central lattice site δn 0 /n of a ground-state BEC in an optical lattice are sub-Poissonian (squeezed), such that δn 0 /n 1/2 0 < 1 [9,11]. For fast solitons, the soliton position uncertainty is small [figure 8 (b)], soliton cores between different realizations overlap, so there are predictable times when the soliton at high probability is not in the central lattice site. During these times number statistics in the site remain squeezed. The crossing of the soliton at x = 0 can be recognized as a sharp super-Poissonian peak in the dynamics of the squeezing parameter. For slower solitons, the position uncertainty of the soliton is so large that there are large atom number fluctuations in the central lattice site due to fluctuations in the crossing times of the solitons at x = 0.