Dark soliton collisions in superfluid Fermi gases

In this work dark soliton collisions in a one-dimensional superfluid Fermi gas are studied across the BEC-BCS crossover by means of a recently developed finite-temperature effective field theory [S. N. Klimin, J. Tempere, G. Lombardi, J. T. Devreese, Eur. Phys. J. B 88, 122 (2015)] . The evolution of two counter-propagating solitons is simulated numerically based on the theory's nonlinear equation of motion for the pair field. The resulting collisions are observed to introduce a spatial shift into the trajectories of the solitons. The magnitude of this shift is calculated and studied in different conditions of temperature and spin-imbalance. When moving away from the BEC-regime, the collisions are found to become inelastic, emitting the lost energy in the form of small-amplitude density oscillations. This inelasticity is quantified and its behavior analyzed and compared to the results of other works. The dispersion relation of the density oscillations is calculated and is demonstrated to show a good agreement with the spectrum of collective excitations of the superfluid.


I. INTRODUCTION
Solitons are among the most fascinating nonlinear phenomena in physics. Arising from an interplay between dispersive and nonlinear effects of the underlying medium, these solitary waves retain their shape while propagating at a constant velocity. Solitons emerge in a wide variety of physical systems, including optical fibers, classical fluids and plasmas. More recently, they have also become a subject of interest in the ultracold atoms community. Due to the fact that ultracold quantum gases are well-controlled and highly tunable nonlinear systems, they form an ideal environment for studying the properties and dynamics of solitons. In particular, solitons in ultracold atom clouds most often manifest themselves as dark solitons, which are characterized by a localized density dip and a jump in the phase profile of the order parameter.
Dark solitons have been theoretically and experimentally studied in Bose-Einstein condensates [1][2][3][4][5] and superfluid Fermi gases [6][7][8][9]. In both of these systems, they are subject to an instability mechanism called the snake instability [9][10][11][12], which makes the soliton decay into vortices if the radial width of the atom cloud is too large. Since this inhibits the creation and observation of stable solitary waves in three-dimensional (3D) quantum gases, the preferred set-ups to study dark solitons are elongated quasi-one-dimensional (1D) clouds.
Solitons are often portrayed as particle-like excitations. The main reason for this is that they preserve their identity not only while propagating, but also when interacting with each other: when two solitons collide, they re-emerge again as two solitonic waves [13,14]. In the case of Bose-Einstein condensates, collisions between dark solitons have been thoroughly studied both theoretically [15,16] and experimentally [17,18].
For fermionic superfluids on the other hand, this topic has been investigated far less extensively. Dark soliton collisions have been theoretically analyzed in quasi-1D Fermi superfluids by means of a generalized nonlinear Schrödinger equation [19] and through numerical simulations of the time-dependend Bogoliuobov-de Gennes (TDBdG) equations [20], but the methods that were used in these works imposed limitations on either the domain or the amount of the resulting data. The goal of the present paper is to further extend the study of dark soliton collisions in superfluid Fermi gases by using a recently developed effective field theory (EFT) that is capable of describing Fermi superfluids across the BEC-BCS crossover regime in a wide temperature domain [21]. This theory is based on the assumption that the order parameter changes slowly in both space and time, corresponding to the condition [22] that the pair field should vary over a spatial region larger than the pair correlation length [23,24]. The consequent limitations and validity of the EFT are discussed in Sec. II B. The theory has already been successfully employed in the description of both stable dark solitons and the snake instability mechanism in different regimes of temperature and population imbalance [25][26][27].
In this work, we use the EFT equation of motion that governs the dynamics of the order parameter to numerically simulate the collision of two solitons in a 1D Fermi superfluid and study the properties of the reemerging solitons across the BEC-BCS crossover. We demonstrate how the collision introduces a phase shift into the space-time trajectories of the solitons and analyze the effects of temperature and spin-imbalance on this quantity. We observe that the soliton interactions become inelastic when moving away from the BECregime, the lost energy being converted into small-amplitude density ripples that emanate from the point of collision. We determine the dispersion of these ripples and compare it to the dispersion of the collective excitations of the superfluid. Wherever possible, we qualitatively compare our results and their validity to those of the works mentioned above.
The (quasi-)1D setting in which we investigate the collisions ensures that the solitons are stable with respect to the snake instability [27] and is indeed the regime of interest for most of the current theoretical and experimental studies on this subject. The 1D propagating soliton solutions which constitute the initial states for the numerical simulations are obtained from analytic expressions for the phase and amplitude profile of the order parameter, which were derived in [25,26] for the case of a dark soliton on a uniform background. While traditionally ultracold gases are studied in set-ups with harmonic trapping potentials, the recent realization of box-like optical traps [28] provides an incentive to investigate uniform superfluids and the opportunity to experimentally test the predictions of the present work.
The remainder of the article is organized as follows: in Sec. II we give a brief overview of the theoretical model employed to examine ultracold fermionic systems. In Sec. III, we investigate and discuss the properties of dark soliton collisions in a 1D Fermi superfluid using numerical simulations. Finally, the conclusions of our study are presented in Sec. IV.

