{'fraction_non_alphanumeric': 0.04002952532364297, 'fraction_numerical': 0.01436520554167613, 'mean_word_length': 4.758541768840935, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stability of matter-wave solitons in a density-dependent gauge theory

We consider the linear stability of chiral matter-wave solitons described by a density-dependent gauge theory. By studying the associated Bogoliubov-de Gennes equations both numerically and analytically, we find that the stability problem effectively reduces to that of the standard Gross-Pitaevskii equation, proving that the solitons are stable to linear perturbations. In addition, we formulate the stability problem in the framework of the Vakhitov-Kolokolov criterion and provide supplementary numerical simulations which illustrate the absence of instabilities when the soliton is initially perturbed.


I. INTRODUCTION
The stability of solitons in the presence of perturbations represents one of the fundamental problems in the study of solitary waves. Only solitons which are generally robust to perturbations are suitable for study in a physical setting, and by extension, implementing into potential applications in both science and industry [1,2]. Over the years, this topic has been studied extensively for various nonlinear models, with analysis generally falling into two frameworks: the study of small (linear) perturbations of the soliton envelope via a linear stability analysis [3][4][5][6], or the study of additional perturbative terms in the model through a perturbation theory or variational analysis [7][8][9][10][11][12]. The motivation here is twofold: to establish the intrinsic stability of the soliton in a given model, but also to consider the effect of physically relevant perturbations which may influence or potentially damage the soliton. The latter point is of particular importance, as realistic systems are generally described by non-integrable models, in which solitons can potentially be unstable.
In the setting of nonlinear optics, described by the generalized nonlinear Schrödinger equation, understanding the stability of solitons in the presence of perturbations has been a fundamental area of research in the design of soliton-based optical communications [13,14]. Most notable is the Vakhitov-Kolokolov criterion [15,16], which connects the linear stability of bright solitons to two key properties: the number of negative eigenvalues in the spectral problem, and the behaviour of the power integral with respect to the propagation constant. In addition to this pioneering work, several studies have addressed the addition of perturbative terms in the model, such as but not limited to, the excitation of internal modes in non-Kerr media [17][18][19], the effects of third-order dispersion and self-steepening [10,20], and more recently for PTsymmetric potentials [21][22][23][24], which describe media with complex refractive indices.
In ultracold atomic gases, the linear stability framework is more commonly referred to as the Bogoliubov-de Gennes equations and plays a historic role in understanding the superfluid properties of the gas [25][26][27]. Following the first experimental realisation of these ultracold gases, considerable work was centered around studying the response of the trapped condensate to small perturbations [28][29][30][31][32], highlighting the collective nature of the low-lying excitations. The second generation of studies considered the case of both dark and bright solitons [33][34][35][36][37][38], which in turn has lead to a number of key works, such as the interactions between trapped bright solitons [39,40], generation of soliton trains [38,41,42], the reflection properties of bright solitons [43,44], and understanding the snaking instability of dark solitons [45,46].
As the design of these ultracold systems become more involved, it is important to understand how the dynamics and stability of the condensate is modified, as retaining the coherent property of the gas is of vital importance for the purposes of interferometry and future atomtronic based technologies [47][48][49]. In particular, the ability to simulate artificial gauge potentials in these systems [50,51] can lead to the emergence of vortices [52] and spin-orbit coupling [53], offering new means to control and manipulate the gas. Recently, proposals have appeared which center around the engineering of density-dependent gauge potentials which feature a back-action between the matter-field and the gauge potential [54,55]. The condensate dynamics in these models can be greatly modified, including collective modes which violate Kohn's theorem [56][57][58][59], unconventional vortex dynamics [60,61], and the emergence of chiral solitons which feature non-integrable collision dynamics [62]. Very recently, experiments have appeared in which a dynamical gauge theory was realised in trapped ion systems [63] and a density-dependent synthetic gauge field in a Bose-Einstein condensate loaded into a two-dimensional lattice [64].
In this paper, we study the linear stability of onedimensional chiral matter-wave solitons. We begin, by briefly reviewing the physical model, in which a densitydependent gauge potential is optically engineered in a ultracold atomic gas. Then, in Section III, we derive the Bogoliubov-de Gennes equations and study the stability spectrum both numerically and analytically in Sec. IV and Sec. V respectively. Finally, in Sec. VI, we reformulate the stability problem using the Vakhitov We consider the model studied in Ref. [54], in which a harmonically-trapped two-level Bose-Einstein condensate is optically dressed by an external laser field. The Hamiltonian describing the system can be written aŝ with the interaction matrix (2) The light-matter interactions are parametrised by the Rabi frequency Ω and phase φ (r) of the laser field, with the mean-field interactions controlled by the scattering parameters g ii = 4π 2 a ii /m, with a ii corresponding to the scattering lengths for collisions between atoms in state i and i . The harmonic potential appearing in Eq. (1) is chosen to be tightly confined in the radial plane r ⊥ (y, z), but free along the axial x-axis such that the condensate dynamics is effectively one-dimensional.
By treating the mean-field interactions as a small perturbation to the laser coupling Ω g ii |Ψ i | 2 , it can be shown in the dressed state picture [50,54], that the effective detuning induced by the interacting gas can give rise to a density-dependent gauge potential where A (0) is the single-particle vector potential and a 1 controls the strength of the density-dependent gauge potential with ± indices labelling the dressed states. In the following we choose without loss of generality, one of the dressed states and drop the ± indices. By reducing to the one-dimensional picture, the dynamics is described by the Gross-Pitaevskii equation which contains a probability current of the form in addition to the standard cubic nonlinearity. The strengths of the scattering parameters are given by a 1 = k l (g 11 − g 22 ) / 16πl 2 ⊥ Ω and g 1D = (g 11 + g 22 + 2g 12 ) / 8πl 2 ⊥ , which are scaled by the harmonic length l ⊥ = /mω ⊥ and laser phase φ l = k l x. The interacting gauge theory described by Eq. (4) represents a novel nonlinear model, in which the condensate dynamics is influenced by a back-action between the matter-field and the gauge potential.
Rather than working with Eq. (4) directly, we will instead consider the system with which is arrived at by using the nonlinear transformation In the literature, Eq. (6) is often referred to as a 'chiral nonlinear Schrödinger equation', which was originally studied in the context of one-dimensional anyons [65]. Compared to the standard Gross-Pitaevskii equation, this model is generally non-integrable [66,67], and posses chiral soliton solutions which arise due to the breakdown of Galilean invariance [68,69]. As solitons in nonintegrable models can be unstable to perturbations [5,6], we are naturally concerned with the stability of the soliton solutions in our model.

A. Conservation laws
The principle conservation laws underlying the chiral model are given by the integral expressions [68] and which quantify the number of atoms (power integral), momentum, and energy of the condensate respectively. An important dynamical feature of the model is highlighted by Eq. (10), which shows that the gauge field contributes to the momentum of the condensate in a similar manner for a particle travelling in an electromagnetic field. Note, that the Hamiltonian density defined by the integrand of Eq. (11) excludes the probability current, but does correctly reduce to Eq. (6) provided Hamilton's equations are also transformed via Eq. (8).

B. Chiral solitons
In order to derive and subsequently study the stability of the soliton solutions of Eq. (6), it will prove advantageous to work in the moving frame of the condensate as opposed to the stationary frame. To this end, we introduce the Gailiean transformation where the stationary coordinates (x, t) and moving coordinates (x , t ) are related by the translations, x → x−vt and t → t, with frame velocity v. The dynamics of the condensate in the moving frame is then described by the equation of motion The introduction of the renormalised scattering parameterg 1D = g 1D − 2a 1 v, highlights that Eq. (6) is not Galilean invariant, with the strength of the mean-field interactions dependent on both the magnitude and direction that the condensate is moving. For the remainder of this paper, we will explicitly drop the prime notation in the coordinates for brevity and work exclusively in the moving frame unless otherwise stated.
The bright soliton solutions of Eq. (13) then admit the standard form Φ (x, t) = ϕ S (x)e −iµt/ , with envelope width b = −2 2 /mg 1D , and chemical potential (propagation constant) µ = −mg 2 1D /8 2 . In this example, each of these quantities is constrained by normalising ϕ S to unity, providedg 1D < 0. Due to the breakdown of Galilean invariance, both the width and chemical potential of the soliton will depend on the direction of motion. The soliton solution described by Eq. (14) is therefore chiral, such that under appropriate conditions the soliton can either be stable or unstable in a given direction [54]. In recent work, we have shown how the non-integrability of the model can lead to interesting collision dynamics between pairs of chiral solitons, featuring inelastic trajectories and population transfer, in addition to soliton fission and the formation of two-bounce resonance states [62].

III. BOGOLIUBOV-DE GENNES EQUATIONS
To study the linear stability of the chiral soliton, we proceed in the standard way by introducing the condensate wave function [18,26] in which the soliton envelope ϕ S is perturbed by smallamplitude excitations, u (x) and v (x), with frequency ω. Substituting Eq. (15) into Eq. (13) and linearising to first order in u (x) and v (x), leads to the zeroth-order equation and the pair of Bogoliubov-de Gennes equations [56] which feature the standard linear operators for κ = {1, 3}, in addition to the current operator Together, Eqs. (16) and (17) describe the perturbation dynamics of the soliton, with the stability properties determined by the nature of the eigenvalues, or 'stability spectrum', of the Bogoliubov-de Gennes equations. Note, that the linearised operatorL is not self-adjoint (see Eq. (34)), even in the standard case with a 1 = 0.
A key feature of these equations is highlighted by the property that in the moving frame, the current operator does not explicitly couple to the envelope of the soliton, but does couple to the excitations around it. In turn, this leads to the zeroth-order equation for the stationary soliton being described by the integrable Gross-Pitaevskii equation, despite Eq. (6) being generally non-integrable. Therefore, it is expected that the spectrum of excitations around the chiral soliton will be similar to that of integrable models, with only their form modified slightly due to the coupling of the current operator, as per Eq. (19). This proposition will be a key underlying point in the analysis to follow.
We can also conclude several additional properties of the excitations by studying the matrix which is obtained by substituting the soliton solution into Eq. (17). Appearing in Eq. (20) are two potential functions and which are a standard reflectionless potential [70,71] and a gain-loss distribution for the excitations respectively [21][22][23][24]. Together, they form a modified (hyperbolic) Scarf-II potential with the following properties [21]: These properties highlight that the soliton acts as a complex effective potential for the excitations, with bound states and scattering states supported for the attractive potential V (x), and the gains and losses of the excitations balanced by the symmetry of the imaginary potential ∞ −∞ dx W (x) = 0. As this potential is not self-adjoint, the stability spectrum for the excitations could potentially contain complex eigenvalues, in addition to a PTsymmetry breaking point featuring exceptional points [72].

IV. NUMERICAL RESULTS
With the Bogoliubov de-Gennes equations derived and their properties reviewed, we can now proceed in solving for the stability spectrum of the soliton. To achieve this, we first consider a numerical solution, in which Eq. (17) is discretised with periodic boundary conditions and subsequently solved using a sparse eigenvalue solver. The resulting eigenvalues and eigenvectors are shown in Fig As expected, we find that the eigenspectrum retains several characteristic features commonly encountered in integrable models [5]. The eigenvalues consist of the union of two sets: a continuous spectrum with two symmetric branches each gapped from the origin by | ω| = |µ|, and a discrete spectrum with one pair of eigenvalues located at ω = 0 and another pair displaced from the origin by a small imaginary component. At first glance, this pair of imaginary eigenvalues would indicate an instability mode where the soliton state can collapse. However, this component arises instead due to the discreteness of the numerical analysis and subsequently vanishes in the continuum limit (see Appendix A). The eigenvalues of Eq. (17) are therefore entirely real with a four-fold degenerate eigenvalue at ω = 0, in an identical manner as for the Gross-Pitaevskii equation. In reference to the PT -symmetry breaking proposed earlier, we find that the eigenvalues are consistently real in the parameter spacẽ g 1D < 0, indicating the absence of a PT -symmetry breaking point for the regime in which the soliton solution is valid. Therefore we may conclude that in free space, the chiral soliton is stable to linear perturbations.
The eigenvectors of Eq. (17) are also consistent with that of integrable models, with the discrete spectrum corresponding to localised real-valued solutions in the vicinity of the soliton (Fig. 2.(a-d)), while the continuous states are complex-valued and generally oscillatory at x = ±∞ (Fig. 2.(e-f)). In fact, the discrete states pictured are exactly the same as for the Bogoliubov-de Gennes equations associated with the Gross-Pitaevskii equation, with the first and second states taking the form of the envelope of the soliton and its derivative respectively. This similarity, or rather, the invariance of the form of the discrete states in the presence of the current  operator, will become clear from the analytical results in the next section.
In obtaining these numerical results, the proposition that the stability spectrum of the chiral soliton is similar to that of integrable models is validated. As such, the interacting gauge theory represents a non-integrable model in which the soliton solutions are stable to linear perturbations. The absence of instability modes in the model will be revisited in the numerical simulations presented in Sec. VI A.

V. ANALYTICAL RESULTS
For several integrable models, the spectrum of excitations around the soliton solution can be derived analytically [4,73,74]. This is generally achieved either through a connection to the squared eigenfunctions of the associated eigenvalue problem [75], or in some cases, by direct methods. As our model is generally non-integrable, it is unclear whether the former method would be applicable. However, using the numerical results obtained previously as a basis, we can at the very least obtain expressions for the discrete spectrum using traditional methods.

A. Discrete Spectrum -Bound States
As was demonstrated previously, the discrete spectrum of the excitations correspond to a set of bound states which are a zero-eigenvalue solution of Eq. (17). For this reason, the Bogoliubov-de Gennes equations can be written as and The above set of differential equations can be solved analytically using either a hypergeometric series approach [76][77][78], or by operator methods [79] (see Appendix B for details on the former). Then, by denoting ζ n = u v T for n ∈ Z + , we find the pair of zero-eigenvalue solutions and where χ = x/b, and c 1 and c 4 are arbitrary constants to be determined. Due to the four-fold degeneracy of the zero-eigenvalue, two linearly-independent solutions, ζ 3 = 2 (ζ 3,− + ζ 3,+ ) and ζ 4 = 2 (ζ 4,− + ζ 4,+ ), are also required for completeness, where and are a set of generalized eigenvectors which satisfy the eigenrelations,Lζ 3,± = ±2ζ 2 andLζ 4,± = ±2ζ 1 , up to a constant [7,74]. It is straightforward to verify that the solutions form a linearly independent set by calculating the Wronskian and likewise for the remaining solutions. Together, these states compose the discrete spectrum of Eq. (17), corresponding each to small variations of the soliton with respect to its four free parameters: phase, position, velocity, and chemical potential. As an example, by expanding the soliton solution around a small variation of the velocity δv, such that b → b + δb, one finds to first-order that sech χ e iδvx = sech χ which is captured by ζ 3 with reference to Eq. (15). Surprisingly, this set of discrete states is identical to that of the Bogoliubov-de Gennes equations for the Gross-Pitaevskii equation, except for ζ 3 which features currentdependent terms due to the width of the soliton being defined in part by its velocity. This similarity becomes clear when obtaining the solutions, as one finds that the soliton envelope is a zero eigenvalue solution of the operator L 1 ϕ S = 0 and the linearised current operator J 0 ϕ S = 0. Therefore, the eigenvalue problem for both ζ 1 and ζ 2 effectively reduces to the standard case, which naturally leads to the form of the solutions pictured in Fig. 2.
Defining the inner product between any pair of solutions as together with the normalisation constraint defined by Eq. (9), leads to the values of |c 1 | 2 = 1/ (2b) and |c 4 | 2 = 3/ (2b) for the arbitrary constants. Note, that since ζ 3 and ζ 4 are not square-integrable solutions, we are not required to compute inner-products for the remaining constants as these states are not physical.
With the form of the solutions now reduced, we are now in a position to compare our analytical results to the numerical ones obtained earlier. In Fig. 2, we find exact agreement between the analytical (dots) and numerical (solid-line) results for the degenerate bound states ζ 1 and ζ 2 . This reinforces the statement that the discrete eigenvalues in the chiral model are four-fold degenerate at ω = 0, with the imaginary component in the numerical analysis attributed to numerical artifacts. Therefore we may conclude, that due to the consistent results obtained from both methods, that the chiral soliton is stable to linear perturbations.
In addition to the standard eigenvalue problem, we can also solve the corresponding adjoint problem using the same analytical techniques. Although this system does not have any physical relevance, the equivalence of both eigenvalue problems will be a key property which we will exploit when deriving the Vakihitov-Kolokolov criterion in the next section.
For the sub-space of square-integrable solutions ofL with the inner-product defined by Eq. (33), we write the whose right-eigenvectors are the adjoint of the lefteigenvectors ofL. The adjoint of the linearised currentoperator appearing in Eq. (34) is given by Together, Eqs. (34) and (35) explicitly highlight thatL, as mentioned previously, is not self-adjoint.

C. Continuous States
For the continuous states, we are unfortunately unable to derive a closed-form expression due to the complexity of the eigenvalue problem. Although we may be able to obtain these solutions using a power series method [80,81], the complicated nature of the calculation presents little benefit for the knowledge gained, since the continuous states are irrelevant for addressing the stability of the soliton. However, we can still obtain an expression for the eigenvalues using an asymptotic approach as follows.
The key point to note, is that for large distances away from the centre of the soliton, the continuous states are generally oscillatory such that they can then be written as plane waves of the form with wave-number q. Then, by substituting Eq. (42) into the asymptotically reduced form of Eq. (17), one finds the continuous-eigenvalue expression which as expected, is simply a free-particle dispersion relation gapped by the chemical potential of the soliton. In Fig. 1(a), we compare both the numerical (red solidline) and analytic (red dots) values for the continuousstate band edge (q = 0) which as shown, is in exact agreement.

VI. VAKHITOV-KOLOKOLOV CRITERION
As the chiral solitons present in our model are solutions to a generalized Gross-Pitaevskii equation with a real positive envelope, we can also establish their stability properties using the Vakhitov-Kolokolov criterion [5,15]. For the standard Gross-Pitaevskii, the stability criterion requires: (i) The eigenspectrum of the operators L 1 and L 3 , should contain at most only a single negative eigenvalue.
The first condition is straightforward to verify by noting that both operators have a positive continuous spectrum ω = [−µ, ∞), and a discrete spectrum defined by the eigenrelations, and The second condition can also be proven easily by direct integration, provided the particle number (power) is correctly posed [3]. This criterion directly follows from the definiteness of the operators in the eigenvalue problem and therefore requires us to extend the Vakhitov-Kolokolov analysis to include the linearised current operator as follows.
As the eigenvalues ofL andL † are conjugate to each other, we can without loss of generality consider the spectral properties of either system. To this end, it will prove advantageous to work in the adjoint picture, as the following stability analysis is simpler while exploiting the property that the operators L 1 and L 3 are self-adjoint.
We start, by restating the adjoint eigenvalue problem and Taking the inner-product of Eq. (46) with ϕ S (x) leads to the expression As L 1 is self-adjoint with L 1 ϕ S = 0, the above expression is true ∀ω * , provided the orthogonality condition is satisfied. Therefore, for the non-zero eigenvalues in the stability analysis, we may restrict ourselves to the function space where the inverse operators L −1 1 and L −1 3 are definable. Returning to Eqs. (46) and (47), we can now proceed in constructing the stability criterion for the chiral soliton by combining both equations into the fourth-order equation in the function space S. Multiplying Eq. (51) by L −1 1 , taking the inner product with respect to ν(x), and completing the square leads to the expression with α = ν, L −1 1 ν , β = ν, L −1 1 J † 0 ν , and γ = ν, L 3 ν . The condition of stability is now set by requiring that the right-hand side of Eq. (52) be non-negative, such that ω * by extension is real. Otherwise for negative values, ω * would be imaginary, thereby indicating an instability. As both L 1 and L 3 (and their inverses) are known to be positive definite in the space S [75], the standard term γ/α will not need to be considered in our analysis.
Instead, the stability of the soliton will be resolved by studying the definiteness of L −1 1 J † 0 . Despite being a trivial reduction of the problem, the non-negativeness of Eq. (52) can be guaranteed if J † 0 is nilpotent in the domain considered, i.e. J † 0 ν = 0, such that all the eigenvalues of J † 0 are zero. This can be proven by directly solving the eigenvalue problem with eigenvalue λ and eigenfunction ν. The general solution of Eq. (53) can be readily found, where = mb 2 /a 1 , and C is an arbitrary constant. For a ring domain of length (circumference) L with periodic boundary conditions ν (−L/2) = ν (L/2), the eigenvalues form a continuous spectrum with σ = 0, ±1, ±2, . . . . Then in the combined limit where L → ∞ and σ → ±∞, the eigenvalues of Eq. (53) coalesce at λ = 0, highlighting that for the ring domain considered in our numerics, the adjoint current operator J † 0 is nilpotent. Therefore, the stability condition, Eq. (52), reduces to the form encountered for the standard Gross-Pitaevskii equation, from which the Vakhitov-Kolokolov criterion is known to be satisfied [5,15].

A. Numerics
In addition to the results obtained analytically, we also consider a set of numerical simulations which illustrate the stability of the chiral soliton under the influence of a perturbation. To achieve this, we follow the standard numerical scheme in which the initial number of atoms (power) of the soliton differs from the exact solution and observe whether the soliton collapses or retains its shape [5,6]. As such, we define the perturbed soliton state as and show two examples of the pertubation dynamics in Fig. 3, each for a different sign of the perturbation parameter ∆ϕ. In both cases, the soliton maintains its shape over the course of the simulation and does not collapse, disperse, or oscillate due to the excitation of an internal mode [6,19]. Instead, the soliton emits a small (nonvisible) amount of radiation and decays to the stable lowamplitude solution, in a similar manner to solitons of the Gross-Pitaevskii equation. However, as the initial width of the soliton changes due to the perturbation, with g 1D and a 1 fixed, the velocity of the perturbed soliton will differ from the frame velocity set by Eq. (12). This results in the soliton drifting in the moving frame, with the direction controlled by the sign of the perturbation. This effect is not captured by the stability spectrum of the Bogoliubov-de Gennes, due to it being a higher-order effect which is neglected in the linearisation. Instead, we are required to consider an alternative framework to explain the presence of the soliton drift.

B. Variational equations
To quantify the drift of the soliton due to the action of the perturbation, a variational approach can be employed to derive a set of coupled equations which effectively describe the soliton dynamics [82]. Although we will not be able to solve these equations analytically, their numerical solutions will provide sufficient illustrations of the perturbation dynamics to compare with the results pictured in Fig. 3.
We begin by writing the Lagrangian density [69], which is written in the stationary frame of the soliton, with stationary coordinates (x, t). To accurately describe the perturbation dynamics, we choose a general variational ansatz of the form [71,83] ψ (x, t) ≡ a sech ( in which the width and position of the soliton envelope can vary dynamically through the spatially-varying phase Here, a(t), b(t), ξ(t), k(t), w(t), and φ(t) are timedependent variational parameters corresponding to the amplitude, width, centre-of-mass coordinate, velocity, curvature, and absolute phase of the soliton. As the form of Eq. (58) explicitly assumes that the shape and particle number of the soliton is conserved, the interplay of radiation will therefore be absent in the analysis. However, as we will demonstrate, provided the magnitude of the perturbation is kept small, this discrepancy will not have significant implications. Substituting Eq. (58) into Eq. (57) and minimizing the corresponding action functional leads to the set of coupled differential equations and which collectively describe the motion of the soliton. The source of the drift is now clear from the coupling between Eqs. (60) and (61); that a time-dependent variation of the soliton's width, induced by a perturbation, can lead to a change in the soliton's centre-of-mass proportional to the strength of the gauge field. The trajectory of the soliton is then set by Eq. (62), which in a consistent manner to Eq. (10), contains an additional contribution from the gauge field. Therefore, for either an increase or decrease in the particle number, it is expected that the soliton will drift in the moving frame. To illustrate the above reasoning, we solve the set of differential equations numerically using a fourth-order Runge-Kutta method and plot the predicted soliton trajectories (red solid-line) in Fig. 3. Both the direction and magnitude of the drift is captured correctly by the variational equations and therefore validates that the drift of the soliton arises due to how the initial state is prepared. In addition, these results show that in this weak perturbation regime, the emission of radiation from the soliton plays no significant role in the dynamics. However, its absence in the variational description does lead to inconsistencies, as demonstrated by the presence of small-amplitude oscillations in the predicted trajectories which persist indefinitely.
To conclude, although we cannot strictly say the soliton is stable due to the presence of the drift, we stress that it is a manageable feature which does not destroy or damage the envelope of the soliton. Therefore, we may view the soliton as effectively stable, with the absence of the traditional instability mechanisms consistent with the Bogoliubov-de Gennes analysis.

VII. CONCLUSION
In this paper, we have demonstrated the linear stability of chiral matter-wave solitons in an interacting gauge theory. Despite being described by a non-integrable model, we found that the stability spectrum of the soliton reduces to the standard integrable case, with entirely real eigenvalues and the absence of instability modes. This was then further understood by studying the Vakhitov-Kolokolov criterion, which highlighted that the linearised current operator was nilpotent in the numerical domain and therefore does not contribute to dynamical instabilities.
The drift of the soliton due to the presence of a perturbation represents an interesting property of the chiral model. By generalising the study to a broader class of perturbations [84], several questions are inspired not only in regards to the stability, but also to features which could be exploited in order to control the soliton. For example, could a perturbation be designed which when applied continuously, enables the soliton to accelerate or decelerate with minimal radiation losses? These questions, together with the linear stability properties concluded in this work, offers a promising candidate for practical transport dynamics in atomtronic systems [49,85], where retaining the coherent properties of the gas can be an important factor.
In Sec. IV, we claim that the imaginary component of the discrete eigenvalues is a numerical artifact which vanishes in the continuum limit. This feature is a common occurrence in the study of spectral stability, and arises from the numerical model being ill-conditioned; that the soliton, which is strictly speaking a solution in free space, is discretised and truncated in the numerical picture.
To resolve this discrepancy, we define the numerical domain of the soliton as [−L/2, L/2], with length L and spacing ∆x, provided L > b. Then, in the continuum limit, where L → ∞ and ∆x → 0, it is expected that the numerical eigenvalue problem will become wellconditioned and match the analytical results. To demonstrate this, we numerically solve Eq. (17) in each limit independently for fixedg 1D and plot the behaviour of the eigenvalues in Fig. 4(a) and Fig. 4(b) respectively.
In both cases, we find that the eigenspectrum converges to the exact values, with the q = 0 continuous state approaching | ω| = |µ| (as per Eq. (43)) when the domain length is increased, and both pairs of discrete eigenvalues converging linearly to ω = 0 when the domain spacing decreases. Note, that the discrete eigenvalues do not need to be considered in the former case as they are localised to the width of the soliton and will therefore be invariant to variations in L, provided L > b by an order of magnitude. For these reasons, we may conclude that for the fixed chemical potential in this example, a reasonably accurate solution for the eigenspectrum can be obtained with a modest domain length of L = 100 and spacing ∆x = 10 −3 , corresponding each to a error of ≈ 10 −3 .

Appendix B: Analytical methods
In this appendix, we show how to construct the zeroeigenvalue solutions of the Bogoliubov-de Gennes equations. As this method uses standard techniques which are well documented in the literature [76][77][78], we present the following calculations solely for the sake of completeness.
We return to the zero-eigenvalue problem described by Eqs. (23) and (24), which contain a homogeneous and an inhomogeneous eigenvalue problem for the eigenvector components u(x) and v(x) respectively. The general solution will therefore be composed of a complementary solution ζ c for the reduced homogeneous system and a particular solution ζ p to be solved for successively.

Complementary solution
By introducing the soliton width as a scaling parameter, we write the homogeneous system as with χ = x/b, ∈ Z + 0 , and dimensionless eigenvalue E. The potential function appearing in Eq. (B2) is commonly referred to as a modified Pöschl-Teller potential, and has been studied in the context of reflectionless scattering [70,71] and supersymmetry [79,86]. As this potential is attractive, and converges to zero when x → ±∞, the corresponding eigenspectrum will consist of two sets: a bound-state spectrum for E < 0, and scattering states with E > 0. For our purposes, we will only be concerned with the bound-state spectrum.

(B8)
The pair of asymptotic forms described above converge provided the ratio of Γ-functions vanishes. Therefore we require β − γ + 1 = −n and 1 − α = −n, for n ∈ Z + 0 provided E ≤ 0. The resulting expression for the eigenvalues can then be determined iteratively, and takes the form with the constraint 0 ≤ n ≤ − 1.

Particular solution
With the complementary solution derived, we may now proceed in solving for the particular solution of Eq. (24) using the method of variation of parameters. Labelling the pair of fundamental solutions of the homogeneous problem as v 1 = sech χ 1 − χ tanh χ − cosh 2 χ/3 , and v 2 = tanh χ sech χ, we write the Wronskian relation The particular solution can then be obtained by direct integration with h(x) = −2J 0 u c corresponding to the inhomogeneous part of Eq. (24). As the expressions derived from Eq. (B19) are often quite cumbersome, we refer the reader to the solutions presented in Sec. V.