The decay and collisions of dark solitons in superfluid Fermi gases

We study soliton collisions and the decay of solitons into sound in superfluid Fermi gases across the Bose-Einstein condensate to Bardeen-Cooper-Schrieffer (BEC-BCS) crossover by performing numerical simulations of the time-dependent Bogoliubov-de Gennes equations. This decay process occurs when the solitons are accelerated to the bulk pair-breaking speed by an external potential. A similar decay process may occur when solitons are accelerated by an inelastic collision with another soliton. We find that soliton collisions become increasingly inelastic as we move from the BEC to BCS regimes, and the excess energy is converted into sound. We interpret this effect as being due to evolution of Andreev bound states localized within the soliton.


I. INTRODUCTION
Since the realization of Bose-Einstein condensation in dilute gases, experimental and theoretical research has shown that solitons play a key role in their dynamics (see Ref. [1] for recent reviews). Solitons have been generated in a wide variety of contexts such as, for example, phase-imprinting in elongated Bose-Einstein condensates (BECs) [2,3], transport in optical lattices [4,5], transport past single defects or through disordered potentials [6][7][8], BEC collisions [9][10][11] and interferometry [12,13]. Consequently, solitons might be expected to play an equally important role in the dynamics of degenerate Fermi gases. Recent work has predicted the existence of black [14] and grey solitons [15][16][17] across the Bose-Einstein condensate to Bardeen-Cooper-Schrieffer (BEC-BCS) crossover, and that they may perform stable oscillations in a trap [15]. However, if solitons in Fermi gases are to be as ubiquitous as solitons in Bose gases, they must be robust in the presence of defects, rapidly-varying potentials, superfluid flow, sound and other solitons. Can they be easily produced and are they long-lived? If not, what excitations will be observed?
In this paper, we address some of these questions by solving the time-dependent Bogoliubov-de Gennes (TDBdG) equations across the BEC-BCS crossover. We confirm a previous hypothesis that soliton solutions do not exist with a speed v above the bulk pair-breaking speed v pb [16]. This means that on the BCS side of unitarity, where v pb is smaller than the sound speed c, the maximum v is v pb , in accordance with the Landau criterion. In contrast, in the BEC regime the maximum v is c, as predicted by the Gross-Pitaevskii equation [18,19]. If a soliton in the BCS regime is accelerated to v pb by an external potential, it abruptly disappears and its energy is converted into sound. We stress that this decay process is distinct from the snake instability [20][21][22][23], by which a soliton can decay into quantized vortices in a three-dimensional superfluid. In our calculations, in fact, we assume that the dynamics are restricted to longitudinal motion only, the gas remaining uniform in the transverse directions. This is equivalent to assuming that, in a real experiment, the superfluid is sufficiently tightly-confined in the transverse directions to suppress the snake instability. Moreover, the predicted decay process occurs at zero temperature and is distinct from other decay mechanisms considered previously [24][25][26].
We analyse this decay process by looking at the low-lying eigenstates of the Bogoliubov quasi-particle spectrum. The lowest state is found to be localized in the vicinity of the soliton, with an energy below the gap for the continuous spectrum of the extended Bogoliubov states. The origin of this localized state (Andreev state [27][28][29]) is the fact that the energy cost for creating a fermionic excitation near a minimum of the order parameter is reduced with respect to the bulk value. The Andreev states associated with dark solitons at rest have already been discussed in Ref. [14]. Here we find that Andreev states also exist for solitons with finite v. However, as v approaches v pb , the excitation gap in the continuous spectrum approaches zero and the Andreev bound state is lost. We also study the soliton energy E s as a function of v (the energy dispersion), and find that the soliton still has a finite E s and phase jump across it at the maximum v. The energy dispersion shows interesting structure as v approaches v pb , in particular a local minimum. We explain that, in a real experiment, this local minimum would cause the soliton to decay at the value of v in the minimum rather than at the maximum v where the energy dispersion truncates. Moreover, the presence of a harmonic trap causes a further reduction in the value of v at which the soliton decays. This finite-size effect is reduced as we weaken the trap.
The existence of a maximum soliton speed, or minimum soliton energy, has important consequences for the physics of soliton collisions across the crossover. In the BEC limit, due to the integrability of the one-dimensional Gross-Pitaevskii equation for a homogeneous gas, soliton collisions are elastic and hence solitons always emerge from a collision with the same velocity and energy at which they collided [30,31]. However, we find that soliton collisions become increasingly inelastic as we reduce 1/k f a. As for the decay of a single soliton, fermionic quasi-particles are found to play a key role also in these inelastic collisions. This is consistent with the fact that the same collisions between solitons are found to be almost completely elastic [32] when described by means of a nonlinear Schrödinger equation, which is often used to investigate the dynamics of superfluid fermions in the BEC-BCS crossover, but including only bosonic degrees of freedom. Since, counter-intuitively, fast solitons have less energy than slow solitons [16,18,19], the inelastic collisions cause the solitons to accelerate. The energy lost by the solitons is dissipated as sound. We interpret this effect as being due to the evolution of the Andreev bound states as the soliton changes shape during the collision; this evolution eventually causes an energy transfer from the soliton to the continuum of bulk excitations. If the solitons lose so much energy that their energy after the collision would be less than the minimum soliton energy, they are destroyed by the collision.

