Algebraic Bethe ansatz for thesℓ(2)Gaudin model with boundary

Algebraic Bethe ansatz for the sℓ(2) Gaudin model with boundary N. Cirilo António, ∗ N. Manojlović, † E. Ragoucy ‡ and I. Salom § ∗ Centro de Análise Funcional e Aplicações arXiv:1412.1396v2 [nlin.SI] 10 Feb 2015 Instituto Superior Técnico, Universidade de Lisboa Av. Rovisco Pais, 1049-001 Lisboa, Portugal † Grupode Física Matemática da Universidade de Lisboa Av. Prof. Gama Pinto 2, PT-1649-003 Lisboa, Portugal † Departamento de Matemática, F. C. T., Universidade do Algarve Campus de Gambelas, PT-8005-139 Faro, Portugal ‡ Laboratoire d’Annecy-le-Vieux de Physique Théorique LAPTh CNRS et Univerité de Savoie, UMR 5108, B.P. 110 74941 Annecy-le-Vieux Cedex, France § Instituteof Physics, University of Belgrade P.O. Box 57, 11080 Belgrade, Serbia Abstract Following Sklyanin’s proposal in the periodic case, we derive the generating function of the Gaudin Hamiltonians with boundary terms. Our derivation is based on the quasi- classical expansion of the linear combination of the transfer matrix of the XXX Heisenberg spin chain and the central element, the so-called Sklyanin determinant. The corresponding Gaudin Hamiltonians with boundary terms are obtained as the residues of the generating function. By defining the appropriate Bethe vectors which yield strikingly simple off shell action of the generating function, we fully implement the algebraic Bethe ansatz, obtaining the spectrum of the generating function and the corresponding Bethe equations. ∗ E-mail address: nantonio@math.ist.utl.pt † E-mail address: nmanoj@ualg.pt ‡ E-mail address: eric.ragoucy@lapth.cnrs.fr § E-mail address: isalom@ipb.ac.rs  I INTRODUCTION I Introduction A model of interacting spins in a chain was first considered by Gaudin [1, 2]. Gaudin derived these models as a quasi-classical limit of the quantum chains. Sklyanin studied the rational sℓ(2) model in the framework of the quantum inverse scattering method using the sℓ(2) invari- ant classical r-matrix [3]. A generalization of these results to all cases when skew-symmetric r-matrix satisfies the classical Yang-Baxter equation [4] was relatively straightforward [5, 6]. Therefore, considerable attention has been devoted to Gaudin models corresponding to the the classical r-matrices of simple Lie algebras [7–12] and Lie superalgebras [13–17]. Hikami, Kulish and Wadati showed that the quasi-classical expansion of the transfer ma- trix of the periodic chain, calculated at the special values of the spectral parameter, yields the Gaudin Hamiltonians [18, 19]. Hikami showed how the quasi-classical expansion of the transfer matrix, calculated at the special values of the spectral parameter, yields the Gaudin Hamiltonians in the case of non-periodic boundary conditions [20]. Then the ABA was ap- plied to open Gaudin model in the context of the the Vertex-IRF correspondence [21–23]. Also, results were obtained for the open Gaudin models based on Lie superalgebras [24]. An ap- proach to study the open Gaudin models based on the classical reflection equation [25–27] and the non-unitary r-matrices was developed recently, see [28, 29] and the references therein. For a review of the open Gaudin model see [30]. Progress in applying Bethe ansatz to the Heisen- berg spin chain with non-periodic boundary conditions compatible with the integrability of the quantum systems [31–41] had recent impact on the study of the corresponding Gaudin model [41, 42]. The so-called T − Q approach to implementation of Bethe ansatz [35, 36] was used to obtain the eigenvalues of the associated Gaudin Hamiltonians and the corresponding Bethe ansatz equations [42]. In [41] the off shell action of the generating function of the Gaudin Hamiltonians on the Bethe vectors was obtained through the so-called quasi-classical limit. Here we derive the generating function of the Gaudin Hamiltonians with boundary terms following Sklyanin’s approach in the periodic case [3]. Our derivation is based on the quasi- classical expansion of the linear combination of the transfer matrix of the inhomogeneous XXX Heisenberg spin chain and the central element, the so-called Sklyanin determinant. The es- sential step in this derivation is the expansion of the monodromy matrix in powers of the quasi-classical parameter. Moreover, we show how the representation of the relevant Lax ma- trix in terms of local spin operators yields the partial fraction decomposition of the generating function. Consequently, the Gaudin Hamiltonians with the boundary terms are obtained from the residues of the generating function at poles. We derive the relevant linear bracket for the Gaudin Lax operator and certain classical r-matrix, obtained form the sℓ(2) invariant classi- cal r-matrix and the corresponding K-matrix. The local realisation of the Lax matrix together with the linear bracket provide the necessary structure for the implementation of the algebraic Bethe ansatz. In this framework, the Bethe vectors, defined as the symmetric functions of its arguments, have a remarkable property that the off shell action of the generating function on them is strikingly simple. Actually, it is as simple as it can be since it practically coincide with the corresponding formula in the case when the boundary matrix is diagonal [20]. The off -2-  II Sℓ(2) GAUDIN MODEL shell action of the generating function of the Gaudin Hamiltonians with the boundary terms yields the spectrum of the system and the corresponding Bethe equations. As usual, when the Bethe equations are imposed on the parameters of the Bethe vectors, the unwanted terms in the action of the generating function are annihilated. However, more compact form of the Bethe vector ϕ M (µ1 , µ2 , . . . , µ M ), for an arbitrary posi- tive integer M, requires further studies. As it is evident form the formulas for the Bethe vector ϕ4 (µ1 , µ2 , µ3 , µ4 ) given in the Appendix B, the problem lies in the definition the scalar coeffi- (m) cients c M (µ1 , . . . µm ; µm+1 , . . . , µ M ), with m = 1, 2, . . . , M. Some of them are straightforward ( M) to obtain but, in particular, the coefficient c M (µ1 , µ2 . . . , µ M ) still represents a challenge, at least in the present form. This paper is organized as follows. In Section 2 we review the SL(2)-invariant Yang R- matrix and provide fundamental tools for the study of the inhomogeneous XXX Heisenberg spin chain and the corresponding Gaudin model. Moreover, we outline Sklyanin’s derivation of the rational sℓ(2) Gaudin model. The general solutions of the reflection equation and the dual reflection equation are given in Section 3. As one of the main results of the paper, the generating function of the Gaudin Hamiltonians with boundary terms is derived in Section 4, using the quasi-classical expansion of the linear combination of the transfer matrix of the in- homogeneous XXX spin chain and the so-called Sklyanin determinant. The relevant algebraic structure, including the classical reflection equation, is given in Section 5. The implementation of the algebraic Bethe ansatz is presented in Section 6, including the definition of the Bethe vectors and the formulae of the off shell action of the generating function of the Gaudin Hamil- tonians. Our conclusions are presented in the Section 7. Finally, in Appendix A are given some basic definitions for the convenience of the reader. II sℓ(2) Gaudin model The XXX Heisenberg spin chain is related to the SL(2)-invariant Yang R-matrix [43]   λ+η 0 0 0  0 λ η 0  R ( λ ) = λ1 + η P =   0 , (II.1) η λ 0  0 0 0 λ+η where λ is a spectral parameter, η is a quasi-classical parameter, 1 is the identity operator and we use P for the permutation in C2 ⊗ C2 . The Yang R-matrix satisfies the Yang-Baxter equation [43–46] in the space C2 ⊗ C2 ⊗ C2 R12 (λ − µ) R13 (λ) R23 (µ) = R23 (µ) R13 (λ) R12 (λ − µ), (II.2) we use the standard notation of the quantum inverse scattering method to denote spaces on which corresponding R-matrices Rij , ij = 12, 13, 23 act non-trivially and suppress the depen- dence on the quasi-classical parameter η [45, 46]. The Yang R-matrix also satisfies other relevant properties such as -3-  II Sℓ(2) GAUDIN MODEL unitarity R12 (λ) R21 (−λ) = (η 2 − λ2 )1; parity invariance R21 (λ) = R12 (λ); temporal invariance t ( λ ) = R ( λ ); R12 12 crossing symmetry R(λ) = J1 Rt2 (−λ − η )J1 , where t2 denotes the transpose in the second space and the entries of the two-by-two matrix J are J ab = (−1) a−1 δa,3−b . Here we study the inhomogeneous XXX spin chain with N sites, characterised by the local space Vm = C2s+1 and inhomogeneous parameter αm . For simplicity, we start by considering the periodic boundary conditions. The Hilbert space of the system is N H = ⊗ Vm = (C2s+1 )⊗ N . (II.3) m =1 Following [3] we introduce the Lax operator [41]   η  1 λ + ηS3m − ηSm L0m (λ) = 1 + ~σ0 · ~Sm = + . (II.4) λ λ ηSm λ − ηS3m Notice that L (λ) is a two-by-two matrix in the auxiliary space V0 = C2 . It obeys  η 2 c2,m  L0m (λ)L0m (η − λ) = 1 + 10 , (II.5) λ(η − λ) where c2,m is the value of the Casimir operator on the space Vm [41]. When the quantum space is also a spin 21 representation, the Lax operator becomes the R-matrix, L0m (λ) = λ1 R0m (λ − η/2). Due to the commutation relations (A.1), it is straightforward to check that the Lax operator satisfies the RLL-relations R00′ (λ − µ)L0m (λ)L0′ m (µ) = L0′ m (µ)L0m (λ) R00′ (λ − µ). (II.6) The so-called monodromy matrix T (λ) = L0N (λ − α N ) · · · L01 (λ − α1 ) (II.7) is used to describe the system. For simplicity we have omitted the dependence on the quasi- classical parameter η and the inhomogeneous parameters {α j , j = 1, . . . , N }. Notice that T (λ) is a two-by-two matrix acting in the auxiliary space V0 = C2 , whose entries are operators acting in H. From RLL-relations (II.6) it follows that the monodromy matrix satisfies the RTT- relations R00′ (λ − µ) T0 (λ) T0′ (µ) = T0′ (µ) T0 (λ) R00′ (λ − µ). (II.8) The periodic boundary conditions and the RTT-relations (II.8) imply that the transfer matrix t(λ) = tr0 T (λ), (II.9) -4-  II Sℓ(2) GAUDIN MODEL commute at different values of the spectral parameter, [t(µ), t(ν)] = 0, (II.10) here we have omitted the nonessential arguments. The RTT-relations admit a central element − ∆ [ T (λ)] = tr00′ P00 ′ T0 ( λ − η/2) T0′ ( λ + η/2) , (II.11) where − 1 − P00′ 1 P00 ′ = =− R00′ (−η ) . (II.12) 2 2η A straightforward calculation shows that h   i ∆ T (µ) , T (ν) = 0. (II.13) As the first step toward the study of the Gaudin model we consider the expansion of the monodromy matrix (II.7) with respect to the quasi-classical parameter η   N ~ ~σ0 · Sm η 2 N 1 0 ~Sm · ~Sn T (λ) = 1 + η ∑ + ∑ (λ − αm )(λ − αn ) m =1 λ − α m 2 n,m =1 n6=m       η2 N N ı~σ0 · ~Sn × ~Sm ı~σ0 · ~Sm × ~Sn N ∑  + O(η 3 ). 2 m∑ + + ∑ (II.14) =1 n > m ( λ − α m )( λ − α n ) n<m ( λ − α m )( λ − α n ) If the Gaudin Lax matrix is defined by [3] N ~σ0 · ~Sm L0 ( λ ) = ∑ λ − αm (II.15) m =1 and the quasi-classical property of the Yang R-matrix [3] 1 P R(λ) = 1 − ηr(λ), where r (λ) = − (II.16) λ λ is taken into account, then substitution of the expansion (II.14) into the RTT-relations (II.8) yields the so-called Sklyanin linear bracket [3] [ L1 (λ), L2 (µ)] = [r12 (λ − µ), L1 (λ) + L2 (µ)] . (II.17) Using the expansion (II.14) it is evident that N N ~Sm · ~Sn t( λ) = 2 + η 2 ∑ ∑ + O(η 3 ). (II.18) m =1 n 6 = m ( λ − α m )( λ − α n ) -5-  II Sℓ(2) GAUDIN MODEL The same expansion (II.14) leads to ! η2 N ~ σ · ~ S ~ σ · ~ S ∆ [ T (λ)] = 1 + η trL(λ) + tr00′ P00 0 m ′ m − ′ ∑ 2 − 0 2 m =1 ( λ − α m ) ( λ − α m )2       η2 N N 10 Sm · Sn ~ ~ 10 Sm · Sn ′ ~ ~ −   + tr00′ P00 ′ ∑ ∑ + 2 m =1 n 6 = m (λ − αm )(λ − αn ) (λ − αm )(λ − αn )       2 N N ı~σ0 · ~Sn × ~Sm N ı~σ 0 · ~Sm × ~Sn η − ∑  + tr00′ P00 ′ ∑ + ∑ 2 m =1 n > m ( λ − α m )( λ − α n ) n < m ( λ − α m )( λ − α n )       2 N N ı~σ ′ · ~ S n × ~Sm N ı~σ ′ · ~ S m × ~Sn η − ∑ 0 0  + tr00′ P00 ′ ∑ + ∑ 2 m =1 n > m ( λ − α m )( λ − α n ) n < m ( λ − α m )( λ − α n ) − + η 2 tr00′ P00 3 ′ L0 ( λ ) L0′ ( λ ) + O( η ), (II.19) where L(λ) is given in (II.15). The final expression for the expansion of ∆ [ T (λ)] is obtained after taking all the traces ! N N ~ S · ~ S 1 ∆ [ T (λ)] = 1 + η 2 ∑ ∑ m n − tr L2 (λ) + O(η 3 ). (II.20) m =1 n 6 = m ( λ − α m )( λ − α n ) 2 To obtain the generation function of the Gaudin Hamiltonians notice that (II.18) and (II.20) yield η2 t(λ) − ∆ [ T (λ)] = 1 + tr L2 (λ) + O(η 3 ). (II.21) 2 Therefore 1 τ (λ) = tr L2 (λ) (II.22) 2 commute for different values of the spectral parameter, [τ (λ), τ (µ)] = 0. (II.23) Moreover, from (II.15) it is straightforward to obtain the expansion N 2Hm N ~Sm · ~Sm N 2Hm N s m ( s m + 1) τ (λ) = ∑ λ − α m ∑ ( λ − α m )2 = + ∑ λ − α m ∑ ( λ − α m )2 , + (II.24) m =1 m =1 m =1 m =1 and the Gaudin Hamiltonians, in the periodic case, are N ~Sm · ~Sn Hm = ∑ α − αn . (II.25) n6=m m This shows that τ (λ) is the generating function of Gaudin Hamiltonians when the periodic boundary conditions are imposed [3]. -6-  IV Sℓ(2) GAUDIN MODEL WITH BOUNDARY TERMS III Reflection equation A way to introduce non-periodic boundary conditions which are compatible with the integra- bility of the bulk model, was developed in [27]. Boundary conditions on the left and right sites of the system are encoded in the left and right reflection matrices K − and K + . The compatibil- ity condition between the bulk and the boundary of the system takes the form of the so-called reflection equation. It is written in the following form for the left reflection matrix acting on the space C2 at the first site K − (λ) ∈ End(C2 ) R12 (λ − µ)K1− (λ) R21 (λ + µ)K2− (µ) = K2− (µ) R12 (λ + µ)K1− (λ) R21 (λ − µ). (III.1) Due to the properties of the Yang R-matrix the dual reflection equation can be presented in the following form R12 (µ − λ)K1+ (λ) R21 (−λ − µ − 2η )K2+ (µ) = K2+ (µ) R12 (−λ − µ − 2η )K1+ (λ) R21 (µ − λ). (III.2) One can then verify that the mapping K + (λ) = K − (−λ − η ) (III.3) is a bijection between solutions of the reflection equation and the dual reflection equation. Af- ter substitution of (III.3) into the dual reflection equation (III.2) one gets the reflection equation (III.1) with shifted arguments. The general, spectral parameter dependent solutions of the reflection equation (III.1) can be written as follows [47]   ξ − λ ψλ K − (λ) = . (III.4) φλ ξ + λ It is straightforward to check the following useful identities

