Position-dependent fusion and non-fusion of colliding dark solitons in harmonically trapped Bose–Einstein condensates

Abstract

Dark-soliton collisions exhibit two distinct regimes, fusion and non-fusion, characterized by different density profiles at the collision instant. For symmetric soliton pairs on a homogeneous background, the two regimes are separated by the critical velocity v_c=c_s/2. In a harmonically trapped Bose–Einstein condensate, however, both the background density and the soliton velocities vary along the trajectories, making the incoming local conditions position dependent. Combining the exact homogeneous two-soliton solution with a local-density description of the trapped background, we obtain {\rho }_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}=n(x_c^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}})[2A_{\mathrm{a}\mathrm{v}\mathrm{g},\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}(x_c^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}})-1{]}^2, where x_c^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}} is the encounter position predicted from the isolated-soliton trajectories and A_{\mathrm{a}\mathrm{v}\mathrm{g},\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}=|v_1-v_2|/(2c_s) is the normalized local relative speed evaluated at that position. The expression predicts the minimum density of the merged dip in the fusion regime and the height of the central bump in the non-fusion regime, with the local threshold at A_{\mathrm{a}\mathrm{v}\mathrm{g},\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}=1/2. Gross–Pitaevskii simulations of symmetric and asymmetric soliton pairs in both head-on and catch-up geometries provide a direct determination of the collision outcome and of the collision-center density {\rho }_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}}. After normalization by n(x_c^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}), the measured densities collapse onto the curve (2A_{\mathrm{a}\mathrm{v}\mathrm{g},\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}-1{)}^2 over both collision regimes.

Graphical abstract

Keywords

dark solitons / Bose–Einstein condensates / soliton collisions / Gross–Pitaevskii equation / local-density approximation

1 Introduction

Dark solitons are localized density depletions on a uniform background, accompanied by a characteristic phase jump across the notch; they represent prototypical nonlinear excitations of the defocusing nonlinear Schrödinger equation [1–3]. They were first explored experimentally in nonlinear optical systems [4, 5] and have since become a prominent class of excitations in repulsively interacting Bose–Einstein condensates. In this context, dark solitons have been generated and tracked using phase imprinting [6, 7], density and phase engineering [8], and the release of a condensate from a double-well potential into a harmonic trap [9]. Recent theoretical proposals have further explored impurity-assisted control schemes for manipulating dark-soliton motion and controlled collisions [10]. Related studies have also addressed interaction-quench-induced nonlinear excitations and ring-dark-soliton dynamics in multicomponent condensates [11–14]. Dark-soliton families in nonlinear and linear lattices have likewise been reported [15, 16].

In the limit of weak trapping and slowly varying background density, a dark soliton behaves as an effective quasiparticle whose position follows an oscillatory trajectory in a harmonic trap [17–20]. This trajectory-level description has been used to capture single-soliton oscillations, adiabatic motion in inhomogeneous backgrounds, and interaction-induced phase or position shifts [3, 21–24]. In the present work, we do not aim to introduce a new single-soliton trajectory theory. Instead, we use the local trajectory information as a kinematic input for a different question: what density structure is formed at the instant when two dark solitons collide in a spatially inhomogeneous condensate?

Dark-soliton collisions are also well understood for the homogeneous, integrable nonlinear Schrödinger equation. In that setting, solitons pass through each other elastically, acquiring only phase or position shifts, and the exact two-soliton solution gives the full density profile during the collision [2, 3, 25, 26]. Beyond the integrable limit, collision and fusion have also been studied in cubic–quintic nonlinear Schrödinger systems with external potentials [27]. For symmetric dark-soliton pairs on a uniform background, this solution yields two distinct density profiles at the collision instant. When the normalized soliton speed is sufficiently large, the two density notches merge into a single composite dip, which we refer to as fusion. When the speed is smaller, the two notches remain separated by a central density bump, which we refer to as non-fusion. The corresponding symmetric homogeneous threshold is commonly expressed as the critical velocity v_c=c_s/2, where c_s is the sound speed [3, 22, 28].

Although this homogeneous velocity-threshold picture is useful, it is incomplete for trapped condensates. In a harmonic trap, both the background density and the soliton velocities vary along the trajectories. As a result, the collision is not characterized solely by the initial soliton parameters or by an individual soliton velocity. The encounter position selects the local background density, sound speed, and incoming relative velocity that enter the collision criterion. Therefore, a pair of solitons generated with fixed imprinting phases can sample different local dynamical conditions depending on where they meet in the trap.

A complementary way to state this problem is to distinguish between trajectory-level and density-level descriptions. Effective-particle theories describe where the solitons move and how their trajectories are shifted by the trap or by mutual interaction. However, they do not directly determine the density structure formed at the instant of overlap. Conversely, the exact homogeneous two-soliton solution gives the collision-core density, but only for a spatially uniform background. A local criterion is therefore needed to connect these two descriptions: the local motion fixes the relative velocity at impact, while the local background fixes the density scale on which the collision core is formed.

