Exact fermionic matrix elements from the Bethe-ansatz solution of the Richardson-Gaudin model

We introduce an effective Hamiltonian description of mesoscopic superconducting systems based on the Bethe-ansatz solution of the Richardson--Gaudin model. Our central result is a compact determinant formula for the matrix elements of fermionic creation and annihilation operators, expressed entirely in terms of eigenvalue-based variables. This formula allows single-electron tunneling terms to be represented directly in the eigenbasis of the Richardson--Gaudin model, extending the eigenvalue-based framework beyond overlaps and previously known form factors. As a benchmark, we apply the method to a superconducting Anderson impurity model comprising a quantum dot coupled to an ultrasmall superconducting grain. The resulting low-energy description, formulated naturally in the canonical ensemble, captures even-odd effects, subgap excitations, and quantum-dot physics. The same framework gives access to observables not previously computed for such grain-dot systems in the canonical ensemble, such as the fixed-particle-number optical conductivity and flux-tunable spectra.