Excitations in the Yang-Gaudin Bose gas

# Excitations in the Yang-Gaudin Bose gas

Neil J Robinson Condensed Matter Physics and Materials Science Division,
Brookhaven National Laboratory, Upton, NY 11973-5000, USA
E-mail: nrobinson@bnl.gov
Robert M Konik Condensed Matter Physics and Materials Science Division,
Brookhaven National Laboratory, Upton, NY 11973-5000, USA
E-mail: rmk@bnl.gov
August 25, 2019
###### Abstract

We study the excitation spectrum of two-component delta-function interacting bosons confined to a single spatial dimension, the Yang-Gaudin Bose gas. We show that there are pronounced finite-size effects in the dispersion relations of excitations, perhaps best illustrated by the spinon single particle dispersion which exhibits a gap at and a finite-momentum roton-like minimum. Such features occur at energies far above the finite volume excitation gap, vanish slowly as for fixed spinon number, and can persist to the thermodynamic limit at fixed spinon density. Features such as the gap also persist to multi-particle excitation continua. Our results show that excitations in the finite system can behave in a qualitatively different manner to analogous excitations in the thermodynamic limit.

The Yang-Gaudin Bose gas is also host to multi-spinon bound states, known as -strings. We study these excitations both in the thermodynamic limit under the string hypothesis and in finite size systems where string deviations are taken into account. In the zero-temperature limit we present a simple relation between the length -string dressed energies and the dressed energy . We solve the Yang-Yang-Takahashi equations numerically and compare to the analytical solution obtained under the strong couple expansion, revealing that the length -string dressed energy is Lorentzian over a wide range of real string centers in the vicinity of . We then examine the finite size effects present in the dispersion of the two-spinon bound states by numerically solving the Bethe ansatz equations with string deviations.

Keywords: Quantum integrability (Bethe Ansatz), Yang-Gaudin model, Bound States, Finite-size Effects

## 1 Introduction

### 1.1 Background

Quantum integrable models provide firm ground from which one can gain understanding of the physics of strongly correlated systems. They include paradigmatic examples of magnetism (the spin- Heisenberg model [1, 2]), interacting electrons on the lattice (the Hubbard model [3, 4]), and interacting Bose and Fermi gases (the Lieb-Liniger [5][7] and Yang-Gaudin [8][11] models). As a result of their exact solutions, integrable models can be used to study the non-perturbative effects of interactions on quasi-particle excitations, hopefully paving the way for insights in more generic cases.

As a starting point for attacking problems with strong correlations, it is clearly desirable to understand the nature and characteristics of quasi-particle excitations in quantum integrable models. By now there is a well-trodden path for such studies: the logarithmic Bethe ansatz equations suggest a natural definition of excitations in terms of the Bethe roots (and their defining integer quantum numbers) [1]. One can then extract the dispersion relation for single-particle excitations by direct numerical computation in finite-size systems, perturbative calculations in the weak or strong coupling limits, or by working directly in the thermodynamic limit and performing manipulations of systems of coupled integral equations. The quasi-particle excitations can be combined to describe the continua of multi-particle excitations. Furthermore, the Bethe ansatz equations also provide a description of emergent multi-particle excitations, such as bound states, which are characterized by complex Bethe roots [1, 2, 4, 7]. The presence of such bound states is model dependent (for example, there are no multi-particle bound states in the repulsive Lieb-Liniger model [5][7], although they are present in the attractive limit [7]). There may be multiple types of bound states, such as the well known - and strings in the one-dimensional Hubbard model [4].

Except for some special cases, numerical computations are usually a necessity – especially in the case of multicomponent systems, where the system of Bethe ansatz equations is nested (e.g., there are additional sets of auxiliary Bethe roots), see for example Ref. [12]. Often such calculations are performed on a finite-size systems, so it is useful to understand how the properties of solutions change with the system size (so-called finite-size effects). The study of finite-size effects can be a useful theoretical tool in the study of critical phenomena, allowing one to extract the central charge [13, 14], critical exponents [15][18] and the operator content [19, 20] of conformal field theories. They also have extensive applications in integrable quantum field theories [21] and have recently received attention for the computation of threshold singularities in integrable lattice models [22][27]. Of course, understanding the properties of finite-size systems is also useful when comparing the result of theoretical calculations to experiments on intrinsically finite-size systems.

