# Skyrmion-induced subgap states in -wave superconductors

###### Abstract

In -wave systems, it has been theoretically shown that a ferromagnetic film hosting a skyrmion can induce a bound state embedded in the opposite-spin continuum. In this work, we consider a case of skyrmion-induced state in a -wave superconductor. We find that the skyrmion induces a bound state that generally resides *within* the spectral gap and is isolated from all other states, in contrast to the case of conventional superconductors. To this end, we derive an approximate expression for the -matrix, through which we calculate the spin-polarized local density of states which is observable in scanning tunneling microscopy measurements. We find the unique spectroscopic features of the skyrmion-induced bound state and discuss how our predictions could be employed as novel experimental probes for -wave superconducting states.

###### pacs:

74.70.Pq, 74.78.Na,74.78.FkIntroduction — Topology has played a significant role in our understanding of robust features of condensed-matter systems. Topological protection guarantees nontrivial materials properties and enables promising potential applications within electronics and technology. Materials with suitable structure provide an excellent ground to produce low-energy excitations analogous to concepts originating in particle physics. One such example is the magnetic skyrmion, a topological defect in a magnetic field which manifests as a vortex-like spin configuration bogdanov:1989:1 ; rossler:2006:1 . The magnetic configuration is then characterized by a topological invariant given by

(1) |

where is a unit vector aligned with the local magnetic field. takes on integer values and is denoted the topological charge of the skyrmion. Configurations with different topological charges are separated from each other by a finite energy barrier, making skyrmions robust excitations.

In recent years, significant experimental progress has been made within the field muhlbauer:2009:1 ; munzer:2010:1 ; yu:2010:1 ; yu:2011:1 ; heinze:2011:1 ; seki:2012:1 ; neubauer:2009:1 ; ritz:2013:1 ; romming:2015:1 . Notably, skyrmions can easily moved by applying spin currents kiselev:2011:1 ; zhang:2015:1 . Employing spin-polarized scanning tunneling microscopy, Romming et al. romming:2013:1 demonstrated a controlled method of creating and destroying individual skyrmions. As a consequence, skyrmions are of special interest from a technological perspective due to their properties being potentially suitable for use in electronics.

Concurrently with the expanded experimental possibilities, there has been a rise in interest towards skyrmion-superconductor heterostructures pershoguba:2016:1 ; yang:2016:1 ; hals:2016:1 , not least because the interplay between skyrmions and superconductivity is expected to give rise to topological systems. In -wave superconductors, skyrmions with have been theoretically shown to give rise to Yu-Shiba-Rusinov-like states with long-range wavefunctions pershoguba:2016:1 . These states are within the spectral gap of the bulk states with parallel spin-polarization, but generally still reside within the bulk with anti-parallel spin-polarization, making them resonance peaks in the density of states. Also, Skyrmions with even charge have been argued to host Majorana zero-energy states yang:2016:1 on two-dimensional (2D) -wave substrates.

In this letter, we extend the study of skyrmions on superconductors (SCs) to -wave SCs. These are of particular interest due to the rich physics stemming from anisotropic pairing that could itself support topological superconductivity. While this type of unconventional superconductivity has not been conclusively shown to exist in nature, there are nevertheless some possible candidate materials under consideration, most notably SrRuO mackenzie:2003:1 . Hence it is interesting to consider what sort of properties such a system would be expected to have. Using standard techniques to analyze the spectra of superconducting states, we find that: i) skyrmions with bind localized subgap states, which is in stark contrast to the results previously obtained for -wave systems; ii) depending on the type of -wave coupling, states bound to Bloch and Néel skyrmions show qualitatively different behaviour.

System – The object of interest in this paper is a magnet-SC heterostructure, the magnetic texture being cylindrically symmetric. The simplest case is that of a constant Zeeman field. Consider at first a system consisting of a two-dimensional -wave SC with a constant background magnetic field. The Hamiltonian density of this system is

(2) |

where is the kinetic energy, is the background magnetic field, and is the superconducting pairing function. This Bogoliubov-de Gennes Hamiltonian acts on the Nambu spinor , where the operator creates an electron of spin and momentum . The superconducting triplet pairing is encoded into the vector . In this paper, we will consider two different types of pairing vectors: out-of-plane and in-plane . These result in two different Hamiltonians

