Thermodynamic limit and boundary energy of the spin-1 Heisenberg chain with non-diagonal boundary fields

The thermodynamic limit and boundary energy of the isotropic spin-1 Heisenberg chain with non-diagonal boundary fields are studied. The finite size scaling properties of the inhomogeneous term in the $ T-Q $ relation at the ground state are calculated by the density matrix renormalization group. Based on our findings, the boundary energy of the system in the thermodynamic limit can be obtained from Bethe ansatz equations of a related model with parallel boundary fields. These results can be generalized to the $SU(2)$ symmetric high spin Heisenberg model directly.


Introduction
The study of quantum integrable models is an interesting subject in the fields of cold atoms, quantum field theory, condensed matter physics and statistic mechanics [1,2,3,4,5]. The spin-1/2 Heisenberg model can effectively quantify the spin-exchanging interaction and plays an important role in the quantum magnetism and many-body theory. By using the Bethe ansatz method, the one-dimensional (1D) spin-1/2 Heisenberg model can be solved exactly [6]. The typical spin-exchanging couplings in the 1D spin-1 system are characterized by the bilinear biquadratic model, where the Hamiltonian reads (1.1) Here S k (S x k , S y k , S z k ) is the spin-1 operator at site k, N is the number of sites, and the periodic boundary condition gives S N +1 = S 1 . If J 2 /J 1 = 1, the system (1.1) has the SU (3) symmetry and is integrable. If J 2 /J 1 = −1, the SU (2) symmetry exists, and the system is known as the Zamalodchikov-Fateev (ZF) model [7]. The Bethe ansatz solution and thermodynamic properties of the ZF model are studied by Takhtajan [8] and Babujian [9,10]. If J 2 = 0, the system is no longer integrable. Starting from the nonlinear sigma model, Haldane conjectures that the excitation of the system has a gap [11,12]. If J 2 /J 1 = 1/3, the Hamiltonian (1.1) degenerates into a projector operator that is in fact the projection onto the sum of the spin-0 and spin-1 subspaces (up to a constant) and the ground state is the famous valence bond solid state [13,14]. If J 1 = 0, by using the Temperley-Lieb algebra, the system can be mapped into the XXZ spin chain and is also integrable [15,16,17].
Besides the periodic boundary condition, the integrable open one is also an interesting subject, which means that the system has magnetic impurity or the boundary magnetic fields [18,19]. In the past few decades, the exact results of high spin models with periodic [7,8,9,10,20,21,22,23,24,25] and parallel boundary fields [26,27,28,29] have been extensively studied. It is emphasized that the integrable boundary reflection matrix can have non-diagonal elements, which means that the boundary fields are unparallel. Then the U (1) symmetry is broken and it is very hard to study the exact solution of the system. It is known that the integrable systems without U (1) symmetry have many applications in the open string theory and the stochastic process of nonequilibrium statistics. Therefore, many interesting works of high spin models with non-diagonal boundary reflections have been done [30,31,32,33,34,35].
Many attentions have been paid for quantum integrable models without U (1) symmetry during past decades [36,37,38,39,40,41,42,43,44,45,46,47,48,49]. Recently, a systematic method, i.e., the off-diagonal Bethe ansatz (ODBA) is proposed to solve the models with or without U (1) symmetry [50]. Eigenvalues and eigenstates of several typical integrable models are obtained, where eigenvalues are given in terms of some homogeneous/inhomogeneous T − Q relation [50,51,52,53]. The next task is to derive the physical quantities in the thermodynamic limit, which is very complicated because the related Bethe ansatz equations (BAEs) are inhomogeneous and the traditional thermodynamic Bethe ansatz can not be employed. In order to overcome this difficulty, an effective method is to study the finite size scaling effects of the inhomogeneous term in the T − Q relation. With the help of this idea, the thermodynamic limit, surface energy and elementary excitations of spin-1/2 XXZ spin chain with arbitrary boundary fields are studied [54]. The boundary energy of the SU (3) symmetric spin-1 chain with generic integrable open boundaries is also obtained [55]. However, the corresponding thermodynamic properties of the SU (2) symmetric spin-1 Heisenberg model are still missing.
In this paper, we study the thermodynamic limit and boundary energy of the spin-1 isotropic Heisenberg spin chain with non-diagonal boundary reflections. The finite size scaling analysis of the contribution of the inhomogeneous term in the T − Q relation (namely, the third term in (2.17) below) to the ground state energy is studied as follows. We first introduce a very function Λ hom (u) which is given in terms of a reduced T − Q relation 2 (see (3.1) and (3.2) below) [51,52,53] and the associated BAEs are homogeneous ones (see (3.3) below). For any finite N , Λ hom (u) is actually not an eigenvalue of the transfer matrix with generic off-diagonal boundary K-matrices. Since that the function is given by a homogeneous T − Q relation, we can apply the conventional thermodynamic Bethe ansatz [2] to investigate its thermodynamic limit. Then, comparing with the result of its thermodynamic limit and that of the density matrix renormalization group (DMRG) numerical [56,57,58] studies, we conclude that Λ hom (u), in the limit N → ∞, really gives the correct boundary energy. Moreover, we find that most Bethe roots of the reduced BAEs at the ground state in the thermodynamic limit form 2-strings, associated with certain boundary strings and the 2 The function Λ hom (u) can be simulated by eigenvalue of the transfer matrix with parallel boundary fields of the strengthes: p → p/ 1 + α 2 − ; q → q/ 1 + α 2 + . rearrangement of the Fermi sea. The different structures of Bethe roots in different regimes of model parameters are given explicitly. Based on them, we obtain the boundary energy induced by the boundary magnetic fields. We also check the analytic results by the numerical extrapolation, and find that the analytical results and the numerical ones coincide with each other very well. The results given in this paper can be generalized to the SU (2) symmetric spin-s Heisenberg model directly.
This paper is organized as follows. Section 2 serves as an introduction to the notations for the spin-1 Heisenberg model with non-diagonal boundary fields. The ODBA exact solution is also briefly reviewed. In Section 3, we focus on the contribution of the inhomogeneous term in the T − Q relation to the ground state energy. In Section 4, by using the patterns of Bethe roots of the reduced BAEs, we study the boundary energy of the model in the thermodynamic limit. We summarize the results and give some discussions in Section 5.

