Integrability and weak diffraction in a two-particle Bose-Hubbard model

# Integrability and weak diffraction in a two-particle Bose-Hubbard model

Daniel Braak Experimental Physics VI, Center for Electronic Correlations and Magnetism, University of Augsburg, 86135 Augsburg, Germany    J. M. Zhang Theoretical Physics III, Center for Electronic Correlations and Magnetism, University of Augsburg, 86135 Augsburg, Germany Max Planck Institute for the Physics of Complex Systems, Nöthnitzer Str. 38, 01187 Dresden, Germany    Marcus Kollar Theoretical Physics III, Center for Electronic Correlations and Magnetism, University of Augsburg, 86135 Augsburg, Germany
###### Abstract

A recently introduced one-dimensional two-particle Bose-Hubbard model with a single impurity [J. M. Zhang et al., Phys. Rev. Lett. 109, 116405 (2012)] is studied on finite lattices. The model possesses a discrete reflection symmetry and we demonstrate that all eigenstates odd under this symmetry can be obtained with a generalized Bethe ansatz if periodic boundary conditions are imposed. Furthermore, we provide numerical evidence that this holds true for open boundary conditions as well. The model exhibits backscattering at the impurity site—which usually destroys integrability—yet there exists an integrable subspace. We investigate the non-integrable even sector numerically and find a class of states which have almost the Bethe ansatz form. These weakly diffractive states correspond to a weak violation of the non-local Yang-Baxter relation which is satisfied in the odd sector. We bring up a method based on the Prony algorithm to check whether a numerically obtained wave function is in the Bethe form or not, and if so, to extract parameters from it. This technique is applicable to a wide variety of other lattice models.

###### pacs:
03.65.Ge, 03.75.-b, 71.10.Fd

## I Introduction

Exactly solvable models are always appreciated, not only because of the beauty they embody but also because of the pivotal role they can play. Among the exactly solvable models known so far, many are solved by the Bethe ansatz for the wave function bethe (). A prerequisite for the consistency of this ansatz is the absence of diffraction sutherland (), which was first explicitly pointed out by McGuire mcguire (), and is formally characterized by the Yang-Baxter equation yang (); baxter ().

Recently, in an investigation motivated by atomic physics, the present authors discovered that a two-particle lattice model, whose continuum counterpart was previously claimed to be diffractive by McGuire mcguire (), has eigenstates in the Bethe ansatz form although they do not span the whole Hilbert space prl2012 (); pra2013 (). The model is very simple—it consists of two interacting identical bosons moving on a one-dimensional infinite lattice with a single site defect. The model possesses a reflection symmetry (parity) and the states odd under parity conform to the Bethe ansatz. Its original purpose was to investigate how the competition or cooperation between the interaction and the defect potential will affect the formation of bound states in the system, as interesting consequences are known in a similar object, i.e., the negative hydrogen ion H rau01 (). It turns out that the model can exhibit two kinds of exotic bound state, i.e., the bound state in the continuum wigner () and the bound state at threshold zero (); mattis (); nieto ().

It was convenient to consider the model on the infinite line because the bound state can be unambiguously discerned from the extended states belonging to the continuous spectrum. The Yang-Baxter relation is then sufficient to prove the validity of the Bethe ansatz in the odd subspace pra2013 (). However, the situation is different on a finite lattice zvyagin (). The boundary conditions impose now restrictions on the allowed momenta and as the impurity potential generates backscattered waves which in all known cases inevitably spoil integrability, it was not clear whether meaningful Bethe ansatz equations can be derived to determine the discrete spectrum and if so, whether their solutions span the full subspace with odd parity.

We shall address these questions in the following. Both can be answered in the affirmative for odd lattice size and periodic boundary conditions. However, the framework of the Bethe ansatz has to be generalized to incorporate the non-trivial backscattering at the impurity. The standard approach (see e.g. andrei ()) defines the -matrix as a relation between local amplitudes, whereas in our case the amplitudes have to be defined in a non-local way.

This paper is organized as follows. In Sec. II we define the model and outline the formalism of the generalized Bethe ansatz. We show that it allows for a natural inclusion of the reflection (parity) symmetry and thus explains at the same time the integrability of the odd sector and the failure of the ansatz in the even sector. The representation of periodic boundary conditions in the formalism is discussed in some detail. We derive the Bethe ansatz equations for this case. The proof that their solutions span indeed the whole space is given in the appendix.