II. MODEL & METHOD
A. Effective field theory The system under consideration is an ultracold Fermi gas in which particles of opposite pseudo-spin interact via an s-wave contact potential. In the context of a recently developed finite-temperature effective field theory [21], this system can be described across the BEC-BCS crossover in terms of a superfluid order parameter Ψ(r, t) (representing the bosonic pair field), under the assumption that Ψ(r, t) varies slowly in both space and time. In the natural units ofh = 1, 2m = 1, E F = 1, the Euclidean-time action functional for the system is then given by where β is the inverse temperature, and the Hamiltonian density H is given by The thermodynamic potential Ω s reads while the coefficients C, D and E are defined as The functions f j (β, , ζ) in the above expressions are defined by with the fermionic Matsubara frequencies ω n = (2n + 1)π/β. In this treatment, the chemical potentials of the two pseudo-spin species µ ↑ and µ ↓ are combined into the average chemical potential µ = (µ ↑ + µ ↓ )/2 and the imbalance chemical potential ζ = (µ ↑ − µ ↓ )/2, the latter determining the difference between the number of particles in each spin-population. The quantity ξ k = k 2 2m − µ is the dispersion relation for a free fermion, E k = (ξ k + |Ψ| 2 ) 1/2 is the Bogoliubov excitation energy, and a s is the s-wave scattering length that determines the strength and the sign of the contact interaction. It is important to note that, while both the coefficients C and E in the action functional are kept constant and equal to the value they assume in the uniform system case, the coefficient D and the thermodynamic potential Ω s fully depend upon the order parameter [25]. The regularized real-time Lagrangian density that follows from (1) reads where |Ψ ∞ | is the value of the order parameter for a uniform system (i.e. the superfluid gap). Consequently, the subtraction of the term Ω s (|Ψ ∞ |) means that in the following calculations all energy values are considered as energy differences with respect to the energy of the uniform system. Since this background value |Ψ ∞ | and the average chemical potential µ essentially serve as input parameters for the coefficients (3)-(6), there is a certain freedom of choice for the values of these quantities which allows to tune the quantitative accuracy of the EFT. In this work, we assign to these parameters the mean-field values that are obtained by simultaneously solving the saddle-point gap and number equations [29]. At unitarity, this yields the values |Ψ ∞ | = 0.69 E F , µ = 0.59 E F and T c = 0.5 T F for respectively the gap, chemical potential and critical temperature of the system. These background values could be further improved upon by, for example, including fluctuations around the saddle point, using the values of quantum Monte-Carlo simulations [30][31][32] or even using the values derived from experimental measurements [33].
From the Lagrangian density (8), the equation of motion for the pair field Ψ is found to be: where the coefficientsD and A are defined as This equation is a type of non-linear Schrödinger equation which is closely related to both the Gross-Pitaevskii (GP) equation for Bose-Einstein condensates and the Ginzburg-Landau (GL) equation for Fermi superfluids [34]. The first term on the right-hand side can be identified as a kinetic energy term, while the non-linear term represents a system-inherent potential for the field. The ratioD/C can be interpreted as a renormalization factor for the mass of the fermion pairs [35] and the coefficient A determines the uniform background value of the system, since A(Ψ) Ψ = 0 is nothing but the aforementioned gap equation. It has been verified that in the deep BEC-limit (1/k F a S 1) this equation of motion correctly tends to the GP equation for bosons with a mass M = 2m and an s-wave boson-boson scattering length a B = 2 a s [36]. In the most general case, equation (9) describes a 3D fermionic superfluid. To limit ourselves to the study of (quasi-)1D systems with a uniform background, we will consider only solutions of the form Ψ(r, t) = Ψ(x, t), thus imposing the one-dimensionality of the system directly onto the pair field. The 1D equation of motion then becomes: The expressions for the coefficients C,D, E and A remain the same as in (3)- (6) and (10), butD and A now only depend upon one spatial parameter through their dependence on Ψ(x, t). For the case of a 1D dark soliton that propagates with a constant velocity v s on a uniform background, equation (11) can be solved analytically and an exact solution Ψ s (x − v s t) can be found [25,26]. Specifically, we write the complex order parameter as Ψ(x, t) = |Ψ ∞ | a(x, t) e iθ(x,t) and solve the resulting equations for the amplitude modulation a(x) and the phase field θ(x). This results in an expression for the phase profile and a relation for the position (i.e. the distance from the soliton center) in function of the amplitude where X(a) and Y (a) are defined as The relative bulk amplitude is given by a ∞ = 1, while the relative amplitude at the center of the soliton a 0 = a(x = 0) is determined as the solution of X(a 0 ) − v 2 s Y (a 0 ) = 0. This solution for the pair field exhibits the characteristic phase jump and amplitude dip at the position of the dark soliton, both of which are intrinsically connected to the soliton velocity v s : the lower the magnitude of v s , the larger the phase jump across the soliton position and the deeper the density dip. A stationary soliton, which is often referred to as a black soliton, possesses the highest possible phase jump |∆θ| = π and a maximal depth (zero pair density) at its center. Finite-velocity solitons, often called gray solitons, will become less and less deep as their velocity increases, up until a critical velocity v c (which coincides with the sound velocity c s of the superfluid) for which both the density dip and phase jump disappear completely and a soliton solution no longer exists. This relation between the soliton velocity, depth and phase jump is illustrated in Figure 1.
To simulate a dark soliton collision in a 1D superfluid in the context of the above model, an initial state Ψ 0 (x, t) consisting of two counter-propagating solitons is numerically evolved in time, using the EFT equation of motion (9). This initial state is constructed by combining two analytical solutions Ψ s (x − x 1 − v 1 t) and Ψ s (x − x 2 − v 2 t) for dark solitons with respective velocities v 1 and v 2 and respective initial positions x 1 and x 2 . While in general the superposition of two solutions of a nonlinear equation is not necessarily a solution of the equation as well, it has been demonstrated that, due to their localized character, the superposition principle holds asymptotically for solitary waves at large spatial separation [37]. Indeed, sufficiently far from the position of each soliton, the modulus and phase of the order parameter will assume constant values, so the two solutions can be connected to each other if the distance between x 1 and x 2 is sufficiently large. The subsequent numerical time evolution of the initial state is carried out by discretizing the space-time grid and applying a finite-difference fourth order Runge-Kutta (RK4) algorithm. A detailed explanation of this numerical procedure is given in Appendix A.

