Non-integrable dynamics of matter-wave solitons in a density-dependent gauge theory

We study interactions between bright matter-wave solitons which acquire chiral transport dynamics due to an optically-induced density-dependent gauge potential. Through numerical simulations, we find that the collision dynamics feature several non-integrable phenomena, from inelastic collisions including population transfer and radiation losses to short-lived bound states and soliton fission. An effective quasi-particle model for the interaction between the solitons is derived by means of a variational approximation, which demonstrates that the inelastic nature of the collision arises from a coupling of the gauge field to velocities of the solitons. In addition, we derive a set of interaction potentials which show that the influence of the gauge field appears as a short-range potential, that can give rise to both attractive and repulsive interactions.


I. INTRODUCTION
One of the defining properties of solitons in systems such as the nonlinear Schrödinger (NLS) and Kortewegde Vries equations is that they pass through and emerge from collisions with other solitons unperturbed, with the exception of a phase shift arising from the nonlinear interaction. [1,2]. All of the dynamical quantities of solitons, such as their velocities and masses, are conserved during the collision. In such systems, the elastic nature of these collisions is a consequence of the integrability of the model, which heavily restricts the allowed dynamics due to the existence of an infinite set of conservation laws.
In non-integrable systems, stable solitons may exist too, but their collisions are, generally, inelastic and can lead to trajectories which are chaotic [1,3]. In this case, the defining feature is the existence of short-lived bound states in which the number of collision events depends fractally on the initial conditions. Generally, this mechanism arises from the excitation of an internal mode of the solitons in non-integrable models, either with [4][5][6][7][8] or without radiation losses [9,10], or through the presence of a weak perturbation [11][12][13]. Solitons can merge or fracture into new products through fission and fusion processes [14][15][16], which has also been studied in the context of three-soliton and soliton-breather collisions [17]. These effects highlight a strong contrast to soliton dynamics in integrable systems, which are not only interesting from a fundamental point of view, but offer insight into the description of realistic systems where the influence of perturbations can be consequential.
Ultracold atomic gases represent an attractive platform to study nonlinear physics due to their unprecedented experimental controllability. The ability to precisely engineer both the dimensionality and interactions in these systems has lead to the realization of isolated bright matter-wave solitons [18,19], as well as solitons trains [20]. The second generation of experiments addressed controlled collisions with potential barriers [21,22], as well as understanding both the role and origin of the relative phase for the stability of bright soliton states [23][24][25].
The underlying integrability of the focusing nonlinear Schrödinger equation leads to a hierarchy of analytical higher-order soliton solutions provided by the inverse scattering transform. This directly led to the development of a classical particle model [26], describing the dynamics of the bright solitons, from which regions of chaotic behaviour have been eventually predicted for trapped bright solitons in non-integrable settings [27,28].
Due to their inherent coherence, bright solitons represent a useful tool for investigating interferometry in the quantum realm. Recent experimental progress in this direction has seen the first realization of a 85 Rb matterwave Interferometer, as well as theoretical proposals for precision measurement using the Sagnac effect [29] and the creation of Bell states using quantum bright solitons [30], as well taking advantage of the interaction of solitons with nonlinear splitters [31].
The ability to simulate artificial gauge theories with ultracold gases offers a new opportunity to understand the interplay of effective magnetism in these systems with nonlinear effects [32,33]. This has led to the realization of vortex states [34] as well as spin-orbit coupling [35]. Interest was recently focused on schemes for generating gauge potentials with an effective back-action between the matter and the gauge potential [36,37], which leads to a number of novel phenomena, including the violation of the Kohn's theorem [38,39], and unconventional vortex dynamics [40,41]. Very recently, the first experimen-arXiv:1709.08037v2 [cond-mat.quant-gas] 22 Feb 2018 tal realization of a dynamical gauge theory in a trapped ion system was shown [42].
In this paper, we study the nonlinear dynamics of two interacting one-dimensional chiral matter-wave solitons. We begin by reviewing how these solitons can be engineered in ultracold gases using optical techniques which induce an effective density-dependent gauge potential in the atomic cloud. The resulting equation of motion for the gas, which takes the form of a chiral nonlinear Schrödinger equation, is then solved numerically in Sec. III. Following this in Sec. IV, we develop a variational approach to further understand the soliton dynamics in this system, both in linear and asymptotic limits, before concluding in Sec. V.