In this work, we formulate such a density-resolved local criterion for two-dark-soliton collisions in a harmonically trapped Bose–Einstein condensate. We take the measured collision-center density {\rho }_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}} as the central observable and compare it with the forward prediction {\rho }_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}} constructed from the incoming local conditions. Starting from the exact homogeneous two-soliton solution, we obtain a closed-form expression for the collision-center density and extend it to asymmetric collisions in terms of the normalized relative-speed parameter. We then generalize the expression to a trapped condensate by evaluating the density, sound speed, and incoming soliton velocities at the encounter position predicted by the isolated-soliton trajectories, while the collision outcome and {\rho }_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}} are obtained from Gross–Pitaevskii simulations. The resulting criterion contains the homogeneous critical velocity v_c=c_s/2 as a special symmetric limit, but it also provides a quantitative density observable across the fusion, non-fusion, and critical regimes.

The main contribution of the present work is not the use of the local-density approximation itself, but the identification of a measurable local criterion for the collision core. This criterion turns the conventional velocity-threshold picture into a density-resolved relation: in the fusion regime it gives the minimum density of the composite dip, whereas in the non-fusion regime it gives the height of the central bump between two separated notches. It also predicts a position-selected collision outcome that has no counterpart in a homogeneous system: the same pair of solitons can undergo fusion or non-fusion at different encounter positions in the trapped condensate.

We validate this description by Gross–Pitaevskii simulations on both uniform and trapped backgrounds. The simulations include symmetric and asymmetric soliton pairs, on- and off-center collisions, and both head-on and catch-up geometries. In particular, we show that trapped collision data collapse onto a single curve after normalization by the local background density, demonstrating that the normalized collision-center density is controlled by the local relative-speed parameter rather than by the specific preparation protocol. This collapse establishes the proposed expression as a local organizing principle for the density structure formed at impact.

The remainder of this paper is organized as follows. Section 2 establishes the analytical mapping from imprinting parameters to soliton trajectories, which provides the local velocities and background densities required for the collision criterion. Section 3 derives and validates the collision-center density formula on a homogeneous background. Section 4 extends this description to a harmonic trap, presents numerical verifications, and examines the position dependence of the collision outcome. Section 5 summarizes the findings and outlines possible extensions.

2 Soliton dynamics and parameter mapping in a harmonic trap

We consider a quasi-one-dimensional Bose–Einstein condensate in a highly elongated harmonic trap with {\omega }_x\ll {\omega }_{\perp } [29, 30], described within the mean-field approximation by the dimensionless Gross–Pitaevskii equation [31–33]

(1)\mathrm{i}{\displaystyle \frac{\mathrm{\partial }\psi }{\mathrm{\partial }t}}=-{\displaystyle \frac{1}{2}}{\displaystyle \frac{{\mathrm{\partial }}^2\psi }{\mathrm{\partial }{x}^2}}+{\displaystyle \frac{1}{2}}{\mathrm{\Omega }}^2{x}^2\psi +|\psi {|}^2\psi ,

where \mathrm{\Omega }={\omega }_x/{\omega }_{\perp }\ll 1. Length and time are measured in units of the transverse oscillator length a_{\perp }=\sqrt{\hslash /(m{\omega }_{\perp })} and the inverse transverse frequency {\omega }_{\perp }^{-1}, respectively, and the wave function is scaled so that the nonlinear coefficient is unity. We choose the dimensionless chemical potential \mu =1, yielding a ground-state peak density n(0)\simeq 1.

The ground state is obtained numerically by imaginary-time propagation [34–36]; within the Thomas–Fermi approximation it follows the profile n(x)=1-{\omega }_0^2{x}^2 for , where R_{\mathrm{T}\mathrm{F}} is the Thomas–Fermi radius defining the cloud boundary and {\omega }_0=\mathrm{\Omega }/\sqrt{2} is the soliton oscillation frequency [17]. The numerical ground state is used as the background for soliton imprinting and for extracting local densities; the Thomas–Fermi expression is used in analytical predictions.

This section summarizes the map from imprinting parameters (x_0,\theta ) to the soliton velocity and depth at any point along its trajectory, providing the analytical input for the collision analysis of Sections 3 and 4. The map combines the homogeneous soliton relations, applied at the imprinting point through a local-density approximation, with energy conservation along the subsequent trajectory.

On a uniform background of density n_0, a dark soliton generated by a phase step \theta has velocity v=\sqrt{n_0}\cos (\theta /2) and depth D=\sqrt{n_0}\sin (\theta /2), related by v^2+D^2=n_0 [2, 3, 37]. In a trapped condensate, the imprinting process acts on a spatial scale much smaller than the length scale over which the background density varies. Within this local-density approximation [18, 19], the homogeneous relations apply at the imprinting point with n_0\to n(x_0):

(2)v_0=s\sqrt{n(x_0)}\cos (\theta /2),s=\pm 1,

(3)D=\sqrt{n(x_0)}\sin (\theta /2).