The partial solvability of this model contradicts common lore about the possibility to use the Bethe ansatz for systems with impurities which are neither located at the edges sklyanin () nor have fine-tuned features which effectively remove the backscattering anjo (). In the present case, our primary evidence comes from numerical data obtained by exact diagonalization. This is the subject of Sec. III. Here we study both open and periodic boundary conditions with three different methods to corroborate the fact that the odd sector is always integrable while the even sector is not. In particular, the Bethe ansatz checking method in Sec. III.2 allows us to confirm that all the eigenstates in the odd subspace have the Bethe form. On the other hand, the even sector exhibits eigenstates which are composed of only a few Fourier components and thus resemble Bethe states. These “weakly diffractive” states are presented in Sec. III.3. We show in Sec. IV that some of them correspond to a weak violation of the generalized Yang-Baxter equation, while the ordinary Yang-Baxter equation is strongly violated. Sec. V summarizes the results and lays out some directions for future work.

## Ii The model on a finite lattice

The model describes two spinless interacting bosons hopping on a lattice with a single site defect. On a one-dimensional lattice with an odd number sites, the Hamiltonian reads

 ^H=M′∑i=−M′[−(^a†i^ai+1+h.c.)+U2^a†i^a†i^ai^ai]+V^a†0^a0. (1)

We measure energy in units of the hopping integral; the parameters and are the strengths of on-site interaction and the defect potential at site 0, respectively. A crucial property of (1) is its invariance under the reflection .

The Hamiltonian (1) conserves particle number and we work in the two-particle subspace. To set up the Bethe ansatz, we use the first quantized form of (1); the wave function reads with being two integers, and Bose symmetry requires . The Hamiltonian acts on as

 ^Hψ(x1,x2) =∑Δ=±1−[ψ(x1+Δ,x2)+ψ(x1,x2+Δ)] +[V(δx1,0+δx2,0)+Uδx1,x2]ψ(x1,x2). (2)

The group generated by is and a state will be called even [odd] under parity if it satisfies [].

### ii.1 Generalized Bethe ansatz

#### ii.1.1 Non-local Yang-Baxter equation

The one-dimensional elastic scattering between the particles () leaves the momenta pairwise invariant but the impurity potential generates reflected waves with opposite momentum. Therefore only the absolute values of the momenta could serve as good quantum numbers. We try to generalize the Bethe ansatz, accounting for this simple type of diffraction through the introduction of an internal quantum number for each particle which denotes the sign of the corresponding momentum. In this way the scattering phases in a system of spinless bosons are replaced by non-commutative -matrices. The generalized ansatz reads,

 ψ(x1,x2)=S∑R,σ1,σ2ARσ1σ2ei(k1σ1x1+k2σ2x2)χR(x1,x2). (3)

The wavefunction is composed from plane waves defined in each region with denoting the set in the - plane with

 [ijk]↔xi≤xj≤xk (4)

and . is the characteristic function of region and symmetrizes the wave function with respect to and . The pseudo-spin index accounts for the possible backscattering at the impurity site. Due to the interactions at the boundaries of the regions , the amplitudes are in general different in the six regions given by (4). For example, the amplitudes and are related by the -dependent scattering phase (, )

 A[012]σ1σ2=σ1s1−σ2s2+iU2σ1s1−σ2s2−iU2A[021]σ1σ2. (5)

The scattering phases (5) are unimodular for real momenta . On the other hand, the amplitudes are related to as

 (6)

The matrix appearing in (6) maps the amplitudes with particle 1 on the right of the impurity to the amplitudes corresponding to particle 1 being on the left of the impurity. This entails that it is not unitary but preserves instead the particle current across the impurity:

 |A[102]++|2−|A[102]−+|2=|A[012]++|2−|A[012]−+|2. (7)

This non-unitarity is the main obstacle to derive the Bethe ansatz equations even if the ansatz (3) should be consistent - which is not the case in general, the reason being the partial transmission/reflection generated by the impurity as exemplified in (6). It turns out that the parity symmetry of the model is able to restore integrability to a certain degree in the two-particle sector. To see that, it is mandatory to define the amplitude vectors in such a way that all -matrices become unitary.

All scattering processes are necessarily unitary if considered not as a transfer from “right to left” as in (6) but from “past to future”: One relates the set to the set , i.e. the incoming to the outgoing waves,

 (8)

