%% 
%% Copyright 2007-2020 Elsevier Ltd
%% 
%% This file is part of the 'Elsarticle Bundle'.
%% ---------------------------------------------
%% 
%% It may be distributed under the conditions of the LaTeX Project Public
%% License, either version 1.2 of this license or (at your option) any
%% later version.  The latest version of this license is in
%%    http://www.latex-project.org/lppl.txt
%% and version 1.2 or later is part of all distributions of LaTeX
%% version 1999/12/01 or later.
%% 
%% The list of all files belonging to the 'Elsarticle Bundle' is
%% given in the file `manifest.txt'.
%% 

%% Template article for Elsevier's document class `elsarticle'
%% with numbered style bibliographic references
%% SP 2008/03/01
%%
%% 
%%
%% $Id: elsarticle-template-num.tex 190 2020-11-23 11:12:32Z rishi $
%%
%%
\documentclass[preprint,12pt]{elsarticle}

%% Use the option review to obtain double line spacing
%% \documentclass[authoryear,preprint,review,12pt]{elsarticle}

%% Use the options 1p,twocolumn; 3p; 3p,twocolumn; 5p; or 5p,twocolumn
%% for a journal layout:
%% \documentclass[final,1p,times]{elsarticle}
%% \documentclass[final,1p,times,twocolumn]{elsarticle}
%% \documentclass[final,3p,times]{elsarticle}
%% \documentclass[final,3p,times,twocolumn]{elsarticle}
%% \documentclass[final,5p,times]{elsarticle}
%% \documentclass[final,5p,times,twocolumn]{elsarticle}

%% For including figures, graphicx.sty has been loaded in
%% elsarticle.cls. If you prefer to use the old commands
%% please give \usepackage{epsfig}

%% The amssymb package provides various useful mathematical symbols
%\usepackage{amssymb}
%% The amsthm package provides extended theorem environments

%\usepackage{amsthm}


\usepackage{amsfonts,latexsym,amsthm,amssymb}
\usepackage{amsmath,bbm}  % Better maths support & more symbols
\usepackage{mathtools}
\usepackage[mathscr]{eucal}


% \usepackage{mathpazo} % math & rm
% \usepackage[scaled]{helvet} % ss
% \usepackage{courier} % tt
% \usepackage[T1]{fontenc}

%\usepackage[english]{babel}
%\usepackage[english]{layout}
%\usepackage{cite}

\newcommand{\so}{{\scriptscriptstyle \rm I}}
\newcommand{\st}{{\scriptscriptstyle \rm I\hspace{-.1ex}I}}
\newcommand{\sth}{{\scriptscriptstyle \rm I\hspace{-.1ex}I\hspace{-.1ex}I}}
\newcommand{\Res}{\operatorname{Res}}
\newcommand{\tr}{\operatorname{tr}}


%% The lineno packages adds line numbers. Start line numbering with
%% \begin{linenumbers}, end it with \end{linenumbers}. Or switch it on
%% for the whole article with \linenumbers.
%% \usepackage{lineno}

\journal{Nuclear Physics B}

\begin{document}

\begin{frontmatter}

%% Title, authors and addresses

%% use the tnoteref command within \title for footnotes;
%% use the tnotetext command for theassociated footnote;
%% use the fnref command within \author or \address for footnotes;
%% use the fntext command for theassociated footnote;
%% use the corref command within \author for corresponding author footnotes;
%% use the cortext command for theassociated footnote;
%% use the ead command for the email address,
%% and the form \ead[url] for the home page:
%% \title{Title\tnoteref{label1}}
%% \tnotetext[label1]{}
%% \author{Name\corref{cor1}\fnref{label2}}
%% \ead{email address}
%% \ead[url]{home page}
%% \fntext[label2]{}
%% \cortext[cor1]{}
%% \affiliation{organization={},
%%             addressline={},
%%             city={},
%%             postcode={},
%%             state={},
%%             country={}}
%% \fntext[label3]{}

\title{Rational $so(3)$ Gaudin model with general boundary terms}

%% use optional labels to link authors explicitly to addresses:
%% \author[label1,label2]{}
%% \affiliation[label1]{organization={},
%%             addressline={},
%%             city={},
%%             postcode={},
%%             state={},
%%             country={}}
%%
%% \affiliation[label2]{organization={},
%%             addressline={},
%%             city={},
%%             postcode={},
%%             state={},
%%             country={}}