K − (−λ)K − (λ) = ξ 2 − λ2 (1 + φψ) 1 = det K − (λ) 1, (III.5) − − − K (−λ) = tr K (λ) − K (λ). (III.6) IV sℓ(2) Gaudin model with boundary terms With the aim of describing the inhomogeneous XXX spin chain with non-periodic boundary condition it is instructive to recall some properties of the Lax operator (II.4). The identity (II.5) can be rewritten in the form [41]  η 2 s m ( s m + 1)  L0m (λ − αm )L0m (−λ + αm + η ) = 1 + 10 . (IV.1) (λ − αm )(−λ + αm + η ) It follows from the equation above and the RLL-relations (II.6) that the RTT-relations (II.8) can be recast as follows e0′ (µ) R00′ (λ + µ) T0 (λ) = T0 (λ) R00′ (λ + µ) T T e0′ (µ), (IV.2) e0 (λ) T T e0′ (µ) R00′ (µ − λ) = R00′ (µ − λ) T e0 (λ), e0′ (µ) T (IV.3) -7-  IV Sℓ(2) GAUDIN MODEL WITH BOUNDARY TERMS where e(λ) = L01 (λ + α1 + η ) · · · L0N (λ + α N + η ). T (IV.4) The Sklyanin monodromy matrix T (λ) of the inhomogeneous XXX spin chain with non-periodic e0 (λ) (IV.4) and a reflection matrix K − (λ) boundary consists of the two matrices T (λ) (II.7) and T (III.4), T0 (λ) = T0 (λ)K0− (λ) Te0 (λ). (IV.5) Using the RTT-relations (II.8), (IV.2), (IV.3) and the reflection equation (III.1) it is straightfor- ward to show that the exchange relations of the monodromy matrix T (λ) in V0 ⊗ V0′ are [41] R00′ (λ − µ)T0 (λ) R0′ 0 (λ + µ)T0′ (µ) = T0′ (µ) R00′ (λ + µ)T0 (λ) R0′ 0 (λ − µ), (IV.6) The open chain transfer matrix is given by the trace of T (λ) over the auxiliary space V0 with an extra reflection matrix K + (λ) [27],
t(λ) = tr0 K + (λ)T (λ) . (IV.7) The reflection matrix K + (λ) (III.3) is the corresponding solution of the dual reflection equation (III.2). The commutativity of the transfer matrix for different values of the spectral parameter [t(λ), t(µ)] = 0, (IV.8) is guaranteed by the dual reflection equation (III.2) and the exchange relations (IV.6) of the monodromy matrix T (λ). The exchange relations (IV.6) admit a central element − ∆ [T (λ)] = tr00′ P00 ′ T0 ( λ − η/2) R00′ (2λ )T0′ ( λ + η/2). (IV.9) For the study of the open Gaudin model we impose  
lim K + (λ)K − (λ) = ξ 2 − λ2 (1 + φψ) 1. (IV.10) η →0 In particular, this implies that the parameters of the reflection matrices on the left and on the right end of the chain are the same. In general this not the case in the study of the open spin chain. However, this condition is essential for the Gaudin model. Then we will write K − ( λ ) ≡ K ( λ ), (IV.11) so that   + −1 ψ K (λ) = K (−λ − η ) = K (−λ) − η M with M= . (IV.12) φ 1 Remark that the matrix M obeys M2 = (1 + ψφ)1. -8-  IV Sℓ(2) GAUDIN MODEL WITH BOUNDARY TERMS e(λ) as intro- The expansion of T (λ) is given in (II.14). It is easy to get the expansion for T duced in (IV.4) and then, the one for T (λ). Using the relation (IV.12), we deduce the expansion of t(λ) (VI.24) in powers of η:
t(λ) = 2 ξ 2 − λ2 (1 + φψ) − 2ηλ ( 1 + φψ) − η 2 tr0 ( M0 ( L0 (λ)K0 (λ) − K0 (λ) L0 (−λ))) ! 2 2 2
N ~Sm · ~Sn ~Sm · ~Sn +η ξ − λ (1 + φψ) ∑ + (λ − αm )(λ − αn ) (λ + αm )(λ + αn ) m,n =1 n6=m − η 2 tr0 L0 (λ)K0 (λ) L0 (−λ)K0 (−λ) + O(η 3 ). (IV.13) Our next step is to obtain the expansion of ∆ [T (λ)] (IV.9) in powers of η. We follow the analogous steps as for the periodic case, and after some tedious but straightforward calcula- tions we get
∆ [T (λ)] = λ tr20 K0 (λ) − tr0 K02 (λ) + 2ηλtr0 K0 (λ) tr0 ( L0 (λ)K0 (λ) − K0 (λ) L0 (−λ))     η − 2ηλ tr0 L0 (λ)K02 (λ) − tr0 L0 (−λ)K02 (λ) − tr0 K02 (λ) 2 ! N ~ Sm · Sn~ ~ ~ Sm · Sn + η2 λ ∑ + tr0 K0 (−λ)K0 (λ) m,n =1 (λ − αm )(λ − αn ) (λ + αm )(λ + αn ) n6=m 2 − 2η λ tr0 L0 (λ)K0 (λ) L0 (−λ)K0 (−λ) 
− η 2 tr0 L0 (λ)K0 (λ) − K0 (λ) L0 (−λ) K0 (λ)  2 + η 2 λ tr0 { L0 (λ)K0 (λ) − K0 (λ) L0 (−λ)} n  o − η 2 λ tr0 L0 (λ)K0 (λ) − K0 (λ) L0 (−λ) L0 (λ)K0 (λ) − K0 (λ) L0 (−λ) η2 λ + tr0 M02 + O(η 3 ). (IV.14) 4 Using the relations (III.5) and (III.6) the first term of the expansion above simplifies and the second and third term together turn out to be propositional to the trace of L(λ) (II.15) and therefore vanish, -9-  IV Sℓ(2) GAUDIN MODEL WITH BOUNDARY TERMS

∆ [T (λ)] = 2λ ξ 2 − λ2 (1 + φψ) − η ξ 2 + λ2 (1 + φψ) ! 2 2 2
N ~Sm · ~Sn ~Sm · ~Sn + 2η λ ξ − λ (1 + φψ) ∑ + (λ − αm )(λ − αn ) (λ + αm )(λ + αn ) m,n =1 n6=m − 2η 2 λ tr0 L0 (λ)K0 (λ) L0 (−λ)K0 (−λ) − η 2 tr0 (( L0 (λ)K0 (λ) − K0 (λ) L0 (−λ)) K0 (λ)) n  + η 2 λ tr0 tr0′ { L0′ (λ)K0′ (λ) − K0′ (λ) L0′ (−λ)} − L0 (λ)K0 (λ) + K0 (λ) L0 (−λ) ×  o × L0 (λ)K0 (λ) − K0 (λ) L0 (−λ) η2 λ + (1 + φψ) + O(η 3 ). (IV.15) 2 In order to simplify some formulae we introduce the following notation L0 (λ) = L0 (λ) − K0 (λ) L0 (−λ)K0−1 (λ), (IV.16) Using the formulas (IV.13) and (IV.15) we calculate the expansion in powers of η of the difference