The matrix in (8) is unitary for real as it should be. However, it collects amplitutes belonging to different regions into a single vector, in contrast to (6). A consistent formalism has to express the scattering between particle 1 and 2 with non-local amplitudes as well which makes it convenient to split the regions and into two subregions (see Fig. 1),

 [102A]↔x1≤0≤x2, |x1|≤|x2|, ↔x1≤0≤x2, |x1|≥|x2|, ↔x2≤0≤x1, |x1|≤|x2|, (9) ↔x2≤0≤x1, |x1|≥|x2|.

The interaction part of the Hamiltonian vanishes on the boundary between and , so the corresponding amplitudes have to be the same. This will be guaranteed if the -matrices are chosen correctly. We collect the local amplitudes into the following eight four-vectors,

 A1=(A[012]−−,A[201A]−+,A[102A]+−,A[210]++)T, A2=(A[012]+−,A[201A]++,A[102A]−−,A[210]−+)T, A3=(A[021]+−,A[201B]++,A[102B]−−,A[120]−+)T, A4=(A[021]++,A[201B]+−,A[102B]−+,A[120]−−)T, A5=(A[012]++,A[201A]+−,A[102A]−+,A[210]−−)T, (10) A6=(A[012]−+,A[201A]−−,A[102A]++,A[210]+−)T, A7=(A[021]−+,A[201B]−−,A[102B]++,A[120]+−)T, A8=(A[021]−−,A[201B]−+,A[102B]+−,A[120]++)T.

The wavefunction reads in terms of the components , (, ),

 ψ (x1,x2)= S{ (∑k[Ak8χAk8+Ak1χAk1])e−i(k1|x1|+k2|x2|) + (∑k[Ak2χAk2+Ak3χAk3])eik1|x1|−k2|x2| + (∑k[Ak4χAk4+Ak5χAk5])ei(k1|x1|+k2|x2|) (11) + (∑k[Ak6χAk6+Ak7χAk7])e−ik1|x1|+k2|x2|},

where denotes the characteristic function of the region belonging to . Note that the exponential factor is invariant under and the same for all components of . Moreover, we have for all the relation

 ^P[χAkj](x1,x2)=χA5−kj(x1,x2). (12)

It follows that

 ^P[ψ](x1,x2;{Akj})=ψ(x1,x2;{A5−kj}), (13)

which means that acts on by transforming each amplitude vector with the matrix , . The projection of onto the even (odd) subspace is done by projecting the amplitudes onto the two-dimensional eigenspace of with eigenvalue (). The -matrices for connect “adjacent” amplitudes ( because is adjacent to ). correspond to an exchange of particle 1 with 2 and [] to scattering of particle 1 [2] at the impurity, see Fig. 2.

The determination of the -matrices from (II) is straightforward pra2013 (). Defining

 α=s1−s2+iU2s1−s2−iU2,β=s1+s2+iU2s1+s2−iU2, (14)

we find e.g. for ,

 ^S(1,8)=⎛⎜ ⎜ ⎜⎝α−100001000010000α−1⎞⎟ ⎟ ⎟⎠, (15)

and for ,

 ^S(2,1)=1s1+iV2⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝−iV20s100−iV20s1s10−iV200s10−iV2⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠. (16)

Each amplitude can be obtained from any other by following a path connecting them on the circle depicted in Fig. 2. The ansatz (3), (11) will be consistent if the amplitudes are not overdetermined by the -matrices. The amplitude is given by in two ways,

 A5 = ^S(5,4)^S(4,3)^S(3,2)^S(2,1)A1, (17a) A5 = ^S(5,6)^S(6,7)^S(7,8)^S(8,1)A1. (17b)

Because of the relations

 ^S(8,1) =^S(5,4),^S(7,8)=^S(4,3), ^S(6,7) =^S(3,2),^S(5,6)=^S(2,1), (18)

we obtain the analogue of the Yang-Baxter equation (YBE),

 ^S(8,1)^S(7,8)^S(6,7)^S(5,6) =^S(5,6)^S(6,7)^S(7,8)^S(8,1). (19)

This generalized consistency relation is quadrilinear instead of trilinear in the -matrices. The reason is the existence of two different -matrices describing the scattering between particle 1 and 2: Although both and correspond to the exchange of the two particles, they are not identical; if we define and together with , , we obtain the generalized YBE in the following, more familiar form,

 ^S10^S′12^S20^S12=^S12^S20^S′12^S10. (20)

