Journal of the Physical Society of Japan FULL PAPERS Efficient Numerical Self-Consistent Mean-Field Approach for Fermionic Many-Body Systems by Polynomial Expansion on 2 Spectral Density 1 0 2 Yuki Nagai, Yukihiro Ota∗, and Masahiko Machida n a CCSE, Japan Atomic Energy Agency, 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8587, J 7 Japan 1 CREST(JST), 4-1-8 Honcho, Kawaguchi, Saitama, 332-0012, Japan ] n o We propose an efficient numerical algorithm to solve Bogoliubov de Gennes equations c - r self-consistently for inhomogeneous superconducting systems with a reformulated poly- p u nomial expansion scheme. This proposed method is applied to typical issues such as a s . t vortex under randomly distributed impurities and a normal conducting junction sand- a m wiched between superconductors. With various technical remarks, we show that its - d efficiency becomes remarkable in large-scale parallel performance. n o KEYWORDS: mean-fieldtheory,superconductivity,Bogoliubov-deGennesequations,parallelcomputation c [ 3 v 9 3 1. Introduction 9 4 The mean-field approach is one of the most convenient and efficient ways to clarify . 5 0 relationshipbetweenmicroscopicdescriptionsandmacroscopicphenomenaincondensed 1 1 matter physics. So far, this approach has been applied to various fermionic interacting : v many-body systems. Recently, inhomogeneous superconducting systems such as nano- i X scale superconducting wires and artificial/intrinsic Josephson junctions have been its r a intriguing application examples.1–4) The Bogoliubov-de Gennes (BdG) equation5) is a theoretical starting-point in these systems. Its self-consistent solution gives information such as non-trivial quasi-particle excitation-spectra and inhomogeneous superconduct- inggap.However, suchamean-fieldBdGapproachhasbeenregardedtobenotpractical for a long time, since its calculations require huge computation resources beyond the contemporary standard level of present computers. Very recently, a few groups6–8) have proposed a highly efficient numerical method to ∗ Present address: Advanced Research Institute, RIKEN, 2-1 Hirosawa,Wako-shi, Saitama 351-0198, Japan 1/16 J. Phys. Soc. Jpn. FULL PAPERS solve the BdG equations by using the kernel polynomial expansion.9) In these papers, the key idea is to expand Green’s function with a set of the Chebyshev polynomials. Thepolynomialexpansiondrasticallyreducesthecomputationcostandhasanexcellent parallel efficiency in contrast to the conventional way, i.e., direct full diagonalization of the BdG-Hamiltonian. In the history of condensed matter physics, such a Chebyshev- based expansion of Green’s function goes back to seminal studies by Tanaka et al.10,11) TheseauthorscalculatedGreen’sfunctionsforagenericHamiltonianwithuseofvarious kinds of orthonormal polynomials. Sota et al.12) developed this idea as an oscillation- free Fourier expansion scheme. Afterwards, several papers have been published in terms of the polynomials expansion of Green’s function.9,13) Thus, the polynomial-based ex- pansion has been an attractive numerical method in large-scale fermionic many-body systems. Specifically, the application to inhomogeneous superconducting systems is an important target since spatial profile of the superconducting gap has key information to predict their physical properties. The aim of this paper is to develop a fast and tractable method to self-consistently calculate the BdG equations for various inhomogeneous superconducting systems. We reformulate the polynomial expansion scheme in a more comprehensible manner and performfirst self-consistent calculations ontypical inhomogeneous superconducting sys- tems. Atfirst, weclaimthatthepolynomialexpansion isappliedtonotGreen’sfunction itself but the spectral density of Green’s function. By expanding the Dirac’s delta func- tion on the spectral density, we obtain a full mean-field calculation scheme. The present tool, spectral-density polynomial expansion allows a straightforward numerical calcula- tion with the mean-fields, i.e., superconducting gap and/or general Hartree-Fock terms. Next,weclaimthatthemean-fielditselfconverges fasterthanthelocaldensity ofstates. We focus ons-wave superconductivity to demonstrate its efficiency. We remark that our approach can be effective for more general cases including d-wave superconductivity in magnetic fields14) and multi-orbital superconductivity with spin-orbit coupling. This paper is organized as follows. The polynomial-expansion reformula is given in Sec. II. The present style is theoretically more complete and easier to do programing on parallel cluster machines. In Sec. III., we demonstrate three typical examples; an s-wave disordered superconductor in the magnetic field, i.e., a vortex formation under randomly distributed impurities, a nano-scale superconductor-normal-superconductor (SNS) junction, i.e., superconducting proximity effects, and a numerical challenge, i.e., a reproduction of full temperature dependence of the gap amplitude from zero to T . c 2/16 J. Phys. Soc. Jpn. FULL PAPERS Someconvenient technical remarksforpracticalsimulations aregiveninSec. IV.Section V is devoted to the summary. 2. Formulation We give a full formula based on the orthonormal-polynomial expansion to solve the BdG equations. The essence is that the spectral density of Green’s functions can be expanded by orthonormal polynomials. Afterwards, we show the explicit expressions with use of the Chebyshev polynomials. The present formalism can be applied to any fermionic quadratic Hamiltonian with mean-fields. 2.1 Hamiltonian We start with the Hamiltonian associated with the BdG equations. Covaciet al.6) proposed a way to calculate the eigenvalues and the eigenvectors of this Hamiltonian without full diagonalization. We describe their formula in a more general way. Let us consider a Hamiltonian for a fermion system given as H = Ψ† ˆΨ/2. The H column vector Ψ is composed of N fermionic annihilation c and creation operators i c† (i = 1, 2,..., N), Ψ = ( c , c† )T, where c = (c , c ,..., c )T and c† = i { i} { i} { i} 1 2 N { i} (c†, c†,..., c† )T.TherowvectorΨ† isalsodefinedasΨ† = ( c† T, c T).ThesymbolT 1 2 N { i} { i} † meanstransposition.Thefermioniccanonicalanti-commutationrelationleads[c ,c ] = i j + † δ . The subscription i in c or c indicates a quantum index depending on spatial site, ij i i spin, orbital, etc. The “Hamiltonian” matrix ˆ is a 2N 2N Hermite matrix given as H × Aˆ Bˆ ˆ =  , (1) H Bˆ† AˆT  −  where Aˆ and Bˆ are N N complex matrices. These matrices have the relation given as × Aˆ† = Aˆ, BˆT = Bˆ, (2) − because of the hermitian property of H and the fermionic canonical anti-commutation relations. When we consider a superconductor, ˆ corresponds to the mean-field H Bardeen-Cooper-Schrieffer (BCS) Hamiltonian and Bˆ contains the superconducting gap. 3/16 J. Phys. Soc. Jpn. FULL PAPERS 2.2 BdG equations The BdG equations are regarded as th eigenvalue equation with respect to ˆ ex- H pressed as ˆf = ǫ f (γ = 1,2,...,2N), (3a) (γ) γ (γ) H u (γ) f =  . (3b) (γ) v (γ)   The column vectors u and v are N-component complex vectors. To solve the BdG (γ) (γ) equations is equivalent to diagonalization of ˆ with a unitary matrix Uˆ (Ref. 15), H Uˆ† ˆUˆ = Dˆ, Dˆ = diag(ǫ ,ǫ ,...,ǫ ). (4) 1 2 2N H Theeigenvaluesǫ ’sarenotindependent ofeachother,i.e.,ǫ = ǫ (i = 1, 2,...,N). c i i+N − The matrix elements of Uˆ lead as U = u U = v . (5) iγ (γ),i i+Nγ (γ),i 2.3 Spectral density ˆ Here, we concentrate on a key quantity, i.e., spectral density (or discontinuity) d(ω), which is a 2N 2N matrix, to solve the BdG equations using orthonormal polynomials. × ˆ All essential physical observables are described by bilinear forms with respect to d(ω). Now, let us define the Green’s function as Gˆ(z) = (z ˆ)−1, which is a 2N 2N −H × ˆ ˆ complex matrix. With the use of the unitary matrix U, each component of G(z) is expressed as 2N 1 G (z) = U U∗ (1 α,β 2N). (6) αβ αγ βγz ǫ ≤ ≤ X γ γ=1 − If we set z = iω with the Matsubara frequency ω = (2n + 1)/β, the above formula n n corresponds to the temperature Green’s function.16) The retarded and the advanced Green’s functions are, respectively, defined as GˆR(ω) = lim Gˆ(ω +iη), (7a) η→0+ GˆA(ω) = lim Gˆ(ω iη). (7b) η→0+ − The spectral density16) is given as a difference between the retarded and the advanced Green’s functions, dˆ(ω) GˆR(ω) GˆA(ω), whose matrix elements are expressed as ≡ − 2N dˆ(ω) = 2πi U U∗ δ(ω ǫ ). (8) h iαβ − X αγ βγ − γ γ=1 4/16 J. Phys. Soc. Jpn. FULL PAPERS ˆ In order to obtain physical observables (e.g., density of states) from d(ω), we introduce the following useful 2N-component unit-vectors e(i) and h(i) (1 i N), which are, ≤ ≤ respectively, defined as [e(i)] = δ , [h(i)] = δ . (9) γ i,γ γ i+N,γ Specifically, employing e(i) and h(i), the elements of the column vectors u and v are γ γ rewritten as u = [e(i)TUˆ] , (10) (γ),i γ v∗ = [Uˆ†h(i)] . (11) (γ),i γ Furthermore, the local density of states with respect to the site i is given as 1 N(ω,i) = e(i)Tdˆ(ω)e(i), (12) −2πi with use of dˆ(ω) and e(i), since N(ω,i) is defined as 2N N(ω,i) = u 2δ(ω ǫ ), (13) (γ),i γ | | − X γ=1 N N = u 2δ(ω ǫ )+ v 2δ(ω +ǫ ). (14) (j),i j (j),i j | | − | | X X j=1 j=1 A typical self-consistent BdG calculation for a superconductor requires two types of mean-fields c†c and c c . These mean-fields can be expressed as, h i ji h i ji 1 ∞ c†c = dωf(ω)e(j)Tdˆ(ω)e(i), (15a) h i ji −2πi Z −∞ 1 ∞ c c = dωf(ω)e(j)Tdˆ(ω)h(i), (15b) h i ji −2πi Z −∞ with f(x) = 1/(eβx +1). Here, β is the inverse temperature β = 1/T. 2.4 Orthonormal polynomial expansion In this paper, we focus on an orthonormal polynomial φ (x) with interval [-1,1] n (n = 0,1,...). In principle, various orthonormal polynomials are applicable to solve the BdG equations (as shown in Ref.11). A set of polynomials φ (x), which is assumed to n be a real function with respect to x, fulfills the relations ∞ W(x) δ(x x′) = φ (x)φ (x′), (16a) n n − w X n n=0 1 w δ = φ (x)φ (x)W(x)dx. (16b) n n,m Z n m −1 5/16 J. Phys. Soc. Jpn. FULL PAPERS A recurrence formula is generally given as17) φ (x) = (a +b x)φ (x) c φ (x). (16c) n+1 n n n n n−1 − In order to confine the eigenvalue range inside the interval, we rescale the energy scale ˆ of by the following manner H ˆ bIˆ ǫ b ˆ = H− , ξ = γ − , (17) γ K a a where a = (E E )/2 and b = (E +E )/2 with E ǫ E . It seems max min max min min γ max − ≤ ≤ that these relations require a tough computing task to obtain E and E . However, max min their rough estimations are practically enough. We employ a convenient criterion to determine E and E from physical information. See Sec. 4.2 for detail. max min Now, let us derive a formula using the polynomial expansion. At first, we define a matrix form by polynomial functions as 2N φ (ˆ) = U U∗ φ (ξ ), (18) h n K iαβ X αγ βγ n γ γ=1 where φ (ξ ) is well-defined in the interval ξ [ 1,1]. This implies that ω-integrals n γ γ ∈ − in Eqs. (15a) and (15b) are also bound in the finite energy range. According to a similar manner to Eq. (17), these integral intervals become [ 1,1], i.e. ω = ax+b with − ˆ x [ 1,1]. Substituting the right hand side of Eq. (16a) for the definition of d(ω), we ∈ − have ∞ 2πi W(ω) pTdˆ(ω)q = φ (ω)pTq , (19) n n − a w X n n=0 for arbitrary 2N-component real vectors p and q. A sequence of the vector q ( n ≡ φ (ˆ)q) is recursively generated by n K q = (a +b ˆ)q c q (n 2), (20a) n+1 n n n n n−1 K − ≥ q = φ (ˆ)q, q = φ (ˆ)q. (20b) 1 1 0 0 K K The coefficients of Eq. (20a) are the same as the ones in Eq. (16c). Accordingly, the mean-fields (15a) and (15b) are expressed as ∞ c†c = e(j)Te (i)Tn, (21a) h i ji n w X n n=0 ∞ c c = e(j)Th (i)Tn, (21b) i j n h i w X n n=0 6/16 J. Phys. Soc. Jpn. FULL PAPERS where 1 = dxf(ax+b)W(x)φ (x), (22) Tn Z n −1 e (i) = φ (ˆ)e(i), h (i) = φ (ˆ)h(i). (23) n n n n K K Here, it should be noted that does not depend on the index i. Therefore, the calcula- n T tion of can be done before any self-consistent calculations. The essential mathemat- n T ical relations are Eqs. (16a)-(16c). Many useful formulae about orthogonal polynomials applicable to physical problems are found, for example, in Refs.18,19) Hereafter, as shown by by Covaci et al.,6) we use the Chebyshev polynomials,17) i.e., φ (x) = cos[narccos(x)], (24a) n 1 π W(x) = , w = (1+δ ), x [ 1,1]. (24b) n n0 √1 x2 2 ∈ − − The coefficients in the recursive formula (16c) are a = 0, b = 2, and c = 1. The n n n vector form of the formula associated with Eq. (20a) is given as q = 2ˆq q (n 2), (24c) n+1 n n−1 K − ≥ with q = q and q = ˆq. At the zero temperature (β ), we find that 0 1 K → ∞ = π arccos( b/a), (25a) 0 T − − sin[narccos( b/a)] = − . (25b) n6=0 T − n Equation(24c)allowstheevaluationoftherighthandsideofEq.(19)withoutanydirect diagonalization of ˆ. Specifically, we find that the calculations of the mean-fields (21a) H and (21b) do not require any heavy computation in contrast to matrix diagonalization. Finally, we remark that the approach shown in this section is applicable to solve the other equations in condensed matter physics such as the Kohn-Sham equation in the density functional theory with real-space formalism. 3. Results: Numerical Demonstrations We demonstrate three examples of self-consistent calculations in inhomogeneous su- perconducting systems by using the above spectral-density polynomial expansion. In this section, we focus on a sigle-band s-wave superconductor. Therefore, the quantum index i(j) introduced in the previous section indicates a single spatial site. Then, the spin index ( ) is separately written. The label (i, j) simply means a spatial site in ↑ ↓ real 2D space. Together with the BdG equations, the gap equation for s-wave supercon- 7/16 J. Phys. Soc. Jpn. FULL PAPERS ductivity is given as ∆ = V c c . In numerical calculations, its right hand side is ij ij i,↓ j,↑ h i expressed by using Eq. (21b) with truncation in the infinite summation. The maximum number of the summation is written as n . For simplicity, we consider a simple square c lattice tight-binding model only with nearest neighbor hopping. The hopping magni- tude is t. In all simulations, we adopt n = 1000. We find that this number is enough c from convergent tendency except for the calculations of the temperature dependence of ∆ . ij 3.1 Vortex lattice Solution The first example is a vortex solution and related issues. The electromagnetic re- sponse of a superconductor under the magnetic field has been examined on the basis of the BdG equations.20–22) In type II superconductors, the low-lying vortex-core excita- tions can be studied by taking account of the presence of vortex lattice. Here, we investigate a two-dimensional N N site system with periodic boundary x y × condition. This system has a vortex square lattice. Moreover, we introduce randomly distributed impurities at the zero temperature. The parameters are set as follows: V = ij 2.2tδ , the chemical potential µ = 1.5t, the system size N N = 64 64, and ij x y − − × × the impurity potential V = t. We successfully calculate the gap amplitude as shown imp in Fig. 1. We note that a calculation with 20 times iterations takes about 40 minutes Cross section at y = 32 0.2 ||∆∆|| 0.2 0.22 0.15 0.15 0.2 0.18 0.16 0.1 0.14 00 0..01.182 0.05 ∆|| 0.1 00..0046 0 0.02 50 60 0.05 40 10 20x 30 40 50 60 10 20 30 y 0 0 10 20 30 40 50 60 70 (a) (b) x Fig. 1. (Color online) (a) Spatial modulation of the gap amplitude in a two-dimensional vortex lattice system with ten impurities (b) Cross section of the spatial modulation at y =32. by a desktop computer with 8 CPU cores (Intel Xeon X5550 2.67 GHz 2). Since the × Hamiltonian is a sparse matrix, the calculation needs little computational memories. The separate calculations of ∆ is performed on each CPU core. The communication in ij this case, which is a one-to-all communication, is needed only when updating ∆ . We ij 8/16 J. Phys. Soc. Jpn. FULL PAPERS 000...111 SSSCCC NNNooorrrmmmaaalll SSSCCC 6)6)6) |(Mean-field potential)| =1=1=1 000...000888 0.12 YYY at at at 0000....00013692 0000...000369 d potential)| (d potential)| (d potential)| ( 000000......000000464646 0 elelel 28 n-fin-fin-fi 21 eaeaea 000...000222 12 2X4 36 48 7 14 Y |(M|(M|(M 000 111888 222111 222444 222777 333000 333333 333666 333999 444222 444555 444888 (a) (b) XXX Fig. 2. (Color online) (a) Spatial modulation of mean-field potential in a two-dimensional SNS junction. (b) Superconducting proximity effect in the vicinity of the normal region at Y =16. confirm that these advantage points are still true in more general cases such as d-wave and multi-orbital superconductors.23) 3.2 SNS Junction System The next example is the proximity effects in an SNS Josephson junction, in which the gap symmetry of the superconducting electrodes is assumed to be s-wave. Joseph- son junctions and superconducting weak links have attracted great attention in vari- ous research fields such as superconducting device engineering24) including Josephson qubits,25) superconducting transport studies in superconducting wires,26) fundamental studies on unconventional nano-superconductors,7,27–29) and so on. Numerical simula- tions for the BdG equations has been often carried out to study Josephson effects.30–32) Now, let us show a spatial modulation of the mean-field, i.e., the superconducting gap in the SNS junction with spatial size N N = 64 32. The normal conducting x y × × part is inside the region (x,y);27 x 37, 1 y N . We set V = 0 inside y ij { ≤ ≤ ≤ ≤ } this region. Otherwise, V = 2.2tδ . The chemical potential µ = 1.5t and the ij ij − − temperature T = 0. The hopping between the superconducting and normal regions t′ = 0.8t. We note that the averaged gap amplitude outside the normal region (i.e., in a uniform superconductor) is 0.198 in the present parameter set. It indicates that the coherence length is around 5 sites (Ref.3). The mean-field superconducting-gap c c i,↑ i,↓ h i distribution is shown in Figs.2(a) and (b). In this case, the self-consistent iteration number is 20, which is confirmed to be enough through its convergence check. The mean-field takes a non-zero value even in the normal region. We find that the spatial modulation of the mean-field is well characterized by the length scale compatible with 9/16 J. Phys. Soc. Jpn. FULL PAPERS 0.25 30 iteration steps 300 iteration steps 0.2 BCS gap function 0.15 ∆/t 0.1 0.05 0 0 0.05 0.1 0.15 0.2 T/t Fig. 3. (Coloronline)Temperaturedependence ofthe gapamplitude inthe two-dimensionals-wave superconductor. the coherence length ( 5 sites). ∼ 3.3 Temperature Dependence of Gap Amplitude We demonstrate temperature dependence of the superconducting gap amplitude. We consider a two-dimensional N N square-plate s-wave superconductor. We set × V = 2.2tδ , the chemical potential µ = 1.5t, and the system size N N = 28 28. ij ij − − × × The 30 and 300 times iterations are adopted to evaluate the temperature dependent su- perconducting gap, respectively. As shown in Fig. 3, we successfully obtain the temper- ature dependence of the gap amplitude,33) although the convergence tendency depends onthetemperature range.These obtainedtemperature dependences arefittedby anan- alytical formula of the BCS gap function ∆(T) = ∆(0)tanh(1.74 ∆(0)/(1.764T) 1). − p The slow convergence appears around T , since the superconducting gap almost van- c ishes near the critical point. However, this is not a fault of the expansion scheme. It is known that the diagonalization scheme also shows similar tendency. A calculation with 300 times iterations takes about 3 minutes when using a parallel cluster with 112 cores. 4. Technical remarks Inthissection,werefertotechnicalremarkstodescribetheadvantagesofthepresent polynomial expansion scheme. These points are useful when actually performing large- scale numerical calculations. 10/16