II. THE THEORETICAL MODEL
We consider a Bose-Einstein condensate of N two-level interacting atoms, in which two internal states of the atoms (labelled |1 and |2 ) are resonantly coupled by an external laser field. The Hamiltonian describing the interacting trapped gas as well as the light-matter interaction can be written aŝ whereĤ lm = Ω 2 0 e −iφ e iφ 0 (2) describes the optical coupling of the two internal states of the atoms, with strength Ω and laser phase φ . In order to obtain an equation of motion for the many-particle system, the state of the system is defined as a Hartree product |Ψ = ⊗ N l=1 |χ (0) l , where |χ (0) l defines one of the single-particle eigenstates of Eq. (2). In this work we assume that the gas is harmonically trapped in the (x, y) plane (described by the vectorial r ⊥ coordinates in Eq. (2)), but free along the axial z direction. The mean-field interactions appearing in Eq.
where |Φ i | 2 denotes the population of the state i.
Provided that the gas is sufficiently dilute, we can diagonalize the Hamiltonian by treating the mean-field interactions g ii |Φ i | 2 as a small perturbation to the laser coupling Ω. The eigenvectors of Eq. (1) can then be written in the dressed state basis {+, −} as where |χ (0) ± = |1 ± e iφ |2 / √ 2 denotes the unperturbed dressed states. The associated eigenvalues are given by g|Φ ± | 2 ± Ω/2, with the dressed scattering parameter g = (g 11 + g 22 + 2g 12 )/4. From the interacting dressed states, we can write the state vector of the system as |ξ = i=+,− Φ i (r, t) |χ i . Then, the effective Hamiltonian is written aŝ Equation (4) introduces the geometric phase A ± = i χ ± |∇χ ± . Accompanying this is a scalar geometric phase, whose leading-order effect is inducing an energy offset, which may be dropped. Then, using the definition given by Eq. (3), to lowest-order the density-dependent geometric phase appears as A ± = A (0) +a 1 |Φ ± (r)| 2 , with the single-particle vector potential, A (0) = − 2 ∇φ l (r), while a 1 = ∇φ l (r)(g 11 − g 22 )/8Ω defines the strength of the density-dependent gauge potential. The equation of motion governing the evolution of the wave-function amplitude Φ + (r, t) is found from minimization of the system's energy functional, E = Ψ|(i ∂ t −Ĥ ± )|Ψ . After dropping ± subscripts, the resulting mean-field densitydependent Gross-Pitaevskii equation is obtained as defines the current nonlinearity appearing in Eq. (5). The current-coupled nonlinear Schrödinger equation, captured by Eq. (5), describes a novel nonlinear gauge theory where there is an effective back-action between the matter-field and the gauge potential [36]. This feedback ingredient of the system is somewhat similar to the local field effect, which affects a "soft" optical lattice trapping the condensate, that gives rise to various consequences, such as formation of bright solitons in the absence of contact interactions between atoms [43].