This relation is formally equivalent to Sklyanin’s condition for integrable boundary terms sklyanin (); the difference lies in the physical interpretation: in the systems considered by Sklyanin, the particles are only reflected but not transmitted by the boundary potential (the amplitudes are local). Due to the parity invariance of (1), the -matrices are block-diagonal in the eigenbasis of . One finds for the even sector,

 ^S+10= ^S+20= 1s2−iV2⎛⎝iV2s2s2iV2⎞⎠,^S+12=(α−1001), (21)

and for the odd sector,

 ^S−10= −1s1−iV2⎛⎝−iV2s1s1−iV2⎞⎠,^S′−12=(100β), ^S−20= 1s2−iV2⎛⎝iV2s2s2iV2⎞⎠,^S−12=(100α−1). (22)

One checks now that (20) is violated for in the even sector but satisfied identically in in the odd sector. The ansatz (3) is therefore consistent for parity-odd states but fails in the even sector (see Sec. IV).

#### ii.1.2 Bethe ansatz equations

We confine the analysis from now on to the odd sector. The momenta are quantized by the boundary conditions (BC). For periodic BC we have . Open BC correspond to . We shall now demonstrate how periodic BC can be implemented within our non-local framework, leaving open BC for future study (they are treated numerically in Sec. III).

Because of the non-local nature of the amplitudes , a given particle (say particle 1) moves in both directions simultaneously if the regions are connected via the -matrices . It is therefore convenient to treat each component separately, before projection onto the odd subspace. For example, the amplitude describes particle 1 with negative momentum to the right of the impurity, . maps it to to the left of site 0. corresponds to crossing the border between and , leading to , of course without changing the value of amplitude . Now particle 1 is close to the left border of the system, periodicity for values means proportionality of to with . We have

 ψ∝A33ei(−k1)(−M′−x)=A33′ei(−k1)(M′−x+1), (23)

or . Now the configuration of can be recognized as that of , which is mapped by back onto . Similarly, describes a right-moving particle with positive momentum , which is transferred by to , and in turn by to . This time the exchange between the two particles leads to a non-trivial phase factor, because the border between and is crossed. The regions and with are related by periodicity as

 ψ∝A13eik1(M′+x)=A13′eik1(−M′−1+x), (24)

meaning . corresponds to , mapped by onto the original amplitude . The phase factor picked up during transfer from to (see the red arrows in Fig. 2) is independent from the sign of the momenta. Collecting all amplitudes and projecting onto the odd subspace, we conclude that the odd part of must be an eigenvector of

 ^M1=1s1+iV2⎛⎝s1β−1iV2β−1iV2α−1s1α−1⎞⎠ (25)

with eigenvalue . In an analogous way, by transporting particle 2 around the ring, one obtains that the odd part of is an eigenvector of

 ^M2=1s2+iV2⎛⎝s2β−1−iV2αβ−1−iV2s2α⎞⎠ (26)

with eigenvalue . Indeed, , so both matrices are simultaneously diagonalizable. Moreover, they are unitary for real . Let us denote their eigenbasis as . The corresponding eigenvalues of and read

 λ(1)±= s1(s21−s22+U24)±i[(U2−V2)s21s22+V24(s21+s22+U24)2]1/2s1(s21−s22−U24−UV2)+i[V2(s21−s22−U24)+Us21], (27a) λ(2)±= s2(s22−s21+U24)±i[(U2−V2)s21s22+V24(s21+s22+U24)2]1/2s2(s22−s21−U24−UV2)+i[V2(s22−s21−U24)+Us22]. (27b)

One notes that the eigenvalue of obtains from the corresponding eigenvalue of by exchange of and , as expected. If one writes

 λ(j)±=aj±ibcj+idj, (28)

we have

 a2j+b2=c2j+d2j, (29)

which guarantees unimodularity of for real . The Bethe ansatz equations read

 e−ik1M=λ(1)±(k1,k2),e−ik2M=λ(2)±(k1,k2), (30)

which are coupled equations for and . In general both and yield solutions obtained from and , respectively.

An interesting detail can be read off from (27) immediately: if and () the eigenvalues () become 1 for all . This means that all states in the subspaces () have momenta quantized in multiples of as in the non-interacting system. The effects of the impurity and the interaction (being equal but with opposite sign) compensate each other in one of the invariant subspaces and the particles behave in this subspace as if they were free bosons on a ring with sites.

