Quadratic operator relations and Bethe equations for spin-1/2 Richardson-Gaudin models

In this work we demonstrate how one can, in a generic approach, derive a set of $N$ simple quadratic Bethe equations for integrable Richardson-Gaudin (RG) models built out of $N$ spins-1/2. These equations depend only on the $N$ eigenvalues of the various conserved charges so that any solution of these equations defines, indirectly through the corresponding set of eigenvalues, one particular eigenstate. The proposed construction covers the full class of integrable RG models of the XYZ (including the subclasses of XXZ and XXX models) type realised in terms of spins-1/2, coupled with one another through $\sigma_i^x \sigma_j^x$, $\sigma_i^y \sigma_j^y$, $\sigma_i^z \sigma_j^z$ terms, including, as well, magnetic field-like terms linear in the Pauli matrices. The approach exclusively requires integrability, defined here only by the requirement that $N$ conserved charges $R_i$ (with $i = 1,2 \dots N$) such that $\left[R_i,R_j\right] =0 \ (\forall \ i,j)$ exist . The result is therefore valid, and equally simple, for models with or without $U(1)$ symmetry, with or without a properly defined pseudo-vacuum as well as for models with non-skew symmetric couplings.