  In this work we present numerical results for physical quantities in the
steady-state obtained after a variety of product-states initial conditions are
evolved unitarily, driven by the dynamics of quantum integrable models of the
rational (XXX) Richardson-Gaudin family, which includes notably Tavis-Cummings
models. The problem of interest here is one where a completely inhomogeneous
ensemble of two-level systems (spins-1/2) are coupled to a single bosonic mode.
  The long-time averaged magnetisation along the z-axis as well as the bosonic
occupation are evaluated in the diagonal ensemble by performing the complete
sum over the full Hilbert space for small system sizes. These numerically exact
results are independent of any particular choice of Hamiltonian and therefore
describe general results valid for any member of this class of quantum
integrable models built out of the same underlying conserved quantities.
  The collection of numerical results obtained can be qualitatively understood
by a relaxation process for which, at infinitely strong coupling, every initial
state will relax to a common state where each spin is in a maximally coherent
superposition of its $\left|\uparrow\right>$ and $\left|\downarrow\right>$
states, i.e. they are in-plane polarised, and consequently the bosonic mode is
also in a maximally coherent superposition of different occupation number
states. This bosonic coherence being a feature of a superradiant state, we
shall loosely use the term superradiant steady-state to describe it.
  A finite value of the coupling between the spins and the bosonic mode then
leads to a long-time limit steady-state whose properties are qualitatively
captured by a simple "dynamical" vision in which the coupling strength $V$
plays the role of a time $t_V$ at which this "relaxation process" towards the
common strong coupling superradiant steady-state is interrupted.