Although (30) has formally the same structure as the Bethe ansatz equations in the two-particle sector of known integrable systems, the right-hand sides are much more complicated, leading to more types of possible solutions. It turns out that there are not just real solutions and strings andrei () but two additional types, which we shall discuss in detail in the next section, presenting a numerical analysis of the eigenstates.

## Iii Numerical analysis in three aspects

In this section we study the model via exact diagonalization on finite lattices, for both open and periodic BC. We shall use three methods to analyze the numerical data. These will corroborate the results of Sec. II.1 and provide numerical evidence that the odd sector is integrable regardless of the boundary condition or the parity of the lattice size.

### iii.1 Spectral graph

The first evidence that the odd-parity subspace is integrable while the even one is not comes from the behavior of the energy spectrum of the model as a function of (or as well).

In Fig. 3, we plot the spectrum of the model as a function of in the odd and even parity subspaces and with open or periodic boundary conditions. We see that regardless of the boundary condition, levels in the odd subspace cross each other while levels in the even subspace repel each other. These two different behaviors strongly hint on the integrability of the odd subspace and the nonintegrability of the even subspace stepanov ().

### iii.2 Bethe-form checking

By full exact diagonalization of the Hamiltonian, we can obtain numerical values of all the eigenvalues and eigenstates of the model. A natural question is then whether it is possible to demonstrate that a wave function, whose value at each site is known, can or cannot be decomposed into the form of (3). To answer this question, first we note that, if the wave function has this form, its values on a line which lies entirely in a particular region (say, the line , in region ), is a superposition of at most four exponentials (if the four wave vectors and are all different, as is demonstrated below). Therefore, first we need to check whether the wave function evaluated on this line is a superposition of four exponentials.

#### iii.2.1 Prony’s algorithm

Fortunately, there exists a beautiful algorithm, namely Prony’s algorithm prony (), which can be used to check whether a function is a superposition of several exponentials and if it is, to extract the exponents. The logic of this algorithm is as follows. Let a function defined on the integers be a superposition of four exponentials,

 gn=w1ec1n+w2ec2n+w3ec3n+w4ec4n, (31)

where the exponents ’s are assumed to be all different, and the ’s are all constant coefficients. Here we take four exponentials merely because of the context—the algorithm works for an arbitrary number of exponentials. Because the ’s are all different, the Vandermonde matrix

 K=⎛⎜ ⎜ ⎜ ⎜⎝1ec1e2c1e3c11ec2e2c2e3c21ec3e2c3e3c31ec4e2c4e3c4⎞⎟ ⎟ ⎟ ⎟⎠ (32)

is non-singular and the linear equation

 ⎛⎜ ⎜ ⎜ ⎜⎝1ec1e2c1e3c11ec2e2c2e3c21ec3e2c3e3c31ec4e2c4e3c4⎞⎟ ⎟ ⎟ ⎟⎠⎛⎜ ⎜ ⎜⎝r0r1r2r3⎞⎟ ⎟ ⎟⎠=⎛⎜ ⎜ ⎜ ⎜⎝e4c1e4c2e4c3e4c4⎞⎟ ⎟ ⎟ ⎟⎠ (33)

has a unique solution for the ’s. It is readily seen that this fact implies the linear iterative relation

 gn+4=r3gn+3+r2gn+2+r1gn+1+r0gn, (34)

which holds regardless of the values of the ’s.

Now consider such a linear equation for the ’s (note that the matrix on the left hand side is a Hankel matrix),

 ⎛⎜ ⎜ ⎜⎝gngn+1gn+2gn+3gn+1gn+2gn+3gn+4gn+2gn+3gn+4gn+5gn+3gn+4gn+5gn+6⎞⎟ ⎟ ⎟⎠⎛⎜ ⎜ ⎜⎝z0z1z2z3⎞⎟ ⎟ ⎟⎠=⎛⎜ ⎜ ⎜⎝gn+4gn+5gn+6gn+7⎞⎟ ⎟ ⎟⎠. (35)

If the value of the function is known at each integer , the matrix on the left hand side and the vector on the right hand side can be constructed, and the linear equation is well-posed. Now the point is that, in view of (34), the solution would be equal to if has the form (31). That is, it is a constant vector independent of . This is a necessary condition for the function to have the form (31).