II. METHODOLOGY
We consider a three-dimensional superfluid Fermi gas with equal populations of the two spin components. We model its dynamics across the BEC-BCS crossover by solving the TDBdG equations [15,33] , in which m is the atomic mass, U is the external potential and µ is the chemical potential. The order parameter is calculated as ∆ (r, t) = −g η u η v * η , in which g is given by [34]. Here a is the 3D s-wave scattering length characterizing the interaction between atoms of different spins, while E f = 2 k 2 f /2m and k f = 3π 2 n 1/3 are the Fermi energy and momentum of an ideal Fermi gas of density n, respectively. The cut-off energy E c is introduced in order to remove the ultraviolet divergences in the BdG equations with contact potentials. The density of the gas is n(r, t) = 2 η |v η (r, t)| 2 . We impose that the potential U has no y or z dependence, and hence we may write the functions u η (r, t) and v η (r, t) as u η (x, t)e i(ky y+kzz) and v η (x, t)e i(ky y+kzz) respectively, in which k y and k z are quantized according to k y = 2πα y /L ⊥ and k z = 2πα z /L ⊥ , where α y and α z are integers and L ⊥ is the width of the box in the y-and z-directions. As initial states at t = 0, we find stationary solutions of Eq. (1) [14,15]. We also search for solutions of Eq. (1) satisfying ∆ (x, t) = ∆ (x − vt), given that U = 0, by solving the equation where ξ = x − vt,Ĥ ξ = − 2 /2m ∂ 2 /∂ξ 2 − k 2 y − k 2 z − µ and ǫ η is the energy of level η. We use a generalized secant (Broyden's) method to find self-consistent solutions [16]. Hence we find travelling soliton solutions in a homogeneous gas which are stationary in the frame of the soliton, thus generalising to the BEC-BCS crossover the solutions for solitons propagating in a homogeneous Bose gas [18,19,35]. This technique enables us to determine the soliton properties more accurately in the homogeneous gas and identify effects due to the trapping potential by comparison with the time-dependent simulations.

III. SOLITON DECAY
Some of us have shown in a previous publication [16] that the energy dispersion of the soliton in the BCS regime truncates below the speed of sound c. This result suggested that, more generally, soliton solutions do not exist for speeds greater than v pb or c, whichever is the smaller. To support this hypothesis, and to show how the effect would manifest itself in a real experiment, we now present time-dependent simulations of the Bogoliubov-de Gennes equations [15,33] in which we accelerate the soliton above the critical velocity.
Firstly, for comparison, in Fig. 1 we present a simulation of a stable soliton oscillation in a trap at 1/k f a = −0.5. As in a previous publication [15], we consider a soliton oscillating in a 40 K superfluid, contained in the harmonic trapping potential U (x) = mω 2 x x 2 /2, with ω x = 2π × 50 rad s −1 , L ⊥ = 3.3 µm and a peak density n p = 1.  cloud near the center of the trap is shown, the low density tails of the cloud are outside of the field-of-view. The soliton begins at rest at a distance X 0 = 3.3 µm from the trap centre, and is consequently accelerated by the harmonic potential. As this happens, the density profile of the soliton becomes shallower and the phase jump across the soliton reduces from π. However, as the soliton passes the center of the trap and reaches its maximum speed, it is still localized in the density profile and the phase jump is well-defined. The soliton then begins to decelerate as it climbs the trap potential, and the phase jump increases towards π. If the simulation were allowed to continue the soliton would perform a complete oscillation in the harmonic potential.

Figures 1(c) and (d)
show the corresponding density profile and phase for the larger X 0 of 5.4 µm. The soliton begins to accelerate, and initially the soliton remains localized in the density profile with a clear phase jump. However, when the soliton reaches a critical speed it rapidly spreads out in the density profile and the phase jump disappears. At the end of the simulation we see only low-amplitude modulations of the density profile. By measuring their speed we identify these modulations as sound.
To gain a deeper physical understanding of this effect, we now study the soliton in the homogeneous gas using the BdG equations (2). The magnitude of the order parameter |∆(x)| is plotted in Fig. 2(a) for 1/k f a = −0.5 and v = 0 (solid blue curve), 0.2c (solid green curve with circles), 0.3c (dashed red curve) and 0.38c (dot-dashed black curve).
As v approaches the pair-breaking velocity v pb = 0.41c, the soliton becomes very shallow and broad. This is reflected in the profile of the lowest Andreev state, shown in Fig. 2 2 is tightly localised at the soliton position [14] (blue solid curve). Also note that there are many nodes in the function. As v increases to 0.2c (green solid curve with circles) and 0.3c (red dashed curve), the function widens and the nodes becomes gentle oscillations.
When v reaches 0.38c, the function is broad and smooth with almost no discernable oscillations (black dot-dashed curve).
As v approaches the pair-breaking velocity v pb = 0.41c for 1/k f a = −0.5, the excitation gap in the continuous spectrum of the extended Bogoliubov states, shown by the black solid line in the Fig. 2 inset, approaches zero. Below the continuous spectrum, there is the single Andreev bound state, which is plotted for various v in Fig. 2(b). As the excitation gap closes, the lower boundary of the continuous spectrum approaches the energy of the Andreev bound state, shown by the dashed line in the Fig. 2 inset. Note that the energy of the Andreev bound state is almost independent of velocity. Eventually, the localized Andreev bound state is lost and becomes an extended state.
The loss of the Andreev state can be understood by a simple analytic argument. We treat Eq. (2) semiclassically by considering the situation where the order parameter ∆ varies slowly over scales of the order of k −1 f [27]. This is true for the soliton near its maximum velocity, as shown in Fig. 2(a). This assumption allows us to decouple Eq. (2) into two separate second-order differential equations for the amplitudes u 0 = f (z) exp(−ik f z) and v 0 = g(z) exp(−ik f z), by neglecting terms proportional to ∂ 2 z f compared to terms in k f ∂ z f , and the same for g. Each of the resulting equations has two independent solutions of the form exp  calculations suggest that the soliton solution ceases to exist at a velocity slightly below v pb where the excitation gap is still finite.
In Fig. 3(a) we plot the energy of the soliton E s as a function of v (the energy dispersion) for various values of 1/k f a. In the BEC limit, the energy is given by the well-known Gross-Pitaevskii result [18,19,35] where α = 3/2. For 1/k f a = 1 (black dotted curve), our numerical data closely follows this result. At unitarity (red dashed curve), the energy initially decreases smoothly with increasing v, as for 1/k f a = 1, but then reaches a local minimum, subsequently briefly increases, before decreasing again rapidly. Finally the energy dispersion truncates at a small but finite value of E s [36]. In a previous publication, some of us predicted that at unitarity α = 2 in Eq.
(3) [16]. This is approximately true for small v. The presence of a local minimum in the energy dispersion is a crucial point, as we now explain. In general, the soliton energy E s is a function of µ and v 2 [15]. In Fig 1(c), the soliton begins at a distance X 0 from the trap center with v = 0. The soliton may move closer to the center of the trap, to larger µ, maintaining a constant E s , by increasing v 2 . However, once the soliton reaches the local minimum in the energy dispersion, it can no longer move to larger µ by changing v 2 . Hence, in a real experiment, the soliton must decay when v reaches the position of the local minimum, not when it reaches the point where the curve truncates.
Following this line of reasoning, we may reach another important conclusion. If the energy of the soliton were to decrease monotonically to zero with increasing |v|, as in the Gross-Pitaevskii equation, decay of the kind shown in Fig. 1(c) would be impossible. This is because the soliton could always maintain a constant E s , however much µ increased, by increasing v closer to the maximum v where E s goes to zero. Hence, the existence of the decay process shown in Fig. 1(c) implies a truncation in the energy dispersion, or at least a non-monotonic energy dispersion.
When we reduce 1/k f a from 0 to −0.2 [green solid curve in Fig. 3(a)], we again find a local minimum in the energy dispersion, but the curve truncates at a much larger value of E s . We now obtain α ≃ 2.2 in Eq. (3). The minimum value of E s is even larger for 1/k f a = −0.5 (blue dot-dashed curve), and α has increased to 2.8, but now the energy dispersion does not contain a minimum. Instead, we see a broad bump in the curve as the soliton approaches the  The non-monotonic slope of the energy dispersion in the BCS regime can be tentatively explained by considering the Friedel oscillations. These arise due to scattering of fermions at the Fermi energy by the soliton, and appear as ripples in the density or order parameter profile either side of the central minimum of the soliton [as shown in Fig. 1(a), for example]. For Friedel oscillations to appear, the atomic wave length at the Fermi energy, λ ≈ 2π/k f , should be larger than the size of the soliton. When the velocity of the soliton increases, its size grows and can exceed λ, causing the loss of the Friedel oscillations. For example, we find that the size of the bound state l b and λ are comparable at v = 0.3c for 1/k f a = −0.5. For v = 0, 0.2c and 0.3c, the Friedel oscillations are almost identical in the profile of the order parameter [blue solid curve, green solid curve with circles and red dashed curve in Fig. 2(a), respectively], although the depth of the soliton varies dramatically. In contrast, if v increases to 0.38c, the Friedel oscillations disappear and the order parameter acquires a smooth broad profile [black dot-dashed curve in Fig. 2(a)]. Due to the disappearance of the Friedel oscillations, we would expect different behaviour in the energy dispersion, as observed in Fig. 3(a). Notice that the existance of a small range of v where E s increases with v implies that, if a soliton were created by phase imprinting in this range of v, the snake instability would be absent [37]. Figure 3(b) shows the equivalent plot to Fig. 3(a) for the phase jump across the soliton ∆φ. In the BEC limit, the Gross-Pitaevskii equation predicts that ∆φ = −2 arccos (v/c). For 1/k f a = 1 (black dotted curve), our numerical data again agrees well with the Gross-Pitaevskii prediction. At unitarity (red dashed curve), ∆φ varies linearly with v for small v, as some of us predicted analytically in a previous publication [16]. However, as v approaches the maximum velocity we see a departure from this linear behaviour: in a similar manner to the energy dispersion, the ∆φ(v) curve forms a local minimum, then briefly rises before decreasing rapidly and terminating. This plot makes the important point that the soliton still has a finite phase jump when the energy dispersion terminates. This is expected because, as shown in Fig. 3(a), the energy dispersion truncates at finite energy. In the case of a soliton moving with a nonzero velocity v, a finite phase step implies finite energy by virtue of the kinetic energy associated with a superfluid current. When we decrease 1/k f a from 0 to −0.2 we find a similar variation of ∆φ with v (green solid curve), except that the curve now truncates at a larger ∆φ. As for the energy dispersion, at 1/k f a = −0.5 (blue dot-dashed curve) we find a broad bump instead of the local minimum.
Using the above results we plot the maximum soliton velocity v m in the homogeneous gas as a function of 1/k f a as red squares in Fig. 4. As explained earlier, in a real or numerical experiment, we predict that the measured v m is given by the position of the local minimum in the energy dispersion, not by the the point where the curve truncates, and hence we define v m to be the position of that minimum. energy dispersion, but we take v m to be the velocity at the "broad bump" for consistency. We compare this data to the bulk sound speed c (dashed curve) and bulk pair-breaking velocity v pb (solid curve) for a homogeneous gas, given by mv 2 pb = ∆ 2 + µ 2 − µ. For 1/k f a = 1, the maximum soliton velocity is given by c, as expected from the Gross-Pitaevskii equation [18]. The pair-breaking speed is much larger and hence not relevant because the pairs are strongly bound as molecules. However, for 1/k f a ≤ 0, the pair-breaking speed determines the Landau critical velocity, and hence limits v m . Thus, at unitarity we find a v m which is slightly smaller than v pb . As we decrease 1/k f a we find that v m closely follows the v pb curve. In fact, the data points lie slightly below v pb , and the deviation becomes larger as 1/k f a decreases.
On the same graph we also include a number of points (black crosses) showing the critical velocity predicted by the time-dependent simulations of a trapped superfluid. The points lie slightly below the predictions of the timeindependent calculations, the relative deviation being less at unitarity than at smaller values of 1/k f a. This small deviation is a finite-size effect, which probably occurs because the pair size for small 1/k f a becomes comparable to the size of the cloud. Hence it would be exceedingly difficult to probe soliton dynamics near the bulk pair-breaking speed in a real experiment. We include one point (green plus sign) in Fig. 4 for ω x = 2π × 25 rad s −1 and 1/k f a = −0.5 to illustrate that the critical velocity in the time-dependent simulations approaches the prediction of the time-independent calculations as the trap frequency is reduced.

IV. SOLITON COLLISIONS
We place two black solitons in a trapped Fermi superfluid, with a displacement of ±X 0 from the trap centre, and then evolve in time. The solitons are accelerated by the harmonic potential and collide at the trap centre. We may increase the speed of the collision by increasing X 0 . Figure 5 shows the evolution of the density profile with time for different X 0 and 1/k f a, with the other parameters as in Fig. 1. In the upper (middle) row X 0 = 5.1 (2.1) µm, and 1/k f a decreases from left to right (see caption).
Let us begin with the BEC regime 1/k f a = 1, shown in the left-hand column. For X 0 = 5.1 µm [ Fig. 5(a)], the two solitons pass through each other without changing their form, and subsequently slow down and come to rest at x ≈ ±X 0 , indicating that no energy has been lost by the solitons and hence the collision is elastic. For X 0 = 2.1 µm [ Fig. 5(f)], the solitons behave differently. In this case, the solitons slow down as they approach one another, come to a halt, and then move away without crossing. This is known as a black collision, whereas a crossing of two solitons, as shown in Fig. 5(a), is known as a grey collision. We again note that the two solitons return to their original positions after the collision, indicating that the collision is elastic. This behaviour of grey collisions for high impact speeds, and black collisions for low impact speeds, has already been observed in the GP equation [30]. When we reduce 1/k f a to 0.2, the behaviour is similar. We again observe a grey collision for X 0 = 5.1 µm [ Fig. 5(b)], and a black collision for X 0 = 2.1 µm [ Fig. 5(g)]. However, on close inspection we observe that, after the collision, the solitons come to rest at a value of |x| > X 0 , indicating that the solitons have a greater speed after the collision than before. This is particularly noticeable for the slower collision [ Fig. 5(g)]. Since, counter-intuitively, fast solitons have less energy than slow solitons (as shown in Fig 3), we conclude that the solitons have lost energy due to an inelastic collision. We also observe small ripples in the density which emanate from the collision, showing that the energy lost by the solitons is converted into sound.
We interpret this effect as being due to the evolution of fermionic quasiparticles localised in the solitons. In the bottom row of Fig. 5, we plot the evolution of the lowest Andreev state (specifically |v 0 (x, t)| 2 ) during the slow soliton collisions shown in the middle row of Fig. 5. As the solitons approach, the Andreev states see a double-well potential, and subsequently a deep single-well (if the collision is grey), causing them to oscillate and breathe in a non-adiabatic way. In the BEC regime this effect is negligible, because the Andreev states make a small contribution to the overall density ( u 0 dx dy dz ≫ v 0 dx dy dz, for example), and because the solitons repel each other at slow speeds. Consequently the slow collision causes no disruption to the lowest Andreev state for 1/k f a = 1.0 [ Fig. 5(k)].
However, for 1/k f a = 0.2 we observe oscillations in the amplitude of the lowest Andreev state [ Fig. 5(l)]. At one instant, shortly after the collision, |v 0 (x, t)| 2 reduces almost to zero. This means that particles are being transferred between different eigenstates of the Bogoliubov spectrum, with a coupling between the localized Andreev state and the states in the continuum. These transitions are associated with density and phase oscillations which eventually cause an emission of sound and loss of energy from the soliton. This effect is enhanced for low impact-collisions, because the fermionic quasiparticles have more time to move in response to the change in potential. The key role played by fermionic quasiparticles in this process is also confirmed by the fact that the same collisions are found to be elastic when a purely bosonic density functional theory (i.e., a nonlinear Schrödinger equation for the pairing field) is used in place of the TDBdG equations [32].
As we further reduce 1/k f a, the collisions become increasingly inelastic. At unitarity [Figs. 5(c) and (h)], the collisions create a great deal of sound and the solitons subsequently come to rest at much larger |x| than X 0 , particularly for X 0 = 2.1 µm [Fig. 5(h)]. The inelastic collision is associated with very pronounced oscillations in |v 0 (x, t)| 2 [ Fig. 5(m)]. We also observe that the collision for X 0 = 2.1 µm is no longer black, due to the increasing mass of the solitons [15]. When 1/k f a reaches −0.35, the solitons are destroyed by the collision for X 0 = 5.1 µm [ Fig. 5(d)], but survive if X 0 = 2.1 µm [ Fig. 5(i)]. Solitons can be destroyed by the inelastic collisions if their residual energy after the collision is less than the minimum soliton energy [38]. For 1/k f a = −0.5 and X 0 = 5.1 µm, the minimum energy (or maximum speed) is reached before the solitons even collide [ Fig. 5(e)]. For X 0 = 2.1 µm, the solitons collide but then are immediately destroyed [ Fig. 5(j)]. Consequently the lowest Andreev state is no longer localised after the collision [ Fig. 5(o)].
We quantify how inelastically the solitons collide in Fig. 4. Here we plot the soliton speed immediately before and after the soliton collisions presented in the middle row of Fig. 5 as pairs of downward and upward-pointing triangles respectively. For 1/k f a = 1, the collision is elastic, and hence the downward and upward-pointing triangles lie on top of each other. For 1/k f a = 0.2, the triangles are a small distance apart, indicating a slightly inelastic collision. This distance increases as 1/k f a is reduced to 0, showing that the collision is becoming increasingly inelastic. However, the speed of the soliton after the collision is still far below the maximum speed observed for a single soliton, shown by the black crosses. However, when 1/k f a is reduced to −0.35 the final soliton speed is very close to the critical value. This explains why collisions for larger X 0 or smaller 1/k f a destroy the solitons.

V. CONCLUSIONS
We have shown that in a two-component superfluid Fermi gas soliton solutions cease to exist as the soliton approaches the bulk pair-breaking velocity v pb on the BCS side of the resonance. At this point, the variation of soliton energy with velocity (the energy dispersion) truncates at a finite value of energy and phase jump across the soliton. Close to v pb , the energy dispersion shows interesting structure, and sometimes a local minimum, which can reduce the maximum soliton velocity as measured in a real or numerical experiment. The presence of a harmonic trap causes a further reduction in the maximum observed soliton velocity. If a soliton is accelerated to its maximum velocity by an external potential, it abruptly disappears and its energy is converted into sound. Solitons may also be accelerated by inelastic collisions with other solitons, and the energy lost by the solitons is emitted as sound. This can destroy the solitons if the energy lost as sound reduces the soliton energy below the truncation point in the energy dispersion. We find that soliton collisions become increasingly inelastic as we move from the BEC to BCS regime.
On the one hand, our results impose some limitations on experiments aimed to observe soliton oscillations and collisions in Fermi gases. Solitons on the BCS side of the resonance must be prepared close to the trap centre and accelerated gently in order to avoid their decay into sound. On the other hand, it is encouraging that solitons in fermionic superfluids behave differently to those in bosonic superfluids, and that this physics is not captured by bosonic models. The solitons in fermionic superfluids are sensitive to the fermionic degrees of freedom and so may be used as a tool to further characterize these gases.