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

Thank you very much for your many helpful points raised about our paper and we really appreciate your support. We believe that these comments and suggestions have significantly improved our manuscript. To these comments we respond as follows. The page numbers and equation numbers refer to revised version, unless specify.

1. Reviewer: The main results derived are for boundary conditions different from those declared from the author. $\quad$The manuscript presents some interesting and apparently new results on the surface energy of the open integrable spin-1 Heisenberg chain while the authors also claim to have derived some interesting numerical results for large chains. The presentation of the material is however strongly misleading if not even incorrect, leading to state analytic results not for the proper boundary conditions. $\quad$They call this the "reduced Bethe equation", but in fact they are just the Bethe equations of the same open spin-1 XXX chain for different boundary conditions. This changes the boundary conditions in "diagonal ones". This changes the u (spectral parameter) asymptotic behavior of the transfer matrix allowing to the eigenvalues of the (1/2,1) transfer matrix to satisfy the homogeneous TQ-equation. Then, the analysis developed in section 3 and 4 on the Bethe roots and the boundary energy in the thermodynamic limit, in fact, are for these "diagonal boundary conditions" and not for the un non-diagonal boundary magnetic fields case. Indeed, in section 4, for all the considered regimes I-IV, they always refer to the Bethe equations associated to the homogeneous TQ-equations. This is very different from what the authors write already in the title and in their abstract, see also for example the statements at the beginning of section 4: “In this section, we study the physical effects induced by the unparallel boundary magnetic fields and compute the boundary energy”. Authors: The motivation of this paper is that we want to calculate the physical quantities such as the boundary energy of the spin-1 Heisenberg chain induced by the unparallel boundary magnetic fields. The exact energy spectrum of the system is characterized by the inhomogeneous $T-Q$ relations (17)-(23) and associated Bethe ansatz equations (24) for any finite $N$. Due to the existence of inhomogeneous terms in Eq.(24) induced by the unparallel boundary fields, the traditional thermodynamic Bethe asnatz method does not work. By using the density matrix renormalization group simulation, we find that the contribution of inhomogeneous terms in the $T-Q$ relations can be neglected at the ground state only in the thermodynamic limit. Namely, one of most important results of this paper is that we show here the inhomogeneous $T-Q$ relation (17)-(23) and Bethe ansatz equations (24) for the spin-1 Heisenberg chain with the unparallel boundary fields, only in the thermodynamic limit, can be simulated by the homogeneous ones (27)-(29) (which are equivalent to those with diagonal boundary fields). The physical picture is that the two boundary fields can not "see" each other only in the thermodynamic limit (i.e., the contributions to the boundary energy of the two boundary fields independently). [We shall note that the eigenstates with diagonal and those with non-diagonal boundary fields are quite different.] From this equivalent diagonal boundary fields (only in thermodynamic limit), we then obtain the analytical expression for the boundary energy. We note that the boundary energy (68) includes the contribution of off-diagonal boundary parameters $\alpha_{\pm}$. All the results are reasonable and consistent with each other.

2. Reviewer: I found moreover, that the manuscript is rather scarce in the citation of the relevant literature and too much centered on the previous production of the authors. Authors: Thank you for your kind advice. We deleted some self-citations and added three new references: [33] J. Cao, W.-L. Yang, K. Shi and Y. Wang, Off-diagonal Bethe ansatz solution of the XXX spin-chain with arbitrary boundary conditions, Nucl. Phys. B 875, 152 (2013). [34] R.I. Nepomechie, An inhomogeneous $T-Q$ equation for the open XXX chain with general boundary terms: completeness and arbitrary spin, J. Phys. A: Math. Theor. 46, 442002 (2013). [35] R.I. Nepomechie and C. Wang, Boundary energy of the open XXX chain with a non-diagonal boundary term, J. Phys. A: Math. Theor. 47, 032001 (2014).