A. One-dimensional reduction
We are interested in studying solitary-wave solutions in the frameworks of the dimensionally reduced form of Eqs. (5) and (6). To do this, we assume the system is in the ground state of the transverse trap, such that one can factorize the wave function as Φ( is the transverse ground-state wave function, and l ⊥ = /mω ⊥ is the transverse harmonic length scale. Equation (5) is then reduced to an effective one-dimensional form, where both the scattering parameters a 1 = k(g 11 − g 22 )/(16πl 2 ⊥ Ω) and g 1D = g/(2πl 2 ⊥ ) have been scaled by the transverse area of the cloud, and the corresponding one-dimensional current nonlinearity is defined as In writing Eq. (7), we have defined the laser phase as φ = kx and subsequently eliminated the zeroth-order vector potential through a momentum boost. We can further simplify Eq. (7) by introducing the nonlinear phase transformation which acts to decouple the vector potential from the canonical momentum appearing in the one-dimensional Gross-Pitaevskii equation. Substituting Eq. (9) into Eq. (7) and (8) leads to the simplified equation, where j (x) = ( /m)Im(ψ * ∂ x ψ) is the gauge-transformed current operator. Equation (10) belongs to the class of derivative or 'chiral' NLS equations [44,45], which was originally studied, in particular, in the context of onedimensional anyons [46]. The model features several key differences from the standard NLS equation in that it is generally non-integrable, does not obey the Galilean invariance, and possesses chiral soliton solutions [45][46][47]. These properties are expected to contribute to unconventional soliton dynamics in the one-dimensional case. An experimental realization of the interacting gauge theory relies, basically, on two conditions, viz., an atomic species possessing long-lived excited states that the adiabatic motion of the atoms requires, as well as the spontaneous-emission rate that is negligible on the time scale of cold-atom experiments. A promising candidate that fulfils these conditions are alkali-earth atoms, which have been recently used to create a spin-orbit-coupled Fermi gas of 173 Yb atoms [48], using a methodology similar to that outlined here. Very recently, an interactioninduced synthetic gauge potential was realised experimentally in a Bose-Einstein condensate loaded into a modulated two-dimensional lattice [49].

B. Chiral solitons
Single-soliton solutions of Eq. (10) can be derived by first boosting into the moving frame via the transformation which consists of a Galilean transformation in which the stationary coordinates (x, t) and moving coordinates (x , t ) are connected by the translations, x → x −vt and t → t, with frame velocity v. The resulting differential equation for the real-valued wave function ϕ(x ) becomes in which the current is contained as j(x) = vϕ 2 . Integrating Eq. (12), and requiring that the wave function converges to ϕ(±∞) = 0 for g 1D = g 1D − 2a 1 v < 0, one finds the single bright-soliton solution, in which we have transformed back into the stationary frame. The chemical potential appearing in Eq. (12) is µ = −mg 2 1D N 2 /8 2 , with the amplitude factor 1/ √ 2b provided normalization N = 1. In contrast to the Gross-Pitaevskii theory used to model non-chiral bright solitons [50], the lack of the Galilean invariance of Eq. (10) results in a change of the soliton width, b = −2 2 /mg 1D , depending on the direction of the soliton's motion. These solitons are therefore chiral, in the sense that the meanfield interactions, and hence the soliton's size depend on the direction in which it is travelling. An illustrative example of this can be demonstrated by the reflection of a chiral soliton off a hard wall, which causes the soliton to disperse [36].

C. Conservation laws
Although the equations of motion defined by Eqs. (7) and (10) are non-integrable, a set of conservation laws in the present system can be derived directly from the Noether's theorem [1]. With the respective Lagrangian density, it is straightforward to show that, at least, three conserved quantities exist [51], given by the integral expressions: and which correspond to the number of particles, momentum, and energy of the system, respectively. The integrands of Eq. (17) introduce the Hamiltonian densities for both the transformed and non-transformed representations, which is expected due to the fact that the underlying Lagrangian is Hermitian. A specific peculiarity of the chiral model arises in Eq. (16), which shows that, as a consequence of the breakdown of the Galilean invariance, the canonical momentum is not conserved, but the quantity mj + a 1 |ψ| 4 is conserved.

III. NUMERICAL SIMULATIONS
The focussing NLS equation is an integrable model, where solitons collide elastically, with the same shape and velocity before and after scattering. The one-dimensional gauge theory defined by Eq. (10) breaks the integrability due to the presence of the current nonlinearity. In this section we numerically solve Eq. (10) for binary chiralsoliton collisions using the known single-soliton solution, given by Eq. (13), in order to understand how the broken integrability manifests itself. The system is prepared initially in the state where δ ∈ [−π, π] is the relative phase difference between the solitons, while ξ 1,2 are initial center-of-mass coordinates of the solitons, and the normalization condition is Since a full parameter scan of chiral-soliton collisions, featuring every degree of freedom, presents a formidable problem, we restrict our analysis to two parameter regimes, each set by a ratio of interaction strengths, which illustrate the essential physics present in the model: • In the first instance, we consider the case of strong chiral interactions, |g 1D | |a 1 (v 1 + v 2 ) |, where effects stemming from the current nonlinearity are made influential by the solitons' high velocities.
• For the second, we treat the case of weak chiral interactions, |g 1D | |a 1 (v 1 + v 2 ) |, where the current nonlinearity is treated as a small perturbation added to the usual mean-field dynamics.
To integrate Eq. (10) numerically, we construct an explicit central-difference scheme for the evolution of the wave function, and compare the results to those produced by a split-step Fourier method, to ensure consistency. The numerical domain is chosen to be at least two orders of magnitude larger than the widths of the solitons to avoid radiation back-reflecting into the solitons (the aliasing effect). It is useful at this stage to point out that the evolution of the soliton's relative phase depends on their separation [52]. Therefore, each result which we present is defined up to a choice of the initial phase difference and separation, although altering these parameters does not yield a qualitative difference. The colourbar limit has been intentionally lowered in (c)-(d) to display the solitons more clearly.

A. Strong interactions
In Fig. 1 we show a set of density plots for the collision of two co-moving chiral solitons in the presence of the strong chiral interactions. In each case, the interaction induced by the current nonlinearity dominates over the usual mean-field effects due to the solitons' high velocities, which effectively reduces the mean-field scattering parameter. Figure 1(a) and Fig. 1(b) highlight extreme manifestations of this regime, where the dynamics are solely influenced by the current nonlinearity, setting g 1D m / 2 = 0, i.e., the mean-field contact interaction is completely disregarded.
Surprisingly, the collisions in these two instances are similar to those produced by the conventional NLS dynamics, with the solitons surviving the collision and retaining the general shape of their envelopes. However, two key differences are visible in the trajectories of the solitons, particularly in the case of a 1 / = 4, Fig. 1 The first are a pair of inelastic trajectories, with velocities of the outgoing solitons differing from their initial velocities, with the left-and right-most solitons reducing and increasing their velocities, respectively. The second difference from the mean-field setting is the appearance of a density node at the interaction centre, which is reminiscent of a repulsive interaction, despite the initial phase difference taken as δ = 0, which is conventionally an attractive interaction (for the solitons colliding in-phase). A third feature is also spotted in the form of population transfer, as illustrated in the Fig. 1(b.i) for the collision pictured in Fig. 1(b). Here, approximately a quarter of the initial mass of the right-most soliton is transferred to the left-most soliton in the course of the collision. Note, that the populations of the solitons are not calculated during the collision at ( /m 2 )t = 4 to ( /m 2 )t = 12, due to the overlap of the soliton envelops.
Each of these effects can be traced to the nonintegrability of the current nonlinearity, which permits the transfer of stored interaction energy into the kinetic energy, and a possibility of a non-trivial shift of the soliton phase difference. The presence of the population transfer is directly linked to the phase shift, as both quantities are mutually conjugate. The energy exchange, or, to a greater extent, the inelasticity of the collision, appears to be minimised when the collision parameters are chosen so that the solitons interact repulsively. This is evident from comparisons between the two figures, where we note that a 1 / = 1 leads to a repulsive interaction with elastic trajectories, while a 1 / = 4 gives rise to a more attractive interaction which features inelastic trajectories.
The lower row of Fig. 1 shows the dynamics where the collision is destructive, causing fission of the solitons. In both cases, three solitons emerge from the collision (the third soliton in Fig. 1(d) with δ = π is located at the leading edge of the other two) with the populations and velocities of each outgoing soliton depending on the initial phase difference. In addition, a modest amount of radiation is ejected during the collision, as seen in the trailing edge in Fig. 1(c), and the interference pattern located between the slowest two solitons in Fig. 1(d). The main difference against the previous case is a larger difference in the initial velocities, which, if coupled with a larger gauge-field strength, a 1 / = 5, produces two soliton envelopes with a greater disparity of widths. As such, in the course of the collision, the solitons effectively interact over a longer period, thus enhancing effects stemming from the interaction.

B. Weak interactions
In the previous subsection it was seen how chiral solitons undergo inelastic collisions. To quantify the elasticity of the collisions, we introduce the coefficient of resti- tution [14] η = which compares the difference of the kinetic energy before and after the collision. For η = 1, the collision is perfectly elastic with conserved masses and velocities of the solitons, while η = 0 indicates an inelastic collision.
Here, m 1,2 and v 1,2 play the role of the masses and velocities of the solitons in our semi-classical description, and are calculated from the respective expectation values, The integration in each case at either the initial or final time is performed locally around each soliton's centre of mass to exclude contributions from radiation and overlap with the other soliton. Due to the occurrence of both the population transfer and changes in the outgoing velocities of the solitons, one cannot distinguish whether the chiral solitons pass through or rebound off each other during the collision. Therefore, to remain consistent, we denote the solitons located in the x < 0 and x > 0 regions as the first and second ones, respectively. By varying the strength of the gauge field, we have performed a detailed parameter scan of the soliton-soliton collisions as a function of the initial phase difference, with the initial soliton velocities fixed. The coefficient of restitution is computed and plotted in Fig. 2 with corresponding examples of the dynamics shown in Fig. 3. For each value of the gauge- field strength, three regimes of the collision dynamics can be identified, depending on the initial phase difference between the solitons. The first is an elastic scattering regime highlighted by a plateau in the restitution data at η = 1, with an example of the dynamics shown in Fig. 3(a). Here, the interaction is notably repulsive, with a distinct node in the density at the interaction centre and the soliton parameters keeping their values after the collision.
Away from this plateau, two distinct regimes of inelastic dynamics are found with η > 1, as illustrated in Fig. 3(b) and Fig. 3(d). Here, the dynamics are also similar to the case of strong interactions, with inelastic trajectories that feature a redistribution of the soliton masses, as well as evolution of the initial phase difference, resulting in shifts of the in-and out-of-phase collision points. Comparing these two plots, one notices that, depending on the direction in which one moves away from the plateau in the parameter space, the soliton mass can be transferred, chiefly, in either the left (Fig. 3(d)) or right-hand ( Fig. 3(b)) outgoing soliton.
The final inelastic regime, indicated by the 'resonance' peak in the restitution data represented by the cross labelled (c) in Fig. 2), features the turning point of this population transfer, where the mass is transferred from one soliton to the other. We show in Fig. 3(c) an example of the dynamics in this regime, at the maxima of the 'resonance' peak. Here, a peculiar soliton state is formed, where there is a strong interplay between emitted radiation, excitation of an oscillatory mode in the right-hand soliton, and a weak left-hand one. This regime appears to be an example of in-phase (fully attractive) dynamics, which features the formation of a metastable (short-lived) bound state.
A feature universal to the restitution data presented in Fig. 2 is that the location of each inelastic regime is cyclically shifted left-wards for an increasing current strength.
Comparing different gauge-potential strengths, one can see that the elastic region shrinks for larger values, which can be explained by enhancement of the non-integrability effects for stronger gauge-field strengths. Although not shown here, the dip in the restitution data initially appears close to δ = 0 at small values of the current strength, and cyclically displaces towards lower δ for increasing current strengths.
To complete the analysis for the weak-chiral regime, we perform a similar parameter scan as before, but now in the case when the relative phase difference is fixed to δ = 0, with the initial velocity of the left-hand soliton allowed to vary. In this case, the coefficient of restitution provides a poor illustration of the underlying dynamics, therefore we, instead, plot the outgoing velocity of the soliton travelling to the right for increasing values of the gauge-field strength, as shown in Fig. 4(a).
Depending on the choice of the initial velocity and gauge-field strength, the strength of the chiral interactions |a 1 (v 1 + v 2 )| may be either small or comparable to the mean-field strength |g 1D |. Therefore, for extreme values of the parameters, it is expected that the dynamics will be generally inelastic in a similar manner to Fig.  (1), whereas for smaller values the dynamics will be, generally, elastic. This reasoning is reflected in the pair of curves corresponding to a 1 / = 1 and a 1 / = 1.25 in Fig.  4(a), in which the soliton velocity does not change significantly after the interaction for small initial velocities. As the velocity increases (and hence the interaction strength increases too), this invariance begins to break, which is particularly notable in the case of a 1 / = 1.25, which exhibits a sinusoidal behaviour, at the velocity exceeding a critical value, v 1 ≈ 0.5.
As the gauge field strength is increased further, as in the case of a 1 / = 1.5, a 'resonance' feature appears in the data where a two-bound resonance state is formed, as shown in Fig. 4(b). As mentioned previously, such states are a common occurrence in non-integrable models [4][5][6][7][8][9][10][11][12][13], perhaps most notably in non-integrable versions of the sine-Gordon equation [8], where, depending on the strength and shape of the interaction potential, an nbound resonance state may emerge. Here the underlying mechanism is the energy exchange between the colliding solitons as a whole and their internal modes, which requires the solitons to collide several times before escaping, thus regaining the energy temporarily transferred into the internal mode. In performing this parameter scan, higher-order bound states, where the solitons collide more than twice, were not observed, as for stronger interaction strengths the appearance of a bound state tends to be suppressed in a similar manner to Fig. 3(c).

