Self-consistent-field ensembles of disordered Hamiltonians: Efficient solver and application to superconducting films

Our general interest is in self-consistent-field (scf) theories of disordered fermions. They generate physically relevant subensembles (“scf ensembles”) within a given Altland-Zirnbauer class. We are motivated to investigate such ensembles (i) by the possibility to discover new fixed points due to (long-range) interactions; (ii) by analytical scf theories that rely on partial self-consistency approximations awaiting a numerical validation; and (iii) by the overall importance of scf theories for the understanding of complex interaction-mediated phenomena in terms of effective single-particle pictures. In this paper we present an efficient, parallelized implementation solving scf problems with spatially local fields by applying a kernel-polynomial approach. Our first application is the Boguliubov-deGennes theory of the attractive-$U$ Hubbard model in the presence of on-site disorder; the sc fields are the particle density $n\left(\mathbf{r}\right)$ and the gap function $\mathrm{\Delta }\left(\mathbf{r}\right)$. For this case, we reach system sizes unprecedented in earlier work. They allow us to study phenomena emerging at scales substantially larger than the lattice constant, such as the interplay of multifractality and interactions or the formation of superconducting islands. For example, we observe that the coherence length exhibits a nonmonotonic behavior with increasing disorder strength already at moderate $U$. With respect to methodology our results are important because we establish that partial self-consistency (“energy-only”) schemes as typically employed in analytical approaches tend to miss qualitative physics such as island formation.

Figure 1

Benchmarking intranode parallelization and code performance. Left: Speedup with the number of cores per process for different system sizes. The performance dips (green and blue: near 16; red: near 7, 14, 21, 25) with the rising number of cores we assign to a hardware issue related to caching [parameters: $L=96$ (blue, center trace), $L=192$ (green, top trace), and $L=288$ (red, bottom trace).] Right: Performance check comparing the matrix-free implementation (orange, bottom trace) with standard mkl_sparse_d_mm of the MKL Sparse BLAS library (blue, top trace). One iteration corresponds to one sparse matrix-vector product. The ratio of the timings of the MKL and matrix-free algorithms is shown in red at $L=192$ and $L=384$.

Figure 2

Left: Disorder averaged gap $\overline{\mathrm{\Delta }\left(\mathbf{r}\right)}$ in the $U\text{−}W$ parameter plane. Parameters: $L=64$; $N_E=500$, $\alpha =1\%$. Right: Density of states for typical samples shown at four characteristic points. The red lines indicate the energies at which the LDoS is investigated in Fig. 3. Parameters: $W=0.5$ (bottom) and 1.5 (top) and $U=1.5$ (left) and 3.0 (right), $L=192$; $N_{\mathcal{C}}=6144$, $\alpha =3\%$.

Figure 3

Distribution of the local density of states (LDoS), $\rho (E,\mathbf{r})$. Left: Spatial distribution for a typical sample at peak energy of DoS ($E=0.11$, cf. Fig. 2). Right: Corresponding distribution function of LDoS, ${\mathcal{P}}_{\mathcal{\text{ld}}}$, at energies $E=0.11,1.0,3.0$ illustrating the flow of the distribution with $E$. In Fig. 2 the corresponding DoS can be found (parameters: $W=1.5$, $U=1.5$; energy resolution 0.013; $N_{\mathcal{C}}=6144$, $\alpha =3\%$).

Figure 4

Distribution of the local gap function, ${\mathcal{P}}_{\mathcal{\text{lg}}}$, with interaction strength $U$ for a typical sample with $L=192$ at weak disorder, $W=0.5$ (left plot), and stronger disorder, $W=2.0$ (right plot). As reference energy the pairing amplitude of the clean system, ${\mathrm{\Delta }}_{\text{BCS}}\left(U\right)$ has been chosen (parameters: $U=0.8,1.5,3.0$; $N_{\mathcal{C}}=8192,3072,1024$, $\alpha =0.1\%$).

Figure 5

The gap autocorrelation function ${\mathrm{\Phi }}_{\text{lg}}\left(q\right)=\overline{{\left|\mathrm{\Delta }\left(q\right)\right|}^2}$ in logarithmic representation for $L=192$ and $U=1.5$ at two values of disorder, $W=0.5$ (left) and $W=2.5$ (right); $N_E\approx 900$–1000, $\alpha =0.1\%$, $N_{\mathcal{C}}=1024$.

