Excitations in the Yang-Gaudin Bose gas

TL;DR

This paper investigates the excitation spectrum of the Yang-Gaudin Bose gas, revealing significant finite-size effects, novel features like a $2k_F$ gap and roton minimum, and analyzing multi-spinon bound states both analytically and numerically.

Contribution

It provides a detailed analysis of finite-size effects on excitations, including the $2k_F$ gap, roton minimum, and multi-spinon bound states, using both analytical and numerical methods.

Findings

Pronounced finite-size effects in excitation dispersion relations.

Persistence of $2k_F$ gap and roton minimum at finite density.

Lorentzian form of $ abla$-string dressed energies near $ ext{lambda}=0$.

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 $2k_F$ and a finite-momentum roton minimum. Such features occur at energies far above the finite volume excitation gap, vanish slowly as $1/L$ for fixed spinon number, and can persist to the thermodynamic limit at fixed spinon density. Features such as the $2k_F$ 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 $\Lambda$-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 $n$ $\Lambda$-string dressed energies $\epsilon_n(\lambda)$ and the dressed energy $\epsilon(k)$. We solve the Yang-Yang-Takahashi equations numerically and compare to the analytical solution obtained under the strong couple expansion, revealing that the length $n$ $\Lambda$-string dressed energy is Lorentzian over a wide range of real string centers $\lambda$ in the vicinity of $\lambda = 0$. 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.

Equations

$$\mathcal{H}=\frac{{\hslash }^2}{2m}\sum _{j=1,2}{\partial }_x{\Psi }_j^{†}(x){\partial }_x{\Psi }_j(x)+c\sum _{j,l=1,2}{\Psi }_j^{†}(x){\Psi }_l^{†}(x){\Psi }_l(x){\Psi }_j(x),$$

$$\Big{[}\Psi_{j}(x),\Psi^{\dagger}_{l}(y)\Big{]}=\delta_{j,l}\delta(x-y).$$

$$2\pi I_j$$

$$2\pi J_{\alpha }$$

$${\varphi }_n(u,v)=\text{rmi}\log \left(\frac{\text{rmi}c+(u-v)n}{\text{rmi}c-(u-v)n}\right)\equiv 2\arctan \left(\frac{(u-v)n}{c}\right),$$