\author[1]{N.~Manojlovi\'c}

\affiliation[1]{organization={Departamento de Matem\'atica, Faculdade de Ci\^encias e Tecnologia},%Department and Organization
            addressline={Universidade do Algarve, Campus de Gambelas}, 
            city={Faro},
            postcode={PT-8005-139}, 
%            state={},
            country={Portugal}}

\author[2]{I.~Salom}

\affiliation[2]{organization={Institute of Physics},%Department and Organization
            addressline={University of Belgrade}, 
            city={Belgrade},
            postcode={P.O. Box 57, 11080}, 
%            state={},
            country={Serbia}}
            
            
\begin{abstract}
%% Text of abstract
We study the $so(3)$ Gaudin model with general boundary K-matrix in the framework of the algebraic Bethe ansatz.  The off-shell action of the generating function of the $so(3)$ Gaudin Hamiltonians is determined. The proof based on the mathematical induction is presented on the algebraic level without any restriction whatsoever on the boundary parameters. The $so(3)$ Gaudin Hamiltonians with general boundary terms are given explicitly as well as their off-shell action on the Bethe states. The correspondence between the Bethe states and the solutions to the generalized $so(3)$ Knizhnik-Zamolodchikov equations is established. In this context, the on-shell norm of the Bethe states is determined as well as their off-shell scalar product.
\end{abstract}

%%Graphical abstract
% \begin{graphicalabstract}
% \includegraphics{grabs}
% \end{graphicalabstract}
% 
% Research highlights
% \begin{highlights}
% \item Research highlight 1
% \item Research highlight 2
% \end{highlights}

%\begin{keyword}
%% keywords here, in the form: keyword \sep keyword

%% PACS codes here, in the form: \PACS code \sep code

%% MSC codes here, in the form: \MSC code \sep code
%% or \MSC[2008] code \sep code (2000 is the default)

%\end{keyword}

\end{frontmatter}

%% \linenumbers

%% main text

\section{Introduction} 
%
%
The systems obtained as the quasi-classical limit of the Heisenberg spin chains \cite{Heisenberg28} were first studied by Gaudin   \cite{Gaudin76, Gaudin83, Gaudin-English}. In the framework of the coordinate as well as the algebraic Bethe ansatz Gaudin has found the spectrum of the generating function of the corresponding Hamiltonians \cite{Gaudin76, Gaudin83, Gaudin-English}. This system has been recasted in the framework of the quantum inverse scattering method \cite{TakhtajanFaddeev79, KulishSklyanin82, Faddeev98} by exploring the so-called Sklyanin linear bracket  using an $s\ell(2)$ invariant, unitary classical r-matrix \cite{Sklyanin:1987ih}. This result enabled further generalisations based on other unitary solutions to the classical Yang-Baxter equation \cite{KulSkl80, BelavinDrinfeld}, prompting the interest in the Gaudin systems based on higher-rank simple Lie algebras \cite{Jurco90, Semenov97, WagnerMacfarlane00} as well as  Lie superalgebras \cite{BrzezinskiMacfarlane94, Kulish:2001hr, Kulish:2001vs, Kulish:2003da, LimaUtiel01}. The relation with the Knizhnik-Zamolodchikov equations of conformal field theory in two dimension \cite{KZ84} and the representation theory of the Kac-Moody algebras \cite{ReshetikhinVarchenko} was further strengthened when the connection between the Bethe states of the Gaudin model and the solution to the Knizhnik-Zamolodchikov equations was established \cite{FeginFrenkelReshetikhin, BabujianFlume, Kulish:2001hr, Kulish:2001vs, Kulish:2003da, KurakLima04}. It is also interesting to note that, in a somewhat more physical approach, the long-range interaction of these  systems was studied in \cite{HikamiKulishWadati92, HikamiKulishWadati92a}. The Gaudin system on an elliptic curve was studied in \cite{SklyaninTakebe96}, while the  $s\ell(2)$ Gaudin with the Jordanian twist was studied in \cite{CiriloAntonio:2005qts, KMSS, Cirilo-Antonio:2010xmg}. On the classical level, the Gaudin model corresponds to the so-called Schlesinger system in the theory of isomonodromic deformation \cite{Schlesinger,Jimbo:1981aa,Jimbo:1981bb,Jimbo:1981cc,Korotkin:1999xx,Manojlovic:2001mv,Dubrovin2007}. 

In our considerations of the quantum Heisenberg spin chains with non-periodic boundary conditions we follow Sklyanin's approach where the boundary conditions are expressed in the form of the left and right reflection matrices \cite{Sklyanin88}. The so-called reflection equation and the dual reflection equation represent the compatibility conditions between the bulk and the boundary of the system at the left and, respectively,  right site of the system. The commutativity of the transfer matrix, in this case, is guaranteed on the one hand by the fact that the matrix form of the exchange relations between the entries of the Sklyanin monodromy matrix is analogous to the reflection equation and on the other hand, by the dual reflection equation \cite{Sklyanin88, FreidelMaillet91, FreidelMaillet91a}.

Renewed interest has emerged in the implementation of algebraic Bethe ansatz on solvable Heisenberg chains with non-periodic boundary conditions \cite{Martins05, Eric13, Belliard13, Lima13, Belliard15, BelliardPimenta15, AvanEtal15, Nepomechie15,Belliard:2021esn}. As for alternative approaches, a review of the coordinate Bethe ansatz in this case is given in \cite{Eric2013RMP}, the Bethe ansatz based on the functional relation between the eigenvalues of the transfer matrix and the quantum determinant, as well as the associated T-Q relation are studied in \cite{LucEricRafael07,China13, China15}, functional relations for the eigenvalues of the transfer matrix based on fusion hierarchy were discussed in \cite{Nepomechie04} and the Vertex-IRF correspondence in \cite{Jimbo95, China03}, while the Jordanian deformation of the open XXX chain is analysed in \cite{Kulish:2009ci}. For the latest results, as well as an excellent review on the application of the separation of variables method on the 6-vertex model and the associate XXZ quantum chains see \cite{Maillet17}.

Our interest in open Heisenberg spin chains was twofold. On the one hand, we were interested in the implementation of the algebraic Bethe ansatz, and on the other hand, we wanted to consider the quasi-classical limit which yields the corresponding Gaudin model \cite{CiriloAntonio:2014ymf, Manojlovic:2017nbh}. As it is well known \cite{Eric13, Belliard13, CiriloAntonio:2014ymf}, due to the symmetry of the R-matrix of the non-periodic XXX spin chain, the accomplishment of the algebraic Bethe ansatz does not imply any restriction on the boundary parameters. However, in the case of the open XXZ chain \cite{CiriloAntonio:2013eng}, the existence of the so-called vacuum vector requires the triangular form of the boundary K-matrix \cite{Lima13, Belliard15, Manojlovic:2017nbh}. As for the quasi-classical limit, Hikami showed, in complete analogy with the periodic case \cite{HikamiKulishWadati92, HikamiKulishWadati92a}, how the expansion of the XXZ transfer matrix, calculated at the special values of the spectral parameter, yields the Gaudin Hamiltonians in the case when both reflection matrices are diagonal \cite{Hikami95}. Similar expansion was done for the  Jordanian deformation of the ration $s\ell(2)$ Gaudin model with generic boundaries \cite{CiriloAntonio:2013rbz}. The algebraic Bethe ansatz was applied to open Gaudin model in the context of the Vertex-IRF correspondence \cite{YangZhangSasakic04, YangZhangSasakic05, YangZhang12}. Also,  results were obtained for the open Gaudin models based on Lie superalgebras \cite{Lima09}. Returning back to the quasi-classical limit, following the Sklyanin proposal for the periodic boundary conditions  \cite{Sklyanin:1987ih, Manojlovic:2015jhi}, we have derived the generating function of the Gaudin Hamiltonians both for the XXX \cite{CiriloAntonio:2014ymf} and the XXZ chain \cite{Manojlovic:2017nbh} as well as for the Jordanian deformation of the XXX Heisenberg spin chain \cite{Manojlovic:2018abc}. Moreover, we have shown \cite{CiriloAntonio:2014bfg} how,  in the context of the quasi-classical limit, the solutions to the classical Yang-Baxter equation \cite{KulSkl80, BelavinDrinfeld} can be combined with the solutions to the classical reflection equation \cite{Sklyanin86, Sklyanin87} to yield solutions to the so-called generalized classical Yang-Baxter equation \cite{Maillet85, Maillet86, BabelonViallet90, AvanTalon90}. These solutions are the non-unitary classical r-matrices \cite{Skrypnyk06, Skrypnyk07, Skrypnyk09, Skrypnyk10, Skrypnyk10a, Skrypnyk13, Skrypnyk16}. In particular, the generic elliptic $s\ell(2)$ non-unitary r-matrix was studied in \cite{Manojlovic:2020xou}. Also, we draw attention to the recent study of the generalized Gaudin and Richardson models based on a class of non-unitary r-matrices \cite{Skrypnyk:2021xxc}.

An approach to the implementation of the algebraic Bethe ansatz for the rational as well as the trigonometric $s\ell(2)$ Gaudin model based on the corresponding Maillet linear bracket was developed in \cite{Salom:2015wok, Salom:2019myz, Manojlovic:2017yxv, Crampe17, Salom:2021gox}. Once a suitable set of generators of the relevant generalized Gaudin algebra is found, the local realization of these generators becomes compact and it naturally leads to the definition of the so-called creation operators. In both the rational and trigonometric case \cite{Salom:2019myz, Salom:2021gox}, these creation operators define Bethe states in such a way that the off-shell action of the generating function of the Gaudin Hamiltonians can be computed explicitly, and  a completely algebraic proof of this action given. 

This paper is centred on the application of the algebraic Bethe ansatz to the rational $so(3)$ Gaudin model with generic classical boundary K-matrix. We recall that this K-matrix can be obtained by the so-called fusion procedure \cite{KulishSklyanin82, Mezincescu:1991ke, Kulish:2009cx}, starting from the $s\ell(2)$ K-matrix \cite{CiriloAntonio:2014ymf, VegaGonzalez, Zamolodchikov94, Zamolodchikov94b}. The outline of this methods in the trigonometric $so(3)$ case was given in \cite{Salom:2020dxt}. Alternatively, one can use the so-called scaling limit \cite{Kulish:2009aa}, to obtain the K-matrix 
from the trigonometric $so(3)$ boundary K-matrix \cite{Salom:2020dxt, Inami:1996np}. The non-unitary $so(3)$ classical r-matrix \eqref{Eq: rK} is then obtained by combining the unitary $so(3)$ invariant classical r-matrix and the K-matrix. This non-unitary classical r-matrix defines the $so(3)$ Maillet linear bracket, for the suitable Lax operator, and provides an algebraic framework for our study of the non-periodic $so(3)$ Gaudin model. As an immediate consequence of the definition of the $so(3)$ Maillet bracket follows the mutual commutativity of the generating function for different values of the spectral parameter. However, as it will be confirmed in the following, the natural set of generators unfortunately turns out not to be adequate for the  implementation of the algebraic Bethe ansatz. Thus we will here propose a new set of generators. Besides the relative simplicity of the local realization of the new generators, their most striking feature will be the compact form of their commutation relations. This is of great significance since it efficiently enables the algebraic proof of the off-shell action of the generating function on the Bethe states. Furthermore, it is important to stress that these results will be obtained without any restriction whatsoever on the boundary parameters. It is only when solving the generalized $so(3)$ Knizhnik-Zamolodchikov equations that the key identity in the proof will require one of four boundary parameters to be set to zero. However, in spite of this constraint we will retain a large improvement in generality over the previous studies: while the formulas that we here provide for the solutions to the generalized $so(3)$ Knizhnik-Zamolodchikov equations, the on-shell norm of the Bethe vectors and the off-shell scalar product of the Bethe vectors do superficially look similar to the analogous formulae in the $s\ell(2)$ case \cite{Salom:2019myz}, only one of the boundary parameters will be fixed here, instead of all four of them (as, for example, in \cite{Salom:2019myz}).

The paper is organised as follows. In Section \ref{Sec: Mbr} we study the $so(3)$ Maillet linear bracket  which provides the algebraic framework for implementation of the Bethe ansatz. In the same section we propose the novel set of generators with simplified commutation relations and introduce Gaudin Hamiltonians. The implementation of the algebraic Bethe ansatz is the principal topic of the Section \ref{Sec: ABA}. There we will obtain the expression for the off-shell action of the generating function $\tau (\lambda)$,  as well as for the off-shell action of the $so(3)$ Gaudin Hamiltonians with general boundary terms  -- and prove these formulas by mathematical induction. The solutions to the generalized  $so(3)$ Knizhnik-Zamolodchikov equations will be given in the Section \ref{Sec: gKZeq}.  Our results will be summarised in the concluding Section \ref{Sec: Con}. Fundamental definitions regarding the $so(3)$ Lie algebra, including the two $so(3)$ invariant operators in $\mathbb{C}^3\otimes \mathbb{C}^3$ which generate the relevant Brauer algebra, are presented in the Appendix \ref{App: Prelim1}. Finally, the cornerstone of our study -- the non-unitary $so(3)$ classical r-matrix 
-- is given in the Appendix \ref{App: Prelim2}.



\section{The $so(3)$ Maillet linear bracket \label{Sec: Mbr} }
%
%
In this section we show how the non-unitary $so(3)$ classical r-matrix \eqref{Eq: rK} helps define the $so(3)$ Maillet linear bracket \eqref{Eq: M-br-rk} for the suitable Lax operator \eqref{Eq: cal-L}. Although this Maillet bracket provides an appropriate algebraic framework for studying the quantum $so(3)$ Gaudin model, yielding the generating function of the $so(3)$ Gaudin Hamiltonians with general boundary terms, it will be shown below that the natural set of generators unfortunately does not provide the most efficient way for implementing the algebraic Bethe ansatz in this case. Thus we will propos a new set off generators of the corresponding generalized $so(3)$ Gaudin algebra.



In our study we use the Lax operator
\begin{equation}
\label{Eq: L-Gdin}
L_0 (\lambda) = \sum _{m=1}^N \frac{\vec{\mathscr{S}}_0 \cdot \vec{S}_m}{\lambda - \alpha _m} 
= \sum _{m=1}^N \frac{1}{\lambda - \alpha _m} \left( \mathscr{S}_{0}^3 \otimes S_{m}^3 + \frac{1}{2} \left(  \mathscr{S}_{0}^+ \otimes S_{m}^- +  \mathscr{S}_{0}^- \otimes S_{m}^+ \right) \right)  \>, 
\end{equation}
where the spin operators $S_m^{\alpha}$, with $\alpha = +, - , 3$ and $m = 1, 2 ,\dots , N$, are introduced in \eqref{Eq: S-local} and the matrices $\mathscr{S} ^3_0$ and $\mathscr{S} ^{\pm}_0$ in the auxiliary space $\mathbb{C}^{3}$ are specified by \eqref{Eq: S-mats} and \eqref{Eq: S-pm-mats}, respectively.
The Lax operator \eqref{Eq: L-Gdin} can also be represented in the following form
\begin{equation}
\label{Eq: L-Gdin-S}
L_0 (\lambda) = \vec{\mathscr{S}}_0 \cdot \vec{S} (\lambda) \>,
\end{equation}
where the generators of the $so(3)$ Gaudin algebra are defined by \cite{Gaudin76, Gaudin83, Gaudin-English}
\begin{equation}
\label{Eq: S-gens}
S^3(\lambda) = \sum _{m=1}^N \frac{S^3_m}{\lambda - \alpha _m} \>, \qquad S^{\pm}(\lambda) = \sum _{m=1}^N \frac{S^{\pm}_m}{\lambda - \alpha _m} \>.
\end{equation}
The so-called Sklyanin linear bracket \cite{Sklyanin:1987ih, Kulish:2001hr, Kulish:2001vs} for the Lax operator \eqref{Eq: L-Gdin} and the r-matrix \eqref{Eq: r-so3} 
\begin{equation}
\label{Eq: Skly-br}
\left[ L _1(\lambda) , L _2 (\mu) \right] = \left[ r_{12}(\lambda - \mu) , L _1(\lambda) + L _2 (\mu) \right] 
\end{equation}
yields nontrivial commutation relations for the generators \eqref{Eq: S-gens}
\begin{align}
\left[ S^3(\lambda) , S^{\pm}(\mu) \right] &= \mp \, \frac{S^{\pm}(\lambda) - S^{\pm}(\mu) }{\lambda - \mu} \>, \nonumber \\ 
%
\left[ S^{+}(\lambda) , S^{-}(\mu) \right] &= ( - 2 ) \, \frac{S^{3}(\lambda) - S^{3}(\mu) }{\lambda - \mu} \>. 
\label{Eq: so3-S-rels}
\end{align}




The Lax operator corresponding to the generalized  $so(3)$ Gaudin algebra is given by
\begin{equation}
\label{Eq: cal-L-def}
\mathcal{L} _0 (\lambda) = L _0 (\lambda) - K_0 (\lambda)  L _0 (- \lambda) K_0^{-1}(\lambda) \>,
\end{equation}
where $L _0 (\lambda)$ is the Lax operator  \eqref{Eq: L-Gdin} and  $K_0 (\lambda)$ is the reflection K-matrix defined in \eqref{Eq: K-mat}. This form of the Lax operator can be obtained by following a relatively general procedure of quasi-classical expansion of the Sklyanin monodromy \cite{CiriloAntonio:2014bfg}. By direct substitution we obtain
\begin{equation}
\label{Eq: cal-L}
\begin{aligned}
\mathcal{L} _0 (\lambda) 
&= \vec{\mathscr{S}}_0 \cdot \vec{S} (\lambda) - \left( K_0 (\lambda) \vec{\mathscr{S}}_0 K_0^{-1} (\lambda) \right) \cdot \vec{S} (-\lambda) 
\\[1ex]
%
&= \left(\begin{array}{ccc}
H (\lambda) & \frac{1}{\sqrt{2}} F (\lambda) & 0 \\[1ex]
\frac{1}{\sqrt{2}} E (\lambda) & 0 & \frac{1}{\sqrt{2}} F (\lambda) \\[1ex]
0 & \frac{1}{\sqrt{2}} E (\lambda) & - H (\lambda)
\end{array}\right) \>.
\end{aligned}
\end{equation}
The Lax operator \eqref{Eq: cal-L} obeys the following $so(3)$ Maillet linear bracket \cite{Maillet86, BabelonViallet90, AvanTalon90, FreidelMaillet91,  CiriloAntonio:2014bfg}
\begin{equation}
\label{Eq: M-br-rk}
\left[ \mathcal{L} _{0} (\lambda) , \mathcal{L} _{0'}(\mu) \right] = \left[ r_{00'}^{K} (\lambda , \mu) , \mathcal{L} _{0} (\lambda) \right] - \left[ r_{0'0}^{K} (\mu , \lambda) , \mathcal{L} _{0'} (\mu) \right] .
\end{equation}
This linear bracket is obviously anti-symmetric and it obeys the Jacobi identity because the $r$-matrix \eqref{Eq: rK-def} satisfies the generalized classical Yang-Baxter equation \eqref{Eq: g-c-YBE}.


The $so(3)$ Maillet bracket \eqref{Eq: M-br-rk} implies the following commutation relations for the generators $E (\lambda)$, $F (\lambda) $ and $H (\lambda)$ \eqref{Eq: cal-L}
\begin{align}
\label{Eq: cr-EE}
\left[ E (\lambda) , E (\mu)  \right] &= \frac{-2 \varphi ^2}{\lambda + \mu}  \left( \frac{\mu^2}{\xi ^2 - (\psi \varphi + \nu ^2) \mu ^2} \, H (\lambda) - \frac{\lambda^2}{\xi ^2 - (\psi \varphi + \nu ^2) \lambda ^2} \, H (\mu) \right)  \>,
\notag \\[1ex]
%
&+ \frac{ 2 \varphi}{\lambda + \mu}  \left( \frac{(\xi + \nu \mu) \, \mu}{\xi ^2 - (\psi \varphi + \nu ^2) \mu ^2} \, E (\lambda) - \frac{(\xi + \nu \lambda) \, \lambda}{\xi ^2 - (\psi \varphi + \nu ^2) \lambda} \, E (\mu) \right)   \>,
\\[2ex]
%
%
\label{Eq: cr-FF}
\left[ F (\lambda) , F (\mu)  \right] &=  \frac{2 \psi ^2}{\lambda + \mu} \left( \frac{\mu^2}{\xi ^2 - (\psi \varphi + \nu ^2) \mu ^2} \, H (\lambda) - \frac{\lambda^2}{\xi ^2 - (\psi \varphi + \nu ^2) \lambda ^2} \, H (\mu) \right)
\notag \\[1ex]
%
&+ \frac{2 \psi}{\lambda + \mu} \left( \frac{(\xi - \nu \mu) \, \mu}{\xi ^2 - (\psi \varphi + \nu ^2) \mu ^2} \, F (\lambda) - \frac{(\xi - \nu \lambda) \, \lambda}{\xi ^2 - (\psi \varphi + \nu ^2) \lambda} \, F (\mu) \right)   \>,
\\[2ex]
% 
%
\label{Eq: cr-HH}
\left[ H (\lambda) , H (\mu)  \right] &= \frac{ - \psi}{\lambda + \mu} \left(  \frac{(\xi + \nu \mu) \, \mu}{\xi ^2 - (\psi \varphi + \nu ^2) \, \mu^2} \, E (\lambda) - \frac{(\xi + \nu \lambda) \, \lambda}{\xi ^2 - (\psi \varphi + \, \nu ^2) \lambda ^2} \, E (\mu) \right)   \notag \\[1ex]
% 
&+ \frac{ - \varphi}{\lambda + \mu}  \left( \frac{(\xi - \nu \mu) \, \mu}{\xi ^2 - (\psi \varphi + \nu ^2) \mu^2} \, F (\lambda) - \frac{(\xi - \nu \lambda) \, \lambda}{\xi ^2 - (\psi \varphi + \nu ^2) \lambda ^2} \, F (\mu) \right) \>,
\end{align} 
and 
\begin{align}
\label{Eq: cr-HE}
\left[ H (\lambda) , E (\mu)  \right] &=   \frac{\varphi}{\lambda + \mu} \left( 
\frac{\varphi \, \mu^2 }{\xi ^2 - (\psi \varphi + \nu ^2) \, \mu^2} \, F (\lambda) - \frac{2 \, (\xi - \nu \lambda) \, \lambda}{\xi ^2 - (\psi \varphi + \, \nu ^2) \lambda ^2} \, H (\mu) \right) 
\notag \\[1ex]
%
& - \frac{1}{(\lambda - \mu)(\lambda + \mu)} 
\bigg( \frac{\left(2 (\xi - \nu \lambda ) \, (\xi + \nu \mu )  - \psi \varphi \, (\lambda + \mu) \mu \right)  \mu}{\xi ^2 - (\psi \varphi + \nu ^2) \mu^2} \, E (\lambda)  
\notag \\[1ex]
%
& - \frac{2 \left( \xi ^2 - \left( \psi \varphi \, \mu + \nu ^2 \lambda \right) \lambda \right) \lambda}{\xi ^2 - (\psi \varphi + \, \nu ^2) \lambda ^2} \, E (\mu) \bigg)  \>,
\\[2.5ex]
%
%
\label{Eq: cr-HF}    
\left[ H (\lambda) , F (\mu)  \right] &=   \frac{ - \psi }{\lambda + \mu} \left(  \frac{\psi \, \mu^2 }{\xi ^2 - (\psi \varphi + \nu ^2) \, \mu^2} \, E (\lambda) + \frac{2 \, (\xi + \nu \lambda) \, \lambda}{\xi ^2 - (\psi \varphi + \, \nu ^2) \lambda ^2} \, H (\mu) 
\right)  
\notag \\[1ex]
%
& + \frac{1}{(\lambda - \mu)(\lambda + \mu)}  \bigg( \frac{\left( 2 ( \xi - \nu \mu) \, ( \xi + \nu \lambda) - \psi \varphi \, (\lambda + \mu) \mu \right) \mu}{\xi ^2 - (\psi \varphi + \nu ^2) \mu^2} \, F (\lambda) 
\notag \\[1ex]
%
&- \frac{2 \left( \xi ^2 - \left( \psi \varphi \, \mu + \nu ^2 \lambda \right) \lambda \right) \lambda}{\xi ^2 - (\psi \varphi + \, \nu ^2) \lambda ^2} \, F (\mu) \bigg)  \>,
\\[2.5ex]
%
%
\label{Eq: cr-FE}
\left[ F (\lambda) , E (\mu)  \right]  &=  \frac{-2}{\lambda + \mu} \left( \frac{\varphi \, (\xi + \nu \mu ) \, \mu}{\xi ^2 - (\psi \varphi + \nu ^2) \, \mu^2} \, F (\lambda) - \frac{\psi \, (\xi - \nu \lambda) \, \lambda}{\xi ^2 - (\psi \varphi + \nu ^2) \lambda ^2} \, E (\mu) \right)  
\notag \\[1ex]
%
& + \frac{2}{(\lambda - \mu)(\lambda + \mu)} \bigg( \frac{\left(2 (\xi - \nu \lambda) \, ( \xi + \nu \mu) - \psi \varphi \, (\lambda + \mu) \mu \right) \mu}{\xi ^2 - (\psi \varphi + \nu ^2) \mu^2} H (\lambda) 
\notag \\[1ex]
%
& - \frac{\left(2 (\xi - \nu \lambda) \, ( \xi + \nu \mu) - \psi \varphi \, (\lambda + \mu) \lambda \right)\lambda}{\xi ^2 - (\psi \varphi + \nu ^2) \lambda ^2} H (\mu) \bigg)  \>. 
\end{align} 



Moreover, the Maillet linear bracket \eqref{Eq: M-br-rk} yields the expression for the generating function of the $so(3)$ Gaudin Hamiltonians with general boundary terms in terms of the Lax operator \eqref{Eq: cal-L-def}
\begin{equation}
\label{Eq: tau}
\tau (\lambda) = \frac{1}{2} \, \mathrm{tr}_0 \left( \mathcal{L}_0 ^2(\lambda) \right) \>.
\end{equation}
Namely, using the Maillet bracket \eqref{Eq: M-br-rk}, it is straightforward to check that the operator $\tau (\lambda)$ commutes for different values of the spectral parameter,
\begin{equation}
\label{Eq: tau-tau}
[ \tau (\lambda) , \tau (\mu) ] = 0 \>.
\end{equation}
From \eqref{Eq: cal-L} it follows that
\begin{equation}
\label{Eq: tau-gen-inital}
\tau (\lambda) = H^2 (\lambda) + \frac{1}{2} \left(  E (\lambda) F (\lambda) + F (\lambda) E(\lambda) \right) \>.
\end{equation}


Our aim here is to obtain the spectrum and the corresponding states of the generating function $\tau (\lambda)$ by algebraic methods. To this end we would have to use the relations \eqref{Eq: cr-EE} -- \eqref{Eq: cr-FE}. As it is evident from the formulae above, these relations do not seem to be suitable to efficiently address this problem. Therefore, we propose a new set of generators
\begin{align}
\label{Eq: e-scr}
\mathscr{E} (\lambda) =& \frac{1}{2\psi \sqrt{\psi \varphi + \nu^2}} \ \bigg( \psi^2  E(\lambda) - \left( \psi \varphi + 2 \nu \left( \nu - \sqrt{\psi \varphi + \nu^2}\right) \right) F(\lambda) 
\notag \\[1ex]
%
&+ 2 \psi \left( \nu - \sqrt{\psi \varphi + \nu^2} \right) H(\lambda) \bigg)
\>, 
\\[2ex]
%
%
\label{Eq: f-scr}
\mathscr{F} (\lambda) &= \frac{1}{2\psi \sqrt{\psi \varphi + \nu^2}} \ \bigg( - \psi^2  E(\lambda) + \left( \psi \varphi + 2 \nu \left( \nu + \sqrt{\psi \varphi + \nu^2}\right) \right) F(\lambda) 
\notag \\[1ex]
%
&- 2 \psi \left( \nu + \sqrt{\psi \varphi + \nu^2} \right) H(\lambda) \bigg)
\>, \\[2ex]
%
%
\label{Eq: h-scr}
\mathscr{H} (\lambda) &= \frac{1}{2\sqrt{\psi \varphi + \nu^2}} \ \left( \psi E(\lambda) + \varphi F(\lambda) + 2\nu H(\lambda) \right)
 \>.
\end{align}
The commutation relations we obtain for the new generators are substantially simpler than the initial relations \eqref{Eq: cr-HE}--\eqref{Eq: cr-HH}. In particular, 
\begin{equation}
\label{Eq: cr-zero}
%&
\left[ \mathscr{E} (\lambda) , \mathscr{E} (\mu)  \right] = \left[\mathscr{F} (\lambda) , \mathscr{F} (\mu)  \right] = \left[ \mathscr{H} (\lambda) , \mathscr{H} (\mu)  \right] = 0 \>,
\end{equation}
and the three non-trivial relations are
\begin{align}
\label{Eq: cr-HE-scr}
&\left[ \mathscr{H} (\lambda) , \mathscr{E} (\mu)  \right] = \frac{-2}{\lambda ^2 - \mu ^2} \left( \mu \ \frac{\xi - \lambda \sqrt{\psi \varphi + \nu ^2}}{\xi - \mu \sqrt{\psi \varphi + \nu ^2}} \ \mathscr{E} (\lambda) - \lambda \ \mathscr{E} (\mu) \right) \>,
\\[2ex]
%
%
\label{Eq: cr-HF-scr} 
&\left[  \mathscr{H} (\lambda) , \mathscr{F} (\mu)  \right] = \frac{2}{\lambda ^2 - \mu ^2} \left( \mu \ \frac{\xi + \lambda \sqrt{\psi \varphi + \nu ^2}}{\xi + \mu \sqrt{\psi \varphi + \nu ^2}} \ \mathscr{F} (\lambda) - \lambda  \, \mathscr{F} (\mu) \right)  \>,
\\[2ex]
%
%
\label{Eq: cr-EF-scr}
&\left[ \mathscr{F} (\lambda) , \mathscr{E} (\mu)  \right] = \frac{4}{\lambda ^2 - \mu ^2} \left(\mu \ \frac{\xi - \lambda \sqrt{\psi \varphi + \nu ^2}}{\xi - \mu \sqrt{\psi \varphi + \nu ^2}}  \ \mathscr{H} (\lambda) - \lambda  \ \frac{\xi + \mu \sqrt{\psi \varphi + \nu ^2}}{\xi + \lambda \sqrt{\psi \varphi + \nu ^2}}  \ \mathscr{H} (\mu) \right)  \>.
\end{align} 




Furthermore, the local realization of the new generators is:
\begin{align}
\label{Eq: e-scr-l}
\mathscr{E} (\lambda) &= \frac{\lambda}{\sqrt{\psi \varphi + \nu^2}} \ \sum _{m=1}^N \frac{\xi  - \alpha_m \sqrt{\psi \varphi + \nu^2}}{\xi - \lambda \sqrt{\psi \varphi + \nu ^2}} 
\notag \\[1ex]
%
&\times  \frac{2 ( \nu - \sqrt{\psi \varphi + \nu^2} ) S^3_m + \psi S^{+}_m  - \frac{\psi \varphi + 2 \nu \, ( \nu - \sqrt{\psi \varphi + \nu^2})}{\psi}   S^{-} _m}{(\lambda - \alpha_m)(\lambda + \alpha_m)}  
  \>, \\[2ex]
%
%
\label{Eq: f-scr-l}
\mathscr{F} (\lambda) &= \frac{-\lambda}{\sqrt{\psi \varphi + \nu^2} } \ \sum _{m=1}^N \frac{\xi  + \alpha_m \sqrt{\psi \varphi + \nu^2}}{\xi + \lambda \sqrt{\psi \varphi + \nu ^2}} 
\notag \\[1ex]
%
&\times \frac{2 ( \nu + \sqrt{\psi \varphi + \nu^2} ) S^3_m + \psi S^{+}_m  - \frac{\psi \varphi + 2 \nu \, ( \nu + \sqrt{\psi \varphi + \nu^2})}{\psi}   S^{-} _m}{(\lambda - \alpha_m)(\lambda + \alpha_m)}  
\>, 
\\[2ex]
%
%
\label{Eq: h-scr-l}
\mathscr{H} (\lambda) &= \frac{\lambda}{\sqrt{\psi \varphi + \nu^2}} \ \sum _{m=1}^N \frac{2\nu S^3_m + \psi S^{+}_m  + \varphi  S^{-} _m}{(\lambda - \alpha_m)(\lambda + \alpha_m)}
\>.
\end{align}

A straightforward but somewhat lengthy calculation shows that the generating function $\tau (\lambda)$ \eqref{Eq: tau-gen-inital} has exactly the same form when expressed in terms of the new generators
\begin{equation}
\label{Eq: tau-gen}
\tau (\lambda) = \mathscr{H}^2 (\lambda) + \frac{1}{2} \left(  \mathscr{E} (\lambda) \mathscr{F} (\lambda) + \mathscr{F} (\lambda) \mathscr{E}(\lambda) \right) \>.
\end{equation}

The explicit expressions for the $so(3)$ Gaudin Hamiltonians with general boundary terms are derived by substituting the the local realization of the new generators \eqref{Eq: e-scr-l} -- \eqref{Eq: h-scr-l} in the right-hand-side of \eqref{Eq: tau-gen}
\begin{align}
\label{Eq: Ham-m}
\mathrm{H}_m =& (\pm) \ \underset{\lambda = \pm \alpha_m}{\Res} \tau (\lambda) = \frac{1}{\xi ^2 - \left( \psi \varphi + \nu ^2 \right) \alpha_m^2}
\notag \\[1.5ex]
%
&\times \bigg( \frac{\xi ^2 + \left( \psi \varphi - \nu ^2 \right) \alpha_m^2}{\alpha_m} \ (S^3_m ) ^2 - \frac{\alpha_m}{2} \left( \psi ^2  \ (S_m^{+})^2 + \varphi ^2 \ (S_m^{-})^2 \right.
\notag \\[1.5ex]
%
& \left.+ 2 \psi \nu \left( S_m^+ S_m^3 +  S_m^3 S_m^+ \right) + 2 \varphi  \nu \left( S_m^- S_m^3 + S_m^3 S_m^- \right) \right)
\notag \\[1.5ex]
%
&+ \frac{\xi ^2 + \nu^2 \alpha_m^2}{2 \alpha_m}  \left( S_m^{+} S_m^{-} + S_m^{-} S_m^{+} \right)
\bigg) 
\notag \\[1.5ex]
%
&+ \frac{\alpha_m }{\xi ^2 - \left( \psi \varphi + \nu ^2 \right) \alpha_m^2} \ \sum _{n \neq m}^{N} \bigg( \frac{4 \left( \xi ^2 - \psi \varphi \, \alpha_m \alpha_n - \nu ^2 \alpha_m^2 \right)}{\alpha_m^2 -\alpha_n^2}  \ S_m^3 S_n^3  
\notag \\[1.5ex]
%
&- \frac{\alpha_m}{\alpha_m + \alpha_n} \big(  
\psi ^2  \ S_m^+ S_n^+ + \varphi ^2 \ S_m^- S_n^- 
\notag \\[1.5ex]
%
&  + 2 \psi \nu \left( S_m^+ S_n^3 +  S_m^3 S_n^+ \right) + 2 \varphi  \nu \left( S_m^- S_n^3 + S_m^3 S_n^- \right) \big) 
\notag \\[1.5ex]
%
&+ \frac{2 \left( \xi ^2 - \left( \psi \varphi + \nu ^2 \right)  \alpha_m \alpha_n \right) - \psi \varphi \, \alpha_m \, (\alpha_m -\alpha_n)}{\alpha_m^2 -\alpha_n^2} \left( S_m^- S_n^+ + S_m^+ S_n^- \right) \bigg)
\notag \\[1.5ex]
%   
%\left. 
%
& + \frac{\xi \cdot \alpha_m} {\xi ^2 - \left( \psi \varphi + \nu ^2 \right) \alpha_m^2} \sum _{n \neq m}^{N} \frac{1}{\alpha_m +\alpha_n} \Big( 2 \psi \left( S_m^+ S_n^3 -  S_m^3 S_n^+ \right) 
\notag \\[1.5ex]
%
&+ 2 \nu \left( S_m^- S_n^+ - S_m^+ S_n^- \right) + 2 \varphi \left( S_m^3 S_n^- - S_m^- S_n^3\right) \Big) \>.
\end{align}

Besides the formula above for the $so(3)$ Gaudin Hamiltonians with general boundary terms, our main result in this section is the new form of generators of the generalized $so(3)$ Gaudin algebra \eqref{Eq: e-scr-l} -- \eqref{Eq: h-scr-l}. Due to their strikingly simple  commutation relations \eqref{Eq: cr-zero} -- \eqref{Eq: cr-EF-scr} they now provide a suitable framework for applying the algebraic Bethe ansatz without any restrictions on boundary parameters.




%
%
\section{Implementation of the algebraic Bethe ansatz \label{Sec: ABA}}
%
%
Before we can proceed to find Bethe vectors and determine the off-shell action of the generating function $\tau(\lambda)$, we have to establish several intermediary results.

In the Hilbert space $\mathcal{H}$ \eqref{Eq: H-space} of the system we have to define the so-called vacuum vector $\Omega _+ \in \mathcal{H}$ together with the appropriate action of the generators \eqref{Eq: e-scr-l} -- \eqref{Eq: h-scr-l} on it. To this purpose we observe that in
every local space $V_ m =  \mathbb{C}^{3}$, $m \in \{ 1 , \ldots , N \}$ there exists a vector $ \omega_m \in V_ m$ given by
\begin{equation}
\label{Eq: om}
\omega_m = \left(\begin{array}{c} 
\psi ^2 \\
-\sqrt{2} \, \psi \left( \nu - \sqrt{\psi \varphi  + \nu ^2} \right) \\
\left( \nu - \sqrt{\psi \varphi + \nu^2} \right) ^2
\end{array}\right)
\in \mathbb{C}^{3} = V_ m \>,
\end{equation}
where the parameters  $\nu$, $\psi$ and $\varphi$ are the parameters of the boundary  K-matrix \eqref{Eq: K-mat}. Then it is easy to check that
\begin{align}
\label{Eq: S+-om}
&\left( 2 \left( \nu - \sqrt{\psi \varphi + \nu^2} \right) S^3_m + \psi \ S^{+}_m  - \frac{\psi \varphi + 2 \nu \, ( \nu - \sqrt{\psi \varphi + \nu^2})}{\psi}  \ S^{-}_m \right) \omega_m = 0 \> , \\[2ex]
%
&\left( 2\nu S^3_m + \psi S^{+}_m  + \varphi  S^{-} _m \right) \omega_m = 2 \sqrt{\psi \varphi + \nu ^2} \ \omega_m \>.
% 
\end{align}
Therefore the vacuum vector $\Omega _+$, defined as
\begin{equation}
\label{Eq: Omega+}
\Omega _+ = \omega _1 \otimes \cdots \otimes \omega _N \in \mathcal{H}
\end{equation} 
has the desired properties. Namely, it is annihilated by the generator $\mathscr{E} (\lambda) $ \eqref{Eq: e-scr-l} and, at the same time, it is an eigenvector of the generator $\mathscr{H} (\lambda)$ \eqref{Eq: h-scr-l}, that is
\begin{equation}
\label{Eq: Vac}
\mathscr{E} (\lambda) \ \Omega _+ = 0 \quad \text{and} \quad \mathscr{H} (\lambda) \ \Omega _+ = \rho (\lambda) \ \Omega _+ \quad \text{with} \quad \rho (\lambda) = \sum _{m=1}^N \frac{2\lambda}{\lambda ^2 - \alpha _m^2} \>.
\end{equation}



Our next aim is to rewrite the formula for $\tau (\lambda)$ \eqref{Eq: tau-gen} in a more suitable way so that the action of the generating function $\tau (\lambda)$ on the vacuum vector $\Omega _+$ \eqref{Eq: Omega+} becomes more transparent. With this aim, we first note that the commutation relations \eqref{Eq: cr-HE-scr} -- \eqref{Eq: cr-EF-scr} imply
\begin{align}
\label{Eq: cr-HF-scr-la}
\left[  \mathscr{H} (\lambda) , \mathscr{F} (\lambda)  \right] &= \frac{-\xi}{\lambda \left(\xi + \lambda \sqrt{\psi \varphi + \nu ^2}\right)} \ \mathscr{F} (\lambda) + \ \mathscr{F}^{\prime} (\lambda)  \>,
\\[2.5ex]
%
%
\label{Eq: cr-HE-scr-la}
\left[  \mathscr{H} (\lambda) , \mathscr{E} (\lambda)  \right] &= \frac{\xi}{\lambda \left(\xi - \lambda \sqrt{\psi \varphi + \nu ^2}\right)} \ \mathscr{E} (\lambda) - \ \mathscr{E}^{\prime} (\lambda)  \>,
\\[2.5ex]
%
%
\label{Eq: cr-EF-scr-lam}
\left[ \mathscr{F} (\lambda) , \mathscr{E} (\lambda)  \right] &= 2 \left( \frac{-1}{\lambda} \ \frac{ \xi^2 + (\psi \varphi + \nu ^2) \lambda^2}{ \xi^2 - (\psi \varphi + \nu ^2) \lambda^2}  \ \mathscr{H} (\lambda)  +  \ \mathscr{H}^{\prime} (\lambda) \right) \>,
\end{align}
where prime denotes derivative with respect to parameter.
Therefore we can express the generating function $\tau (\lambda)$ \eqref{Eq: tau-gen} as follows
\begin{equation}
\label{Eq: tau-g}
\tau (\lambda) = \mathscr{H}^2 (\lambda) + \frac{1}{\lambda} \ \frac{ \xi^2 + (\psi \varphi + \nu ^2) \lambda^2}{ \xi^2 - (\psi \varphi + \nu ^2) \lambda^2}  \ \mathscr{H} (\lambda) - \ \mathscr{H}^{\prime} (\lambda) + \mathscr{F} (\lambda) \mathscr{E}(\lambda) \>.
\end{equation}
Taking into account \eqref{Eq: Vac} and \eqref{Eq: tau-g}, it is evident that  the vacuum vector $\Omega _+$ \eqref{Eq: Omega+} is an eigenvector of the generating function 
\begin{equation}
\label{Eq: tau-Om}
\tau (\lambda) \, \Omega _+ = \chi_0 (\lambda) \, \Omega _+  \quad \text{with} \quad \chi_0 (\lambda) = \rho ^2 (\lambda) + \frac{ \xi^2 + (\psi \varphi + \nu ^2) \lambda^2}{ \xi^2 - (\psi \varphi + \nu ^2) \lambda^2}  \ \frac{\rho (\lambda)}{\lambda}   -  \rho^{\prime} (\lambda) \>.
\end{equation}
 

In our approach, one fo the essential steps in the implementation of algebraic Bethe ansatz is to find the commutation relation between the generating function $\tau (\lambda)$ \eqref{Eq: tau-g} and the generator $\mathscr{F} (\mu)$ \eqref{Eq: f-scr}. To this end, we will also need the following auxiliary result which follows from \eqref{Eq: cr-HF-scr}
\begin{equation}
\begin{split}
\label{Eq: cr-HprimF-scr-la}
\left[  \mathscr{H}^{\prime} (\lambda) , \mathscr{F} (\mu)  \right] &= \frac{2}{\lambda ^2 - \mu^2} \left( \frac{\xi \ (\lambda - \mu)}{(\lambda + \mu) \left( \xi + \mu \sqrt{\psi \varphi + \nu^2} \right)} \ \mathscr{F} (\lambda) \right.
\notag \\[1ex]
%
&+\left. \frac{\lambda ^2 + \mu^2}{\lambda ^2 - \mu^2} \left( \mathscr{F} (\mu) - \mathscr{F} (\lambda) \right)  
  + \mu \ \frac{\xi + \lambda \sqrt{\psi \varphi + \nu ^2}}{\xi + \mu \sqrt{\psi \varphi + \nu ^2}} \ \mathscr{F}^{\prime} (\lambda) 
\right)
\>.
\end{split}
\end{equation} 

Now we can compute the commutator by a straightforward calculation, based on the formulae \eqref{Eq: tau-g}, \eqref{Eq: cr-HF-scr}, \eqref{Eq: cr-EF-scr}  and \eqref{Eq: cr-HprimF-scr-la}:
\begin{equation}
\label{Eq: tau-f}
\begin{split}
&\left[ \tau (\lambda) , \mathscr{F} (\mu) \right] = - \frac{4}{\lambda ^2 - \mu^2} \ \mathscr{F} (\mu) \left( \lambda \mathscr{H}(\lambda) + \frac{(\psi \varphi + \nu ^2) \lambda^2}{\xi^2 - (\psi \varphi + \nu ^2) \lambda^2} \right) \\[1ex]
%
&+ \frac{4}{\lambda ^2 - \mu^2} \ \frac{\lambda}{\mu} \ \frac{\xi - \mu \sqrt{\psi \varphi + \nu ^2}}{\xi - \lambda \sqrt{\psi \varphi + \nu ^2}}
 \ \mathscr{F} (\lambda) \left( \mu \mathscr{H}(\mu) + \frac{(\psi \varphi + \nu ^2) \mu^2}{\xi^2 - (\psi \varphi + \nu ^2) \mu^2} \right) \>.
\end{split} 
\end{equation}


The relative simplicity of the right hand side of the equation above has encouraged us to seek the commutator between the operator $\tau (\lambda)$ and the product $\mathscr{F} (\mu _1) \mathscr{F} (\mu _2)$ as the next step. In this case, an analogous direct calculation based on the previous formulae, leads to
\begin{equation}
\label{Eq: tau-ff}
\begin{split}
[ \tau (\lambda), & \mathscr{F} (\mu _1) \mathscr{F} (\mu _2)] =\\[1ex]
%
&- \frac{4}{\lambda ^2 - \mu_1^2} \ \mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \left( \lambda \mathscr{H}(\lambda) + \frac{(\psi \varphi + \nu ^2) \lambda^2}{\xi^2 - (\psi \varphi + \nu ^2) \lambda^2} - \frac{\lambda ^2}{\lambda ^2 - \mu_2^2} \right) \\[1ex]
%
&- \frac{4}{\lambda ^2 - \mu_2^2} \ \mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \left(\lambda \mathscr{H}(\lambda) + \frac{(\psi \varphi + \nu ^2) \lambda^2}{\xi^2 - (\psi \varphi + \nu ^2) \lambda^2} - \frac{\lambda ^2}{\lambda ^2 - \mu_1^2} \right) \\[1ex]
%
&+ \frac{4}{\lambda ^2 - \mu _1^2} \ \frac{\lambda}{\mu _1} \ \frac{\xi - \mu _1 \sqrt{\psi \varphi + \nu ^2}}{\xi - \lambda \sqrt{\psi \varphi + \nu ^2}} \ \mathscr{F} (\lambda) \mathscr{F} (\mu _2) \\[1ex]
%
&\times \left( \mu _1 \mathscr{H}(\mu _1) + \frac{(\psi \varphi + \nu ^2) \mu _1^2}{\xi^2 - (\psi \varphi + \nu ^2) \mu _1^2} - \frac{2 \mu _1^2}{\mu _1^2 - \mu_2^2} \right) \\[1ex]
%
&+ \frac{4}{\lambda ^2 - \mu _2^2} \ \frac{\lambda}{\mu _2} \ \frac{\xi - \mu _2 \sqrt{\psi \varphi + \nu ^2}}{\xi - \lambda \sqrt{\psi \varphi + \nu ^2}} \ \mathscr{F} (\mu _1) \mathscr{F} (\lambda)\\[1ex]
%
&\times \left( \mu _2 \mathscr{H}(\mu _2) + \frac{(\psi \varphi + \nu ^2) \mu _2^2}{\xi^2 - (\psi \varphi + \nu ^2) \mu _2^2} - \frac{2 \mu _2^2}{\mu _2^2 - \mu_1^2} \right) \>.
\end{split} 
\end{equation}



Evidently, the right hand side of \eqref{Eq: tau-ff} has extra lines and every line has extra terms in comparison with \eqref{Eq: tau-f}. While certain pattern is already visible, we will explicitly compute one more step befor conjecturing the general case. The commutation relation between the generating function $\tau (\lambda)$ and the product $\mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \mathscr{F} (\mu _3)$ is obtained in a similar manner, using the previous results,
\begin{equation}
\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!
\label{Eq: tau-fff}
\begin{split}
&\left[ \tau (\lambda) , \mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \mathscr{F} (\mu _3) \right] =  \\[1ex]
%
& - \frac{4}{\lambda ^2 - \mu_1^2} \ \mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \mathscr{F} (\mu _3) \left( \lambda \mathscr{H}(\lambda) + \frac{(\psi \varphi + \nu ^2) \lambda^2}{\xi^2 - (\psi \varphi + \nu ^2) \lambda^2} - \frac{\lambda ^2}{\lambda ^2 - \mu_2^2} - \frac{\lambda ^2}{\lambda ^2 - \mu_3^2} \right) \\[1ex]
%
&- \frac{4}{\lambda ^2 - \mu_2^2} \ \mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \mathscr{F} (\mu _3) \left(\lambda \mathscr{H}(\lambda) + \frac{(\psi \varphi + \nu ^2) \lambda^2}{\xi^2 - (\psi \varphi + \nu ^2) \lambda^2} - \frac{\lambda ^2}{\lambda ^2 - \mu_1^2}  - \frac{\lambda ^2}{\lambda ^2 - \mu_3^2}  \right) \\[1ex]
%
&- \frac{4}{\lambda ^2 - \mu_3^2} \ \mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \mathscr{F} (\mu _3) \left(\lambda \mathscr{H}(\lambda) + \frac{(\psi \varphi + \nu ^2) \lambda^2}{\xi^2 - (\psi \varphi + \nu ^2) \lambda^2} - \frac{\lambda ^2}{\lambda ^2 - \mu_1^2}  - \frac{\lambda ^2}{\lambda ^2 - \mu_2^2}  \right) \\[1ex]
%
&+ \frac{4}{\lambda ^2 - \mu _1^2} \ \frac{\lambda}{\mu _1} \ \frac{\xi - \mu _1 \sqrt{\psi \varphi + \nu ^2}}{\xi - \lambda \sqrt{\psi \varphi + \nu ^2}} \ \mathscr{F} (\lambda) \mathscr{F} (\mu _2) \mathscr{F} (\mu _3)\\[1ex]
%
&\qquad \times \left( \mu _1 \mathscr{H}(\mu _1) + \frac{(\psi \varphi + \nu ^2) \mu _1^2}{\xi^2 - (\psi \varphi + \nu ^2) \mu _1^2} - \frac{2 \mu _1^2}{\mu _1^2 - \mu_2^2} - \frac{2 \mu _1^2}{\mu _1^2 - \mu_3^2}\right) \\[1ex]
%
&+ \frac{4}{\lambda ^2 - \mu _2^2} \ \frac{\lambda}{\mu _2} \ \frac{\xi - \mu _2 \sqrt{\psi \varphi + \nu ^2}}{\xi - \lambda \sqrt{\psi \varphi + \nu ^2}} \ \mathscr{F} (\mu _1) \mathscr{F} (\lambda)  \mathscr{F} (\mu _3)\\[1ex]
%
&\qquad \times \left( \mu _2 \mathscr{H}(\mu _2) + \frac{(\psi \varphi + \nu ^2) \mu _2^2}{\xi^2 - (\psi \varphi + \nu ^2) \mu _2^2} - \frac{2 \mu _2^2}{\mu _2^2 - \mu_1^2}  - \frac{2 \mu _2^2}{\mu _2^2 - \mu_3^2} \right) 
\\[1ex]
%
&+ \frac{4}{\lambda ^2 - \mu _3^2} \ \frac{\lambda}{\mu _3} \ \frac{\xi - \mu _3 \sqrt{\psi \varphi + \nu ^2}}{\xi - \lambda \sqrt{\psi \varphi + \nu ^2}} \ \mathscr{F} (\mu _1) \mathscr{F} (\mu _2 )  \mathscr{F} (\lambda )\\[1ex]
%
&\qquad \times \left( \mu _3 \mathscr{H}(\mu _3) + \frac{(\psi \varphi + \nu ^2) \mu _3^2}{\xi^2 - (\psi \varphi + \nu ^2) \mu _3^2} - \frac{2 \mu _3^2}{\mu _3^2 - \mu_1^2}  - \frac{2 \mu _3^2}{\mu _3^2 - \mu_2^2} \right) \>.
\end{split} 
\end{equation}

In the general case, we conjecture validity of the following relation:
\begin{equation}
\label{Eq: tau-Mf}
\begin{split}
&\left[ \tau (\lambda) , \mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \cdots \mathscr{F} (\mu _M) \right] = - \mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \cdots \mathscr{F} (\mu _M) \sum _{j=1}^M \frac{4}{\lambda ^2 - \mu_j^2} \times 
\\[1.5ex]
%
&\times \left( \lambda \mathscr{H}(\lambda) + \frac{(\psi \varphi + \nu ^2) \lambda^2}{\xi^2 - (\psi \varphi + \nu ^2) \lambda^2} - \sum _{k \neq j}^M \frac{\lambda ^2}{\lambda ^2 - \mu_k^2} \right) \\[1.5ex]
%
&+ \sum _{j=1}^M \frac{4}{\lambda ^2 - \mu _j^2} \ \frac{\lambda}{\mu _j} \ \frac{\xi - \mu _j \sqrt{\psi \varphi + \nu ^2}}{\xi - \lambda \sqrt{\psi \varphi + \nu ^2}} \mathscr{F} (\lambda) \mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \cdots \widehat{\mathscr{F} (\mu _j)} \cdots
\mathscr{F} (\mu _M)\\[1.5ex]
%
&\times \left( \mu _j \mathscr{H}(\mu _j) + \frac{(\psi \varphi + \nu ^2) \mu _j^2}{\xi^2 - (\psi \varphi + \nu ^2) \mu _j^2} - \sum _{k \neq j}^M \frac{2 \mu _j^2}{\mu _j^2 - \mu_k^2} \right) 
 \>,
\end{split} 
\end{equation}
where the notation $\widehat{\mathscr{F} (\mu _j)}$ means that the operator $\mathscr{F} (\mu _j)$ is omitted.


The proof of the formula above is by mathematical induction. Our initial hypothesis is that  the equation \eqref{Eq: tau-Mf} is valid for some natural number $M$. Thus, we have to show that the analogous equation is valid for $M+1$. To this end, we write  
\begin{align}
\label{Eq: tau-M+1f-s1}
[ \tau (\lambda) , &\mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \cdots \mathscr{F} (\mu _M)\mathscr{F} (\mu _{M+1}) ] = 
\notag \\[1ex]
%
&=\left[ \tau (\lambda) , \mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \cdots \mathscr{F} (\mu _M) \right] \mathscr{F} (\mu _{M+1}) 
\notag \\[1ex]
%
&+ \mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \cdots \mathscr{F} (\mu _M)\mathscr{F} (\mu _{M+1}) \left[ \tau (\lambda) , \mathscr{F} (\mu _{M+1}) \right] 
\\[1.5ex]
%
\label{Eq: tau-M+1f-s2}
&=  \left[ \left[   \tau (\lambda) , \mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \cdots \mathscr{F} (\mu _M) \right] , \mathscr{F} (\mu _{M+1}) \right] 
\notag \\[1ex]
%
&+ \mathscr{F} (\mu _{M+1}) \left[ \tau (\lambda) , \mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \cdots \mathscr{F} (\mu _M) \right]
\notag \\[1ex]
%
&+ \mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \cdots \mathscr{F} (\mu _M)\mathscr{F} (\mu _{M+1}) \left[ \tau (\lambda) , \mathscr{F} (\mu _{M+1}) \right] \>.
\end{align}
It follows from \eqref{Eq: tau-Mf} that the first term on the right hand side of \eqref{Eq: tau-M+1f-s2} yields two type of terms. In the second term of \eqref{Eq: tau-M+1f-s2} we can just substitute the right hand side of \eqref{Eq: tau-Mf}. Finally, in the last term of  \eqref{Eq: tau-M+1f-s2} we use \eqref{Eq: tau-f}. In this way we obtain
\small
\begin{align}
\label{Eq: tau-M+1f-s3}
&\left[ \tau (\lambda) , \mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \cdots \mathscr{F} (\mu _M) \mathscr{F} (\mu _{M+1}) \right] = 
\notag \\[1ex]
%
&- \mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \cdots \mathscr{F} (\mu _M) \sum _{j=1}^M \frac{4\lambda}{\lambda ^2 - \mu_j^2} \left[
\mathscr{H} (\lambda) , \mathscr{F} (\mu _{M+1}) \right] 
\notag \\[1ex]
%
&+ \sum _{j=1}^M \frac{4\lambda}{\lambda ^2 - \mu _j^2}  \ \frac{\xi - \mu _j \sqrt{\psi \varphi + \nu ^2}}{\xi - \lambda \sqrt{\psi \varphi + \nu ^2}} \times
\notag \\[1ex] 
%
&\qquad \times \mathscr{F} (\lambda) \mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \cdots \widehat{\mathscr{F} (\mu _j)} \cdots
\mathscr{F} (\mu _M) \left[ \mathscr{H} (\mu _j) , \mathscr{F} (\mu _{M+1}) \right] 
\notag \\[1ex]
%
&- \mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \cdots \mathscr{F} (\mu _M) \mathscr{F} (\mu _{M+1}) \times 
\notag \\[1ex]
%
&\qquad \times  \sum _{j=1}^M \frac{4}{\lambda ^2 - \mu_j^2}  \left( \lambda \mathscr{H}(\lambda) + \frac{(\psi \varphi + \nu ^2) \lambda^2}{\xi^2 - (\psi \varphi + \nu ^2) \lambda^2} - \sum _{k \neq j}^M \frac{\lambda ^2}{\lambda ^2 - \mu_k^2} \right) 
\notag \\[1ex]
%
&+ \sum _{j=1}^M \frac{4}{\lambda ^2 - \mu _j^2} \ \frac{\lambda}{\mu _j} \ \frac{\xi - \mu _j \sqrt{\psi \varphi + \nu ^2}}{\xi - \lambda \sqrt{\psi \varphi + \nu ^2}} \ \mathscr{F} (\lambda) \mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \cdots \widehat{\mathscr{F} (\mu _j)} \cdots
\mathscr{F} (\mu _M)  \mathscr{F} (\mu _{M+1}) 
\times  
\notag \\[1ex]
%
&\qquad \times \left( \mu _j \mathscr{H}(\mu _j) + \frac{(\psi \varphi + \nu ^2) \mu _j^2}{\xi^2 - (\psi \varphi + \nu ^2) \mu _j^2} - \sum _{k \neq j}^M \frac{2 \mu _j^2}{\mu _j^2 - \mu_k^2} \right) 
\notag \\[1ex]
%
&+ \mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \cdots \mathscr{F} (\mu _M) \Bigg( - \frac{4}{\lambda ^2 - \mu _{M+1}^2} \ \mathscr{F} (\mu _{M+1}) \left( \lambda \mathscr{H}(\lambda) + \frac{(\psi \varphi + \nu ^2) \lambda^2}{\xi^2 - (\psi \varphi + \nu ^2) \lambda^2} \right) 
\notag \\[1ex]
%
& + \frac{4}{\lambda ^2 - \mu _{M+1}^2} \ \frac{\lambda}{\mu _{M+1}} \ \frac{\xi - \mu _{M+1} \sqrt{\psi \varphi + \nu ^2}}{\xi - \lambda \sqrt{\psi \varphi + \nu ^2}}
 \ \mathscr{F} (\lambda) 
\notag \times \\[1ex]
%
&\qquad \times \left( \mu _{M+1} \mathscr{H}(\mu _{M+1}) + \frac{(\psi \varphi + \nu ^2) \mu _{M+1}^2}{\xi^2 - (\psi \varphi + \nu ^2) \mu _{M+1}^2} \right) \Bigg) \>.
\end{align}
\normalsize

In the first two terms on the right hand side of \eqref{Eq: tau-M+1f-s3} we used the equation \eqref{Eq: cr-HF-scr} and the remaining  terms we rewrite in a more appropriate order
\small
\begin{align}
\label{Eq: tau-M+1f-s4}
&\left[ \tau (\lambda) , \mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \cdots \mathscr{F} (\mu _M) \mathscr{F} (\mu _{M+1}) \right] = 
- \mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \cdots \mathscr{F} (\mu _M) \sum _{j=1}^M \frac{4\lambda}{\lambda ^2 - \mu_j^2} \times
\notag \\[1ex]
% \left[\mathscr{H} (\lambda) , \mathscr{F} (\mu _{M+1}) \right] 
&\qquad \times \frac{2}{\lambda ^2 - \mu_{M+1} ^2} \left( \mu_{M+1} \ \frac{\xi + \lambda \sqrt{\psi \varphi + \nu ^2}}{\xi + \mu _{M+1} \sqrt{\psi \varphi + \nu ^2}} \ \mathscr{F} (\lambda) - \lambda  \, \mathscr{F} (\mu _{M+1}) \right)
\notag \\[1ex]
%
&+ \sum _{j=1}^M \frac{4\lambda}{\lambda ^2 - \mu _j^2}  \ \frac{\xi - \mu _j \sqrt{\psi \varphi + \nu ^2}}{\xi - \lambda \sqrt{\psi \varphi + \nu ^2}} \ \mathscr{F} (\lambda) \mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \cdots \widehat{\mathscr{F} (\mu _j)} \cdots
\mathscr{F} (\mu _M)  \times
\notag \\[1ex]
%
&\qquad \times \frac{2}{\mu _{j} ^2 - \mu_{M+1} ^2} \left( \mu_{M+1} \ \frac{\xi + \mu _{j} \sqrt{\psi \varphi + \nu ^2}}{\xi + \mu _{M+1} \sqrt{\psi \varphi + \nu ^2}} \ \mathscr{F} (\mu _{j}) - \mu _{j}  \, \mathscr{F} (\mu _{M+1}) \right)
\notag \\[1ex]
%
&- \mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \cdots \mathscr{F} (\mu _M) \mathscr{F} (\mu _{M+1})\times
\notag \\[1ex]
%
&\qquad \times  \sum _{j=1}^M \frac{4}{\lambda ^2 - \mu_j^2}  \left( \lambda \mathscr{H}(\lambda) + \frac{(\psi \varphi + \nu ^2) \lambda^2}{\xi^2 - (\psi \varphi + \nu ^2) \lambda^2} - \sum _{k \neq j}^M \frac{\lambda ^2}{\lambda ^2 - \mu_k^2} \right) 
\notag \\[1ex]
%
&- \mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \cdots \mathscr{F} (\mu _M) \mathscr{F} (\mu _{M+1}) \times
\notag \\[1ex]
%
&\qquad \times \frac{4}{\lambda ^2 - \mu _{M+1}^2}
\left( \lambda \mathscr{H}(\lambda) + \frac{(\psi \varphi + \nu ^2) \lambda^2}{\xi^2 - (\psi \varphi + \nu ^2) \lambda^2} \right)
\notag \\[1ex]
%
&+ \sum _{j=1}^M \frac{4}{\lambda ^2 - \mu _j^2} \ \frac{\lambda}{\mu _j} \ \frac{\xi - \mu _j \sqrt{\psi \varphi + \nu ^2}}{\xi - \lambda \sqrt{\psi \varphi + \nu ^2}} \ \mathscr{F} (\lambda) \mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \cdots \widehat{\mathscr{F} (\mu _j)} \cdots
\mathscr{F} (\mu _M)  \mathscr{F} (\mu _{M+1}) 
\times  
\notag \\[1ex]
%
&\qquad \times \left( \mu _j \mathscr{H}(\mu _j) + \frac{(\psi \varphi + \nu ^2) \mu _j^2}{\xi^2 - (\psi \varphi + \nu ^2) \mu _j^2} - \sum _{k \neq j}^M \frac{2 \mu _j^2}{\mu _j^2 - \mu_k^2} \right) 
\notag \\[1ex]
%
&+ \frac{4}{\lambda ^2 - \mu _{M+1}^2} \ \frac{\lambda}{\mu _{M+1}} \ \frac{\xi - \mu _{M+1} \sqrt{\psi \varphi + \nu ^2}}{\xi - \lambda \sqrt{\psi \varphi + \nu ^2}} \ \mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \cdots \mathscr{F} (\mu _M) \mathscr{F} (\lambda) \times
\notag \\[1ex]
%
&\qquad \times \left( \mu _{M+1} \mathscr{H}(\mu _{M+1}) + \frac{(\psi \varphi + \nu ^2) \mu _{M+1}^2}{\xi^2 - (\psi \varphi + \nu ^2) \mu _{M+1}^2} \right) \>.
\end{align}
\normalsize

Now it is just a question of reordering the terms in a more suitable manner
\begin{align}
\label{Eq: tau-M+1f-s5}
&\left[ \tau (\lambda) , \mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \cdots \mathscr{F} (\mu _M) \mathscr{F} (\mu _{M+1}) \right] = 
\notag \\[1ex]
%
&- \mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \cdots \mathscr{F} (\mu _M) \mathscr{F} (\mu _{M+1}) \times
\notag \\[1ex]
%
&\qquad \times   \sum _{j=1}^M \frac{4}{\lambda ^2 - \mu_j^2}  \left( \lambda \mathscr{H}(\lambda) + \frac{(\psi \varphi + \nu ^2) \lambda^2}{\xi^2 - (\psi \varphi + \nu ^2) \lambda^2} - \sum _{k \neq j}^{M+1} \frac{\lambda ^2}{\lambda ^2 - \mu_k^2} \right) 
\notag \\[1ex]
%
&- \mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \cdots \mathscr{F} (\mu _M) \mathscr{F} (\mu _{M+1}) \times
\notag \\[1ex]
%
&\qquad \times  \frac{4}{\lambda ^2 - \mu _{M+1}^2}
\left( \lambda \mathscr{H}(\lambda) + \frac{(\psi \varphi + \nu ^2) \lambda^2}{\xi^2 - (\psi \varphi + \nu ^2) \lambda^2} - \sum _{k=1}^{M} \frac{\lambda ^2}{\lambda ^2 - \mu_k^2} 
\right)
\notag \\[1ex]
%
&+ \sum _{j=1}^M \frac{4}{\lambda ^2 - \mu _j^2} \ \frac{\lambda}{\mu _j} \ \frac{\xi - \mu _j \sqrt{\psi \varphi + \nu ^2}}{\xi - \lambda \sqrt{\psi \varphi + \nu ^2}} \ \mathscr{F} (\lambda) \mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \cdots \widehat{\mathscr{F} (\mu _j)} \cdots\notag \\[1ex]
%
& \cdots \mathscr{F} (\mu _M)  \mathscr{F} (\mu _{M+1}) 
\left( \mu _j \mathscr{H}(\mu _j) + \frac{(\psi \varphi + \nu ^2) \mu _j^2}{\xi^2 - (\psi \varphi + \nu ^2) \mu _j^2} - \sum _{k \neq j}^{M+1} \frac{2 \mu _j^2}{\mu _j^2 - \mu_k^2} \right) 
\notag \\[1ex]
%
&+ \frac{4}{\lambda ^2 - \mu _{M+1}^2} \ \frac{\lambda}{\mu _{M+1}} \ \frac{\xi - \mu _{M+1} \sqrt{\psi \varphi + \nu ^2}}{\xi - \lambda \sqrt{\psi \varphi + \nu ^2}} \ \mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \cdots \mathscr{F} (\mu _M) \mathscr{F} (\lambda) \times
\notag \\[1ex]
%
&\qquad\times \left( \mu _{M+1} \mathscr{H}(\mu _{M+1}) + \frac{(\psi \varphi + \nu ^2) \mu _{M+1}^2}{\xi^2 - (\psi \varphi + \nu ^2) \mu _{M+1}^2} 
- \sum _{k=1}^M \frac{2 \mu _{M+1}^2}{\mu _{M+1}^2 -\mu _{k}^2} \right) \>.
\end{align}
To obtain all the terms in the last sum (in the last line above) we had to use a generally valid, purely algebraic identity:
\begin{equation}
\label{Eq: ident}
\begin{split}
&\frac{-1}{(\lambda ^2 - \mu _j^2)(\lambda ^2 - \mu _{M+1}^2)} \ \frac{\xi^2 - (\psi \varphi + \nu ^2) \lambda ^2}{\xi^2 - (\psi \varphi + \nu ^2) \mu _{M+1}^2} + \frac{1}{(\lambda ^2 - \mu _j^2)(\mu _{j}^2 - \mu _{M+1}^2)} \times 
 \\[1ex]
%
&\times \frac{\xi^2 - (\psi \varphi + \nu ^2) \mu _{j}^2}{\xi^2 - (\psi \varphi + \nu ^2) \mu _{M+1}^2} = \frac{-1}{(\lambda ^2 - \mu _{M+1}^2 )(\mu _{M+1}^2 - \mu _j^2)} \>,
\end{split}
\end{equation}
which is valid for every $j = 1 , 2, \ldots , M$.

Finally, to complete the proof we can simply combine together similar terms and obtain the desired result
\begin{align}
\label{Eq: tau-M+1f-s6}
&\left[ \tau (\lambda) , \mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \cdots \mathscr{F} (\mu _M) \mathscr{F} (\mu _{M+1}) \right] = 
\notag \\[1ex]
%
&- \mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \cdots \mathscr{F} (\mu _M) \mathscr{F} (\mu _{M+1}) \times
\notag \\[1ex]
%
&\qquad \times    \sum _{j=1}^{M+1} \frac{4}{\lambda ^2 - \mu_j^2}  \left( \lambda \mathscr{H}(\lambda) + \frac{(\psi \varphi + \nu ^2) \lambda^2}{\xi^2 - (\psi \varphi + \nu ^2) \lambda^2} - \sum _{k \neq j}^{M+1} \frac{\lambda ^2}{\lambda ^2 - \mu_k^2} \right) 
\notag \\[1ex]
%
&+ \sum _{j=1}^{M+1} \frac{4}{\lambda ^2 - \mu _j^2} \ \frac{\lambda}{\mu _j} \ \frac{\xi - \mu _j \sqrt{\psi \varphi + \nu ^2}}{\xi - \lambda \sqrt{\psi \varphi + \nu ^2}} \ \mathscr{F} (\lambda) \mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \cdots \widehat{\mathscr{F} (\mu _j)} \cdots
\notag \\[1ex]
%
&\cdots \mathscr{F} (\mu _M)  \mathscr{F} (\mu _{M+1}) 
\left( \mu _j \mathscr{H}(\mu _j) + \frac{(\psi \varphi + \nu ^2) \mu _j^2}{\xi^2 - (\psi \varphi + \nu ^2) \mu _j^2} - \sum _{k \neq j}^{M+1} \frac{2 \mu _j^2}{\mu _j^2 - \mu_k^2} \right) \>.
\end{align}
This completes the proof by mathematical induction of the formula \eqref{Eq: tau-Mf}. It should be stressed that the right hand side is in an algebraically closed form. Thus the off-shell action of the generating function of the $so(3)$ Gaudin Hamiltonians with general boundary terms becomes a simple corollary of this result. 





Namely, from \eqref{Eq: tau-Mf} it follows that, for an arbitrary natural number $M$, the off-shel action of the generating function $\tau (\lambda)$ on the Bethe vectors 
\begin{equation}
\label{Eq: Phi-M}
\Phi _M ( \mu _{1}, \mu _{2}, \ldots , \mu _{M}) = \mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \cdots \mathscr{F} (\mu _M) \Omega _+
\> ,
\end{equation}
is given by
\begin{equation}
\label{Eq: tau-Phi}
\begin{split}
&\tau (\lambda) \Phi _M ( \mu _{1}, \mu _{2}, \ldots , \mu _{M}) = \chi _M (\lambda , \mu _{1}, \mu _{2}, \ldots , \mu _{M}) \ \Phi _M ( \mu _{1}, \mu _{2}, \ldots , \mu _{M}) 
\\[1.5ex]
%
&+ \sum _{j=1}^{M} \frac{4\lambda }{\lambda ^2 - \mu _j^2} \frac{\xi - \mu _j \sqrt{\psi \varphi + \nu ^2}}{\xi - \lambda \sqrt{\psi \varphi + \nu ^2}} 
\left( \rho (\mu _j) + \frac{(\psi \varphi + \nu ^2) \mu _j}{\xi^2 - (\psi \varphi + \nu ^2) \mu _j^2} - \sum _{k \neq j}^{M} \frac{2 \mu _j}{\mu _j^2 - \mu_k^2} \right) \times \\[1.5ex]
%
&\qquad \times \Phi _M ( \lambda, \mu _{1}, \ldots , \widehat{\mu} _{j}, \dots , \mu _{M}) 
\>,
\end{split}
\end{equation}
where the eigenvalue $\chi _M (\lambda , \mu _{1}, \mu _{2}, \ldots , \mu _{M})$  is given by
\begin{equation}
\label{Eq: eig-val}
\begin{split}
&\chi _M (\lambda , \mu _{1}, \mu _{2}, \ldots , \mu _{M}) = \\
%
&\qquad \chi_0 (\lambda) - \sum _{j=1}^{M} \frac{4\lambda}{\lambda ^2 - \mu_j^2}  \left( \rho (\lambda) + \frac{(\psi \varphi + \nu ^2) \lambda}{\xi^2 - (\psi \varphi + \nu ^2) \lambda^2} - \sum _{k \neq j}^{M} \frac{\lambda}{\lambda ^2 - \mu_k^2} \right) \> .
\end{split}
\end{equation}
The unwanted terms on the right hand side of \eqref{Eq: tau-Phi} are annihilated once the Bethe equations
\begin{equation}
\label{Eq: BE}
\rho (\mu _j) + \frac{(\psi \varphi + \nu ^2) \mu _j}{\xi^2 - (\psi \varphi + \nu ^2) \mu _j^2} - \sum _{k \neq j}^{M} \frac{2 \mu _j}{\mu _j^2 - \mu_k^2} = 0 \>, \qquad j= 1, 2, \ldots , M \>,
\end{equation}
are imposed on the parameters $\mu _{1}, \mu _{2}, \ldots , \mu _{M}$.


The off-shell action of the $so(3)$ Gaudin Hamiltonians with general boundary terms \eqref{Eq: Ham-m} on the Bethe vectors \eqref{Eq: Phi-M} is obtained by taking the residue, at $\lambda = \alpha_m$, of the equation \eqref{Eq: tau-Phi} 
\begin{equation}
\label{Eq: HamPhi}
\begin{split}
&\mathrm{H}_m \ \Phi _M ( \mu _{1}, \mu _{2}, \ldots , \mu _{M}) = \mathbf{\mathcal{E}} _{m,M} \ \Phi _M ( \mu _{1}, \mu _{2}, \ldots , \mu _{M}) \\[1.5ex]
%
&+ \sum _{j=1}^{M} \frac{4\alpha_m }{\alpha_m ^2 - \mu _j^2} \frac{\xi - \mu _j \sqrt{\psi \varphi + \nu ^2}}{\xi - \alpha_m \sqrt{\psi \varphi + \nu ^2}} \times 
\\[1.5ex]
%
&\qquad \times  
\left( \rho (\mu _j) + \frac{(\psi \varphi + \nu ^2) \mu _j}{\xi^2 - (\psi \varphi + \nu ^2) \mu _j^2} - \sum _{k \neq j}^{M} \frac{2 \mu _j}{\mu _j^2 - \mu_k^2} \right)\times 
\\[1.5ex]
%
&\qquad\times 
\left( \frac{-2 ( \nu + \sqrt{\psi \varphi + \nu^2} ) S^3_m - \psi S^{+}_m  + \frac{\psi \varphi + 2 \nu \, ( \nu + \sqrt{\psi \varphi + \nu^2})}{\psi}   S^{-} _m}{2\sqrt{\psi \varphi + \nu^2}} \right)
\times 
\\[1.5ex]
%
&\qquad\times \Phi _{M-1} ( \mu _{1}, \ldots , \widehat{\mu} _{j}, \dots , \mu _{M}) \>,
\end{split}
\end{equation}
where 
\begin{equation}
\mathrm{H}_m = \underset{\lambda = \alpha_m}{\Res} \tau (\lambda)
\end{equation}
and the eigenvalues $\mathbf{\mathcal{E}} _{m,M}$ of the $so(3)$ Gaudin Hamiltonians are the residues of the eigenvalues $\chi _M (\lambda , \mu _{1}, \mu _{2}, \ldots , \mu _{M})$ \eqref{Eq: eig-val} of the generating function $\tau (\lambda)$ at $\lambda = \alpha_m$, 
\begin{equation}
\label{Eq: Reschi}
\begin{split}
\mathbf{\mathcal{E}} _{m,M} &= \underset{\lambda = \alpha_m}{\Res} \chi _M (\lambda , \mu _{1}, \mu _{2}, \ldots , \mu _{M})  \\
%
&=\frac{2 \xi ^2}{\left( \xi ^2 - \left( \psi \varphi + \nu ^2 \right) \alpha_m^2 \right) \alpha_m} + \sum _{n \neq m}^{N} \frac{4\alpha_m}{\alpha_m^2 - \alpha_n^2}
- \sum _{j =1}^{M} \frac{4\alpha_m}{\alpha_m^2 - \mu _j^2} \>,
\end{split}
\end{equation}
and
\begin{align}
\label{Eq: ResPhiM}
&\underset{\lambda = \alpha_m}{\Res} \Phi _M ( \lambda, \mu _{1}, \ldots , \widehat{\mu} _{j}, \dots , \mu _{M}) = \underset{\lambda = \alpha_m}{\Res} \left( \mathscr{F} (\lambda) \right) \cdot \Phi _{M-1} ( \mu _{1}, \ldots , \widehat{\mu} _{j}, \dots , \mu _{M})
\notag \\[1.5ex]
%
&= \left( \frac{-2 ( \nu + \sqrt{\psi \varphi + \nu^2} ) S^3_m - \psi S^{+}_m  + \frac{\psi \varphi + 2 \nu \, ( \nu + \sqrt{\psi \varphi + \nu^2})}{\psi}   S^{-} _m}{2\sqrt{\psi \varphi + \nu^2}} \right)  \cdot 
\notag \\[1.5ex]
%
&\qquad\cdot \Phi _{M-1} ( \mu _{1}, \ldots , \widehat{\mu} _{j}, \dots , \mu _{M})
\>, 
\end{align}
where the notation $\widehat{\mu} _{j}$ means that the argument $\mu _{j}$ is omitted. 

As a closing remark for this section we must underline the complete gerality of these results: the formulae for the off-shell action of the generating function $\tau (\lambda)$ \eqref{Eq: tau-Phi} and the $so(3)$ Gaudin Hamiltonians on the Bethe vectors \eqref{Eq: Phi-M} are obtained  for an arbitrary natural number $M$ and without any restriction whatsoever on all four boundary parameters. In this sense we can say that these formulae are as general as they can possibly be. In the next section we will establish a correspondence between the  Bethe vectors \eqref{Eq: Phi-M} established here and the solutions to the generalized $so(3)$ Knizhnik-Zamolodchikov equations.



%
%
%
%
\section{Generalized $so(3)$ Knizhnik-Zamolodchikov equations \label{Sec: gKZeq} }
%
%
In this section we study solutions to the generalized $so(3)$ Knizhnik-Zamolodchikov equations. To proceed further, we now have to set the parameter $\xi$ to zero, i.e. $\xi = 0$. Consequently, the local realization of the generators \eqref{Eq: e-scr-l} -- \eqref{Eq: h-scr-l} simplifies to
\begin{align}
\label{Eq: e-scr-til}
\widetilde{\mathscr{E}} (\lambda) &=  \sum _{m=1}^N \frac{\alpha_m}{\sqrt{\psi \varphi + \nu ^2}} \  \frac{2 ( \nu - \sqrt{\psi \varphi + \nu^2} ) S^3_m + \psi S^{+}_m  - \frac{\psi \varphi + 2 \nu \, ( \nu - \sqrt{\psi \varphi + \nu^2})}{\psi}   S^{-} _m}{\lambda ^2 - \alpha_m^2} 
\>, \\[2ex]
%
%
\label{Eq: f-scr-til}
\widetilde{\mathscr{F}} (\lambda) &=  \sum _{m=1}^N \frac{- \alpha_m}{\sqrt{\psi \varphi + \nu ^2}} \ \frac{2 ( \nu + \sqrt{\psi \varphi + \nu^2} ) S^3_m + \psi S^{+}_m  - \frac{\psi \varphi + 2 \nu \, ( \nu + \sqrt{\psi \varphi + \nu^2})}{\psi}   S^{-} _m}{\lambda ^2 - \alpha_m^2}  
  \>, \\[2ex]
%
%
\label{Eq: h-scr-loc}
\mathscr{H} (\lambda) &= \frac{\lambda}{\sqrt{\psi \varphi + \nu^2}} \ \sum _{m=1}^N \frac{2\nu S^3_m + \psi S^{+}_m  + \varphi  S^{-} _m}{\lambda ^2 - \alpha_m^2} 
\>.
\end{align}
These generators have the following non-trivial commutation relations 
\begin{align}
\label{Eq: cr-HE-scr-til}
&\left[ \mathscr{H} (\lambda) , \widetilde{\mathscr{E}} (\mu)  \right] = \frac{-2\lambda}{\lambda ^2 - \mu ^2} \left(  \widetilde{\mathscr{E}} (\lambda) -  \widetilde{\mathscr{E}} (\mu) \right) \>,
\\[2ex]
%
%
\label{Eq: cr-HF-scr-til} 
&\left[  \mathscr{H} (\lambda) , \widetilde{\mathscr{F}} (\mu)  \right] = \frac{2\lambda}{\lambda ^2 - \mu ^2} \left( \widetilde{ \mathscr{F}} (\lambda) - \widetilde{\mathscr{F}} (\mu) \right)  \>,
\\[2ex]
%
%
\label{Eq: cr-EF-scr-til}
&\left[ \widetilde{\mathscr{E}} (\lambda) , \widetilde{\mathscr{F}} (\mu)  \right] = \frac{-4}{\lambda ^2 - \mu ^2} \left( \lambda \, \mathscr{H} (\lambda) - \mu \, \mathscr{H} (\mu) \right)  \>.
\end{align} 
Thus, in this case, the first Bethe vector is here defined by
\begin{equation}
\label{Eq: BV-Phi-1-til}
\widetilde{\Phi} _1 ( \mu ) = \Phi _1 ( \mu ) \Big{|}_{\xi = 0} = \widetilde{\mathscr{F}} (\mu) \Omega_{+} \>, 
\end{equation}
and, in general case, for the Bethe vectors $\widetilde{\Phi} _M ( \mu _{1}, \mu _{2}, \ldots , \mu _{M}) $ we have
\begin{equation}
\label{Eq: Phi-M-til}
\widetilde{\Phi} _M ( \mu _{1}, \mu _{2}, \ldots , \mu _{M}) = \Phi _M ( \mu _{1}, \mu _{2}, \ldots , \mu _{M}) \Big{|}_{\xi = 0} 
= \widetilde{\mathscr{F}} (\mu _1) \widetilde{\mathscr{F}} (\mu _2) \cdots \widetilde{\mathscr{F}} (\mu _M) \Omega _+
\> .
\end{equation}


It is useful to explicitly write $\mathrm{H}_m$ for $\xi = 0$:
\begin{align}
\label{Eq: Ham-m-til}
&\widetilde{\mathrm{H}}_m = \mathrm{H}_m  \Big{|}_{\xi = 0} = \frac{1}{2 \left( \psi \varphi + \nu ^2 \right) \alpha_m}
\Big( 2 \left( \nu ^2  -  \psi \varphi \right) (S^3_m ) ^2 
\notag \\[1.5ex]
%
&+ \psi ^2  \ (S_m^{+})^2 + \varphi ^2 \ (S_m^{-})^2 + 2 \psi \nu \left( S_m^+ S_m^3 +  S_m^3 S_m^+ \right) 
\notag \\[1.5ex]
%
& + 2 \varphi  \nu \left( S_m^- S_m^3 + S_m^3 S_m^- \right) - \nu^2  \left( S_m^{+} S_m^{-} + S_m^{-} S_m^{+} \right) \Big) 
\notag \\[1.5ex]
%
& + \frac{1}{\psi \varphi + \nu ^2 } \ \sum _{n \neq m}^{N} \bigg( \frac{4 \left(\psi \varphi \, \alpha_n + \nu ^2 \alpha_m \right)}{\alpha_m^2 -\alpha_n^2}  \ S_m^3 S_n^3 + \frac{1}{\alpha_m + \alpha_n} \times 
\notag \\[1.5ex]
%
&\quad\times \left( \psi ^2  \ S_m^+ S_n^+ + \varphi ^2 \ S_m^- S_n^- + 2 \psi \nu \left( S_m^+ S_n^3 +  S_m^3 S_n^+ \right) + 2 \varphi  \nu \left( S_m^- S_n^3 + S_m^3 S_n^- \right) \right)
\notag \\[1.5ex]
%
&+ \frac{2 \nu ^2 \alpha_n  + \psi \varphi \, (\alpha_m +\alpha_n)}{\alpha_m^2 -\alpha_n^2}
 \left( S_m^- S_n^+ + S_m^+ S_n^- \right) \bigg) 
\>.
\end{align}
Therefore the  off-shell action of these Hamiltonians reads 
\begin{equation}
\label{Eq: HamPhi-till}
\begin{split}
&\widetilde{\mathrm{H}}_m \, \widetilde{\Phi} _M ( \mu _{1}, \mu _{2}, \ldots , \mu _{M}) = \\
%
&\qquad \widetilde{\mathbf{\mathcal{E}}} _{m,M} \ \widetilde{\Phi} _M ( \mu _{1}, \mu _{2}, \ldots , \mu _{M}) + \sum _{j=1}^{M} \frac{(-2) \mu _j }{\alpha_m ^2- \mu _j^2} \ \beta _{M} (\mu _j) \ \widetilde{\Phi} _{M-1}^{(j,m)} \>,
\end{split}
\end{equation}
where
\begin{align}
\label{Eq: EmM-til}
\widetilde{\mathbf{\mathcal{E}}} _{m,M} &= \mathbf{\mathcal{E}} _{m,M} \Big{|}_{\xi = 0} = 
\sum _{n \neq m}^{N} \frac{4\alpha_m}{\alpha_m^2 - \alpha_n^2}
- \sum _{j =1}^{M} \frac{4\alpha_m}{\alpha_m^2 - \mu _j^2} 
\>,
\\[2ex]
%
\label{Eq: beta}
\beta _{M} (\mu _j) &= -2 \left( \rho (\mu _j) - \frac{1}{\mu _j} - \sum _{k \neq j}^{M} \frac{2 \mu _j}{\mu _j^2 - \mu_k^2} \right) \> ,
\end{align}
and 
\begin{equation}
\begin{split}
\label{Eq: Phi-jm-til}
&\widetilde{\Phi} _{M-1}^{(j,m)} = \left( \frac{-2 ( \nu + \sqrt{\psi \varphi + \nu^2} ) S^3_m - \psi S^{+}_m  + \frac{\psi \varphi + 2 \nu \, ( \nu + \sqrt{\psi \varphi + \nu^2})}{\psi}   S^{-} _m}{2\sqrt{\psi \varphi + \nu^2}} \right)\cdot  \\
%
&\qquad \cdot \widetilde{\Phi} _{M-1} ( \mu _{1}, \ldots , \widehat{\mu} _{j}, \dots , \mu _{M}) \> .
\end{split}
\end{equation}





Our main objective in this section is to show how to each Bethe vector \eqref{Eq: Phi-M-til} we can relate a solution to the generalized $so(3)$ Knizhnik-Zamolodchikov equations  
\begin{equation}
\label{Eq: KZe}
\kappa \,  \partial _{\alpha_m} \Psi ( \alpha _1 , \alpha _2 , \ldots , \alpha _N ) = \widetilde{\mathrm{H}}_m \ \Psi ( \alpha _1 , \alpha _2 , \ldots , \alpha _N ) \> . 
\end{equation}
Within our approach  \cite{BabujianFlume, Hikami95, Kulish:2001hr, Kulish:2001vs, Kulish:2003da, Salom:2019myz, Salom:2021gox} this correspondence is defined by a closed contour integration with respect to the variables $\mu _1, \mu _2, \ldots , \mu _M$
\begin{equation}
\label{Eq: Sol-KZe}
\Psi ( \alpha _1 , \alpha _2 , \ldots , \alpha _N ) = \oint \oint \cdots \oint \  \Upsilon \left( \overrightarrow{\mu} ; \overrightarrow{\alpha} \right) 
\cdot \widetilde{ \Phi} _M \left( \overrightarrow{\mu} ; \overrightarrow{\alpha} \right) \, \mathrm{d} \mu _1 \, \mathrm{d} \mu _2 \cdots \mathrm{d} \mu _M \> . 
\end{equation}
The scalar function $\Upsilon \left( \overrightarrow{\mu} ; \overrightarrow{\alpha}\right) $ is defined by 
\begin{equation}
\label{Eq: Psi-f}
\Upsilon \left( \overrightarrow{\mu} ; \overrightarrow{\alpha} \right) = \exp \left( \frac{S ( \overrightarrow{\mu} ; \overrightarrow{\alpha} )}{\kappa} \right) \>,
\end{equation}
with the constant $\kappa$ and the function $S ( \overrightarrow{\mu} ; \overrightarrow{\alpha} )$ specified by
\begin{align}
\label{Eq: S-fun}
S \left( \overrightarrow{\mu} ; \overrightarrow{\alpha} \right)  =&  \sum _{m=1}^N \left( \sum_{n \neq m}^N  \ln \left(\alpha _n^2 - \alpha _m^2 \right) 
- \sum _{j=1}^M 2 \, \ln \left( \mu _j^2 - \alpha _m^2  \right) \right) \\
%
&+ \sum _{j=1}^M  \left(  \ln \left( \mu _j^2 \right) +  \sum_{k \neq j}^M \ln \left( \mu _j^2  -  \mu _k^2 \right) \right) \>.
\end{align}
It is straightforward to check that the function $\Upsilon \left( \overrightarrow{\mu} ; \overrightarrow{\alpha}\right) $ satisfies the system 
\begin{align}
\label{Eq: sys-1}
\kappa \,  \partial _{\alpha_m}  \Upsilon  &= \widetilde{\mathbf{\mathcal{E}}} _{m, M}  \, \Upsilon  \>, \\[2ex]
%
\label{Eq: sys-2}   
\kappa \,  \partial _{\mu _j}  \Upsilon  &=  \beta _{M} (\mu _j) \, \Upsilon  \>,
\end{align}
where $\widetilde{\mathbf{\mathcal{E}}} _{m, M}$  and $  \beta _{M} (\mu _j)$ are defined in \eqref{Eq: EmM-til} and \eqref{Eq: beta}, respectively.


The crucial identity in our approach is
\begin{equation}
\label{Eq: KZ-Identity}
\partial _{\alpha _m} \widetilde{\Phi }_M = \sum _{j=1}^{M} \partial _{\mu_j} \left( \frac{2 \mu _j  }{\alpha_m ^2- \mu _j^2} \ \widetilde{\Phi} _{M-1}^{(j,m)} \right) \> .
\end{equation}
It takes a few rather simple steps to confirm that the function $\Psi ( \alpha _1 , \alpha _2 , \ldots , \alpha _N )$ \eqref{Eq: Sol-KZe} is a solution to the to the generalized $so(3)$ Knizhnik-Zamolodchikov equations \eqref{Eq: KZe}. As the first step, using the Leibniz rule, we calculate 
\begin{equation}
\label{Eq: Proof-1}
\kappa \, \partial _{\alpha_m} \left( \Upsilon \cdot \widetilde{\Phi }_M \right) = \left( \kappa \,  \partial _{\alpha_m} \Upsilon \right) \cdot \widetilde{\Phi} _M + \Upsilon \cdot \left( \kappa \,  \partial _{\alpha_m} \widetilde{\Phi} _M \right) \> .
\end{equation}
Then we use the equation \eqref{Eq: sys-1} in the first term on the right hand side of the equation above and the identity \eqref{Eq: KZ-Identity} in the second term 
\begin{equation}
\label{Eq: Proof-2}
\kappa \, \partial _{\alpha_m} \left( \Upsilon \cdot \widetilde{\Phi} _M \right) = \mathbf{\mathcal{E}} _{m, M}  \left(  \Upsilon \cdot \widetilde{\Phi }_M \right) + \Upsilon \cdot \kappa \ \sum _{j=1}^{M} \partial _{\mu_j} \left( \frac{2 \mu _j}{\alpha_m ^2- \mu _j^2} \ \widetilde{\Phi} _{M-1}^{(j,m)} \right) \> .
\end{equation}
In the following step we use the equation \eqref{Eq: HamPhi-till} in the first term and the Leibniz rule in the second term on the right hand side of the equation above
\begin{equation}
\label{Eq: Proof-3}
\begin{split}
&\kappa \, \partial _{\alpha_m} \left( \Upsilon \cdot \widetilde{\Phi} _M \right) = \widetilde{\mathrm{H}}_m  \left( \Upsilon \cdot \widetilde{\Phi} _M \right) + \sum _{j=1}^{M} \frac{2 \mu _j }{\alpha_m ^2- \mu _j^2} \ \beta _{M} (\mu _j) \cdot \Upsilon \cdot \widetilde{\Phi} _{M-1}^{(j,m)} 
\\[1ex]
%
&+  \kappa \ \sum _{j=1}^{M} \partial _{\mu_j} \left( \frac{2 \mu _j }{\alpha_m ^2- \mu _j^2} \  \Upsilon \cdot \widetilde{\Phi} _{M-1}^{(j,m)} \right) 
- \kappa \ \sum _{j=1}^{M} \left( \partial _{\mu_j} \Upsilon \right) \frac{2 \mu _j }{\alpha_m ^2- \mu _j^2} \ \widetilde{\Phi }_{M-1}^{(j,m)} \>.
\end{split}
\end{equation}
Now it remains to rewrite the last terms on the right hand side using the equation \eqref{Eq: sys-2} 
\begin{equation}
\label{Eq: Proof-4}
\begin{split}
&\kappa \, \partial _{\alpha_m} \left( \Upsilon \cdot \widetilde{\Phi }_M \right) = \widetilde{\mathrm{H}}_m  \left( \Upsilon \cdot \widetilde{\Phi} _M \right) + \sum _{j=1}^{M} \frac{2 \mu _j}{\alpha_m ^2- \mu _j^2} \ \beta _{M} (\mu _j) \cdot \Upsilon \cdot \widetilde{\Phi} _{M-1}^{(j,m)} 
\\[1ex]
%
&+  \kappa \ \sum _{j=1}^{M} \partial _{\mu_j} \left( \frac{2 \mu _j }{\alpha_m ^2- \mu _j^2} \  \Upsilon \cdot \widetilde{\Phi} _{M-1}^{(j,m)} \right) 
- \sum _{j=1}^{M} \frac{2 \mu _j }{\alpha_m ^2- \mu _j^2}  \ \beta _{M} (\mu _j) \cdot \Upsilon \cdot \widetilde{\Phi} _{M-1}^{(j,m)} \>.
\end{split}
\end{equation}
As the final step of this demonstration, we simplify the second and the last term in the equation above in order to obtain the desired result
\begin{equation}
\label{Eq: Proof-5}
\kappa \, \partial _{\alpha_m} \left( \Upsilon \cdot \widetilde{\Phi} _M \right) = \widetilde{\mathrm{H}}_m  \left( \Upsilon \cdot \widetilde{\Phi} _M \right) +  \kappa \ \sum _{j=1}^{M} \partial _{\mu_j} \left( \frac{2 \mu _j }{\alpha_m ^2- \mu _j^2} \  \Upsilon \cdot \widetilde{\Phi} _{M-1}^{(j,m)} \right) \>.
\end{equation}
This shows that the function $\Psi ( \alpha _1 , \alpha _2 , \ldots , \alpha _N )$ \eqref{Eq: Sol-KZe} is a solution to the to the generalized $so(3)$ Knizhnik-Zamolodchikov equations \eqref{Eq: KZe} since the terms in the sum will not contribute to the closed contour integrals with respect to the variables $\mu _j$, $j=1, 2, \ldots , M$.




In the final part of this section we determine the on-shell norm as well as the off-shell scalar products of the Bethe vectors \eqref{Eq: Phi-M-til}. In particular, the on-shell norm of the Bethe vector \eqref{Eq: BV-Phi-1-til} is obtained to be
\begin{equation}
\label{Eq: norm-BV-Phi-til-1}
\begin{split}
\big{\|} \widetilde{\Phi} _1 (\mu ) \big{\|}^2 &= \lim _{\nu \to \mu} \big{\langle}  \Omega _+ , \widetilde{\mathscr{E}} (\nu) \widetilde{\mathscr{F}} (\mu) \Omega _+  \big{\rangle}  = -2 \left( \rho ^{\prime}(\mu) + \frac{\rho (\mu) }{\mu}\right) \\[2ex]
&= \frac{\partial \beta _1 (\mu)}{\partial \mu} \Big{|}_{\beta _1 (\mu) = 0} = \frac{\partial ^2 S (\mu)}{\partial \mu ^2} \Big{|}_{\beta _1 (\mu) = 0}
\>.
\end{split}
\end{equation}
Similarly, the norm of the Bethe vector
\begin{equation}
\label{Eq: BV-Phi-til-2}
\widetilde{\Phi} _2 ( \mu _1 , \mu _2 ) = \widetilde{\mathscr{F}} (\mu _1)  \widetilde{\mathscr{F}} (\mu _2) \Omega_{+} \>, 
\end{equation}
when the Bethe equations are imposed on the parameters $\mu _1$ and $\mu _2$, is given by
\begin{equation}
\label{Eq: norm-BV-Phi-til-2}
\begin{split}
&\big{\|} \widetilde{\Phi} _2 (\mu _1 , \mu_2 ) \big{\|}^2 = \underset{{\nu_2 \to \mu_2}}{\lim _{\nu_1 \to \mu_1}}  \big{\langle} \Omega _+ ,  \ \widetilde{\mathscr{E}} (\nu_1) \widetilde{\mathscr{E}} (\nu_2)  \widetilde{\mathscr{F}} (\mu_2) \widetilde{\mathscr{F}} (\mu_1) \Omega _+  \big{\rangle}  
\\[1.5ex]
%
&= 4  \ \rho ^{\prime} (\mu_1) \, \rho ^{\prime} (\mu_2) + 4 \ \rho ^{\prime} (\mu_1) \left( \frac{1}{\mu _2^2} + \frac{2}{\mu _2^2 - \mu _1^2} + \frac{4 \, \mu _1^2}{\left( \mu _2^2 - \mu _1^2 \right) ^2} \right)\\[1.5ex]
%
&\qquad + 4 \ \rho ^{\prime} (\mu_2) \left( \frac{1}{\mu _1^2} + \frac{2}{\mu _1^2 - \mu _2^2} + \frac{4 \, \mu _2^2}{\left( \mu _2^2 - \mu _1^2 \right) ^2} \right) + 12 \ \frac{\left( \mu_1^2 + \mu_2^2 \right)  ^2}{\mu _1^2 \mu _2^2 \left( \mu _2^2 - \mu _1^2 \right) ^2} \\[2ex]
%
%
&= \det \left(\begin{array}{cc}  
\frac{\partial \beta _2 (\mu_1)}{\partial \mu_1} & \frac{\partial \beta _2 (\mu_2)}{\partial \mu_1}  \\[1.5ex]
\frac{\partial \beta _2 (\mu_1)}{\partial \mu_2} & \frac{\partial \beta _2 (\mu_2)}{\partial \mu_2} 
 \end{array} \right) \Bigg{|}_{\genfrac{}{}{0pt}{2}{\beta _2 (\mu_1) = 0}{\beta _2 (\mu_2) = 0}} 
 %
 = \det \left(\begin{array}{cc}  
\frac{\partial ^2 S}{\partial \mu_1 ^2} & \frac{\partial ^2 S}{\partial \mu_1 \, \partial \mu_2}  \\[1.5ex]
\frac{\partial ^2 S}{\partial \mu_2 \, \partial \mu_1} & \frac{\partial ^2 S}{\partial \mu_2^2 } 
 \end{array} \right) \Bigg{|}_{\genfrac{}{}{0pt}{2}{\beta _2 (\mu_1) = 0}{\beta _2 (\mu_2) = 0}}.
\end{split}
\end{equation}
In the general case, for an arbitrary positive integer $M$, the norm of the Bethe vector\break $\widetilde{\Phi} _M ( \mu _{1}, \mu _{2}, \ldots , \mu _{M})$ \eqref{Eq: Phi-M-til}, when the Bethe equations 
\begin{equation}
\label{Eq: beta-M}
\beta _{M} (\mu _j) = -2 \left( \rho (\mu _j) - \frac{1}{\mu _j} - \sum _{k \neq j}^{M} \frac{2 \mu _j}{\mu _j^2 - \mu_k^2} \right) = 0
\> , \quad j= 1, 2, \ldots , M \> ,
\end{equation}
are imposed on the parameter $\mu _1, \ldots , \mu _M$, is obtained to be
\begin{equation}
\label{Eq: norm-BV-Phi-til-M}
\big{\|} \widetilde{\Phi} _M (\mu _1 , \mu _2 ,  \ldots ,  \mu _M ) \big{\|}^2 = \det \left(\begin{array}{cccc}  
\frac{\partial^2 S}{\partial \mu_1 ^2} & \frac{\partial^2 S}{\partial \mu_1 \partial \mu_2} & \hdots &  \frac{\partial^2 S}{\partial \mu_1 \partial \mu_M} \\[1ex]
\vdots & \ddots & &\vdots  \\[1ex]
 \frac{\partial^2 S}{\partial \mu_M \partial \mu_1} & \frac{\partial^2 S}{\partial \mu_M \partial \mu_2} & \hdots & \frac{\partial^2 S}{\partial \mu_M ^2}
 \end{array} \right) \vline_{\genfrac{}{}{0pt}{2}{\beta _M (\mu_1) = 0}{\genfrac{}{}{0pt}{2}{\vdots}{\beta _M (\mu_M) = 0}}} \> .
\end{equation}



Finally, we also calculate the off-shell scalar products of the Bethe vectors $\widetilde{\Phi} _M  (\mu _1 , \mu _2 ,  \ldots ,  \mu _M)$ \eqref{Eq: Phi-M-til}. As our first step, we observe that in the case when $M=1$ the scalar product is 
\begin{equation}
\label{Eq: ScPr--BV-Phi-til-1}
\big{\langle}  \widetilde{\Phi}  _1 (\mu ) , \widetilde{\Phi}  _1 (\nu ) \big{\rangle}  
= 4 \left( - \frac{\mu \, \rho (\mu) - \nu \, \rho (\nu)}{\mu ^2 - \nu ^2} \right) \>.
\end{equation}
For $M=2$, a straightforward calculation yields
\begin{equation}
\label{Eq: ScPr--BV-Phi-til-2}
\big{\langle}  \widetilde{\Phi} _2 (\mu _1 , \mu_2  ) , \widetilde{\Phi} _2 ( \nu _1 ,  \nu _2 ) \big{\rangle} = 4^2 \sum _{\sigma \in \mathscr{S}_2} \det \mathscr{M} ^{\sigma} = 16 \left( \det \mathscr{M} ^{1} + \det \mathscr{M} ^{2} \right) \>,
\end{equation}
where $\mathscr{S}_2$ is the symmetric group of degree two and the two-by-two matrices $\mathscr{M} ^{1}$ and $\mathscr{M} ^{2}$ are given by
\begin{equation}
\label{Eq: M1}
\begin{aligned}
\mathscr{M} ^{1}_{11}  &= - \frac{ \mu _1 \, \rho (\mu _1) - \nu _1 \, \rho (\nu _1)}{\mu _1^2 - \nu _1^2}  
- \frac{ \mu _2^2 + \nu _2 ^2}{\left( \mu _1^2 - \mu _2^2 \right) \left( \nu _1^2 - \nu _2^2 \right)} \>,
%
\\[2ex]
%
\mathscr{M} ^{1}_{12}  &= - \frac{\mu _2^2 +  \nu _2^2}{\left(\mu _1^2 - \mu _2^2 \right) \left( \nu _1^2 - \nu _2^2 \right)} \>,
\\[2ex]
%
\mathscr{M} ^{1}_{22} &=  - \frac{\mu _2 \, \rho (\mu _2) - \nu _2 \, \rho (\nu _2)}{\mu _2^2 - \nu _2^2} 
- \frac{\mu _1^2 + \nu _1^2}{\left( \mu _2 ^2 - \mu _1 ^2 \right) \left( \nu _2 ^2 - \nu _1^2 \right)} \>,
%
\\[2ex]
%
\mathscr{M} ^{1}_{21}  &= - \frac{\mu _1^2 + \nu _1^2}{\left( \mu _2^2 - \mu _1^2 \right) \left( \nu _2^2 - \nu _1^2 \right)} \> ,
\end{aligned}
\end{equation}
and 
\begin{equation}
\label{Eq: M2}
\begin{aligned}
\mathscr{M} ^{2}_{11}  &= - \frac{\mu _1 \, \rho( \mu _1) -  \nu _2 \,  \rho (\nu _2)}{\mu _1^2 - \nu _2^2}  
- \frac{\mu _2^2 + \nu _1^2}{\left( \mu _1^2 - \mu _2^2 \right) \left( \nu _2 ^2 - \nu _1^2 \right)} \>,
%
\\[2ex]
%
 \mathscr{M} ^{2}_{12}  &= - \frac{\mu _2^2  +  \nu _1^2}{\left( \mu _1^2 - \mu _2^2 \right) \left( \nu _2 ^2 - \nu _1^2 \right)}  \>,
\\[2ex]
%
\mathscr{M} ^{2}_{22}  &=  - \frac{\mu _2 \, \rho (\mu _2) - \nu _1 \, \rho (\nu _1)}{\mu _2 ^2 - \nu _1^2} 
- \frac{\mu _1^2 + \nu _2^2}{\left( \mu _2^2 - \mu _1^2 \right) \left( \nu _1 ^2 - \nu _2^2 \right)} \>,
%
\\[2ex]
%
\mathscr{M} ^{2}_{21}  &= -\frac{\mu _1^2  + \nu _2^2}{\left( \mu _2 ^2 - \mu _1^2 \right) \left( \nu _1^2  - \nu _2^2 \right)} \> .
\end{aligned}
\end{equation}
In general case, for an arbitrary positive integer $M$, we have
\begin{equation}
\label{Eq: ScPr--BV-Phi-til-M}
\big{\langle} \widetilde{\Phi} _M (\mu _1 , \mu_2 , \ldots , \mu_M ) , \ \widetilde{\Phi} _M ( \nu _1 , \nu_2 , \ldots , \nu_M ) \big{\rangle}
= 4^M \sum _{\sigma \in \mathscr{S}_M} \det \mathscr{M} ^{\sigma} \>,
\end{equation}
where $S_M$ is the symmetric group of degree M and the matrix entries of the $M\times M$ matrix $\mathcal{M}^{\sigma}$ are given by
\begin{align}
\label{Eq: com-M}
\mathscr{M}^{\sigma}_{jj}  &= - \frac{\mu _j \, \rho (\mu_j) - \nu _{\sigma(j)} \, \rho (\nu _{\sigma(j)}) }{\mu_j^2 - \nu _{\sigma(j)}^2} - \sum _{k\neq j} \frac{\mu _k^2 +\nu_{\sigma(k)}^2}{(\mu_j^2 - \mu _k^2)(\nu_{\sigma(j)}^2 - \nu _{\sigma(k)}^2)} \> ,\\[2ex]
%
\mathscr{M}^{\sigma}_{jk} &=  - \frac{\mu _k^2 +\nu_{\sigma(k)}^2}{(\mu _j^2 - \mu _k^2)(\nu _{\sigma(j)}^2 - \nu _{\sigma(k)}^2)} \ , \quad \text{for} \quad j , k = 1, 2, \ldots , M \>.
\end{align}

The off-shell scalar products of the Bethe vectors $\widetilde{\Phi} _M  (\mu _1 , \mu _2 ,  \ldots ,  \mu _M)$ \eqref{Eq: Phi-M-til} in the $M=1$ case \eqref{Eq: ScPr--BV-Phi-til-1} and in the $M=2$ case \eqref{Eq: ScPr--BV-Phi-til-2} were derived by a direct, straightforward calculations.  The formula \eqref{Eq: ScPr--BV-Phi-til-M} was obtained by symbolic computer calculations for $M = 3, 4, 5$, for some values of $N$. In the general case, the proof of these formulae by induction would be very difficult, since it would require some highly non-trivial relations between a certain type of determinants of a different order. In this sense, the general formula \eqref{Eq: ScPr--BV-Phi-til-M}, strictly speaking, remains a conjecture.




%
%
%
%
\section{Conclusions \label{Sec: Con}}
%
%
The cornerstone of our study of the non-periodic $so(3)$ Gaudin model was  the $so(3)$ Maillet linear bracket \eqref{Eq: M-br-rk} for the suitable Lax operator \eqref{Eq: cal-L} and the non-unitary $so(3)$ classical r-matrix \eqref{Eq: rK} constructed from the generic boundary K-matrix \eqref{Eq: K-mat}. Based on Maillet bracket we obtained the generating function \eqref{Eq: tau} of the $so(3)$ Gaudin Hamiltonians with general boundary terms. However it turned out that the natural set of generators was not the most efficient choice for implementing the algebraic Bethe ansatz due to the cumbersome commutation relations \eqref{Eq: cr-EE} - \eqref{Eq: cr-FE}. For this reason we proposed a new set of generators: \eqref{Eq: e-scr} - \eqref{Eq: h-scr}. Not only that their commutation relations had a strikingly compact form \eqref{Eq: cr-zero} - \eqref{Eq: cr-EF-scr} but also their local realization was fairly simple \eqref{Eq: e-scr-l} - \eqref{Eq: h-scr-l}, yielding the explicit expression for the Gaudin Hamiltonians with all four boundary parameters \eqref{Eq: Ham-m}. 

There were several preceding objectives which we had to address before attempting to find the off-shell action of the generating function $\tau(\lambda)$ \eqref{Eq: tau-gen}. In the first place, we had to define 
the so-called vacuum vector $\Omega _+$ \eqref{Eq: om} - \eqref{Eq: Omega+}. Then we had to confirm the action of the generators on the vacuum vectors \eqref{Eq: Vac} and to show that the vector $\Omega _+$ is the eigenvector of the generating function $\tau(\lambda)$ \eqref{Eq: tau-Om}. 
As our next step, we have calculated the commutation relations between the generating function $\tau(\lambda)$ and the remaining generator $\mathscr{F} (\mu)$ of the generalized $so(3)$ Gaudin algebra \eqref{Eq: tau-f}. The idea of using the generator $\mathscr{F} (\mu)$ as the so-called creation operator prompted us to calculate the commutation relations between the generating function $\tau(\lambda)$ and the product $\mathscr{F} (\mu _1) \mathscr{F} (\mu _2)$ \eqref{Eq: tau-ff} as well as the product $\mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \mathscr{F} (\mu _3)$ \eqref{Eq: tau-fff}. Hence, we have conjectured the formula \eqref{Eq: tau-Mf},  in the general case, for the commutator between  the generating function $\tau(\lambda)$ and the product $\mathscr{F} (\mu _1) \mathscr{F} (\mu _2) \cdots \mathscr{F} (\mu _M) $, for an arbitrary natural number $M$. The proof of the formula \eqref{Eq: tau-Mf} based on the mathematical induction was presented in \eqref{Eq: tau-M+1f-s1} - \eqref{Eq: tau-M+1f-s6}.
Once we have accordingly defined  the Bethe vectors \eqref{Eq: Phi-M}, the off-shell action \eqref{Eq: tau-Phi} of the generating function $\tau(\lambda)$,  including the formulae for the eigenvalues \eqref{Eq: eig-val} and the Bethe equations \eqref{Eq: BE}, followed from \eqref{Eq: tau-Mf}. 
Moreover, the off-shell action \eqref{Eq: HamPhi} of the Gaudin Hamiltonians \eqref{Eq: Ham-m} on the Bethe vectors \eqref{Eq: Phi-M} was obtained by taking the the residue, at $\lambda = \alpha_m$, of the left and the right hand side of \eqref{Eq: tau-Phi}.  It should be stressed that the formulae of the off-shell action \eqref{Eq: tau-Phi} and \eqref{Eq: HamPhi} have been obtained  without any restriction whatsoever on any of the four boundary parameters and therefore  we can say that these formulae are as general as they can possibly be. 

Next, we found the solutions to the generalized $so(3)$ Knizhnik-Zamolod\-chikov equations \eqref{Eq: KZe}. In spite that the key identity \eqref{Eq: KZ-Identity} in the proof required the parameter $\xi$ to be set to zero, the formulae we obtained for the solutions to the generalized $so(3)$ Knizhnik-Zamolodchikov equations \eqref{Eq: Sol-KZe}, the on-shell norm of the Bethe vectors \eqref{Eq: norm-BV-Phi-til-M} and the off-shell scalar product of the Bethe vectors \eqref{Eq: ScPr--BV-Phi-til-M} -- all possess higher degree of generality than the analogous formulae in the $s\ell(2)$ case \cite{Salom:2019myz} (here we have fixed only one of the four boundary parameters instead of all four).

In our future research we hope to address the remaining open problem of correlation functions for the $so(3)$ Gaudin model with general boundary, following Sklyanin approach in the periodic $s\ell (2)$ case \cite{Sklyanin:1997zz}.


\bigskip


\noindent
\textbf{Acknowledgments}

\medskip

\noindent
The authors thank to Rodrigo Pimenta for helpful discussions. IS is supported by the Ministry of Education, Science and Technological Development (MPNTR) of the Republic of Serbia, and by the Science Fund of the Republic of Serbia, Program DIASPORA, No. 6427195, SQ2020.



%
% 
%
%
\appendix
%
%


%
%
%
%
\section{Preliminaries \label{App: Prelim1}}
%
%
Some essential definitions regarding the $so(3)$ Lie algebra and its fundamental representation are given in the Appendix \ref{App: Prelim1}. Namely, 
we consider the spin one operators 
\begin{equation}
\label{Eq: S-mats}
\mathscr{S}^x = \frac{1}{\sqrt{2}}\left(\begin{array}{ccc}0 & 1 & 0 \\1 & 0 & 1 \\0 & 1 & 0\end{array}\right), \
\mathscr{S}^y=\frac{\imath}{\sqrt{2}}\left(\begin{array}{ccc}0 & -1 & 0 \\1 & 0 & -1 \\0 & 1 & 0\end{array}\right), \
\mathscr{S}^z=\left(\begin{array}{ccc}1 & 0 & 0 \\0 & 0 & 0 \\0 & 0 & -1\end{array}\right)
\end{equation}
acting in the space $V^{(1)}=\mathbb{C}^3$ with the commutation relations
$$
[ \mathscr{S}^x , \mathscr{S}^y ]=\imath \mathscr{S}^z, \quad [ \mathscr{S}^z , \mathscr{S}^x ]=\imath \mathscr{S}^y, \quad [ \mathscr{S}^y , \mathscr{S}^z ]=\imath S^x
$$
and the Casimir element
\begin{equation}
\label{C2}
    c_2 = \vec{\mathscr{S}} \cdot \vec{\mathscr{S}}  = (\mathscr{S}^x) ^2 + (\mathscr{S}^y) ^2 + (\mathscr{S}^z) ^2 = 2 \mathbbm{1}.
\end{equation}
Introducing raising and lowering operators 
\begin{equation}
\label{Eq: S-pm-mats}
\mathscr{S}^{+} = \mathscr{S}^x + \imath \mathscr{S}^y = 
\left(\begin{array}{ccc}0 & \sqrt{2} & 0 \\0 & 0 & \sqrt{2} \\0 & 0 & 0\end{array}\right), \quad
\mathscr{S}^{-}= \mathscr{S}^x - \imath \mathscr{S}^y =
\left(\begin{array}{ccc}0 & 0 & 0 \\\sqrt{2} & 0 & 0 \\0 & \sqrt{2} & 0\end{array}\right),
\end{equation}
the relations above can also be written as
\begin{equation}
\label{crspin1}
[ \mathscr{S}^z , \mathscr{S}^{\pm}] = \pm \mathscr{S}^{\pm}, \quad [ \mathscr{S}^+ , \mathscr{S}^-] = 2 \mathscr{S}^z , 
\end{equation}
and
\begin{equation}
\label{C2-a}
c_2 = (\mathscr{S}^z) ^2 + \frac{1}{2} (\mathscr{S}^+ \mathscr{S}^- + \mathscr{S}^- \mathscr{S}^+) = (\mathscr{S}^z) ^2 + \mathscr{S}^z + \mathscr{S}^- \mathscr{S}^+.
\end{equation}



 It is useful to notice that the tensor Casimir operator can be expressed as follows
\begin{align}
\label{Eq: c2tensor-12}
    c_2^{\otimes}(1,2)    &=  \vec{\mathscr{S}}_{1} \cdot \vec{\mathscr{S}}_{2}  = \mathcal{P} - 3 \mathcal{K} .
\end{align}
The permutation operator
\begin{equation}
\label{Pmat}
 \mathcal{P}
    =  \left(\begin{array}{rrr|rrr|rrr}1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ \hline
  0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ \hline
  0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1\end{array}\right),
\end{equation}
the rank 1 projector  
\begin{equation}
\label{Ktensor}
 \mathcal{K}  =  \frac{1}{3} \left(\begin{array}{rrr|rrr|rrr}
  0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 1 & 0 & -1 & 0 & 1 & 0 & 0 \\ \hline
  0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & -1 & 0 & 1 & 0 & -1 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \hline
  0 & 0 & 1 & 0 & -1 & 0 & 1 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\end{array}\right) 
\end{equation}
and the identity operator $\mathbbm{1}$ satisfy the relations
\begin{equation}
\label{Br}
\mathcal{P}^2 =\mathbbm{1} \>, \quad   \mathcal{K}^2 =  \mathcal{K} \>, \quad  \mathcal{P} \mathcal{K}= \mathcal{K} \mathcal{P}=  \mathcal{K} \>.
\end{equation}
and therefore define the representation of the Brauer algebra in  $\mathbb{C}^3\otimes \mathbb{C}^3$. Moreover these are the three invariant operators acting on $\mathbb{C}^3\otimes \mathbb{C}^3$ 
\begin{equation}
\label{inv}
[ \mathscr{S}^{\alpha} \otimes \mathbbm{1} + \mathbbm{1} \otimes \mathscr{S}^{\alpha} ,\mathcal{P}] = 0 \>, \
[ \mathscr{S}^{\alpha} \otimes \mathbbm{1} + \mathbbm{1} \otimes \mathscr{S}^{\alpha} ,\mathcal{K}] = 0 \>, 
\end{equation}
here $\alpha = x,y,z.$ 



In our study of the $so(3)$ Gaudin model with $N$ sites, characterised by the local space $V_ m = \mathbb{C}^{3}$ together with the corresponding inhomogeneous parameter $\alpha _m$, the Hilbert space is given by 
\begin{equation}
\label{Eq: H-space}
\mathcal{H} = \underset {m=1}{\overset {N}{\otimes}}  V_m = (\mathbb{C}^{3} ) ^{\otimes N} \>.
\end{equation}
The local spin operators 
\begin{equation}
\label{Eq: S-local}
S_m^{\alpha} = \mathbbm{1} \otimes \cdots \otimes \underbrace{S^{\alpha}} _m \otimes \cdots \otimes \mathbbm{1} \>,
\end{equation}
with $\alpha = x, y, z$ and $m = 1, 2 ,\dots , N$, are given by the matrices \eqref{Eq: S-pm-mats} and \eqref{Eq: S-mats} in every local space $V_ m = \mathbb{C}^{3}$. Evidently, they satisfy the usual commutation relations
\begin{equation}
\label{Eq: crsl2}
[S_m^3, S_n^{\pm}] = \pm S_m^{\pm} \, \delta_{mn} \>, \quad [S_m^+,S_n^-] = 2 S_m^3 \, \delta_{mn} \>. 
\end{equation}





%
%
%
%
\section{The non-unitary $so(3)$ classical r-matrix \label{App: Prelim2}}
%
%
The cornerstone of our study presented in this paper is the non-unitary $so(3)$ classical r-matrix \eqref{Eq: rK}. Here, in the Appendix \ref{App: Prelim2}, we recount how the r-matrix \eqref{Eq: rK} can be obtained starting from the unitary, $so(3)$ invariant classical r-matrix
\begin{equation}
\label{Eq: r-so3}
r (\lambda) = - \frac{\vec{\mathscr{S}}_{1} \cdot \vec{\mathscr{S}}_{2}}{\lambda} = - \frac{\mathcal{P} - 3 \mathcal{K}}{\lambda} \>,
\end{equation}
where we have used the notation introduced in the Appendix \ref{App: Prelim1}. In particular,  the classical r-matrix \eqref{Eq: r-so3} can be obtained as a quasi-classical limit of the $SO(3)$ quantum R-matrix \cite{ZZ78, ZZ79, KulishSklyanin82}. Evidently, this  classical r-matrix satisfies the classical Yang-Baxter equation \cite{BelavinDrinfeld}
\begin{equation}
\label{Eq: cl-YBE}
 [r_{12}(\lambda - \mu) , r_{13} (\lambda - \nu)] +  [r_{12}(\lambda - \mu) , r_{23}(\mu - \nu)] + [r_{13} (\lambda - \nu) , r_{23} (\mu - \nu) ] = 0 ,
\end{equation}
and has the unitarity property 
\begin{equation}
\label{Rq: unit-r}
r_{21}( - \lambda ) = - r_{12} ( \lambda) .
\end{equation}


We also consider the following reflection matrix
\begin{equation}
\label{Eq: K-mat}
K(\lambda) = %\frac{1}{(\xi - \lambda \sqrt{\psi \varphi  + \nu ^2} )^2}
\left(\begin{array}{ccc}
(\xi  - \nu \lambda  )^2 & - \sqrt{2} \, \psi  \, \lambda  (\xi  - \nu \lambda) & \psi ^2 \lambda ^2 \\[1ex] 
- \sqrt{2} \, \varphi \lambda (\xi  - \nu \lambda) & \xi ^2 + ( \psi \varphi  - \nu ^2) \, \lambda^2  &  - \sqrt{2} \, \psi \lambda (\xi  + \nu \lambda) \, \\[1ex] 
\varphi ^2 \lambda ^2  & - \sqrt{2} \, \varphi \lambda (\xi  + \nu \lambda) & (\xi  + \nu \lambda)^2 
\end{array}\right) \>,
\end{equation}
here $\xi , \nu , \psi , \varphi$ are arbitrary parameters. As it is well known, this K-matrix can be obtained by the so-called fusion procedure \cite{KulishSklyanin82, Mezincescu:1991ke, Kulish:2009cx}, starting from the $s\ell(2)$ K-matrix \cite{CiriloAntonio:2014ymf, VegaGonzalez, Zamolodchikov94, Zamolodchikov94b}.  This method is outlined, in the trigonometric $so(3)$, in \cite{Salom:2020dxt}. Alternatively, the so-called scaling limit \cite{Kulish:2009aa} can be used to obtain the K-matrix \eqref{Eq: K-mat} from the trigonometric $so(3)$ boundary K-matrix \cite{Salom:2020dxt, Inami:1996np}. Evidently, this K-matrix satisfies the classical reflection equation 
\begin{equation}
\label{Eq: c-RE}
\begin{split}
r_{12}(\lambda - \mu) K _1(\lambda) K_2(\mu) + K _1(\lambda) r_{21}(\lambda + \mu) K_2(\mu) = \\
= K_2(\mu) r_{12}(\lambda + \mu) K _1(\lambda) +  K_2(\mu) K _1(\lambda) r_{21}(\lambda - \mu) .
\end{split}
\end{equation}
It is worth mentioning that, while in the context of Heisenberg's open spin chain, one should also consider the dual reflection equation, this is not the case in the Gaudin model. Namely, as a consequence of long-range Gaudin model interactions, the ``two ends of the chain'' cannot have the same interpretation as in the case of Heisenberg's spin chain. In the Gaudin case, boundary parameters must be fixed in a way that the reflection equation and its dual effectively degenerate into a single equation. 

Therefore, it follows that the corresponding non-unitary classical r-matrix, given by \cite{CiriloAntonio:2014bfg, Skrypnyk06, Skrypnyk07, Skrypnyk09, Skrypnyk10, Skrypnyk10a, Skrypnyk13, Skrypnyk16}
\begin{equation}
\label{Eq: rK-def}
r_{12}^{K} (\lambda , \mu) = r_{12} (\lambda - \mu) - K_{2} (\mu) r_{12} (\lambda + \mu) K_{2}^{-1} (\mu) \>,
\end{equation}
satisfies the generalized classical Yang-Baxter equation \cite{CiriloAntonio:2014bfg, Maillet85, Maillet86, BabelonViallet90, AvanTalon90}
\begin{equation}
\label{Eq: g-c-YBE}
\left[ r^K_{32}(\nu , \mu) , r^K_{13} (\lambda , \nu ) \right] + \left[  r^K_{12}(\lambda , \mu ) ,  r^K_{13} (\lambda ,\nu )  \right]  + \left[  r^K_{12}(\lambda , \mu ) ,   r^K_{23} ( \mu , \nu ) \right] = 0. 
\end{equation}
The explicit form of this non-unitary $so(3)$ classical r-matrix is the following 
\begin{equation}
\label{Eq: rK}
r_{12}^{K} (\lambda , \mu) =  - \left( \frac{\vec{\mathscr{S}}_{1} \cdot \vec{\mathscr{S}}_{2}}{\lambda - \mu} - \frac{\vec{\mathscr{S}}_{1} \cdot  \left( K_{2} (\mu) \vec{\mathscr{S}}_{2}K_{2}^{-1}(\mu) \right)}{\lambda + \mu} \right) .
\end{equation}







%\label{}

%% The Appendices part is started with the command \appendix;
%% appendix sections are then done as normal sections
%% \appendix

%% \section{}
%% \label{}

%% If you have bibdatabase file and want bibtex to generate the
%% bibitems, please use
%%
%%  \bibliographystyle{elsarticle-num} 
%%  \bibliography{<your bibdatabase>}

%% else use the following coding to input the bibitems directly in the
%% TeX file.

\begin{thebibliography}{00}

%% \bibitem{label}
%% Text of bibliographic item

\bibitem{Heisenberg28} W. Heisenberg, \textsl{Zur Theorie der Ferromagnetismus}, Zeitschrift f{\"u}r Physik 49, 619--636 (1928).
\bibitem{Gaudin76} M. ~Gaudin, \textsl{Diagonalisation d'une classe d'hamiltoniens de spin}, J. Physique 37 (1976) 1087--1098.
\bibitem{Gaudin83} M. ~Gaudin, \textsl{La fonction d'onde de {B}ethe}, Masson, Paris, 1983.
\bibitem{Gaudin-English} M. ~Gaudin, \textsl{The Bethe Wavefunction}, Cambridge University Press, 2014.
\bibitem{TakhtajanFaddeev79} L. ~A. ~Takhtajan and ~L. ~D. ~Faddeev, \textsl{The quantum method for the inverse problem and the $XYZ$ Heisenberg model}, (in Russian) Uspekhi Mat. Nauk \textbf{34} No. 5 (1979) 13--63; translation in Russian Math. Surveys \textbf{34} No.5 (1979) 11--68.
\bibitem{KulishSklyanin82} P. ~P. ~Kulish and ~E. ~K. ~Sklyanin, \textsl{Quantum spectral transform method. Recent developments}, Lect. Notes Phys. \textbf{151} (1982), 61--119.
\bibitem{Faddeev98} L. ~D. ~Faddeev,  \textsl{How the algebraic {B}ethe {A}nsatz works for integrable models}, In \textrm{Quantum symmetries / Symetries quantiques}, Proceedings of the \textrm{Les Houches summer school, Session LXIV}.
Eds. A. Connes, K. Gawedzki and J. Zinn-Justin. North-Holland, 1998, 149--219; \texttt{hep-th/9605187}.
%\cite{Sklyanin:1987ih}
\bibitem{Sklyanin:1987ih} E.~K.~Sklyanin, \emph{Separation of variables in the Gaudin model}, Zap. Nauchn. Semin. \textbf{164} (1987), 151-169;
doi:10.1007/BF01840429
%113 citations counted in INSPIRE as of 08 Sep 2021  

\bibitem{KulSkl80} P. ~P. ~Kulish and ~E. ~K. ~Sklyanin, \emph{Solutions of the Yang-Baxter equation}, (Russian) Zap. Nauchn. Sem. Leningrad. Otdel. Mat. Inst. Steklov. (LOMI) 95 (1980), 129--160; translation in  J. Soviet Math. Vol. 19 (1982), 1596--1620.
\bibitem{BelavinDrinfeld} A. ~A. ~Belavin and  V. ~G. ~Drinfeld. \textsl{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.

\bibitem{Jurco90} B. ~Jur\v{c}o, \textsl{Classical Yang-Baxter equations and quantum integrable systems (Gaudin models)}, in Quantum groups (Clausthal, 1989), Lecture Notes in Phys. Volume 370 (1990) 219--227.
\bibitem{Semenov97} M. ~A. ~Semenov-Tian-Shansky, \textsl{Quantum and classical integrable systems}, in Integrability of Nonlinear Systems, Lecture Notes in Physics Volume 495 (1997) 314-377.
%\bibitem{Jurco89} B. ~Jur\v{c}o, \textsl{Classical {Y}ang-{B}axter equations and quantum integrable systems}, J. Math. Phys. Volume 30 (1989) 1289--1293.
\bibitem{WagnerMacfarlane00} F. ~Wagner and A. ~J. ~Macfarlane,  \textsl{Solvable Gaudin models for higher rank symplectic algebras. Quantum groups and integrable systems (Prague, 2000)} Czechoslovak J. Phys. 50 (2000) 1371--1377.
\bibitem{BrzezinskiMacfarlane94} T. ~Brzezinski and A. ~J. ~Macfarlane,  \textsl{On integrable models related to the osp(1,2) Gaudin algebra}, J. Math. Phys. 35 (1994), no. 7, 3261--3272.

%\cite{Kulish:2001hr}
\bibitem{Kulish:2001hr} P.~P.~Kulish and N.~Manojlovic, \emph{Bethe vectors of the osp(1|2) Gaudin model}, Lett. Math. Phys. \textbf{55} (2001), 77-95; doi:10.1023/A:1010950003268
%\cite{Kulish:2001vs}
\bibitem{Kulish:2001vs} P.~P.~Kulish and N.~Manojlovic, \emph{Creation operators and Bethe vectors of the osp(1|2) Gaudin model}, J. Math. Phys. \textbf{42} (2001), 4757-4778; doi:10.1063/1.1398584
%19 citations counted in INSPIRE as of 08 Sep 2021
%\cite{Kulish:2003da}
\bibitem{Kulish:2003da} P.~P.~Kulish and N.~Manojlovic, \emph{Trigonometric $osp(1 | 2)$ Gaudin model}, J. Math. Phys. \textbf{44} (2003), 676-700; doi:10.1063/1.1531250
%16 citations counted in INSPIRE as of 08 Sep 2021
\bibitem{LimaUtiel01} A.~Lima-Santos and W.~Utiel, \textsl{Off-shell {B}ethe ansatz equation for $osp(2|1)$ {G}audin magnets}, Nucl. Phys. B  600 (2001) 512--530.


\bibitem{KZ84} V. G. Knizhnik and A. B. Zamolodchikov, \emph{Current algebras and Wess-Zumino model in two dimensions}, Nucl. Phys. B 247 (1984) 83--103.
\bibitem{ReshetikhinVarchenko} N. ~Reshetikhin and A. ~Varchenko, \emph{Quasiclassical asymptotics of solutions to the KZ equations},  Geometry, topology \& physics for Raul Bott,  Conference Proceedings Lecture Notes Geometry Topology VI, Int. Press, Cambridge, MA., pp. 293--322, (1995). 
\bibitem{FeginFrenkelReshetikhin} B. ~Fegin, E. ~Frenkel and N. ~Reshetikhin, \emph{Gaudin Model, Bethe Ansatz and Critical Level},  Commun. Math. Phys. \textbf{166 }(1994) 27--62.
\bibitem{BabujianFlume} H. M. ~Babujian and R. ~Flume, \emph{Off-shell Bethe ansatz equations for Gaudin magnets and\break solutions of {K}nizhnik-{Z}amolodchikov equations}, Modern Physics Letters A, \textbf{Vol 9} No. 22 (1994) 2029--2039.
\bibitem{KurakLima04} V.~Kurak and A.~Lima-Santos, \textsl{$sl(2|1)^{(2)}$ Gaudin magnet and its associated Knizhnik-Zamolodchikov equation}, Nuclear Physics B 701 (2004) 497--515.


\bibitem{HikamiKulishWadati92} K. ~Hikami, P. ~P. ~Kulish and M. ~Wadati, \textsl{Integrable Spin Systems with Long-Range Interaction}, Chaos, Solitons \& Fractals  Vol. 2 No. 5 (1992) 543--550.
\bibitem{HikamiKulishWadati92a} K. ~Hikami, P. ~P. ~Kulish and M. ~Wadati, \textsl{Construction of Integrable Spin Systems with Long-Range Interaction}, J. Phys. Soc. Japan Vol. 61 No. 9 (1992) 3071--3076.

\bibitem{SklyaninTakebe96} E.~K. ~Sklyanin and T. ~Takebe, \textsl{Algebraic Bethe ansatz for the XYZ Gaudin model}, Phys. Lett. A 219 (1996) 217-225.

%\cite{CiriloAntonio:2005qts}
\bibitem{CiriloAntonio:2005qts}
N.~Cirilo Ant\'onio and N.~Manojlovi\'c, \emph{$sl_2$ Gaudin model with jordanian twist}, J. Math. Phys. \textbf{46} (2005) no.10, 102701;
doi:10.1063/1.2036932
\bibitem{KMSS}  P.~P.~Kulish,  N. ~Manojlovi\'c, M. ~Samsonov and A. ~Stolin, \emph{Bethe Ansatz for deformed Gaudin model}, (AGMF Tartu08 Workshop Proceedings) in the Proceedings of the Estonian Academy of Sciences \textbf{Vol 59} No. 4 (2010) 326--331.
%\cite{Cirilo-Antonio:2010xmg}
\bibitem{Cirilo-Antonio:2010xmg} 
N.~Cirilo Antonio, N.~Manojlovic and A.~Stolin, \emph{Algebraic Bethe Ansatz for deformed Gaudin model},
J. Math. Phys. \textbf{52} (2011), 103501; doi:10.1063/1.3644345
%[arXiv:1002.4951 [nlin.SI]].
%3 citations counted in INSPIRE as of 13 Sep 2021


\bibitem{Schlesinger} L. ~Schlesinger, \"Uber eine Klasse von Differentialsystemen beliebiger Ordnung mit festen kritischen Punkten, J. Reine Angew. Math. 141 (1912) 96-145.
\bibitem{Jimbo:1981aa}
M.~Jimbo, T.~Miwa and a.~K.~Ueno, \emph{Monodromy preserving deformation of linear ordinary differential equations with rational coefficients: I. General theory and $\tau$-function}, Physica D \textbf{2} (1981) 306-352. 
\bibitem{Jimbo:1981bb}
M.~Jimbo and T.~Miwa, \emph{Monodromy perserving deformation of linear ordinary differential equations with rational coefficients. II},
Physica D \textbf{2} (1981) 407-448.
\bibitem{Jimbo:1981cc}
M.~Jimbo and T.~Miwa, \emph{Monodromy preserving deformation of linear ordinary differential equations with rational coefficients. III},
Physica D \textbf{4} (1981), 26-46.
%\cite{Korotkin:1999xx}
\bibitem{Korotkin:1999xx}
D.~Korotkin, N.~Manojlovic and H.~Samtleben, \emph{Schlesinger transformations for elliptic isomonodromic deformations}, J. Math. Phys. \textbf{41} (2000), 3125-3141; doi:10.1063/1.533296
%\cite{Manojlovic:2001mv}
\bibitem{Manojlovic:2001mv}
N.~Manojlovic and H.~Samtleben, \emph{Schlesinger transformations and quantum R matrices}, Commun. Math. Phys. \textbf{230} (2002), 517-537;
doi:10.1007/s00220-002-0716-1
%2 citations counted in INSPIRE as of 13 Sep 2021
\bibitem{Dubrovin2007}
B.~Dubrovin and M.~ Mazzocco, \emph{Canonical Structure and Symmetries of the Schlesinger \break Equations}, Commun. Math. Phys. 271, 289--373 (2007);  doi:10.1007/s00220-006-0165-3


\bibitem{Sklyanin88} E. ~K. ~Sklyanin, \textsl{Boundary conditions for integrable quantum systems}, J. Phys. A: Math. Gen. \textbf{21} (1988) 2375--2389.
\bibitem{FreidelMaillet91}  L.~Freidel and  J.M.~Maillet, \emph{Quadratic algebras and integrable systems}, Physics  Letters B \textbf{262} (1991) 278--284.
\bibitem{FreidelMaillet91a} L. ~Freidel  and J.M.~Maillet,  \textsl{On classical and quantum integrable field theories associated to Kac-Moody current algebras}, Phys. Lett. \textbf{B 263 }(1991) 403--410.


\bibitem{Martins05} C. ~S. ~Melo,  G. ~A. ~P. ~Ribeiro and M. ~J. ~Martins,  \textsl{Bethe ansatz for the $XXX-S$ chain with non-diagonal open boundaries}, Nuclear Phys. B 711, no. 3 (2005) 565--603.
\bibitem{Eric13} S. ~Belliard, N. ~Cramp\'e and E. ~Ragoucy, \textsl{Algebraic Bethe ansatz for open XXX model with triangular boundary matrices}, Lett. Math. Phys. 103 No. 5 (2013) 493--506. 
\bibitem{Belliard13} S. ~Belliard and N. ~Cramp\'e \textsl{Heisenberg XXX model with general boundaries: eigenvectors from algebraic Bethe ansatz}, SIGMA Symmetry Integrability Geom. Methods Appl. 9 (2013), Paper 072, 12 pp.
\bibitem{Lima13} R. ~A. ~Pimenta and A. ~Lima-Santos, \textsl{Algebraic Bethe ansatz for the six vertex model with upper triangular K-matrices}, J. Phys. A 46 No. 45 (2013)  455002, 13 pp.
\bibitem{Belliard15} S. ~Belliard, \emph{Modified algebraic Bethe ansatz for XXZ chain on the segment -  I: Triangular cases}, Nuclear Physics \textbf{B 892} (2015) 1--20.
\bibitem{BelliardPimenta15} S. ~Belliard and R. ~A. ~Pimenta \emph{Modified algebraic Bethe ansatz for XXZ chain on the segment - II - general cases}, Nuclear Physics \textbf{B 894} (2015) 527--552.
\bibitem{AvanEtal15} J. ~Avan, S. ~Belliard, N. ~Grosjean and R. ~A. ~Pimenta, \emph{Modified algebraic Bethe ansatz for XXZ chain on the segment - III - Proof}, Nuclear Physics \textbf{B 899} (2015) 229--246.
\bibitem{Nepomechie15} A. ~M. ~Gainutdinov and R. ~I. ~Nepomechie, \textsl{Algebraic Bethe ansatz for the quantum group invariant open XXZ chain at roots of unity}, Nuclear Physics B 909 (2016) 796--839. 
\bibitem{Eric2013RMP} E. ~Ragoucy, \textsl{Coordinate Bethe ans\"atze for non-diagonal boundaries}, Rev. Math. Phys. 25 (2013), no. 10, 1343007.
\bibitem{LucEricRafael07} L. ~Frappat, R. ~I. ~Nepomechie, and E. ~Ragoucy, \textsl{A complete Bethe ansatz solution for the open spin-s XXZ chain with general integrable boundary terms}, Journal of Statistical Mechanics: Theory and Experiment,  0709 (2007)  P09009.
\bibitem{China13} J. ~Cao, W.- L. ~Yang, K. ~Shi and Y. ~Wang, \textsl{Off-diagonal Bethe ansatz solution of the XXX spin chain with arbitrary boundary conditions}, Nuclear Physics B 875 (2013) 152--165.
\bibitem{China15} X. ~Zhang, Y.-Y. ~Li, J. ~Cao, W.-L. ~Yang, K. ~Shi and Y. ~Wang, \textsl{Bethe states of the XXZ spin-$\textstyle{\frac{1}{2}}$ chain with arbitrary boundary fields} Nuclear Physics B 893 (2015) 70--88. 
\bibitem{Nepomechie04} R. ~I. ~Nepomechie, \textsl{Bethe ansatz solution of the open XXZ chain with nondiagonal boundary terms}, Journal of Physics A: Mathematical and Theoretical 37 (2004), no. 2, 433--440. 
\bibitem{Jimbo95} M. ~Jimbo, R. ~Kedem, T. ~Kojima, H. ~Konno and T. ~Miwa, \textsl{XXZ chain with a boundary}, Nuclear Physics B 441 (1995) 437--470. 
\bibitem{China03} J. ~Cao, H. ~Lin, K. ~Shi and Y. ~Wang, \textsl{Exact solutions and elementary excitations in the XXZ spin chain with unparallel boundary fields}, Nucl. Phys. B 663 (2003)  487--519.
%\cite{Kulish:2009ci}
\bibitem{Kulish:2009ci}
P.~P.~Kulish, N.~Manojlovic and Z.~Nagy, \emph{Jordanian deformation of the open XXX-spin chain}, Theor. Math. Phys. \textbf{163} (2010), 644-652;
doi:10.1007/s11232-010-0047-x
%[arXiv:0911.5592 [nlin.SI]].
%7 citations counted in INSPIRE as of 13 Sep 2021
\bibitem{Maillet17} J.M. ~Maillet, G. ~Niccoli and B. Pezelier, \textsl{Transfer matrix spectrum for cyclic representations of the 6-vertex reflection algebra I}, SciPost Phys. 2, 009 (2017).



%\cite{CiriloAntonio:2014ymf}
\bibitem{CiriloAntonio:2014ymf}
N.~Cirilo Ant\'onio, N.~Manojlovi\'c and I.~Salom, \emph{Algebraic Bethe ansatz for the XXX chain with triangular boundaries and Gaudin model},
Nucl. Phys. B \textbf{889} (2014), 87-108; doi:10.1016/j.nuclphysb.2014.10.014
%16 citations counted in INSPIRE as of 09 Sep 2021
%\cite{Manojlovic:2017nbh}
\bibitem{Manojlovic:2017nbh} N.~Manojlovi\'c and I.~Salom, \emph{Algebraic Bethe ansatz for the XXZ Heisenberg spin chain with triangular boundaries and the corresponding Gaudin model}, Nucl. Phys. B \textbf{923} (2017), 73-106; doi:10.1016/j.nuclphysb.2017.07.017
%8 citations counted in INSPIRE as of 09 Sep 2021
%\cite{CiriloAntonio:2013eng}
\bibitem{CiriloAntonio:2013eng}
N.~Cirilo Ant\'onio, N.~Manojlovi\'c and Z.~Nagy, \emph{Trigonometric $sl(2)$ Gaudin model with boundary terms},
Rev. Math. Phys. \textbf{25} no.10 (2013) 1343004;
doi:10.1142/S0129055X13430046
%11 citations counted in INSPIRE as of 22 Sep 2021
\bibitem{Hikami95} K. ~Hikami, \textsl{Gaudin magnet with boundary and generalized Knizhnik-Zamolodchikov equation}, J. Phys. A Math. Gen. \textbf{28} (1995) 4997--5007.
%\cite{CiriloAntonio:2013rbz}
\bibitem{CiriloAntonio:2013rbz}
N.~Cirilo Ant\'onio, N.~Manojlovi\'c and Z.~Nagy, \emph{Jordanian deformation of the open $s\ell(2)$ Gaudin model},
Theor. Math. Phys. \textbf{179} (2014), 462-471; doi:10.1007/s11232-014-0155-0
%[arXiv:1304.6918 [nlin.SI]].
%2 citations counted in INSPIRE as of 13 Sep 2021
\bibitem{YangZhangSasakic04} W. ~L. Yang,  R. ~Sasaki and Y. ~Z. ~Zhang, \textsl{$\mathbb{Z}_n$ elliptic Gaudin model with open boundaries}, JHEP 09 (2004) 046.
\bibitem{YangZhangSasakic05} W. ~L. Yang,  R. ~Sasaki and Y. ~Z. ~Zhang, \textsl{$A_{n-1}$ Gaudin model with open boundaries}, Nuclear Physics B 729 (2005) 594--610.
\bibitem{YangZhang12} K. ~Hao, W.-L. ~Yang, H. ~Fan, S. ~Y. ~Liu, K. ~Wu, Z. ~Y. ~Yang and Y. ~Z. ~Zhang, \textsl{Determinant representations for scalar products of the XXZ Gaudin model with general boundary terms}, Nuclear Physics \textbf{B 862} (2012) 835--849.
\bibitem{Lima09}A.~Lima-Santos, \textsl{The $sl(2|1)^{(2)}$ Gaudin magnet with diagonal boundary terms}, J. Stat. Mech. (2009) P07025.

%\cite{Manojlovic:2015jhi}
\bibitem{Manojlovic:2015jhi} N. ~Manojlovi\'c,  Z. ~Nagy and I. ~Salom, \emph{Derivation of the trigonometric Gaudin Hamiltonians}, Proceedings of the 8th Mathematical Physics meeting: Summer School and Conference on Modern Mathematical Physics, 24 - 31 August 2014, Belgrade, Serbia, SFIN \textbf{XXVIII} Series A: Conferences No. A1(2015) 127--135; ISBN: 978-86-82441-43-4.
%0 citations counted in INSPIRE as of 20 Sep 2021
\bibitem{Manojlovic:2018abc}  N. ~Manojlovi\'c, N. ~Cirilo Ant\'onio  and I. ~Salom, \emph{Quasi-classical limit of the open Jordanian XXX spin chain}, Proceedings of the 9th Mathematical Physics meeting: Summer School and Conference on Modern Mathematical Physics, 18 - 23 September 2017, Belgrade,\break Serbia, SFIN \textbf{XXXI} Series A: Conferences No. A1 (2018) 259--266; ISBN: 978-86-82441-48-9.


%\cite{CiriloAntonio:2014bfg}
\bibitem{CiriloAntonio:2014bfg}N.~Cirilo Ant\'onio, N.~Manojlovi\'c, E.~Ragoucy and I.~Salom, \emph{Algebraic Bethe ansatz for the $s\ell$(2) Gaudin model with boundary}, Nucl. Phys. B \textbf{893} (2015), 305-331; doi:10.1016/j.nuclphysb.2015.02.011
%11 citations counted in INSPIRE as of 08 Sep 2021

\bibitem{Sklyanin86} E. ~K. ~Sklyanin, \textsl{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.
\bibitem{Sklyanin87} E. ~K. ~Sklyanin, \textsl{Boundary conditions for integrable systems}, in the Proceedings of the VIIIth international congress on mathematical physics (Marseille, 1986), World Sci. Publishing, Singapore, (1987) 402--408.
\bibitem{Maillet85} J.M. ~Maillet, \textsl{Kac-Moody algebra and extended Yang-Baxter relations in the $O(N)$ non-linear $\sigma$-model}, Physics Letters B \textbf{162} (1985) 137--142.
\bibitem{Maillet86} J.M. Maillet, \textsl{New integrable canonical structures in two-dimensional models}, Nuclear Physics \textbf{B269} (1986) 54--76.
\bibitem{BabelonViallet90} O.~Babelon and C.~Viallet, \emph{Hamiltonian structures and Lax equations} Physics  Letters B  {\bf 237}, 411 (1990).
\bibitem{AvanTalon90}  J. ~Avan and  M. ~Talon, \emph{Rational and trigonometric constant non-antisymmetric r-matrices},  Physics  Letters B \textbf{241}, (1990) 77--82.



\bibitem{Skrypnyk06}  T.~Skrypnyk, \emph{Generalized quantum Gaudin spin chains, involutive automorphisms and twisted classical r-matrices}, J. Math. Phys. \textbf{47}, 033511 (2006); https://doi.org/10.1063/1.2179052
\bibitem{Skrypnyk07}  T.~Skrypnyk, \emph{Generalized Gaudin spin chains, non-skew-symmetric r-matrices and reflection equation algebras}, J. Math. Phys. 48, 113521 (2007); https://doi.org/10.1063/1.2816256
\bibitem{Skrypnyk09} T. ~Skrypnyk, \textsl{Non-skew-symmetric classical r-matrix, algebraic Bethe ansatz, and Bardeen-Cooper-Schrieffer-type integrable systems}, J. Math. Phys. 50 (2009) 033540, 28 pages.
\bibitem{Skrypnyk10} T.~Skrypnyk, \emph{Generalized Knizhnik-Zamolodchikov equations, off-shell Bethe ansatz and non-skew-symmetric classical r-matrices},  Nuclear Physics \textbf{B 824} (2010) 436-451;
\bibitem{Skrypnyk10a} T.~Skrypnyk,  \emph{Isomonodromic deformations, generalized Knizhnik-Zamolodchikov equations and non-skew-symmetric classical r-matrices}, Journal of Mathematical Physics \textbf{51}, 083516 (2010);
\bibitem{Skrypnyk13}T. ~Skrypnyk,  \textsl{"Z2-graded'' Gaudin models and analytical Bethe ansatz}, Nuclear Physics B 870 (2013), no. 3, 495--529.
\bibitem{Skrypnyk16}T. ~Skrypnyk, \textsl{"Generalized" algebraic Bethe ansatz, Gaudin-type models and Zp-graded classical r-matrices}, Nuclear Physics B 913 (2016) 327--356. 



%\cite{Manojlovic:2020xou}
\bibitem{Manojlovic:2020xou}
N.~Manojlovic, I.~Salom and N.~Cirilo Ant\'onio, \emph{XYZ Gaudin model with boundary terms}, Proceedings of the 10th Mathematical Physics meeting: Summer School and Conference on Modern Mathematical Physics, 9 - 14 September 2019, Belgrade, Serbia,
SFIN  \textbf{XXXIII} Series A: Conferences, No. A1 (2020)  143 -- 160; ISBN 978-86-82441-51-9.
%0 citations counted in INSPIRE as of 15 Sep 2021


%\cite{Skrypnyk:2021xxc}
\bibitem{Skrypnyk:2021xxc}
T.~Skrypnyk and N.~Manojlovi\'c, \emph{Twisted rational r -matrices and algebraic Bethe ansatz: Application to generalized Gaudin and Richardson models}, Nucl. Phys. B \textbf{967} (2021), 115424; doi:10.1016/j.nuclphysb.2021.115424
%0 citations counted in INSPIRE as of 09 Sep 2021



%\cite{Salom:2015wok}    
\bibitem{Salom:2015wok}
I.~Salom and N.~Manojlovi\'c, \emph{Creation operators of the non-periodic $sl(2)$ Gaudin model}, Proceedings of the 8th Mathematical Physics meeting: Summer School and Conference on Modern Mathematical Physics, 24 - 31 August 2014, Belgrade, Serbia, SFIN \textbf{XXVIII} Series A: Conferences No. A1 (2015) 149--155; ISBN: 978-86-82441-43-4.
%0 citations counted in INSPIRE as of 21 Sep 2021
%\cite{Salom:2019myz}
\bibitem{Salom:2019myz} I.~Salom, N.~Manojlovi\'c and N.~Cirilo Ant\'onio, \emph{Generalized $s\ell$ (2) Gaudin algebra and corresponding Knizhnik\textendash{}Zamolodchikov equation}, Nucl. Phys. B \textbf{939} (2019), 358-371; doi:10.1016/j.nuclphysb.2018.12.025
%4 citations counted in INSPIRE as of 08 Sep 2021
%\cite{Manojlovic:2017yxv} 
\bibitem{Manojlovic:2017yxv}
N.~Manojlovi\'c and I.~Salom, \emph{Algebraic Bethe Ansatz for the Trigonometric s$\ell$(2) Gaudin Model with Triangular Boundary}, 
Symmetry \textbf{12} no.3 (2020) 352; doi:10.3390/sym12030352
%[arXiv:1709.06419 [nlin.SI]].
%7 citations counted in INSPIRE as of 21 Sep 2021
%\cite{Salom:2021gox}
\bibitem{Salom:2021gox} I.~Salom and N.~Manojlovi\'c, \emph{Bethe states and Knizhnik-Zamolodchikov equations of the trigonometric Gaudin model with triangular boundary}, Nucl. Phys. B \textbf{969} (2021), 115462; doi:10.1016/j.nuclphysb.2021.115462
%0 citations counted in INSPIRE as of 08 Sep 2021

%\cite{Mezincescu:1991ke}
\bibitem{Mezincescu:1991ke} L.~Mezincescu and R.~I.~Nepomechie, \emph{Fusion procedure for open chains}, J. Phys. A \textbf{25} (1992) 2533-2544.
%92 citations counted in INSPIRE as of 23 Sep 2021

%\cite{Kulish:2009cx}
\bibitem{Kulish:2009cx} P.~P.~Kulish, N.~Manojlovic and Z.~Nagy, \emph{Symmetries of spin systems and Birman-Wenzl-Murakami algebra}, J. Math. Phys. \textbf{51} (2010), 043516; doi:10.1063/1.3366259
%[arXiv:0910.4036 [nlin.SI]].
%8 citations counted in INSPIRE as of 23 Sep 2021


\bibitem{VegaGonzalez} H. ~J. ~de Vega and A. ~Gonz\'alez ~Ruiz, \textsl{Boundary $K$-matrices for the $XYZ$, $XXZ$, $XXX$ spin chains}, J. Phys. A: Math. Gen. \textbf{27} (1994), 6129--6137.
\bibitem{Zamolodchikov94} S. ~Ghoshal and A. ~B. ~Zamolodchikov, \textsl{Boundary S-matrix and boundary state in two-dimensional integrable quantum field theory}, International Journal of Modern Physics A \textbf{09}, 3841 (1994) 3841--3885.
\bibitem{Zamolodchikov94b} S. ~Ghoshal and A. ~B. ~Zamolodchikov, \textsl{Errata: Boundary S-matrix and boundary state in two-dimensional integrable quantum field theory}, International Journal of Modern Physics A \textbf{09}, 4353 (1994) 4353. 


%\cite{Salom:2020dxt}
\bibitem{Salom:2020dxt}
I.~Salom, N.~Manojlovic and N.~Cirilo Ant\'onio, \emph{The spin 1 XXZ Gaudin model with boundary}, Proceedings of the 10th Mathematical Physics meeting: Summer School and Conference on Modern Mathematical Physics, 9 - 14 September 2019, Belgrade, Serbia,
SFIN \textbf{XXXIII} Series A: Conferences, No. A1 (2020) 277 -- 286; ISBN 978-86-82441-51-9.
%0 citations counted in INSPIRE as of 09 Sep 2021

\bibitem{Kulish:2009aa} P.~P.~Kulish, \emph{Twist Deformations of Quantum Integrable Spin Chains}, Lect. Notes Phys. \textbf{774} (2009) 165--188.

%\cite{Inami:1996np}
\bibitem{Inami:1996np}
T.~Inami, S.~Odake and Y.~Z.~Zhang, \emph{Reflection K matrices of the 19 vertex model and XXZ spin 1 chain with general boundary terms},
Nucl. Phys. B \textbf{470} (1996), 419-434; doi:10.1016/0550-3213(96)00133-2
%[arXiv:hep-th/9601049 [hep-th]].
%25 citations counted in INSPIRE as of 23 Sep 2021



\bibitem{ZZ78} Alexander B. Zamolodchikov and Alexey B. Zamolodchikov, \emph{Relativistic factorized S-matrix in two dimensions having $O(N)$ isotopic symmetry} Nuclear Phys. B 133 (1978) 525--535.
\bibitem{ZZ79} Alexander B. Zamolodchikov and Alexey B. Zamolodchikov, \emph{Factorized S-matrices in two dimensions as the exact solutions of certain relativistic quantum field theory models} Annals of Physics Vol. 120 (1979) 253--291.



%\cite{Sklyanin:1997zz}
\bibitem{Sklyanin:1997zz}
E.~K.~Sklyanin, \emph{Generating function of correlators in the sl(2) Gaudin model}, Lett. Math. Phys. \textbf{47} (1999), 275-292
doi:10.1023/A:1007585716273
%31 citations counted in INSPIRE as of 08 Sep 2021    


\end{thebibliography}
\end{document}
\endinput
%%
%% End of file `elsarticle-template-num.tex'.