Here, s=\pm 1 specifies the initial propagation direction of the soliton, while the depth D is set by the local background at the imprinting point and is conserved along the subsequent trajectory [17, 18]; the velocity, by contrast, varies along the trajectory, with v_0 denoting its value at the imprinting point.

The soliton subsequently behaves as a quasiparticle oscillating in the harmonic potential with frequency {\omega }_0, following x(t)=A_{\mathrm{o}\mathrm{s}\mathrm{c}}\sin ({\omega }_{0}t+\phi ) [17, 20, 23, 24]. The velocity at any position x along the trajectory is determined by energy conservation,

(4)v^2(x)=v_0^2+{\omega }_0^2(x_0^2-x^2),

and the oscillation amplitude A_{\mathrm{o}\mathrm{s}\mathrm{c}} follows from the turning-point condition v(A_{\mathrm{o}\mathrm{s}\mathrm{c}})=0,

(5)A_{\mathrm{o}\mathrm{s}\mathrm{c}}=\sqrt{x_0^2+{\displaystyle \frac{n(x_0){\cos }^2(\theta /2)}{{\omega }_0^2}}}.

For two solitons generated at (x_{0,j},{\theta }_j), j=1,2, each follows its isolated-soliton trajectory [Eqs. (4) and (5)] before appreciable wave-function overlap. These isolated trajectories define a predicted encounter point (t_c^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}},x_c^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}) through x_1(t_c^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}})=x_2(t_c^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}})=x_c^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}, which provides the local inputs used below. The incoming velocity of soliton j at the predicted encounter position follows from Eq. (4):

(6)v_j(x_c^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}})=\mathrm{s}\mathrm{g}\mathrm{n}({\dot{x}}_j)\sqrt{v_{0,j}^2+{\omega }_0^2(x_{0,j}^2-(x_c^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}{)}^2)}.

Throughout this work, dark solitons are initialized using the product ansatz [38]

(7)\psi (x,0)={\psi }_0(x){\displaystyle \prod _{j=1}^2\left[\cos {\displaystyle \frac{{\theta }_j}{2}}-i{s}_j\sin {\displaystyle \frac{{\theta }_j}{2}}\tanh (\sin {\displaystyle \frac{{\theta }_j}{2}}(x-x_{0,j}))\right],}

where {\psi }_0(x) is the background condensate wave function and s_j=\pm 1 specifies the initial propagation direction of soliton j. This construction provides a clean two-soliton initial state with the phase structure of imprinted dark solitons, while avoiding the explicit modeling of density waves that may accompany experimental phase imprinting [6, 7, 38]. We impose an initial separation |x_{0,1}-x_{0,2}|\ge 5/\sin ({\theta }_{min}/2), where {\theta }_{min}=min({\theta }_1,{\theta }_2) corresponds to the broader of the two solitons and therefore sets the relevant width scale. For this separation, the interaction-induced corrections to the individual soliton parameters remain below 1%; we verified this by computing the inverse-scattering eigenvalues of the box-type two-soliton profile following Refs. [39–42]. Because the inter-soliton interaction remains weak before appreciable overlap, the isolated-soliton trajectories provide a reliable baseline for evaluating the incoming velocities at x_c^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}, although the nonlinear interaction during overlap may shift the GPE-identified collision center away from this predicted point.

Figure 1 confirms that the position-corrected map [Eqs. (2)–(5)] reproduces both the oscillation amplitude and the conserved depth of the imprinted soliton to within 0.5\mathrm{\%} over the full range of imprinting positions. This single-soliton mapping is used to determine the local velocity and background density that enter the collision criterion. In this sense, Eqs. (2)–(6) provide the kinematic input for the density-resolved description developed in Sections 3 and 4.

3 Two-soliton collision density on a uniform background

On a uniform background, the density at the collision instant follows from the exact two-soliton solution of the defocusing nonlinear Schrödinger equation. In this section, we evaluate this solution at the collision center to obtain a closed-form prediction for the collision-center density, denoted by {\rho }_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}, valid across both the fusion and non-fusion regimes.

3.1 Collision-center density formula

We begin with the exact two-dark-soliton solution of the defocusing nonlinear Schrödinger equation on a uniform background of constant density n_0, following the inverse-scattering theory of the defocusing NLS equation [25, 26] and the explicit notation of Frantzeskakis [3] [Eqs. (69)–(71) therein]. This solution takes its most compact form for symmetric pairs with identical depths and opposing velocities, v_1=-v_2=v. In this case, the collision-center density depends on a single dimensionless speed parameter A=v/c_s, with c_s=\sqrt{n_0} the homogeneous sound speed.

For asymmetric collision configurations where |v_1|\ne |v_2|, the general two-dark-soliton solution shows that the collision-center density depends on the relative velocity of the two solitons (see Appendix A). We therefore introduce

(8)A_{\mathrm{a}\mathrm{v}\mathrm{g}}={\displaystyle \frac{|v_1-v_2|}{2c_s}}

as the generalized normalized relative-speed parameter.

The resulting predicted collision-center density can then be written as

(9){\rho }_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}=n_0(2A_{\mathrm{a}\mathrm{v}\mathrm{g}}-1{)}^2.

