As another manifestation of the integrability of the model (REF ), it admits a Lax representation, i.e., a pair of matrix-valued functions \(L, A\,:\, \mathbb {C}^{2N}\rightarrow \mathrm {Mat}(2N, \mathbb {C})\)  such that the Hamiltonian dynamics takes the form \(\frac{d}{dt} L=[L, A]\) . An immediate corollary is that the quantities \(\mathrm {tr}(L^k)\) , as well as the eigenvalues of \(L\) , are constants of motion, which means that \(L\)  remains isospectral for all \(t\) . Originally, Inozemtsev constructed in [1]} a Lax pair of size \(3N\times 3N\)  (see also [2]}); other Lax pairs of smaller size have since been found [3]}, [4]}. We will use the Lax matrix of size \(2N\times 2N\)  from [4]}. To write it down, we need the functions \({\sigma }_\alpha (x)\)  and \(v_\alpha (x) := \sum _{r=0}^{3} g_r {\sigma }_{2 \alpha }^r(x)\)  whose definition and basic properties are given in the Appendix.
We have:
\(L & = \sum _{i=1}^N \left(p_iE_{i,i}-p_iE_{i+N, i+N}+v_\alpha (q_i)E_{i, i+N}+v_\alpha (-q_i)E_{i+N,i}\right)\\& + g\sum _{i\ne j}^N \left(\sigma _{\alpha }(q_{ij})E_{i,j}+\sigma _{\alpha }(q_{ij}^+)E_{i, j+N}+\sigma _{\alpha }(-q_{ij}^+)E_{i+N, j}+\sigma _{\alpha }(-q_{ij})E_{i+N, j+N}\right)\,, \nonumber \) 