$$\displaystyle I_{j},J_{\alpha}\in\left\{\begin{array}[]{lcl}\mathbbm{Z},&&\mathrm{if~{}}(N+M)\in 2\mathbbm{Z}+1,\\ \mathbbm{Z}+\frac{1}{2},&&\mathrm{if~{}}(N+M)\in 2\mathbbm{Z}.\end{array}\right.$$

$$P(\{I\},\{J\})$$

$$E(\{I\},\{J\})$$

$$I_{\mathrm{GS}}=\left\{-\frac{N-1}{2},-\frac{N-1}{2}+1,\dots ,\frac{N-1}{2}\right\},$$

$$I_s=\left\{-\frac{N}{2},\dots ,\frac{N}{2}-1\right\},|J_s|\le \frac{N}{2}.$$

$$I_{h\bar{h}}=\left\{-\frac{N-1}{2},\dots ,-\frac{N-1}{2}+j-1,-\frac{N-1}{2}+j+1,\dots ,\frac{N-1}{2},I_N\right\},$$

$$\mathrm{with}{I}_N>\frac{N-1}{2}\mathrm{and}1\le j\le N.$$

$$I_{sh}=\left\{-\frac{N}{2},\dots ,-\frac{N}{2}+j-1,-\frac{N}{2}+j+1,\dots ,\frac{N}{2}\right\},$$

$$J_{sh}\le \frac{N}{2},\mathrm{and}1\le j<N,$$

$${\Lambda }_{\alpha }^{n,a}={\omega }_{\alpha }^n+i\frac{c}{2}(n+1-2a)+i{\zeta }_{\alpha }^{n,a}.$$

$$2\pi I_j$$

$$2\pi J$$

$$I_0\equiv \left\{-\frac{N}{2},-\frac{N}{2}+1,\dots ,\frac{N}{2}-1\right\},J_0=-\frac{N}{2},$$

$$0$$

$$\displaystyle\phi_{1}(k_{j},k_{l})=\phi_{1}(k_{j}^{(0)},k_{l}^{(0)})+K(k_{j}^{(0)},k_{l}^{(0)})\Big{[}\Big{(}k_{j}-k_{j}^{(0)}\Big{)}-\Big{(}k_{l}-k_{l}^{(0)}\Big{)}\Big{]},$$

$$K(u,v)\equiv {\varphi }_1^{\prime }(u,v)=\frac{2c}{c^2+(u-v{)}^2}.$$

$$\sum_{l}\phi_{1}\Big{(}k_{j}^{(0)},k_{l}^{(0)}\Big{)}-{\sum_{l}}^{\prime}\phi_{1}\Big{(}k_{j}^{(0)},k_{l}^{(0)}\Big{)}=O(L^{-1}),$$

$$0=\Big{(}k_{j}-k_{j}^{(0)}\Big{)}2\pi L\rho\Big{(}k_{j}^{(0)}\Big{)}-\phi_{2}(k_{j},\Lambda)+\pi-\sum_{l}\Big{(}k_{l}-k_{l}^{(0)}\Big{)}K\Big{(}k_{j}^{(0)},k_{l}^{(0)}\Big{)},$$

$$\displaystyle 1+\frac{1}{L}\sum_{l}K\Big{(}k_{j}^{(0)},k_{l}^{(0)}\Big{)}=2\pi\rho\Big{(}k_{j}^{(0)}\Big{)}.$$

$$F_{s}\Big{(}k_{j}^{(0)}|\Lambda\Big{)}=\frac{k_{j}-k_{j}^{(0)}}{k_{j}^{(0)}-k_{j-1}^{(0)}},$$

$$\displaystyle 0=2\pi F_{s}\Big{(}k_{j}^{(0)}|\Lambda\Big{)}-\int^{q_{R}}_{-q_{L}}{\rm d}q\,F_{s}(q|\Lambda)K\Big{(}k_{j}^{(0)},q\Big{)}-\phi_{2}\Big{(}k_{j}^{(0)},\Lambda\Big{)}+\pi,$$

$$2\pi J=\sum _{l=1}^N{\varphi }_2(\Lambda -k_l).$$

$$2\pi J=L{\int }_{-q_L}^{q_R}\mathrm{d}k\rho (k){\varphi }_2(\Lambda ,k)-{\int }_{-q_L}^{q_R}\mathrm{d}k{F}_s(k|\Lambda ){\varphi }_2^{\prime }(\Lambda ,k).$$

$$E_s(\Lambda )={\int }_{-q_L}^{q_R}\mathrm{d}k2k{F}_s(k|\Lambda ),P_s(\Lambda )={\int }_{-q_L}^{q_R}\mathrm{d}k{F}_s(k|\Lambda ),$$

$$\frac{N}{L}={\int }_{-q_L}^{q_R}\mathrm{d}k\rho (k),\frac{P}{L}={\int }_{-q_L}^{q_R}\mathrm{d}kk\rho (k),$$

$$\rho (k)=\frac{1}{2\pi }+\frac{1}{2{\pi }^{2}c}(q_R+q_L)+O\left(c^{-2}\right)\equiv \frac{1}{2\pi }+\frac{1}{\pi c}\varrho +O\left(c^{-2}\right),$$

$$\displaystyle q_{R}=\pi\varrho\Bigg{(}1+\frac{2\varrho}{c}\Bigg{)}-\frac{\pi}{L}\left(\frac{1}{2}+\frac{J}{N}\right)\Bigg{(}1-\frac{4\varrho}{c}\Bigg{)}+O\left(c^{-2}\right),$$

$$\displaystyle q_{L}=\pi\varrho\Bigg{(}1+\frac{2\varrho}{c}\Bigg{)}+\frac{\pi}{L}\left(\frac{1}{2}+\frac{J}{N}\right)\Bigg{(}1-\frac{4\varrho}{c}\Bigg{)}+O\left(c^{-2}\right),$$

$$F_s(k|\Lambda )=\frac{1}{2\pi }{\varphi }_2(k,\Lambda )-\frac{1}{2}+\frac{1}{2{\pi }^2}\left[G\left(\frac{2(q_R-\Lambda )}{c}\right)-G\left(\frac{2(q_L+\Lambda )}{c}\right)\right]-\frac{q_R+q_L}{2\pi c}+\dots ,$$

Full text

The Yang-Gaudin Bose gas is also host to multi-spinon bound states, known as $\Lambda$-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 $n$ $\Lambda$-string dressed energies $\epsilon_{n}(\lambda)$ and the dressed energy $\epsilon(k)$. We solve the Yang-Yang-Takahashi equations numerically and compare to the analytical solution obtained under the strong couple expansion, revealing that the length $n$ $\Lambda$-string dressed energy is Lorentzian over a wide range of real string centers $\lambda$ in the vicinity of $\lambda=0$. 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. 1 Introduction

2. 1.1 Background

3. 1.2 This work

4. 1.3 The Yang-Gaudin Bose gas

5. 1.4 Low-energy excitations of the Yang-Gaudin Bose gas

6. 2 The spinon single-particle dispersion

7. 2.1 Strong coupling expansion

8. 2.1.1 Integral equation for the shift function.

9. 2.1.2 The spin rapidity $\Lambda$.

10. 2.1.3 The dispersion relation.

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

12. 2.3 Roton-like excitations in the thermodynamic limit

13. 2.3.1 A non-rigorous argument.

14. 2.3.2 Numerical supporting evidence.

15. 2.3.3 Origin of the roton-like minima.

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

17. 2.5 Extending beyond $2k_{F}$

18. 3 Two particle excitations

19. 3.1 The holon-antiholon ($h\bar{h}$) continuum

20. 3.1.1 Strong coupling expansion.

21. 3.2 The spinon-holon ($sh$) continuum

22. 4 Bound states of spinons: $\Lambda$-strings

23. 4.1 The Bethe-Takahashi equations

24. 4.2 The Yang-Yang-Takahashi equations

25. 4.3 Zero temperature limit of the YYT equations

26. 4.3.1 Strong coupling expansion in the zero temperature limit.

27. 4.4 Numerical solution of the YYT equations

28. 4.4.1 Results.

29. 4.5 A single 2-string in a finite system

30. 4.5.1 Small systems: string deviations and finite-size effects.

31. 4.5.2 Numerical results.

32. 4.5.3 Strong coupling expansion.

33. 5 Discussion

34. A The continuum limit

35. B $\Lambda$-string dressed energies at zero temperature

36. C Strong coupling expansion for the string dressed energies

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-$1/2$ 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 $k$-$\Lambda$ and $\Lambda$ 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 $(s)$ dispersion, and the two-particle continua for holon-antiholon ($h\bar{h}$) and spinon-holon ($sh$) 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 ($M/L\to 0$, with $M$ 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 ($M,N,L\to\infty$ with $N/L$, $M/L$ fixed, $N$ 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 ($M=1$, $N/L$ 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 $\Lambda$-strings. We present a simple relation between the dressed energies of the $\Lambda$-strings and antiholon dressed energy $\epsilon(k)$, 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-$1/2$ fermionic fields. Here we will only discuss two-component bosons. is described by the Hamiltonian density

[TABLE]

where $j,l=1,2$ label the two different boson species, $m$ is the boson mass, and $c$ characterizes the interaction strength. The bosonic fields obey canonical commutation relations,

[TABLE]

Herein we set $\hbar=2m=1$. This model is integrable, and may be solved via the nested Bethe ansatz [8]–[11]. The $N$ particle eigenstates, with $M$ particles of the second species, are described by two sets of quantum numbers: the momenta $\{k\}\equiv\{k_{1},\ldots,k_{N}\}$ and the spin rapidities $\{\Lambda\}\equiv\{\Lambda_{1},\ldots,\Lambda_{M}\}$. 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. (9) and (10). The quantum numbers satisfy the Bethe ansatz equations, which read in their logarithmic form [8]–[10]

[TABLE]

Here we have defined the scattering phase

[TABLE]

and the sets of ‘integers’ $I_{j},\ J_{\alpha}$ which obey

[TABLE]

The eigenstate $|\{k\};\{\Lambda\}\rangle$ associated with the sets of integers $\{I\}$, $\{J\}$ has momentum $P$ and energy $E$ given by

[TABLE]

The zero-temperature ground state of the Yang-Gaudin Bose gas is fully polarized: it is found in the sector with $M=0$, 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 $N$ even (and herein we take $N$ to be even), the $N$-particle ground state is described by the set of momenta integers

[TABLE]

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 (8) provide a natural definition for excitations. Restricting our attention to states containing $N$ particles, and starting from the absolute ground state (11), we can construct the following types of excitations [12]:

The spinon ($s$) excitation. This corresponds to the presence of a single spin rapidity $\Lambda$ in the Bethe ansatz equations (3),(4), characterized by the integer $J_{s}$. According to the rules (8), the momenta integers shift and we have the following configuration

[TABLE]

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 $N=8$ particles below:

I$$-9$$-8$$-7$$-6$$-5$$-4$$-3$$-2$$-1[math]1$$2$$3$$4$$5$$6$$7$$8

J$$-4$$-3$$-2$$-1[math]1$$2$$3$$4$$J_{s} 2. 2.

The holon-antiholon ($h\bar{h}$) excitation. Here we start from the ground state configuration of integers (11) and remove the $j$th integer [leaving a single hole at $-(N-1)/2+j$]. We then add an integer $I_{N}$ outside the Fermi sea:

[TABLE]

We illustrate this configuration (with $I_{h}=j-(N-1)/2$ the position of the hole in the sea of integers) for $N=8$ particles below:

I_{\alpha}^{1}$$-\frac{15}{2}$$-\frac{13}{2}$$-\frac{11}{2}$$-\frac{9}{2}$$-\frac{7}{2}$$-\frac{5}{2}$$-\frac{3}{2}$$-\frac{1}{2}$$\frac{1}{2}$$\frac{3}{2}$$\frac{5}{2}$$\frac{7}{2}$$\frac{9}{2}$$\frac{11}{2}$$\frac{13}{2}$$\frac{15}{2}$$I_{h}$$I_{N} 3. 3.

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

[TABLE]

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

I$$-\frac{15}{2}$$-\frac{13}{2}$$-\frac{11}{2}$$-\frac{9}{2}$$-\frac{7}{2}$$-\frac{5}{2}$$-\frac{3}{2}$$-\frac{1}{2}$$\frac{1}{2}$$\frac{3}{2}$$\frac{5}{2}$$\frac{7}{2}$$\frac{9}{2}$$\frac{11}{2}$$\frac{13}{2}$$\frac{15}{2}$$I_{h}^{\prime}$$I_{N}

J$$-4$$-3$$-2$$-1[math]1$$2$$3$$4$$J_{sh}

where $I_{h}^{\prime}=-N/2+j$ is the position of the hole in the Fermi sea of integers. 4. 4.

Spinon bound states ($\Lambda$-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 $n$-string $\Lambda^{(n)}_{\alpha}\equiv\{\Lambda^{n,1}_{\alpha},\ldots,\Lambda^{n,n}_{\alpha}\}$ consists of $n$ spin rapidities which share the same real part $\omega_{\alpha}^{n}$

[TABLE]

$\zeta_{\alpha}^{n,a}\in\mathbbm{C}$ 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 $L$, and in the thermodynamic limit can be set to $\zeta_{\alpha}^{n,a}=0$; 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 $N$ particles, of which there is a single $M=1$ 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

[TABLE]

where the ‘integers’ $I_{j},\ J$ satisfy Eq. (8) and $|J|\leq N/2$, which follows from the bounding of $\phi_{2}(u)$. The momentum and energy of the eigenstates are as previously described in Eqs. (9) and (10), respectively.

Herein we focus on the case with $N\in 2\mathbbm{Z}$. We choose conventions where the ground state in the sector with $M=1$ is described by the integers.

[TABLE]

and we exclude $J=N/2$ from future discussions, to avoid double counting. Diagrammatically, for $N=8$ particles this is:

I_{0}$$-9$$-8$$-7$$-6$$-5$$-4$$-3$$-2$$-1[math]1$$2$$3$$4$$5$$6$$7$$8

J$$-4$$-3$$-2$$-1[math]1$$2$$3$$J_{0}

Notice that when $J>-N/2$, it follows from Eq. (9) that the state has finite-momentum and the left/right Fermi points of the set of momenta $\{k_{j}\}$ do not coincide.

2.1 Strong coupling expansion

To compute the spin wave dispersion under a strong coupling ($1/c$) 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 $\Lambda$, or integer $J$) in the system. We compute the energy of the spin wave state $J>-N/2$ above the ground state configuration (18) with the set of integers $I_{0}$ 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 $I$ in the ground state configuration $I_{0}$ and vary the integer $J$ from its ground state value $J_{0}=-N/2$ to $J>-N/2$ (the ground state value of the spin rapidity is $\Lambda_{0}=-\infty$, 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 $\{k_{j}^{(0)}\}$ and the excited state momenta by $\{k_{j}\}$ with the accompanying finite spin rapidity $\Lambda$. We have

[TABLE]

Here the factor of $\pi$ arises from $\phi_{2}(k_{j}^{(0)},\Lambda_{0})$ with $\Lambda_{0}=-\infty$.

We now use that the difference between the roots in the presence of the finite spin rapidity and those in the ground state are $k_{j}-k_{j}^{(0)}=O(L^{-1})$. We can then expand the scattering phase as

[TABLE]

where we neglect terms $O(L^{-2})$ and the derivative of the phase is defined as

[TABLE]

Keeping track of sums which pass to the continuum with symmetric (${\sum}^{\prime}$) or non-symmetric ($\sum$) limits (see A), and using

[TABLE]

we arrive at

[TABLE]

where we have kept terms to order $O(L^{0})$ and we used the following identity for the root distribution

[TABLE]

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

[TABLE]

which can be interpreted physically as measuring the effect of the finite spin rapidity $\Lambda$ on the ground state momenta $k^{(0)}_{j}$. We pass to the continuum ($L\to\infty$ with $N/L$ fixed, see A for details) and obtain the integral equation

[TABLE]

where we’ve used that the derivative term in the expansion of $\phi_{2}(k_{j},\Lambda)$ about $k_{j}^{(0)}$ is $O(L^{-1})$, thus allowing us to drop it.

2.1.2 The spin rapidity $\Lambda$.

Let us recap how the spin rapidity $\Lambda$ is quantized. For a given integer $J$, the second Bethe equation (17) reads

[TABLE]

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

[TABLE]

The second term is sub-leading in $L$ 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 $E_{s}(\Lambda)$ and momentum $P_{s}(\Lambda)$ of the state (defined with respect to the ground state) with spin rapidity $\Lambda$ are given by

[TABLE]

where $F_{s}(k|\Lambda)$ satisfies the integral equation (26), with the root distribution $\rho(k)$ determined from Eq. (24) and $\Lambda$ the solution of Eq. (28). The Fermi momenta $q_{L},q_{R}$ can be determined from the following relations:

[TABLE]

where the momentum $P=-\pi(N+2J)/L$ is given by Eq. (9) and implicitly depends on the spin rapidity $\Lambda$.

We now compute the root distribution $\rho(k)$, the Fermi momenta $q_{R,L}$ and the shift function $F(k|\Lambda)$ under a $1/c$ expansion. We find

[TABLE]

where $\varrho=N/L$ is the average particle density and we define the function $G(x)=x\arctan(x)-\frac{1}{2}\log(1+x^{2})$.

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

At first glance, the spin wave dispersion is rather surprising. Let’s first limit our attention to the well-studied region with $|k|\ll\varrho$. There we recognize the usual quadratic spin wave dispersion $E(k)\sim k^{2}/2m_{*}$ which can be understood from general symmetry considerations [73, 74]. The dependence of the effective mass $m_{*}$ 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

[TABLE]

There is a simple picture for the mass $m_{*}\propto N$ being proportional to the number of particles $N$ in the strong coupling $c\to\infty$ limit. Consider a single flipped spin: due to the hardcore repulsion in the $c=\infty$ 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 $N$ [69].

Moving our attention away from the well-studied $|k|\ll\varrho$ region, the spin wave dispersion becomes non-monotonic with a gapped roton-like minima at close to $|k|\sim 2k_{F}$. The excitations about the roton-like minima have the dispersion

[TABLE]

where $p_{r}$ is the momentum of the roton-like minima, $m_{r}$ is the effective mass for excitations about the minima, $\Delta_{r}$ is the energy gap (herein the roton gap, $\Delta_{r}\approx 1.5$ in Fig. 1), and the ellipses refer to terms higher order in $(p-p_{r})$. Notice that the roton gap is many times larger than $\delta E_{YG}$, the finite volume excitation gap in the Yang-Gaudin Bose gas,

[TABLE]

which follows from the strong coupling expansion of the spinon effective mass $m_{*}$ [69], Eq. (35). For the data presented in Fig. 1 the finite volume excitation gap is $\delta E_{YG}\approx 0.0085$.

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 $k=-2k_{F}$. The problem that we have considered, $N$ particles with a single impurity $M=1$, 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 $L\to\infty$. 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 $L\sim 100$ systems and the infinite volume limit is qualitatively different to the finite volume.

To begin, we present the spinon dispersion at fixed density $\varrho\equiv N/L=2$ for a number of system sizes $L=50\to 500$ in Fig. 2. We see that the roton-like minima shifts towards $2k_{F}$ with increasing system size, and the roton gap $\Delta_{r}$ is suppressed. In fact, the roton/$2k_{F}$ gap is suppressed as a power law in the system size: $\Delta_{r}\approx\Delta_{2k_{F}}\propto 1/L$. The value of the pre-factor can be computed directly from the strong coupling expansion (see the previous section) using that $k=-2k_{F}$ corresponds to $J=N/2$ and spin rapidity $\Lambda=\infty$. For $c$ large but finite, the shift function (34) then becomes

[TABLE]

and hence the energy (29) at $k=-2k_{F}$ is

[TABLE]

Here we see both the scaling of the gap $\Delta_{2k_{F}}\sim 1/L$ and the origin of the large finite-size effects: the pre-factor $4\pi^{2}\varrho\gg 1$ is large. The large pre-factor arises from the term $(q_{R}^{2}-q_{L}^{2})$ when the Fermi sea is imbalanced $q_{R}\neq q_{L}$, see Eqs. (32) and (33). Notice that there is not a strong $c$-dependence of the pre-factor; this is consistent with the roton/$2k_{F}$ 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 $2k_{F}$ gap (and the roton-like minima) vanish in the limit $L\to\infty$ when the spinon number $M=1$ is fixed. A natural question to ask is whether these features still vanish in the thermodynamic limit with fixed spinon density ($L\to\infty$ with $N/L,\ M/L$ 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 $2k_{F}$ 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 $N/L,M/L$. Let us begin by writing the single spinon Bethe ansatz equations, defined by the integers $I^{(1)}_{j}$ and $J^{(1)}$ for $N$ particles with $M=1$ on the length $L$ ring:

[TABLE]

We next consider the case with double the number of spinons, $M=2$. We want to consider two close-to-identical spinons, so we choose neighboring values for the spin rapidity integers: $J^{(2)}_{1}=J_{2}^{(2)}-1$. As we have doubled the number of spinons, we must also double $L$ and $N$ to continue working at fixed particle and spinon density. Recall that this means the range of the integers $I_{j}^{(2)}$ and $J^{(2)}$ is also doubled in comparison to the case with $M=1$. The Bethe ansatz equations in this case read

[TABLE]

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

[TABLE]

Taking the continuum limit, the first term on the right hand side of Eq. (45) acquires a factor of $2L$ [cf. Eq. (28)] which implies that $\Lambda_{2}^{(2)}-\Lambda_{1}^{(2)}=O(L^{-1})$. Accordingly, we parameterize the spin rapidities by

[TABLE]

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

[TABLE]

We can simplify this equation in two further manners. Firstly, as we are working at fixed particle density $N/L$, the Fermi surface for the momenta quantum numbers is unchanged (up to corrections of order $1/L$) 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

[TABLE]

Secondly, we can approximate $J_{1}^{(2)}$ by $2J^{(1)}$ up to an additive factor of unity. As the right hand side will be proportional to $L$ in the continuum limit, single factors of one are unimportant. Thus, we can approximate Eq. (44) by

[TABLE]

which is the same equation as for the spin rapidity in the $M=1$ case, see Eq. (41). In other words, the spinon rapidity is (approximately) the same for the case with $M=1$ and $M=2$ provided the spinon density is identical. Working through the same arguments for Eq. (42), we find that the $k_{2j}^{(2)},k_{2j+1}^{(2)}\approx k_{j}^{(1)}$ and hence the energy and momentum of a two spinon states is approximately described by $E^{(2)}\approx 2E^{(1)}$, $P^{(2)}\approx 2P^{(1)}$, 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 $N/L$ and $M/L$ 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 $M=2,3$ “almost identical” spinons (that is, we choose the spin rapidity quantum numbers to be $J_{1}=J_{2}-1$ for $M=2$ and $J_{1}=J_{2}-1=J_{3}-2$ for $M=3$) at constant particle density $N/L$ and spinon density $M/L$. 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 $M$, when working at fixed particle and spinon density. This supports the idea that the roton-like minima and $2k_{F}$ gap persist to the thermodynamic limit when the spinon density is finite (while maintaining the limit $M\ll N$).444In Sec. 4.2, we will see that finite temperature in the thermodynamic limit also opens a $2k_{F}$ gap in the dressed spinon energy. The spin wave excitations in the Yang-Gaudin Bose gas at fixed spinon density $M/L$ 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 $2k_{F}$ 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. (16) and (17). To begin, it is useful to consider the ground state with $J_{0}=-N/2$, corresponding to a spin rapidity of $\Lambda=-\infty$. Then, the Bethe equation (16) becomes

[TABLE]

where we define the shifted ‘integers’ $\tilde{I}_{j}^{(0)}=I_{j}^{(0)}+1/2$. The ground state is realized as the symmetric Fermi sea of the shifted integers $\tilde{I}^{(0)}=\left\{-\frac{N-1}{2},\ldots,\frac{N-1}{2}\right\}$, and the solution of Eq. (50) 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 $k\sim-2k_{F}$, which corresponds to solutions with large spin rapidity $\Lambda\gg 1$. As a result, the Bethe equation (16) can be approximated by

[TABLE]

The solutions to Eq. (51) are formally equivalent to the scenario where the right-most shifted integer of the ground state configuration $\tilde{I}^{(0)}$ is scattered across the Fermi sea: $\tilde{I}^{(0)}_{N}=\frac{N-1}{2}\to-\frac{N+1}{2}$. The momentum $-2k_{F}$ and energy gap $\Delta_{2k_{F}}$ thus coincide with that of the holon-antiholon $-2k_{F}$ excitation in the (one-component) Lieb-Liniger model (see also the following section) [5]–[7].

At a formal level we see that the $2k_{F}$ gap can be pictured in a similar manner to the $-2k_{F}$ 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 $\varrho_{s}=1/L$ and vary the total particle density $\varrho$. We see in Fig. 4 that the $2k_{F}$ gap $\Delta_{2k_{F}}$ is proportional to the total particle density $\varrho$, as can be expected from the strong coupling expansion (39). If instead we fix the total particle density and increase the spinon density $\varrho_{s}$ (by consider the “almost identical” spinon configuration discussed in the previous sections), we find similar: $\Delta_{2k_{F}}$ is proportional to the spinon density. This is also consistent with the strong coupling expansion for $M=1$ ($\varrho_{s}=1/L$), where $\varrho_{s}$ enters through the imbalance of the Fermi points, $q_{R}-q_{L}$. Combining these results, we see that the $2k_{F}$ 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 $N$.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 $m_{\ast}$ [69] (see Eq. (35))

[TABLE]

see also Eq. (37) and the accompanying discussion.

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 $I_{FS}$ and moving it to $I_{FS}+1$ (creating, effectively, a $h\bar{h}$ excitation, see the following section). As a result, the finite volume excitation gap of the (one-component) Lieb-Liniger model is

[TABLE]

as can be found from, e.g., Eq. (63). We note that the Yang-Gaudin model with $M=0$ 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. (53).

Let us now briefly highlight an important point. Setting $c=\infty$ for clarity, we have

[TABLE]

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. (54) 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 $\delta E_{LL}$ in the Yang-Gaudin Bose gas.

2.5 Extending beyond $2k_{F}$

We have considered the dispersion of excitations described by the configuration of integers (12), in which the momenta integers $\{I\}$ are described by the ground state configuration and the integer $J_{s}$ is varied between its bounds $|J_{s}|\leq N/2$. With fixed momenta integers $\{I\}$, the bounding of $J_{s}$ implies the spinon excitation has bounded momentum $0\leq|k|\leq 2k_{F}$. To realize spinon-like states with higher momenta, it is necessary to modify the configuration of the momenta integers $\{I\}$, creating multi-particle excitations. To extend to the momentum range $2k_{F}\leq|k|\leq 4k_{F}$, 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 $sh\bar{h}$ excitation. For $N=8$ particles, such a configuration of integers is shown below:

I$$-9$$-8$$-7$$-6$$-5$$-4$$-3$$-2$$-1[math]1$$2$$3$$4$$5$$6$$7$$8$$-$$+

J$$-4$$-3$$-2$$-1[math]1$$2$$3$$4$$J_{s}

In Fig. 5 we present the dispersion relation for the spinon ($0\leq|k|\leq 2k_{F}$) and the spinon-like ($2k_{F}<|k|\leq 4k_{F}$) states computed from the numerical solution of the Bethe ansatz equations (16), (17) for $N=100$ bosons on the length $L=50$ ring with interaction parameter $c=50$. 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 $4k_{F}$. 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 $\sim k^{2}$ 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 $E>2\Delta_{r}$ 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 $h\bar{h}$ and $sh$ excitations. In particular, we will focus on whether there exists a finite energy gap at $|k|=2k_{F}$ in the thermodynamic limit.

3.1 The holon-antiholon ($h\bar{h}$) continuum

The $h\bar{h}$ continuum constructed on top of the absolute ground state (11) 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 $\Lambda=-\infty$, we recover the Bethe ansatz equations for the Lieb-Liniger model and so the $h\bar{h}$ continuum above this state is also equivalent to that in the Lieb-Liniger model. We consider the problem of removing a particle with momentum $-q_{0}<k_{h}<q_{0}$ from within the Fermi sea and replacing it with a particle outside the Fermi sea $|k_{p}|>q_{0}$. The two-particle continuum of excitations is characterized in terms of the energy $E(k_{p},k_{h})$ and momentum $P(k_{p},k_{h})$ given by [7]

[TABLE]

Here we define $F_{h\bar{h}}(k|k_{p},k_{h})=(k_{j}-\bar{k}_{j})/(k_{j}-k_{j-1})$ the shift function defined with ground state momenta $k_{j}$ and the excited state (with the $h\bar{h}$ excitation) momenta $\bar{k}_{j}$. The shift function obeys the following integral equation [7]

[TABLE]

3.1.1 Strong coupling expansion.

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

[TABLE]

For a configuration of integers $\{I_{1},\ldots,I_{N}\}$, the strong coupling expansion for the momenta is

[TABLE]

which is consistent with our expression for the Fermi momentum $q_{0}$, Eq. (59).

As the shift function is functionally independent of $k$ to order $c^{-2}$, we can directly integrate Eqs. (55) and (56) to obtain

[TABLE]

Denoting the hole in the Fermi sea of integers by $-\frac{N-1}{2}\leq I_{h}\leq\frac{N-1}{2}$ and the excited integer by $I_{p}>\frac{N-1}{2}$, we have the simple relations

[TABLE]

This realises the usual holon-antiholon continuum of the Lieb-Liniger model [7], shown in Fig. 6. From Eq. (63), the sole effect of finite interaction strength [to order $O(c^{-3})$] is to rescale the energy compared to the $c=\infty$ limit; for $\varrho=2$, $c=50$ (as shown in Fig. 6) the energy is rescaled by $(1-2\varrho/c+4\varrho^{2}/c^{2})\approx 0.85$, e.g., $15\%$. As is well known in the Lieb-Liniger model [7], there is a gap at $k=2k_{F}$ due to the finite volume excitation gap – the lowest energy $2k_{F}$ excitation is umklapp-like: an integer is removed at the left Fermi point $-k_{F}$ and inserted immediately to the right of the Fermi point at $k_{F}$. The state is thus separated from the ground state by the finite volume excitation gap $\delta E_{LL}$ given in Eq. (53) (cf. Eq. (37)). With increasing system size $L$, there is no qualitative change in the $h\bar{h}$ continuum, unlike in the spinon dispersion (cf. Fig. 2). We note that the $2k_{F}$ gap in the $h\bar{h}$ continuum is comparable in size to the roton gap for the spinon excitations $\delta E_{LL}\approx\Delta_{r}$, cf. Eq. (53) and Eq. (39)

3.2 The spinon-holon ($sh$) continuum

We now turn our attention to the spinon-holon two-particle continuum. We compute the energy and momentum of the states with the configurations of integers described in Eqs. (14) along similar lines to the previous computation. The holon corresponds to a hole in the symmetric Fermi sea with momentum $k_{h}$; the position of the hole in the Fermi sea of integers is characterized by $0\leq j\leq N$. The total momentum of the $N$ particle state characterized by integers $j$ and $J$ is

[TABLE]

An example of the $sh$ excitation continuum for $N=100$ particles at density $\varrho=2$ is shown in Fig. 7. There, we have computed the continuum from the numerically exact solution of Eq. (3) and Eq. (4). As the spinon excitation is, in fact, a special case of the state under consideration (that with $j=N$), the lower bound of the continuum is given exactly by the spinon dispersion, as is expected, and hence the finite-size effects for the two are identical: the $2k_{F}$ gap in the $sh$ continuum vanishes in the thermodynamic limit with $M=1$.

It is worth reminding the reader that the $2k_{F}$ gap is not a result of the finite volume excitation gap: the gap is much larger than the finite volume excitation gap $\delta E_{YG}$ (37) (that is, many small momentum $k\sim 0$ spin wave excitations have energy below $\Delta_{2k_{F}}$). Instead, the non-zero excitation energy at $2k_{F}$ is a result of interactions between the spinons and holons/antiholons.

4 Bound states of spinons: $\Lambda$-strings

So far, we have focused on elementary excitations of the Yang-Gaudin Bose gas, the spinons and the holons. In the few cases where we have considered more than one spinon, $M>1$, we have been interested in solutions where both spin rapidities are real. However, this is not generically the case: in sectors with $M>1$ the spin rapidities can take complex values, describing multi-spinon bound states.

The complex solutions to the Bethe ansatz equations are organized in regular patterns in the complex plane, known as ‘strings’. Such solutions were already known by Bethe in his original solution of the Heisenberg model [1] and they have been well-studied in spin-$1/2$ quantum magnets (see, for example, Ref. [82] and the literature review of Ref. [83]), the Hubbard model [4], and the attractive Lieb-Liniger model [2]. They take the form

[TABLE]

where $\alpha$ labels the different possible $n$-spinon $\Lambda$-strings, whilst $a$ labels the rapidities contained within the string. The rapidities are distributed symmetrically about the real line, $\Lambda_{\alpha}^{n,a}=(\Lambda_{\alpha}^{n,n+1-a})^{*}$ [84], and have the real center $\omega^{n}_{\alpha}$. $\zeta_{\alpha}^{n,a}\in\mathbbm{C}$ is known as the string deviation – in the thermodynamic limit it is assumed to be exponentially small in the system size. The validity of the string hypothesis in the thermodynamic limit of the Yang-Gaudin Bose for computing thermodynamic observables has been confirmed [85], in keeping with general arguments about the structure of the Bethe ansatz equations [86]. We will consider this case in detail and comment on string deviations at the end of this section.

4.1 The Bethe-Takahashi equations

In the thermodynamic limit, the string hypothesis $\zeta_{\alpha}^{n,a}=0$ is believed to be valid; spin rapidities are organized in regular patterns in the complex plane described by

[TABLE]

which we characterize through their real centers $\lambda^{n}_{\alpha}$ and length $n$. Inserting (67) into the Bethe ansatz equations (3), (4), we arrive at the Bethe-Takahashi equations [85]

[TABLE]

where $\bar{I}_{j}\in(\mathbb{Z}+\frac{1}{2})$, $\bar{J}^{n}_{\alpha}$ are $N_{n}$ sets of $n$ numbers in $\mathbb{Z}$ ($\mathbb{Z}+\frac{1}{2}$) if $N_{n}$ is even (odd), and we define a convenient reparameterisation of the scattering phase

[TABLE]

4.2 The Yang-Yang-Takahashi equations

Now the thermodynamic limit can be taken. Introducing chemical potentials for each of the particle species ($\mu_{1}$, $\mu_{2}$) and finite temperature $T$ (we set $k_{B}=1$), one can derive the system of Yang-Yang-Takahashi (YYT) equations [87] for the dressed energy $\epsilon(k)$ and the length $n$ $\Lambda$-string dressed energies $\epsilon_{n}(\lambda)$. A detailed derivation can be found in Ref. [85]. The YYT equations read

[TABLE]

where $\mu=(\mu_{1}+\mu_{2})/2$, $\Omega=(\mu_{1}-\mu_{2})/2$, and we denote the convolution by

[TABLE]

We also define the functions

[TABLE]

The YYT equations (71)–(73) are supplemented by the asymptotic condition

[TABLE]

From this, the large momentum asymptotic form of these functions can be derived [85]

[TABLE]

For $n=1$ (i.e., the spinon) Eq. (77) simplifies to

[TABLE]

For the case with $SU(2)$ symmetry preserved, $\Omega=0$, we see that finite temperature $T$ opens a gap

[TABLE]

at large spin rapidity (e.g., $k\to 2k_{F}$). This is consistent with the picture presented in Sec. 2.3.3 that finite spinon density (which indeed is induced by finite temperature) leads to a $2k_{F}$ gap in the spinon dispersion.777 More generally, with $\Omega=0$ the length $n$ $\Lambda$-string has a gap of $\lim_{\lambda\to\infty}\epsilon_{n}(\lambda)|_{\Omega=0}=T\ln[n(2+n)].$

The Gibbs free energy per unit length $g$ can be expressed solely in terms of the dressed energy $\epsilon(k)$

[TABLE]

The density of particles $N/L$ and the spin flip density $M/L$ can be computed by taking derivatives of the Gibbs free energy $g$ with respect to the relevant chemical potentials

[TABLE]

A detailed study of the finite-temperature thermodynamics of the Yang-Gaudin Bose gas, including properties such as the polarization of the gas, the specific heat, and local pair correlations, can be found in Ref. [85].

4.3 Zero temperature limit of the YYT equations

Let us now consider the zero temperature limit of the YYT equations, on which the remainder of this work will be focused. We begin by defining the momentum $k_{F}$ for which the dressed energy changes sign $\epsilon(k_{F})=0$ (which may be zero in the case of vanishing particle density). Next, we use that the right hand sides of Eqs. (72) and (73) are positive, which implies that $\epsilon_{n}(\lambda)\geq 0$. It then follows that the YYT equations can be reduced to the form

[TABLE]

Notice that Eq. (82) has reduced to a Fredholm equation of the second kind, and so it may be solved numerically using standard routines, see [88] for one such example.

As is shown in B, the solutions to Eqs. (83) and (84) can be expressed solely in terms of the dressed energy $\epsilon(k)$. The length $n$ string dressed energies satisfy

[TABLE]

where we define $a_{0}(k)=\delta(k)$ and

[TABLE]

4.3.1 Strong coupling expansion in the zero temperature limit.

To serve as a check of our numerical solution of the YYT equations, we also consider Eqs. (82)–(84) in the strong coupling limit. As previously stated, the dressed energy $\epsilon(k)$ is given by a Fredholm equation of the second kind; it has the formal solution

[TABLE]

As this is expressed as a perturbative series, this provides a natural starting point for computing the corrections to the dressed energy $\epsilon(k)$. At zeroth order (e.g., $\pi c=\infty$) have the expected result

[TABLE]

where $k_{0}=\sqrt{\mu+\Omega}$ is the momentum at which the dispersion relation $\epsilon(k)|_{c=\infty}$ changes sign. Proceeding to order $1/(\pi c)$, we find

[TABLE]

At small momentum and with a finite particle density (e.g., $k/c\ll 1$ and $k_{F}/c\ll 1$) this simplifies to

[TABLE]

where we used $k_{F}=k_{0}+O(c^{-1})$. Equation (94) can also be obtained by a straight forward perturbative expansion of Eq. (82) by insertion of $a_{2}(x)=1/\pi c+O(c^{-3})$, valid provided $k_{F}/c\ll 1$. Notice that this also implies that the effective mass $m_{\rm eff}$ (cf. $\epsilon(k)=E_{0}+k^{2}/2m_{\rm eff}$) is unchanged to $O(c^{-2})$.

Let us now turn our attention to the length $n$-string dressed energies $\epsilon_{n}(\lambda)$ determined from Eqs. (83) and (84). Using

[TABLE]

we arrive at the zeroth order expansion

[TABLE]

That is, the string dressed energies are equal to their asymptotic values at zeroth order. The $c=\infty$ point is rather simple as the function $f(k)$ becomes infinitely wide, and hence there is no $\lambda$ dependence to the string energies $\epsilon_{n}(\lambda)$. Physically, this is reasonable: at $c=\infty$ spin flip excitations (e.g., spinons) cannot move through the gas as they are trapped by the hardcore interaction, experiencing Pauli exclusion. Indeed, this is seen as a divergence in the spinon effective mass as $c\to\infty$ [69, 70].

The next order correction to the length $n$-string dressed energy is considered in detail in C. For $\Lambda$-strings of length $n\geq 2$, we find a remarkably simple expression for the dressed energy

[TABLE]

whilst the spinon dressed energy is

[TABLE]

where $\psi(x)$ is the digamma function.

4.4 Numerical solution of the YYT equations

As we want to compute the dispersion relation (e.g., the dressed energy) of $\Lambda$-strings, we need to numerically solve the system of coupled integral equations (71–73). To do so, we follow one of the procedures outlined in Ref. [85]:

We begin by introducing a function cutoff $n_{\rm max}$, replacing $\epsilon_{n>n_{\rm max}}$ with their asymptotic value (76). 2. 2.

We introduce a momentum cutoff $\Delta_{i}$, reducing the range of the convolution integral (74) to $[-\Delta_{n},\Delta_{n}]$. This is quite reasonable as $f(k)\sim e^{-\pi k/c}$ approaches zero exponentially for large $k$, whilst it seems reasonable that $\epsilon(k)\sim k^{2}$ and the string dressed energies $\epsilon_{n}(\lambda)$ approaches a constant, cf. (77). 3. 3.

The solution is now estimated by the set of functions $\{\epsilon,\epsilon_{1},\ldots,\epsilon_{n_{\rm max}}\}$ where integrals are evaluated on the truncated range $[-\Delta_{i},\Delta_{i}]$ with $N_{i}$ evaluation points within this interval. 4. 4.

Starting from the non-interacting $c=0$ solution for $\epsilon(k)$ and the asymptotic values (77) for $\epsilon_{n}(k)$, the solution is then found iteratively. 5. 5.

We check convergence by varying $n_{\rm max}$, $\Delta_{i}$ and the discretization grid $N_{i}$.

Alternative numerical approaches, such as using a Möbius transformation to map the infinite line to the interval $(-1,1)$ (see, for example, Ref. [89]) are tricky to use here as the Kernel $f(k-q)$ becomes close to singular.

4.4.1 Results.

Having discussed both the strong coupling expansion, and our numerical approach to solving the YYT equations, we turn our attention to the results for the dressed energies $\epsilon(k)$ and $\epsilon_{n}(\lambda)$.

Let us begin with the dressed energy $\epsilon(k)$, from which many thermodynamic properties of the model follow, cf. Eq. (80). In Fig. 8 we present a comparison between the numerical solution of Eq. (82) and the strong coupling expansion (94), as well as the hardcore $c=\infty$ result (92), for the case with $\mu=2$, $\Omega=1$ and $c=10$. At small momenta we see that the dressed energy is significantly renormalized to lower energies from the hardcore limit, with the first order expansion (94) agreeing extremely well with the full numerical solution. For larger momenta (not shown in the plot) the dressed energy $\epsilon(k)$ remains parabolic and close to the $c=\infty$ result.

We see that the strong coupling result $\epsilon(k)=k^{2}-\tilde{k}_{F}^{2}+O(c^{-2})$ (where $\tilde{k}_{F}=k_{0}[1+2k_{0}/(3\pi c)]\approx k_{F}$) is a very good approximation for the whole range $k\in[-k_{F},k_{F}]$. As a result, the zero-temperature Gibbs free energy density (80) can be evaluated as

[TABLE]

Hence we can evaluate, within the $1/c$ expansion, the number density for the different species

[TABLE]

For the parameters in Fig. 8 this gives $N/L=0.6324+O(c^{-2})$, which is in excellent agreement with the numerically exact result $N/L=0.6326$.

Now, let us turn our attention to the dressed energies $\epsilon_{n}(\lambda)$, corresponding to the length $n$ $\Lambda$-strings with real center $\lambda$. In Fig. 9 we present the dispersion relations as a function of the real string center $\lambda$ for the length $n=1$ and $n=2$ $\Lambda$-strings, whilst in Fig. 10 we present similar for $n=3-6$. In all cases, we consider the same parameter set as previously, see Fig. 8. The dispersion relations obey the asymptotic condition $\lim_{\lambda\to 0}\epsilon_{n}(\lambda)=2n\Omega$,888We note that the spinon dispersion in Sec. 2 is computed for a single spinon in a finite system. Taking the thermodynamic limit $L\to\infty$ with $M=1$ fixed results in $\Omega=0$. Then we have $\lim_{\lambda\to\infty}\epsilon_{1}(\lambda)\to 0$, which is consistent with the previously presented results in finite-size systems. and they are Lorentzian in shape close to $\lambda=0$ with amplitude decaying as $1/n$. This is consistent with the strong coupling expansion (97), although we note that Eq. (97) significantly underestimates the amplitude of Lorentzian for the data presented in Figs. 9 and (10), which is not too surprising as the strong coupling expansion parameter $\bar{k}_{F}=\pi k_{F}/c>0.5$ is sizeable.

4.5 A single 2-string in a finite system

So far we have studied the thermodynamic limit with finite density of both particle species, which corresponds to a finite density of string excitations of all lengths [cf. Eqs. (71-73)]. We now turn our attention to a slightly different scenario, closely related to the cases consider in Secs. 2 and 3: a state containing a single $2$-string on top of the ground state configuration of integers $\{I_{j}\}$. We begin by considering small finite systems, in which one needs to take string deviations into account, before considering the large $L$ limit via a strong coupling expansion (see, for example, Secs. 2 and 3).

4.5.1 Small systems: string deviations and finite-size effects.

In a finite size system, the string solutions (67) are deviated (or deformed) away from their perfect ordering $\Lambda_{\alpha}^{n,a}\to\Lambda_{\alpha}^{n,a}+i\zeta_{\alpha}^{n,a}$, where $\zeta_{\alpha}^{n,a}\in\mathbb{C}$. These so-called string deviations have been well studied in integrable lattice models, such as the Heisenberg chain [2, 83, 90], the Hubbard model [4, 91], and the spin-1 Babujan-Takhtajan chain [92]. In continuum integrable models, they have been extensively studied in the attractive limit of the Lieb-Liniger model [93]. In all of these cases, the string deviations typically vanish exponentially in the system size $L$, and this is expected to be the case in the Yang-Gaudin Bose gas as well [85].

In this section we will numerically solve the Bethe ansatz equations for small numbers of particles and spinons without assuming the string hypothesis. We limit ourselves to small numbers of particles as this is numerically challenging: we need to solve a system of coupled complex non-linear algebraic equations (there are no good general methods for doing this). In particular, we will focus on the case of the length $n=2$ $\Lambda$-string.

Our aim is to solve the Bethe ansatz equations for $N$ particles with the spin rapidities describing a deviated length $2$-string, parameterized by

[TABLE]

where $\lambda,\delta$ are real parameters characterizing the string center and the deviation in the complex plane, respectively.

We proceed along similar lines to Refs. [83, 92], by first carefully treating the branch cuts of the arctangent. Working with the branch cut of the logarithm defined by $-\pi<{\rm Im}\ln z\leq\pi$, we have

[TABLE]

Taking the sum of the two arctangents with complex conjugate arguments, we have [83, 92]

[TABLE]

with

[TABLE]

Here $\Theta(x)$ is the Heaviside function with $\Theta(0)=1$, ${\rm sgn}(a)$ is the signum function with ${\rm sgn}(0)=0$, and we note that the function has the limit $\lim_{b\to 0}\xi(a,b)={\rm sgn}(a)\pi/2$.

With these definitions at hand, we can derive the allowed integers $J_{1},J_{2}$ that described the $2$-string.999When the string hypothesis is assumed, a $2$-string is described by a single integer ($\bar{J}^{n}_{\alpha}$ in Eq. (69)) which fixes the real string center $\lambda$. In a finite system, two integers $J_{1},\,J_{2}$ enter the Bethe equations (reflecting that this is a state with two flipped spins) and the string is parameterized in terms of the real string center $\lambda$ and the string deviation $\delta$. We will see that the two integers, $J_{1},\,J_{2}$ are not independent and, indeed, if $\delta=0$ are fixed to be the same. Inserting (102) into the Bethe ansatz equation (4) and then taking the difference between the two equations with integers $J_{1},\ J_{2}$ we have

[TABLE]

Matching the real part of the left- and right-hand sides and using ${\rm Re}[\arctan(z)-\arctan(z^{*})]=0$ if ${\rm Re}(z)\neq 0$, we find the following condition on the integers $J_{1},J_{2}$

[TABLE]

We see that there are two types of string deviations: those which narrow the string in the complex plane $\delta<0$ (described by $J_{1}=J_{2}$), and those which expand it $\delta>0$ (described by $J_{2}=J_{1}+1$). Analogous behavior is seen in other integrable models, see for example Refs. [83, 92]. It is worth noting that once the spin rapidities become complex, we lose the “Pauli-like exclusion principle” for the integers $J_{1},J_{2}$.

Now that we understand the allowed integers $J_{1},J_{2}$ for the spin rapidities, we solve the Bethe ansatz equations obtained by inserting (102) into Eqs. (3) and (4)

[TABLE]

The final equation is equivalent to (106) (it is obtained instead from the non-logarithmic Bethe ansatz equations) and is used for numerical convenience [92]. We compute the energy and momentum of the $n=2$ $\Lambda$-string relative to the absolute ground state according to

[TABLE]

where $E_{0}$ ($P_{0}$) is the ground state energy (momentum).

4.5.2 Numerical results.

We now proceed to solve Eqs. (110)–(112) for the case with the ground state configuration of integers $I_{j}$ (11). We vary the integers $J_{1},\ J_{2}$ [subject to the conditions (109)] to obtain the set of states containing a single $2$-string and we compute the dispersion relation through Eqs. (113). Figure 11(a) shows the dispersion relation as a function of the real string center $\lambda$ for $|\lambda|\leq 1$ (corresponding to $|P_{2s}|\lesssim k_{F}$), whilst Fig. 11(b) presents the dispersion relation as a function of the $2$-string momentum, $P_{2s}$.

The structure of the dispersion relation shows two distinct regions of behavior separated by a ‘kink’ at $|P_{2s}|\sim k_{F}$. This can easily be understood in the following manner. Consider the $2$-string with real string center $\lambda$ in the set of momenta $\{k\}_{j}$ which is bounded by $k_{R,L}\approx\pm k_{F}$. In the large system limit, we assume that string deviations are small (i.e., $c\delta\ll|k_{j}-\lambda|$) and expand the scattering phase between the string center and the momenta as

[TABLE]

Inserting (114) into (110), we see that the real string center causes an effective shift of the integers $I_{j}$ according to

[TABLE]

In short, the real string center causes all integers defining momenta to its left to shift to the left, while all those to the right of it shift to the right. This is very much like introducing a hole into the sea of momenta. The kink in the dispersion occurs at $P_{2s}=k_{R,L}\approx\pm k_{F}$, where the real string center exits/enters the Fermi sea of momenta. Accordingly, in the intermediate region $|P_{2s}|\lesssim k_{F}$) the dispersion in Fig. 11 can be pictured as describing an excitation similar to a holon coupled to a $2$-string.

We see that there are significant finite size effects in the $2$-string dispersion for small numbers of particle $N=8-16$. Qualitatively, we see that there are two different regions: firstly, for $|P_{2s}|<k_{F}$ the dispersion converges rapidly with increasing particle number. Secondly, the outer region $|P_{2s}|>k_{F}$ exhibits an upturn that slowly diminishes with increasing system size. Finite size scaling of the $2$-string excitation energies, see Fig. 12, is consistent with $\lim_{L\to\infty}E_{2s}(2k_{F})=0$, as expected in the thermodynamic limit with $M=2$ (cf. Eq. (77) with $\Omega=0$). The finite size scaling also shows that the excitation with momentum $P_{2s}=\pm k_{F}$ remains of finite energy with $L\to\infty$.

4.5.3 Strong coupling expansion.

To confirm the structure of dispersion relation for a single 2-string observed in the previous section, we perform a strong coupling expansion along similar lines to Secs. 2 and 3. In this section we will consider large $L$, adopting the string hypothesis as a working assumption (e.g., we set $\delta\to 0$). The calculation proceeds along similar lines to the spinon case (cf. Sec. 2), resulting in the following equations

[TABLE]

where $\Phi(u,v)=\phi_{1}(u,v)+\pi{\rm sgn}(u-v)$, $F_{2s}(k|\lambda)$ is the shift function that characterizes how the presence of the $2$-string with real center $\lambda$ modifies the distribution of momenta $\{k_{j}\}$ from its ground state configuration $\{k^{(0)}_{j}\}$, and $J_{1,2}\in\mathbb{Z}+\frac{1}{2}$ with $J_{2}-J_{1}=0,1$ as discussed in the previous section. $q_{R,L}$ denotes the right/left edges of the momenta distribution, which can be determined from the relations (30) with $P=-2\pi(J_{1}+J_{2})/L$, which follows from Eq. (9) with $\sum_{j}I_{j}=0$.

The energy and momentum of the $2$-string state with real center $\lambda$ are given by [cf. Eqs. (29)]

[TABLE]

We present the $c=\infty$ result in Fig. 13, comparing to the exact result for the case of $N=16$ particles at unit density with $c=50$. We see excellent agreement between the strong coupling result and the exact result, further confirming the unusual structure of the $2$-string dispersion. This also implies that corrections to the dispersion due to string deviations are rather small, as we take the $\delta\to 0$ limit in the strong coupling expansion. We also present $c=\infty$ results for two larger systems, which suggests that in the $N\to\infty$ limit, the $2$-string has a flat (or very close to flat) dispersion when the real string center is located outside the Fermi sea of momenta. This is not entirely surprising: the mass of a single spinon, $m_{*}$, also diverges in the $c\to\infty$ limit [69, 70] (see also Sec. 2).

5 Discussion

In this work, we have studied the excitation spectrum of the Yang-Gaudin gas, with a particular focus on how finite system size $L$ affects the dispersion relation for elementary excitations. Using both a strong coupling $1/c$ expansion and exact numerical solution of the Bethe ansatz equations, we have shown that there are prominent finite-size effects in this system, with qualitative and quantitative features changing at energy scales much greater than the finite volume excitation gap. The origin of this is the drastic change in the finite volume excitation gap between the one-component and two-component cases, arising from the new low-energy spin degrees of freedom.

We started by considering the spinon (spin wave) dispersion above the fully polarized ground state of the two-component gas. It was seen that the spinon dispersion is non-monotonic: at small momentum $k\ll\varrho$, conventional spin wave excitations $E(k)\sim k^{2}/2m_{*}$ with effective mass $m_{*}$ are observed [69, 70], whilst at large momentum the energy of the spinon drops, and a roton-like minima occurs close to momentum $|k|=2k_{F}$. This minima has finite energy, denoted the roton gap $\Delta_{r}$, which decreases with increasing system size at fixed density $\Delta_{r}\sim 1/L$ for a single spinon $M=1$. For reasonably large systems $L=50-200$, the roton gap is many times greater than the finite volume excitation gap $\Delta_{r}\gg\delta E_{YG}$. The spinon dispersion also defines the lower bound of the spinon-holon continuum, and accordingly there is a pronounced gap at $k=2k_{F}$ gap which vanishes as $1/L$ in the system size.

The origin of the roton-like minima in the spin wave dispersion is the finite interactions between the excitations in the Yang-Gaudin Bose gas. The $2k_{F}$ gap varies as the product of the spinon density and the particle density. As a result, taking the thermodynamic limit with fixed spinon density ($N$, $M$, $L\to\infty$ with $N/L$, $M/L$ fixed) results in a spinon dispersion with a roton-like minima and $2k_{F}$ gap. Numerical evidence for this was presented in Fig. 3. The Yang-Gaudin Bose gas thus constitutes an example of roton-like excitations occurring in an exactly solvable microscopic model, a problem which has received attention for over 75 years [75]–[78].

The holon-antiholon continuum coincides with that of the (one-component) Lieb-Liniger model. As in that case, there is an energy gap at momentum $2k_{F}$ which arises from the finite volume excitation gap of the one-component gas, $\delta E_{LL}$, and hence this gap similarly vanishes as $1/L$. However, in the context of the two-component Yang-Gaudin model at strong coupling, we should bear in mind that the energy scale $\delta E_{LL}$ is much greater than the finite volume excitation gap $\delta E_{YG}$ set by the spinon excitations. In other words, unlike the (one-component) Lieb-Liniger model where $\delta E_{LL}$ is a very small energy scale, the $2k_{F}$ ‘gap’ in the holon-anitholon continuum is actually a large energy scale in the Yang-Gaudin gas (comparable to the roton gap $\Delta_{r}$) and there are many spinon excitations with energy less than $\delta E_{LL}$, see Eq. (54).

What are the implications of these finite-size effects in the two-component Bose gas? Consider first the finite-temperature properties of the model in a finite-size system. At low temperatures $T\ll\Delta_{r}$ the system is unaware of the gapped roton-like excitations and hence the low temperature properties are uninfluenced. Raising the temperature, at $T\sim\Delta_{r}$ the system starts to sample the roton-like minimum and the additional gapped excitations can lead to properties of the system changing – new excitations with non-zero average momentum are activated and the entropic properties of the system change. It is tempting to suggest that non-monotonicities in physical observables (such as the local density-density correlator) with increasing temperature observed in Ref. [94] may be due to such excitations becoming thermally activated, but this requires further investigation.

The roton-like minima and anomalously large $2k_{F}$ gaps in excitation continua may also have profound influences on the non-equilibrium dynamics of the system. Consider the spinon dispersion – there are now two dominant velocities for excitations: the traditional ($|k|<\varrho$) spin wave velocity $v_{s}$ and the velocity associated with excitations about the roton-like minima $v_{r}$. These are generically different in finite-size systems, and their behavior when varying the interaction strength $c$ and the density $N/L$ is different. Accordingly, if a non-equilibrium initial state projects on to the low-lying spinon excitations, one may directly observe (via imaging the local density) the propagation of both types of excitations [56]. This may provide a root to observing these excitations in experiments on cold atomic gases (see, e.g., Refs. [95, 96]). The role such roton-like excitations play in the dynamics of a distinguishable impurity immersed in the Bose gas (see, for example, Ref. [55], where unusual dynamics are observed) is yet to be understood.

In the second part of work, our attention was turned to the bound state excitations, so-called $\Lambda$-strings, in the Yang-Gaudin Bose gas. These are analogous to well-known excitations in the Heisenberg XXZ model [2] and the Hubbard model [4]. Focusing on the thermodynamic limit, we solved the Yang-Yang-Takahashi equations numerically and compared results to the strong coupling expansion. We found that the dispersion of $\Lambda$-strings of length $n$, $\epsilon_{n}(\lambda)$ can be well approximated by a Lorentzian for a wide range of real string centers $\lambda$ about $\lambda=0$. In the zero temperature limit, this can be understood from the simple relationship (86) between the string energies and the dressed energy $\epsilon(k)$. The relation (86) can be solved under a strong coupling expansion (97), yielding a Lorentzian dressed energy for length $n\geq 2$ $\Lambda$-strings.

We finished by commenting on finite-size effects for the length $n=2$ $\Lambda$-string. To do so, we numerically solved the Bethe ansatz equations for a single $2$-string by treating the string deviations (which are neglected in the thermodynamic limit under the string hypothesis). As with the other single- and multi-particle excitations, we found that there are significant finite size effects in the dispersion of the $2$-strings when system sizes are small; parts of the dispersion relation (those corresponding to the real string center being located outside of the Fermi sea of momenta) converge slowly with increasing system size. Comparison with the strong coupling expansion suggests that corrections to the 2-string dispersion due to the string deviations are very small.

Bound states in integrable models are generally difficult to observe with equilibrium probes, such as inelastic neutron scattering in spin chain materials. However, recently it was realized that bound states can be visible in the non-equilibrium dynamics arising from a local quantum quench [97]. The non-equilibrium setting can also allow access to phases which would be thermodynamically unstable in equilibrium, where bound states dominate the physics [89, 98]. Whilst there have been some studies of local quenches in the two-component Bose gas [55, 56], they were in a limit where bound states are absent. Attacking non-equilibrium problems with bound states in multi-component integrable quantum gases remains an outstanding problem.

We are grateful to Jean-Sébastien Caux, Tamás Pálmai, and Eoin Quinn for useful conversations, and Alexei Tsvelik for enthusiastic discussions and encouragement. This work was supported by the Condensed Matter Physics and Materials Science Division, Brookhaven National Laboratory, in turn funded by the U.S. Department of Energy, Office of Basic Energy Sciences, under Contract DE-SC0012704. We acknowledge the hospitality and support of the program “Mathematical Aspects of Quantum Integrable Models in and out of Equilibrium” at the Isaac Newton Institute for Mathematical Sciences, University of Cambridge (under EPSRC grant EP/K032208/1), where this work was initiated.

Appendix A The continuum limit

Let us briefly consider taking the continuum limit ($L\to\infty$, $N\to\infty$ with $N/L$ fixed) of Eq. (19). There, we have the momenta $\tilde{k}_{j}$ for a finite momentum state (see the discussion preceding Eq. (19)) and the momenta $k_{j}$ of the ground state (forming a symmetric Fermi sea with zero total momentum). To take the continuum limit, we use the following identities

[TABLE]

where $q_{L}$ ($q_{R}$) is the momentum of the left (right) Fermi point for the configuration of integers $\tilde{I}_{0}$ and $\tilde{J}>-N/2$ (e.g., the finite momentum state), whilst $q_{0}$ is the Fermi momentum of the ground state. We have also introduced the root distributions

[TABLE]

which are identical to leading order in $L$. Notice that this means in the $L\to\infty$ limit, the finite momentum of the state with $\tilde{k}_{j}$ is realized through the non-symmetric bounds of the integral $q_{L}$, $q_{R}$, which implicitly depend upon $\tilde{I}_{0}$ and $\tilde{J}$.

Appendix B $\Lambda$-string dressed energies at zero temperature

Let us now consider the string dressed energies which satisfy Eqs. (83) and (84). The integral featuring the dressed energy on the right hand side of (83) can be considered as a known function (as the dressed energy can be determined using standard methods for Fredholm equations of the second kind [88]). Let us denote this function by

[TABLE]

Motivated by the limit (85) and the zeroth order contribution to the strong coupling expansion (96), we write the string dressed energy in the form

[TABLE]

where $\lim_{n\to\infty}\bar{\epsilon}_{n}(\lambda)/n=0$. The YYT equations (83), (84) now read

[TABLE]

We proceed by applying a Fourier transform ${\cal F}$ to the above system of equations. The convolution theorem allows us to decompose the Fourier transform of the convolution into the product of the Fourier transforms

[TABLE]

where we define the Fourier transform as

[TABLE]

The YYT equations are now a set of functional equations for the Fourier transform of the dressed energies ${\cal F}(\bar{\epsilon}_{n})=\tilde{\epsilon}_{n}(\omega)$:

[TABLE]

Here $\tilde{f}(\omega)=1/[2\cosh(\pi c\omega)]$. We next define the vector $\vec{\varepsilon}(\omega)=(\tilde{\epsilon}_{1}(\omega),\tilde{\epsilon}_{2}(\omega),\ldots)^{\mathrm{T}}$ this set of equations can be written as

[TABLE]

Next we truncate the infinite hierarchy of equations at some level $n=n_{\rm max}\gg 1$, and we replace $\epsilon_{n_{\rm max}+1}$ with its asymptotic value. This leads to an additional contribution to the matrix on the RHS of Eq. (135) in the lower right corner. Rewriting Eq. (135), we have

[TABLE]

where $\tilde{g}(\omega)=\tilde{h}(\omega)/\tilde{f}(\omega)$ and $\tilde{d}(\omega)=-1/\tilde{f}(\omega)$.

The solution of the set of equations (126) and (127) is now obtained by inverting the matrix on the LHS of Eq. (145)

[TABLE]

Due to the particularly simple structure of the matrix being inverted, the matrix elements of $R$ are known and have the simple form [99]

[TABLE]

Ignoring the contribution from the lower right element on the right hand side (which will be justified as we will take $n_{\rm max}\to\infty$), we have

[TABLE]

We can now take the $n_{\rm max}\to\infty$ limit, which gives

[TABLE]

We now find the dressed energy for the length $n$ string by application of the inverse Fourier transform ${\cal F}^{-1}$ and remembering that $\epsilon_{n}(\lambda)=2n\Omega+\bar{\epsilon}_{n}(\lambda)$

[TABLE]

that is

[TABLE]

Substituting in the Fourier transform $\tilde{h}(q)$ gives

[TABLE]

For $n>1$, computing the $\omega$ integral gives

[TABLE]

and hence we have

[TABLE]

where $a_{n}(k)$ is defined in Eq. (75). We see that this expression satisfies the limit conditions $\lim_{n\to\infty}\epsilon_{n}(\lambda)/n=2\Omega$ and $\lim_{\lambda\to\infty}\epsilon_{n}(\lambda)=2n\Omega$.

For $n=1$ we have

[TABLE]

and hence

[TABLE]

Appendix C Strong coupling expansion for the string dressed energies

We begin with Eq. (86), which was derived in the previous appendix (162). The strong coupling solution of this equation is a little more involved, as one cannot take the naive strong coupling expansion of $f(k)$ (due to it featuring in the convolution, so its argument can become arbitrarily large). Instead, we compute $h(q)$ using the strong coupling expression for the dressed energy $\epsilon(k)$ without further approximation. That is, we compute

[TABLE]

where we use the strong coupling expression for $k_{F}$. The result is

[TABLE]

Here we use the shorthand notation $\bar{q}=\pi q/c$, $\bar{k}_{F}=\pi k_{F}/c$ and we define the polylogarithm function

[TABLE]

At this point, we perform an expansion in $\bar{k}_{F}\ll 1$, which is reasonable for large $c$ and finite particle density. Taking care with the factors of $c/\pi$ which appear in the above expression and expanding the accompanying expression to appropriate orders in $\bar{k}_{F}$, we eventually arrive at the simple expression

[TABLE]

Now that we have a simple form for $h(q)$, we can compute the convolution

[TABLE]

using the convolution theorem. Fourier transforming the two functions gives

[TABLE]

Multiplying the two Fourier transforms and subsequently taking the inverse Fourier transform yields the convolution

[TABLE]

where $\psi(x)$ is the digamma function. Now, there is a rather beautiful simplification when we consider the full term appearing in Eq. (86):

[TABLE]

which follows from the property of the digamma function $\psi(1+z)=\psi(z)+1/z$. So, we arrive at the simple expression

[TABLE]

For the special case of $n=1$, this simplification does not occur and we have, explicitly,

[TABLE]