2λt(λ) − ∆ [T (λ)] = 2λ ξ 2 − λ2 (1 + φψ) + η ξ 2 − 3λ2 (1 + φψ) − 2η 2 λ tr0 ( M0 L0 (λ)K0 (λ))
+ η 2 tr0 L0 (λ)K02 (λ) − η 2 λ tr0 ((tr0′ (L0′ (λ)K0′ (λ)) 10 − L0 (λ)K0 (λ)) L0 (λ)K0 (λ)) η2 λ − (1 + φψ) + O(η 3 ). (IV.17) 2 Actually the third and the fourth term in the expression above vanish
tr0 L0 (λ)K02 (λ) − 2λ tr0 ( M0 L0 (λ)K0 (λ)) = tr0 ((L0 (λ)K0 (λ)) (K0 (λ) − 2λM0 ))
= tr0 (L0 (λ)K0 (λ)K0 (−λ)) = ξ 2 − λ2 (1 + φψ) tr0 L0 = 0, (IV.18) due to the fact that the tr0 L0 is equal to zero. Therefore the expansion (IV.17) reads

2λt(λ) − ∆ [T (λ)] = 2λ ξ 2 − λ2 (1 + φψ) + η ξ 2 − 3λ2 (1 + φψ) − η 2 λ tr0 ((tr0′ (L0′ (λ)K0′ (λ)) 10 − L0 (λ)K0 (λ)) L0 (λ)K0 (λ)) η2 λ − (1 + φψ) + O(η 3 ). (IV.19) 2 It is important to notice that using the following identity tr0′ (L0′ (λ)K0′ (λ)) 10 − L0 (λ)K0 (λ) = −K0 (−λ)L0 (λ), (IV.20) the third term in (IV.19) can be simplified
tr0 K0 (−λ)L0 (λ)L0 (λ)K0 (λ) = ξ 2 − λ2 (1 + φψ) tr0 L20 (λ). (IV.21) -10-  IV Sℓ(2) GAUDIN MODEL WITH BOUNDARY TERMS Finally, the expansion (IV.19) reads

2λt(λ) − ∆ [T (λ)] = 2λ ξ 2 − λ2 (1 + φψ) + η ξ 2 − 3λ2 (1 + φψ)
+ η 2 λ ξ 2 − λ2 (1 + φψ) tr0 L20 (λ) η2 λ − (1 + φψ) + O(η 3 ). (IV.22) 2 This shows that τ (λ) = tr0 L20 (λ) (IV.23) commute for different values of the spectral parameter, [τ (λ), τ (µ)] = 0. (IV.24) and therefore can be considered to be the generating function of Gaudin Hamiltonians with boundary terms. The multiplicative factor in (VI.31), which is equal to the determinant of K (λ), will be useful in the partial fraction decomposition of the generating function. With the aim of obtaining the Gaudin Hamiltonians with the boundary terms from the generating function (VI.31), it is instructive to study the representation of L0 (λ) (IV.16) in terms of the local spin operators     −1 ~Sm N ~ ~σ0 · Sm K 0 ( λ )~ σ 0 K 0 ( λ ) · L0 ( λ ) = ∑  + , (IV.25) m =1 λ − α m λ + α m noticing that    ~ σ · K −1 ( λ)~ S K ( λ ) N ~σ0 · ~Sm 0 m m m L0 ( λ ) = ∑  + . (IV.26) m =1 λ − αm λ + αm Now it is straightforward to obtain the expression for the generating function (VI.31) in terms of the local operators     N ~Sm · ~Sn Sm · Kn (λ)Sn Kn (λ) + Kn (λ)Sn Kn (λ) · ~Sm ~ − 1 ~ − 1 ~ τ (λ) = 2 ∑ + m,n =1 (λ − αm )(λ − αn ) (λ − αm )(λ + αn )     Km−1 ( λ)~ Sm Km (λ) · Kn−1 (λ)~Sn Kn (λ) + . (IV.27) (λ + αm )(λ + αn ) It is important to notice that (IV.27) simplifies further N ~Sm · ~Sn ~Sm · ~Sn τ (λ) = 2 ∑ + (λ − αm )(λ − αn ) (λ + αm )(λ + αn ) m,n =1      ~Sm · Kn−1 (λ)~Sn Kn (λ) + Kn−1 (λ)~Sn Kn (λ) · ~Sm + . (IV.28) (λ − αm )(λ + αn ) -11-  V LINEAR BRACKET RELATIONS The Gaudin Hamiltonians with the boundary terms are obtained from the residues of the generating function (IV.28) at poles λ = ±αm : Resλ=αm τ (λ) = 4 Hm and em Resλ=−αm τ (λ) = (−4) H (IV.29) where     N ~Sm · ~Sn N ~Sm · Kn−1 (αm )~Sn Kn (αm ) + Kn−1 (αm )~Sn Kn (αm ) · ~Sm Hm = ∑ α − αn n∑ + , (IV.30) n6=m m =1 2( α m + α n ) and     ~Sm · ~Sn N N ~Sm · Kn−1 (−αm )~Sn Kn (−αm ) + Kn−1 (−αm )~Sn Kn (−αm ) · ~Sm em = ∑ α − αn n∑ H + . n6=m m =1 2( α m + α n ) (IV.31) The above Hamiltonians can be expressed in somewhat a more symmetric form     N ~ ~ − 1 ~ ~ ~ − 1 K m ( α m ) Sm K m ( α m ) · Sn + Sn · K m ( α m ) Sm K m ( α m ) Sm · ~Sn N Hm = ∑ +∑ , (IV.32) α − α n n =1 n6=m m 2( α m + α n ) and     N ~Sm · ~Sn N Km (−αm )~Sm Km −1 (− α ) · ~ m S n + ~ S n · K m (− α )~S m m mK −1 (− α ) m em = ∑ α − αn n∑ H + . n6=m m =1 2( α m + α n ) (IV.33) The next step is to study the quasi-classical limit of the exchange relations (IV.6) with the aim of deriving relevant algebraic structure for the Lax operator (IV.16). V Linear bracket relations The implementation of the algebraic Bethe ansatz requires the commutation relations between the entries of the Lax operator (IV.16). Our aim is to derive these relations as the quasi-classical limit of (IV.6). As the first step in this direction we observe that using (II.16) the reflection equation (III.1) can be expressed as

1 − ηr12 (λ − µ) K1 (λ) 1 − ηr21 (λ + µ) K2 (µ) =