3. Reviewer: However, one should stress that the thermodynamic convergence of the ground state boundary energy density of the non-diagonal and diagonal cases does not imply that in general one can just consider the "reduced" Bethe equations to develop the thermodynamic analysis of the relevant physical quantities in the non-diagonal case. This has been stressed already by the second referee and I agree with the first one on the relevance of having an analytic derivation of the $\beta$ as functions of the boundary parameters. For example, this type of finite size analysis can be important for relevant physical quantities such as the correlation functions where a modification of order 1/N in the ground state Bethe root density can produce a modification of the correlation functions; e.g. this is clearly the case if one thinks to the exact results for the XXZ spin 1/2 chain. Authors: Very good suggestion! The finite size correction of quantum integrable system is an interesting issue. Based on the logarithmic and integral forms of Bethe ansatz equations, some physical quantities such as the long-range correlation function, elementary excitation, dressed energy, ground state energy correction, susceptibility and Drude weight can be calculated. In the references [A. Klumper, M. T. Batchelor and P. A. Pearce, Central charges of the 6- and 19-vertex models with twisted boundary conditions, J. Phys. A 24, 3111 (1991); J Suzuki, Spinons in magnetic chains of arbitrary spins at finite temperatures, J. Phys. A: Math. Gen. 32, 2341 (1999)], the results of some models with periodic boundary condition have been obtained. We think that this method is also valid for the open boundary conditions, if the Bethe ansatz equations are homogeneous. Because the starting points of this method are the energy spectrum and Bethe ansatz equations, and the processes of treating are the same. Obviously, the finial physical properties could be different. $\quad$Focus on the present model (15). Starting from the reduced Bethe ansatz equation (29), we can analytically compute the finite size corrections and obtain the approximate solutions. The more accurate solutions should be obtained from the Bethe ansatz equations (24). However, due to the existence of inhomogeneous terms in Eq.(24), it is hard to analytically calculate the finite size correction. We think that this difficulty could be overcome by using the recently proposed $t-W$ relation [Phys. Rev. B 102, 085115 (2020)] or the fusion hierarchy at the special points of $u={\theta_k,k=1,\cdots,N}$ [Phys. Rev. B 103, L220401 (2021)]. The main idea is that the eigenvalue of transfer matrix $t(u)$ can also be characterized by its zero points, instead of the Bethe roots. The zero points satisfy the homogeneous Bethe ansatz equations. Based on them, we can analytically calculate the finite size correction.

4. Reviewer: The numerical analysis developed is mainly stated without presenting in the manuscript enough data. $\quad$It would be nice and useful to have more detailed data of their DMRG analysis and a figure, for each chosen values of $p$ and $q$ for which they have done their computations. In particular, they should plot the DMRG results for the chain of N=4, 14, …,194 and from which the asymptotic values are extrapolated. Authors: Very good suggestion! According to your suggestion, we have added the Figure 4 to show the DMRG analysis. The corresponding description in the manuscript is also revised. $\quad$In order to check the correction of analytical expression of boundary energy (68), we perform the numerical simulation with DMRG algorithm. The ground state energy $E_g(N)$ of the model (15) with finite system size $N$ can be calculated by using the DMRG method. Then we consider the physical quantity

   $$ e_b(N) =E_g(N)-Ne_g,\tag{1} $$
   where $e_g=-1$ is the ground state energy density of the system with periodic boundary condition [Please see Eq.(55) in the manuscript]. Obviously, in the thermodynamic limit, the value of $e_b(N\rightarrow\infty)$ gives the boundary energy. $\quad$In the following Figs.1-4 (please see attachment), we show the numerical results. The quantity $e_b(N)$ is calculated with the system size $N = 10(n-1)+4$ and $n=1,2,\cdots, 20$. The boundary parameters $p$ and $q$ take the values as those given by Fig.3 in the manuscript. The red points are the numerical values of $e_b(N)$, the blue solid lines are the fitting curves, and the red solid lines are the extrapolated boundary energies. From the data, we find that the $e_b(N)$ and $N$ satisfy the power law relation, i.e., $e_b(N)=aN^b+c$. Due to the fact that $b<0$, in the thermodynamic limit, the asymptotic value $c$ determines the boundary energy. Comparing $c$ with the analytical result (68), we find that all the results agree with each other very well. $~$

