\documentclass[11pt]{article}
\usepackage{geometry} % see geometry.pdf on how to lay out the page. There's lots.
\geometry{a4paper} % or letter or a5paper or ... etc
%\geometry{margins=2in} % for example, change the margins to 2 inches all round
\usepackage[english]{babel}
\usepackage[english]{layout}
\usepackage{cite}
%\usepackage{showkeys}

% Palatino for rm and math | Helvetica for ss | Courier for tt
\usepackage{mathpazo} % math & rm
\linespread{1.05}        % Palatino needs more leading (space between lines)
\usepackage[scaled]{helvet} % ss
\usepackage{courier} % tt
\normalfont
\usepackage[T1]{fontenc}

%% LaTeX - Article customise
%%% PACKAGES
%\usepackage{booktabs} % for much better looking tables
\usepackage{array} % for better arrays (eg matrices) in maths
\usepackage{paralist} % very flexible & customisable lists (eg. enumerate/itemize, etc.)
\usepackage{verbatim} % adds environment for commenting out blocks of text & for better verbatim
%\usepackage{subfigure} % make it possible to include more than one captioned figure/table in a single float
% These packages are all incorporated in the memoir class to one degree or another...
%% LaTeX Preamble - Common packages

\usepackage[utf8]{inputenc} % Any characters can be typed directly from the keyboard, eg éçñ
\usepackage{textcomp} % provide lots of new symbols
\usepackage{graphicx}  % Add graphics capabilities
\usepackage{epstopdf} % to include .eps graphics files with pdfLaTeX
\usepackage{flafter}  % Don't place floats before their definition
%\usepackage{topcapt}   % Define \topcation for placing captions above tables (not in gwTeX)
%\usepackage{natbib} % use author/date bibliographic citations

\usepackage{amsfonts,latexsym,amsthm,amssymb}
\usepackage{amsmath,amssymb,bbm,amsfonts}  % Better maths support & more symbols
\usepackage{bm}  % Define \bm{} to use bold math fonts
\usepackage{mathtools}
\usepackage[mathscr]{eucal}

\usepackage{memhfixc}  % remove conflict between the memoir class & hyperref
%\usepackage[activate]{pdfcprot}  % Turn on margin kerning (not in gwTeX)
%\usepackage{pdfsync}  % enable tex source and pdf output syncronicity
%%% HEADERS & FOOTERS
\usepackage{fancyhdr} % This should be set AFTER setting up the page geometry
\pagestyle{fancy} % options: empty , plain , fancy
%\renewcommand{\headrulewidth}{.3pt} % customise the layout...
%\lhead{}\chead{}\rhead{}
%\lfoot{preliminary draft}
\cfoot{-\thepage-}
%\rfoot{\today}

%%% SECTION TITLE APPEARANCE
\usepackage{sectsty}
\allsectionsfont{\sffamily\mdseries\upshape} % (See the fntguide.pdf for font help)
% (This matches ConTeXt defaults)

%%% colours

\usepackage{color}
\newcommand{\red}{\color{red}}
\newcommand{\blue}{\color{blue}}
\newcommand{\gr}{\color{green}}
\newcommand{\mage}{\color{magenta}}



\def\vph{\varphi}
\def\eps{\varepsilon}
\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}}


%%%%%%

%\numberwithin{equation}{section}

%\renewcommand{\thesection}{\Roman{section}}
%\renewcommand{\thesubsection}{\thesection.\roman{subsection}}

%\newtheorem{theorem}{Theorem}[section]
%\newtheorem{proposition}{Proposition}[section]
%\newtheorem{remark}{Remark}[section]
%\newtheorem{corollary}{Corollary}[section]
%\newtheorem{lemma}{Lemma}[section]

%%%%%%%%%%%%%

\begin{document}