B. Validity of the model
There are in general two approximations that determine the validity and limitations of the effective field theory outlined above. The main assumption of the model is that the order parameter changes slowly in both space and time, which corresponds to the condition that the pair field should vary over a spatial region much larger than the pair correlation length (also referred to as the Pippard correlation length in the context of superconducting systems). A detailed study of the limitations imposed by this condition was carried out in [26]. By comparing the soliton width to the pair correlation length for the whole relevant region of the {(k F a s ) −1 , β, ζ} space, it was revealed that the reliability of the theory is not guaranteed on the BCS-side of the resonance at low temperatures. Hence, results found in this regime should be treated with additional caution. A second remark that must be made is that in the full version of the effective field theory, as it was introduced in [21], the action functional (1) contains additional terms that result in terms with second order time derivatives of the pair field in the equation of motion (9) (with a form similar to the terms for the Here, vc is the critical velocity above which a soliton solution no longer exists. The position is given in units of k −1 F . spatial derivatives). While it would in principle be better to use the full second order time derivative (SOTD) model, we found that numerical SOTD simulations of soliton-soliton collisions exhibit inherent instabilities in a large area of the parameter domain. These instabilities are most likely caused by the fast deformation of the condensate at the moment of the collision. We therefore omit the second order time derivatives and make use of the first order time derivative (FOTD) version of the theory. It has been verified that for low-velocity solitons these omitted terms are of less importance and the deviations between the SOTD and FOTD calculations are usually small [36].

