2D CFT and efficient Bethe ansatz for exactly solvable Richardson-Gaudin models

This work inaugurates a series of complementary studies on Richardson-Gaudin integrable models. We begin by reviewing the foundations of classical and quantum integrability, recalling the algebraic Bethe ansatz solution of the Richardson (reduced BCS) and Gaudin (central spin) models, and presenting a proof of their integrability based on the Knizhnik-Zamolodchikov equations and their generalizations to perturbed affine conformal blocks. Building on this foundation, we then describe an alternative CFT-based formulation. In this approach, the Bethe ansatz equations for these exactly solvable models are embedded within two-dimensional Virasoro CFT via irregular, degenerate conformal blocks. To probe new formulations within the Richardson-Gaudin class, we develop a high-performance numerical solver. The Bethe roots are encoded in the Baxter polynomial, with initial estimates obtained from a secular matrix eigenproblem and subsequently refined using a deflation-assisted hybrid Newton-Raphson/Laguerre algorithm. The solver proves effective in practical applications: when applied to picket-fence, harmonic oscillator, and hydrogen-like spectra, it accurately reproduces known rapidity trajectories and reveals consistent merging and branching patterns of arcs in the complex rapidity plane. We also explain how to generalize our computational approach to finite temperatures, allowing us to calculate temperature-dependent pairing energies and other thermodynamic observables directly within the discrete Richardson model. We propose an application of the solver to Gaudin-type Bethe equations, which emerge in the classical (large central charge) limit of Virasoro conformal blocks. We conclude by outlining future directions: direct minimization of the Yang-Yang function as an alternative root-finding strategy; revisiting time-dependent extensions; and ... .