C. Bound states
To further investigate the inelastic dynamics of the density-dependent gauge theory, we consider a set of symmetric collisions (see Fig. 5), where two interacting solitons form a molecule-like bound state. In a similar manner to results obtained in weakly perturbed cubic-quintic NLS [7,11] and sine-Gordon systems [12], a weak current nonlinearity is found to support short-lived bound states, where the solitons collide several times before escaping. As before, the underlying mechanism here is the transfer of a part of the energy of the interaction between the solitons into the kinetic energy and redistribution of the soliton masses. However, the magnitude of the interaction energy is now comparable to or larger than the kinetic energy of the solitons, requiring them to collide several times in order to gain enough kinetic energy for escap-ing the attractive interaction. Compared to the standard Gross-Pitaevskii dynamics shown in Fig. 5(a) with a 1 / = 0, where the solitons are perpetually trapped with a fixed oscillation amplitude and frequency, a modest current strength can begin to destabilize the bound state, such as in Fig. 5(c) for a 1 / = 0.25, where the solitons collide four times before escaping, and in Fig. 5(d) for a 1 / = 0.5, where they collide twice before escaping.
Interestingly, despite the interaction being initially symmetric, effects stemming from the chiral dynamics result in a left-handedness in the post-collision behaviour. For example, in Fig. 5(b) with a 1 / = 0.125, the first collision at ( /m 2 )t = 10 is noticeably attractive due to the presence of the anti-node at the interaction centre, but every subsequent collision becomes increasingly repulsive with the amplitude of the anti-node decreasing and its position shifting towards the left. In addition, a density node fills the vacancy left by the anti-node at each interaction centre, with some manifestation of the population transfer. This effect is seen to be most profound in Fig. 5(d) for a 1 / = 0.5, where ∼ 60% of the outgoing mass in captured in the left soliton. Returning to the two-bound resonance state in Fig. 4, the existence of higher-order bound states appears unlikely due to the fact that each subsequent interaction becomes more repulsive, hence more elastic than the previous one. These dynamics highlight the role that the lack of the Galilean invariance plays in the interacting gauge theory. In particular, the current operator appearing in Eq. (10) induces the population transfer between the two solitons, resulting in the suppression and, ultimately, breakdown of the bound state.

