Dark-soliton-like excitations in the Yang-Gaudin gas of attractively interacting fermions

Yrast states are the lowest energy states at given non-zero momentum and provide a natural extension of the concept of dark solitons to strongly-interacting one-dimensional quantum gases. Here we study the yrast states of the balanced spin-$\frac{1}{2}$ Fermi gas with attractive delta-function interactions in one dimension with the exactly solvable Yang-Gaudin model. The corresponding Bethe-ansatz equations are solved for finite particle number and in the thermodynamic limit. Properties corresponding to the soliton-like nature of the yrast excitations are calculated including the missing particle number, phase step, and inertial and physical masses. The inertial to physical mass ratio, which is related to the frequency of oscillations in a trapped gas, is found to be unity in the limits of strong and weak attraction and falls to $\approx 0.78$ in the crossover regime. This result is contrasted by one-dimensional mean field theory, which predicts a divergent mass ratio in the weakly attractive limit. By means of an exact mapping our results also predict the existence and properties of dark-soliton-like excitations in the super Tonks-Girardeau gas. The prospects for experimental observations are briefly discussed.


Introduction
Solitons are ubiquitous nonlinear wave phenomena that appear in many different fields of physics [1]. They often appear in the context of fluid dynamics or the nonlinear mean-field theories of quantum systems [2]. A very interesting open question is to what degree solitary wave phenomena can survive in low-dimensional and strongly correlated quantum systems, where mean-field theories or fluid-dynamic-like descriptions are inappropriate due to the dominance of quantum effects.
An interesting feature of dark solitons is that they exhibit oscillatory motion around the position of maximum density in a trapped gas. The frequency of oscillation Ω is universally related to the ratio of the inertial mass m I and the physical mass m P under the condition of small oscillations [16,17,24]  was predicted from the Gross-Pitaevskii equation in the nonlinear Thomas-Fermi limit [24,25], and measured experimentally in [6,26]. In the BEC to BCS crossover in a Fermi superfluid, the ratio was predicted to increase from 2 in the BEC regime to 3 in the unitary regime and grow even larger towards the BCS regime in [16,17]. Unfortunately, the dark solitons observed in experiments decayed before their oscillation frequency could be measured by means of the snaking instability due to the three-dimensional (3D) nature of the trapped gas [20,27,28]. Instead, a solitonic vortex [29][30][31] was observed [32,33], that owes its large mass ratio m I /m P (in the hundreds) to hydrodynamic effects [33,34].
A natural way to stabilize the dark soliton is to reduce the dimensionality by confining the gas to a tight waveguide-like trap [35]. In order to suppress the snaking instability for the crossover Fermi superfluid, a transverse confinement length scale comparable to the mean particle spacing is required [27,28], which means entering a 1D regime. An analytically solvable 1D mean-field model for moving dark solitons in the BCS regime predicts that the mass ratio m I /m P diverges to +¥ in the BCS limit of weak attractive interactions [36] for a quasi onedimensional channel (whereas  -¥ m m I P for a strict 1D version of their model). These results are consistent with the predictions from 3D mean-field theory [16,17], where m I /m P increases towards the BCS regime. As we will show in this work, however, exact solutions of the Yang-Gaudin model predict that m I /m P approaches unity in the weakly coupled BCS limit for a 1D Fermi gas. The validity of mean-field theory is generally questionable in strongly correlated regimes where no small parameter exists 1 , and in the strongly confined 1D regime where phase fluctuations prevent true long-range superfluid order according to the Mermin-Wagner Hohenberg theorem [37,38]. Nevertheless, the use of BCS-like mean-field theory in the weakly attractive 1D Fermi gas had previously been supported by the comparison of bulk properties with exact solutions of the Yang-Gaudin model [39].
Here we approach the tightly confined 1D Fermi gas from the exactly solvable Yang-Gaudin model [40,41] of spin- 1 2 fermions with attractive δ-function interactions. Using a purely 1D model with contact interactions is well justified for an ultra-cold Fermi gas in the low density regime, where the relevant energy scale, e.g. the Fermi energy of the 1D gas, is small compared to the spacing of transverse excitation levels of the wave-guide trap [42][43][44][45]. While the two-particle problem in a wave-guide trap always has a bound state [46], the many-body problem of spin-1 2 fermions knows two distinct regimes that are separated by a confinement-induced resonance, where the 1D scattering length a 1D passes through zero. The regime of < a 0 1D is well described by the Yang-Gaudin model of spin- 1 2 fermions with attractive contact interactions while positive scattering lengths > a 0 1D lead to a gas of repulsively interacting bosonic dimers described by the Lieb-Liniger model [47,48]. The situation is sketched in figure 1 with reference to the relevant dimensionless coupling parameter g = a n 2 1D 0 of the Yang-Gaudin model. Here, = n N L 0 is the particle number density, N the total number of fermions and L the length of the system.
The Yang-Gaudin model with attractive interactions (g < 0) has two-particle bound states, dimers, with length scale~| | a 1D . The coupling constant g = a n 2 1D 0 thus provides the ratio of the mean particle spacingn 0 1 to the dimer size and g | | 1  is a regime of Cooper-pair-like large dimers. An interesting crossover regime appears at g » | | 1, where the length scales are comparable, while the dimers are tightly bound for g  -¥. This is the position of the confinement-induced resonance, where the dimers are tightly bound but also the dimer-dimer interactions diverge [45]. The many-body physics becomes that of a gas of impenetrable bosons, known as a Tonks-Girardeau gas [49,50]. In this regime, the energy spectrum becomes that of N 2 spinless fermions with mass m 2 , where m is the mass of the original fermionic atoms. In fact, there is a one-to-one mapping between the impenetrable Bose gas (here referring to dimers) and non-interacting spinless fermions (of mass m 2 ) [50]. In this sense both limits of the attractive Yang-Gaudin model are non-interacting fermions and the interesting correlated physics is expected around g » | | 1. There is a close relation between the attractively interacting Fermi gas with a balanced population of spin-up and -down fermions and the repulsively interacting Bose gas. The 1D Bose gas with repulsive δ-function interactions is known as the Lieb-Liniger model and its ground state energy and excitation spectrum were solved Fermi gas with attractive interactions in a one-dimensional wave guide. The physical regimes of the many-body physics are determined by the dimensionless coupling parameter g = = =  a n gm n c n 2 1D 0 0 2 0 .
for in [47,48]. An important finding was a gapless branch of excitations with lower energy than the gapless branch of phonon excitations (similar to the Bogoliubov phonons of the 3D Bose gas). Lieb called these type-II elementary excitations [48]. They are part of the set of yrast states, which refers to the lowest energy state at given momentum 2 . By now there is mounting evidence for the association of the yrast excitations of the Lieb-Liniger model to dark solitons known to exist in the 3D BEC and found mathematically in the nonlinear Schrödinger equation [51][52][53][54][55][56][57][58][59][60]. References [52,53] first identified the energy-momentum dispersion relation of dark solitons with Lieb's type-II branch in the weak coupling limit, which was later extended to multi-soliton solutions [51,54,55]. Reference [56] further strengthened the connection by establishing how the localized dark solitons of mean-field theories (and experimental observations) emerge from a superposition of the translationally invariant yrast states. Syrwid and Sacha went on to demonstrate that localized density depressions typical of dark solitons emerge at random positions from the translationally invariant yrast states under a single-shot measurement procedure, where the position of the atoms is detected one-by-one [57].
The association between dark solitons and yrast excitations in the Lieb-Liniger model opens the door to a deeper investigation of solitonic properties of yrast states in the more strongly interacting regimes of 1D quantum gases, outside the validity of mean-field theory. Astrakharchik and Pitaevskii have studied the (type-II) dispersion relation of yrast excitations in the Lieb-Liniger model throughout the crossover from a weakly interacting Bose gas to the strongly correlated Tonks-Girardeau regime and deduced a number of quantities of interest [61]. These include the inertial and physical masses. One important finding is that the effective particle number m P /m diverges, and thus becomes macroscopic, in the limit of small coupling constant γ, while corresponding to a single boson hole in the Tonks-Girardeau regime of large γ. The effective particle number < m m 0 P (called N s in [61]) is closely related to the number of atoms < N 0 d in the dark soliton, i.e. the number of atoms that have to be removed from the background density to create an observed density depression. Another result of [61] is that the inertial-to-physical mass ratio m m I P , which relates to the oscillation frequency of the soliton in a harmonically trapped gas, smoothly decreases from the value of 2 in the weak coupling limit of the 1D Bose gas (where it agrees with the mean-field prediction), to 1 in the Tonks-Girardeau regime.
In this work we extend the study of yrast excitations to the spin-1 2 Fermi gas with attractive interactions within the Yang-Gaudin model. We identify yrast excitations in the finite system as dimer-hole excitations in the strongly attractive regime and clarify their continuous connection to single fermion holes in the non-interacting limit by solving the Bethe-ansatz equations directly and with the aid of a string hypothesis for fermion pairs in the strongly attractive regime. We then solve the integral equations for dimer holes, which govern the yrast dispersion relation in the thermodynamic limit. From these, a number of physical quantities including the missing particle number, the phase step, and the inertial and physical masses are calculated. We find that the missing particle number varies from one to two missing fermions as we tune from the Tonks-Girardeau gas of dimers to the BCS limit, while both the physical and the inertial masses correspondingly change fromm 2 to -m. Interestingly, the mass ratio m I /m P is 1 in both limits and reaches a minimum value of about 0.78 in the crossover regime g » -1.3 (see figure 9). This contrasts sharply with the results of 1D mean field theory [36], where the mass ratio diverges in the BCS limit.
There is a close relation between the attractive Fermi gas and the super Tonks-Girardeau gas, a highly excited gaseous phase of attractive bosons in 1D [62][63][64][65][66]. Even though attractive bosons form cluster-like bound states or bright solitons in the ground state, a metastable gaseous phase can be supported by kinetic energy pressure in a similar way that Fermi pressure stabilizes the attractive Fermi gas. Due to a one-to-one correspondence of the Bethe-ansatz solutions of the Yang-Gaudin model to the Lieb-Liniger model established by Chen et al [66], our results carry over to the super Tonks-Girardeau gas where they predict the existence of dark soliton-like features, which is quite remarkable. Due to the stability properties of yrast states we may expect the life time of these states to be comparable to that of the 'ground-state' super Tonks-Girardeau gas.
This paper is organized as follows: section 2 introduces the Yang-Gaudin model and the corresponding Bethe ansatz equations for finite particle number. The structure of the ground and yrast exctited states, as well as the methods for solving the equations, are discussed in both the strongly and weakly interacting regimes, and the excitation spectra are computed. Section 3 then deals with the thermodynamic limit dispersion relations. After discussing the ground state thermodynamic-limit equations in section 3.1, the equations for the yrast dispersion relations are given in section 3.2, the missing particle number and phase step are calculated in section 3.3 and the inertial and physical masses in section 3.4. Conclusions are drawn in section 4 and three appendices provide technical details regarding the numerical solution of the finite-system Bethe ansatz equations (appendix A), an outline of the derivation of the excited-state thermodynamic limit integral equations (appendix B), and details pertaining to the numerical computation of the missing particle number (appendix C), respectively.