(V.1) = K2 (µ) 1 − ηr12 (λ + µ) K1 (λ) 1 − ηr21 (λ − µ) The conditions obtained from the zero and first order in η are identically satisfied for the matrix K (λ). In particular, it obeys the classical reflection equation [25, 26]: r12 (λ − µ)K1 (λ)K2 (µ) + K1 (λ)r21 (λ + µ)K2 (µ) = (V.2) = K2 (µ)r12 (λ + µ)K1 (λ) + K2 (µ)K1 (λ)r21 (λ − µ). -12-  V LINEAR BRACKET RELATIONS The terms of the second order in η in (V.1) are r12 (λ − µ)K1 (λ)r21 (λ + µ)K2 (µ) = K2 (µ)r12 (λ + µ)K1 (λ)r21 (λ − µ). (V.3) This equation is also satisfied by the the K-matrix (III.4) and the classical r-matrix (II.16). In addition, the classical r-matrix (II.16) has the unitarity property r21 (−λ) = −r12 (λ), (V.4) and satisfies the classical Yang-Baxter equation [r13 (λ), r23 (µ)] + [r12 (λ − µ), r13 (λ) + r23 (µ)] = 0. (V.5) Now we can proceed to the derivation of the relevant linear bracket relations of the Lax opera- tor (IV.16). The desired relations can be obtained by writing the equation (IV.6) in the following form, using (II.16), (1 − ηr00′ (λ − µ)) T0 (λ) (1 − ηr0′ 0 (λ + µ)) T0′ (µ) = (V.6) = T0′ (µ) (1 − ηr00′ (λ + µ)) T0 (λ) (1 − ηr0′ 0 (λ − µ)) and substituting the expansion of T (λ) (IV.5) in powers of η η 2 d2 T ( λ ) T (λ) = K (λ) + η L(λ)K (λ) + |η =0 + O(η 3 ). (V.7) 2 dη 2 The zero and first orders in η are identically satisfied for the matrix K (λ) defined in (III.4). The relations we seek follow from the terms of the second order in η in (V.6). When the terms containing the second order derivatives of T are eliminated and the equation (V.3) is used to eliminate the other two terms, there are ten terms remaining. Using twice the classical reflection equation (V.2) and the unitarity property (V.4) one obtains (L0 (λ)L0′ (µ) − L0′ (µ)L0 (λ)) K0 (λ)K0′ (µ) = (r00′ (λ − µ)L0 (λ) − L0 (λ)r00′ (λ − µ)) × × K0 (λ)K0′ (µ) + (L0 (λ)K0′ (µ)r00′ (λ + µ) − K0′ (µ)r00′ (λ + µ)L0 (λ)) K0 (λ) − (r0′ 0 (µ − λ)L0′ (µ) − L0′ (µ)r0′ 0 (µ − λ)) K0 (λ)K0′ (µ) + (K0 (λ)r0′ 0 (µ + λ)L0′ (µ) − L0′ (µ)K0 (λ)r0′ 0 (µ + λ)) K0′ (µ). (V.8) Multiplying both sides of the equation (V.8) from the right by K0−1 (λ)K0−′ 1 (µ), (V.8) can be express as h i −1 [L0 (λ), L0 (µ)] = r00 (λ − µ) − K0 (µ)r00 (λ + µ)K0′ (µ), L0 (λ) ′ ′ ′ ′ h i (V.9) − r0′ 0 (µ − λ) − K0 (λ)r0′ 0 (µ + λ)K0−1 (λ), L0′ (µ) . Defining K −1 r00 ′ ( λ, µ ) = r00′ ( λ − µ ) − K0′ ( µ )r00′ ( λ + µ ) K0′ ( µ ), (V.10) -13-  VI ALGEBRAIC BETHE ANSATZ (V.9) can be written as h i h i K K [L0 (λ), L0′ (µ)] = r00 ′ ( λ, µ ), L0 ( λ ) − r0′ 0 ( µ, λ ), L0′ ( µ ) . (V.11) The commutator (V.11) is obviously anti-symmetric. It obeys the Jacobi identity because the r-matrix (V.10) satisfies the classical YB equation K K K K K [r32 (λ3 , λ2 ), r13 (λ1 , λ3 )] + [r12 (λ1 , λ2 ), r13 (λ1 , λ3 ) + r23 (λ2 , λ3 )] = 0. (V.12) The commutator (V.11) can also be recasted as an (r, s) Maillet algebra [48]. In the following we study the algebraic Bethe ansatz based on the linear bracket (V.11). VI Algebraic Bethe Ansatz Our preliminary step in the implementation of the algebraic Bethe ansatz for the open Gaudin model is to bring the boundary K-matrix to the upper, or lower, triangular form. As it was pointed out in (III.4), the general form of the K-matrix (IV.11) is   e ξ − λ ψλ e( λ ) = K . (VI.1) e φλ ξ + λ It is straightforward to check that the matrix   −1 − ν e φ U= e , (VI.2) φ −1 − ν q with ν = 1 + φ eψ e , which does not depend on the spectral parameter λ, brings the K-matrix to the upper triangular form by the similarity transformation   −1 e ξ − λν λψ K (λ) = U K ( λ )U = , (VI.3) 0 ξ + λν e + ψ. where ψ = φ e Evidently, the inverse matrix is   −1 1 ξ + λν −λψ K (λ) = 2 . (VI.4) ξ − λ 2 ν2 0 ξ − λν Direct substitution of the formulas above into (IV.25),   ! H (λ) F (λ) N ~σ0 · ~Sm K0 (λ)~σ0 K0−1 (λ) · ~Sm L0 ( λ ) = E(λ) − H (λ) = ∑ λ − αm + λ + αm , (VI.5) m =1 -14-  VI ALGEBRAIC BETHE ANSATZ yields the following local realisation for the entries of the Lax matrix N  + +  Sm (ξ + λν)Sm E(λ) = ∑ λ − αm + (ξ − λν)(λ + αm ) , (VI.6) m =1 N  − − − λ2 ψ2 S + − 2λψ( ξ − λν) S3  Sm (ξ − λν)2 Sm m m F (λ) = ∑ λ − αm + (ξ + λν)(ξ − λν)(λ + αm ) , (VI.7) m =1 N  + + ( ξ − λν) S3  S3m λψ Sm m H (λ) = ∑ λ − αm + (ξ − λν)(λ + αm ) . (VI.8) m =1 The linear bracket (V.11) based on the r-matrix r00 K ( λ, µ ) (V.10), corresponding to (VI.3), (VI.4) ′ and the classical r-matrix (II.16), defines the Lie algebra relevant for the open Gaudin model [E(λ), E(µ)] = 0, (VI.9)   2 ξ − λν [ H (λ), E(µ)] = λ E(µ) −µ E(λ) , (VI.10) λ2 − µ2 ξ − µν   2ψµ 4 ξ − µν ξ + λν [E(λ), F (µ)] = E(λ) + 2 λ H (µ) − µ H (λ) , (VI.11) (λ + µ)(ξ + µν) λ − µ2 ξ − λν ξ + µν   −ψ λ µ [ H (λ), H (µ)] = E(µ) − E(λ) , (VI.12) λ + µ ξ − λν ξ − µν     ψ 2λ ψµ2 2 ξ + λν [ H (λ), F (µ)] = H (µ) − 2 E(λ) − 2 λ F (µ) − µ F (λ) , λ+µ ξ − λν ξ − µ 2 ν2 λ − µ2 ξ + µν (VI.13)     2ψ λ µ 2ψ2 λ2 µ2 [ F (λ), F (µ)] = F (µ) − F (λ) − H (µ) − 2 H (λ) . λ + µ ξ + λν ξ + µν λ + µ ξ 2 − λ 2 ν2 ξ − µ 2 ν2 (VI.14) Our next step is to introduce the new generators e(λ), h(λ) and f (λ) as the following linear combinations of the original generators   −ξ + λν 1 ψλ 1 e( λ) = E ( λ ), h( λ ) = H (λ) − E(λ) , f (λ) = ((ξ + λν) F (λ) + ψλH (λ)) . λ λ 2ξ λ (VI.15) The key observation is that in the new basis we have [e(λ), e(µ)] = [h(λ), h(µ)] = [ f (λ), f (µ)] = 0. (VI.16) -15-  VI ALGEBRAIC BETHE ANSATZ Therefore there are only three nontrivial relations 2 [h(λ), e(µ)] = (e(µ) − e(λ)) , (VI.17) λ2 − µ2 −2 2ψν
[h(λ), f (µ)] = ( f (µ) − f (λ)) − 2 µ2 h ( µ ) − λ2 h ( λ ) λ2 −µ 2 (λ − µ )ξ 2 ψ 2
− 2 2 2 µ2 e ( µ ) − λ2 e ( λ ) , (VI.18) (λ − µ )ξ 2ψν 2 2
4 2 2 2 2 2 2
[e(λ), f (µ)] = µ e ( µ ) − λ e ( λ ) − ( ξ − µ ν ) h ( µ ) − ( ξ − λ ν ) h ( λ ) . ( λ2 − µ2 ) ξ λ2 − µ2 (VI.19) In the Hilbert space H (II.3), in every Vm = C2s+1 there exists a vector ωm ∈ Vm such that + S3m ωm = sm ωm and Sm ωm = 0. (VI.20) We define a vector Ω+ to be Ω + = ω1 ⊗ · · · ⊗ ω N ∈ H . (VI.21) From the definitions above, the formulas (VI.6) - (VI.8) and (VI.15) it is straightforward to obtain the action of the generators e(λ) and h(λ) on the vector Ω+ e(λ)Ω+ = 0, and h( λ) Ω + = ρ( λ) Ω + , (VI.22) with   1 N sm sm N 2sm λ m∑ ∑ ρ( λ) = + = 2 − α2 . (VI.23) =1 λ − αm λ + αm m =1 λ m The generating function of the Gaudin Hamiltonians (VI.31) in terms of the entries of the Lax matrix is given by τ (λ) = tr0 L20 (λ) = 2H 2 (λ) + 2F (λ) E(λ) + [ E(λ), F (λ)] . (VI.24) From (VI.11) we have that the last term is ξ 2 + λ 2 ν2 ψ [E(λ), F (λ)] = 2 H (λ) − 2H ′ (λ) + E ( λ ), (VI.25) ( ξ 2 − λ 2 ν2 ) λ ξ + λν and therefore the final expression is     2 ξ 2 + λ 2 ν2 ′ ψ τ (λ) = 2 H (λ) + 2 H (λ) − H (λ) + 2F (λ) + E ( λ ). (VI.26) ( ξ − λ 2 ν2 ) λ ξ + λν Our aim is to implement the algebraic Bethe ansatz based on the Lie algebra (VI.16) - (VI.19). To this end we need to obtain the expression for the generating function τ (λ) in terms -16-  VI ALGEBRAIC BETHE ANSATZ of the generators e(λ), h(λ) and f (λ). The first step is to invert the relations (VI.15) −λ E(λ) = e ( λ ), (VI.27) ξ − λν   ψλ H ( λ) = λ h( λ) − e( λ) , (VI.28) 2ξ (ξ − λν)   λ ψ2 λ2 F (λ) = f (λ) − ψλ h(λ) + e( λ) . (VI.29) ξ + λν 2ξ (ξ − λν) In particular, we have   2 2 2 ψλ ψ2 λ2 2 H (λ) = λ h (λ) − (2h(λ)e(λ) − [h(λ), e(λ)]) + 2 e (λ) 2ξ (ξ − λν) 4ξ (ξ − λν)2     2 2 ψλ e′ ( λ) ψ2 λ2 2 = λ h (λ) − 2h(λ)e(λ) + + 2 e (λ) . (VI.30) 2ξ (ξ − λν) λ 4ξ (ξ − λν)2 Substituting (VI.27) – (VI.30) into (VI.26) we obtain the desired expression for the generating function   2 2 2ν2 h′ ( λ) τ (λ) = 2λ h (λ) + 2 h( λ) − ξ − λ 2 ν2 λ   (VI.31) 2λ2 ψλ2 ν ψ2 λ2 ψν − 2 f (λ) + h( λ) + e( λ) − e ( λ ). ξ − λ 2 ν2 ξ 4ξ 2 ξ An important initial observation in the implementation of the algebraic Bethe ansatz is that the vector Ω+ (VI.21) is an eigenvector of the generating function τ (λ), to show this we use the expression above, (VI.22) and (VI.23)   2 2 2ν2 ρ(λ) ρ′ ( λ) τ (λ)Ω+ = χ0 (λ)Ω+ = 2λ ρ (λ) + 2 − Ω+ , (VI.32) ξ − λ 2 ν2 λ using (VI.23) the eigenvalue χ0 (λ) can be expressed as !! N N N s m ( s m + 1) sm ν2 2sn χ0 (λ) = 8λ2 ∑ ( λ2 − α2 )2 ∑ λ2 − α2 + 2 2 2 + ∑ . (VI.33) m =1 m m =1 m ξ −λ ν α − α2n 2 n6=m m An essential step in the algebraic Bethe ansatz is a definition of the corresponding Bethe vectors. In this case, they are symmetric functions of their arguments and are such that the off shell action of the generating function of the Gaudin Hamiltonians is as simple as possible. With this aim we attempt to show that the Bethe vector ϕ1 (µ) has the form ϕ1 (µ) = ( f (µ) + c1 (µ)) Ω+ , (VI.34) -17-  VI ALGEBRAIC BETHE ANSATZ where c1 (µ) is given by ψν
1 − µ2 ρ ( µ ) . c1 ( µ ) = − (VI.35) ξ Evidently, the action of the generating function of the Gaudin Hamiltonians reads τ (λ) ϕ1 (µ) = [ τ (λ), f (µ)] Ω+ + χ0 (λ) ϕ1 (µ). (VI.36) A straightforward calculation show that the commutator in the first term of (VI.36) is given by   8λ2 ν2 [τ (λ), f (µ)] Ω+ = − 2 ρ( λ) + 2 ϕ1 ( µ ) λ − µ2 ξ − λ 2 ν2   8λ2 (ξ 2 − µ2 ν2 ) ν2 + 2 ρ( µ) + 2 ϕ1 ( λ ). (VI.37) (λ − µ2 )(ξ 2 − λ2 ν2 ) ξ − µ 2 ν2 Therefore the action of the generating function τ (λ) on ϕ1 (µ) is given by   8λ2 (ξ 2 − µ2 ν2 ) ν2 τ (λ) ϕ1 (µ) = χ1 (λ, µ) ϕ1 (µ) + 2 ρ( µ) + 2 ϕ1 ( λ ), (VI.38) (λ − µ2 )(ξ 2 − λ2 ν2 ) ξ − µ 2 ν2 with   8λ2 ν2 χ1 (λ, µ) = χ0 (λ) − 2 ρ( λ) + 2 . (VI.39) λ − µ2 ξ − λ 2 ν2 The unwanted term in (VI.38) vanishes when the following Bethe equation is imposed on the parameter µ, ν2 ρ( µ) + 2 = 0. (VI.40) ξ − µ 2 ν2 Thus we have shown that ϕ1 (µ) (VI.34) is the desired Bethe vector of the generating function τ (λ) corresponding to the eigenvalue χ1 (λ, µ). We seek the Bethe vector ϕ2 (µ1 , µ2 ) as the following symmetric function ( 1) ( 1) ( 2) ϕ 2 ( µ1 , µ2 ) = f ( µ1 ) f ( µ2 ) Ω + + c2 ( µ2 ; µ1 ) f ( µ1 ) Ω + + c2 ( µ1 ; µ2 ) f ( µ2 ) Ω + + c2 ( µ1 , µ2 ) Ω + , (VI.41) ( 1) ( 2) where the scalar coefficients c2 (µ1 ; µ2 ) and c2 (µ1 , µ2 ) are   ( 1) ψν 2 2µ21 c2 ( µ1 ; µ2 ) = − 1 − µ1 ρ ( µ1 ) + 2 , (VI.42) ξ µ1 − µ22   ( 2) ψ2 (ξ 2 − 3µ22 ν2 )ρ(µ1 ) − (ξ 2 − 3µ21 ν2 )ρ(µ2 ) 2 2 2 2 c2 ( µ1 , µ2 ) = − 2 + ( ξ − ( µ1 + µ2 ) ν ) ρ ( µ1 ) ρ ( µ2 ) . ν µ21 − µ22 (VI.43) One way to obtain the action of τ (λ) on ϕ2 (µ1 , µ2 ) is to write   ( 1) τ (λ) ϕ2 (µ1 , µ2 ) = [[ τ (λ), f (µ1 )] , f (µ2 )] Ω+ + f (µ2 ) + c2 (µ2 ; µ1 ) [ τ (λ), f (µ1 )] Ω+   (VI.44) ( 1) + f (µ1 ) + c2 (µ1 ; µ2 ) [τ (λ), f (µ2 )] Ω+ + χ0 (λ) ϕ2 (µ1 , µ2 ). -18-  VI ALGEBRAIC BETHE ANSATZ Then to substitute (VI.37) in the second and third term above and use the relation   ψν 2µ22 ( 1) f ( µ1 ) + c2 ( µ1 ; µ2 ) ϕ 1 ( µ2 ) = ϕ 2 ( µ1 , µ2 ) − ϕ 1 ( µ1 ) ξ µ21 − µ22   ( 2) ψν µ21 c1 (µ2 ) − µ22 c1 (µ1 ) − c2 ( µ1 , µ2 ) − c1 ( µ1 ) c1 ( µ2 ) + 2 Ω+ , ξ µ21 − µ22 (VI.45) which follows from the definition (VI.41). A straightforward calculation shows that the off shell action of the generating function τ (λ) on ϕ2 (µ1 , µ2 ) is given by 8λ2 (ξ 2 − µ2i ν2 ) 2 τ (λ) ϕ2 (µ1 , µ2 ) = χ2 (λ, µ1 , µ2 ) ϕ2 (µ1 , µ2 ) + ∑ 2 2 2 2 2 × i=1 ( λ − µ i )( ξ − λ ν ) ! (VI.46) ν2 2 × ρ( µi ) + 2 − 2 ϕ2 (λ, µ3−i ), ξ − µ2i ν2 µi − µ23−i with the eigenvalue ! 2 8λ2 ν2 1 χ2 (λ, µ1 , µ2 ) = χ0 (λ) − ∑ 2 2 ρ( λ) + 2 2 2 − 2 . (VI.47) i=1 λ − µ i ξ −λ ν λ − µ23−i The two unwanted terms in the action above (VI.46) vanish when the Bethe equations are imposed on the parameters µ1 and µ2 , ν2 2 ρ( µi ) + 2 − 2 = 0, (VI.48) 2 ξ − µi ν 2 µi − µ23−i with i = 1, 2. Therefore ϕ2 (µ1 , µ2 ) is the Bethe vector of the generating function of the Gaudin Hamiltonians with the eigenvalue χ2 (λ, µ1 , µ2 ). As our next step we propose the Bethe vector ϕ3 (µ1 , µ2 , µ3 ) in the form of the following symmetric function of its arguments ( 1) ( 1) ϕ 3 ( µ1 , µ2 , µ3 ) = f ( µ1 ) f ( µ2 ) f ( µ3 ) Ω + + c3 ( µ1 ; µ2 , µ3 ) f ( µ2 ) f ( µ3 ) Ω + + c3 ( µ2 ; µ3 , µ1 ) f ( µ3 ) f ( µ1 ) Ω + ( 1) ( 2) ( 2) + c3 ( µ3 ; µ1 , µ2 ) f ( µ1 ) f ( µ2 ) Ω + + c3 ( µ1 , µ2 ; µ3 ) f ( µ3 ) Ω + + c3 ( µ2 , µ3 ; µ1 ) f ( µ1 ) Ω + ( 2) ( 3) + c3 ( µ3 , µ1 ; µ2 ) f ( µ2 ) Ω + + c3 ( µ1 , µ2 , µ3 ) Ω + , (VI.49) where the three scalar coefficients above are given by -19-  VI ALGEBRAIC BETHE ANSATZ   ( 1) ψν 2µ2 2µ2 c3 ( µ1 ; µ2 , µ3 ) =− 1 − µ21 ρ(µ1 ) + 2 1 2+ 2 1 2 , (VI.50) ξ µ1 − µ2 µ1 − µ3      ( 2) ψ2 ξ 2 − 3µ22 ν2 2 ξ 2 − 3µ21 ν2 2 c3 ( µ1 , µ2 ; µ3 ) =− 2 ρ ( µ1 ) − 2 − ρ ( µ2 ) − 2 ν µ21 − µ22 µ1 − µ23 µ21 − µ22 µ2 − µ23    ψ2 2 2 − 2 (ξ 2 − (µ21 + µ22 )ν2 ) ρ(µ1 ) − 2 ρ ( µ 2 ) − , (VI.51) ν µ1 − µ23 µ22 − µ23

