N-independent Localized Krylov Bogoliubov-de Gennes Method: Ultra-fast Numerical Approach to Large-scale Inhomogeneous Superconductors

We propose the ultra-fast numerical approach to large-scale inhomogeneous superconductors, which we call the Localized Krylov Bogoliubov-de Gennes method (LK-BdG). In the LK-BdG method, the computational complexity of the local Green's function, which is used to calculate the local density of states and the mean-fields, does $not$ depend on the system size $N$. The calculation cost of self-consistent calculations is ${\cal O}(N)$, which enables us to open a new avenue for treating extremely large systems with millions of lattice sites. To show the power of the LK-BdG method, we demonstrate a self-consistent calculation on the 143806-site Penrose quasicrystal lattice with a vortex and a calculation on 1016064-site two-dimensional nearest-neighbor square-lattice tight-binding model with many vortices. We also demonstrate that it takes less than 30 seconds with one CPU core to calculate the local density of states with whole energy range in 100-millions-site tight-binding model.

We propose the ultra-fast numerical approach to large-scale inhomogeneous superconductors, which we call the Localized Krylov Bogoliubov-de Gennes method (LK-BdG). In the LK-BdG method, the computational complexity of the local Green's function, which is used to calculate the local density of states and the mean-fields, does not depend on the system size N . The calculation cost of self-consistent calculations is O(N ), which enables us to open a new avenue for treating extremely large systems with millions of lattice sites. To show the power of the LK-BdG method, we demonstrate a self-consistent calculation on the 143806-site Penrose quasicrystal lattice with a vortex and a calculation on 1016064-site two-dimensional nearest-neighbor square-lattice tightbinding model with many vortices. We also demonstrate that it takes less than 30 seconds with one CPU core to calculate the local density of states with whole energy range in 100-millions-site tight-binding model.
Introduction. The mean-field approach through the Bogoliubov-de Gennes (BdG) equations is one of the most convenient and efficient ways to describe inhomogeneous superconductivity. In the past two decades, studies about quasiparticle excitations in superconducting systems with junctions or vortices become more important since these systems can be stages for a topological quantum computing with the use of the Majorana quasiparticles [1][2][3]. Recently, the Majorana zero modes have been observed in many systems such as the iron-based superconductor FeTe x Se 1−x [4,5]. Since in the systems with junctions or vortices one needs to use a real-space formulation, the numerical simulation becomes computationally involved. Although there are alternative approaches to inhomogeneous superconductivity like quasiclassical Eilenberger theory or Ginzburg-Landau methods, these methods can not treat discretized quantum modes like Majorana zero modes. In addition to topological materials, there are many interesting systems to be solved such as high-T c superconductors for which the superconducting coherence length is of the order of the Fermi wavelength, or nanoscale superconductivity for which superconducting coherence length is comparable to the system size. Therefore, the need for a fully quantummechanical approach has become imperative.
It is very hard to diagonalize the Hamiltonian in large inhomogeneous systems, since the computational complexity to diagonalize the Hamiltonian matrix is O(N 3 ). Here, N is the matrix size. In the last decade, various kinds of numerical approaches to solve the BdG equations for inhomogeneous systems have been developed [6][7][8][9]. The computational complexities of these approaches are O(N 2 ) for self-consistent calculations and O(N ) for calculating the local quantities like the local density of states (LDOS), respectively. However, if one wants to treat inhomogeneous systems with internal degrees of freedom (e.g. the iron-based superconductors are multiband systems), the computational complexity becomes huge even with the use of supercomputing systems with thousands of CPU cores. Therefore, it is still hard to treat large realistic inhomogeneous systems.
In this letter, by focusing on the fact that the oneparticle local Green's function is constructed locally in real space, we propose the ultra-fast numerical approach to large-scale inhomogeneous superconductors, which we call the Localized Krylov Bogoliubov-de Gennes method (LK-BdG). We show that vectors in the Krylov subspace to calculate the Green's function are localized. The computational complexities in the LK-BdG are O(N ) for selfconsistent calculations and O(1) for calculating the local quantities, respectively. To show the power of the LK-BdG method, we demonstrate a self-consistent calculation on the s-wave superconducting 143806-site Penrose quasicrystal lattice with a vortex and a calculation on 1016064-site two-dimensional s-wave nearest-neighbor square-lattice tight-binding model with many vortices. Finally, the summary is given.
Model. The Bogoliubov-de Gennes equations describe the behavior of electrons and holes in superconductors, which are coupled to mean-fields. A general BdG Hamiltonian is given as H = Ψ †Ĥ Ψ/2. The column vector Ψ is composed of N fermionic annihilation c i and creation operators c † i (i = 1, 2, · · · , N ), . The subscription i in c i or c † i indicates a quantum index depending on spatial site, spin, orbital, etc. For simplicity, we regard i as a spatial site index. The Hamiltonian matrixĤ is a 2N × 2N Hermitian matrix given aŝ arXiv:2001.02362v1 [cond-mat.supr-con] 8 Jan 2020 Here,Ĥ N is a Hamiltonian matrix in normal states and ∆ is a superconducting order parameter. Without diagonalizing the BdG Hamiltonian directly, we can calculate physical observables and mean-fields with the use of the one-particle Green's functionĜ(z) = (zÎ −Ĥ) −1 . The important quantity is the difference of the retarded and advanced Green's function matrices determined asd(ω) =Ĝ R (ω) −Ĝ A (ω). For example, the LDOS at site i and the mean-field c i c j are expressed [6,7]. Here, e(i) and h(i) are 2N -component unit-vector defined as Localized Krylov subspace. We focus on the fact that the vectors e(i) and h(i) are localized in real space. We introduce the order-m Krylov subspace generated by the Hamiltonian matrixĤ and an unit vector b (= e(i) or h(i)) given as We call K m (Ĥ, b) the localized Krylov subspace, since m vectors are localized as follows. The element of the second vector is expressed as If the Hamiltonian matrixĤ is sparse, the number of finite value elements of [Ĥe(i)] k is a few. For example, in the case of the two-dimensional square-lattice tightbinding model with nearest neighbor hopping, there are only four elements in the normal state Hamiltonian matrix. The element of the third vector is expressed as [Ĥ 2 e(i)] k = lĤ klĤli , where the index l for the summation is restricted due to the sparseness of the Hamiltonian. Thus, the number of the finite value elements of m-th vector is ∼ m d , where d is the dimension of the system. With the use of the above discussion, we can reduce the computational complexity of the matrix-vector product: where the number of the indices in the summation does not depend on the size of the system N , as shown in Fig. 1. The first vector h 0 , defined in Eq. (2), is localized in hole-space. Although the amplitude of the vectors spreads in electron-and hole-spaces due to the superconducting order parameter matrix∆, the amplitudes of h 1 , h 2 , h 3 are still localized (See, Fig. 1). We should note that this property is useless for finding eigenvalues and eigenvectors with the use of the Krylov-subspace based methods such as Lanczos or Arnoldi methods, since the eigenvectors are usually not localized so that we should take large m in the Krylov subspace where the elements of the m-th vector are all finite. Chebyshev polynomial method. In the Chebyshev polynomial method, we generate vectors by a recurrence formula: with q 0 = b, q 1 =Kb andK = (Ĥ − bÎ)/a [6,7]. n generated vectors are in the Krylov subspace K n+1 (Ĥ, b). The mean-field c i c j is expressed as , w n = (1 + δ n0 )π/2 and f (x) = 1/(1 + exp(−x/T )). We found that the mean-field c i c j can be calculated with good enough accuracy in the localized Krylov subspace whose order m is much smaller than the dimension of the matrix if the index j is not far from the index i in real space. Therefore, the Chebyshev polynomial method with the localized Krylov subspace to calculate elements of the matrixd(ω) does not depend on the matrix dimension N . We note that Furukawa and Motome have shown that the free energy in the Monte Carlo simulations can be calculated with the use of the similar localized Krylov subspace in fermion systems coupled with classical degree of freedom [12]. According to their paper, we can introduce the truncated localized Krylov subspace, which accelerates BdG simulations more effectively as discussed later.
Lanczos method for a Green's function. Another Krylov-subspace-based method to calculate a Green's function is the Lanczos method [13]. The diagonal elements of Green's function [Ĝ(z)] ii (i ≤ N ) can be calculated by the continued fraction expansion expressed as )), where the coefficients a n and b 2 n are calculated by a recurrence formula: j n+1 =Ĥj n − a n j n − b 2 n j n−1 , with a n = j T nĤ j n /j T n j n and b 2 n = j T n j n /j T n−1 j n−1 , supplemented by b 2 0 = 0, j −1 = 0 and j 0 = e(i). In the Lanczos method for a Green's function, the n-th order Krylov subspace is given as K n (Ĥ, e(i)), which is localized. We note that the Lanczos method with the localized Krylov subspace has been used in the field of the order-N first-principles calculations [14,15]. In this letter, we adopt the Lanczos method to calculate the LDOS, since this quantity is calculated by the diagonal element of the Green's function.
Truncated Localized Krylov subspace. According to the paper written by Furukawa and Motome [12], we introduce the truncated matrix-vector product expressed as We can reduce the range where calculations are restricted, by introducing a threshold . The details are discussed in their paper [12]. In this letter, we adopt = 10 −6 . Demonstration I: quasicrystal with the Penrose lattice. We demonstrate a self-consistent calculation of quasicrystalline superconductors to show the power of the LK-BdG method. Recently, the superconductivity of quasicrystals was discovered in Al-Zn-Mg quasicrystalline alloys [16]. Sakai and Arita have studied possible superconductivity on a Penrose-tiling structure, which is a prototype of quasicrystalline structures [10,11]. They found that there exists a superconducting state with spatially extended Cooper pairs in the attractive Hubbard model. This Penrose-tiling structure that we call Penrose lattice is one of good demonstrations for the LK-BdG, since the inhomogeneous superconducting order parameter naturally appears. We consider the tight-binding model on the Penrose lattice proposed in the previous papers [10,11] (Also see the supplemental materials[17]). We introduce a vortex located at a center. We calculate a s-wave onsite superconducting order parameter c i c i with the interaction U = −3t and the chemical potential µ = −1t at the zero temperature. The Chebyshev cutoff parameter is n c = 200, which is enough to obtain the mean-field. The renormalized parameters a and b for the Hamiltonian matrix are a = 10t and b = 0, respectively. In Fig. 2, we show the self-consistent solution on the 143806-site Penrose lattice (Also see Fig. S1 in the supplemental materials[17]). By comparing with 143806-site and 21106-site systems, we show that a center region can be regarded as a bulk since the LDOS around a center does not depend on the lattice size as shown in Fig. 2(d). While the level spacing of vortex bound states in conventional systems is characterized by ∆ 2 /E F with the Fermi energy E F [20], the energy level in the Penrose lattice is close to zero. The size of the vortex core is small as shown in Fig. 2(a). In the conventional theory for the vortex bound states, a minimum energy level in a small vortex core is high due to a quantum confinement of quasiparticles. Note that this energy level depends on the position of a vortex center [18]. Because there is no Fermi wave length due to the absence of the translational symmetry, this suggests that interesting vortex physics exists in quasicrystalline superconductors.
Let us show the system size dependence of the computational complexity for self-consistent calculations. Figure 3 shows that the computational complexity is O(N ). For example, the elapsed time for one iteration step in 375971-site (143806-site) Penrose lattice with 240 CPU cores on the supercomputing system SGI ICE X at the Japan Atomic Energy Agency [19] is about 1265 (440) seconds.
Demonstration II: Ultra-large 2D tight-binding model. We demonstrate that the LK-BdG method can treat ultra-large 2D tight-binding model. We consider 2D N x × N y square-lattice nearest-neighbor tight-binding model with many vortices at zero temperature. We calculate the s-wave order parameter with U = −2.4t, µ = −1.5t, n c = 200, a = 10t, b = 0. We consider the Peierls phase to introduce vortices perpendicular to the system [7,21]. The symmetric gauge A(r) = (1/2)H × r with H = (0, 0, H z ) is used. Here, we consider m vortices in the system (H z = mφ 0 /(N x N y )) and the type II limit (the magnetic penetration depth λ → ∞). As the initial guess of the superconducting order parameter, we introduce vortices with the Penrose tiling pattern, whose phase singularities are located at vertices of the Penrose tiling. Although this vortex configuration is not a true ground state, we can obtain a similar vortex configuration as a metastable state if the vortex- vortex distance is long enough in the type II-limit superconductor. After 120 iteration steps with solving gap equations [22], we obtain the superconducting gap distribution in the 1008 × 1008 2D tight-binding model. As shown in Fig. 4, the LK-BdG method can calculate superconducting mean-fields and the LDOS in large systems with many vortices.
To compare with a previous method, we measure the elapsed time for calculating the LDOS with a single vortex located at a center. For simplicity, we consider a single vortex whose coherence length is 10 (the unit is the lattice spacing). We consider 4000 energy meshes from −2.5t to 5.5t. The number of the Lanczos iterations is 400 and the smearing factor for the LDOS is 5 × 10 −2 . The calculations with a single CPU core are done on a laptop PC (MacBook Pro (13-inch, 2018, Four Thunderbolt 3 ports) with 2.7GHz Intel Core i7 CPU with 4 cores). As shown in Fig. 5, the computational complexity of our method with the localized Krylov subspace does not depend on N , where N is the dimension of the Hamiltonian matrix. We confirm that the elapsed time in the 5000 × 5000 lattice system, whose matrix dimension is 5×10 7 , is only about 26 seconds on the laptop PC. On the other hand, the complexity of the previous Lanczos-based method is O(N ) as shown in Fig. 5. Summary. We proposed the LK-BdG, the ultra-fast numerical approach to large-scale inhomogeneous superconductors, by focusing on the fact that the vectors in the Krylov subspace for the Green's function is localized. In a self-consistent calculation, the computational complexity is O(N ). We also showed that the computational complexity to calculate the local density of states does not depend on the system size N . The LK-BdG method enables us to open a new avenue for treating extremely large systems with millions of lattice sites.

PENROSE TILING
We show the Penrose lattice that we used in Fig. S1. We regard each vertex of the rhombuses as a site, and put an electron hopping t between two sites connected by the edge of the rhombuses. The Hamiltonian is expressed as where n iσ = c † iσ c iσ . For simplicity, we neglect the Hartree mean-fields. The site-dependent superconducting order parameter is The self-consistent solution of the superconducting mean-fields is shown in Fig. S2.