(3) |

where is for out-of-plane -vector (pairing of antiparallel spins) and for in-plane -vector (pairing of parallel spins); and are Pauli matrices in particle-hole and spin space, respectively. Throughout this work we will treat both types of -wave pairing in parallel.

When the system is coupled to a skyrmion, the magnetic-field term in the Hamiltonian acquires a spatial dependence . We will here consider the two different configurations of the magnetic field seen in Fig. 1, known as Néel and Bloch skyrmions, respectively. The fundamental difference between the two cases lies in the direction of rotation for the magnetization vector as a function of radius. The magnetic textures of each is described by the following vectors:

(4) |

where we model the position-dependent angle as

(5) |

In the above, can be viewed as the radius of the skyrmions. Calculating the topological charge for the two configurations above as per Eq. (1), we find that they are equal: for both configurations, . This indicates that the two skyrmions are topologically equivalent, and indeed it is possible to transform between the two through a unitary transformation . However, since the transformation between skyrmions does not leave the in-plane -wave Hamiltonian invariant, one might expect to see differences between the two types of skyrmions in that model. In order to ascertain the effect of the skyrmions on the -wave system, we will proceed to calculate the local density of states (LDOS), as it contains the relevant spectral properties of the system and further is amenable to experimental analysis.

Green’s function and the multipole expansion — For generic spatial dependence, an explicit analytic solution of the system including the exact magnetic texture is challenging. As a first approximation we assume that the skyrmion is small compared to the superconducting coherence length. We then perform a multipole expansion around the origin to obtain corrections to a desired order. This gives us an effective potential and allows us to employ the -matrix formalism to find an approximative solution to the full Green’s function of the system, the validity of which is tied to that of the multipole expansion. We restrict ourselves to the second order of the calculation, corresponding to the monopole/anapole term for the Néel/Bloch-type skyrmions, respectively. The expansion for both skyrmions contain a constant magnetic field term, which for the purposes of this work can be treated as part of the unperturbed background; consequently, Eq. (3) will be used as a starting point of the expansion. The remainder of the multipole terms will then be treated as a scattering potential. We further note that, in real systems, skyrmions can generally be moved around by perturbations unless they are pinned down by static terms fukuyama:1978:1 ; lee:1979:1 . This can be modelled by adding a pointlike scalar potential to the skyrmionic terms. Taking this into account, we will consider a change to the Hamiltonian of the form , where

(6) |

Correspondingly, in momentum space, the perturbative potential can be written

(7) |

The magnitudes of and can be obtained as the respective moments of the expansion. As this calculation makes no reference to the superconductivity, and the skyrmions are generally related by a simple rotation, we can in both cases simply use the magnitudes obtained from the Néel magnetic texture:

(8) |

The above relations fix the values of , as a function of the background magnetic field and the skyrmion radius . Using the truncated multipole expansion we can then find an approximation for the -matrix of the skyrmion through use of the Lippmann-Schwinger equation

(9) |

The equation can be solved analytically in the approximation where the incoming and outgoing momenta are close to the Fermi level supp . The explicit solution for the -matrix allows us to calculate the full Green’s function of the system. Inserting our result for the -matrix up to second order in the multipole expansion supp results in an expression for the Green’s function

(10) |

where we used the fact that the integrals over the two momenta can in each case be separated into two different integrals. Hence finding reduces to finding the value of the integrals

(11) | ||||

(12) |

where is the spatial Green’s function of the system without a skyrmion. The integrals are analytically tractable, and the solutions for both types of -wave pairing are presented in the supplemental material. We hence have an analytic expression for the full Green’s function in terms of these integrals. Moving on, we can use to calculate the spin-polarized local density of states (SPLDOS) through use of the formula

(13) |

where corresponds to spin up (down) electrons. The SPLDOS is useful in that it can be probed by spin-polarized STM and hence provides a direct way of comparing theory to experiment. Further, the sum of the terms for spin up and down directly yields the LDOS measured in typical STM experiments.