( 3) ψ3 4ξ 4 + ξ 2 + µ21 ν2 4µ21 − 5(µ22 + µ23 ) ν2 c3 ( µ1 , µ2 , µ3 ) =− 3 ρ ( µ1 ) ν ξ (µ21 − µ22 )(µ21 − µ23 )



! 4ξ 4 + ξ 2 + µ22 ν2 4µ22 − 5(µ23 + µ21 ) ν2 4ξ 4 + ξ 2 + µ23 ν2 4µ23 − 5(µ21 + µ22 ) ν2 + ρ ( µ2 ) + ρ ( µ3 ) (µ22 − µ23 )(µ22 − µ21 ) (µ23 − µ21 )(µ23 − µ22 )



ψ3 ξ 2 ν2 µ41 + µ42 − µ21 µ22 + 2µ23 µ21 + µ22 − 5µ43 − 2ξ 4 − µ21 µ22 ν4 µ21 + µ22 − 2µ23 − 3 ρ ( µ1 ) ρ ( µ2 ) ν ξ (µ21 − µ23 )(µ22 − µ23 )



ξ 2 ν2 µ42 + µ43 − µ22 µ23 + 2µ21 µ22 + µ23 − 5µ41 − 2ξ 4 − µ22 µ23 ν4 µ22 + µ23 − 2µ21 + ρ ( µ2 ) ρ ( µ3 ) (µ22 − µ21 )(µ23 − µ21 )



! ξ 2 ν2 µ43 + µ41 − µ23 µ21 + 2µ22 µ23 + µ21 − 5µ42 − 2ξ 4 − µ23 µ21 ν4 µ23 + µ21 − 2µ22 + ρ ( µ3 ) ρ ( µ1 ) (µ23 − µ22 )(µ21 − µ22 ) ψ3

