  Motivated by recent experiments [Nat. Phys. $\textbf{16}$, 1227 (2020)], we
present here a theoretical study of the DC Josephson effect in a system
comprising two magnetic impurities coupled to their respective superconducting
electrodes and which exhibit Yu-Shiba-Rusinov (YSR) states. We make use of a
mean-field Anderson model with broken spin symmetry to compute the supercurrent
in this system for an arbitrary range of parameters (coupling between the
impurities, orientation of the impurity spins, etc.). We predict a variety of
physical phenomena such as (i) the occurrence of multiple $0-\pi$ transitions
in the regime of weak coupling that can be induced by changing the energy of
the YSR states or the temperature; (ii) the critical current strongly depends
on the relative orientation of the impurity spins and it is maximized when the
spins are either parallel or antiparallel, depending on the ground state of the
impurities; and (iii) upon increasing the coupling between impurities, triplet
superconductivity is generated in the system and it is manifested in a highly
nonsinusoidal current-phase relation. In principle, these predictions can be
tested experimentally with the existing realization of this system and the main
lessons of this work are of great relevance for the field of superconducting
spintronics.