\section*{Reply to the referee's comments}

We thank the referee for the careful perusal of our paper and for the positive overall assessment.

In particular, referee's questions 2 and 3 initially prompted us to consider introducing large extensions of the paper, in order to provide fully detailed answers in the text. However, we realized that entirely covering these particular aspects of the problem here would divert the reader from the main results of this work, so we eventually opted to add only short clarifications in the present paper (while shortly postponing their more detailed treatment for another publication).

While we find the full discussion of 2 and 3 unsuitable for inserting in the present paper, we nevertheless provide it below, preceded by concisely addressing all of the issues brought up by the referee:

\begin{enumerate}
\item {\it The computer-output typesetting is very poor. Many formulas are too wide. The number of equations should be consistent with the style of NPB.}

We accept the referee's criticism, noting that, style-wise, we adhered to the "Your Paper Your Way" NPB policy, allowing the authors to choose arbitrary style during the peer review stage. However, for the revised version we now used Elsevier's template/style, took care of the width of formulas, and changed equation numbering - hopefully to good effect.


\item {\it Eqs.(B.4) and (B.5) give the reflection matrix and corresponding classical reflection equation, respectively. Does there exist the dual reflection matrix which satisfies the dual classical reflection equation?}

The short answer is: 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. Thus, 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. 

We briefly clarify this in Appendix B. But, to fully elucidate this matter requires considering the Gaudin model as a quasi-classical limit of the spin chain, which we show below.

\item {\it Is the construction (II.6) of the Lax operator of the model with boundary reflection unique? }

The Lax operator (II.6) can be obtained by following a relatively general procedure of quasi-classically expanding the Sklyanin monodromy (detailed below) and, in this sense, its form is fixed by construction. On the other hand, strictly speaking, we are not aware of any formal proof of its uniqueness (i.e. whether there is any freedom to nontrivially modify II.6 while still retaining effectively the same model).

In the text below (II.6) we have now mentioned that this form can be, in a straightforward manner, obtained by quasi-classical expansion of the Sklyanin monodromy of the corresponding open Heisenberg spin chain.


\item {\it Does the identity (III.21) give the constraints of Bethe roots? Is it included in the Bethe ansatz equations (III.26) or satisfied naturally?}

The identity (III.21) is purely algebraical (i.e. satisfied by an arbitrary set of complex values $\lambda$, $\mu_j$, $\xi$, $\psi$, $\phi$, $\nu$) and it has nothing to do with Bethe equations or their solutions. We now stressed this in the text.


\item {\it Is the scalar product Eq.(IV.35) a conjecture? Can it be proved by the mathematical induction? }

The off-shell scalar products of the Bethe vectors $\widetilde{\Phi} _M  (\mu _1 , \mu _2 ,  \ldots ,  \mu _M)$ (IV.8) in the $M=1$ case, the formula (IV.31), and in the $M=2$ case, the formulae (IV.32-34) were derived by a direct, straightforward calculations.  The formulae (IV.35-37) were 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 form (IV.35), strictly speaking, remains a conjecture - which we have now clearly stated in the text.

\end{enumerate}

Bellow, we address questions 2 and 3 in full detail, putting them in a more general context. We use the same notation as in the manuscript. At the end of the text we also include the list of references relevant for this discussion. 

\vspace{2cm}

We first notice that the $O(3)$ invariant R-matrix was found in \cite{ZZ78,ZZ79} (see also \cite{KulSkl80})
\begin{equation}
\label{ZZ-R}
 R (\lambda, \eta) = \lambda\left( \lambda + \eta \right ) \mathbbm{1} + 2 \eta
\left( \lambda +  \eta \right) \mathcal{P} - 6 \lambda \eta  \mathcal{K} \>,
\end{equation}
where $\lambda$ is a spectral parameter and  $\eta$ is a quasi-classical parameter. This R-matrix satisfies the Yang-Baxter equation 
\begin{equation}
\label{Eq: YBE}
R_{12}(\lambda - \mu) R_{13}(\lambda) R_{23}(\mu ) 
=  R_{23}(\mu ) R_{13}(\lambda) R_{12}(\lambda - \mu)  ,
\end{equation}
here the standard notation is used to denote spaces $V_j, j=1, 2, 3$ on which corresponding $R$-matrices $R_{ij}, ij = 12, 13, 23$ act non-trivially. In the present case $V_1 = V_2 = V_3 =\mathbb{C}^3$.

The R-matrix \eqref{ZZ-R} admits the spectral decomposition
\begin{align}
\label{spec-R}
\check{R}  (\lambda, \eta) &= \mathcal{P} R (\lambda, \eta) = 2 \eta \left( \lambda +  \eta \right) \mathbbm{1} + \lambda \left( \lambda + \eta \right ) \mathcal{P} - 6 \lambda \eta  \mathcal{K}    \notag \\
 &=  ( \lambda - \eta ) ( \lambda - 2 \eta) P_0  
- ( \lambda + \eta) ( \lambda - 2 \eta ) P_1 + ( \lambda + \eta ) ( \lambda + 2 \eta ) P_2 ,  
\end{align}
 using the projectors on the irreducible representation components in the Clebsch-Gordan decomposition
$V ^{(1)}\otimes V ^{(1)} = V ^{(2)}\oplus V ^{(1)}\oplus V ^{(0)}$
\begin{equation}
\label{projs}
P_2 = \frac{1}{2} \left( \mathbbm{1}  + \mathcal{P} - 2 \mathcal{K} \right)  , \quad
P_1 = \frac{1}{2} \left( \mathbbm{1}  - \mathcal{P} \right) , \quad
P_0 = \mathcal{K} .
\end{equation}

The R-matrix \eqref{ZZ-R} has some important properties: regularity, unitarity, PT-symmetry, and crossing symmetry.
The regularity condition at $\lambda = 0$ reads
\begin{equation}
\label{Eq: reg}
R(0, \eta) = 2 \eta ^2  \mathcal{P}.
\end{equation}
The unitarity relation is
\begin{equation}
\label{Eq: uni}
R_{12}( \lambda) R_{21}(-\lambda) = \rho ( \lambda) \mathbbm{1},
\end{equation}
with $\rho ( \lambda) = \lambda ^4 - 5 \lambda^2 \eta ^2 + 4 \eta^4$.
The so-called PT-symmetry states 
\begin{equation}
\label{Eq: PT-s}
R _{12}^t (\lambda) =  R _{21}(\lambda) \>.
\end{equation}
One particular consequence of \eqref{Eq: uni} and \eqref{Eq: PT-s} is the fact that $\rho ( \lambda)$ is necessarily an even function. Finally, the R-matrix \eqref{ZZ-R} has the following crossing symmetry property:
\begin{equation}
\label{Eq: crossing-s}
R ( \lambda) =  \left(\mathcal{J}\otimes \mathbbm{1}\right) R ^{t_2}( -\lambda -\eta ) \left(\mathcal{J}^{-1} \otimes \mathbbm{1}\right),
\end{equation}%
 where $t_2$ denotes the transpose in the second space and the matrix $\mathcal{J}$ is given by 
\begin{equation}
\label{J-mat}
\mathcal{J} = \left(\begin{array}{rrr}0 & 0 & -1 \\0 & 1 & 0 \\- 1 & 0 & 0\end{array}\right) \>.
\end{equation}
An important consequence of the above relations is the following identity, the so-called crossing unitarity,
\begin{equation}
\label{Eq: crossing-u}
R _{12}^{t_2} (\lambda) R _{12}^{t_1} (-\lambda-2\eta)  = \rho ( -\lambda-\eta ) \mathbbm{1} \>,
\end{equation}
here  $t_1$ and $t_2$ denote the transpose in the first and in the second space, respectively.
%\begin{equation}
%\label{M-matrix}
%M = \mathcal{J}^t \mathcal{J} .
%\end{equation}
%Moreover, $M^t = M$ and 
%\begin{equation}
%\label{Msymmetry}
%[M\otimes M, R (\lambda)] = 0 .
%\end{equation}


A way to introduce non-periodic boundary conditions which are compatible with the integrability of the bulk model, was developed in \cite{Sklyanin88, KMN3}. Boundary conditions on the left and right sites of the chain are encoded in the left and right reflection matrices $K^-$ and $K^+$. The compatibility condition between the bulk and the boundary of the system takes the form of the so-called reflection equation. It is written in the following form for the left reflection matrix acting on the space $\mathbb{C}^3$ at the first site $K^-(\lambda) \in \mathrm{End} (\mathbb{C}^3)$
\begin{equation}
\label{Eq: RE}
R_{12}(\lambda - \mu) K^-_1(\lambda) R_{21}(\lambda + \mu) K^-_2(\mu)=
K^-_2(\mu) R_{12}(\lambda + \mu) K^-_1(\lambda) R_{21}(\lambda - \mu) \>.
\end{equation}
Due to the properties of the R-matrix \eqref{Eq: uni}--\eqref{Eq: crossing-u} the dual reflection equation can be written in the following form
\begin{equation}
\label{Eq: d-RE}
R_{12}(-\lambda + \mu) K^+_1(\lambda) R_{21}(-\lambda - \mu -2\eta) K^+_2(\mu)=
K^+_2(\mu) R_{12}(-\lambda - \mu -2\eta) K^+_1(\lambda) R_{21}(-\lambda + \mu) \>.
\end{equation}
Evidently, 
\begin{equation}
\label{Eq: Kpl}
K^+ (\lambda) = K^- (-\lambda - \eta) 
\end{equation}
is a bijection between solutions to the reflection equation and the dual reflection equation, in the sense that the substitution of \eqref{Eq: Kpl} into the dual reflection equation \eqref{Eq: d-RE} yields the reflection equation \eqref{Eq: RE} with shifted arguments.







The left reflection matrix, which we have obtained by the fusion procedure \cite{KulishSklyanin82, Mezincescu:1991ke, Kulish:2009cx} starting from the $s\ell(2)$ K-matrix \cite{CiriloAntonio:2014ymf, VegaGonzalez, Zamolodchikov94, Zamolodchikov94b},
\begin{equation}
\label{Eq: K-min}
 K ^{-}(\lambda , \eta) = \big{(} k ^-_{ij} (\lambda , \eta) \big{)} \>, \qquad \text{where} \qquad i , j = 1 , 2 , 3 ,
\end{equation}
whose matrix elements are given below
\begin{equation}
\begin{aligned}
\label{Eq: K-min-com}
k ^-_{11} (\lambda , \eta)  &= 4 \xi ^{- \, 2} -  8 \xi^{-} \nu^{-} \lambda + (2\lambda - \eta) \left(2 \lambda \nu ^{- \, 2} + \eta (\nu ^{- \, 2} + \psi^{-} \varphi^{-}) \right)  \>, \\[1ex]
k ^-_{12} (\lambda , \eta)  &= -2 \sqrt{2} \psi^{-} \lambda \left(2 \xi^{-} - \nu^{-} \, (2\lambda - \eta) \right)  \>, \\[1ex]
k ^-_{13} (\lambda , \eta) &= 2 \psi^{- \, 2} \lambda \, (2\lambda - \eta)  \>, \\[2ex]
%
k ^-_{21} (\lambda , \eta)  &= -2 \sqrt{2} \varphi^{-} \lambda \left(2 \xi^{-} - \nu^{-} \, ( 2\lambda - \eta) \right) 
\>, \\[1ex]
k ^-_{22} (\lambda , \eta)  &= 4 \xi ^{- \, 2} - (2\lambda - \eta) \left( 2 \lambda \, ( \nu^{- \, 2} - \psi^{-} \varphi^{-} ) - \eta \, (\nu^{- \, 2} + \psi^{-} \varphi^{-}) \right) \>, \\[1ex]
k ^-_{23} (\lambda , \eta) &=  -2 \sqrt{2} \psi^{-}  \lambda \left( 2 \xi^{-}  + \nu^{-}  \, ( 2 \lambda - \eta) \right) \>, \\[2ex]
%
k ^-_{31} (\lambda , \eta)  &= \varphi^{- \, 2} \lambda \, (2\lambda - \eta) \>, \\[1ex]
k ^-_{32} (\lambda , \eta)  &= -2 \sqrt{2} \varphi^{-} \lambda \left( 2 \xi^{-} + \nu^{-} \, (2 \lambda - \eta) \right) \>, \\[1ex]
k ^-_{33} (\lambda , \eta)  &= 4 \xi ^{- \, 2} + 8 \xi^{-} \nu^{-} \lambda + (2\lambda - \eta) \left(2\lambda \nu^{- \, 2} + \eta \, (\nu^{- \, 2} + \psi^{-} \phi^{-}) \right)  \>, 
\end{aligned}
\end{equation}
satisfies the reflection equation \eqref{Eq: RE}. In particular, this K-matrix has the following important property
\begin{equation}
\label{Eq: Kmin-inv}
K ^{-}( - \lambda , \eta) \cdot  K ^{-}(\lambda , \eta) =  d ^- (\lambda , \eta ) \, d ^- (- \lambda , \eta ) \mathbbm{1}  \>, \\[2ex]
\end{equation}
%{%\small 
%\tiny
%\begin{equation}
%\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!
%%\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!
%%\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!
%\label{Eq: K-min}
%%\begin{split} 
%%&
%K ^{-}(\lambda) =  
%%\frac{1}{4 \xi^{- \, 2} - 8 \xi^{-} \lambda \sqrt{\nu^{-\, 2} + \psi^{-} \varphi^{-}} + (\nu^{- \, 2} + \psi^{-} \varphi^{-}) (4 \lambda ^2 - \eta ^2)} \times \\[2ex]
%%
%%&\times 
%\left(\begin{array}{ccc}
%4 \xi ^{- \, 2} -  8 \xi^{-} \nu^{-} \lambda + (2\lambda - \eta) \left(2 \lambda \nu ^{- \, 2} + \eta (\nu ^{- \, 2} + \psi^{-} \varphi^{-}) \right) & -2 \sqrt{2} \psi^{-} \lambda \left(2 \xi^{-} - \nu^{-} \, (2\lambda - \eta) \right) & 2 \psi^{- \, 2} \lambda \, (2\lambda - \eta)  \\[1ex]
%%
%-2 \sqrt{2} \varphi^{-} \lambda \left(2 \xi^{-} - \nu^{-} \, ( 2\lambda - \eta) \right) & 4 \xi ^{- \, 2} - (2\lambda - \eta) \left( 2 \lambda \, ( \nu^{- \, 2} - \psi^{-} \varphi^{-} ) - \eta \, (\nu^{- \, 2} + \psi^{-} \varphi^{-}) \right)  & -2 \sqrt{2} \psi^{-}  \lambda \left( 2 \xi^{-}  + \nu^{-}  \, ( 2 \lambda - \eta) \right) \\[1ex]
%%
%2 \varphi^{- \, 2} \lambda \, (2\lambda - \eta) & -2 \sqrt{2} \varphi^{-} \lambda \left( 2 \xi^{-} + \nu^{-} \, (2 \lambda - \eta) \right)  & 4 \xi ^{- \, 2} + 8 \xi^{-} \nu^{-} \lambda + (2\lambda - \eta) \left(2\lambda \nu^{- \, 2} + \eta \, (\nu^{- \, 2} + \psi^{-} \phi^{-}) \right)  
%\end{array}\right) 
%%\end{split}
%\end{equation}}%
where the function $d(\lambda , \eta )$ is given by
\begin{equation}
\label{Eq: d-min-fun}
d ^- (\lambda , \eta ) = 4 \xi^{- \, 2} - 8 \xi^{-} \lambda \sqrt{\nu^{-\, 2} + \psi^{-} \varphi^{-}} 
+ (4 \lambda ^2 - \eta ^2) (\nu^{- \, 2} + \psi^{-} \varphi^{-}) \>.
\end{equation}
%Notice that the multiplicative factor in \eqref{Eq: K-min} is not essential, we have normalised  the left  reflection matrix $K^-$ so that the following properties are valid 
%\begin{align}
%\label{Eq: Kmin-inv}
%&K ^{-}( - \lambda) \cdot  K ^{-}(\lambda) =  \mathbbm{1}  \>, \\[2ex]
%\label{Eq:}
%&K ^{-}(- \lambda) = \frac{1}{\det K^{-} (\lambda)} \left( {K ^{-}}^2 (\lambda) - \tr K ^{-} (\lambda) \cdot K ^{-} (\lambda) 
%+ \left( \sum_{j =1}^3 \det \widehat{K}_{jj}(\lambda)  \right) \cdot \mathbbm{1} \right) \>,
%\end{align}
%where $\widehat{K}_{jj}(\lambda), j= 1, 2, 3$ are the principal minors of the matrix \eqref{Eq: K-min}. 
As it was shown above \eqref{Eq: Kpl}, the corresponding general solution to the dual reflection equation \eqref{Eq: d-RE} is given by
\begin{equation}
\label{Eq: Kpl-gen}
K^+ (\lambda, \eta, \xi^{+} , \nu^{+} ,  \psi^{+} , \varphi^{+} ) = K^- (-\lambda - \eta,  \eta, \xi^{+} , \nu^{+} ,  \psi^{+} , \varphi^{+}) \>.
\end{equation}
 

However, in order to study the Gaudin model one has to identify the parameters of the reflection matrices at the left and at the right end of the chain, i.e.
\begin{equation}
\label{Eq: Iden-par}
\begin{aligned}
\xi ^-  = \xi^+ = \xi \>,    \qquad  \nu ^- = \nu ^+ = \nu  \>, \qquad  \psi^- = \psi^+= \psi  \>, \qquad  \varphi^-= \varphi^+= \varphi \>.
\end{aligned}
\end{equation}
Usually this is done by imposing the following condition \cite{CAMN1, CAMN2, CiriloAntonio:2014ymf, MS01}
\begin{equation}
\label{normalizationKpl}
\lim_{\eta \to 0}\left(  K^+(\lambda , \eta) \cdot K^{-} (\lambda , \eta) \right)  =  16 \ d (\lambda ) \, d (- \lambda) \, \mathbbm{1} \>,
\end{equation}
here
\begin{equation}
\label{Eq: d-min}
d (\lambda ) =  \left( \xi - \lambda \sqrt{\nu^2 + \psi \varphi} \right) ^2 \>.
%+ (4 \lambda ^2 - \eta ^2) (\nu^{- \, 2} + \psi^{-} \varphi^{-}) 
\end{equation}
This identification is essential due to the to the long range interaction of  the Gaudin model \cite{HikamiKulishWadati92, HikamiKulishWadati92a, Hikami95}.
%
%Due to the normalisation \eqref{Eq: Kmin-inv} of the matrix $K ^{-}(\lambda)$ \eqref{Eq: K-min} we have imposed on the right hand side of the equation above we have the identity matrix. Notice that for some other normalisation of the K-matrix a multiplying scalar function may appear on the right hand side. A particular consequence of the condition \eqref{normalizationKpl} is that the parameters of the reflection matrices on the left and on the right end of the chain are the same. In general this not the case in the study of the open spin chain. However, due , this condition is essential. 
Thus, in the following we will use
\begin{equation}
\label{K-min-notation}
K^-(\lambda, \eta)\equiv K(\lambda, \eta) \>,
\end{equation}
and
\begin{equation}
\label{Eq: def-Kplus}
K^+(\lambda, \eta)= K(-\lambda-\eta, \eta) %= K(-\lambda) -\eta\, M \ \mathbf{with} \ M= \left(\begin{array}{cc} -1 & \psi  \\ \phi  &  1 \end{array}\right)
\>. 
\end{equation}
%Remark that the matrix $M$ obeys $M^2=(1+ \psi \phi) \mathbbm{1}$.

The next step in the derivation of the Gaudin model is the so-called the quasi-classical limit \cite{CiriloAntonio:2014ymf, MS01}. For the step-by-step derivation of the classical Reflection Equation in the $s\ell(2)$ case see\cite{CAMRS} (for the trigonometric case see \cite{MS02}). Here we write the R-matrix as follows
\begin{equation}
\label{Eq: R-exp}
 \frac{1}{\lambda^2}  R (\lambda, \eta)  = \mathbbm{1} + \frac{\eta}{\lambda} \mathbbm{1} - 2 \eta r (\lambda)  + 2 \left( \frac{\eta}{\lambda} \right) ^2 \mathcal{P} \>, 
%\label{Eq:K-exp}    
%K(\lambda, \eta) &=  4 K(\lambda) + 2 \eta \lambda M - \eta ^2 \left( \nu^2 + \psi \varphi \right)  \mathbbm{1} \>,
\end{equation}
where $ r (\lambda)$ is the classical r-matrix 
\begin{equation}
\label{Eq: r-cl}
r (\lambda) = - \frac{\mathcal{P} - 3 \mathcal{K}}{\lambda} \>,
\end{equation}
and the K-matrix in the following form
\begin{equation}
\label{Eq: K-exp}    
K(\lambda, \eta) =  4 K(\lambda) + 2 \eta \lambda M - \eta ^2 \left( \nu^2 + \psi \varphi \right)  \mathbbm{1} \>,
\end{equation}
where $K(\lambda)$ is the classical K-matrix 
\begin{equation}
\label{Eq: K-cl}
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}
and the matrix $M$ is given by
\begin{equation}
\label{Eq: M-mat}
M = \left(\begin{array}{ccc}
\psi \varphi & -\sqrt{2} \nu \psi & - \psi^2 \\[1ex]
-\sqrt{2} \nu \varphi & 2\nu ^2 & \sqrt{2} \nu \varphi \\[1ex]
-\varphi^2 & \sqrt{2} \nu \varphi &  \psi \varphi 
\end{array}\right) \>.
\end{equation}
In a direct substitution of the formulae \eqref{Eq: R-exp} and \eqref{Eq: K-exp} into the Reflection Equation \eqref{Eq: RE}, the terms in the zero order in 
$\eta$ lead to some obvious identity satisfied by the classical K-matrix \eqref{Eq: K-cl}. However, the first order terms in $\eta$, 
besides some obvious identities for the matrices $K(\lambda)$ and $M$, yield 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}