Interest in quantum integrable models has recently undergone a resurgence (see, for example, the recent reviews [28][37]), thanks to ground-breaking progress in the field of ultra-cold atomic gases [38]. Experiments have achieved both unprecedented levels of isolation from the environment and control of the Hamiltonian, which have lead to extremely accurate realizations of oft-studied theoretical models [38]. Thanks to this, integrable quantum systems are now routinely studied in the laboratory111To be more precise, the experiments on cold atomic gases are weakly non-integrable, in the sense that integrability is (weakly) broken by small inhomogeneities, the presence of a trapping potential, boundary conditions, and so forth. However, it has been understood theoretically that such weak breaking of integrability can lead to dynamics that remain approximately integrable for very long periods of time, see for example, Refs. [35],[39][54]., with both their equilibrium properties [38, 59, 60] and their non-equilibrium dynamics [61][63] receiving a great deal of scrutiny. In the latter case, the ground-breaking experiments of Kinoshita, Wenger and Weiss [64] highlighted the dramatic consequences of integrability on non-equilibrium dynamics; integrable systems driven out of equilibrium do not thermalize, but instead equilibrate to a generalized Gibbs ensemble whose form is fixed by the initial expectation values of local and quasi-local conservation laws [65, 66].

### 1.2 This work

In this work we study properties of the excitations in the Yang-Gaudin Bose gas, with a particular focus on finite-size effects and the multi-particle bound states. In the first part of this paper, we will focus our attention on the spinon dispersion, and the two-particle continua for holon-antiholon () and spinon-holon () excitations. These will be studied using both the exact numerical solution of the Bethe ansatz equations and the strong coupling expansion. Whilst these excitations have previously received some attention [12], only small finite systems were considered and there was no study of finite-size effects. We will show that the finite-size effects in this system are large and quite surprising: excitations in the finite system can behave in a qualitatively different manner to analogous excitations in the thermodynamic limit. This is particularly well illustrated by the spinon dispersion, which exhibits a pronounced finite-momentum roton-like minima in small systems which vanishes in the thermodynamic limit (, with the number of spinons). We study how such features evolve with system size and present numerical evidence that the finite-momentum roton-like minima can persist to the thermodynamic limit of the two-component gas provided the spinon density does not vanish ( with , fixed, being the number of particles).

Our interest in understanding the simple few-particle excitations of the Yang-Gaudin Bose gas stems from recent works on the non-equilibrium dynamics of a distinguishable impurity in the Bose gas [55, 56]. The dynamics of an impurity exhibit some interesting features: an initially localized impurity can undergo arrested expansion or may move through the gas in a snaking motion [55]. Furthermore, for initial states containing a superposition of spinon excitations the spreading shows clear signs of a ‘double light cone’, which is related to the presence of a roton-like minima in the spinon dispersion [56]. As the impurity limit of the Yang-Gaudin Bose gas (, finite) is perhaps the simplest limit to consider and already exhibits unusual non-equilibrium behavior, it is necessary to develop a good understanding of the excitations of the model. There has also been interesting recent work which relates the non-relativistic limit of various integrable relativistic quantum field theories to multi-component Lieb-Liniger and Yang-Gaudin models [57, 58].

In the remainder of the paper, we turn our attention to the multi-spinon bound states present in the Yang-Gaudin Bose gas. Starting in the thermodynamic limit, we numerically solve the Yang-Yang-Takahashi equations to compute the dispersion relation for bound states, so-called -strings. We present a simple relation between the dressed energies of the -strings and antiholon dressed energy , which we solve under a strong coupling expansion to reveal a particularly simple closed form. We then consider finite size effects for the bound state excitations in small systems, taking into account string deviations.

We finish with a discussion of our results, including implications for comparison between theory and experiments in cold atomic gases.

### 1.3 The Yang-Gaudin Bose gas

The Yang-Gaudin Bose gas222This is often also called the two-component Lieb-Liniger model or the spinor (two-component) Bose gas. The Yang-Gaudin model can also refer to the same Hamiltonian with spin- fermionic fields. Here we will only discuss two-component bosons. is described by the Hamiltonian density

 H=ℏ22m∑j=1,2∂xΨ†j(x)∂xΨj(x)+c∑j,l=1,2Ψ†j(x)Ψ†l(x)Ψl(x)Ψj(x), (1)

where label the two different boson species, is the boson mass, and characterizes the interaction strength. The bosonic fields obey canonical commutation relations,

 [Ψj(x),Ψ†l(y)]=δj,lδ(x−y). (2)