Non-diagonal boundary Spin-1 Heisenberg model
The spin-1 Heisenberg model with non-diagonal boundary fields is related to the 19-vertex where the non-vanishing elements are u is the spectral parameter, and η is the crossing parameter. Here we are dealing with the isotropic model, and η can be scaled out. Throughout this paper, we adopt the standard notations. For any matrix A ∈ End(V), A j is an embedding operator in the tensor space V ⊗ V ⊗ · · · , which acts as A on the j-th space and as identity on the other factor spaces.
acts as an identity on the factor spaces except for the i-th and j-th ones. The R-matrix R 12 (u) satisfies the quantum Yang-Baxter equation (QYBE) [59,60] Besides, the R-matrix (2.1) also enjoys the properties Initial condition : R 12 (0) = 2η 2 P 12 , where P 12 is the permutation operator and P 12 is the projector in the total spin-0 channel. The most general off-diagonal boundary reflection on one side of the chain is quantified by the reflection matrix obtained in [61,62] where the matrix elements are p − , α − and φ − are the boundary parameters which measure the strength and direction of the boundary field. The reflection matrix K − (u) satisfies the reflection equation (RE) The most general off-diagonal boundary reflection at the other side is quantified by the dual reflection matrix where p + , α + and φ + are the boundary parameters characterizing the strength and direction of the corresponding boundary field. The dual reflection matrix K + (u) satisfies the dual RE From the R-matrix (2.1), we construct the single row monodromy matrices T 0 (u) and T 0 (u) as where {θ k , k = 1, · · · , N } are the inhomogeneous parameters, and the subscript 0 means the auxiliary space and 1, · · · , N denote the quantum spaces. The single row monodromy matrices T 0 (u) andT 0 (u) are the 3×3 matrices in the auxiliary space V 0 and their elements act on the quantum space V ⊗N . The transfer matrix of the system reads (2.14) Now, we seek the exact solution of the system (2.14). Let |Ψ be an arbitrary eigenstate of t(u) with the eigenvalue Λ(u), i.e., Using the ODBA method [50] and fusion hierarchy, in the homogeneous limit {θ k = 0}, the eigenvalue Λ(u) can be expressed as the inhomogeneous T − Q relation, The eigenvalue of Hamiltonian (2.14) reads where {u k } should satisfy the BAEs (2.23) and . is naturally reduced to the conventional diagonal one [30] obtained by the algebraic Bethe Ansatz. 3 For the other case with unparallel boundary fields, the parameter c does not vanish.
Thus the corresponding T − Q relation has to include a non-vanishing inhomogeneous term for any finite N .
The present BAEs (2.23) are inhomogeneous, thus it is very hard to investigate the thermodynamic properties of the system by using the traditional thermodynamic Bethe ansatz. In order to overcome this difficulty, we first analyze the contribution of inhomogeneous term in the T − Q relation (2.17).
Define the reduced T − Q relation as It should be emphasized that although the non-diagonal boundary parameters {p ± , α ± } except φ ± are included in the above reduced T − Q relation (3.2), the Λ hom (u) is not the eigenvalue Λ(u) for any finite N but rather that of the transfer matrix with parallel boundary fields of the same strength. In the limit N → ∞ it will give, however, the correct boundary energy (see the following parts of the paper). From the singularity analysis of the reduced T − Q relation (3.2), we obtain the following reduced BAEs where M = 1, · · · , 2N and we have put η = 1, From the Λ hom (u) given by Eq.(3.1), we obtain the reduced energy which is defined as Taking the logarithm of above Eq.(3.5), we obtain 2πI j = W (λ j ; M 1 ) + θ 2p−1 (λ j ) + θ 2p+1 (λ j ) + θ 2q−1 (λ j ) + θ 2q+1 (λ j ), j = 1, · · · , M 1 , (3.6) where [2θ 2 (λ j − λ l ) + 2θ 2 (λ j + λ l ) + θ 4 (λ j − λ l ) + θ 4 (λ j + λ l )] , (3.7) I j is the quantum number, θ n (x) = 2 arctan(2x/n) and M 1 = N/2. The ground state is characterized by the set of quantum numbers Solving the reduced BAEs (3.6) and substituting the values of Bethe roots into Eq.(3.4), we obtain the reduced ground state energy as Now, we are ready to characterize the contribution of inhomogeneous term in the T − Q relation (2.17) at the ground state by the quantity where E hom is the reduced ground state energy given by (3.9) and E g is the actual ground state energy (2.24) of the Hamiltonian (2.14). The ground state energy E g can be obtained by two methods. One is solving the inhomogeneous BAEs (2.23) directly and the other is DMRG [56,57,58]. We have checked that the ground state energy E obtained by these two methods are the same.
In Figure 2(a), we give the values of E inh versus the system size N in the regime I. The red circles are the data calculated from Eq.(3.10) and the blue solid line is the fitted curve.
From the fitted curve, we find that E inh and N satisfy the power law relation E inh = γN β .
Due to the fact that β < 0, the value of E inh tends to zero when the system size N tends to infinity. Therefore, in the thermodynamic limit, the inhomogeneous term in the T − Q relation (2.17) can be neglected at the ground state and E hom = E g . The inset shows the distribution of Bethe roots with N = 10.
where W (λ j ; M 2 ) is given by Eq.
The corresponding reduced ground state energy reads where G(λ j ; M 2 ) is given by Eq.(3.9) with the replacing of M 1 by M 2 .
The procedure in the regime III is similar and reduced ground state energy is (3.14) In the regime IV, substituting the string solutions including four boundary strings into Eq.(3.3) and taking the logarithm, we have where M 3 = N/2 − 2 and the quantum numbers are The reduced ground state energy is In the regime V, the logarithm form of the BAEs are where M 4 = N/2 − 1 and the quantum numbers are {I j } = {1, 2, · · · , M 4 }. We shall note that the quantum number corresponding to the real Bethe root λ 0 is 0. The reduced ground state energy reads Similarly, the reduced ground state energy in the regime VI is Substituting the reduced ground state energies in different regimes into Eq. Due to higher order correction terms, the effective exponents β determined in the paper differ from −1.