− 3 ξ 2ξ 2 − µ21 + µ22 + µ23 ν2 ρ(µ1 )ρ(µ2 )ρ(µ3 ). (VI.52) ν A lengthy but straightforward calculation based on appropriate generalisation of (VI.44) and (VI.45) shows that the action of the generating function τ (λ) on ϕ3 (µ1 , µ2 , µ3 ) is given by 3 8λ2 (ξ 2 − µ2i ν2 ) τ (λ) ϕ3 (µ1 , µ2 , µ3 ) = χ3 (λ, µ1 , µ2 µ3 ) ϕ3 (µ1 , µ2 , µ3 ) + ∑ 2 2 2 2 2 × i=1 ( λ − µ i )( ξ − λ ν ) ! (VI.53) 3 ν2 2 × ρ( µi ) + 2 2 ξ − µi ν 2 − ∑ 2 2 ϕ3 (λ, {µ j } j6=i ), j6 = i µ i − µ j where the eigenvalue is ! 3 3 8λ2 ν2 1 χ3 (λ, µ1 , µ2 , µ3 ) = χ0 (λ) − ∑ 2 2 ρ( λ) + 2 2 2 − ∑ 2 . (VI.54) i=1 λ − µ i ξ −λ ν 2 j6 = i λ − µ j The three unwanted terms in (VI.53) vanish when the Bethe equation are imposed on the pa- rameters µi , 3 ν2 2 ρ( µi ) + 2 2 ξ − µi ν 2 − ∑ 2 2 = 0, (VI.55) j6 = i µ i − µ j -20-  VI ALGEBRAIC BETHE ANSATZ with i = 1, 2, 3. As a symmetric function of its arguments the Bethe vector ϕ4 (µ1 , µ2 , µ3 , µ4 ) is given explic- itly in the Appendix B. It is possible to check that the off shell action of the generating function τ (λ) on the Bethe vector ϕ4 (µ1 , µ2 , µ3 , µ4 ) is given by 4 8λ2 (ξ 2 − µ2i ν2 ) τ (λ) ϕ4 (µ1 , µ2 , µ3 , µ4 ) = χ4 (λ, µ1 , µ2 µ3 ) ϕ4 (µ1 , µ2 , µ3 , µ4 ) + ∑ 2 2 2 2 2 × i=1 ( λ − µ i )( ξ − λ ν ) ! (VI.56) 4 ν2 2 × ρ( µi ) + − ∑ ξ 2 − µ2i ν2 j6=i µ2i − µ2j ϕ4 (λ, {µ j } j6=i ), with the eigenvalue ! 4 4 8λ2 ν2 1 χ4 (λ, µ1 , µ2 , µ3 , µ4 ) = χ0 (λ) − ∑ 2 2 ρ( λ) + 2 2 2 − ∑ 2 . (VI.57) i=1 λ − µ i ξ −λ ν 2 j6 = i λ − µ j The four unwanted terms on the right hand side of (VI.56) vanish when the Bethe equation are imposed on the parameters µi , 4 ν2 2 ρ( µi ) + − ∑ ξ 2 − µi ν2 j6=i µi − µ2j 2 2 = 0, (VI.58) with i = 1, 2, 3, 4. Based on the results presented above we can conclude that the local realisation (VI.6) - (VI.8) of the Lie algebra (VI.15) - (VI.19) yields the spectrum χ M (λ, µ1 , . . . , µ M ) of the generat- ing function of the Gaudin Hamiltonians ! M M 8λ2 ν2 1 χ M (λ, µ1 , . . . , µ M ) = χ0 (λ) − ∑ 2 ξ − λ 2 ν2 ∑ 2 ρ( λ) + 2 − 2 2 , (VI.59) i=1 λ − µ i j6 = i λ − µ j and the corresponding Bethe equations which should be imposed on the parameters µi M ν2 2 ρ( µi ) + − ∑ ξ 2 − µi ν2 j6=i µi − µ2j 2 2 = 0, (VI.60) where i = 1, 2, . . . , M. Moreover, from (IV.29) and (VI.59) it follows that the eigenvalues of the Gaudin Hamiltonians (IV.32) and (IV.33) can be obtained as the residues of χ M (λ, µ1 , . . . , µ M ) at poles λ = ±αm ! N 1 s m ( s m + 1) ν2 2sn ξ 2 − α2m ν2 n∑ Em = Resλ=αm χ M (λ, µ1 , . . . , µ M ) = + αm sm + 4 2αm α2 − α2n 6=m m M 1 − 2αm sm ∑ , 2 i=1 α m − µ2i (VI.61) -21-  VII CONCLUSION and ! N 1 s ( s + 1) ν2 2sn Eem = − Resλ=−αm χ M (λ, µ1 , . . . , µ M ) = m m + αm sm 2 2 2 + ∑ 4 2αm ξ − αm ν 2 α − α2n n6=m m M 1 − 2αm sm ∑ . 2 i=1 α m − µ2i (VI.62) Evidently, the respective eigenvalues (VI.61) and (VI.62) of the Hamiltonians (IV.32) and (IV.33) coincide. When all the spin sm are set to one half, these energies coincide with the expressions obtained in [42] (up to normalisation). The Bethe equations are also equivalent, the corre- spondence between our variables and the one used in [42] being given by (the left hand sides correspond to our variables, the left hand sides to the ones used in [42]): λj θm ξ ξ µj = ; αm = ; = . (VI.63) 1 − ξ ( 1) 1 − ξ ( 1) ν 1 − ξ ( 1) However, explicit and compact form of the relevant Bethe vector ϕ M (µ1 , µ2 , . . . , µ M ), for an arbitrary positive integer M, requires further studies and will be reported elsewhere. As it is evident form the formulas for the Bethe vector ϕ4 (µ1 , µ2 , µ3 , µ4 ) given in the Appendix (m) B, the main problem lies in the definition the scalar coefficients c M (µ1 , . . . µm ; µm+1 , . . . , µ M ), with m = 1, 2, . . . , M. Some of them can be obtained straightforwardly, but, in particular, the ( M) coefficient c M (µ1 , µ2 . . . , µ M ) still represents a challenge, at least in the present form of the Bethe vectors. VII Conclusion Following Sklyanin’s proposal in the periodic case [3], here we have derived the generating function of the Gaudin Hamiltonians with boundary terms. Our derivation is based on the quasi-classical expansion of the linear combination of the transfer matrix of the XXX Heisen- berg spin chain and the central element, the so-called Sklyanin determinant. The correspond- ing Gaudin Hamiltonians with boundary terms are obtained as the residues of the generating function. Then we have studied the appropriate algebraic structure, including the classical reflection equation. Our approach to the algebraic Bethe ansatz is based on the relevant Lax matrix which satisfies certain linear bracket and simultaneously provides the local realisation for the corresponding Lie algebra. By defining the appropriate Bethe vectors we have obtained the strikingly simple off shell action of the generating function of the Gaudin Hamiltonians. Actually, the action of the generating function is as simple as it could possible be since it al- most coincides with the one in the case when the boundary matrix is diagonal [20]. In this way we have implemented the algebraic Bethe ansatz, obtaining the spectrum of the generating function and the corresponding Bethe equations. -22-  A BASIC DEFINITIONS Although the off shell action of the generating function which we have established is very simple, it would be important to obtain more compact formula for the Bethe vector ϕ M (µ1 , µ2 , . . . , µ M ), for an arbitrary positive integer M. In particular, simpler expression for (m) the scalar coefficients c M (µ1 , . . . µm ; µm+1 , . . . , µ M ), with m = 1, 2, . . . , M would be of utmost importance. Such a formula would be crucial for the off shell scalar product of the Bethe vec- tors and these results could lead to the correlations functions of Gaudin model with boundary. Moreover, it would be of considerable interest to establish a relation between Bethe vectors and solutions of the corresponding Knizhnik-Zamolodchikov equations, along the lines it was done in the case when the boundary matrix is diagonal [20]. Acknowledgments We acknowledge useful discussions with Zoltán Nagy. E.R. would like to thank to the staff of the GFM-UL and the Department of Mathematics of the University of the Algarve for warm hospitality while a part of the work was done. I. S. was supported in part by the Ser- bian Ministry of Science and Technological Development under grant number ON 171031. N. M. is thankful to Professor Victor Kac and the staff of the Mathematics Department at MIT for their kind hospitality. N. M. was supported in part by the FCT sabbatical fellowship SFRH/BSAB/1366/2013. A Basic definitions We consider the spin operators Sα with α = +, −, 3, acting in some (spin s) representation space C2s+1 with the commutation relations [S3 , S ± ] = ± S ± , [S+ , S− ] = 2S3 , (A.1) and Casimir operator 1 c2 = (S3 )2 + (S+ S− + S− S+ ) = (S3 )2 + S3 + S− S+ = ~S · ~S. 2 1 In the particular case of spin 2 representation, one recovers the Pauli matrices   1 1 δα3 2δα+ S = σα = α . 2 2 2δα− −δα3 We consider a spin chain with N sites with spin s representations, i.e. a local C2s+1 space at each site and the operators α Sm = 1 ⊗ · · · ⊗ |{z} Sα ⊗ · · · ⊗ 1, (A.2) m with α = +, −, 3 and m = 1, 2, . . . , N. -23-  B BETHE VECTOR ϕ4 (µ1 , µ2 , µ3 , µ4 ) B Bethe vector ϕ4 (µ1 , µ2 , µ3 , µ4 ) Here we present explicit formulas of the Bethe vector ϕ4 (µ1 , µ2 , µ3 , µ4 ). The vector ϕ4 (µ1 , µ2 , µ3 , µ4 ) is an symmetric function of its arguments and is given by ( 1) ϕ 4 ( µ1 , µ2 , µ3 , µ4 ) = f ( µ1 ) f ( µ2 ) f ( µ3 ) f ( µ4 ) Ω + + c4 ( µ4 ; µ1 , µ2 , µ3 ) f ( µ1 ) f ( µ2 ) f ( µ3 ) Ω + ( 1) ( 1) + c4 ( µ3 ; µ1 , µ2 , µ4 ) f ( µ1 ) f ( µ2 ) f ( µ4 ) Ω + + c4 ( µ2 ; µ1 , µ3 , µ4 , ) f ( µ1 ) f ( µ3 ) f ( µ4 ) Ω + ( 1) ( 2) + c4 ( µ1 ; µ2 , µ3 , µ4 ) f ( µ2 ) f ( µ3 ) f ( µ4 ) Ω + + c4 ( µ3 , µ4 ; µ1 , µ2 ) f ( µ1 ) f ( µ2 ) Ω + ( 2) ( 2) + c4 ( µ2 , µ4 ; µ1 , µ3 ) f ( µ1 ) f ( µ3 ) Ω + + c4 ( µ2 , µ3 ; µ1 , µ4 ) f ( µ1 ) f ( µ4 ) Ω + ( 2) ( 2) + c4 ( µ1 , µ4 ; µ2 , µ3 ) f ( µ2 ) f ( µ3 ) Ω + + c4 ( µ1 , µ3 ; µ2 , µ4 ) f ( µ2 ) f ( µ4 ) Ω + ( 2) ( 3) + c4 ( µ1 , µ2 ; µ3 , µ4 ) f ( µ3 ) f ( µ4 ) Ω + + c4 ( µ2 , µ3 , µ4 ; µ1 ) f ( µ1 ) Ω + ( 3) ( 3) + c4 ( µ1 , µ2 , µ4 ; µ2 ) f ( µ2 ) Ω + + c4 ( µ1 , µ2 , µ4 ; µ3 ) f ( µ3 ) Ω + ( 3) ( 4) + c4 ( µ1 , µ2 , µ3 ; µ4 ) f ( µ4 ) Ω + + c4 ( µ1 , µ2 , µ3 , µ4 ) Ω + , (B.1) where the four scalar coefficients are ! 4 ( 1) ψν 2µ2 c4 ( µ1 ; µ2 , µ3 , µ4 ) = − 1 − µ21 ρ(µ1 ) + ∑ 2 1 2 , (B.2) i =2 µ1 − µ i ξ ! !! 4 4 ( 2) ψ2 ξ 2 − 3µ22 ν2 2 ξ 2 − 3µ21 ν2 2 c4 ( µ1 , µ2 ; µ3 , µ4 ) =− 2 ρ ( µ1 ) − ∑ 2 − ρ ( µ2 ) − ∑ 2 ν µ21 − µ22 2 i =3 µ1 − µ i µ21 − µ22 2 j =3 µ2 − µ j ! ! 4 4 ψ2 2 2 2 − (ξ − (µ21 + µ22 )ν2 ) ρ ( µ1 ) − ∑ ρ ( µ2 ) − ∑ 2 , ν2 2 i =3 µ1 − µ2i j =3 µ2 − µ j 2 (B.3)

  ( 3) ψ3 4ξ 4 + ξ 2 + µ21 ν2 4µ21 − 5(µ22 + µ23 ) ν2 2 c4 ( µ1 , µ2 , µ3 ; µ4 ) =− 3 ρ ( µ1 ) − 2 ν ξ (µ21 − µ22 )(µ21 − µ23 ) µ1 − µ24

  4ξ 4 + ξ 2 + µ22 ν2 4µ22 − 5(µ23 + µ21 ) ν2 2 + ρ ( µ2 ) − 2 (µ22 − µ23 )(µ22 − µ21 ) µ2 − µ24

 ! 4ξ 4 + ξ 2 + µ23 ν2 4µ23 − 5(µ21 + µ22 ) ν2 2 + ρ ( µ3 ) − 2 (µ23 − µ21 )(µ23 − µ22 ) µ3 − µ24 -24-  B BETHE VECTOR ϕ4 (µ1 , µ2 , µ3 , µ4 )



ψ3 ξ 2 ν2 µ41 + µ42 − µ21 µ22 + 2µ23 µ21 + µ22 − 5µ43 − 2ξ 4 − µ21 µ22 ν4 µ21 + µ22 − 2µ23 − 3 × ν ξ (µ21 − µ23 )(µ22 − µ23 )    2 2 × ρ ( µ1 ) − 2 ρ ( µ2 ) − 2 µ1 − µ24 µ2 − µ24



ξ 2 ν2 µ42 + µ43 − µ22 µ23 + 2µ21 µ22 + µ23 − 5µ41 − 2ξ 4 − µ22 µ23 ν4 µ22 + µ23 − 2µ21 + × (µ22 − µ21 )(µ23 − µ21 )    2 2 × ρ ( µ2 ) − 2 ρ ( µ3 ) − 2 µ2 − µ24 µ3 − µ24



ξ 2 ν2 µ43 + µ41 − µ23 µ21 + 2µ22 µ23 + µ21 − 5µ42 − 2ξ 4 − µ23 µ21 ν4 µ23 + µ21 − 2µ22 + × (µ23 − µ22 )(µ21 − µ22 )    2 2 × ρ ( µ3 ) − 2 ρ ( µ 1 ) − µ3 − µ24 µ21 − µ24     ψ3

2 2 2 − 3 ξ 2ξ 2 − µ21 + µ22 + µ23 ν2 ρ(µ1 ) − 2 ρ ( µ 2 ) − ρ ( µ 3 ) − . ν µ1 − µ24 µ22 − µ24 µ23 − µ24 (B.4)

( 4) 2ψ49ξ 4 + ξ 2 ν2 27µ21 − 7(µ22 + µ23 + µ24 ) + 3µ21 ν4 8µ21 − 7(µ22 + µ23 + µ24 ) c4 ( µ1 , µ2 , µ3 , µ4 ) =− 4 ρ ( µ1 ) ν (µ21 − µ22 )(µ21 − µ23 )(µ21 − µ24 )

9ξ 4 + ξ 2 ν2 27µ22 − 7(µ21 + µ23 + µ24 ) + 3µ22 ν4 8µ22 − 7(µ22 + µ23 + µ24 ) + ρ ( µ2 ) (µ22 − µ21 )(µ22 − µ23 )(µ22 − µ24 )

9ξ 4 + ξ 2 ν2 27µ23 − 7(µ21 + µ22 + µ24 ) + 3µ23 ν4 8µ23 − 7(µ21 + µ22 + µ24 ) + ρ ( µ3 ) (µ23 − µ21 )(µ23 − µ22 )(µ23 − µ24 )