Figure 6

The nontrivial part of the inverted normalized gap autocorrelation function ${\mathrm{\Phi }}_{\text{lg}}^{-1}\left(\mathbf{q}\right)={\overline{{\left|\mathrm{\Delta }\left(\mathbf{q}\right)\right|}^2}}^{-1}$ evolving with $W$ at fixed $U$. ${\mathrm{\Phi }}_{\text{lg}}$ is shown averaged over equivalent directions $(\pi /a,0)$, $(0,\pi /a)$ (full symbols) and $(\pi /a,\pi /a)$, $(\pi /a,-\pi /a)$ (open). The inset shows a blowup of the small wave-number regime where open and closed symbols collapse, so all traces are isotropic (parameters: $U=1.5$, $W=1.5$ (orange), 2.5 (blue), 3.5 (green), $L=192$; $N_{\mathcal{C}}=1024$, $N_E\approx 600$–1000, $\alpha =0.5\%$).

Figure 7

Variation of ${\mathrm{\Phi }}_{\text{lg}}(\pi /L,0)$ and ${\mathrm{\Phi }}_{\text{lg}}\left(0\right)$ (red) and the correlation length, $\xi$ (blue), with increasing disorder (left, $U=1.5$) and increasing interaction (right, $W=2.5$). ${\mathrm{\Phi }}_{\text{lg}}\left(0\right)$ coincides with ${\mathrm{\Phi }}_{\text{lg}}(\pi /L)$ within the symbol size as portrayed here. The error bars depict the ensemble average error. The uncertainty due to cutoff $\alpha$ for ${\xi }_W$ is discussed in the Appendix [parameters (left): $N_{\mathcal{C}}=1024$, $N_E\approx 600$–1000, $\alpha =0.1\%$; parameters (right): $N_{\mathcal{C}}=16384$, $N_E\approx 500$, $\alpha =3\%]$.

Figure 8

Evolution of islands with disorder increasing from top to bottom, $W=0.5,2.0,3.5$. Different self-consistency schemes are compared. Left column: Full self-consistency. Center column: Energy-only self-consistency with inhomogeneous Hartree shift. Calculation is done with the single-particle (“screened”) potential as it results from the full scf calculation, left. Right column: Energy-only scheme. The energy-only data have been calculated employing full diagonalization (parameters: $U=1.5$; $N_{\mathcal{C}}=1024$, $\alpha =0.5\%$).

Figure 9

Corresponding plot to Fig. 8 in 3D. A representative 2D slice of a sample is shown. The data have been calculated employing full diagonalization (parameters: $L=24$, $W=4.0$, $U=2.5$, $n=0.3$; $\alpha =0.5\%$).

Figure 10

Gap autocorrelation function ${\mathrm{\Phi }}_{\text{lg}}\left(\mathbf{q}\right)=\overline{{\left|\mathrm{\Delta }\left(\mathbf{q}\right)\right|}^2}$ calculated employing two different energy-only self-consistency schemes. ${\mathrm{\Phi }}_{\text{lg}}$ is shown along directions $(\pi /a,0)$ (full symbols) and $(\pi /a,\pi /a)$ (open); traces for four different disorder values are shown, $W=0.05,0.5,1.5,2.5$, from bottom to top. Left: Energy-only self-consistency with screened potential. Right: Energy-only self-consistency (parameters: $U=1.5$, $L=96$; $N_E=1000$, $\alpha =0.1\%$).

Figure 11

Development of $\xi$ with cutoff $\alpha$ for disorder strengths from bottom to top W = 3.5 (black), 0.5 (green), 1.5 (blue), 2.0 (purple), and 2.5 (red), error bars depict the uncertainty stemming from the ensemble average. The dashed lines show a linear fit accounting for the three smallest $\alpha$ values (parameters: $U=1.5L=192$; $N_E\approx 600$–1000).

Figure 12

Replot of the data Fig. (11) to illustrate the (converged) variation of $\xi$ with $W$.