Moreover the substitution of the formula \eqref{Eq: def-Kplus} into the dual reflection equation \eqref{Eq: d-RE} yields
\begin{equation}
\label{Eq: d-RE-new}
\begin{split}
R_{12}(-\lambda + \mu , \eta) K _1( -\lambda -\eta , \eta) R_{21}(-\lambda - \mu -2\eta , \eta) K _2(-\mu-\eta , \eta) = \\
= K _2(-\mu-\eta , \eta) R_{12}(-\lambda - \mu -2\eta , \eta) K _1( -\lambda -\eta , \eta) R_{21}(-\lambda + \mu , \eta) \>.
\end{split}
\end{equation}  
As expected we have obtained the reflection equation \eqref{Eq: RE} with the shifted arguments. To study the quasi-classical limit in this case, we use the formulae \eqref{Eq: R-exp} and \eqref{Eq: K-exp} with the shifted arguments $\lambda \to -\lambda - \eta$ and $\mu \to -\mu - \eta$. By following the same steps we have done in deriving the classical reflection equation \eqref{Eq: c-RE}, we obtain the same equation but with shifted arguments 
\begin{equation}
\label{Eq: c-RE-shif}
\begin{split}
r_{12}(-\lambda + \mu) K _1(-\lambda - \eta) K_2(-\mu - \eta) + K _1(-\lambda - \eta) r_{21}(-\lambda - \mu -2\eta) K_2(-\mu - \eta) = \\
= K_2(-\mu - \eta) r_{12}(-\lambda - \mu -2\eta) K _1(-\lambda - \eta) +  K_2(-\mu - \eta) K _1(-\lambda - \eta) r_{21}(-\lambda + \mu) \>.
\end{split}
\end{equation}