III. RESULTS
In this section we present the results of the numerical simulations of dark soliton collisions in a 1D superfluid Fermi gas, based on the EFT model that was described in the previous section. The initial state of each collision process is constructed as described in Sec. II A. The initial positions of the counter-propagating solitons are always chosen as x 1 = −20 k −1 F and x 2 = 20 k −1 F , while the box size is chosen to be L = 300 k −1 F . The resolution of the spatial grid is taken to be ∆x = 0.02 k −1 F . Figure 10 in Appendix A indicates the corresponding maximal value the time step ∆t can have in each interaction regime for the simulation to remain stable.

A. General collision process
We will start by describing and characterizing the evolution of the different types of dark soliton collisions.
An important feature of these collisions is that the solitons preserve their identity, re-emerging as two solitary waves after the interaction. This is often interpreted as the fact that the two solitons simply propagate through each other. As an initial state, we consider two dark solitons with constant velocities v 1 and v 2 , propagating towards each other through a uniform superfluid. Depending on the magnitude of these initial velocities, two different types of collisions can be identified. The first type of collision, which occurs if the initial velocities of the counter-propagating solitons are small enough, is illustrated in Figure 2. The left panels of Figure 2a show the evolution of the pair density profile of the superfluid, while the right panels show the evolution of the phase profile. We observe that at some time before the soliton centers reach other, both their respective velocities and the density at their centers decrease until they become zero. Accordingly, the phase jump associated with each soliton grows until it reaches a value of −π, at which point it flips to a value of π and causes the soliton velocity to change sign. As a result, the solitons propagate back into the opposite direction without having fully come together. On the space-time diagram 2b that tracks the motion of the soliton centers, it appears as if the two solitons are simply reflected off each other in a particle-like collision. However, as can be observed from the evolution of the phase and density, the solitons actually exchange energy during the collision, adjusting their velocity, depth and phase jump accordingly. As a result, the soliton that was initially propagating to the right will still be found propagating to the right after the collision and vice versa, supporting the notion that the solitons move through one another. This type of collision is sometimes referred to as a black collision [15], because for each soliton there will be some point in time for which the density at its center vanishes completely. The second type of collision, which is sometimes   referred to as a gray collision [15], is illustrated in a similar way in Figure 3. If the initial velocities of the solitons are sufficiently large, they will never become zero and the density at the soliton centers will never fully vanish during the collision. Instead, the solitons come together to form a single composite structure, at which point their phase jumps cancel each other out. After a limited amount of time, the solitons separate again and continue propagating through the superfluid, conforming to the picture that the solitons move through one another. It is interesting to relate this collision process to the general fact that creating a single density dip in a superfluid with a uniform phase background results in the creation of one or more pairs of finite-velocity solitons, a method that can be used for the density-engineering of solitons in experiments [16].
In  In the more general case of asymmetrical collisions, the studied quantities depend on both of the individual soliton velocities instead of only the relative soliton velocity, but their general behavior still agrees with the results discussed below.  showing an initial rapid decrease and then leveling out slowly. To explain this behavior, we also study the quantity ∆x min , which represents the smallest distance that occurs between the soliton centers during each collision. For a black collision, in which the solitons turn around before their centers reach each other, this value will be a finite number, while for a gray collision, in which the solitons merge completely, this quantity will be zero by definition. As can be seen in the inset of the figure, the behavior of ∆x min in function of both v r and (k F a s ) −1 seems to be related to that of the curves in the main graph. Slow moving solitons with a large depth will reach their turning point while they are still at a relatively large distance from each other, and subsequently acquire a larger shift with respect to their free trajectory. When the solitons come closer to merging completely however, this effect diminishes and the magnitude of the shift slowly levels out. The inset also indicates that the transition from a black to a gray collision consistently occurs when the relative velocity of the solitons is equal to the critical soliton velocity v c . This was already demonstrated for the case of Bose-Einstein condensates [38], but is now observed to be valid across the whole BEC-BCS crossover. Figure 5 shows the magnitude of ∆s for a symmetrical collision with relative velocity v r = 0.3 v F across the BEC-BCS crossover for varying degrees of population-imbalance. The spin-imbalance parameter η is defined as η = ζ/|Ψ ∞ | where |Ψ ∞ | represents the superfluid gap of the homogeneous unpolarized gas [39]. For the regular spin-balanced case (η = 0), ∆s decreases when moving away from the deep BEC-regime, reaches a minimum around (k F a s ) −1 = 0.5 and then increases again towards the BCS-regime. The increase at the BCS-side of the resonance could be expected: the pair density at the center of the soliton is much lower in this interaction regime, causing the solitons to turn around sooner and thus acquiring a bigger shift. Additionally, the non-monotonous behavior of the shift around (k F a s ) −1 = 0.5 can be explained by noticing that it is very similar to the behavior of the soliton width across the BEC-BCS crossover [26]. This implies that solitons with a bigger width start feeling each other's influence sooner and as a result acquire a greater spatial shift in the collision. This idea also persists when we consider the effects of spin-imbalance on the value of the shift. It has been demonstrated that dark solitons broaden as the system becomes more imbalanced [26]. In accordance with the above observation, this seems to result in an increase of the magnitude of the spatial shift with η. Using spin-imbalanced systems may therefore be a convenient way to make the spatial shifts of the solitons more clearly observable in experiments. The inset of the figure shows the effect of temperature on the shift. Since increasing the temperature slightly increases the soliton width, this also results in a small increase of ∆s.
Our predictions for the behavior of ∆s across the BEC-BCS crossover in a spin-balanced Fermi superfluid can be qualitatively compared to those of Ref. [19], where the authors also use a superfluid order parameter equation to calculate the relative shift of a collision between two high-speed solitons across the BEC-BCS crossover in a quasi-1D system. Their result displays a local minimum around (k F a s ) −1 = 0.5, a local maximum close to unitarity and a subsequent decrease towards the deeper BCS-regime. While our own results also display a minimum around (k F a s ) −1 = 0.5, we don't find a maximum and subsequent decrease of this quantity at the BCS-side of the interaction domain. In particular, since our observations suggest that the magnitude of the spatial shift depends on both the width and depth of the colliding solitons, and since neither of these quantities exhibits non-monotonous behavior in the BCS-regime, we do not find a reason for the spatial shift to display any non-monotonous behavior in this regime either. On the other hand, the approximations and calculations in Ref. [19] focus mainly on solitons with a high velocity, which happens to be a regime in which our own results might exhibit deviations on the BCS-side, as was discussed in Sec. II B. Consequentially, a definite conclusion can not be drawn from the available data. Blue colors indicate a lower relative density, yellow colors a higher relative density. One can clearly observe how density waves emanate from the point of collision, carrying away a small amount of energy from the colliding solitons. Our numerical simulations demonstrate that dark soliton collisions are not elastic across the whole BEC-BCS crossover. In the case of a (quasi-)1D BEC, soliton interactions are known to be completely elastic [15,17], and indeed we find that no inelasticity occurs in the deep BEC-limit of the interaction domain  which identifies the largest density wave that is produced in each of these collisions. One can clearly see the correlation between the behavior of these quantities: the larger the generated density waves are, the more energy the solitons lose in the collision, and the higher the resulting velocity increase. In general, the inelasticity is small for slow-moving solitons, then increases to reach a maximum value at a fixed value of v r /v c , and finally becomes very small again for high-velocity solitons. Considering how the critical soliton velocity and thus the velocity range become larger when moving towards the BCS-regime, the inelasticity will peak at increasingly higher velocities. Figure 8 shows 1.8 Our results for the inelasticity of soliton collisions across the BEC-BCS crossover can be compared to those of Ref. [19] and [20]. In the former work, in which the authors also make use of a nonlinear macroscopic order parameter equation, soliton collisions are observed to be as good as elastic across the whole BEC-BCS crossover. In the latter work on the other hand, in which a harmonically trapped Fermi superfluid is studied through numerical simulations of the TDBdG equations, inelastic soliton collisions are clearly observed. Due to the fact that the TDBdG method is computationally very demanding, only four values of ∆v were calculated across the BEC-BCS crossover. While this prevents a rigorous comparison between the different formalisms, we can nevertheless observe that the velocity increases found with the TDBdG method are almost one order of magnitude larger than those found in the present work. The authors determine that the inelasticity is caused by fermionic quasiparticles that are localized in the solitons (i.e. Andreev bound states) and identify this as the reason that in Ref. [19] only elastic collisions were observed, since these type of localized quasiparticle states are not present in a bosonic effective field theory. In the context of the currently applied EFT, the existence of a continuum of fermionic single-particle excitations with energy E k does enter the formalism in an implicit way through the background theory, but no dynamical pair-breaking processes are contained within the Lagrangian for the bosonic pair field. This means that, similar to Ref. [19], localized quasiparticles and Andreev bound states are absent in the current description. However, the present results do reveal a small but observable inelasticity for the collisions. This suggests that the reason that no inelastic collisions were observed in Ref. [19] might actually be due to the fact that the work focuses on solitons with a very high velocity, for which, according to the results in Figure 7, the energy loss will become close to negligible. On the other hand, the fact that the observed energy losses are much larger in the TDBdG formalism than in the present case indicates that, while the inelasticity might not be completely due to the localized fermionic quasiparticles, they do provide a large contribution to this process.

