Characterising arbitrary dark solitons in trapped one-dimensional Bose-Einstein condensates^(a)

Quantifying the nature of turbulence in a three-dimensional (3D) Bose-Einstein condensate requires a way of identifying individual phase defects present in the system. Quantum turbulence, characterised by a disordered system of entangled vortices and sound waves with a power law in the energy spectrum, proves to be an experimental challenge. Such an entangled 3D system provides a difficult system to optically image; if the individual vortices are not aligned with the direction of visualisation then they are barely visible within the turbulent cloud. We wish to provide a quantitative method of detecting vortices within 3D condensates. We hence take a step back and begin our journey with the study of one-dimensional (1D) repulsive Bose-Einstein Condensates (BECs) and the dark solitons (the 1D analogue for vortices in a repulsive condensate) within them. Solitons, localised waves characterised by the balance of nonlinearity and dispersion, are ubiquitous in nature: they are present in many nonlinear systems, from optics [1–3], thin films [4,5] and fluids [6–8] to atomic BECs [9–24]. Multiple methods of creating dark solitons [25–27] were proposed shortly after realisation of BECs [28–31] and the subsequent creation of low-dimension condensates [32]. These methods can be broken into two categories: phase imprinting [9,10,16,33] and density engineering [14,17,26], and also a combination of the two [23,27]. Recently, new methods to create solitons [34] through quenching have also been discussed [35,36]. Dark solitons are stable in BECs confined in quasi-1D geometries, at T = 0, and under certain specific forms of trapping potential, such as a harmonic trapping potential [37]. If any of these constraints are broken, the solitons are prone to decay, due to unstable excitation of the dark soliton into vortex rings/pairs (the snake instability) [11,38–43], thermal dissipation [44–46], and net sound emission [37,47–49], respectively. The first generation of dark soliton experiments in BECs showed the possibility of non-dispersive solitary waves propagating through the background condensate [9–11,28,33]. In these particular experiments, propagation was very short-lived. Dark solitons were shown to break down quickly due to a mixture of thermodynamic and dynamical instabilities which occurred if the background density was not strongly enough within the quasi-1D regime [38]. Unlike theoretical work based on the 1D mean-field limit, experimental quasi-1D condensates consist of cigar-like clouds of trapped atoms in which a strong radial confinement prevents excitations in the radial plane. The second generation of experiments was able to keep the condensate both cold and confined enough in the radial directions to sustain dark solitons for long periods (that is, for at least one oscillation of the soliton) [16,17]. Notable experiments include verification of the oscillation frequency of the soliton in a harmonic trap being $\omega/\sqrt{2}$ (where ω is the trapping frequency of the 1D system) for systems consisting of a single soliton [17,50] and the deviation from this prediction for multiple soliton systems [17]. Interesting theoretical models include the study of interactions of multiple dark solitons [51,52] and solitons in two-component systems [53–55]. Analysis of the complex structure present in a condensate —be that a single vortex/soliton or quantum turbulence [56–59]— requires accurate detection of the density. Techniques to image the density of the condensate include dispersive methods [60], absorption [61,62] and phase-contrast imaging [63,64]. Each technique has its own advantages and disadvantages; the choice is made depending on the type of experiment taking place. Imaging happens either in situ (that is, imaging the trapped condensate) or during a time-of-flight (TOF) expansion (releasing the condensate from the trap). Absorption methods, either in situ or after TOF expansion, are inherently destructive —due to the heating of the condensate by the imaging laser and the loss of the condensate in the case of TOF imaging. The viability [65] and search for minimally/non-destructive imaging techniques has been an active topic recently [66–68]. A particularly effective method of imaging the condensate with minimal destruction of the sample called Partial Transfer Absorption Imaging (PTAI) was experimentally demonstrated by Freilich et al. [69] and perfected by Ramanathan et al. [61] and has the distinct advantage of working for almost any optical depth. A small given percentage of atoms are outcoupled and imaged, leaving the original condensate almost unaltered. Using PTAI, the condensate can be imaged up to 50 times before breakdown [70]. The PTAI method was used successfully to visualise solitonic vortices [71] and reconnecting vortex lines [72,73] in cigar-shaped condensates. PTAI can be particularly useful for studying the evolving dynamics of moving solitons within a condensate. Shining light upon the outcoupled cloud and projecting this onto a charge-coupled device camera, experimentalists can thus gain multiple column integrated density profiles of the system [74,75].

We will employ the success of the PTAI method to show that, by taking multiple snapshots of a condensate, there is a link between the averaged density spectra and the dark solitons present. In this paper, we solve the Gross-Pitaevskii equation for a harmonically trapped 1D system with dark solitons present. We will study condensates with single and multiple solitons and introduce the idea of a spectral shift to identify the depth of any soliton present. We utilise the density engineering method in order to study more experimentally valid systems and, finally, will discuss the implications for 3D condensates.