Results — Based on the treatment above, we have calculated the LDOS of the system. Notably, we find that, the -wave SC-skyrmion system can support subgap bound states well separated from the continuum, unlike the -wave SC, as seen in Fig. 2. An example of the LDOS of a system with parameters supporting subgap bound states can be seen in Fig. 2 (a) for the case of an out-of-plane -vector. To illustrate its location within the gap, in Figs. 2 (b) and (c) we have plotted the spin-polarized DOS at the origin for both types of -vector. In the latter two we have also included the DOS in the case of non-zero scalar potential as dashed lines. As is clear from the figures, the non magnetic potential has a quantitative effect on the energies of pre-existing bound states, although the degree and direction of this shift depends on the -vector. We find in total four bound states (including those at negative energies), consistent with results obtained for -wave Yu-Shiba-Rusinov systems kaladzhyan:2016:2 .

Interestingly, the type of skyrmion makes a significant qualitative difference in the case of an in-plane -vector. In the studied regime, the Néel-type skyrmion did not support subgap bound states at all, whereas the Bloch-type skyrmion supports subgap bound states for a wide parameter range note:skyrmiontrans . In a system with out-of-plane -vector on the other hand, the two skyrmions are equivalent and the system can support subgap states regardless of skyrmion type. When subgap states are present, they are generally localized within the skyrmion, i.e., ; the spatial decay of the LDOS is exponential, as illustrated in Fig. 3, where we have plotted the LDOS as a function of radius and energy for the relevant SC-skyrmion combinations. For completeness, we note that while the in-plane -wave SC coupled to a Néel skyrmion does not support bound states for the parameters used in Figs. 2 and 3, it can host bound states for high enough scalar potentials . However, in this case the states are clearly bound to the scalar potential rather than the skyrmion, and in fact the presence of the skyrmion increases the scalar potential needed to form a subgap bound state in this system. It is important to note that these features are rather generic: the SC-skyrmion combinations that support subgap bound states do so for a wide range of parameters.

Thus we propose that it may be possible to distinguish between -wave and different types of -wave superconductivity depending on the effect seen when the system is coupled to skyrmions. The existence of subgap bound states in general indicates that the superconductive pairing is not -wave, and dependence on the type of skyrmion can act as a separator between out-of-plane and in-plane -vectors.

Conclusions and outlook — We have found new bound states that are generated in a -wave superconductor in proximity to a ferromagnetic film hosting a skyrmion with topological charge . We predict sharp features in the non-polarized as well as spin-polarized local density of states that can be measured by tunneling spectroscopy in superconductors. To be general, we considered two types of -wave superconductors, i.e., in-plane and out-of-plane -vector. In contrast to -wave superconductors, which can only host bound states embedded in a continuum, we found that both studied types of -wave systems can support genuine subgap bound states. In the out-of-plane case, both skyrmions are equivalent due to rotational symmetry. However, for the in-plane -vector, we found remarkable qualitative differences between the two types of skyrmions: namely, a Bloch skyrmion can induce a subgap bound state in the superconductor for a wide range of parameters, while we observed no bound states for the Neel skyrmion note:skyrmiontrans . This feature could be used experimentally to investigate the character of the pairing of a given -wave superconductor.

We also found that a scalar potential can have an impact on the structure of the bound states. In general, will shift the energies of any bound states present, depending on its sign potentially bringing them to low energies or gapping the system altogether. High scalar potentials can result in bound states even for an in-plane system with a Néel-type skyrmion; however, in this case it is clear that the state is bound specifically to the scalar potential well – not the skyrmion – and in fact the presence of the skyrmion increases the minimum scalar potential for which a bound state appears.

Our prediction of subgap bound states also opens up further venues for research, in addition to being useful for distinguishing different types of superconductivity. Within the past decade, it has been experimentally shown that under some circumstances, lattices of skyrmions can form spontaneously muhlbauer:2009:1 ; munzer:2010:1 ; heinze:2011:1 . More recently, 2D lattices of hybridized subgap bound states have been found to give rise to interesting topological behavior in both -wave rontynen:2015:1 ; rontynen:2016:1 and -wave kaladzhyan:2016:1 superconductors. Similarly, we may expect the slowly decaying subgap states in a skyrmion lattice to hybridize and potentially induce interesting topological behavior.