Boundary energy
In this section, we study the physical effects induced by the boundary magnetic fields and compute the boundary energy in the thermodynamic limit [18,35,67,68,69]. As mentioned above, we can calculate the boundary energy based on the string hypothesis of the reduced BAEs (3.3), then the numerical analysis allows us to obtain the boundary energy induced by the boundary fields.
The values of Bethe roots at the ground state are determined by the quantum numbers {I j }. Thus we define the counting function as Z(λ j ) = I j 2N . In the thermodynamic limit, the Bethe roots can take the continuous values and we have Z(λ j ) → Z(u). Taking the derivative of Z(u) with respect to u, we obtain where ρ(u) is the density of Bethe roots and ρ h (u) means the density of holes in the real axis. Again, the distribution of Bethe roots in different regimes are different. We should consider them separately. In regime I, from the BAEs (3.6) with the constraint N → ∞ and using Eq.(4.1), we obtain the density of states as where a n (u) = 1 2π We should note that the presence of delta-function in Eq.(4.2) is due to that λ j = 0 is the solution of BAEs (3.6), which should be excluded because it makes the wavefunction vanish identically [70]. Note that two holes λ h 1 and λ h 2 are introduced to ensure the magnetization satisfying Thus the holes are located at the infinities in the real axis.
With the help of Fourier transformatioñ from Eq.(4.2), we obtainρ (ω) =ρ g (ω) +ρ 0 (ω) +ρ 1 (ω) +ρ 2 (ω), (4.6) wherẽ a n (ω) = e − n|ω| 2 ,ρ g (ω) =ã 1 (ω) +ã 3 (ω) 1 + 2ã 2 (ω) +ã 4 (ω) ,ρ 0 (ω) = 1 2Nã (4.7) Then the ground state energy (3.9) can be expressed as where e g is the ground state energy density which is the same as that for the periodic boundary condition [9], and e s is boundary energy e s = 2π − 4 + E 0 + e 1 + e 2 , (4.10) (4.11) In order to show that there exist the stable boundary bound states, we denote the deviation between ρ p (u) and ρ(u) as δρ p (u) = ρ p (u) − ρ(u). From Eqs.(4.2) and (4.13), we obtain (4.14) Taking the Fourier transformation of Eq.(4.14), we have The energy deviation δe p induced by the density deviation δρ p (ω) can be expressed as Because of δe p < 0, the boundary strings are stable. Then we conclude that in this regime, the ground state energy of the system is E g = N e g + e s + δe p . The total spin along the z-direction is S z = − ∞ −∞ δρ p (u)du = 3/4. Next, we consider the regime III where boundary strings are qi and (q − 1)i. Similarly, the energy deviation δe q between this case and that without boundary strings is Due to the fact δe q < 0, we know that the ground state energy is E g = N e g + e s + δe q and the total spin along the z-direction is S z = 3/4.
In the regime IV, we combine the results (4.16) and (4.17), and conclude that the ground state energy with boundary strings pi, (p − 1)i, qi and (q − 1)i equals to E g = N e g + e s + δe p + δe q .
In the regime VI, there are N − 4 2-string, one real Bethe root λ 0 and three boundary strings qi, pi and (p − 1)i. Combining the results (4.16) and (4.21), we obtain the ground state energy as E g = N e g + e s + δe p + δe λq .
After tedious calculation, we find that the boundary energy e b for all the regimes in Figure 1 can be expressed as  Figure 3 as the red points. Specifically, for each red point that is for the given boundary parameters p and q, we first calculate the ground state energy E g (N ) of the model (2.14) with the system size N = 10(n − 1) + 4 and n = 1, 2, · · · , 20 by using the DMRG method. Then we consider the physical quantity where e g = −1 is the ground state energy density of the system with periodic boundary conditions. Obviously, in the thermodynamic limit, the value of e b (N → ∞) gives the boundary energy. In Figure 4, we show how to extrapolate the boundary energy, where the red points are the numerical values of e b (N ), the blue solid line is the fitting curve, and the red solid line is the extrapolated boundary energy. From the fitting curve, we find that the e b (N ) and N satisfy the power law relation, i.e., e b (N ) = aN β + c. Due to the fact that β < 0, the values of e b (N ) tend to the asymptotic value c when the system size N tends to infinity. Therefore, in the thermodynamic limit, the asymptotic value c determines the boundary energy. Repeating this process, we obtain the boundary energies with other values of boundary parameters. As shown in Figure 3, the analytical and numerical results agree with each other very well.

Conclusions
In this paper, we have studied the thermodynamic limit and boundary energy of the isotropic spin-1 Heisenberg chain with generic integrable non-diagonal boundary reflections. It is shown that the contribution of the inhomogeneous term in the associated T − Q relation (2.17) (due to the unparallel boundary fields) at the ground state can be neglected when the system size N tend to infinity. Then we calculate the analytical expression of boundary energy (4.22) in the thermodynamic limit based on the string hypothesis of the reduced