In the zero temperature limit, the mean-field approximation for the wave function $\Psi(x,y,z,t)$ of a condensate provides a quantitative model of the dynamics in the form of the Gross-Pitaevskii equation (GPE),

$$\begin{equation} i\hbar\frac{\partial\Psi}{\partial t}= -\frac{\hbar^2}{2m}\nabla^2\Psi + V(x,y,z)\Psi + g|\Psi|^2\Psi - \mu\Psi,\quad \end{equation} \tag{ 1 }$$

under the normalization $\int|\Psi|^2 \textrm{d}^3r=\mathcal{N}$, where $\mathcal{N}$ is the number of atoms, ℏ is the reduced Planck constant $h/2\pi$, m the mass of the atomic species, μ the chemical potential, $g=4\pi\hbar^2a/m$ and a is the scattering length of the species. We take the trapping potential, V(x, y,z), to be harmonic, of the form $V=m[\omega^2_\perp(y^2+z^2)+\omega^2_xx^2]/2$, where $\omega_\perp$ and $\omega_x$ are the trap's parameters. We can reduce eq. (1) [76] by taking the perpendicular trapping frequencies to be sufficiently large, $\omega_{\perp} \gg \omega_x$, and integrate out the dependence on y, z. We also hence rescale the chemical potential $\mu_{1D}=\mu-\hbar\omega_\perp$ and $g_{1D}=g/2\pi\ell^2_{\perp}$ where $\ell_{\perp}$ is the harmonic oscillator length $\ell_{\perp}=\sqrt{\hbar/m\omega_{\perp}}$. We present our results using the natural units for a homogeneous condensate, which are time $\tau=\hbar/\mu$, length $\xi=\hbar/\sqrt{m\mu}$ and peak density $n_0=\mu/g$. There are two limits of elongated condensates [77]; when $an_0\gg1$ we enter the 3D cigar limit, and when $an_0\ll1$ we are in the 1D mean-field limit. In this paper, we work solely in the 1D mean-field limit. Written in terms of these natural units, the 1D GPE is

$$\begin{equation} i\frac{\partial\Psi'}{\partial t'}=-\frac{1}{2}\frac{\partial^2 \Psi'}{\partial x'^2}+\frac{1}{2}\omega'^2x'^2\Psi'+|\Psi'|^2\Psi'-\Psi',\quad \end{equation} \tag{ 2 }$$

where $t'=t/\tau$, $x'=x/\xi$ and $\omega'=\omega_x\tau$. All results presented are for $\omega'=0.02$. We normalise the condensate $\int \vert \Psi' \vert^2 \textrm{d}x'=\mathcal{N}'$ so that the peak density (the density at the trap minimum) $n_0'$ is unity. For a condensate with $\omega'=0.02$, this corresponds to $\mathcal{N}'\approx 94.2$. We choose $\omega'=0.02$ to ensure we are within the Thomas-Fermi limit $R_{x}'\gg\omega'^{(-1/2)}$. Here $R_x'=R_x/\xi$ and $\omega'^{(-1/2)}$ are the Thomas-Fermi radius and dimensionless harmonic oscillator length of the condensate, respectively. Provided this relation is satisfied, we can model the shape of the 1D condensate by the Thomas-Fermi profile [77,78],

$$\begin{equation} n_{TF}(x) = \left\{\begin{array}{l@{\quad}l} n'_0\left(1-\displaystyle\frac{{x}^{'2}}{R_x^{'2}}\right)\!, & \text{for } |x'| \leq R'_x,\\[8pt] 0, & \text{for } |x'| > R'_x. \end{array}\right. \end{equation} \tag{ 3 }$$

For $\omega'=0.02$, $R'_{x}\approx70$.

As discussed above, phase defects in a repulsive 1D condensate take the form of dark solitons. The analytic expression for a dark soliton of prescribed speed $v'$ in a homogeneous background is [18]

$$\begin{equation} \Psi'_S(x',t') = \sqrt{n'}\Big[B'\tanh{\Big\{B'(x'-x'_0(t'))\Big\}}+iv'\Big],\quad \end{equation} \tag{ 4 }$$

where $B'=\sqrt{1-v'^2}$, and $v'=v/\sqrt{n'}$. The parameter $x'_0(t')$ defines the location of the soliton at time $t'$, where $x'_0(t) = v't' + b$, where b is an arbitrary constant denoting the initial location of the soliton. To set up the initial condition for a numerical simulation of a single soliton in our condensate we multiply $\Psi'_s$ by the ground state $\Psi'$ of a harmonically trapped condensate. The density depletion at the minimum of the resulting condensate and the speed of the soliton are related via