Herein we set . This model is integrable, and may be solved via the nested Bethe ansatz [8][11]. The particle eigenstates, with particles of the second species, are described by two sets of quantum numbers: the momenta and the spin rapidities . The spin rapidities are often known as auxiliary Bethe roots, as they do not directly enter into expressions for the momentum or energy of the eigenstates, see Eqs. (7) and (8). The quantum numbers satisfy the Bethe ansatz equations, which read in their logarithmic form [8][10]

 2πIj = kjL+N∑l=1ϕ1(kj,kl)−M∑β=1ϕ2(kj,Λβ), (3) 2πJα = N∑l=1ϕ2(Λα,kl)−M∑β=1ϕ1(Λα,Λβ). (4)

Here we have defined the scattering phase

 ϕn(u,v)=\rmilog(\rmic+(u−v)n\rmic−(u−v)n)≡2arctan((u−v)nc), (5)

and the sets of ‘integers’ which obey

 Ij,Jα∈{Z,if (N+M)∈2Z+1,Z+12,if (N+M)∈2Z. (6)

The eigenstate associated with the sets of integers , has momentum and energy given by

 P({I},{J}) = ∑jkj=2πL⎛⎝∑jIj−∑βJβ⎞⎠, (7) E({I},{J}) = ∑jk2j. (8)

The zero-temperature ground state of the Yang-Gaudin Bose gas is fully polarized: it is found in the sector with , which follows from general symmetry considerations [67, 68].333A system of multi-component interacting bosons with component-independent repulsive interaction will have a ferromagnetic ground state [67, 68]. With even (and herein we take to be even), the -particle ground state is described by the set of momenta integers

 IGS={−N−12,−N−12+1,…,N−12}, (9)

which forms a ‘Fermi sea’ about the origin. As the ground state is fully polarized, it coincides with the ground state of the Lieb-Liniger model for a single component Bose gas [5][7].

### 1.4 Low-energy excitations of the Yang-Gaudin Bose gas

The Bethe ansatz equations (3),(4) and the integers (6) provide a natural definition for excitations. Restricting our attention to states containing particles, and starting from the absolute ground state (9), we can construct the following types of excitations [12]:

1. The spinon () excitation. This corresponds to the presence of a single spin rapidity in the Bethe ansatz equations (3),(4), characterized by the integer . According to the rules (6), the momenta integers shift and we have the following configuration

 Is={−N2,…,N2−1},|Js|≤N2. (10)

This corresponds to a ‘spin wave’ excitation [69, 70], where the species index plays the role of spin (accordingly, this excitation is sometimes called an isospinon). We show an example configuration for particles below:

\integers\border\intrange

-98 \fillI-43 \intLabelI \render

\integers\intrange

2. The holon-antiholon () excitation. Here we start from the ground state configuration of integers (9) and remove the th integer [leaving a single hole at ]. We then add an integer outside the Fermi sea:

 Ih¯h={−N−12,…,−N−12+j−1,−N−12+j+1,…,N−12,IN}, with  IN>N−12  and  1≤j≤N. (11)

We illustrate this configuration (with the position of the hole in the sea of integers) for particles below:

\halfintegers\border\intrange

3. The spinon-holon () excitation. This two-particle excitation is similar to the combination of the above two excitations. We consider the configuration of integers,

 Ish={−N2,…,−N2+j−1,−N2+j+1,…,N2}, Jsh≤N2,and1≤j

i.e., we consider a symmetric Fermi sea of integers containing a single hole at position accompanied by a single spin rapidity. For particles this configuration has the following diagrammatic depiction:

\halfintegers\border\intrange

\integers\intrange

where is the position of the hole in the Fermi sea of integers.

4. Spinon bound states (-strings). When the system contains more than one flipped spin the associated spin rapidities can become complex. Such solutions are arranged in regular patterns in the complex plane, known as strings [1, 71]. An -string consists of spin rapidities which share the same real part

 Λn,aα=ωnα+ic2(n+1−2a)+iζn,aα. (13)

are known as the ‘string deviations’, which are non-zero in the finite-size system. It is usually assumed that the string deviations are exponentially small in the system size , and in the thermodynamic limit can be set to ; this is known as the ‘string hypothesis’. Rapidities in each string are then spaced evenly in the complex plane.

## 2 The spinon single-particle dispersion

We consider the case with particles, of which there is a single particle of the second kind. In this case, the Bethe ansatz equations (3) and (4) become particularly simple [72]. In their logarithmic form they read

 2πIj = kjL+N∑l=1ϕ1(kj,kl)−ϕ2(kj,Λ), (14) 2πJ = N∑l=1ϕ2(Λ,kl), (15)

where the ‘integers’ satisfy Eq. (6) and , which follows from the bounding of . The momentum and energy of the eigenstates are as previously described in Eqs. (7) and (8), respectively.

Herein we focus on the case with . We choose conventions where the ground state in the sector with is described by the integers.

 I0≡{−N2,−N2+1,…,N2−1},J0=−N2, (16)

and we exclude from future discussions, to avoid double counting. Diagrammatically, for particles this is:

\integers\border\intrange

-98 \fillI-43 \intLabelI_0 \render

\integers\intrange

Notice that when , it follows from Eq. (7) that the state has finite-momentum and the left/right Fermi points of the set of momenta do not coincide.

### 2.1 Strong coupling expansion

To compute the spin wave dispersion under a strong coupling () expansion, we follow the standard prescription described in, e.g., Ref. [7]. Our aim is to compute the energy associated with the presence of a spin wave (e.g., the presence of a spin rapidity , or integer ) in the system. We compute the energy of the spin wave state above the ground state configuration (16) with the set of integers fixed.

#### 2.1.1 Integral equation for the shift function.

In order to compute the spin wave dispersion above the ground state, we fix the integers in the ground state configuration and vary the integer from its ground state value to (the ground state value of the spin rapidity is , which corresponds to acting on a Lieb-Liniger eigenstate with a global spin lowering operator [72]). We proceed by taking the difference of the first Bethe equation for the two cases; we denote the ground state momenta by and the excited state momenta by with the accompanying finite spin rapidity . We have

 0 = (kj−k(0)j)L−ϕ2(kj,Λ)+π+N∑l=1[ϕ1(kj,kl)−ϕ1(k(0)j,k(0)l)]. (17)

Here the factor of arises from with .

We now use that the difference between the roots in the presence of the finite spin rapidity and those in the ground state are . We can then expand the scattering phase as

 ϕ1(kj,kl)=ϕ1(k(0)j,k(0)l)+K(k(0)j,k(0)l)[(kj−k(0)j)−(kl−k(0)l)], (18)

where we neglect terms and the derivative of the phase is defined as

 K(u,v)≡ϕ′1(u,v)=2cc2+(u−v)2. (19)

Keeping track of sums which pass to the continuum with symmetric () or non-symmetric () limits (see A), and using

 ∑lϕ1(k(0)j,k(0)l)−∑l′ϕ1(k(0)j,k(0)l)=O(L−1), (20)

we arrive at

 0=(kj−k(0)j)2πLρ(k(0)j)−ϕ2(kj,Λ)+π−∑l(kl−k(0)l)K(k(0)j,k(0)l), (21)

where we have kept terms to order and we used the following identity for the root distribution

 1+1L∑lK(k(0)j,k(0)l)=2πρ(k(0)j). (22)

We now define the ‘shift function’ (see, e.g., Ref. [7])

 Fs(k(0)j|Λ)=kj−k(0)jk(0)j−k(0)j−1, (23)

which can be interpreted physically as measuring the effect of the finite spin rapidity on the ground state momenta . We pass to the continuum ( with fixed, see A for details) and obtain the integral equation

 0=2πFs(k(0)j|Λ)−∫qR−qLdqFs(q|Λ)K(k(0)j,q)−ϕ2(k(0)j,Λ)+π, (24)

where we’ve used that the derivative term in the expansion of about is , thus allowing us to drop it.

#### 2.1.2 The spin rapidity Λ.

Let us recap how the spin rapidity is quantized. For a given integer , the second Bethe equation (15) reads

 2πJ=N∑l=1ϕ2(Λ−kl). (25)

Passing to the continuum according to Eqs. (115), this becomes

 2πJ=L∫qR−qLdkρ(k)ϕ2(Λ,k)−∫qR−qLdkFs(k|Λ)ϕ′2(Λ,k). (26)

The second term is sub-leading in and we neglect it herein; this will be justified a posteriori by direct comparison to the full numerical solution of the Bethe ansatz equations.

#### 2.1.3 The dispersion relation.

From the difference equation, it follows that the energy and momentum of the state (defined with respect to the ground state) with spin rapidity are given by

 Es(Λ)=∫qR−qLdk2kFs(k|Λ),Ps(Λ)=∫qR−qLdkFs(k|Λ), (27)

where satisfies the integral equation (24), with the root distribution determined from Eq. (22) and the solution of Eq. (26). The Fermi momenta can be determined from the following relations:

 NL=∫qR−qLdkρ(k),PL=∫qR−qLdkkρ(k), (28)

where the momentum is given by Eq. (7) and implicitly depends on the spin rapidity .

We now compute the root distribution , the Fermi momenta and the shift function under a expansion. We find

 ρ(k)=12π+12π2c(qR+qL)+O(c−2)≡12π+1πcϱ+O(c−2), (29) qR=πϱ(1+2ϱc)−πL(12+JN)(1−4ϱc)+O(c−2), (30) qL=πϱ(1+2ϱc)+πL(12+JN)(1−4ϱc)+O(c−2), (31) Fs(k|Λ)=12πϕ2(k,Λ)−12+12π2[G(2(qR−Λ)c)−G(2(qL+Λ)c)]−qR+qL2πc+…, (32)

where is the average particle density and we define the function .

Numerically integrating the set of equations (27) containing the expansion for the shift function (32) we obtain the energy and momentum for each state with spin rapidity . We present the dispersion relation in Fig. 1, where we compare to the exact result for bosons on the length ring with interaction parameter . The comparison shows that there is excellent agreement between the expansion computed above and the exact result. This validates our dropping of terms which are sub-leading in and .

At first glance, the spin wave dispersion is rather surprising. Let’s first limit our attention to the well-studied region with . There we recognize the usual quadratic spin wave dispersion which can be understood from general symmetry considerations [73, 74]. The dependence of the effective mass on interaction strength is understood [69, 70] – at strong coupling, the effective mass diverges due to the ‘fermionization’ (e.g., the hard core repulsion) of the bosons

 12m∗=1N+4π2ϱ3c+O(c−2). (33)

There is a simple picture for the mass being proportional to the number of particles in the strong coupling limit. Consider a single flipped spin: due to the hardcore repulsion in the limit, one must move all the other particles on the ring in order to move the particle with different spin. As a result, a flipped spin acts much like a particle with mass  [69].

Moving our attention away from the well-studied region, the spin wave dispersion becomes non-monotonic with a gapped roton-like minima at close to . The excitations about the roton-like minima have the dispersion

 Er(p)=Δr+(p−pr)22mr+…, (34)

where is the momentum of the roton-like minima, is the effective mass for excitations about the minima, is the energy gap (herein the roton gap, in Fig. 1), and the ellipses refer to terms higher order in . Notice that the roton gap is many times larger than , the finite volume excitation gap in the Yang-Gaudin Bose gas,

 δEYG∼(2πL)2(1ϱL+4π2ϱ3c), (35)

which follows from the strong coupling expansion of the spinon effective mass  [69], Eq. (33). For the data presented in Fig. 1 the finite volume excitation gap is .

### 2.2 Finite-size effects: vanishing of the roton-like minima

Now we turn our attention to how the spinon dispersion varies with the system size. As we have seen in the previous section, see Fig. 1, there is a pronounced finite-momentum roton-like minima in the dispersion at close to . The problem that we have considered, particles with a single impurity , is clearly susceptible to finite-size effects: the density of the impurity species is non-zero in the finite system, whilst it vanishes in the infinite volume limit . Accordingly, we may expect some changes in the dispersion with varying system size, although the usual assumption is that there changes will be small and only quantitative – they arise as a result of the finite volume excitation gap – and no qualitative features change. Here, we will show that this assumption is not valid – there are large finite-size effects present in the spinon dispersion for relatively large systems and the infinite volume limit is qualitatively different to the finite volume.

To begin, we present the spinon dispersion at fixed density for a number of system sizes in Fig. 2. We see that the roton-like minima shifts towards with increasing system size, and the roton gap is suppressed. In fact, the roton/ gap is suppressed as a power law in the system size: . The value of the pre-factor can be computed directly from the strong coupling expansion (see the previous section) using that corresponds to and spin rapidity . For large but finite, the shift function (32) then becomes

 limΛ→∞Fs(k|Λ)=−1−12πc(qR+qL), (36)

and hence the energy (27) at is

 Δ2kF=4π2ϱL[1+3ϱc+O(c−2)]. (37)

Here we see both the scaling of the gap and the origin of the large finite-size effects: the pre-factor is large. The large pre-factor arises from the term when the Fermi sea is imbalanced , see Eqs. (30) and (31). Notice that there is not a strong -dependence of the pre-factor; this is consistent with the roton/ gap being present for weak interactions, as seen in small systems in Ref. [12].

### 2.3 Roton-like excitations in the thermodynamic limit

We have seen that the gap (and the roton-like minima) vanish in the limit when the spinon number is fixed. A natural question to ask is whether these features still vanish in the thermodynamic limit with fixed spinon density ( with fixed)? Physically, such a scenario corresponds to the case with finite population imbalance. In this section, we attempt to address this question. First, we will present a non-rigorous argument which suggests that, indeed, the roton-like excitations and gap should persist in the thermodynamic limit with fixed spinon density. We follow this with supporting numerical data and a discussion of the physical origin of the roton-like minima.

#### 2.3.1 A non-rigorous argument.

Here we will show that the Bethe ansatz equations for two almost-identical spinons approximately reduce to twice the Bethe ansatz equations for a single spinon when working at fixed density . Let us begin by writing the single spinon Bethe ansatz equations, defined by the integers and for particles with on the length ring:

 2πI(1)j = k(1)jL+N∑l=1ϕ1(k(1)j,k(1)l)−ϕ2(k(1)j,Λ(1)), (38) 2πJ(1) = N∑l=1ϕ2(Λ(1),k(1)l). (39)

We next consider the case with double the number of spinons, . We want to consider two close-to-identical spinons, so we choose neighboring values for the spin rapidity integers: . As we have doubled the number of spinons, we must also double and to continue working at fixed particle and spinon density. Recall that this means the range of the integers and is also doubled in comparison to the case with . The Bethe ansatz equations in this case read

 2πI(2)j = 2k(2)jL+2N∑l=1ϕ1(k(2)j,k(2)l)−2∑β=1ϕ2(k(2)j,Λ(2)β), (40) 2πJ(2)α = 2N∑l=1ϕ2(Λ(2)α,k(2)l)−2∑β=1ϕ1(Λ(2)α,Λ(2)β),   α=1,2. (41)

There are now a number of points to note. Firstly, taking the sum and difference of the lower equations, we have

 2π(2J(2)1+1)=2N∑l=12∑α=1ϕ2(Λ(2)α,k(2)l), (42) (43)

Taking the continuum limit, the first term on the right hand side of Eq. (43) acquires a factor of [cf. Eq. (26)] which implies that . Accordingly, we parameterize the spin rapidities by

 Λ(2)α=Λ(2)+(−1)αdΛwithdΛ∼O(L−1). (44)

We can now expand the right hand side of Eq. (42) to give

 2π(2J(2)1+1)=22N∑l=1ϕ2(Λ(2),k(2)l)+O((dΛ)2). (45)

We can simplify this equation in two further manners. Firstly, as we are working at fixed particle density , the Fermi surface for the momenta quantum numbers is unchanged (up to corrections of order ) and hence increasing the particle number simply increases the density of the momenta within the Fermi surface. We can approximate to leading order the sum as

 2N∑l=1ϕ2(Λ(2),k(2)l)≈2N∑l=1ϕ2(Λ(2),k(1)l). (46)

Secondly, we can approximate by up to an additive factor of unity. As the right hand side will be proportional to in the continuum limit, single factors of one are unimportant. Thus, we can approximate Eq. (42) by

 (47)

which is the same equation as for the spin rapidity in the case, see Eq. (39). In other words, the spinon rapidity is (approximately) the same for the case with and provided the spinon density is identical. Working through the same arguments for Eq. (40), we find that the and hence the energy and momentum of a two spinon states is approximately described by , , as would naively be expected for two almost identical spinons. This gives some support to the idea that the roton gap for the single spinon excitation can persist to infinite volume provided both and are fixed.

#### 2.3.2 Numerical supporting evidence.

The above argument, whilst suggestive that the roton gap persists to infinite volume, is not rigorous. To confirm that this argument is essentially valid, we provide supporting evidence from the numerical exact solution of the full Bethe ansatz equations (3), (4) for the case of “almost identical” spinons (that is, we choose the spin rapidity quantum numbers to be for and for ) at constant particle density and spinon density . In Fig. 3 we present the extracted single spinon dispersion from computation of the energy and momentum of the “almost identical” spinon states. We see that there is an approximate collapse of results upon rescaling by the spinon number , when working at fixed particle and spinon density. This supports the idea that the roton-like minima and gap persist to the thermodynamic limit when the spinon density is finite (while maintaining the limit ).444In Sec. 4.2, we will see that finite temperature in the thermodynamic limit also opens a gap in the dressed spinon energy. The spin wave excitations in the Yang-Gaudin Bose gas at fixed spinon density thus constitute an example of roton-like excitations in an exactly solvable microscopic many-body model, a problem which has attracted attention for over 75 years [75][78].

#### 2.3.3 Origin of the roton-like minima.

Having established the presence of a gap and roton-like minima in the spin wave dispersion, it is natural to ask what is their physical origin? At a technical level, the behavior is quite easy to explain – everything can be understood in terms of the integers appearing in the Bethe equations, Eqs. (14) and (15). To begin, it is useful to consider the ground state with , corresponding to a spin rapidity of . Then, the Bethe equation (14) becomes

 2π(I(0)j+12)≡2π~I(0)j=k(0)jL+N∑l=1ϕ1(kj,kl), (48)

where we define the shifted ‘integers’ . The ground state is realized as the symmetric Fermi sea of the shifted integers , and the solution of Eq. (48) is in one-to-one equivalence with the ground state of the (one-component) Lieb-Liniger model [5][7]. Now, let us turn our attention to the gap at , which corresponds to solutions with large spin rapidity . As a result, the Bethe equation (14) can be approximated by

 2π(I(0)j−12)≡2π(~I(0)j−1)≈k(0)jL+N∑l=1ϕ1(kj,kl). (49)

The solutions to Eq. (49) are formally equivalent to the scenario where the right-most shifted integer of the ground state configuration is scattered across the Fermi sea: . The momentum and energy gap thus coincide with that of the holon-antiholon excitation in the (one-component) Lieb-Liniger model (see also the following section) [5][7].

At a formal level we see that the gap can be pictured in a similar manner to the holon-antiholon excitation of the (one-component) Lieb-Liniger model, but a more physical explanation would be nice. In the above, we clearly see that presence of interactions between the constituent excitations of the Yang-Gaudin Bose gas (e.g., spinons and (anti)holons) plays an important role. The role of interactions can be further elucidated by considering how the spinon dispersion varies with density of spinons, the number of particles, and the system size. To begin, we fix the spinon density and vary the total particle density . We see in Fig. 4 that the gap is proportional to the total particle density , as can be expected from the strong coupling expansion (37). If instead we fix the total particle density and increase the spinon density  (by consider the “almost identical” spinon configuration discussed in the previous sections), we find similar: is proportional to the spinon density. This is also consistent with the strong coupling expansion for (), where enters through the imbalance of the Fermi points, . Combining these results, we see that the gap is proportional to the product of the spinon and particle densities, consistent with the gap being induced by spinon-(anti)holon interactions.

### 2.4 Finite-size effects: the finite volume excitation gap in the Yang-Gaudin Bose gas

Before moving on to discuss two particle excitations, it is useful to briefly discuss the finite volume excitation gap in the Yang-Gaudin model and compare to that in the (one-component) Lieb-Liniger model. Henceforth, when we say the finite volume excitation gap we mean the minimum excitation energy above the ground state (for any type of excitation) when there is a fixed number of particles .555The name reflects the fact that excitations are gapless in the infinite volume limit with fixed density.

In the Yang-Gaudin model at strong coupling, the finite volume excitation gap is set by the spinon excitation and can be derived directly from the strong coupling expression for the effective spinon mass  [69] (see Eq. (33))

 δEYG∼(2πL)2(1ϱL+4π2ϱ3c)+O(c−2), (50)

On the other hand, in the (one-component) Lieb-Liniger model there exist only holon/antiholon excitations. Here, the finite volume excitation gap corresponds to taking the momentum integer at the Fermi surface and moving it to (creating, effectively, a excitation, see the following section). As a result, the finite volume excitation gap of the (one-component) Lieb-Liniger model is

 δELL∼(2πL)2(1−4ϱc)ϱL+O(c−2), (51)

as can be found from, e.g., Eq. (61). We note that the Yang-Gaudin model with reduces to the one-component Lieb-Liniger model and as a result, the holon-antiholon minimum excitation energy in the Yang-Gaudin model is also given by Eq. (51).

Let us now briefly highlight an important point. Setting for clarity, we have

 δEYG∣∣c=∞δELL∣∣c=∞=1N2. (52)

That is, in the hard core limit the finite volume excitation gap of the Yang-Gaudin model (set by the spinon excitation) is always much smaller than that of the (one-component) Lieb-Liniger model.666In the Yang-Gaudin Bose gas at strong coupling, Eq. (52) simply states that the finite volume excitation gap for the spinon excitation is always much smaller than the minimal excitation energy for the holon-antiholon excitation. In other words, there will be many low-energy spinon states below the energy scale in the Yang-Gaudin Bose gas.

### 2.5 Extending beyond 2kF

We have considered the dispersion of excitations described by the configuration of integers (10), in which the momenta integers are described by the ground state configuration and the integer is varied between its bounds . With fixed momenta integers , the bounding of implies the spinon excitation has bounded momentum . To realize spinon-like states with higher momenta, it is necessary to modify the configuration of the momenta integers , creating multi-particle excitations. To extend to the momentum range , we remove () the momenta integer at the right Fermi point and add () it immediately to the left of the left Fermi point, creating a three-particle excitation. For particles, such a configuration of integers is shown below:

\integers\border\intrange

\integers\intrange

In Fig. 5 we present the dispersion relation for the spinon () and the spinon-like () states computed from the numerical solution of the Bethe ansatz equations (14), (15) for bosons on the length ring with interaction parameter . We see that the roton-like minima in the spinon dispersion is indeed a local minima, with the dispersion of the spinon-like states increasing, before the appearance of an additional roton-like minima close to . Computing the dispersion for spinon-like states carrying higher momentum (generated by translating the Fermi sea of integers further to the left), we see that the spinon-like states have dispersion with additional slowly-decaying oscillations superimposed on top of this trend.

We note that the spinon-like excitations considered here are absolutely stable, having an infinite lifetime due to the integrability of the model. This prevents the decay of a spinon-like excitation with energy into two roton-like excitations, an instability that occurs in non-integrable models (and is observed experimentally in high-density liquid Helium [79]) that results in the appearance of a Pitaevskii plateau [80] in the dynamical structure factor. Related phenomena, known as quasi-particle breakdown, is also observed in spin systems when one- and two-particle excitations overlap [81].

## 3 Two particle excitations

Having observed large finite-size effects in the single particle dispersion, we now turn our attention to computation of the two-particle excitation continuum for the and excitations. In particular, we will focus on whether there exists a finite energy gap at in the thermodynamic limit.

### 3.1 The holon-antiholon (h¯h) continuum

The continuum constructed on top of the absolute ground state (9) is identical to the particle-hole excitation continuum of the Lieb-Liniger model (see, e.g., Ref. [7]). Notice also that setting the spin rapidity , we recover the Bethe ansatz equations for the Lieb-Liniger model and so the continuum above this state is also equivalent to that in the Lieb-Liniger model. We consider the problem of removing a particle with momentum from within the Fermi sea and replacing it with a particle outside the Fermi sea . The two-particle continuum of excitations is characterized in terms of the energy and momentum given by [7]

 Eh¯h(kp,kh)=k2p−k2h−∫q0−q0dk2k Fh¯h(k|kp,kh), (53) Ph¯h(kp,kh)=kp−kh−∫q0−q0dkFh¯h(k|kp,kh). (54)

Here we define the shift function defined with ground state momenta and the excited state (with the excitation) momenta . The shift function obeys the following integral equation [7]

 Fh¯h(q|kp,kh)−∫q0−q0dk2πK(q,k)Fh¯h(k|kp,kh)=12π[ϕ1(q−kp)−ϕ1(q−kh)]. (55)

#### 3.1.1 Strong coupling expansion.

We can now compute the expansion for the relevant quantities, as done when considering the spinon dispersion. We find

 ρ(k)=12π(1+2ϱc)+O(c−3), (56) q0=πϱ(1−2ϱc+4ϱ2c2)+O(c−3), (57) Fh¯h(k|kp,kh)=1πc(kh−kp)[1+2ϱc]+O(c−3). (58)

For a configuration of integers , the strong coupling expansion for the momenta is

 kj=2πIjL(1−2ϱc+4ϱ2c2)+O(c−3), (59)

which is consistent with our expression for the Fermi momentum , Eq. (57).

As the shift function is functionally independent of to order , we can directly integrate Eqs. (53) and (54) to obtain

 Eh¯h(kp,kh)=k2p−k2h+O(c−3),Ph¯h(kp,kh)=(kp−kh)[1+2ϱc]+O(c−3). (60)

Denoting the hole in the Fermi sea of integers by and the excited integer by , we have the simple relations

 Eh¯h(Ip,Ih)=(2πL)2(1−4ϱc+12ϱ2c2)(I2p−I2h)+O(c−3), (61) Ph¯h(Ip,Ih)=2πL(Ip−Ih)+O(c−3). (62)

This realises the usual holon-antiholon continuum of the Lieb-Liniger model [7], shown in Fig. 6. From Eq. (61), the sole effect of finite interaction strength [to order ] is to rescale the energy compared to the limit; for ,