Yang-Gaudin model for finite fermion number
The Yang-Gaudin model describes a gas of spin-1/2 fermions, confined to a 1D box with periodic boundary conditions and interacting via a two-body δ-function potential. The Hamiltonian takes the form where in the second term the sum is over all pairs counted once. There are N fermions in total,  M N 2 of which are spin-up and the rest are spin-down. Furthermore, m is the mass of each particle, and L is the length of the box. The model is relevant to ultra-cold fermionic atoms in two hyperfine states confined to a 1D wave guide in the low density limit [42][43][44][45]. The 1D interaction constant g can be related to the 3D scattering length a of the atoms, the wave guide trap frequency w^and length scale [49]. The coupling constant diverges and changes sign at a confinement-induced resonance where = a a A 1 . We are here considering the case of attractive interaction where < g 0. It is convenient to write the interaction constant as = where = c a 2 1D and a 1D is the 1D scattering length [67]. Let us also introduce the 1D density, = n N L 0 , and g = c n 0 which is a useful dimensionless parameter both in the finite system and in the thermodynamic limit.

Bethe ansatz equations in exponential form
The Bethe-ansatz solution of the Yang-Gaudin model (2) consists of superpositions of plane waves for the many-body wave function [68]. These are subject to boundary conditions where particles interact and at the box boundaries, as well as fermionic symmetry constraints. The solutions are uniquely determined by a set of variables with the dimension of wave numbers known as rapidities. They have to satisfy the Bethe ansatz equations in exponential form [40,41]: The charge rapidities k j can be thought of as the quasi-momenta of the fermions. They completely determine the total momentum and energy of the system The spin rapidities a m are auxiliary variables and are present due to the spin degree of freedom. The a m ʼs do not contribute to the energy or momentum but must be solved for as they are coupled to the k j ʼs. There are infinitely many different sets of rapidities that solve (3) and (4) and each one corresponds to an eigenstate of the Hamiltonian (2). In this work we are interested in the yrast states, i.e. the states with the lowest energy E tot at given momentum P tot . We also restrict ourselves to balanced populations of spin-up and -down particles, i.e.
= N M 2 . The rapidities for yrast states can be easily identified in the weak and strong interaction limits, where simple analytic solutions to (3) and (4) are known. The yrast solutions for finite interaction strength can then be found by continuity. Examples of rapidities for yrast states are shown in figure 2.
A particular feature due to the periodic boundary conditions is that a new set of rapidities solving (3) and (4) can be generated from an existing one by adding p L 2 (or an integer mutiple) to every k j and a m , as is easily seen from the equations. This changes the momentum to p ¢ = +  P P n 2 tot tot 0 , where = n N L 0 and the energy to . Physically this corresponds to a Galilean boost of the whole system by the umklapp momentum p º n p 2 4 0 F , while the internal structure of the many-body state is unchanged [47]. Here, is the Fermi momentum of the non-interacting spin-1 2 Fermi gas. It is thus sufficient to consider yrast states in the interval <  P p 0 4 F . The Bethe ansatz equations (3) and (4) can be solved directly by numerical methods. However it is only practical to do so for small particle number and in a regime of very weak interactions. The situation is similar to that in the Lieb-Liniger model with attractive interactions, where bound states are formed as well [69]. Technical details regarding the numerical solutions are given in appendix A.

