Effects of electron-electron interactions in the Yu-Shiba-Rusinov lattice model

Magnetic impurities in conventional and non-magnetic impurities in unconventional superconductors give rise to Yu-Shiba-Rusinov bound states, which have been studied extensively both theoretically and experimentally [1, 2, 3, 4, 5, 6]. An important property of these states in two-dimensional systems is that they spatially extent over long distances [5, 6, 7, 8, 9]. Therefore, in the presence of impurity chains and lattices the coherent overlapping impurity states can give rise to a long-range hopping models and topological superconductivity [10, 11, 12]. The one-dimensional magnetic adatom chains are theoretically predicted to support Majorana zero modes localized at the end of the chain [13, 14, 15, 16, 10, 11, 17] and promising signatures have been experimentally observed [18, 19, 20, 21, 22]. These observations have attracted significant interest in the quantum computing community because the non-Abelian braiding statistics of the Majorana zero modes. In particular, they can be utilized in topological quantum computing and the quantum information stored in the Majorana qubits is protected from local sources of noise [23, 24]. Although most of the braiding proposals have been developed for Majorana nanowires [24], different methods for manipulating the quantum information in Yu-Shiba-Rusinov chains are also currently being explored [25, 26, 27].

The case of a two-dimensional Yu-Shiba-Rusinov lattice is theoretically even more interesting [28, 29, 12, 30, 31, 32, 33, 34, 35, 36]. In this case the hybridization of the Yu-Shiba-Rusinov states gives rise a rich phase diagram containing large number of different topological phases [12, 30, 31, 32]. These topological phases are described by an integer-valued Chern number $C$ [37], which determines the number of chiral Majorana edge modes (see Fig. 1) and value of the quantized thermal conductance [38, 39]. So far the quantization of the thermal conductance in superconductors has remained elusive, but recent experiments show promising signatures of the Majorana edge modes in the local density of states (LDOS) data measured via the scanning tunneling spectroscopy techniques [40, 41, 42, 43], and unbiased methods for identifying the Chern number from LDOS are currently under development [44, 45, 46, 47, 48]. Therefore, there are reasons to be optimistic that the topological phase diagram of the Yu-Shiba-Rusinov lattice system can be probed also experimentally in the near future.

The rich topological phase diagram of the Yu-Shiba-Rusinov lattice model [12, 30, 31, 32] calls for a detailed theoretical analysis of the different factors that may play an important role in the selection of the topological phases, which are actually realized in the experiments. So far most of the studies have mainly focused on extrinsic factors and it is known for example that in the presence of strong disorder only the $C=0$ and $C=1$ phases survive [49], out of the dozens of topologically distinct phases appearing in a clean system [12, 32]. Here, we study an important intrinsic factor which is present also in the highest quality samples: the effect of electron-electron interactions. We may expect that they lead to competing effects. On one hand, the interactions favor states with large energy gaps in order to lower the free energy in the many-particle systems. On the other hand, the phases with large Chern numbers are most fragile to the effects of the perturbations because they arise from a complex interplay of the superconductivity, magnetism and large spatial extent of the Yu-Shiba-Rusinov states. In particular, the increase of the energy gap leads to a shorter coherence length decreasing the long-range coupling between the Shiba states which is essential for realizing large Chern numbers. Therefore, the interactions could favor the phases with small Chern numbers, where it is easier to open sizeable energy gaps. But, in order that the interactions could result in a change of the Chern number from a large to a small value the system, in the absence of first order transitions, has to undergo a series of energy gap closings. Thus, the interactions can also lead to lowering of the energy gaps. We show that these competing tendencies in this type of complicated system supporting dozens of topologically distinct phases lead to rather surprising statistical effects. By sampling a large number of realizations of Yu-Shiba-Rusinov lattice models with different model parameters, we show that statistically the interactions have practically no effect on the realized Chern numbers and typical magnitudes of the topological gaps. However, the interactions substantially enhance the likelihood of the rare realizations of large topological gaps in the tails of the energy gap distribution in comparison to the non-interacting case.

The Hamiltonian of the system takes the form

where $\mathcal{H}_{\text{kin}}$ describes the nearest-neighbor hopping

$\mathcal{H}_{R}$ is the Rashba spin-orbit coupling

$\mathcal{H}_{J}$ is the magnetic exchange coupling induced by the impurities

and $\mathcal{H}_{\text{SC}}$ describes the superconducting pairing

Here $c^{\dagger}_{i,s}$ is the creation operator for an electron in site $i$ and spin $s$, $\mathbf{d}_{ij}=\mathbf{r}_{i}-\mathbf{r}_{j}$, and $i\in{\rm imp}$ indicates that the summation is over the impurity sites. We assume that the impurity sites form a square lattice with lattice constant of $a$, as shown in Fig. 1(a).