$$\begin{equation} \Delta n' = 1-v'^2. \end{equation} \tag{ 5 }$$

For clarity, hereafter primes are dropped throughout. In the numerical simulations the GPE is solved via a fourth-order Runge-Kutta method using MATLAB with $\textrm{d}x=0.1$ and $\textrm{d}t=0.01$. At desired times, we compute the density spectrum, defined as the Fourier transform of the density of the condensate, $\tilde{n}(k) = \mathcal{F}(n(x))=\mathcal{F}(|\Psi(x)|^2)$, using inbuilt MATLAB subroutines for the Fast Fourier Transform, where k is the wave number $k=2\pi/x$.

We begin our investigation by first studying the density spectrum of the ground state of a 1D harmonically trapped condensate. Figure 1 compares the density spectrum of a ground-state solution, obtained numerically, to the density spectrum of its Thomas-Fermi approximation. The two spectra are consistent in the region of k-space corresponding to the central region within the Thomas-Fermi width, R_x , and the healing length, ξ, but they deviate at large wave numbers. This is expected: whereas the edges of the density profile of the ground state taper off smoothly, the abrupt cut-off of the Thomas-Fermi profile means that the spectrum has relatively large power at large k.

Zoom In Zoom Out Reset image size

As discussed above, a condensate containing a single soliton is easily obtained by multiplying the ground-state wave function by the expression for a dark soliton, eq. (4), in a homogeneous system (see fig. 2(a)). Upon comparing the spectrum of the ground state and the spectrum of the single-soliton condensate, we notice a shift rightwards towards the smaller lengthscale (larger k) region. We quantify this spectral shift and relate it to the soliton's depth $\Delta n$ (hence its speed v), and proceed in the following way. Consider a single soliton at the condensate centre. As mentioned earlier, a single soliton in a harmonically trapped condensate oscillates around the trap's minimum at frequency $\omega/\sqrt{2}$ with an amplitude depending on the soliton depth $\Delta n$. We let the system evolve in time for t = 500, averaging density spectra taken for every 0.5 time units. The total time of the simulation equates to a few oscillations of the soliton on the trap, a timescale achievable in experiments (using the trapping parameters of the experiment of Weller et al. [17], this corresponds to a timescale of 20–30 ms). We verify that the the resulting time-averaged spectrum does not change in time. Figure 2(b) shows that the addition of a soliton drastically shifts the density spectrum to larger wave numbers (compare the blue curve with the red and yellow curves). For the sake of making comparisons, the density spectra $\tilde{n}(k)$ we present are rescaled by $\tilde{n}_0$, defined as $\tilde{n}$ for $k\rightarrow0$. To quantify the shift of density spectrum to larger wave numbers arising from the presence of solitons, we define the relative spectral shift $\Delta k_i=k_i/k_{i}^{(0)}$; here k_i , which we refer to as the intercept wave number, is the wave number corresponding to the value $\tilde{n}_i = 10^{-5}\tilde{n}_0$ in the presence of the soliton, and $k_i^{(0)}$ is the intercept wave number of the ground state. For the ground state with $\omega=0.02$, $k_i^{(0)}=1.22$. Comparing $\Delta k_i$ obtained for a variety of soliton depths from 0.1 to 0.9 we observe the following power law relation between soliton depth and relative spectral shift of the density:

$$\begin{equation} \Delta k_i \sim \Delta n^\alpha, \end{equation} \tag{ 6 }$$

with $\alpha=0.55$ (see fig. 3). Clearly, the deeper the soliton the larger the spectral shift. We have verified that this result does not depend on the precise definition of n_i . We have also checked that eq. (6) is valid for a variety of harmonically trapped condensates, as long as these condensates are deeply in the Thomas-Fermi regime. Rescaling n₀ from 1, and hence altering the norm of the system, we have verified that the norm of the condensate does not affect the spectral shift. If we look at systems with $\omega\neq0.02$, we find that for $\omega\ll1$, the same relation (albeit with a different power law) is present. When $\omega\rightarrow1$, we see that it begins to falter; the condensate itself begins to have structure on the same lengthscales as the dark solitons.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The strong dependence of the spectral shift on the depth of a single soliton moving within the condensate enables us to determine the depth, and hence speed, of the soliton existing in the system. The next logical step is to assess how multiple solitons affect the density spectrum, whether a spectral shift is still observable, and finally whether it can be related to the number or the depths of the solitons.