This shows that, due to the to the long range interaction of  the Gaudin model \cite{HikamiKulishWadati92, HikamiKulishWadati92a, Hikami95}, the integrability of the system is guaranteed by the classical Yang-Baxter equation 
\begin{equation}
\label{Eq: cl-YBE}
 [r_{12}(\lambda - \mu) , r_{13} (\lambda)] +  [r_{12}(\lambda - \mu) , r_{23}(\mu)] + [r_{13} (\lambda) , r_{23} (\mu) ] = 0 ,
\end{equation}
and the classical reflection equation \eqref{Eq: c-RE}. Moreover, by introducing the non-unitary, classical r-matrix \cite{CAMRS}
\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}
the two equations can be combined into the generalized classical Yang-Baxter equation \cite{CAMRS}
\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}

\vspace{1cm}

To address the third question posed by the referee (and following our ideas presented in \cite{CAMRS}), we first observe that in the study of the corresponding inhomogeneous  Heisenberg spin chain with N sites one can use the following Lax operator
\begin{equation}
\label{Eq: Lax-op-sc}
\mathbb{L}_{0m} (\lambda) = \frac{1}{\lambda (\lambda+\eta)}  R_{0m} (\lambda, \eta) \>.
\end{equation}
Due to the Yang-Baxter equation  \eqref{Eq: YBE} the Lax operator satisfies the RLL-relations
\begin{equation}
\label{Eq: RLL}
R_{00'} (\lambda -\mu) \mathbb{L}_{0m} (\lambda) \mathbb{L}_{0'm} (\mu) = \mathbb{L}_{0'm} (\mu) \mathbb{L}_{0m} (\lambda) R_{00'} (\lambda -\mu) \>,
\end{equation}
as well as the following important identity 
\begin{equation}
\label{Eq: LL}
\mathbb{L}_{0m} (\lambda) \mathbb{L}_{0m} (-\lambda) = \left( 1 - \left(\frac{2\eta}{\lambda} \right) ^2 \right) \mathbbm{1}_{0} \>.
\end{equation} 
Hence the corresponding monodromy operator reads
\begin{equation}
\label{Eq: T-op}
T (\lambda) = \mathbb{L}_{0N} (\lambda - \alpha _N) \cdots \mathbb{L}_{01} (\lambda - \alpha _1) \>.
\end{equation}
For simplicity we have omitted the dependence on the quasi-classical parameter $\eta$ and the inhomogeneous parameters $\{ \alpha _j , j = 1, . . . , N \}$. Notice that $T (\lambda)$ is a three-by-three matrix acting in the auxiliary space $V_0 = \mathbb{C}^3$, whose entries are operators acting in the Hilbert space of the system. From RLL-relations \eqref{Eq: RLL} it follows that the monodromy matrix satisfies the RTT-relations
\begin{equation}
\label{Eq: RTT}
R_{00'} (\lambda -\mu) T_{0} (\lambda) T_{0'} (\mu) = T_{0'} (\mu) T_{0} (\lambda) R_{00'} (\lambda -\mu) \>.
\end{equation}