The many-body interactions are included by means of an local $U$ and nearest-neighbor $V$ density-density interaction [see Fig. 1(b)] of the form

The Hamiltonian is solved at the mean-field level $\mathcal{H}_{\text{int}}\approx\mathcal{H}_{\text{int}}^{MF}$ including all the bilinear contractions of the mean-field

giving rise to hopping, pairing, chemical potential and spin-orbit coupling renormalization [50]. By studying separately the cases where only $U$ is present ($U$-interacting case) and both $U$ and $V$ are present ($UV$-interacting case) we can also get some idea how the range of the interactions affects the results. The self-consistent calculations are performed at fixed filling of the normal state $\nu\in[0,1]$. Since we assume everywhere that $\Delta\ll t$, $\nu$ also approximately describes the filling of the bands in the presence of superconductivity. We compute the Chern number numerically using a gauge-invariant description of Chern number associated with the Berry connection defined on a discretized Brillouin zone [51]. Importantly, because the sign of the Chern number is not important for our analysis, everywhere in the manuscript we denote with $C$ the absolute value of the Chern number.

For concreteness, we show in Fig. 1(c) the electronic structure of a specific example of the Yu-Shiba-Rusinov lattice model. The corresponding Berry curvature distribution $\Omega(\mathbf{k})$ is shown in Fig. 1(d), and by integrating $\Omega(\mathbf{k})$ over the Brillouin zone we obtain a Chern number $C=7$. As a consequence of $C=7$, the spectral function in a semi-infinite geometry features 7 chiral edge modes crossing the energy gap, as shown in Fig. 1(e).

In the limit when $\Delta\ll t$ and the impurity lattice constant satisfies $a\gg 1$, the non-interacting model of Eqs. (1)-(5) supports a rich topological phase diagram as a function of model parameters [12, 30, 32]. In Fig. 2(a) we show a representative phase diagram for $a=8$ as a function of filling $\nu$ and the exchange interaction $J$ for $\alpha=0.1t$ and $\Delta=0.06t$. It contains topological phases with Chern numbers ranging from $0$ to $15$ arising due to the long-range coupling between the Yu-Shiba-Rusinov states [12]. The high Chern number topological phases originate from rather peculiar situations where the longer range couplings between the Shiba states are larger than the shorter range couplings, and therefore they cover only small areas of the parameter space [see Fig. 2(a)]. This also means that the phases with large $C$ are quite fragile against the effects of perturbations and for example in the presence of strong disorder only the $C=0$ and $C=1$ phases survive [49].

The electron-electron interactions can influence the topological phases because they can modify the spatial profiles of the mean fields. We find that in the presence of interactions the superconducting mean field is modulated with the period $a$ of the magnetic impurities. A representative example of a single unit cell of the impurity lattice with the magnetic impurity located at the center of it is shown in Fig. 2(b). In the case of magnetization, the interactions do not only modify the magnitude of the mean field, but also lead to an appearance of magnetization textures where the direction of the magnetization varies within the unit cell [see Fig. 2(c)]. We find that these effects can both enhance and reduce the topological energy gap, and even cause topological phase transitions [see Figs. 2(d)-(f)]. Therefore, we conclude that the interaction effects can play an important role in the determination of the topological phase and the magnitude of the topological gap in the case of individual samples.

The model parameters are rarely known accurately in experiments. Therefore, rather than studying the system for specific model parameters, it is important to establish the statistical effects of the interactions on the existence of the topological phases with high $C$ and for the magnitudes of the corresponding energy gaps. For this purpose we have collected statistics of the gap sizes and Chern numbers for 50000 sets of model parameters in impurity lattices with lattice constant $a=:2,4,6,8$. In the sampling we used common ranges $\nu\in[0.1,0.45]$, $\alpha\in t[0.08,0.12]$, $\Delta\in t[0.048,0.072]$, $U\in t[-1.0,1.0]$ and $V\in t[-0.5,0.5]$, whereas the sampling of $J$ for each $a$ was chosen so that the topological phase diagrams contained a large number of topologically distinct phases. We used $J\in t[0.76,1.2]$ for $a=2$, $J\in t[1.4,2.1]$ for $a=4$, and $J\in t[1.8,2.7]$ for $a=6$ and $a=8$.

In Fig. 3 we show the median energy gap values for different values of $C$ in the non-interacting, $U$-interacting, $UV$-interacting cases for different impurity lattice constants. The median gap values oscillate as a function of $C$ but there is also a general tendency for the energy gaps to become smaller and smaller with increasing values of $C$. Importantly, we find that the median energy gaps have remarkably similar magnitudes and dependence on $C$ in the non-interacting, $U$-interacting and $UV$-interacting cases (see Fig. 3). This suggests that the interactions are not statistically relevant in the determination of the topological phases and magnitudes of the topological gaps within the quite large range of interaction strengths considered. We believe that the interval of the electron-electron interaction strengths is sufficiently large to draw reliable conclusions, because the interactions are screened by the superconductor and therefore larger interaction strengths are unlikely to occur in the realization of the Shiba lattice models.