This expression has distinct density-profile interpretations on the two sides of the threshold A_{\mathrm{a}\mathrm{v}\mathrm{g}}=1/2. For A_{\mathrm{a}\mathrm{v}\mathrm{g}}>1/2 (fusion), the two soliton density notches merge into a single composite dip whose minimum density is {\rho }_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}. For (non-fusion), they remain spatially separated as two density minima enclosing a central bump of height {\rho }_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}. At A_{\mathrm{a}\mathrm{v}\mathrm{g}}=1/2, the collision-center density vanishes ({\rho }_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}=0). The critical velocity v_c=c_s/2, commonly used to characterize the fusion–non-fusion threshold for symmetric dark-soliton collisions on a homogeneous background [1, 3, 22, 28], is recovered as the symmetric homogeneous limit of the more general condition A_{\mathrm{a}\mathrm{v}\mathrm{g}}=1/2.

3.2 Numerical verification on a uniform background

To validate Eq. (9), we integrate the dimensionless GPE [Eq. (1)] in the uniform limit (\mathrm{\Omega }=0) using a second-order Strang-splitting (split-step Fourier) scheme [36]. No random noise is added to the initial state or during the real-time propagation. The computational domain is a periodic grid of total length L=400 with N=8192 points, giving a spatial resolution \mathrm{\Delta }x\approx 0.049. The time step is \mathrm{\Delta }t={10}^{-3}, and the full density profiles are recorded every 20\mathrm{\Delta }t.

For the asymmetric counterpropagating initial states considered in the uniform-background calculations, the phases at the two ends of the Fourier grid generally differ by \mathrm{\Delta }\varphi ={\theta }_2-{\theta }_1. To ensure periodicity, we apply a uniform phase gradient \exp [\mathrm{i}{q}_{\mathrm{b}\mathrm{g}}(x+L/2)], with q_{\mathrm{b}\mathrm{g}}=-\mathrm{\Delta }\varphi /L. This transformation leaves the density and relative soliton velocity unchanged, and the induced uniform background flow is removed by transforming positions back to the zero-background-flow frame.

Figure 2 shows representative collision density profiles on a uniform background for the critical and non-fusion cases. At the critical point (A_{\mathrm{a}\mathrm{v}\mathrm{g}}=1/2), the collision-center density vanishes, giving a single zero-density node [Fig. 2(c)]. In the non-fusion case, two resolved density minima remain separated by a finite central bump [Fig. 2(d)]. These profiles provide direct signatures of the distinct collision behaviors described above.

The density profiles in the bottom panels of Fig. 2 are taken at the GPE-identified collision time t_c^{\mathrm{G}\mathrm{P}\mathrm{E}}. The collision frame is identified directly from the GPE density evolution, with t_c^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}} used only to define a temporal search window. Within this window, t_c^{\mathrm{G}\mathrm{P}\mathrm{E}} is taken as the frame at which the merged dip reaches its minimum density (fusion), or at which the central bump between the two resolved minima reaches its lowest height (non-fusion). The corresponding collision position and collision-center density are denoted by x_c^{\mathrm{G}\mathrm{P}\mathrm{E}} and {\rho }_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}}, respectively.

Using the same extraction procedure for all cases, Table 1 compares the measured collision-center densities with the prediction of Eq. (9). The agreement is quantitative over the fusion, non-fusion, and critical regimes: for all six cases, the absolute difference between {\rho }_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}} and {\rho }_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}} remains at or below 9.0\times {10}^{-4}. This agreement confirms that Eq. (9) accurately describes the collision-center density over both sides of the fusion–non-fusion threshold.

4 Two-soliton collisions in a harmonic trap

To apply the homogeneous result of Eq. (9) to a harmonically trapped condensate, we evaluate the relevant quantities locally at the predicted encounter position x_c^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}. Within the local-density approximation, the background density and sound speed are n(x_c^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}})=1-{\omega }_0^2(x_c^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}{)}^2 and c_s(x_c^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}})=\sqrt{n(x_c^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}})}, respectively, while the incoming soliton velocities are evaluated at the same position. The resulting local relations are

(10a)\begin{array}{c}{\rho }_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}=n(x_c^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}){\left[2A_{\mathrm{a}\mathrm{v}\mathrm{g},\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}(x_c^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}})-1\right]}^2,\end{array}

(10b)\begin{array}{c}{A}_{\mathrm{a}\mathrm{v}\mathrm{g},\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}(x_c^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}})={\displaystyle \frac{|v_1(x_c^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}})-v_2(x_c^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}})|}{2c_s(x_c^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}})}}.\end{array}

As in the homogeneous case, the local fusion–non-fusion threshold is A_{\mathrm{a}\mathrm{v}\mathrm{g},\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}(x_c^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}})=1/2, while Eq. (10a) gives the corresponding collision-center density. To test the local prediction, we perform trapped simulations with \mathrm{\Omega }=0.01, using the same spatial and temporal discretization as in the uniform-background calculations. Two solitons are imprinted on the numerical ground state {\psi }_0(x) using Eq. (7), and {\rho }_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}} is extracted following the procedure described in Section 3.2. The quantities entering Eqs. (10a) and (10b) are evaluated at x_c^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}, whereas {\rho }_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}} is obtained at the collision center identified directly from the Gross–Pitaevskii simulation.

