Variational solution of the superconducting Anderson impurity model and the band-edge singularity phenomena

We propose a set of variational wavefunctions for the sub-gap spin-doublet and spin-singlet eigenstates of the particle-hole symmetric superconducting Anderson impurity model. The wavefunctions include up to two Bogoliubov quasiparticles in the continuum which is necessary to correctly capture the weak-coupling asymptotics in all parameter regimes. The eigenvalue problems reduce to solving transcendental equations. We investigate how the lowest singlet state evolves with increasing charge repulsion $U$, transitioning from a proximitized state (a superposition of empty and doubly occupied impurity orbitals, corresponding to an Andreev bound state) to a local moment that is Kondo screened by Bogoliubov quasiparticles (Yu-Shiba-Rusinov state). This change occurs for $U = 2\Delta$, where $\Delta$ is the BCS gap. At this point, the band-edge effects make the eigenenergy scale in a singular way as $\Gamma^{2/3}$, where $\Gamma$ is the hybridization strength. Away from this special point, regular $\Gamma$-linear behavior is recovered, but only for $\Gamma \lesssim (U/2-\Delta)^2/\Delta$. The singular behavior thus extends over a broad range of parameters, including those relevant for some quantum devices in current use. The singular state is an equal-superposition state with maximal fluctuations between the local impurity charge configurations. Accurately capturing the band-edge singularity requires a continuum model, and it cannot be correctly described by discrete (truncated) models such as the zero-bandwidth approximation or the superconducting atomic limit. We determine the region of parameter space where the second spin-singlet state exists: in addition to the whole $U<2\Delta$ ABS region, it also includes a small part of the $U>2\Delta$ YSR region for finite values of $\Gamma$, as long as some ABS wavefunction component is admixed.