The median gap values of course give only partial information about the distribution of the gap values, and it could be possible that the interactions still affect the shape of this distribution. In Fig. 4 we plot the whole statistical distribution of the gap values in the cases of small, intermediate and large Chern numbers. The distributions are practically identical in the non-interacting, $U$-interacting and $UV$-interacting cases in the case of small gap sizes. Since the distribution function decays approximately exponentially with the increasing energy gap values (see Fig. 4) the small gap values occur much more commonly than the large ones, so that the typical gap sizes are indeed the same in the interacting and non-interacting cases. However, the rare events corresponding to large topological gaps show a drastically different phenomenology. As shown in Fig. 4, for all the Chern numbers interactions substantially increase the probability of the largest topological gaps. In that regime, the non-interacting gap distribution decreases much faster than exponentially when the energy gap values approach the upper bound of the gap values of the corresponding Chern numbers that can be realized in the non-interacting Yu-Shiba-Rusinov model. In the presence of interactions, the number of samples exhibiting a large topological gap becomes much larger than in the non-interacting case, demonstrating that the interactions play an important role in the case of these rare events.

So far, our analysis has focused on a purely mean-field description of electronic interactions. Effects beyond mean-field are well known to be important for two-dimensional systems, and could give rise to several different types of corrections. First, at the mean-field level, the electronic states have effectively an infinite lifetime. However, in interacting systems the presence of interactions can give rise to finite lifetimes even for low-lying modes [52]. Furthermore, quantum fluctuations are not accounted for in the mean-field methodology, resulting in an overestimation of the symmetry breaking effects [53]. It is also worth noting that our analysis treats the magnetic impurity as a classical spin, neglecting quantum fluctuation effects. In reality, such quantum fluctuations can give rise to a crossover between Kondo and Yu-Shiba-Rusinov physics [54], which by construction are not included in our model. Finally, the effects of the interactions beyond the mean-field theory in the considered high-Chern number quasi-flat-band system could also lead to qualitatively new physics such as fractionalization, leading to more exotic non-Abelian phases.

It is interesting to note that, while a variety of methods exist for interacting two-dimensional models, the current problem would provide a great challenge for most of them. Exact diagonalization provides an accurate methodology to approach small clusters. Nonetheless, due to the large number of orbitals per unit cell, including large number of unit cells in an exact diagonalization would not be feasible. It is worth noting that such a limitation also emerges in moiré systems, where supercells contains many sites. In moiré systems, however, projection of interactions on a single or few low-energy bands allows to get a reduced effective model per supercell [55, 56]. In the Yu-Shiba-Rusinov model considered, the different bands are often highly entangled, and in such kind of situations isolating a single low energy band is not feasible. Tensor-network [57] and neural-network quantum state methodologies [58] have become one of the more powerful unbiased methods to treat two-dimensional systems. However, again due to the large size of the supercell, these calculations would become unfeasible in practice. An alternative strategy would be to use a dynamical mean-field theory [59] for the supercell considered, which would be possible for modest supercells using the current state-of-the art methodology. The dynamical mean-field theory would allow to obtain a more accurate picture of the renormalization of electronic structure due to interactions, and more importantly, to account for additional Kondo physics. Finally, GW methodology of the effective model [60] in the normal state can yield a more accurate picture of the renormalization effects on the electronic structure by electronic repulsion. While the study of the interaction effects beyond the mean-field level is certainly an interesting direction for future research, the interactions in our calculations do not lead to spontaneous breaking of symmetries. Thus, the mean-field calculations can provide a reliable description of the many-body states also in low-dimensional systems [61].

To conclude, we have studied the effects of electron-electron interactions in the two-dimensional Yu-Shiba-Rusinov lattice model on mean-field level, and our results show that the interactions can enhance or reduce the topological gap as well as cause topological phase transitions. We have shown that statistically the distributions of the Chern numbers and the topological energy gaps in the non-interacting and interacting cases are approximately the same for the most common realizations where the energy gaps are reasonably small. Interestingly, we find that the probability of rare realizations with large energy gaps, located at the tails of statistical distributions, are strongly enhanced by interactions. Such realization are challenging to find experimentally due to their small likelihood, but they would be the most suitable ones for the applications of topological superconductivity. Our results provide also a starting point for studies of more subtle correlation effects in Shiba lattice models. The effects beyond the mean-field theory are a great challenge for future theoretical research, but this research direction might be worth pursuing because these models may support unprecedented many-particle phenomenology thanks to the interplay of superconductivity, magnetism, high Chern numbers and the large spatial extent of the Yu-Shiba-Rusinov states.