D. Dispersion of the density oscillations
As can be observed in the space-time plot of the pair density profile in Figure 6, the density ripples that emerge from an inelastic soliton collision generate a chirp-like wave pattern in which the frequency of the oscillations increases along the length of the wave train. In order to obtain the corresponding dispersion relation ω(q), we determine the wavelength λ i (measured as the distance between two consecutive peaks) and the velocity v i of each matter wave and apply the relations q i = 2π/λ i and ω i = q i v i . The resulting dispersion of the oscillations can then be compared to the spectrum ω s (q) of the collective excitations of the superfluid, which in the context of the present EFT is determined by calculating the poles of the inverse propagator for the bosonic fluctuations [21]. Up to fourth order in q, this leads to the form where the sound velocity c s and the coefficient λ depend on the system parameters. While the EFT values for c s at zero temperature show a very good agreement with the results of Gaussian pair fluctuation calculations across the whole BEC-BCS crossover [21], the predictions for the coefficient λ show deviations from those found by Kurkjian et al. [40] when moving towards the BCS-regime [36]. It is also important to recall that, in the context of the current work, we have omitted from the EFT equation of motion all second order time derivatives of the pair field, limiting ourselves to only first order temporal variations. As was noted in Sec. II B, this introduces deviations in particular for high-velocity excitations and on the BCSside of the resonance. In Figure 9, the dispersion of the density ripples is compared to the spectrum of at low q. This agreement between the excitation spectrum and density wave dispersion is also observed at higher values of the temperature.

