  In the derivation of the generating function of the Gaudin Hamiltonians with
boundary terms, we follow the same approach used previously in the rational
case, which in turn was based on Sklyanin's method in the periodic case. Our
derivation is centered on the quasi-classical expansion of the linear
combination of the transfer matrix of the XXZ 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.