IV. THE VARIATIONAL ANALYSIS
To gain insight into how the current nonlinearity modifies the interactions between the solitons, we have performed variational calculations to derive an effective particle model for the soliton dynamics. We achieve this by using two similar, but essentially different methods.
In the first instance, we approximate the two-soliton state as a linear superposition of two individual solitons, with the interaction treated as the spatial overlap of the soliton envelopes. This technique has been previously applied to interaction problems in the NLS [53,54] and Gross-Pitaevskii equations [55,56], in addition to several others [57][58][59][60][61]. An advantage of this method is the ability to derive a set of variational equations which describe the motion of the solitons. From this, two key results can be extracted. The first is an effective potential describing the interaction between the solitons, which will provide details into the phase dependence and range of the interactions. Secondly, by numerically solving the variational equations, we will be able to illustrate the dynamics of the particle model and compare it directly to the full numerical solutions that are presented above.
For the second method, we follow the technique out-lined in Refs. [62,63], in which the soliton state is also approximated as a linear superposition, but restricted to the case of two stationary solitons which are wellseparated. In this case, the interaction is accounted for by the spatial overlap of one soliton with the 'weak tail' of the other, and may therefore be regarded as an asymptotic approximation to the full interaction. As such, this method is well suited to the study of bound states, but will also provide a basis to draw comparisons to the first method in a far-field low-velocity limit.