Table 2 summarizes ten trapped-collision configurations. The first three entries form a \theta =0.70\pi sequence with x_{0,1}=-x_{0,2}=-x_0 and x_0\in \{10,35,50\}. On a uniform background, this phase corresponds to A_{\mathrm{a}\mathrm{v}\mathrm{g}}=0.454 and therefore lies in the non-fusion regime [cf. Table 1, row 3]. In the trap, however, A_{\mathrm{a}\mathrm{v}\mathrm{g},\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}(x_c^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}) increases from 0.458 at x_0=\pm 10 to 0.553 at x_0=\pm 50, crossing the local threshold near x_0=\pm 35. Thus, varying the imprinting position at fixed \theta =0.70\pi causes the local collision parameter to cross the fusion–non-fusion threshold.

The next two entries correspond to the same asymmetric soliton pair ({\theta }_1=0.678\pi, {\theta }_2=0.508\pi) in head-on and catch-up configurations. These two cases show how different encounter positions in the inhomogeneous condensate lead the same imprinted soliton pair to sample different local collision conditions. The contrast between these two configurations is illustrated in Fig. 3. In the head-on configuration [Fig. 3(a, c, e)], the predicted encounter occurs near the trap center, where n(x_c^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}})=0.997 and A_{\mathrm{a}\mathrm{v}\mathrm{g},\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}(x_c^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}})=0.598, and the Gross–Pitaevskii simulation produces a fused density profile. In the catch-up configuration [Fig. 3(b, d, f)], the predicted encounter lies in the outer region of the trap, where n(x_c^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}})=0.809 and A_{\mathrm{a}\mathrm{v}\mathrm{g},\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}(x_c^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}})=0.169, and the two density minima remain separated. The two outcomes therefore differ only because the pair samples different local conditions at the two encounter positions.

The quantitative deviations in Table 2 also reveal a clear dependence on the collision geometry. For the eight head-on configurations, Eq. (10a) reproduces the measured collision-center density to within approximately {10}^{-4}. The two catch-up configurations show larger deviations, with \mathrm{\Delta }\rho =+0.0126 and +0.0065, respectively. They also exhibit appreciable shifts of the GPE-identified collision event relative to the isolated-trajectory prediction: (\mathrm{\Delta }x,\mathrm{\Delta }t)=(-2.35,-6.42) for the first catch-up case and (-0.82,-4.36) for the second.

These deviations are consistent with the expected range of validity of the local construction. The off-center catch-up collisions sample a more rapidly varying background, while their smaller relative velocity leads to a longer period of soliton overlap. The mutual interaction can therefore accumulate over a longer interval, producing a larger displacement of the actual GPE collision event from the encounter predicted by the isolated-soliton trajectories. Consequently, the local approximation is most accurate when the background varies slowly over the collision region and the overlap time is short.

To assess the local scaling over a broader range of collision configurations, Fig. 4 compares the normalized collision-center density for all data sets. Dividing Eq. (10a) by the local background density gives

(11)\begin{array}{c}{\displaystyle \frac{{\rho }_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}}{n(x_c^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}})}}={\left[2A_{\mathrm{a}\mathrm{v}\mathrm{g},\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}(x_c^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}})-1\right]}^2.\end{array}

The data set combines the cases in Tables 1 and 2 with a 46-case symmetric grid scan over \theta /\pi =0.50–0.80 and |x_0|=10–50, spanning A_{\mathrm{a}\mathrm{v}\mathrm{g},\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}\simeq 0.17–0.95. Despite the different collision geometries and the larger deviations of the slow catch-up cases, the normalized densities follow the same functional dependence. This collapse supports the local scaling of Eq. (11) across the range of configurations considered. The sampled interval reflects the configurations included in the present simulations rather than a theoretical bound.

5 Conclusion

We have shown that the fusion–non-fusion behavior of colliding dark solitons in an inhomogeneous Bose–Einstein condensate depends on the local collision conditions rather than on an individual soliton velocity alone. These conditions are set by the local background density and the incoming relative velocity at the predicted encounter position x_c^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}. By evaluating the exact homogeneous two-soliton solution at the collision instant and applying a local-density replacement, we obtain

{\rho }_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}=n(x_c^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}})[2A_{\mathrm{a}\mathrm{v}\mathrm{g},\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}(x_c^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}})-1{]}^2. This expression gives the minimum density of the merged dip in the fusion regime and the height of the central bump between two resolved density minima in the non-fusion regime. In the symmetric homogeneous limit, the local threshold A_{\mathrm{a}\mathrm{v}\mathrm{g},\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}=1/2 reduces to the standard critical-velocity condition v_c=c_s/2.