Acknowledgements — This work was supported by the Academy of Finland and the Aalto Centre for Quantum Engineering (K.P., A.W. and T.O), ITS at ETH Zurich and by US DOE BES E3B7, ERC DM-321031 (S.S.P. and A.V.B.).

## References

- (1) A. N. Bogdanov and D. A. Yablonskii, JETP Letters 68, 101 (1989).
- (2) U. K. Rößler, A. N. Bogdanov, and C. Pfleiderer, Nature 442, 797 (2006).
- (3) S. Mühlbauer et al., Science 323, 915 (2009).
- (4) W. Münzer et al., Phys. Rev. B 81, 041203 (2010).
- (5) X. Z. Yu et al., Nature 465, 901 (2010).
- (6) X. Z. Yu et al., Nat Mater 10, 106 (2011).
- (7) S. Heinze et al., Nat Phys 7, 713 (2011).
- (8) S. Seki, X. Z. Yu, S. Ishiwata, and Y. Tokura, Science 336, 198 (2012).
- (9) A. Neubauer et al., Phys. Rev. Lett. 102, 186602 (2009).
- (10) R. Ritz et al., Nature 497, 231 (2013).
- (11) N. Romming, A. Kubetzka, C. Hanneken, K. von Bergmann, and R. Wiesendanger, Phys. Rev. Lett. 114, 177203 (2015).
- (12) N. S. Kiselev, A. N. Bogdanov, R. Schäfer, and U. K. Rößler, J. of Phys. D 44, 392001 (2011).
- (13) X. Zhang et al., Nanotechnology 26, 225701 (2015).
- (14) N. Romming et al., Science 341, 636 (2013).
- (15) S. S. Pershoguba, S. Nakosai, and A. V. Balatsky, Phys. Rev. B 94, 064513 (2016).
- (16) G. Yang, P. Stano, J. Klinovaja, and D. Loss, Phys. Rev. B 93, 224505 (2016).
- (17) K. M. D. Hals, M. Schecter, and M. S. Rudner, Phys. Rev. Lett. 117, 017001 (2016).
- (18) A. P. Mackenzie and Y. Maeno, Rev. Mod. Phys. 75, 657 (2003).
- (19) H. Fukuyama and P. A. Lee, Phys. Rev. B 17, 535 (1978).
- (20) P. A. Lee and T. M. Rice, Phys. Rev. B 19, 3970 (1979).
- (21) See supplemental material for details.
- (22) V. Kaladzhyan, C. Bena, and P. Simon, Phys. Rev. B 93, 214514 (2016).
- (23) Note that the transformation between the two skyrmion types changes the structure of the in-plane SC from the scalar product term used here to a vector product-type term, i.e., . Hence in the latter type of system, the Néel skyrmion contains bound states while the Bloch skyrmion does not.
- (24) J. Röntynen and T. Ojanen, Phys. Rev. Lett. 114, 236803 (2015).
- (25) J. Röntynen and T. Ojanen, Phys. Rev. B 93, 094521 (2016).
- (26) V. Kaladzhyan, J. Röntynen, P. Simon, and T. Ojanen, Phys. Rev. B 94, 060505 (2016).

## Appendix A -matrix

This appendix is dedicated to the derivation of the -matrix for the multipole expansion of the Skyrmion. The derivation done here closely follows that of what was done in the appendix of Ref. pershoguba:2016:1 . To simplify, we leave out the scalar impurity throughout the derivation and only reintroduce it at the very end.

Our starting point is the Lippmann-Schwinger equation for the -matrix which reads

(A.1) |

In the multipole expansion, the potential is , where the matrix is either equal to the Kronecker delta or the Levi-Civita symbol , depending on whether the Skyrmion is of Néel or Bloch type.

For out-of-plane and in-plane -vector respectively, and as elaborated in Appendix B, this is

(A.2) |

where are the projection operators along the axis – , – and . Henceforth we will denote for notational simplicity. We can insert into the Lippmann-Schwinger equation in order to calculate the full Green’s function of the skyrmion-superconductor composite system:

In order to proceed, we also make a few simplifying observations: first, the potential only consists of a momentum-independent term and a linear term. This observation together with the form of Eq. (A.1) suggests that a good ansatz for the -matrix is one with terms that are at most quadratic in momentum. The -matrix can then be written as

(A.3) |

where are matrices to be determined. Note that the follows from general symmetry arguments for the -matrix. The second observation we make is that the scattering processes primarily occur close to the Fermi surface, allowing us to write ( is a unit vector and is the Fermi momentum) for all momenta in the ansatz expression for the -matrix and the potential term in Eq. (A.3). Inserting all this into Eq. (A.1) gives us

(A.4) |

from which we by matching components get a system of equations for the matrix components:

(A.5) | ||||

(A.6) | ||||

(A.7) |

where we have introduced the the two integrals

(A.8) |

We postpone the evaluation of these integrals to the next section. The matrix components can now be obtained by first solving for in Eq. (A.6). This is achieved by multiplying Eq. (A.6) with and summing over . We can then easily express in terms of and solve the original equation for in terms of . Inserting this into Eq. (A.5) is then trivial albeit cumbersome. To solve for it is then only a matter of performing a similar multiplication and summation trick as we did for . We are finally left with

(A.9) | ||||

(A.10) | ||||

(A.11) |

where we have introduced the matrices

(A.12) |

This concludes the derivation of the -matrix in the multipole expansion. We can now add the scalar impurity term by noting that it would only appear together with , so all we need to do is replace in the expression for in (A.12) and we are done. As a final note, we emphasize that the derivation is valid for both in- and off-plane , since they only differ in whose explicit form was not used in the above derivation.

## Appendix B Integrals

In this Appendix we will evaluate some integrals encountered in the main text and in the previous appendix. Specifically, we consider the integrals in Eq. (11) in the main text,

(B.1) | ||||

(B.2) |

as well as two integrals from the previous appendix:

(B.3) | ||||

(B.4) |

We begin with the integrals from the main text. Due to convergence, it is necessary to consider the cases and separately. For the bare Green’s function, inverting and inserting suitable resolutions of identity yields the following integrals at :

(B.5) |

These integrals do not strictly converge as written, ultimately due to the fact that BCS theory is not valid at high energies. They can be calculated, however, by introducing a suitable cutoff. Here, this can simply be effected by assuming the relevant scale is near the Fermi level and hence approximating . This reduces both Green’s functions to simple residue calculations, immediately yielding the expressions in Eq. (A.2).

The Green’s function for can be calculated from

(B.6) |

The terms in the numerator proportional to the momentum can be written as derivatives with respect to the appropriate coordinate. The angular integral then directly gives a Bessel function of the first kind. Factoring the denominator in terms of its zeroes, we find that calculation of the two Green’s functions reduces to solving three integrals:

(B.7) | ||||

(B.8) |

where we have defined the integrals

(B.9) | ||||

(B.10) | ||||

(B.11) |

where is the modified Bessel function of the second kind, and . Here are the zeroes of the denominator in the respective Green’s function (indices , have been suppressed). Specifically, we have for out-of-plane superconductivity

(B.13) |

and for in-plane SC

(B.14) |

The result (B.9) is obtained through residue integration upon representing the Bessel function as an integral; the others can consequently be deduced through recurrence relations. We have in this work assumed that the argument of the square root in is positive in all cases; to first order in , this is equivalent to requiring that the energies lie within the bulk gap.

The integrals in can be solved in a similar manner. We first consider the case . This can be easily solved by factorizing the denominator as above, yielding

(B.15) | ||||

(B.16) |

Consider then the case . The angular integral is most conveniently handled by pulling out a derivative with respect to , upon which again we obtain a number of integrals over Bessel functions:

(B.17) | ||||

(B.18) |

defined in terms of

(B.19) | ||||

(B.20) | ||||

(B.21) |

where is the Struve function of the first kind, and the other parameters are the same as previously.

Finally we turn to the integrals in the previous appendix. In both cases, the angular integral is trivial and can be integrated out. Looking at the remaining integrals over the momentum, it is then immediately evident that

(B.22) |