What is the $XXX_s$ Hamiltonian in terms of $\vec{S}_i \cdot \vec{S}_{i+1}$? Faddeev, Takhtajan, and others united and discovered many integrable models through the Algebraic Bethe Ansatz. For example, the integrable spin-1/2 Heisenberg model
$$H_{1/2} = \sum_{i=1}^L \vec{S}_i \cdot \vec{S}_{i+1}$$
became just the first model of the family of spin-$s$ $XXX_s$ models, which includes the integrable spin-1 $XXX_1$ (Lai-Sutherland model) of
$$H_{1} = \sum_{i=1}^L \vec{S}_i \cdot \vec{S}_{i+1} \pm (\vec{S}_i \cdot \vec{S}_{i+1})^2$$ where I believe both signs appear in the literature (I think the models are unitarily equivalent, perhaps).
I am trying to learn the explicit forms of the $XXX_s$ models in terms of the quantity $\vec{S}_i \cdot \vec{S}_{i+1}$. In googling it, many references cite Faddeev's lecture notes on the Algebraic Bethe Ansatz as giving exactly my desired form.
Faddeev gives the following for the two-site terms appearing in $H_s$, where I've combined and transcribed a few of his equations hopefully correctly. Note that in practice he strips out additive constants and overall multiplicative constants, and I personally doubt the necessity of $\gamma$:
$$(h_{s})_{i,i+1} = \sum_{j=1}^{2s} \left[ \left( \sum_{k=1}^j \frac{1}{k} - \gamma   \right) \prod_{l=0}^{2s} \frac{\vec{S}_i \cdot \vec{S}_{i+1} - \frac{1}{2}(l(l+1)-2s(s+1))}{\frac{1}{2}(j(j+1)-2s(s+1)) - \frac{1}{2}(l(l+1)-2s(s+1))} \right ]$$
First, there needs to be an additional $l \neq j$ on the product for the product to be well-defined. However, I am unable to use the above formula to recover $H_1$ above, even though Faddeev claims it gives the $H_1$ above with a minus sign.

What is the correct formula for the two-site interactions in the $XXX_{s}$ model in terms of $\vec{S}_i \cdot \vec{S}_{i+1}$? It's possible that I'm just summing wrong, so please note if you get the result for $H_1$ (with the minus sign) up to additive constants and overall multiplicative constants.
 A: I believe I've found the answer. Faddeev is using the Lagrange interpolation formula, which returns a degree-$n$ polynomial from data points $(x_0, y_0), (x_1, y_1),...,(x_n, y_n)$ via $$\sum_{k=0}^n y_k \prod_{j=0, j\neq k}^n \frac{x-x_j}{x_j - x_k}.$$
Note that the lower and upper bounds of the sum and product match in the Lagrange interpolation formula. However, Faddeev's bounds don't match, and it's easy to check that $j=0$ is a perfectly safe contribution to the polynomial ($j=-1$, on the other hand, would suffer a poorly behaved digamma function arising earlier in the problem). Thus Faddeev needed the following:
$$(h_{s})_{i,i+1} = \sum_{j=\color{red}{0}}^{2s} \left[ \left( \sum_{k=1}^j \frac{1}{k} - \gamma   \right) \prod_{l=0, \color{red}{l \neq j}}^{2s} \frac{\vec{S}_i \cdot \vec{S}_{i+1} - \frac{1}{2}(l(l+1)-2s(s+1))}{\frac{1}{2}(j(j+1)-2s(s+1)) - \frac{1}{2}(l(l+1)-2s(s+1))} \right ]$$
This suffices to find the Hamiltonian, as it can only be at most a $2s$th degree polynomial in $\vec{S} \cdot \vec{S}$, as higher powers of $\vec{S} \cdot \vec{S}$ are not linearly independent from the $0,1,...,2S$ powers. As an aside, the form above makes it clear that $\gamma$ only adds a constant term.
This yields the following Hamiltonians, after stripping out constant terms and dividing by overall multiplicative constants:
$$H_{1/2} = \sum_i \vec{S}_i \cdot \vec{S}_{i+1}$$
$$H_{1} = \sum_i \vec{S}_i \cdot \vec{S}_{i+1} - (\vec{S}_i \cdot \vec{S}_{i+1})^2$$
$$H_{3/2} = \sum_i \vec{S}_i \cdot \vec{S}_{i+1} - \frac{8}{27}(\vec{S}_i \cdot \vec{S}_{i+1})^2- \frac{16}{27}(\vec{S}_i \cdot \vec{S}_{i+1})^3$$
$$H_{2} = \sum_i \vec{S}_i \cdot \vec{S}_{i+1} + \frac{43}{234}(\vec{S}_i \cdot \vec{S}_{i+1})^2- \frac{5}{117}(\vec{S}_i \cdot \vec{S}_{i+1})^3- \frac{1}{78}(\vec{S}_i \cdot \vec{S}_{i+1})^4$$