A. Linear calculations
The starting point for the variational calculation is the Lagrangian density [47,64] where we have introduced dimensionless variables x =x , t =tm 2 / , ψ =ψ/ √ , g 1D =g 1D 2 /m , and a 1 =ã 1 , with the length scale . Here, and in the remainder of this paper, we drop the tildes appearing in these variables, working in the dimensionless units.
In what follows, we seek stationary solutions determined by the action which gives rise to the Euler-Lagrange equations for each variational parameter q n (t), namely, The variational ansatz is adopted as in which each soliton contains a spatially-varying phase [65,66], Here, a(t), b(t), ξ j (t), v j (t), and φ j (t) are time-dependent variational parameters corresponding to the amplitude, width, centre-of-mass coordinates, velocities, and central phases of the solitons. This ansatz models two bright solitons in which the individual velocities and positions are allowed to evolve independently, with the interaction treated as the (linear)-overlap of the soliton envelopes. The constraint that the solitons have a common width, which in turn fixes the profiles of the soliton envelopes, is a necessary restriction in order to be able to explicitly calculate interaction integrals. Consequently, this restricts the variational analysis to the regime in which This can be achieved by considering collisions with small velocities and by compensating the effects of the gauge field by a mean field with a modest strength, such that |a 1 (v 1 + v 2 )| |g|. Therefore, the above constraint restricts our variational analysis to the weak-chiral regime for which numerical results are presented above. In spite of these restrictions, we will find that one may be quite liberal with the choice of parameters and still achieve sensible results.
Our choice of ansatz arises due to two reasons. First, our model is non-integrable, therefore a closed-formed expression for a two-soliton state via inverse scattering techniques is not available. Secondly, regardless of whether such a solution existed, Eq. (25) should work as a good approximation to the dynamics pictured in Fig.  (3), as the solitons roughly retain their shape during the interaction. However, it must be stressed that this choice of the ansatz does not fully replicate all the features of the interaction and will therefore lead to inconsistencies at short length scales, when the solitons begin to significantly overlap.
An important example of this which must be considered before proceeding, is related to a divergence of the soliton amplitude a(t) at short length scales. This can be illustrated by evaluating Eq. (15) for our variational ansatz, from which one can obtain the following expression for the soliton's amplitude: To perform the integration, we have introduced the change of variables α = (x − ξ 1 )/b and α + = (x − ξ 2 )/b, where = (ξ 1 − ξ 2 )/b is a new variational parameter describing the relative positions of the solitons. Additionally, we have assumed that the magnitude of the velocity is small, such that the phase difference S 1 − S 2 ≈ δ is an approximate function of solely the central phases. Additional details of the calculation are outlined in Appendix. A. Although this integral can be evaluated exactly without needing this approximation [54], the ensuing exact expression is to cumbersome for extracting explicit results from it. Inspecting Eq. (28), one can see that, in the limit of → 0, the value of a 2 (t) rapidly diverges for δ ∈ [π/2, 3π/4] and approaches a singularity at δ = π, modulo 2π. Our approximations therefore lead to an unphysical divergence which is not representative of soliton collisions, thereby requiring us to restrict our studies to the bounded domain |f ( , δ)| ≥ 1, which corresponds to the interval of δ ∈ [0, π/2]. For such values of f which are bounded from below, a(t) is non-divergent, representing the known dynamics more adequately-for instance, at δ = 0, which corresponds to the fact that the solitons' amplitude is increased by 1/ √ 2b when they con-structively interfere. Substituting our ansatz into Eq. (22) and integrating via Eq. (23) leads to the averaged Lagrangians and where we have defined g 2 = g 1D − a 1 (v 1 + v 2 ), and have split the total Lagrangian into the sum of terms which implicitly and explicitly depend on , as denoted by the free and interacting Lagrangians, L = L 0 + L I .
Equations of motion for each variational parameter can now be derived from Eq. (24), which lead to the set of coupled differential equations a : with and Here, the vertical bar notation in Eqs. (33)(34), denotes the full Lagrangian function in Eq. (30), but excluding terms containing a factor ofφ i ,ξ i , orv i . In Eqs. (32) and (36), for i = 1 positive (+) operations are taken, with the converse for i = 2. A variational equation for b, the solitons' width, is not required to proceed and is excluded.
Both Eqs. (35) and (36) are introduced for notational convenience. From the set of variational equations, we can now extract details of how the gauge field affects the soliton dynamics, and derive several important quantities.
Starting with the first variational equation, which can be obtained by varying either φ 1 or φ 2 , one can identify Eq. (31) as a conservation law for 4a 2 bf ( , δ). This is consistent with both Eqs. (15) and (28), which state that the phase and density of the condensate are conjugate variables. In the asymptotic limit of → ∞, this conserved quantity reduces to a 2 ∼ 1/4b, which is the correct amplitude for a two-soliton state.
The equations for the velocities v i highlight the main result of the variational analysis. The first and last terms of Eq. (32) (and the last term of Eq. (35)), correspond to similar terms in the NLS equation. Together they imply that, in the asymptotic limit of v i ∼ξ i , both solitons move at a constant velocity when they are well separated. However, at → 0 the velocities of the solitons are modified due to their interaction, but, once again, they become constant after the solitons have passed through each other. The additional terms ∝ a 1 , which appear in the velocity equation, are new ones, which arise due to the presence of the gauge potential. The first of these is a non-Galilean effect that redefines the soliton velocities in the asymptotic limit as v i ∼ξ + a 1 /6b, which is consistent with the momentum conservation law stated in Eq. (16). The remaining terms which appear in Eq. (35) are responsible for the interaction-induced velocity shift, which, in both the v 1 and v 2 equations, has the same magnitude and sign.
The variational equations for both a and ξ i are not particularly transparent. However, they do highlight the coupling between all of the variational parameters, and will be required when deriving the interaction potentials later in the section.