Gross–Pitaevskii simulations on both uniform and trapped backgrounds confirm the prediction for symmetric and asymmetric soliton pairs and for both head-on and catch-up collisions. After normalization by n(x_c^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}), the measured collision-center densities collapse onto the common curve (2A_{\mathrm{a}\mathrm{v}\mathrm{g},\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}-1{)}^2. This collapse demonstrates that the collision-center density follows a common local scaling set by the background density and the incoming relative velocity.

A direct consequence is the position dependence of the collision outcome in a trapped condensate. The same imprinted soliton pair can exhibit either fusion or non-fusion under different position–velocity configurations because the trap changes the local background density and incoming velocities sampled at the predicted encounter position. By contrast, on a homogeneous background these quantities are spatially uniform, and the collision outcome is independent of the meeting position.

The predicted position dependence and the collapse of the normalized collision-center density should be accessible in quasi-one-dimensional {}^{87}\mathrm{R}\mathrm{b} or {}^{23}\mathrm{N}\mathrm{a} condensates using high-resolution in situ density imaging and established dark-soliton characterization techniques [43, 44]. Measurements of {\rho }_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}} for controlled initial configurations would provide a direct test of the local scaling identified here.

An interesting extension is to connect the density-resolved collision structure with the interaction-induced trajectory shifts accumulated during the collision. In our simulations, slow catch-up collisions exhibit larger deviations from the isolated-soliton trajectory prediction than the head-on cases considered here, consistent with their prolonged overlap. Developing an effective inter-soliton description that incorporates the collision-density structure may provide a bridge between density-resolved and trajectory-resolved descriptions.

6 Appendix A: Derivation of the collision-center density formula

We derive Eq. (9) from the exact two-dark-soliton solution of the homogeneous defocusing NLS equation, following the notation of Frantzeskakis [3] [Eqs. (69)–(71) therein]. For a symmetric pair of solitons with velocities v_1=-v_2=v=A\sqrt{n_0} on a uniform background of density n_0, the two-soliton wave function reads

(A.1)\begin{array}{c}\psi (x,t)={\displaystyle \frac{F(x,t)}{G(x,t)}}{\mathrm{e}}^{-\mathrm{i}\mu t},\end{array}

with

(A.2)\begin{array}{c}\begin{array}{rl}F(x,t)=& 2(n_0-2n_{min})\cosh (2n_{0}ABt)\\ & -2n_{0}A\cosh (2\sqrt{n_0}Bx)+\mathrm{i}\sinh (2n_{0}ABt),\end{array}\end{array}

(A.3)\begin{array}{c}G(x,t)=2\sqrt{n_0}\cosh (2n_{0}ABt)+2\sqrt{n_{min}}\cosh (2\sqrt{n_0}Bx),\end{array}

where B=\sqrt{1-A^2} is the depth parameter and n_{min}=n_0{A}^2 is the density at the bottom of an individual soliton notch. The collision occurs at x=0 and t=0. Evaluating Eqs. (A.2) and (A.3) at this point gives

\begin{array}{rl}& F(0,0)=2n_0-4n_0{A}^2-2n_{0}A=-2n_0(2A-1)(A+1),\\ & G(0,0)=2\sqrt{n_0}(1+A),\end{array}

so that

\psi (0,0)={\displaystyle \frac{F(0,0)}{G(0,0)}}=-\sqrt{n_0}(2A-1),

and

(A.4)\begin{array}{c}{\rho }_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}=|\psi (0,0){|}^2=n_0(2A-1{)}^2,\end{array}

which gives the collision-center density in Eq. (9) for the symmetric case. For asymmetric pairs, we use the general two-dark-soliton solution of the defocusing NLS equation with nonvanishing boundary conditions [45]. Introducing the signed normalized soliton velocities

a_j={\displaystyle \frac{v_j}{c_s}}={\displaystyle \frac{v_j}{\sqrt{n_0}}},j=1,2,

evaluation of the exact two-soliton solution at the interaction center gives

{\displaystyle \frac{|{\psi }_{\mathrm{c}}|}{\sqrt{n_0}}}=\left|1-|a_1-a_2|\right|.

Therefore, the collision-center density is

(A.5)\begin{array}{c}{\rho }_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}=n_0{\left(1-|a_1-a_2|\right)}^2.\end{array}

Since

A_{\mathrm{a}\mathrm{v}\mathrm{g}}={\displaystyle \frac{|v_1-v_2|}{2c_s}}={\displaystyle \frac{|a_1-a_2|}{2}},

Eq. (A.5) becomes

{\rho }_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}=n_0{\left(2A_{\mathrm{a}\mathrm{v}\mathrm{g}}-1\right)}^2,

which recovers Eq. (9) for a general asymmetric pair. Because a_j are signed velocities, the same expression applies to both counterpropagating and copropagating collisions.

References

[1] T. Tsuzuki , Nonlinear waves in the Pitaevskii–Gross equation, J. Low Temp. Phys. 4(4), 441 (1971)

[2] Yu. S. Kivshar and B. Luther-Davies , Dark optical solitons: Physics and applications, Phys. Rep. 298(2), 81 (1998)

