  Very recent experiments have reported the tunneling between Yu-Shiba-Rusinov
(YSR) bound states at the atomic scale. These experiments have been realized
with the help of a scanning tunneling microscope where a superconducting tip is
functionalized with a magnetic impurity and is used to probe another magnetic
impurity deposited on a superconducting substrate. In this way it has become
possible to study for the first time the spin-dependent transport between
individual superconducting bound states. Motivated by these experiments, we
present here a comprehensive theoretical study of the tunneling processes
between YSR bound states in a system in which two magnetic impurities are
coupled to superconducting leads. Our theory is based on a combination of an
Anderson model with broken spin degeneracy to describe the impurities and
nonequilibrium Green's function techniques to compute the current-voltage
characteristics. This combination allows us to describe the spin-dependent
transport for an arbitrary strength of the tunnel coupling between the
impurities. We first focus on the tunnel regime and show that our theory
naturally explains the experimental observations of the appearance of current
peaks in the subgap region due to both the direct and thermal tunneling between
the YSR states in both impurities. Then, we study in detail the case of
junctions with increasing transparency, which has not been experimentally
explored yet, and predict the occurrence of a large variety of (multiple)
Andreev reflections mediated by YSR states that give rise to a very rich
structure in the subgap current. In particular, we predict the occurrence of
multiple Andreev reflections that involve YSR states in different impurities.
These processes have no analogue in single-impurity junctions and they are
manifested as current peaks with negative differential conductance for subgap
voltages.
