The XXZ spin chain Bethe Ansatz equations are a complicated system of rational function equation:

\[ \left(\frac{\lambda_j + i/2}{\lambda_j - i/2} \right)^N = \prod_{l=1, l \neq j}^M \frac{\lambda_j - \lambda_l + i}{\lambda_j - \lambda_l - i} \]

Generically, since there are M equations and M unknowns generically we could say the solutions in $\lambda$ are a discrete set of points.

Is there any other logic to their roots besides the fact that they solve these equations ?

Spin chains have to do with the raising and lowering operators in the 2-dimensional representation of $SU(2)$ : \[\left.\begin{array}{cccc} S^+ : & |\uparrow \rangle & \mapsto & 0 \\ & |\downarrow \rangle & \mapsto & |\uparrow \rangle \\ \hline S^- : & |\uparrow \rangle & \mapsto & |\downarrow \rangle \\ & |\downarrow \rangle & \mapsto & 0 \\ \hline \\ S^z : & |\uparrow \rangle & \mapsto & \tfrac{1}{2}|\uparrow \rangle \\ & |\downarrow \rangle & \mapsto & \tfrac{1}{2}|\downarrow \rangle \end{array}\right. \]

The Hamiltonian acts on a tensor product of $SU(2)$ representations $V^{\otimes M}$.

\[ \mathcal{H} = - \frac{1}{2} \sum_{n=1}^M \bigg[ \sigma_n^x \sigma_{n+1}^x + \sigma_n^y \sigma_{n+1}^y + \Delta \sigma_n^z \sigma_{n+1}^z\bigg] \]

Various values of $\Delta$ have different interpretations. Since $[\mathcal{H}, S^z] = 0$ we can diagonalize in terms of the $\uparrow, \downarrow$ states. The Bethe Ansatz equations are the eigenvalue equations for $\mathcal{H}$ in this basis.

None of this tells you how to *solve* the eigenvalue equations, here they are.