[3] D. J. Frantzeskakis , Dark solitons in atomic Bose–Einstein condensates: from theory to experiments, J. Phys. A Math. Theor. 43(21), 213001 (2010)

[4] D. Krökel , N. J. Halas , G. Giuliani , and D. Grischkowsky , Dark-pulse propagation in optical fibers, Phys. Rev. Lett. 60(1), 29 (1988)

[5] S. A. Gredeskul and Yu. S. Kivshar , Generation of dark solitons in optical fibers, Phys. Rev. Lett. 62(8), 977 (1989)

[6] S. Burger , K. Bongs , S. Dettmer , W. Ertmer , K. Sengstock , A. Sanpera , G. V. Shlyapnikov , and M. Lewenstein , Dark solitons in Bose–Einstein condensates, Phys. Rev. Lett. 83(25), 5198 (1999)

[7] J. Denschlag , J. E. Simsarian , D. L. Feder , C. W. Clark , L. A. Collins , J. Cubizolles , L. Deng , E. W. Hagley , K. Helmerson , W. P. Reinhardt , S. L. Rolston , B. I. Schneider , and W. D. Phillips , Generating solitons by phase engineering of a Bose–Einstein condensate, Science 287(5450), 97 (2000)

[8] A. R. Fritsch , M. Lu , G. H. Reid , A. M. Piñeiro , and I. B. Spielman , Creating solitons with controllable and near-zero velocity in Bose–Einstein condensates, Phys. Rev. A 101(5), 053629 (2020)

[9] A. Weller , J. P. Ronzheimer , C. Gross , J. Esteve , M. K. Oberthaler , D. J. Frantzeskakis , G. Theocharis , and P. G. Kevrekidis , Experimental observation of oscillating and interacting matter wave dark solitons, Phys. Rev. Lett. 101(13), 130401 (2008)

[10] L. Z. Meng , L. C. Zhao , T. Busch , and Y. Zhang , Controlling dark solitons on the healing length scale, J. Phys. At. Mol. Opt. Phys. 57(14), 145302 (2024)

[11] Z. X. Niu and C. Gao , Dynamical nonlinear excitations induced by interaction quench in a two-dimensional box-trapped Bose–Einstein condensate, Chin. Phys. B 33(2), 020314 (2024)

[12] Y. Zhong , H. Triki , and Q. Zhou , Dynamics of ring dark solitons and the following vortices in spin-1 Bose–Einstein condensates, Chin. Phys. Lett. 41(7), 070501 (2024)

[13] P. H. Lu , X. F. Zhang , and C. Q. Dai , Dynamics and formation of vortices collapsed from ring dark solitons in a two-dimensional spin–orbit coupled Bose–Einstein condensate, Front. Phys. (Beijing) 17(4), 42501 (2022)

[14] X. Y. Yang , S. L. Chen , L. X. Wang , X. F. Zhang , and R. F. Dong , Dynamics of solitons in ultracold Bose gases with tunable interactions, Front. Phys. (Beijing) 21(4), 043202 (2026)

[15] X. Zhu , M. R. Belić , D. Mihalache , J. Li , and L. Zeng , Dark solitons in purely quintic nonlinear lattices with defects, Chaos, Solitons & Fractals 200, 116983 (2025)

[16] J. Chen , D. Mihalache , M. R. Belić , W. Qin , D. Zhu , X. Zhu , and L. Zeng , Dark gap soliton families in coupled nonlinear Schrödinger equations with linear lattices, Nonlinear Dyn. 113(9), 10307 (2025)

[17] Th. Busch and J. R. Anglin , Motion of dark solitons in trapped Bose–Einstein condensates, Phys. Rev. Lett. 84(11), 2298 (2000)

[18] V. V. Konotop and L. Pitaevskii , Landau dynamics of a grey soliton in a trapped condensate, Phys. Rev. Lett. 93(24), 240403 (2004)

[19] V. A. Brazhnyi , V. V. Konotop , and L. P. Pitaevskii , Dark solitons as quasiparticles in trapped condensates, Phys. Rev. A 73(5), 053601 (2006)

[20] D. E. Pelinovsky , D. J. Frantzeskakis , and P. G. Kevrekidis , Oscillations of dark solitons in trapped Bose–Einstein condensates, Phys. Rev. E 72(1), 016615 (2005)

[21] G. Theocharis , A. Weller , J. P. Ronzheimer , C. Gross , M. K. Oberthaler , P. G. Kevrekidis , and D. J. Frantzeskakis , Multiple atomic dark solitons in cigar-shaped Bose–Einstein condensates, Phys. Rev. A 81(6), 063604 (2010)

[22] G. Huang , M. G. Velarde , and V. A. Makarov , Dark solitons and their head-on collisions in Bose–Einstein condensates, Phys. Rev. A 64(1), 013617 (2001)

[23] P. Sprenger , M. A. Hoefer , and G. A. El , Hydrodynamic optical soliton tunneling, Phys. Rev. E 97(3), 032218 (2018)