Collision dynamics
In order to illustrate how the gauge field is the mechanism underlying the inelastic scattering in our system, we set out to first simplify and reduce the number of variational equations, so that an effective particle model can be derived. Subsequently, we can numerically solve our system of equations and compare it to the full numerics presented above.
We begin by first reiterating that we consider the case of weak-chiral interactions |a 1 (v 1 + v 2 )| |g|, which feature the solitons moving slowly for a given choice of the gauge-field strength. Equations (33-34) may then be set up as a set of simultaneous equations in which coordinates,φ 1 ,φ 2 ,ξ 1 , andξ 2 can be eliminated, leading to a pair of equationṡ x [ℓ] Together with Eq. (32), they form a set of coupled differential equations for the soliton dynamics, in which details of the interaction are encoded in expressions for γ( , δ) and L I . Once again, in the asymptotic limit, both equations simplify to the single-soliton resultv i ∼ 0, which highlights that both solitons move independently at a constant velocity when they are well separated. A consequence stemming from the elimination of variablesφ 1 andφ 2 in the equations highlights that the phase in this particle model is static, and does not dynamically evolve. Although this is an important feature in our model which has many consequences in the scattering dynamics, we will still be able to obtain qualitative results which do not strongly depend on phase δ, but will not be able to address issues pertaining to the bound-states dynamics pictured in Fig.  5.
We solve the set of differential equations numerically using a fourth-order Runge-Kutta method and compare our results to an example of the full numerics in Fig.  (6). For the chosen set of parameters, the magnitude of the outgoing velocities are in good agreement with the post-collision trajectories showing that the solitons pass through each other. However, position shifts of the solitons are not captured well, with the left-outgoing soliton shifted too much, and the right-outgoing soliton shifted too little. This particular example represents the configuration that has the best agreement for the velocities. Although not shown here, for δ < 0.3π the particle model predicts that the solitons form a perpetual bound state with a centre of mass coordinate that increases linearly with time. Otherwise, for δ > π/2, the dynamics feature a hard-core elastic interaction where the solitons collide, but rebound off each other. Although the dynamics in these two regimes are similar to what we have obtained numerically, in that we can identify regimes where the interaction is repulsive (and therefore elastic) and attractive (supporting bound states), this correlation actually arises from discrepancies in our model.
To explain how these discrepancies appear in our results, we restate the consequences of the various approximations that we have used in the analysis. In effect, all inconsistencies can be traced back to the initial ansatz used in the analysis, see Eq. (25). The first problem is the obvious fact that the ansatz is a linear superposition of two solitons, neglecting the nonlinear deformation which takes place when they overlap significantly. For this reason, important details of the interaction are omitted.
The second discrepancy in the ansatz, arises from the need to fix and equate the soliton widths. For two chiral solitons to interact, their velocities, and therefore, by extension, their widths, must be different. Furthermore, as the solitons' velocities can change after colliding, b(t) is a time-dependent quantity that requires the additional term w (x − ξ i ) 2 in the expression for the phase given by Eq. (26), with w being the chirp. Although this was derived, it was eventually excluded due the complexity in implementing it in the particle model. The net result is that interaction effects in our particle model are isotropic with respect to each soliton's mutual influence, which is clearly not the case in the full numerics.
Another inconsistency is the divergence of the amplitude Eq. (28) for short length scales in the regime of δ > π/2. This artifact enters due to approximating the phase difference as a function of only the central phases, therefore neglecting velocity contributions. Although this was justified by taking the velocities small, the spatially-varying form of this phase is required to obtain sensible results, as was shown in [54]. Due to the divergence of the amplitude, the effective interaction potential between the solitons also diverges in this regime, as we will demonstrate in the following subsection.
The final discrepancy takes place due to the static nature of the phase in our model. As was demonstrated in our simulations, the current nonlinearity introduces the population transfer and shifts in the soliton's central phase at each collision. The absence of these properties in the model results, therefore, in the existence of perpetual bound states, in addition to the lack of changes in the soliton's amplitudes/widths due to the populations transfer.
From this, it is sensible to conclude that the variational analysis presented here is more suited to studying dynamics at the onset of the collision, before the soliton envelopes significantly overlap, but not in the course of the collision proper.

The interaction potential
The set of coupled differential equations can be reformulated into a mechanical system, in order to derive an effective potential describing the interaction between the solitons. We begin by first restricting our analysis at the onset of the collision, before the solitons began to significantly overlap, as said above. In this regime, the second soliton remains approximately stationary and we can fix that ξ 2 (t) = 0 and v 2 (t) = 0. The set of coupled differential equations then greatly simplify, and we may readily integrate Eq. (37) to obtain v 1ξ1 γ = −4 L I |q n + C, where C is an arbitrary integration constant. Substituting Eq. (32) into Eq. (39), and then substituting again to remove a factor ofξ 1 , leads to the mechanical energy equation, where we identify the soliton kinetic energyξ 2 1 /2, total energy C/2f γ 2 , and the effective interaction potential, The structure of Eq. (40), treats the motion of the first soliton as a classical particle moving through the potential landscape of the second. We plot the interaction curves in Fig. 7, in comparison to an asymptotic calculation derived in the next section. The interaction curves are here plotted only for negative values of the separation, as we have considered the situation in which the moving soliton approaches the stationary soliton up to the start of the collision.
For a typical set of parameters which we have used in our simulations, the curves in Fig. 7 show that both repulsive and attractive interactions are supported, up to a choice of the phase difference. The presence of an attractive potential-in particular, for δ = 0 and δ = π/4,-therefore supports the existence of bound states, and the interpretation of the 'resonance' regime in Fig. 2 as an example of the attractive dynamics. The effect of the divergence of the amplitude for δ > π/2 is illustrated by the upper two curves of Fig. 7, as the magnitude of the potential increases rapidly.
In fact, the nature of these curves does not differ drastically from the results expected in the NLS equation, as contributions from terms stemming from the current feature as a short-range attractive potential, hence they do not have significant influence far from the centre of the interaction. This can be compared to the plot with a 1 = 0, from which it is well known that δ = 0, δ = π, and δ = π/2 correspond, respectively, to the attractive, repulsive, and weak intermediate interactions of solitons in the standard NLS equation.