5. Reviewer: In "The model Hamiltonian is generalized from the transfer matrix $t(u)$ as" it is "generated" and not "generalized". Authors: Thank you for your careful reading. We have replaced the word "generalized" with "generated". $~$
6. Reviewer: The definition of $H$ in terms of the transfer matrix in the first line of (15) is apparently singular in $u=0$. Moreover, it produces a term which is the logarithmic derivative of the transfer matrix which is not the standard definition for the open chain. Authors: Thank you very much for pointing out our misprints. We have corrected the misprint of the definition of the Hamiltonian (15). $~$
7. Reviewer: There is some strange asymmetry between the expressions of the boundary magnetic fields in the site 1 and site $N$ with respect to the boundary parameters minus and plus. A part an overall $1/\eta$ in the site $N$ there is also a different behavior with respect to the boundary parameters $\alpha_{-}$ and $\alpha_{+}$. One expects that when these two values are zero one gets back the diagonal case, i.e., the two boundary fields parallel and oriented along the $z$-direction. However, for $\alpha_{-}=0$ the boundary magnetic field in the site 1 become oriented along the $z$-direction while this is not the case for $\alpha_{+}=0$ for the boundary magnetic field in the site $N$. Even more strange, with the current parametrization, the boundary magnetic field in the site $N$ seems to be orientable in the $z$-direction only for specific values of both the boundary parameters, e.g. ($\alpha{+}=0,p_{+}^2=3\eta^2/4$). If there are no mistakes the authors should comment why from a symmetric writing of the transfer matrix with respect to the boundary matrices (and so the boundary parameters) one gets such asymmetric Hamiltonian. Authors: When the non-diagonal boundary parameters $\alpha_\pm=0$, the Hamiltonian reads
   $$ \begin{eqnarray} H&=&\frac{1}{\eta}\sum_{k=1}^{N-1} \left[\vec S_{k} \cdot \vec S_{k+1} -(\vec S_{k} \cdot \vec S_{k+1})^2 \right] + \frac{1}{p_-^2-\frac{1}{4}\eta^2}\left[ 2p_-S_{1}^z -\eta (S_{1}^z )^2 \right] \ && +\frac{1}{(p_+^2-\frac{1}{4}\eta^2)\eta} \left[ -2p_+ \eta S_{N}^z -\frac{2}{3}p_+^2 [(S_{N}^x )^2+(S_{N}^y )^2+(S_{N}^z )^2] +\frac{\eta^2}{2}[ (S_{N}^x )^2+(S_{N}^y )^2-(S_{N}^z )^2 ]\right] \ && +\frac{\eta}{3p_+^2-\frac{3}{4}\eta^2} +\frac{\eta}{p_-^2-\frac{1}{4}\eta^2}+\frac1{\eta}3N+\frac{4}{\eta} \tag{2} \ &=&\frac{1}{\eta}\sum_{k=1}^{N-1} \left[\vec S_{k} \cdot \vec S_{k+1} -(\vec S_{k} \cdot \vec S_{k+1})^2 \right] + \frac{1}{p_-^2-\frac{1}{4}\eta^2}\left[ 2p_-S_{1}^z -\eta (S_{1}^z )^2 \right] \ && +\frac{1}{p_+^2-\frac{1}{4}\eta^2} \left[ -2p_+ S_{N}^z -\eta (S_{N}^z )^2\right] \ && +\frac{\eta}{p_-^2-\frac{1}{4}\eta^2} + \frac{4\eta^2-4p_+^2}{(3p_+^2-\frac{3}{4}\eta^2)\eta}+\frac1{\eta}3N+\frac{4}{\eta},\tag{3} \end{eqnarray} $$
   where we have used the identity $(S^x)^2+(S^y)^2+(S^z)^2=2$. We see that in the fused Hamiltonian $(3)$, the boundary fields are indeed parallel and both of them are along the $z$-direction. Putting $\eta=1$ and changing the signs of boundary parameters at one side, the Hamiltonian $(3)$ can also be written into the symmetric form. The asymmetry of the Hamiltonian $(2)$ is apparent.