[24] S. K. Ivanov and A. M. Kamchatnov , Motion of dark solitons in a non-uniform flow of Bose–Einstein condensate, Chaos 32(11), 113142 (2022)

[25] V. E. Zakharov and A. B. Shabat , Interaction between solitons in a stable medium, Sov. Phys. JETP 37, 823 (1973)

[26] L. D. Faddeev and L. A. Takhtajan, Hamiltonian Methods in the Theory of Solitons, Springer, Berlin and New York, (1987)

[27] L. Zeng , B. A. Malomed , D. Mihalache , J. Li , and X. Zhu , Collisions and fusion of one- and two-dimensional solitons driven by potential troughs in the cubic–quintic nonlinear Schrödinger equations, Chaos 36(1), 013138 (2026)

[28] D. M. Jezek , P. Capuzzi , and H. M. Cataldo , Dark-soliton collisions in a toroidal Bose–Einstein condensate, Phys. Rev. A 93(2), 023601 (2016)

[29] A. Görlitz , J. M. Vogels , A. E. Leanhardt , C. Raman , T. L. Gustavson , J. R. Abo-Shaeer , A. P. Chikkatur , S. Gupta , S. Inouye , T. Rosenband , and W. Ketterle , Realization of Bose–Einstein condensates in lower dimensions, Phys. Rev. Lett. 87(13), 130402 (2001)

[30] M. Olshanii , Atomic scattering in the presence of an external confinement and a gas of impenetrable bosons, Phys. Rev. Lett. 81(5), 938 (1998)

[31] F. Dalfovo , S. Giorgini , L. P. Pitaevskii , and S. Stringari , Theory of Bose–Einstein condensation in trapped gases, Rev. Mod. Phys. 71(3), 463 (1999)

[32] P. G. Kevrekidis, D. J. Frantzeskakis, and R. Carretero-González, The Defocusing Nonlinear Schrödinger Equation: From Dark Solitons to Vortices and Vortex Rings, Volume 92 of Other Titles in Applied Mathematics, Society for Industrial and Applied Mathematics, Philadelphia, PA, (2015)

[33] X. F. Zhang , X. H. Hu , D. S. Wang , X. X. Liu , and W. M. Liu , Dynamics of Bose–Einstein condensates near Feshbach resonance in external potential, Front. Phys. (Beijing) 6(1), 46 (2011)

[34] F. Dalfovo and S. Stringari , Bosons in anisotropic traps: Ground state and vortices, Phys. Rev. A 53(4), 2477 (1996)

[35] M. L. Chiofalo , S. Succi , and M. P. Tosi , Ground state of trapped interacting Bose–Einstein condensates by an explicit imaginary-time algorithm, Phys. Rev. E 62(5), 7438 (2000)

[36] W. Bao , I. L. Chern , and F. Y. Lim , Efficient and spectrally accurate numerical methods for computing ground and first excited states in Bose–Einstein condensates, J. Comput. Phys. 219(2), 836 (2006)

[37] B. Wu , J. Liu , and Q. Niu , Controlled generation of dark solitons with phase imprinting, Phys. Rev. Lett. 88(3), 034101 (2002)

[38] L. D. Carr , J. Brand , S. Burger , and A. Sanpera , Dark-soliton creation in Bose–Einstein condensates, Phys. Rev. A 63(5), 051601 (2001)

[39] G. Biondini and B. Prinari , On the spectrum of the Dirac operator and the existence of discrete eigenvalues for the defocusing nonlinear Schrödinger equation, Stud. Appl. Math. 132(2), 138 (2014)

[40] A. Romero-Ros , G. C. Katsimiga , P. G. Kevrekidis , B. Prinari , G. Biondini , and P. Schmelcher , On-demand generation of dark soliton trains in Bose–Einstein condensates, Phys. Rev. A 103(2), 023329 (2021)

[41] B. Prinari , M. J. Ablowitz , and G. Biondini , Inverse scattering transform for the vector nonlinear Schrödinger equation with nonvanishing boundary conditions, J. Math. Phys. 47(6), 063508 (2006)

[42] J. H. Kang and S. Y. Zhang , Deterministic generation of dark soliton trains in Bose–Einstein condensates, Phys. Rev. A 113(1), 013330 (2026)

[43] C. Becker , S. Stellmer , P. Soltan-Panahi , S. Dörscher , M. Baumert , E. M. Richter , J. Kronjäger , K. Bongs , and K. Sengstock , Oscillations and interactions of dark and dark–bright solitons in Bose–Einstein condensates, Nat. Phys. 4(6), 496 (2008)

[44] H. A. J. Middleton-Spencer , N. G. Parker , L. Galantucci , and C. F. Barenghi , Characterising arbitrary dark solitons in trapped one-dimensional Bose–Einstein condensates, Europhys. Lett. 133(1), 16001 (2021)

[45] H. Blas and M. Zambrano , Quasi-integrability in the modified defocusing non-linear Schrödinger model and dark solitons, J. High Energy Phys. 2016(3), 5 (2016)

Rights & permissions

Higher Education Press