! 9ξ 4 + ξ 2 ν2 27µ24 − 7(µ21 + µ22 + µ23 ) + 3µ24 ν4 8µ24 − 7(µ21 + µ22 + µ23 ) + ρ ( µ4 ) (µ24 − µ21 )(µ24 − µ22 )(µ24 − µ23 ) ψ4  4  4 2 2 4 2 2 2 2 2 2 2 2 2 2  − 3ξ 2 ( µ 1 + µ µ 1 2 + µ 2 ) − 3 ( µ µ 1 3 + µ µ 1 4 + µ µ 2 3 + µ µ 2 4 − 2µ µ 3 4 ) − 18ξ 2 ν2 µ21 µ22 (µ23 + µ24 ) ν4  + ξ 2 ν2 −2(µ61 − 6µ41 µ22 − 6µ21 µ42 + µ62 ) − 4(µ21 + µ22 )2 (µ23 + µ24 ) + 7(µ21 + µ22 )(µ23 + µ24 )2 + 2µ23 µ24 ×
  ×(5(µ21 + µ22 ) − 7(µ23 + µ24 )) + ν4 14µ21 µ22 (µ43 + µ44 ) − (4µ21 µ22 + 7µ23 µ24 )(µ21 + µ22 )(µ23 + µ24 ) + 2ν4 ×

ρ ( µ1 ) ρ ( µ2 ) × 4µ23 µ24 (µ21 + µ22 )2 − 3µ21 µ22 (µ21 − µ22 )2 − µ21 µ22 (3µ21 µ22 + 5µ23 µ24 ) (µ21 − µ23 )(µ21 − µ24 )(µ22 − µ23 )(µ22 − µ24 ) -25-  B BETHE VECTOR ϕ4 (µ1 , µ2 , µ3 , µ4 ) ψ4  4  4 2 2 4 2 2 2 2 2 2 2 2 2 2  − 3ξ 2 ( µ 1 + µ µ 1 3 + µ 3 ) − 3 ( µ µ 1 2 + µ µ 1 4 + µ µ 2 3 + µ µ 3 4 − 2µ µ 2 4 ) − 18ξ 2 ν2 µ21 µ23 (µ22 + µ24 ) ν4  + ξ 2 ν2 −2(µ61 − 6µ41 µ23 − 6µ21 µ43 + µ63 ) − 4(µ21 + µ23 )2 (µ22 + µ24 ) + 7(µ21 + µ23 )(µ22 + µ24 )2 + 2µ22 µ24 ×
  ×(5(µ21 + µ23 ) − 7(µ22 + µ24 )) + ν4 14µ21 µ23 (µ42 + µ44 ) − (4µ21 µ23 + 7µ22 µ24 )(µ21 + µ23 )(µ22 + µ24 ) + 2ν4 ×

ρ ( µ1 ) ρ ( µ3 ) × 4µ22 µ24 (µ21 + µ23 )2 − 3µ21 µ23 (µ21 − µ23 )2 − µ21 µ23 (3µ21 µ23 + 5µ22 µ24 ) (µ21 − µ22 )(µ21 − µ24 )(µ23 − µ22 )(µ23 − µ24 ) ψ4  4  4 2 2 4 2 2 2 2 2 2 2 2 2 2  − 3ξ 2 ( µ 1 + µ µ 1 4 + µ 4 ) − 3 ( µ µ 1 2 + µ µ 1 3 + µ µ 2 4 + µ µ 3 4 − 2µ µ 2 3 ) − 18ξ 2 ν2 µ21 µ24 (µ22 + µ23 ) ν4  + ξ 2 ν2 −2(µ61 − 6µ41 µ24 − 6µ21 µ44 + µ64 ) − 4(µ21 + µ24 )2 (µ22 + µ23 ) + 7(µ21 + µ24 )(µ22 + µ23 )2 + 2µ22 µ23 ×
  ×(5(µ21 + µ24 ) − 7(µ22 + µ23 )) + ν4 14µ21 µ24 (µ42 + µ43 ) − (4µ21 µ24 + 7µ22 µ23 )(µ21 + µ24 )(µ22 + µ23 ) + 2ν4 ×

ρ ( µ1 ) ρ ( µ4 ) × 4µ22 µ23 (µ21 + µ24 )2 − 3µ21 µ24 (µ21 − µ24 )2 − µ21 µ24 (3µ21 µ24 + 5µ22 µ23 ) (µ21 − µ22 )(µ21 − µ23 )(µ24 − µ22 )(µ24 − µ23 ) ψ4  4  4 2 2 4 2 2 2 2 2 2 2 2 2 2  − 3ξ 2 ( µ 2 + µ µ 2 3 + µ 3 ) − 3 ( µ µ 1 2 + µ µ 2 4 + µ µ 1 3 + µ µ 3 4 − 2µ µ 1 4 ) − 18ξ 2 ν2 µ22 µ23 (µ21 + µ24 ) ν4  + ξ 2 ν2 −2(µ62 − 6µ42 µ23 − 6µ22 µ43 + µ63 ) − 4(µ22 + µ23 )2 (µ21 + µ24 ) + 7(µ22 + µ23 )(µ21 + µ24 )2 + 2µ21 µ24 ×
  ×(5(µ22 + µ23 ) − 7(µ21 + µ24 )) + ν4 14µ22 µ23 (µ41 + µ44 ) − (4µ22 µ23 + 7µ21 µ24 )(µ22 + µ23 )(µ21 + µ24 ) + 2ν4 ×

ρ ( µ2 ) ρ ( µ3 ) × 4µ21 µ24 (µ22 + µ23 )2 − 3µ22 µ23 (µ22 − µ23 )2 − µ22 µ23 (3µ22 µ23 + 5µ21 µ24 ) (µ22 − µ21 )(µ22 − µ24 )(µ23 − µ21 )(µ23 − µ24 ) ψ4  4  4 2 2 4 2 2 2 2 2 2 2 2 2 2  − 3ξ 2 ( µ 2 + µ µ 2 4 + µ 4 ) − 3 ( µ µ 1 2 + µ µ 2 3 + µ µ 1 4 + µ µ 3 4 − 2µ µ 1 3 ) − 18ξ 2 ν2 µ22 µ24 (µ21 + µ23 ) ν4  + ξ 2 ν2 −2(µ62 − 6µ42 µ24 − 6µ22 µ44 + µ64 ) − 4(µ22 + µ24 )2 (µ21 + µ23 ) + 7(µ22 + µ24 )(µ21 + µ23 )2 + 2µ21 µ23 ×
  ×(5(µ22 + µ24 ) − 7(µ21 + µ23 )) + ν4 14µ22 µ24 (µ41 + µ43 ) − (4µ22 µ24 + 7µ21 µ23 )(µ22 + µ24 )(µ21 + µ23 ) + 2ν4 ×

ρ ( µ2 ) ρ ( µ4 ) × 4µ21 µ23 (µ22 + µ24 )2 − 3µ22 µ24 (µ22 − µ24 )2 − µ22 µ24 (3µ22 µ24 + 5µ21 µ23 ) (µ22 − µ21 )(µ22 − µ23 )(µ24 − µ21 )(µ24 − µ23 ) ψ4  4  4 2 2 4 2 2 2 2 2 2 2 2 2 2  − 3ξ 2 ( µ 3 + µ µ 3 4 + µ 4 ) − 3 ( µ µ 2 3 + µ µ 1 3 + µ µ 2 4 + µ µ 1 4 − 2µ µ 1 2 ) − 18ξ 2 ν2 µ23 µ24 (µ21 + µ22 ) ν4  + ξ 2 ν2 −2(µ63 − 6µ43 µ24 − 6µ23 µ44 + µ64 ) − 4(µ23 + µ24 )2 (µ21 + µ22 ) + 7(µ23 + µ24 )(µ21 + µ22 )2 + 2µ21 µ22 ×
  ×(5(µ23 + µ24 ) − 7(µ21 + µ22 )) + ν4 14µ23 µ24 (µ41 + µ42 ) − (4µ23 µ24 + 7µ21 µ22 )(µ23 + µ24 )(µ21 + µ22 ) + 2ν4 ×

ρ ( µ3 ) ρ ( µ4 ) × 4µ21 µ22 (µ23 + µ24 )2 − 3µ23 µ24 (µ23 − µ24 )2 − µ23 µ24 (3µ23 µ24 + 5µ21 µ22 ) (µ23 − µ21 )(µ23 − µ22 )(µ24 − µ21 )(µ24 − µ22 ) -26-  B BETHE VECTOR ϕ4 (µ1 , µ2 , µ3 , µ4 ) ψ4  4  4 2 2 2 2 2 2 2 2 2 2  2 2  − 3ξ 3µ 4 − 2µ ( 4 1µ + µ 2 + µ 3 ) + µ µ 1 2 + µ µ 1 3 + µ µ 2 3 − ξ ν 7µ64 − 4µ44 (µ21 + µ22 + µ23 ) + µ24 × ν4  ×(3(µ21 µ22 + µ21 µ23 + µ22 µ23 ) − 2(µ41 + µ42 + µ43 )) + µ41 µ22 + µ41 µ23 + µ21 µ42 + µ21 µ43 + µ42 µ23 + µ22 µ43 − 4µ21 µ22 µ23   − ν4 2µ44 (µ21 µ22 + µ21 µ23 + µ22 µ23 ) − µ24 (µ41 µ22 + µ21 µ42 + µ41 µ23 + µ21 µ43 + µ42 µ23 + µ22 µ43 + 6µ21 µ22 µ23 )
 ρ ( µ1 ) ρ ( µ2 ) ρ ( µ3 ) −ν4 2µ21 µ22 µ23 (µ21 + µ22 + µ23 ) (µ1 − µ24 )(µ22 − µ24 )(µ23 − µ24 ) 2 ψ4  4  4 2 2 2 2 2 2 2 2 2 2  2 2  − 3ξ 3µ 3 − 2µ ( 3 1µ + µ 2 + µ 4 ) + µ µ 1 2 + µ µ 1 4 + µ µ 2 4 − ξ ν 7µ63 − 4µ43 (µ21 + µ22 + µ24 ) + µ23 × ν4  ×(3(µ21 µ22 + µ21 µ24 + µ22 µ24 ) − 2(µ41 + µ42 + µ44 )) + µ41 µ22 + µ41 µ24 + µ21 µ42 + µ21 µ44 + µ42 µ24 + µ22 µ44 − 4µ21 µ22 µ24   − ν4 2µ43 (µ21 µ22 + µ21 µ24 + µ22 µ24 ) − µ23 (µ41 µ22 + µ21 µ42 + µ41 µ24 + µ21 µ44 + µ42 µ24 + µ22 µ44 + 6µ21 µ22 µ24 )
 ρ ( µ1 ) ρ ( µ2 ) ρ ( µ4 ) −ν4 2µ21 µ22 µ24 (µ21 + µ22 + µ24 ) (µ1 − µ23 )(µ22 − µ23 )(µ24 − µ23 ) 2 ψ4  4  4 2 2 2 2 2 2 3 2 2 2  2 2  − 3ξ 3µ 2 − 2µ ( 2 1µ + µ 3 + µ 4 ) + µ µ 1 3 + µ µ 1 4 + µ µ 3 4 − ξ ν 7µ62 − 4µ42 (µ21 + µ23 + µ24 ) + µ22 × ν4  ×(3(µ21 µ23 + µ21 µ24 + µ23 µ24 ) − 2(µ41 + µ43 + µ44 )) + µ41 µ23 + µ41 µ24 + µ21 µ43 + µ21 µ44 + µ43 µ24 + µ23 µ44 − 4µ21 µ23 µ24   − ν4 2µ42 (µ21 µ23 + µ21 µ24 + µ23 µ24 ) − µ22 (µ41 µ23 + µ21 µ43 + µ41 µ24 + µ21 µ44 + µ43 µ24 + µ23 µ44 + 6µ21 µ23 µ24 )
 ρ ( µ1 ) ρ ( µ3 ) ρ ( µ4 ) −ν4 2µ21 µ23 µ24 (µ21 + µ23 + µ24 ) (µ1 − µ22 )(µ23 − µ22 )(µ24 − µ22 ) 2 ψ4  4  4 2 2 2 2 2 2 3 2 2 2  2 2  − 3ξ 3µ 1 − 2µ ( 1 2µ + µ 3 + µ 4 ) + µ µ 2 3 + µ µ 2 4 + µ µ 3 4 − ξ ν 7µ61 − 4µ41 (µ22 + µ23 + µ24 ) + µ21 × ν4  ×(3(µ22 µ23 + µ22 µ24 + µ23 µ24 ) − 2(µ42 + µ43 + µ44 )) + µ42 µ23 + µ42 µ24 + µ22 µ43 + µ22 µ44 + µ43 µ24 + µ23 µ44 − 4µ22 µ23 µ24   − ν4 2µ41 (µ22 µ23 + µ22 µ24 + µ23 µ24 ) − µ21 (µ42 µ23 + µ22 µ43 + µ42 µ24 + µ22 µ44 + µ43 µ24 + µ23 µ44 + 6µ22 µ23 µ24 )
 ρ ( µ2 ) ρ ( µ3 ) ρ ( µ4 ) −ν4 2µ22 µ23 µ24 (µ22 + µ23 + µ24 ) (µ2 − µ21 )(µ23 − µ21 )(µ24 − µ21 ) 2 ψ4 2