Attachment:

Figs.1-4.pdf

Thank you very much for your many helpful points raised about our paper and we really appreciate your support. We believe that these comments and suggestions have significantly improved our manuscript. To these comments we respond as follows. The page numbers and equation numbers refer to revised version, unless specify.

1. Reviewer: In Section 3 this equation is approximated to give a "reduced TQ-relation" which, in fact, are nothing but the (homogeneous) eigenvalue equation obtained by means of the algebraic Bethe ansatz for parallel boundary fields (which can be seen by comparing (27)-(29) with, e.g., the rational limit of the corresponding equations in Ref. [26]). The boundary parameters p and q appearing in the reduced equations are the combinations of the original non-diagonal parameters obtained by rotation of the boundary matrices (see Ref. [29]). The identification of root configurations to the reduced Bethe equations and the calculations of the boundary contribution (68) to the ground state energy as function of the diagonal boundary fields $p$ and $q$, however, is new (to my knowledge). Putting this into the context of the more general model and its solution by means of the ODBA, however, is misleading. Authors: Thank you very much! The motivation of this paper is that we want to calculate the physical quantities such as the boundary energy of the spin-1 Heisenberg chain induced by the unparallel boundary magnetic fields. The exact energy spectrum of the system is characterized by the inhomogeneous $T-Q$ relations (17)-(23) and associated Bethe ansatz equations (24) for any finite $N$. Due to the existence of inhomogeneous terms in Eq.(24) induced by the unparallel boundary fields, the traditional thermodynamic Bethe asnatz method does not work. By using the density matrix renormalization group simulation, we find that the contribution of inhomogeneous terms in the $T-Q$ relations can be neglected at the ground state only in the thermodynamic limit. Namely, one of most important results of this paper is that we show here the inhomogeneous $T-Q$ relation (17)-(23) and Bethe ansatz equations (24) for the spin-1 Heisenberg chain with the unparallel boundary fields, only in the thermodynamic limit, can be simulated by the homogeneous ones (27)-(29) (which are equivalent to those with diagonal boundary fields). The physical picture is that the two boundary fields can not "see" each other only in the thermodynamic limit (i.e., the contributions to the boundary energy of the two boundary fields independently). [We shall note that the eigenstates with diagonal and those with non-diagonal boundary fields are quite different.] From this equivalent diagonal boundary fields (only in thermodynamic limit), we then obtain the analytical expression for the boundary energy. We note that the boundary energy (68) includes the contribution of off-diagonal boundary parameters $\alpha_{\pm}$. All the results are reasonable and consistent with each other.

2. Reviewer: The ODBA solution from Ref.[40] leading to the inhomogeneous TQ-equation (17,18) for the transfer matrix eigenvalues is presented in Section 2. (NB: the general off-diagonal boundary matrices should be contain 3 parameters for strength and direction). Authors: Thank you very much for pointing the fact that the general off-diagonal boundary matrices should be contain 3 parameters for strength and direction. We use the general off-diagonal boundary matrix $K^-(u)$ in (7)-(8) with the boundary parameters ${p_-,\,\alpha_-,\phi_-}$ characterizing the strength and direction of the boundary field (resp. $K^+(u)$ in (10) with the boundary parameters ${p_+,\,\alpha_+,\phi_+}$ ). The eigenvalues $\Lambda(u)$ of the transfer matrix is given by the inhomogeneous $T-Q$ relations (17)-(23) for any finite $N$, which implies that the energy of the ground state may depend on the boundary parameters through the values: $p_{\pm},\,\alpha_{\pm}$ and $\phi_--\phi_+$. However, in the thermodynamic limit, the boundary energy only depends on the boundary parameters through the values: $p=\frac{p_{-}}{\sqrt{1+\alpha^2_{-}}}-\frac{1}{2}$ and $q=-\frac{p_{+}}{\sqrt{1+\alpha^2_{+}}}-\frac{1}{2}$ (see (63), (67) and (68)).

