Novel application of the kernel polynomial method to inhomogeneous superconductivity

Inhomogeneities (surface, interfaces, impurities, etc.) in superconductors give rise to interesting phenomena, like broken time-reversal states, bound states near surfaces, etc. Numerical solutions of the self-consistent Bogoliubov-de Gennes mean field equations become computationally intensive for systems whose translational symmetry is broken. We propose a new method of solving the mean-field equations based on the Kernel Polynomial Method by expanding the Green's functions in terms of Chebyshev polynomials and calculating the order parameters self-consistently. The benefits of this method are multiple: usage of large systems, easy implementation of symmetries, multiple bands. Although we apply this method to a specific example (formation of Andreev states in 2D superconductors), it is applicable to any mean-field calculation.