For our study of the corresponding Gaudin model it is instructive to consider the expansion of the monodromy matrix \eqref{Eq: T-op} with respect to the quasi-classical parameter $\eta$
\begin{equation}
\label{Eq: T-eta-exp}
T_0 (\lambda) = \mathbbm{1}_0 + 2 \eta L _0 (\lambda) + \mathcal{O}(\eta ^2) \>,
\end{equation}
where
\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}
here we us the same notation as in the manuscript. 


Moreover the identity \eqref{Eq: LL} implies that
\begin{equation}
\label{Eq: LL}
\mathbb{L}_{0m} (\lambda - \alpha _m) \mathbb{L}_{0m} (-\lambda + \alpha _m) = \left( 1 - \frac{4\eta ^2 }{(\lambda - \alpha _m)^2 } \right) \mathbbm{1}_{0} \>.
\end{equation} 
Thus the equation above and the RLL-relations \eqref{Eq: RLL} imply that the RTT-relations  \eqref{Eq: RTT} can be recasted as follows
\begin{align}
\label{Eq: RTT-til}
\widetilde{T}_{0'} (\mu) R_{00'} (\lambda + \mu) T_{0} (\lambda) &=  T_{0} (\lambda)  R_{00'} (\lambda + \mu)  \widetilde{T}_{0'} (\mu) \>,\\%[1ex]
\label{Eq: RT-tilT-til}
\widetilde{T}_{0} (\lambda) \widetilde{T}_{0'} (\mu)  R_{00'} (\mu - \lambda)   &=   R_{00'} (\mu - \lambda) \widetilde{T}_{0'} (\mu) \widetilde{T}_{0} (\lambda) \>.
\end{align}
where
\begin{equation}
\label{Eq: T-til-op}
\widetilde{T}_{0} (\lambda) = \mathbb{L}_{01} (\lambda + \alpha _1) \cdots \mathbb{L}_{0N} (\lambda + \alpha _N) \>.
\end{equation}
 
 
Therefore the Sklyanin monodromy $\mathcal{T} (\lambda)$ of the $SO(3)$ inhomogeneous spin chain with non-periodic boundary consists 
of the operators $T (\lambda)$ \eqref{Eq: T-op} and $\widetilde{T}_{0} (\lambda)$ \eqref{Eq: T-til-op} and the reflection matrix $K ^{-}(\lambda)$ \eqref{Eq: K-min}
\begin{equation}
\label{Eq: Skly-T}
\mathcal{T}_0 (\lambda) = T_0 (\lambda) K _0^{-}(\lambda, \eta) \widetilde{T}_{0} (\lambda) \>.
\end{equation}
It follows from the reflection equation \eqref{Eq: RE} and the RTT-relations \eqref{Eq: RTT}, \eqref{Eq: RTT-til} , \eqref{Eq: RT-tilT-til} that the exchange relations of the Sklyanin monodromy $\mathcal{T} (\lambda)$ in $V_0\otimes V_{0'}$ are
\begin{equation}
\label{Eq: REA}
R_{00'} (\lambda -\mu) \mathcal{T}_0 (\lambda) R_{0'0} (\lambda + \mu) \mathcal{T}_{0'} (\mu) = \mathcal{T}_{0'} (\mu)  R_{0'0} (\lambda + \mu) \mathcal{T}_0 (\lambda) R_{0'0} (\lambda -\mu) \>.
\end{equation}
The open Heisenberg spin chain  transfer matrix is given by the trace of $\mathcal{T}_0 (\lambda)$ over the auxiliary space $V_0$ with an extra reflection matrix $K^+(\lambda, \eta)$ \cite{Sklyanin88},
\begin{equation}
\label{Eq: tm-H-ch}
t (\lambda) = \tr _{0} \left( K^+(\lambda, \eta) \, \mathcal{T}_0 (\lambda) \right) \>.
\end{equation}
The commutativity of the transfer matrix for different values of the spectral parameter
\begin{equation}
\label{Eq: t-t-com}
[t (\lambda), t (\mu) ] = 0 \>,
\end{equation}
is guaranteed by the dual reflection equation \eqref{Eq: d-RE} and the exchange relations \eqref{Eq: REA} of the monodromy matrix $\mathcal{T}_0 (\lambda)$. Moreover, a central element of the reflection equation algebra \eqref{Eq: REA} is given by
\begin{equation}
\label{Eq: Delta}
\Delta [\mathcal{T}_0 (\lambda)] = \frac{1}{6\eta^2} \tr _{00'} R _{00'}(-\eta,\eta) \mathcal{T}_0 (\lambda - \eta/2) R _{00'}(2\lambda,\eta) \mathcal{T}_0 (\lambda + \eta/2) \>.
\end{equation}

The expansion of the Sklyanin monodromy $\mathcal{T} (\lambda)$ in powers of the quasi-classical parameter $\eta$ reads
\begin{equation}
\label{Eq: Skly-T-exp}
\mathcal{T}_0 (\lambda) = 4 K_0(\lambda) + 2 \eta \lambda M_0 + 8 \eta \mathcal{L} _0 (\lambda) K_0 (\lambda) + \frac{\eta ^2}{2} \frac{d ^2 \mathcal{T}_0 (\lambda)}{d\lambda ^2} \Big{|}_{\lambda=0} + \mathcal{O}(\eta ^3) \>,
\end{equation}
where
\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}
By substituting \eqref{Eq: R-exp} and \eqref{Eq: Skly-T-exp} into \eqref{Eq: REA} we can confirm that the zero and first orders in $\eta$ are identically satisfied. The relations we seek follow from the terms of the second order in $\eta$. When the terms containing the second order derivatives of $\mathcal{T} (\lambda)$ are eliminated then, after some calculations, it can be shown that the remaining terms can be recasted in the following form
\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}
here the non-unitary classical r-matrix $r_{00'}^{K} (\lambda , \mu)$ is given in \eqref{Eq: rK-def}. This linear bracket is obviously anti-symmetric and it obeys the Jacobi identity because the $r$-matrix $r_{00'}^{K} (\lambda , \mu)$ \eqref{Eq: rK-def} satisfies the generalized classical Yang-Baxter equation \eqref{Eq: g-c-YBE}.