Thank you very much for your many helpful points raised about our paper and we really appreciate your support. We believe that these comments and suggestions have significantly improved our manuscript. To these comments we respond as follows. The page numbers and equation numbers refer to revised version, unless specify.

$~\textbf{1. Reviewer:}$ After (1) in line 28 the authors write "the $SU(2)$ symmetry survives". I suggest to drop "survives" as model (1) is everywhere $SU(2)$ invariant, just at the points $J_2/J_1=1$ and at $J_1=0$ the model has higher $SU(3)$ invariance.
$\quad~\textbf{Authors:}$ You are right. In the Hamiltonian (1), the spin-exchanging interaction in the bulk is $J_1\vec{S}_k\cdot \vec{S}_{k+1}+J_2(\vec{S}_k\cdot \vec{S}_{k+1})^2$, which has $SU(2)$ symmetries with arbitrary coupling constants. If $J_2/J_1=1$ or $J_1=0$, the model has the $SU(3)$ invariance. So we change the word "survives" by "exists".

$~\textbf{2. Reviewer:}$ In line 33 "$J_2/J_1=1/3$, the Hamiltonian (1) degenerates into a projector operator" should be made more explicit: the "projector operator" is in fact the projection onto the sum of the spin-0 and spin-1 subspaces.
$\quad~\textbf{Authors:}$ According to your suggestion, we adjust the explanation of projector operator. We modify the sentence in line 33 "the Hamiltonian (1) degenerates into a projector operator" into "the Hamiltonian (1) degenerates into a projector operator that is in fact the projection onto the sum of the spin-0 and spin-1 subspaces".

$~\textbf{3. Reviewer:}$ In line 82 the authors write "$\eta$ is the crossing parameter", but here we are dealing with the isotropic model and there is non-trivial crossing parameter. This parameter can be scaled out as is done in line 123. This should be indicated already in line 82.
$\quad~\textbf{Authors:}$ Thank you very much for your suggestion. In line 82, we have indicated that $\eta$ can be scaled out.

$~\textbf{4. Reviewer:}$ In eq(6) and line 85 the operator $\textbf{P}_{12}^{(0)}$ is called antisymmetry but also "projector in the total spin-0 channel". What is really meant? The anti-symmetric states are the spin-1 states!
$\quad~\textbf{Authors:}$ Thank you for your kind advice. At the point of $u=-\eta$, the $R$-matrix (2) is a singular matrix which only has one non-zero eigenvalue. Therefore, when the spectral parameter $u$ equals to $-\eta$, the $R$-matrix (2) degenerates into a one-dimensional projector operator $\textbf{P}_{12}^{(0)}=|\psi_0\rangle\langle\psi_0|$ and the corresponding basis vector is $|\psi_0\rangle=\frac{1}{\sqrt 3}(|1,-1\rangle-|0,0\rangle+|-1,1\rangle)$, where $|1\rangle$, $|0\rangle$ and $|-1\rangle$ are the eigenstates of spin-1 operator $S^z$ with the eigenvalues 1, 0 and -1, respectively. The operator $\textbf{P}_{12}^{(0)}$ projects the tensor space $\mathbb{V} \otimes \mathbb{V}$ into spin-0 subspace. We modified "Antisymmetry" to "Fusion condition" in Eq.(6).

$~\textbf{5. Reviewer:}$ In line 107 "hierarchy fusion" should be called "fusion hierarchy".
$\quad~\textbf{Authors:}$ Thank you for your careful reading. We have modified "hierarchy fusion" as "fusion hierarchy".

