We present a numerical approach for computing the properties of quantum dots (QDs) coupled to superconducting (SC) regions with finite charging energy Ec (e.g., SC grains, thin epitaxial SC layers on nanowires, etc.). It is based on the matrix product operator (MPO) representation of the Hamiltonian in terms of small 9x9 matrices. The low-lying excitations can be computed using the density matrix renormalization group (DMRG) and the (imaginary) time evolution using the time-dependent variational principle (TDVP), providing information on the dynamic response of the system (e.g., dynamical charge susceptibility of the QD). The method treats pairing interaction, electron-electron repulsion and the Kondo exchange interaction induced by the QD-SC hybridization on equal footing, and it is essentially exact.

For the QD-SC problem, we examined the transition from the regime dominated by Yu-Shiba-Rusinov (YSR) physics to the regime governed by the Coulomb blockade, with the nature of the low-lying excitations continuously evolving between the two limits [1]. For Ec>D, the method also allows investigating the case of odd-occupancy of the SC island [1], where a new type of subgap excitation is found with no counterpart in the Ec<D range. These predictions have been confirmed in recent experiments [2].

[1] Pavešić, L., Bauernfeind D., Žitko, R. (2021). Yu-Shiba-Rusinov states in superconducting islands with finite charging energy, arXiv:2101.10168.

[2] Estrada Saldaña, J. C., Vekris, A., Pavešić, L., Krogstrup, P., Žitko, R., Grove-Rasmussen, K., Nygård, J. (2021). Bias asymmetric subgap states mimicking Majorana signatures. arXiv, 2101.10794.