Moreover, the generating function of the Gaudin model 
\begin{equation}
\label{Eq: tau}
\tau (\lambda) = \frac{1}{2} \, \mathrm{tr}_0 \left( \mathcal{L}_0 ^2(\lambda) \right) \>.
\end{equation}
can be obtained by the $\eta$ expansion of a linear combination of the transfer matrix $t (\lambda)$ \eqref{Eq: tm-H-ch} of the corresponding
Heisenberg spin and the central element $\Delta [\mathcal{T}_0 (\lambda)]$ \eqref{Eq: Delta} of the reflection equation algebra \eqref{Eq: REA}.

As we have shown above, our derivation of the Gaudin model with boundary terms can be resumed to the following five formulae. The non-unitary classical r-matrix $r_{00'}^{K} (\lambda , \mu)$ \eqref{Eq: rK-def}, together with the the generalized classical Yang-Baxter equation \eqref{Eq: g-c-YBE},
the Lax operator $\mathcal{L} _{0} (\lambda)$ \eqref{Eq: cal-L-def} of corresponding generalized Gaudin algebra, jointly with the linear bracket \eqref{Eq: M-br-rk}, and finally the generating function $\tau (\lambda)$  \eqref{Eq: tau}. In this sense we can say the this is the unique structure which defines the system at hand.  


\vspace{3cm}
%
%
\begin{thebibliography}{99}
\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.
\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{Maillet86} J.-M. Maillet, \textsl{New integrable canonical structures in two-dimensional models}, Nuclear Physics \textbf{B269} (1986) 54--76.
\bibitem{Sklyanin88} E. ~K. ~Sklyanin, \textsl{Boundary conditions for integrable quantum systems}, J. Phys. A: Math. Gen. \textbf{21} (1988) 2375--2389.
\bibitem{KMN3} P.~P.~Kulish,  N. ~Manojlovi\'c and Z. ~Nagy, \emph{Jordanian deformation of the open XXX spin chain}, (in Russian) Teoreticheskaya i Matematicheskaya Fizika \textbf{Vol. 163} No. 2 (2010) 288-298; translation in Theoretical and  Mathematical Physics \textbf{Vol. 163} No. 2 (2010) 644-652;  {\blue doi: 10.1007/s11232-010-0047-x;}  {\red \texttt{arXiv:0911.5592}.} 
\bibitem{KulishSklyanin82} P. ~P. ~Kulish and ~E. ~K. ~Sklyanin, \textsl{Quantum spectral transform method. Recent developments}, Lect. Notes Phys. \textbf{151} (1982), 61--119.
%\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.~Manojlovi\'c 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{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
\bibitem{CAMN1} N. ~Cirilo Ant\'onio, N. ~Manojlovi\'c and Z. ~Nagy, \emph{Trigonometric $s\ell(2)$ Gaudin model with boundary terms}, Reviews in Mathematical Physics \textbf{Vol. 25} No. 10 (2013) 1343004 (14 pages); {\blue doi: 10.1142/S0129055X13430046;}  {\red \texttt{arXiv:1303.2481}.}
\bibitem{CAMN2} N. ~Cirilo Ant\'onio, N. ~Manojlovi\'c and Z. ~Nagy, \emph{Jordanian deformation of the open $s\ell(2)$ Gaudin model}, Teoreticheskaya i Matematicheskaya Fizika, \textbf{Vol. 179}, No. 1 (2014) 9--101; translation in Theoretical and Mathematical Physics, \textbf{Vol. 179}, No. 1 (2014) 462--471;  {\blue doi: 10.1007/s11232-014-0155-0;} {\red \texttt{arXiv:1304.6918}.}
\bibitem{MS01} N.~Manojlovi\'c and I.~Salom, \ \emph{Algebraic Bethe ansatz for the XXZ Heisenberg spin chain with triangular boundaries and the corresponding Gaudin model}, Nuclear Physics \textbf{B 923} (2017) 73-106; {\blue doi:10.1016/j.nuclphysb.2017.07.017;} {\red \texttt{arXiv:1705.02235}.}

\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{Hikami95} K.~Hikami, \textsl{Gaudin magnet with boundary and generalized Knizhnik-Zamolodchikov equation}, J. Phys. A Math. Gen. \textbf{28} (1995) 4997--5007.
 
\bibitem{CAMRS} 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}, Nuclear Physics \textbf{B 893} (2015) 305-331; {\blue doi:10.1016/j.nuclphysb.2015.02.011;} {\red \texttt{arXiv:1412.1396}.}
\bibitem{MS02} N.~Manojlovi\'c and I.~Salom, \emph{Algebraic Bethe ansatz for the trigonometric $s\ell(2)$ Gaudin model with triangular boundary}, Symmetry 12 (2020) 352; {\blue doi:10.3390/sym12030352}; {\red \texttt{arXiv:1709.06419}.}
\end{thebibliography}









\end{document}