Weak interactions
In the case of weak interactions g | | 1  , an analytical result is available for the rapidities of the ground state [70]. Define the set of integers q m as 2 . An illustration of the rapidities for the case of a non-degenerate ground state with odd M is shown in figure 2 (left, empty symbols). For < c 0 pairs of k j ʼs are complex conjugates and there is a limited range where the exponential equations can be solved in practice.
The energy for weak interactions can be obtained by explicit evaluation of the sum (5) [71]. In terms of the dimensionless coupling parameter g = c n 0 we obtain (assuming odd M) This approximation is compared to the exact results in the bottom panel of figure 3. For c=0 this becomes simply the expression for the energy of a gas of non-interacting spin-1 2 fermions. Yrast excitations in the small | | c limit are obtained by creating a single-fermion hole, which has further implications. When a single k j is moved from To summarize, the real parts of the k j ʼs are determined by the free fermion limit, the splitting in the imaginary axis of both the k j ʼs and a m ʼs is the same as for the ground state, and the real parts of the a m ʼs are always an average of their corresponding k j pairs, again as in the ground state. In the strongly attractive case (right), this excitation is obtained by creating a hole in the 'Fermi-sphere' of dimers, i.e. moving the set of two charge and one spin rapidities from quasi-momentum p -L to the outside of the sphere, at p L 3 . In the weekly interacting regime (left), the corresponding yrast excitation is obtained by taking a single fermion at p -L 2 out of the Fermi sphere. Note that this state is a nondegenerate singlet excitation while the related triplet excitations have higher energy for nonzero c.