B. Asymptotic calculations
In the previous section, the interaction of two chiral solitons was explored using a full variational approach. It  is further useful to consider an asymptotic approach for well-separated solitons. Such an analysis can be accomplished using the methodology presented in Refs. [62,63], where tails of the solitons are used to produce the interaction potential (see also [67]). We begin by writing the one-dimensional Hamiltonian density for the system, which, when minimised, reduces to Eq. (7). To derive an effective interaction potential, we again restrict our analysis to the regime of weak-chiral interactions in which the second soliton can be taken as stationary with respect to the first. In addition, since |a 1 (v 1 + v 2 )| |g|, we note that the dominant contributions to the interaction potential in Eq. (41) arise from terms not containing a factor of v 1 . Therefore we can further impose that the first soliton is also stationary such that v 1 = v 2 = 0. Then provided the solitons' centres of masses are well separated by a distance L = ξ 1 −ξ 2 , we may approximate the two-soliton state in the vicinity of the first soliton as in which ψ 1 represents the envelope of the first soliton, and ψ 2 is the exponential tail of the second soliton. The reciprocal approximation is valid in the region around the second soliton.
Next, we substitute Eq. (43) and its counterpart pertaining to the second soliton into Eq. (42), and retain terms which are linear with respect to the small tails, either ψ 1 or ψ 2 . The resulting expression can be recast in a compact form as where each derivative is evaluated at point ψ = ψ 1 .
To proceed, we note that the Hamiltonian can be defined as H = +∞ −∞ dxH, with the variational (alias Fréchet) derivative, As we consider only the case of exact solutions for which δH/δψ = 0, we may recast Eq. (44), using Eq. (45) and integration by parts, to obtain the expression which, after using Eq. (44), simplifies to 2}.
(47) Rather than taking the integration limits in the domain (−∞, ∞), we here divide the integration domain at an arbitrary point z 0 located between the solitons, and introduce the symmetric contribution {1 2} to account for the contribution from the second soliton. To proceed, we require expressions for single-soliton states of Eq. (43), which are given by where the soliton's width b = −4/g is chosen such that n +∞ −∞ dx|ψ n | 2 = 1, with each soliton containing half the number of atoms. In keeping with the linearisation procedure used in deriving Eq. (46), we can instead simplify these expressions in the vicinity of z 0 , with asymptotic forms To calculate the variational derivative in Eq. (45), we use the full expressions in Eq. (48) to evaluate the upper limit at x = ∞, together with the asymptotic forms in Eq. (49) to evaluate the lower limit at z 0 . To obtain a contribution from the current to the effective potential, which is independent of the choice of arbitrary point z 0 , we must go to the next order in the expression for ψ 2 in the second term in expression (47), taking Substituting these expressions, we obtain an effective interaction potential which is convenient to extract information from. The first term in Eq. (50) comes from the NLS equation, being attractive/repulsive for the correct choice of δ. The second term, which appears due to the current, may also be attractive or repulsive, but it is out of phase with the first term. However this term yields a shorter interaction range, compared to the NLS term, hence it does not contribute significantly to the interaction potential far from the centre of the interaction. It is relevant to compare the interaction potentials produced by both the linear-superposition and asymptotic models. Because the asymptotic calculation neglects some terms in the underlying Hamiltonian density, compared to the linear-superposition ansatz (25), it is necessary to approximate the result produced by the linear superposition, using the asymptotic forms, so that they can be fairly compared. This is detailed in Appendix B. Figure 7 shows a comparison between these interaction potentials, from which we can see that both curves share the same qualitative features, with both repulsive and attractive interactions supported in a similar manner to the standard NLS equation.

V. CONCLUSION
We have demonstrated how the current nonlinearity introduces non-integrable effects in the collision dynamics of bright matter-wave solitons. Using the variational approximation, we have derived an effective particle model for the soliton dynamics, which helps to explain both inelastic scattering and the attractive/repulsive nature of the interactions. We have also derived effective potentials for the interaction between the solitons. We showed that the particle model is valid as long as the current nonlinearity is weak, similar to the situation in integrable models, where inherent symmetries of the system can be exploited [54]. This fact implies that essential results may be produced by collisions between slowly moving solitons. For stronger interactions, the particle model breaks down due to the non-integrability. We observe, in particular, how the strong current-induced nonlinearity can destabilize bound states of solitons and also induce soliton fission, breaking two colliding solitons into several ones after the collision.
These concepts constitute a rich spectrum of dynamics which are interesting from the point of view of the fundamental nonlinear dynamics. In addition, the chiral properties of the quantum gas studied here may provide novel applications to atomtronics [68] and quantum transport. In such scenarios, careful consideration of the collision dynamics is needed. by Res (f (z), z 1 ) = lim Collating the results, we find I 1 = 2b sinh( ) . (A6)

Integral example II
A more involved example is In this case, we consider the following contour integral: C g(z)dz = C z 2 sinh(z + ) cosh(z) cosh 2 (z + ) dz, in which the contour path C is identical to the one used in the previous example. In contrast, the integrand now contains a simple pole and a second-order one. Again, using the residue theorem, we can write b C g(z)dz = −2πiI 2 + π 2 I 3 = 2πib k=1,2 Res (g(z), z k ) , where the task of evaluating I 2 now depends on computing the residues of g(z) and knowledge of the integral tanh(α + )dα cosh(α) cosh(α + ) = − 2b sinh( ) + 2b cosh( ) sinh 2 ( ) . (A10) The residues in this instance are given by Res (g(z), z 1 ) = lim z→z1 (z − z 1 ) g (z) = π 2 cosh( ) 4 sinh 2 ( ) , Res (g(z), Collating the results, we find For the remaining interaction integrals, one can proceed using the same methodology, provided that the correct contour shift is applied to the integral, namely, To directly compare the interaction potentials derived from the variational analysis, it is necessary to reduce the potential obtained from the linear-superpositon ansatz (25) to a form which captures the same approximations which were used in the asymptotic version, based on ansatz (43). This can be done by introducing asymptotic forms for each contribution to the interaction potential and selectively dropping terms which are either short-range ones, or contain an explicit dependence on the soliton velocities.

(B3)
In writing Eq. (B3), we have neglected term ∝ e 3 and the term ∝ a 1 d ( , δ) /6 in Eq. (B1), as they, being small in the asymptotic limit, do not essentially affect the result.