IV. CONCLUSIONS
In this paper we studied the collisions of dark solitons in 1D superfluid Fermi gases by means of a finitetemperature effective field theory already employed to describe the properties of single dark solitons. Using numerical simulations based on the EFT equation of motion, we demonstrated that the collisions introduce a spatial shift into the soliton trajectories and we studied the effects of spin-imbalance and temperature on its magnitude. The fact that the presence of spin-imbalance increases the magnitude of the spatial shift could provide a convenient way to make this quantity easier to observe in experiments. When moving away from the deep BEC-regime, soliton collisions were found to become inelastic, resulting in an increase of the solitons' velocities after the collision. The inelasticity peaks in the middle of the soliton velocity range, but is close to negligible for very slow and very fast moving solitons. This could explain the absence of inelastic collisions in Ref. [19], where only high-velocity solitons were studied. The fact that the observed changes in velocity are much smaller than in Ref. [20] corroborates the fact that the main contribution to the inelasticity in Ref. [20] is due to the presence of localized fermionic quasi-particles, which are absent in the current theory. The energy that is lost in the inelastic collisions was observed to be converted into trains of small-amplitude density oscillations, whose dispersion was demonstrated to show a good agreement with the spectrum of collective excitations of the superfluid. The numerical simulations carried out in this work were based on analytical solutions for 1D stable solitons in a uniform superfluid, as derived in [25,26]. Even though ultracold gases are traditionally studied in set-ups with harmonic trapping potentials, the recent realization of box-like optical traps [28] can provide the opportunity to test the predictions of this work in experiment. While the effects of a single soliton collision is most likely too small to be reliably measured, the (necessarily) non-uniform boundaries of the trapping potential can be used to make the solitons turn around and undergo multiple repeated collisions in the uniform region. Subsequently, the accumulated effect of these consecutive collisions might be observable in experiments.
in (A.2) by central finite difference formulas. Using the notations Ψ m,n = Ψ(x m , t n ), w m,n = |Ψ(x m , t n )| 2 , D m,n =D[Ψ(x m , t n )] and A m,n = A[Ψ(x m , t n )], the discretized form of (A.2) is given by f (Ψ m,n ) = i C 2mD m,n Ψ m+1,n − 2 Ψ m,n + Ψ m−1,n ∆x 2 − i Ψ m,ñ D m,n A m,n + E m w m+1,n − 2 w m,n + w m−1,n ∆x 2 (A.3) Since far from the solitons we expect the superfluid to assume its uniform bulk value, we require the field to remain constant and its derivatives to be zero at the boundary points of the chosen spatial grid. If now, for a certain time step t n , we know the values Ψ m,n for all positions x m , the explicit RK4 method allows us to calculate for every position the value Ψ m,n+1 of the next time step by using the following algorithm [41]: For the finite-difference algorithm to be stable, the ratio ∆x/∆t must be sufficiently large, meaning that the smaller we choose ∆x, the smaller the value of ∆t has to be to prevent the simulation from diverging. This results in a trade-off between the preferred resolution of the spatial grid and the amount of time steps that will be needed to reach the desired point in time. Figure 10 shows the minimal value which ∆x/∆t must have across the interaction domain in order for the algorithm to be stable. One can observe that the condition for stability becomes stricter when moving from the BEC-regime to the BCS-regime. We find that for the grid step ∆x = 0.02 k −1 F that is used in this work, the maximal stable values of the time step ∆t are sufficiently small to ensure that the convergence error of the finite-difference algorithm is much smaller than the error coming from the finite spatial step size. Therefore, the numerical errors on all calculated quantities are determined by the spatial resolution of our grid.