Therefore, by constructing and solving the linear equation in (35) at different , and observing whether the solution varies with , one can determine whether a function has the form (31) or not. If it is, which means the values of the ’s have been determined, one has at least two options to determine the values of the ’s. One can either solve them inversely by using Eq. (33) as we will do below, or by diagonalizing the transfer matrix in the following equation

 ⎛⎜ ⎜ ⎜⎝gn+4gn+3gn+2gn+1⎞⎟ ⎟ ⎟⎠=⎛⎜ ⎜ ⎜⎝r3r2r1r0100001000010⎞⎟ ⎟ ⎟⎠⎛⎜ ⎜ ⎜⎝gn+3gn+2gn+1gn⎞⎟ ⎟ ⎟⎠, (36)

which is a reformulation of the iterative relation (34).

Now let us come back to our problem. Suppose we take the line with and , and consider the wave function on this line as a function of . The range of depends on the boundary condition. In the case of open BC, we have to confine ourself to region and can take values from 0 up to ; while in the case of periodic BC, we can advance into region and can take values from 0 up to . We can then use Prony’s algorithm above to check whether the wave function has the Bethe ansatz form or not.

In our specific problem, the ’s are a permutation of . It is tedious but straightforward to solve the ’s in Eq. (33) as

 r0 = −1, (37a) r1 = eik1+e−ik1+eik2+e−ik2=−E(k1,k2), (37b) r2 = −(eik1+e−ik1)(eik2+e−ik2)−2, (37c) r3 = r1. (37d)

The fact that regardless of the values of is a stringent condition for the wave function to have the Bethe form. Besides this condition, is another very stringent one. Once these two conditions are verified, it is likely that the wave function has indeed the Bethe form. In Fig. 4, we show and as functions of in different subspaces and under different boundary conditions. The lattice has sites and the state in each panel is sampled unbiased as the one whose energy is the closest to 2. We see that, for the odd-parity states, and are indeed independent of the position and take the values expected. On the contrary, for the even-parity states, they both fluctuate significantly, which rules out the possibility that the corresponding wave function has the Bethe form.

By Eqs. (37b) and (37c), we can solve from the values of the values of (), from which in turn we can solve the wave vectors . Finally, we can plug into (3), solve for the amplitudes and compare the result with the one obtained by exact diagonalization. The inner product between them should be unity if they are the same.

We have carried out this algorithm on lattices with sizes between 15 and 31 and with both kinds of boundary condition. We used the Multiple Precision Toolbox for MATLAB matlab () to solve for the eigenstates and eigenvalues to high precisions (say, with 40-60 digits) and then analyzed the eigenvectors one by one by using Prony’s algorithm. In this way, we confirmed that all the odd-parity eigenstates have the form (3) while none of the even-parity eigenstates fit to it.

Here some remarks are necessary. First, we need high precisions for the eigenstates because - as we shall see in the following - some eigenstates decay exponentially in at least one direction. If the precision is not high enough, the true values of the ’s in (35) can be overwhelmed by the noise by orders of magnitude. Second, Prony’s algorithm is very sensitive to any perturbation to the wave function. Specifically, for the weakly diffractive states in Figs. 7 and 8 below, we get ’s as random as in Figs. 4(b) and 4(d). Therefore, the algorithm can only tell whether a state is exactly in the Bethe form or not. It cannot distinguish weakly diffractive states from strongly diffractive states.

#### iii.2.2 Four main categories of states

Previously it was shown that the model on the infinite line has three continuum bands with different nature prl2012 (); pra2013 (). It turns out that in correspondence there are three main categories of -pair in the odd sector, which are shown in the three panels of Fig. 5 separately and explained below (because the wave function is real, the ’s are real and thus must be real simultaneously or complex conjugate to each other):

(i) are both real, and . In this case, are both real and we can always choose them as [see Fig. 5(a)]. The realness of means that the two particles are both delocalized on the lattice and they are not bound to each other. This is the feature of the first-band states prl2012 (); pra2013 ().

(ii) are both real, and . In this case, one wave vector is real while the other is complex (depending on the sign of , it can be purely imaginary or have a real part of ). We observe that [see Fig. 5(b)], the magnitude of the imaginary part can be well approximated by

 νimp=asinh(|V|/2)>0. (38)

Here is the inverse of the localization length of the defect mode induced by the impurity. This kind of indicates that one particle is localized by the impurity potential while the other is not, in accordance with the picture for the second band prl2012 (); pra2013 ().

(iii) are both complex, and . As a result, are complex and can be chosen in the form with and . It is observed that [see Fig. 5(c)], generally the relation

 ν=asi