− 4 ξ 3ξ 2 − µ21 + µ22 + µ23 + µ24 ν2 ρ(µ1 )ρ(µ2 )ρ(µ3 )ρ(µ4 ) (B.5) ν -27- REFERENCES REFERENCES References [1] M. Gaudin, Diagonalisation d’une classe d’hamiltoniens de spin, J. Physique 37 (1976) 1087–1098. [2] M. Gaudin, La fonction d’onde de Bethe, chapter 13 Masson, Paris, 1983. [3] E. K. Sklyanin, Separation of variables in the Gaudin model, Zap. Nauchn. Sem. Leningrad. Otdel. Mat. Inst. Steklov. (LOMI) 164 (1987) 151–169; translation in J. Soviet Math. 47 (1989) 2473–2488. [4] A. A. Belavin and V. G. Drinfeld. Solutions of the classical Yang-Baxter equation for simple Lie algebras (in Russian), Funktsional. Anal. i Prilozhen. 16 (1982), no. 3, 1–29; translation in Funct. Anal. Appl. 16 (1982) no. 3, 159-180. [5] E. K. Sklyanin and T. Takebe, Algebraic Bethe ansatz for the XYZ Gaudin model, Phys. Lett. A 219 (1996) 217-225. [6] M. A. Semenov-Tian-Shansky, Quantum and classical integrable systems, in Integrability of Nonlinear Systems, Lecture Notes in Physics Volume 495 (1997) 314-377. [7] B. Jurčo, Classical Yang-Baxter equations and quantum integrable systems, J. Math. Phys. Volume 30 (1989) 1289–1293. [8] B. Jurčo, Classical Yang-Baxter equations and quantum integrable systems (Gaudin mod- els), in Quantum groups (Clausthal, 1989), Lecture Notes in Phys. Volume 370 (1990) 219– 227. [9] H. M. Babujian and R. Flume, Off-shell Bethe ansatz equation for Gaudin magnets and solutions of Knizhnik-Zamolodchikov equations, Mod. Phys. Lett. A 9 (1994) 2029–2039. [10] B. Feigin, E. Frenkel, and N. Reshetikhin, Gaudin model, Bethe ansatz and correlation func- tions at the critical level, Commun. Math. Phys. 166 (1994) 27–62. [11] N. Reshetikhin and A. Varchenko, Quasiclassical asymptotics of solutions to the KZ equa- tions, in Geometry, Topology and Physics, Conf. Proc. Lecture Notes Geom., pages 293– 273. Internat. Press, Cambridge, MA, 1995. [12] F. Wagner and A. J. Macfarlane, Solvable Gaudin models for higher rank symplectic algebras. Quantum groups and integrable systems (Prague, 2000) Czechoslovak J. Phys. 50 (2000) 1371–1377. [13] T. Brzezinski and A. J. Macfarlane, On integrable models related to the osp(1,2) Gaudin algebra, J. Math. Phys. 35 (1994), no. 7, 3261–3272. [14] P. P. Kulish and N. Manojlović, Creation operators and Bethe vectors of the osp(1|2) Gaudin model, J. Math. Phys. 42 no. 10 (2001) 4757–4778. -28- REFERENCES REFERENCES [15] P. P. Kulish and N. Manojlović, Trigonometric osp(1|2) Gaudin model, J. Math.Phys. 44 no. 2 (2003) 676–700. [16] A. Lima-Santos and W. Utiel, Off-shell Bethe ansatz equation for osp(2|1) Gaudin mag- nets, Nucl. Phys. B 600 (2001) 512–530. [17] V. Kurak and A. Lima-Santos, sl (2|1)(2) Gaudin magnet and its associated Knizhnik- Zamolodchikov equation, Nuclear Physics B 701 (2004) 497–515. [18] K. Hikami, P. P. Kulish and M. Wadati, Integrable Spin Systems with Long-Range Interaction, Chaos, Solitons & Fractals Vol. 2 No. 5 (1992) 543–550. [19] K. Hikami, P. P. Kulish and M. Wadati, Construction of Integrable Spin Systems with Long-Range Interaction, J. Phys. Soc. Japan Vol. 61 No. 9 (1992) 3071–3076. [20] K. Hikami, Gaudin magnet with boundary and generalized Knizhnik-Zamolodchikov equation, J. Phys. A Math. Gen. 28 (1995) 4997–5007. [21] W. L. Yang, R. Sasaki and Y. Z. Zhang, Z n elliptic Gaudin model with open boundaries, JHEP 09 (2004) 046. [22] W. L. Yang, R. Sasaki and Y. Z. Zhang, An−1 Gaudin model with open boundaries, Nuclear Physics B 729 (2005) 594–610. [23] K. Hao, W.-L. Yang, H. Fan, S. Y. Liu, K. Wu, Z. Y. Yang and Y. Z. Zhang, Determi- nant representations for scalar products of the XXZ Gaudin model with general boundary terms, Nuclear Physics B 862 (2012) 835–849. [24] A. Lima-Santos, The sl (2|1)(2) Gaudin magnet with diagonal boundary terms, J. Stat. Mech. (2009) P07025. [25] E. K. Sklyanin, Boundary conditions for integrable equations, (Russian) Funktsional. Anal. i Prilozhen. 21 (1987) 86–87; translation in Functional Analysis and Its Applications Volume 21, Issue 2 (1987) 164–166. [26] E. K. Sklyanin, Boundary conditions for integrable systems, in the Proceedings of the VIIIth international congress on mathematical physics (Marseille, 1986), World Sci. Pub- lishing, Singapore, (1987) 402–408. [27] E. K. Sklyanin, Boundary conditions for integrable quantum systems, J. Phys. A: Math. Gen. 21 (1988) 2375–2389. [28] T. Skrypnyk, Non-skew-symmetric classical r-matrix, algebraic Bethe ansatz, and Bardeen-Cooper-Schrieffer-type integrable systems, J. Math. Phys. 50 (2009) 033540, 28 pages. -29- REFERENCES REFERENCES [29] T. Skrypnyk, "Z2-graded” Gaudin models and analytical Bethe ansatz, Nuclear Physics B 870 (2013), no. 3, 495–529. [30] N. Cirilo António, N. Manojlović and Z. Nagy, Trigonometric sℓ(2) Gaudin model with boundary terms, Reviews in Mathematical Physics Vol. 25 No. 10 (2013) 1343004 (14 pages). [31] J. Cao, H. Lin, K. Shi and Y. Wang, Exact solutions and elementary excitations in the XXZ spin chain with unparallel boundary fields, Nucl. Phys. B 663 (2003) 487–519. [32] R. I. Nepomechie, Bethe ansatz solution of the open XXZ chain with nondiagonal bound- ary terms, J. Phys. A 37 (2004), no. 2, 433–440. [33] D. Arnaudon, A. Doikou, L. Frappat, E. Ragoucy and N. Crampé, Analytical Bethe ansatz in gl ( N ) spin chains, Czechoslovak J. Phys. 56 (2006), no. 2, 141–148. [34] C. S. Melo, G. A. P. Ribeiro and M. J. Martins, Bethe ansatz for the XXX − S chain with non-diagonal open boundaries, Nuclear Phys. B 711, no. 3 (2005) 565–603. [35] J. Cao, W.-L. Yang, K. Shi and Y. Wang, Off-diagonal Bethe ansatz solution of the XXX spin chain with arbitrary boundary conditions, Nuclear Physics B 875 (2013) 152–165. [36] J. Cao, W.-L. Yang, K. Shi and Y. Wang, Off-diagonal Bethe ansatz solutions of the anisotropic spin-1/2 chains with arbitrary boundary fields, Nuclear Physics B 877 (2013) 152–175. [37] E. Ragoucy, Coordinate Bethe ansätze for non-diagonal boundaries, Rev. Math. Phys. 25 (2013), no. 10, 1343007. [38] S. Belliard, N. Crampé and E. Ragoucy, Algebraic Bethe ansatz for open XXX model with triangular boundary matrices, Lett. Math. Phys. 103 No. 5 (2013) 493–506. [39] S. Belliard and N. Crampé Heisenberg XXX model with general boundaries: eigenvec- tors from algebraic Bethe ansatz, SIGMA Symmetry Integrability Geom. Methods Appl. 9 (2013), Paper 072, 12 pp. [40] R. A. Pimenta and A. Lima-Santos, Algebraic Bethe ansatz for the six vertex model with upper triangular K-matrices, J. Phys. A 46 No. 45 (2013) 455002, 13 pp. [41] N. Cirilo António, N. Manojlović and I. Salom, Algebraic Bethe ansatz for the XXX chain with triangular boundaries and Gaudin model, Nuclear Physics B 889 (2014) 87-108. [42] K. Hao, J. Cao, T. Yang and W.-L. Yang, Exact solution of the XXX Gaudin model with the generic open boundaries, arXiv:1408.3012. [43] C. N. Yang, Some exact results for the many-body problem in one dimension with repul- sive delta-function interaction, Phys. Rev. Lett. 19 (1967) 1312-1315. -30- REFERENCES REFERENCES [44] R. J. Baxter, Exactly solved models in statistical mechanics, Academic Press, London (1982). [45] L. D. Faddeev, How the algebraic Bethe Ansatz works for integrable models, In Quan- tum symmetries / Symetries quantiques, Proceedings of the Les Houches summer school, Session LXIV. Eds. A. Connes, K. Gawedzki and J. Zinn-Justin. North-Holland, 1998, 149– 219; hep-th/9605187. [46] P. P. Kulish Twist deformations of quantum integrable spin chains, Lect. Notes Phys. Vol. 774 (2009) 165–188. [47] H. J. de Vega and A. González Ruiz, Boundary K -matrices for the XYZ, XXZ, XXX spin chains, J. Phys. A: Math. Gen. 27 (1994), 6129–6137. [48] J.-M. Maillet, New integrable canonical structures in two-dimensional models, Nuclear Physics B269 (1986) 54–76. -31-