For two solitons of equal depth, the power law relationship between the soliton's depth and the relative spectral shift (eq. (6)) holds true. The combination of two solitons of two different speeds, as presented in fig. 4(a), exhibits a more complicated spectral signature. Figure 4(b) shows results for two solitons. We immediately see that $\Delta k_i$ depends mainly on the deepest soliton (the shallower soliton having only a minor effect). Although direct determination of the depths (and hence speed) is not possible from the spectra, the concentric nature of the results presented in fig. 4(b) shows that if we know the shift, we can easily narrow down the depths to a range of results.

Zoom In Zoom Out Reset image size

We have established that in a two-soliton system the spectral shift is mostly affected by the deepest soliton. However, the methods of creating solitons, as discussed in the introduction, can often lead to a train of solitons. To make better contact with experiments, in this section we describe the spectral shift caused by a relatively large number, N, of solitons. We choose N = 9, with the first soliton of depth $\Delta n_1$ and the other eight of the same depth $\Delta n_{N-1}$, as shown by fig. 5(a). Figure 5 further supports the previous finding that the deepest soliton contributes the most to the shift of spectrum, hence a higher value of $\Delta k_i$. The spectral shifts shown in the figure display the same dependence on solitons' depth displayed in fig. 4: that is $\Delta k_i$ only begins to change when $\Delta n_{N-1}$ is comparable to $\Delta n_{1}$. More precisely, $\Delta k_i$ is noticeably affected only when $\Delta n_{N-1}$ is roughly half the value of $\Delta n_{1}$ (marked with vertical lines in fig. 5(b)).

Zoom In Zoom Out Reset image size

All results described in the previous sections refer to solitons imprinted into the ground state. In many experiments, because of the method used to generate them, solitons coexist with sound waves. The following density engineering method allows us to mimic this more realistic situation numerically. We apply a Gaussian potential of width σ and amplitude A for a time $\mathcal{T}$, before instantaneously removing it. This pulsed Gaussian potential has the form

$$\begin{equation} V(x,t) = \left\{\begin{array}{l@{\quad}l} \displaystyle\frac{1}{2}\omega^2x^2 + A\exp{({-x^2}/{\sigma^2})},& \text{for } t \leq \mathcal{T},\\[8pt] \displaystyle\frac{1}{2}\omega^2x^2, & \text{for } t > \mathcal{T}. \end{array}\right. \end{equation} \tag{ 7 }$$

This mimics the creation of solitons via density engineering methods and causes a density dip to form at the centre of the condensate. When removed, both solitons and sound waves are generated during the collapse inwards of the two sides, as shown in fig. 6(a). We measure the depth of the two main solitons created in order to imprint them into a ground state for comparison (fig. 6(b)). We choose $\sigma=2$, $\mathcal{T}=5$ and vary the amplitude, A, of the central Gaussian pulse. The larger the amplitude A, the deeper the leading solitons created (see fig. 7(a)). In fig. 7(b), we compare the relative spectral shifts $\Delta k_i$ resulting from solitons which are imprinted in the ground state (red line) and density engineered (blue line). We see immediately that although the spectral shifts $\Delta k_i$ in the two systems are not equal, the relationship between the soliton depths and the intercept wave numbers still holds, regardless of the presence of sound waves.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

When $A\rightarrow0$, no solitons are created. Notice the rapid change of the spectral shift for very small amplitudes A (fig. 7(b)); these small spectral shifts correspond to systems with no detectable soliton. The relation between the amplitude and the intercept wave number k_i begins to level out roughly at $A=0.2$, corresponding to two solitons of depth $\Delta n = 0.07$. We conclude that sound waves can shift the density spectrum, but the shift is small compared by the larger shift induced by solitons.

We have presented a method for accurately ascertaining the depth (and hence the speed) of a single soliton in a harmonically trapped condensate from the density spectrum alone. We have also shown that, in a system with multiple dark solitons, the spectral shift is mainly determined by the deepest soliton.

The analysis of the spectral shift which we have presented here for 1D systems may potentially be applied to 3D turbulent systems. A spectral shift of the momentum of the condensate has indeed already been reported in an experiment with a turbulent 3D condensate [79]. While 1D phase defects are solitons whose width depends on its speed, the phase defects present in 3D system are vortices. Since multiply charged vortices are unstable in most cases, the width of the vortex cores is constant (although, in a harmonically trapped condensate, this width increases near the edge). Therefore it may be possible to relate the measured spectral shift to the number or the length of vortices present in the system, thus providing a quantitative measure of the intensity of the turbulence in the condensate.

NGP, LG and CFB acknowledge support by the Engineering and Physical Sciences Research Council (Grant No. EP/R005192/1).