$~\textbf{6. Reviewer:}$ On page 7 and later the terminology "boundary string" is somewhat confusing. It always refers to a single Bethe root.
$\quad~\textbf{Authors:}$ The boundary strings mean that the values of Bethe roots are determined by the boundary parameters. They could be the single value or the bound pairs. For example, the boundary strings in regime II are {$pi, (p-1)i, -pi, -(p-1)i$}, which are clearly related with the boundary fields and distribute symmetrically around the zero point. Due to the symmetry, the boundary strings {$pi, (p-1)i$} and the {$-pi, -(p-1)i$} have the same contribution to the boundary energy.

$~\textbf{7. Reviewer:}$ The word "subfigure" should be replaced by "inset".
$\quad~\textbf{Authors:}$ We have replaced the word "subfigure" with "inset".

$~\textbf{8. Reviewer:}$ Line 213. "Meanwhile" should be replaced by a different word simply by saying "Note that two holes $\lambda_1^h$ and $\lambda_2^h$ are introduced $\cdots$"
$\quad~\textbf{Authors:}$ Thank you for your kind advice. We have replaced "Meanwhile, two holes $ \lambda_1^h $ and $ \lambda_2^h $ should be introduced" with "Note that two holes $\lambda_1^h$ and $\lambda_2^h$ are introduced".

$~\textbf{9. Reviewer:}$ spin long $z$-direction $\rightarrow$ spin along the $z$-direction.
$\quad~\textbf{Authors:}$ Thank you for your careful reading. We have replaced all the "spin long $z$-direction" with "spin along the $z$-direction".

$~\textbf{10. Reviewer:}$ Typo in caption to Fig.3: (68) pand the red points $\rightarrow$ (68) and the red points.
$\quad~\textbf{Authors:}$ Thank you for your careful reading. We have corrected the typo.

$~\textbf{11. Reviewer:}$ The authors should comment on the possibility of evaluating the finite size correction term const.$N^{\beta}$ by analytic means. For periodic boundary conditions the model of interest has been investigated long ago by non-linear integral equations.
A. Klumper, M. T. Batchelor, P. A. Pearce: Central charges of the 6- and 19-vertex models with twisted boundary conditions. J. Phys. A 24, 3111-3133 (1991).
J Suzuki: Spinons in magnetic chains of arbitrary spins at finite temperatures, J. Phys. A: Math. Gen. 32 2341 (1999).
Is a similar treatment conceivable for open boundary conditions?
$\quad~\textbf{Authors:}$ Very good suggestion! The finite size correction of quantum integrable system is an interesting issue. Based on the logarithmic and integral forms of Bethe ansatz equations, some physical quantities such as the long-range correlation function, elementary excitation, dressed energy, ground state energy correction, susceptibility and Drude weight can be calculated. In the references mentioned by the referee, the results of some models with periodic boundary condition have been obtained. We think that this method is also valid for the open boundary conditions, if the Bethe ansatz equations are homogeneous. Because the starting points of this method are the energy spectrum and Bethe ansatz equations, and the processes of treating are the same. Obviously, the finial physical properties could be different.
$\quad$Focus on the present model (15). Starting from the reduced Bethe ansatz equation (29), we can analytically compute the finite size corrections and obtain the approximate solutions. The more accurate solutions should be obtained from the Bethe ansatz equations (24). However, due to the existence of inhomogeneous terms in Eq.(24), it is hard to analytically calculate the finite size correction. We think that this difficulty could be overcome by using the recently proposed $t-W$ relation [Phys. Rev. B 102, 085115 (2020)] or the fusion hierarchy at the special points of $u=\{\theta_k,k=1,\cdots,N\}$ [Phys. Rev. B 103, L220401 (2021)]. The main idea is that the eigenvalue of transfer matrix $t(u)$ can also be characterized by its zero points, instead of the Bethe roots. The zero points satisfy the homogeneous Bethe ansatz equations. Based on them, we can analytically calculate the finite size correction.