String hypothesis
For stronger attractive interaction, as the value of c becomes increasingly negative, the real part of pairs of k j values becomes equal to the value of a m while the imaginary parts become exponentially close to i c 2 . This causes divergent terms in (3) and (4) through vanishing denominators and provides a complication for numerical root finding. In this case, approximate equations can be derived using the so-called string-hypothesis, where the imaginary parts are replaced by the limiting values [70]. Specifically for the ground state, M 2 of the k j ʼs are taken to be a  i m c 2 , and the remaining -N M 2 real k j ʼs are left as unknown variables. The paired charge rapidities can be eliminated from the equations by substitution, with all remaining variables (the a m ʼs and the unpaired k j ʼs) being real and distinct. It is then possible to take the logarithm of the Bethe-ansatz equations. Assuming that all fermions are paired, which is adequate for the balanced ground state and spin-conserving yrast excitations, this leads to [68]: is the two-body phase-shift function of the δ-function potential. The ℓ m ʼs are quantum numbers that specify the state. These are either integers or half-integers, depending on the parity of N and M, and they must be distinct.
Since we know the relation between the a m ʼs and the k j ʼs associated with them, it is easy to write the momentum and energy as . The quantum numbers for the balanced ground state are with the term in brackets only included for even M. In order to solve the logarithmic string hypothesis equations it is easier to start from the  -¥ c limit where it suffices to guess a = pℓ m L m in order to converge to the solution.
The structure of the ground state rapidities in the strongly attractive regime (  -¥ c ) is illustrated in figure 2 on the right by empty symbols. Indeed, it can be seen from the logarithmic equation (9) that the a m ʼs approach p ℓ L m where the ℓ m ʼs are given by (13) (and the corresponding charge rapidities approach . The structure of the ground state is now that of a Tonks-Girardeau gas of dimers: the distribution of the a m rapidities is that of a Fermi sphere with M spinless fermions, each corresponding to a bound pair of elementary fermions. The ground state energy can be obtained from (12) using the expressions for the rapidities in the case of strong interactions. We obtain the energy for a Tonks-Girardeau gas of N 2 bosonic dimers of mass m 2 and binding energy - c m 4 2 2 , as This expression is plotted in the top panel of figure 3 which shows the ground state energy as a function of coupling strength across the range g -< < 100 0 for a system with = = N M 14, 7. In particular, the bottom panel allows us to compare the performance of the string hypothesis equations to the exact exponential equations in the region where < c 0 and g » | | 1 or less. We see that the energy of the ground state is captured very well indeed, which validates the use of the string hypothesis in cases when all particles of opposite spins are paired up into dimers.
In order to identify yrast states in the strongly attractive regime it is important to note that breaking a dimer would cost a large amount of energy ( c m 4 2 2 ). Thus the yrast states of lowest energy for given momentum will be the ones where all fermions are paired, and the pair rapidities a m are rearranged as permitted by (9). The yrast excitations are then obtained by moving a dimer from the inside of the Fermi sphere to the edge and thus creating a hole, as illustrated in figure 2 on the right by the filled symbols. Compared to particle-like or mixed particle-hole-like excitations where a dimer rapidity occupies one of the free slots outside the Fermi sphere, the hole excitation always has lower energy at given momentum in direct analogy to the Lieb-Liniger and Tonks-Girardeau models [48]. For the rapidities this means omitting one element from the set of ground-state ℓ m quantum numbers of (13) and adding + ℓ M 1 . This operation shifts one of the a m values and with it, the corresponding pair of charge rapidities with values a  i m c 2 . It is not yet obvious from the asymptotic representations shown in figure 2, whether and how the dimer hole excitations are connected to the single fermion holes we have discussed in the context of weak attractive interactions. For both to be the yrast excitations of the respective regime, it is necessary that they are continuously connected through a change in the interaction parameter. That this is indeed the case is demonstrated in figure 4 where the real and imaginary values of the rapidities are shown as a function of γ for N=6 particles. Note that the real parts of the charge rapidities converge to the values of the spin rapidities and the imaginary parts become linear in c for strongly negative c, as required by the string hypothesis. In addition, we observe that the pair of k j ʼs that are separated to form the single fermion hole at weak interactions always merge to form a dimer at more negative values of c. Another observation is that the spacing of rapidities on the real axis changes from the single-particle quantization condition p L 2 at c=0 to p L at  -¥ c , which is the appropriate momentum quantization for dimers with mass m 2 . The dispersion relation for the yrast excitation energy = -E E E tot GS versus momentum is shown in figure 5. It is seen that the excitation energy of each yrast state decreases as the interaction strength changes from noninteracting to strongly attractive within a well defined interval. The limits of the interval are indicated by the full parabolic lines drawn from (15) and (16) below, respectively. The excitation at momentum deserves special attention. This value is half of the umklapp momentum p 4 F discussed in section 2.1 and is also the momentum value obtained by adding the unit value p L 2 to every singleparticle wavenumber of half the fermions, e.g. to only the spin-up fermions. In the noninteracting case, the corresponding yrast excitation is obtained by taking a spin-up fermion from the left side of the Fermi sphere and attaching it to the right side, which is equivalent to shifting the Fermi sphere by one unit wavenumber while the spin-down Fermi sphere remains fixed. As can be expected, the energy of this excitation (for g = 0) lies on the parabola P mM 2 2 , which we refer to as the half-system translation parabola. It is shown as a dashed-dotted line in figure 5. For finite attractive interaction strength g < 0, the energy of this excitation decreases and eventually moves to the full system translation parabola P mN 2 2 , which is shown as a dashed line. This is to be expected, since we know that the system dimerizes for g  -¥, and an umklapp excitation of a Tonks-Girardeau gas of = M N 2 dimers can be obtained at momentum p 2 F with energy´( Since we know the ground-and excited-state rapidities in both limits, it is possible to find the asymptotic dispersion relations for weak and strong interactions

=
is the Fermi momentum of the free system. These dispersions are added to figure 5 as solid lines.
Moreover, figure 5 compares the exact and string hypothesis dispersions at g = -2, roughly the lowest c value that can be tackled by the exact exponential equations at the given number of particles. The small difference visible between the dispersions is due to the fact that we have M complex k j pairs, which directly contribute to the energy and their imaginary parts in the exact and approximate equations are different.

Connection to the super Tonks-Girardeau gas
The gaseous phase of attractive bosons known as the super Tonks-Girardeau gas can be described by the regular Bethe ansatz equations for the Lieb-Liniger model of bosons in 1D under the additional restriction that all rapidities are real valued such that bound states of bosons are excluded [65]. These equations map one-to-one to the equations of the attractive Yang-Gaudin model (3) and (4) as well as (9)  , and dimensionless coupling parameter . Thus all results of this work pertaining to the dimerized attractive spin- 1 2 Fermi gas apply to the super Tonks-Girardeau gas as well. An important physical difference between the models is that bound states of fermion dimers are excluded by the Pauli exclusion principle [45], while bound states of bosons exist in a lower energy sector. The yrast states of the spin-  In the finite-system case we saw that the string hypothesis works when | | c L 1  and the exact Bethe ansatz equations can be replaced by the logarithmic form (9) for fully dimerized states. In the thermodynamic limit, this condition is fulfilled for any value of the interaction strength γ since  ¥ L , and thus we only need to consider the logarithmic equations. These coupled algebraic equations become Fredholm integral equations of the second kind, first derived by Gaudin [70]. The ground state integral equations can be found in [72] and there is some earlier work examining excitations in the Yang-Gaudin model [70,[73][74][75], often focusing on spinchanging excitations. In the following we present the integral equations that are appropriate for the ground state and yrast excitations. These are solved numerically using the Fie package for MATLAB [76].

Ground state
The balanced (fully dimerized) ground state properties can be found by introducing the density distribution function a ( ) r for the dimer rapidities, and rewriting (9) into the integral equation [70,72] ò p a p a a a a = -¢ where b is the boundary of the rapidity distribution. We will follow the common convention [72] to redefine the reference point for energy measurements by subtracting the binding energy - Mc 2 m 2 2 2 , which diverges in the strongly attractive limit. This does not affect the dispersion relation, which reports the energy difference between excited and ground state. The ground state energy and momentum then become In order to solve the integral equation numerically we scale the equations, introducing the dimensionless rapidity x and the parameter λ: We thus obtain with the normalization condition where = M N 2 has been used. The parameter λ has a monotonic relationship to the dimensionless interaction strength γ, as is shown in figure 6. For numerical calculations it is practical to set λ first and then obtain the value of γ from (23).
The ground state energy now can be written as is defined by x gx x 2 d. Knowledge of the ground state energy provides us with access to the chemical potential from the binding energy of the dimer that has to be broken in order to remove a single fermion. The chemical potential for the (bosonic) dimers is m 2 .

Yrast dispersion relation
The yrast excited states are characterized by a single hole in the set of dimer rapidities but the rapidities are shifted from their ground state positions. An integral equation can then be derived for a function a ( ) J , which quantifies the shift [48,70]. Appendix B provides some further details on this derivation. The integral equation has a similar structure to that of the ground state but contains the hole rapidity Q as an additional parameter with -< < b Q b: ò p a a a a a p q a = ¢ The energy and momentum of the excitation relative to the ground state are found from The equations can be made dimensionless by a similar rescaling as for the ground state with the additional dimensionless quantities The integral equation becomes   Figure 7 shows example dispersions for this excitation branch. Limiting expressions for weak and strong interactions are easily obtained from the finite particle number results (15) and (16), respectively, by taking the limit  ¥ N . Note that the translation parabolas µP Nm 2 move to the P-axis in the thermodynamic limit as  ¥ N . Thus both the umklapp and half-umklapp excitations have zero energy and the dispersion relation The yrast dispersion relation is gapless and, as demanded from general principles, the slope at the origin is the speed of sound as in the Lieb-Liniger model [48]. A dispersion relation for phonon-like 'particle' excitations with the same slope at the origin but higher energy than the yrast excitations can also be obtained from the Bethe ansatz [70]. In addition there is a triplet of elementary excitations that involve breaking a single dimer pair. The corresponding dispersion relations are gapped for g < 0 and also have higher energy than the yrast excitations.
An explicit form of the dispersion relation can be obtained in the strongly attractive regime by solving the integral equations to leading order in l -1 . Neglecting -¢ ( ) x x 2 compared to l 2 in the kernel, we find p l = --( ) g 2 1 , a constant independent of x. Expanding the inhomogeneous term in the excited-state equation to first order about zero, naturally leads to a linear ansatz for h(x): The excitation energy and momentum can also be obtained, giving Eliminating q and substituting for λ in terms of γ, the dispersion relation becomes This result agrees with (16) for g  -¥ and  ¥ N , as it should.

Missing particle number and phase step
The dispersion relations discussed in the previous section give us access to a number of quantities that define the character of soliton-like excitations, and that can be derived assuming the existence of a localized quasiparticle. In the following we will consider m = ( ) E E P , as a function of the two variables μ and P. Note that both the density n 0 and g = c n 0 depend on the chemical potential μ while parameters of the Hamiltonian, in particular c, do not.
Associating the yrast excitations with a density dip in the inhomogeneous density n s (x) with background density n 0 allows us to define the particle number . The quantity N d can then be obtained from thermodynamic relations [78, equation (5)] and written as An alternative derivation of (43) based on the Hellman-Feynman theorem in a co-moving reference frame will be published elsewhere [79]. The quantity N d differs in general from [17,61] but gives identical results for v s =0 as can be seen directly from (43). For dark solitons as well as for the yrast excitations considered here, < N 0 d and thus is also referred to as the missing particle number, i.e. the number of particles that make up the difference between the soliton's density profile and a constant background density. Strictly speaking, we are using the term only by analogy to classical solitons, e.g. in the nonlinear Schrödinger equation, since the physical meaning of N d for yrast excitations in general is not fully understood. Likewise it is possible to assign to the states on the yrast dispersion a phase step f D . As a well-known characteristic of dark solitons in the nonlinear Schrödinger equation [3], the phase step becomes relevant in a system with periodic boundary conditions, where it determines the size of the backflow current: the entire superfluid flows at . This counter-flow is induced in order to counteract the phase jump of the soliton and to ensure that the total phase is well defined and continuous across the periodic boundary. f D can be calculated from [17] which states that the total (or canonical) momentum P in a periodic box is the sum of = P mv N s s d , the physical momentum of the soliton, and the counter-flow momentum of the entire superfluid. The phase step f D is thus well defined, even though 1D quantum gases are thought to possess at most a fluctuating condensate phase [80]. With trivial adaption to the current case of the Yang-Gaudin model they read 3 : where the functions g(x) and h(x) were introduced in sections 3.1 and 3.2 and still depend on the interaction parameter γ (and the momentum in the case of h(x)). The excellent agreement with the direct evaluation of (43) is seen by the overlapping curves in figure 8.
In the strongly interacting limit, using the approximate analytical solutions of the Bethe ansatz integral equations and Campbell's formulae, we find the following limiting expressions:

Quasiparticle dynamics in a harmonic trap
Konotop and Pitaevskii have argued that the motion of solitons on a slowly varying background is governed by the principles of Landau's quasiparticle dynamics, which leads to simple classical equations of motion [24]. Here we apply these ideas to yrast excitations in the inhomogeneous 1D Fermi gas. The crucial point is to assume that localized quasiparticles can be formed from the yrast excitations, e.g. by wave-packet formation, which move with a velocity m = ¶ ¶ ( ) v E P P , s . There is good reason to expect this to happen given the evidence connecting yrast excitations to dark solitons in the closely related Lieb-Liniger model as discussed in the introduction.
Suppose now the background variations happen on a length scale R TF that is large compared to the quasiparticle length scale ℓ s and the Fermi wavelength x = n 4 F 0 . Since the quasiparticle velocity v s is bounded by the speed of sound , we can also suppose that the acceleration of the quasiparticle moving over the slowly varying background is small enough to suppress the radiation of energy. In the local density approximation, the chemical potential m m varies locally with position and adjusts to the value of the external potential. Requiring the energy of the quasiparticle m = ( ) E E P , to be conserved during the motion then leads to a Hamiltonian equation of motion for the position X(t) of the quasiparticle where we have identified the inertial mass and the physical mass [34] (see also [16,17,32]) At zero velocity it is closely related to the missing particle number Since both mass parameters are negative and approximately constant at small velocities, (54) is solved by harmonic oscillations with a frequency Ω given by (1). The inertial and the physical masses at v s =0 can be calculated from the dispersion relation and are shown in figure 9 along with their ratio, m m I P . The limiting value of 1 for m I /m P in the weakly and strongly interacting regimes means that quasiparticles will oscillate in the harmonically trapped gas with the same frequency as a single fermion would, or, indeed, as the centre-of-mass coordinate of an interacting cloud would if dipole oscillations were initiated 4 . This is not unexpected, since both limits correspond to that of non-interacting fermions: g = 0 is a free gas of spin-1 2 fermions and g  -¥ is a Tonks-Girardeau gas of dimers, which maps to a gas of spinless fermions by the Bose-Fermi mapping [50]. The yrast excitations become free single particle holes in these limits, single fermion holes for g = 0 and single dimer holes for g  -¥, and thus should oscillate with the trap frequency.
In the strongly attractive regime we can obtain expressions for the masses and their ratio from the analytical solutions of the integral equations in the limit g  -¥: Due to the separability of centre-of-mass and relative motion in a harmonic trapping potential, this is true regardless of the strength of the inter-particle interaction. . Note that the strongly attractive regime is the one where the super Tonks-Girardeau gas has the longest life time [64].
Our numerical results of figure 9 indicate that the physical mass monotonically decreases from a value of -m in the weakly interacting tom 2 in the strongly attractive regime. The inertial mass, even though it takes the same limiting values, appears to have a faint maximum, reaching values ofm 0.990 around g » -0.7. An interesting regime is the crossover around g » -1 where the mass ratio takes a minimum value of » m m 0.78 I P around g » -1.3 while it stays below unity for the whole interaction range. The negative physical mass and missing particle number N d confirm that the yrast states correspond to holelike or dark-soliton-like excitations even though the dimer-dimer coupling is negative (see section 2.4). This may seem surprising since attractive interactions usually give rise to bright solitons. For the attractive Yang-Gaudin gas (and by equivalence the super Tonks-Girardeau gas), however, the compressibility is positive in the whole interaction regime [see equation (27)], which implies that it is stable against collapse into bright solitons. The positive compressibility is a consequence of Fermi pressure for the Yang-Gaudin gas, and a corresponding kinetic energy pressure for the super Tonks-Girardeau gas. A manifestation of the effectively attractive interaction between dimers is the fact that the physical and inertial mass parameters are reduced compared to the value of the bare dimer hole massm 2 .

Discussion and conclusions
In this work we have identified the yrast excitations of the spin-1 2 Fermi gas with attractive δ interactions in the crossover from the dimerized regime of strong attraction to the non-interacting limit. We have further examined properties of the yrast dispersion relation in the thermodynamic limit obtained from the exact Bethe ansatz solutions of the Yang-Gaudin model. These predictions can also be directly applied to the super Tonks-Girardeau gas, on account of its mathematical equivalence to the fully dimerized solutions of the Yang-Gaudin model. The yrast dispersion relations very much resemble those of dark solitons in the nonlinear Schrödinger equation and the yrast dispersion of the Lieb-Liniger model of the 1D Bose gas. Indeed, by the nature of the yrast states being the lowest energy state at given momentum, we can expect these excitations to be robust and long lived. Yet there are important quantitative differences between the results of the current work and the repulsive Bose gas. For example, the missing particle number for yrast excitations in the Bose gas can be macroscopically large as g --N 2 d 1 2 diverges for small values of the interaction parameter g  0 of the Lieb-Liniger model [61], whereas in the Yang-Gaudin gas > -N 2 d indicates that hole excitations are microscopic and smaller than a single bare dimer. Also, in contrast to the Bose gas where the effective mass ratio m I /m P is always larger than  unity [61], the ratio falls below one for the Yang-Gaudin gas, which means that yrast quasiparticles would oscillate with a frequency that is slightly faster than the trap frequency. These interesting effects are challenging to observe in experiments, since a 1D Fermi gas would have to be realized and manipulated at the level of single atoms. Important steps have already been taken by realizing one-dimensional Fermi gases [82-84] and singleparticle-level detection with high-numerical aperture imaging [85,86]. A possible route towards measuring the oscillation frequency of the yrast quasiparticles in the Yang-Gaudin regime might be provided by the connection of the Yang-Gaudin with the Lieb-Liniger regimes in the onedimensional gas across the confinement-induced resonance [42][43][44] following loosely the procedure applied in [32]: Yrast excitations could be excited by phase imprinting [5] or combined phase and density imprinting [87] and detected through a visible density depletion in the weakly interacting Lieb-Liniger regime for tightly bound dimers where dark solitons become macroscopic. A rapid ramp of interaction strength through the confinement-induced resonance would then adiabatically transform the dark solitons into the yrast states of the Yang-Gaudin model, thus allowing them to be probed. A similar procedure is feasible with a one-dimensional Bose gas with the aim of probing the yrast excitations of the super Tonks-Girardeau gas. A ramp through the confinement-induced resonance into the super Tonks-Girardeau regime has already been demonstrated [64].
In this appendix we provide all the details on how the Bethe ansatz equations were solved numerically. We used the MATLAB environment, and in particular the fsolve function, included in the optimization toolbox.
In order to solve any of the logarithmic equations (where all variables are necessarily real and distinct), we always begin from the strong-coupling limit, taken as g = | | 100. In this regime, one can use p L 2 -multiples of the quantum numbers as an initial guess, and the solver easily picks up the correct solution. Now, usually this guess is sufficient at any γ (with the logarithmic equations), but we find it is better to follow the solutions from g = | | 100 down to whatever value one needs. This is done with an adapting step of g 10, since in the strongcoupling regime the solutions change very little and large steps can be taken, while in the weak-coupling regime the converse is true. At each consecutive γ-step, the guess is taken as the solution at the previous step. This following-in-γ procedure greatly improves efficiency and accuracy for large systems (~-N 100 1000) but is optional for small systems (Ñ 10). Now, solving the exact exponential Bethe ansatz equations with complex rapidities is a somewhat more involved task. Here there are two possibilities: either the structure of the solution changes with γ (for example, distinct real rapidities may merge into a complex-conjugate pair), or it does not. In the latter case, the situation is simpler, so we start there. The simplest of all cases is the ground state, due to the symmetry of the solution. As such, we can track the solution somewhat further than for excited states.
For the ground state we begin from the known approximate solutions at g » 0 and track in γ using a fixed step of 0.01. Once again, we use a linear interpolation based on the previous two points as a guess for the next. Whenever the solver fails to converge to a valid solution, we call the solver a maximum of five times more, at each attempt adding a vector of (real) random numbers to the initially used guess of order 10 −4 . If all five attempts are exhausted unsuccessfully, the program aborts.
For the single fermion branch, the structure of the solution changes as a function of γ. As always, we begin from g = -0.01 where solutions are approximately known, use a fixed step of 0.01 (unless we are in the vicinity of a merging-see below) and linear interpolation for predictive guessing. At each step, the solver is called once and the solutions are inspected. The Bethe ansatz equations must have distinct roots, so that no two k j ʼs and no two a m ʼs are equal. At points where two real rapidities merge into a complex-conjugate pair, the solver often fails to automatically split up these roots along the imaginary axis, so this must be tested for and corrected.
Thus, having made an initial attempt to solve at a given γ value, we test whether any of the k j ʼs have been returned equal (in this case they will also necessarily be real). If so, we modify the guess by splitting up these rapidities by  in 0.1 0 . Having modified the guess to help the solver pick up the correct solution, we call it again and reduce step size to 0.001 in a γ-interval of 0.01 immediately following the merging point. The step is reduced because usually variables change very rapidly in the vicinity of a merging and in order for our linear guess to be effective, the γ-step must be smaller. Once we perform 10 of these reduced steps, the step size is returned to normal.
Regardless of whether a merging takes place at any given γ value, if the solver fails to converge, we make a maximum of ten further attempts to solve by adding small (real) random numbers of ordern 10 4 0 to the initial guess (which may have already been 'manually' modified due to a merging). If all ten attempts fail, the program is aborted.

Appendix B. Thermodynamic limit equations: outline of derivation
Since the ground state thermodynamic limit equations of the attractive Yang-Gaudin model are available in the literature (e.g. [72]), we will not explain how these are to be obtained from the finite-system discrete equations. However, we will indicate how one can obtain the excited-state equations starting from the ground state.
In the case of single fermion holes (which in the thermodynamic limit are dimer holes), we begin from the most convenient ground state, that of a system with + M 1 dimers. This chosen initial ground state is then perturbed by explicitly removing a dimer to create a hole excitation. This perturbation to the system shifts all other rapidities, and this shifted ground state will be denoted by an asterisk. The explicit removal of a dimer, together with the shift of all other rapidities, will yield the excited state of interest. This excited state must be compared to the initial ground state from which we constructed it, which is not be the ground state of interest (which has M dimers). We must then correct for the 'wrong' ground state explicitly, which is a much simpler task then beginning from the correct but inconvenient ground state.
In our notation, the equations for the single fermion hole branch are derived by following the recipe: This means we must take the finite-system Bethe ansatz equations 5 of + M 1 dimers, remove one charge rapidity (denoted here by a 1 ) with the other rapidities shifted from their true ground state values (denoted by the asterisk) and compare the resulting system to the true ground state of + M 1 dimers. Finally, we add a ground state correction that ensures that both the excited and the ground states have = N M 2 fermions (in practice, the ground state correction simply contributes a m 2 term to the excitation energy since the energy cost of adding/removing two fermions to the ground state is precisely m 2 ). Now let us outline in more detail how the equations are derived: the shift in a m due to the perturbation is written as da L m . We write down the finite-system Bethe ansatz equations for the chosen initial ground state, shifting all the rapidities. Then we explicitly subtract any interaction terms that are associated with the rapidities that must be removed, which are themselves not shifted. Next, we expand the interaction terms (in the case of the δ-function potential, the two-body phase-shifts are given by the θ function of (10)) as a first order Taylor expansion about the unshifted rapidities, with the first order terms proportional to the shifts.
We now subtract the Bethe ansatz equations of the unperturbed, conveniently chosen ground state. At this point we may pass to the continuous limit by using the ground state quasi-momenta density function that appear in the ground state thermodynamic limit equations with the following rule: Some of the terms appearing under the integrals in the equations can be simplified by substituting for them from the ground state thermodynamic limit integral equations. Finally, we define a new function that is the product of L, the ground state distribution, and the shift: in particular a da a = ( ) ( ) J L r . It remains to compute the energy and momentum. This is done according to our recipe, by writing down the finite-system total P and E of the excited state (shifted and with an explicit rapidity removed), subtracting those of the convenient ground state, and adding the ground state correction. The resulting quantities can be expressed through the 'shift function' of the integrals equations (i.e. a ( ) J ), the explicit rapidity removed and the chemical potential. . Taking the derivative of (C.6):