Multi-atom quasiparticle scattering interference for superconductor energy-gap symmetry determination

Complete theoretical understanding of the most complex superconductors requires a detailed knowledge of the symmetry of the superconducting energy-gap Δkα\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathrm{{\Delta}}}_{\mathbf{k}}^\alpha$$\end{document}, for all momenta k on the Fermi surface of every band α. While there are a variety of techniques for determining ∣Δkα∣\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|{\mathrm{{\Delta}}}_{\mathbf{k}}^\alpha |$$\end{document}, no general method existed to measure the signed values of Δkα\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathrm{{\Delta}}}_{\mathbf{k}}^\alpha$$\end{document}. Recently, however, a technique based on phase-resolved visualization of superconducting quasiparticle interference (QPI) patterns, centered on a single non-magnetic impurity atom, was introduced. In principle, energy-resolved and phase-resolved Fourier analysis of these images identifies wavevectors connecting all k-space regions where Δkα\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathrm{{\Delta}}}_{\mathbf{k}}^\alpha$$\end{document} has the same or opposite sign. But use of a single isolated impurity atom, from whose precise location the spatial phase of the scattering interference pattern must be measured, is technically difficult. Here we introduce a generalization of this approach for use with multiple impurity atoms, and demonstrate its validity by comparing the Δkα\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathrm{{\Delta}}}_{\mathbf{k}}^\alpha$$\end{document} it generates to the Δkα\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathrm{{\Delta}}}_{\mathbf{k}}^\alpha$$\end{document} determined from single-atom scattering in FeSe where s± energy-gap symmetry is established. Finally, to exemplify utility, we use the multi-atom technique on LiFeAs and find scattering interference between the hole-like and electron-like pockets as predicted for Δkα\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathrm{{\Delta}}}_{\mathbf{k}}^\alpha$$\end{document} of opposite sign.


INTRODUCTION
The macroscopic quantum condensate of electron pairs in a superconductor is represented by its order-parameter Δ α k / hc αy k c αy Àk i, where c αy k is the creation operator for an electron with momentum k on band α. But electron pair formation can occur through a wide variety of different mechanisms and in states with many possible symmetries 1 . Thus, it is the symmetry properties of Δ α k that are critical for identification of the Cooper pairing mechanism 1 and, moreover, for understanding the macroscopic phenomenology 1 . While macroscopic techniques can reveal energy-gap symmetry for single-band systems 2,3 , no general technique existed to determine the relative signs of Δ α k and Δ β k 0 between k α and k β for all Fermi surface (FS) momenta in an arbitrary superconductor.
In 2015, a conceptually simple and powerful technique for determining Δ α k symmetry was introduced 4 , by Hirschfeld, Eremin, Altenfeld, and Mazin (HAEM). It is based on interference of weakly scattered quasiparticles at a single, non-magnetic, impurity atom. Given a superconductor Hamiltonian where H k is the normal-state Hamiltonian and Δ k the superconducting energy gap, a non-magnetic impurity atom is modeled as a weak point-like potential scatterer, with Hamiltonian H imp ¼ V 0 c y r c r centered at the origin of coordinates r = 0. Effects of scattering are then represented by a T-matrix derived from the local Green's function When the impurity potential is con- with the T-matrix given by where V imp is the impurity matrix.
From G k;k 0 E ð Þ, the perturbations to the local density-of-states δN(r, E) are predicted surrounding the impurity atom, and its Fourier transform can be determined directly from Δ k as which is a purely real quantity because, in the theoretical calculation, the single impurity is exactly at the origin of coordinates. The authors of ref. 4 realized that the particle-hole symmetry of Eq. (2) for scattering interference wavevector q ¼ k β f À k α i , depends on the relative sign of the energy-gaps Δ α ki and Δ β k f at these two momenta. Consequently, the experimentally accessible energy-antisymmetrized function ρ − (q, E) of phase-resolved Bogoliubov scattering interference amplitudes can be used to determine the relative sign of the superconducting energy-gaps connected by q ¼ k β f À k α i . In the simplest case with two isotropic gaps Δ α and Δ β on distinct bands, it was demonstrated that where E þ ¼ E þ i0 þ , so that the functional form of ρ À q; E ð Þ is very different when the product Δ α Δ β is positive or negative. An elementary implication of Eq. (4) is that, when the order parameter has opposite signs on the two bands so that Δ α Δ β < 0, ρ − (q, E) does not change sign and exhibits pronounced maxima or minima near E ≈ Δ α,β whereas if the order parameter has the same sign so that Δ α Δ β > 0, ρ − (q, E) exhibits weak maxima or minima near E ≈ Δ α,β with a sign of change in between. More generally, especially with multiple bands and anisotropic gaps, HAEM requires that ρ − (q, E) be predicted in detail for a specific H k and Δ k in Eq. (1) and then compared with quasiparticle interference imaging 5 in which the scanning tunneling microscope (STM) differential electron tunneling conductance, g r; E ð Þ / δN r; E ð Þ is visualized. This single-atom phase-resolved HAEM method has been established experimentally 6,7 . For example, in the case of the multiband s ± superconductor FeSe, the complete energy and wavevector dependence of ρ − (q, E) was used to determine that the k-space structure including relative sign of Δ α k and Δ β k , for all k α and all k β on two different bands. But this result required that the impurity atom be highly isolated from other impurities and centered precisely at the origin of coordinates, with respect to which the ReδN(q, E) of Eq. (3) is then properly defined. This was critical because, an error of on the order of~1% of a crystal unit cell in the coordinate of the origin (at the impurity atom) produces significant errors in ReδN(q, E) and ImδN(q, E) (Supplementary Note 1 and Fig. S1). Moreover, single impurity atom-based measurements limit the k-space resolution because the field of view (FOV) is typically restricted in size, making them unsuitable for superconductors with large impurity-atom densities. This provides the motivation for a variety of approaches to Δ α k determination beyond single-atom HAEM. One is to study Bogoliubov bound-states at individual impurity atoms [8][9][10] , although this has proven problematic because the elementary HAEM concept (Eq. (3)) is only valid in the weak scattering range, i.e. well below the scattering strength sufficient to generate Bogoliubov bound states 11 . Another approach is to use sparse blind deconvolution 12 to analyze images of scattering interference at multiple atoms, yielding the phase-resolved real space structure of δN(r, E) although not the ρ − (q, E) of Eq. (3). Overall, therefore, widespread application of the HAEM technique (Eq. (3)) as a general tool for Δ α k determination remains a challenge. Here, we introduce a practical technique for determining ρ − (q, E) of Eq. (3) from multiple impurity atoms in a large FOV. To understand this approach, consider the key issue of phase analysis as depicted in Fig. 1, a schematic simulation of Friedel oscillations Intensity (I 0 ) Re ( ) Im ( ) Intensity (I 0 ) Intensity (I 0 ) Intensity(I 0 ) Intensity (I 0 ) e. f. g. . e Imaginary part of Fourier transform ReδN MA (q) calculated using multi-atom technique of Eq. (7). f ReδN(q) from δN(r) in a for ϑ ¼ 0 and ϑ ¼ π, integrated azimuthally from b. Its strong random fluctuations versus |q| are due to summing the Friedel oscillations in δN(r) of a with random phases due to the random locations R i . g ReδN MA (q) from δN(r) in a integrated azimuthally from d. ReδN MA (q) is now orders of magnitude more intense than in f, and the phase of the Friedel oscillations in δN(r) of a is now very well defined because the effects of random locations R i are removed by using Eq. (7). Note that, now, changing the oscillation phase ϑ ¼ 0 to ϑ ¼ π surrounding all R i in δN(r) produces the correct evolution of ReδN MA (q).
atoms at random locations R i . The Fourier transform components of this δN(r) are shown in the top two panels of Fig. 1b. Obviously, the ReδN(q) required for the HAEM technique in Eq. (3), is weak, does not have a clear sign, and is indistinguishable from ImδN(q). Such effects occur because the spatial phases of all the individual Friedel oscillations at R i are being added at random. The consequence is most obvious in the azimuthally integrated ReδN(q) shown in Fig. 1f where the phase information of single-atom Friedel oscillation is completely scrambled and the HAEM technique of Eq. (3) thereby rendered inoperable. This problem could be mitigated if the Fourier transform of the scattering interference pattern surrounding each R i were evaluated as if it were at the zero of coordinates. In this regard, consider the Fourier transform of a scattering interference surrounding a single impurity atom at This "shift theorem" shows how the correctly phase-resolved Fourier transform of a δN i (r) oscillation centered on an atom located at R i = (x i , y i ), can be determined using where δN(q) is the Fourier transform using the same arbitrary origin as determines the R i . Thus we may define a multi-atom phase-preserving algorithm for QPI The consequences of Eq. (7) are illustrated in Fig. 1d, e (Supplementary Note 3 and Supplementary Fig. S4). The real part ReδN MA (q) now becomes well-defined and the overall magnitude is also strongly enhanced compared to ReδN(q). Moreover, the azimuthally integrated ReδN MA (q) plotted in Fig. 1g shows that the sign of ReδN MA (q) changes for ϑ ¼ 0 and ϑ ¼ π as expected.
Here it is essential that the impurity atom coordinates R i be determined accurately so that the phase is well-defined. We therefore employ a picometer-scale transformation 13  showing the scattering between hole-pocket α and electron-pocket ε with scattering vector p 1 , which is the subject of study. The delta pocket is predicted in LDA calculations, hence it is shown dim in the image, but is now reported not to exist in reality 6,38 . Due to the orbital content, the scattering between quasi-parallel Fermi surfaces would be strongly suppressed in this orbitally selective material 6 . c ρ À MA ðq; E ¼ 1:05 meVÞ calculated using Eq. (9) for a FOV containing 17 Fe-vacancies. The circle denotes the region where the α ! ε scattering occurs and we integrate the ρ À MA ðq; EÞ over the q in this region. Black crosses denote the Bragg peaks. d The integrated ρ À MA ðEÞ (dots, black) from our MAHAEM analysis of FeSe compare to the theoretical predictions from an accurate band-and gap-structure model of FeSe for s ++ (solid, pink) and s ± (solid, black) superconducting energy gap symmetry, and to measured ρ ÀExp single ðEÞ (crosses, black) from single impurity analysis as reported in ref. 6 . Clearly, the single atom ρ ÀExp single ðEÞ and the MAHAEM ρ À MA ðEÞ are in good agreement. Note that the 2D plots may show both red and blue colors due to non-ideal nature of real experimental data. However, subsequent to the integration over relevant q-space region, the ρ À MA ðEÞ is well defined as demonstrated here.
recorded g(r,E) to register all the scattering interference oscillations precisely to the crystal lattice (Supplementary Note 2). Equation 7 then allows to correctly define the quantities in Eq. (3) for arbitrarily large numbers of scattering atoms. By using the analog of Eq. (6) for g r; E ð Þ / δNðr; EÞ, ρ − (q, E) for each impurity atom is determined from while from Eq. (7) the sum over these ρ À i q; E ð Þ yields This procedure adds all the individual ρ À i q; E ð Þ signals from every impurity atom at R i in-phase, while effectively averaging out the random phase variations due to both locating the origin and the contributions of all other scatterers (Supplementary Fig. S5). We designate this procedure multi-atom HAEM (MAHAEM).

RESULTS AND DISCUSSIONS
Multi-atom quasiparticle interference for Δ α k determination Determination of the magnitude of superconducting energy gaps jΔ α k j has long been achieved [16][17][18][19][20][21][22][23] using quasiparticle scattering interference (QPI). MAHAEM is the most recent advance of the QPI technique, and to test it we consider FeSe where the single impurity atom HAEM technique for determining Δ α k was established experimentally 6 . We measure the differential tunneling conductance gðr; EÞ dI=dVðr; EÞ in a 30 nm FOV at T = 280 mK, followed by determination of R i = (x i , y i ) for 17 scattering sites (Supplementary Note 3), some of which are shown in the FOV in Fig. 2a (Supplementary Fig. S2 shows all the sites). These sites are well-known Fe-atom vacancies identified by their crystal locations, and are non-magnetic 6 ; their empirical identicality is confirmed by high-resolution electronic structure imaging. We then use Eq. (9) to calculate ρ À MA ðq; EÞ. Figure 2b shows the FeSe FS with the hole-pocket α around Γ-point and electron pockets ε (δ) around X(Y) points. Scattering between α and ε at wavevector p 1 was studied. A representative layer ρ À MA ðq; E ¼ 1:05 meVÞ is shown in Fig. 2c, where the scattering feature at vector p 1 is marked with a circle. We then sum over the encircled q-region to get ρ À MA ðEÞ for this scattering feature which is shown as black dots in Fig. 2d. Results from our MAHAEM measurements agree very well with the experimental results using a single impurity atom ρ ÀExp Single (black crosses) and the theoretically predicted curve for ρ ÀTh s ± (solid, black) in FeSe. This demonstrates the validity and utility of the multi-atom HAEM technique.
Next we consider LiFeAs, a complex iron-based superconductor that is a focus of contemporary physics interest [24][25][26] , particularly the relative sign of Δ α k between all five bands. Figure 3b shows the FS of LiFeAs calculated using a tight-binding fit 27 18 . The electron-hole scattering near q eh appears as a "horn"-shaped feature which is enclosed by a circle. d Theoretical prediction for single atom ρ − (E) integrated over the circular region shown in Fig. 3c for both s ± (black) and s ++ (pink) symmetry.
around Γ-point and two electron pockets e 1 and e 2 around Xpoint. The hole pockets around Γ−point on the FS revealed by spectroscopic imaging STM (SI-STM) 18 and confirmed by angle resolved photoemission spectroscopy (ARPES) 29,30 , are much smaller as compared to most other Fe-based superconductors. Local density approximation (LDA) and dynamical mean field theory (DMFT) calculations have attributed the small size of hole pocket to stronger electron-electron correlation in this material. The superconducting energy-gaps Δ α k are substantially anisotropic 18 . Theoretically, in the case of Δ α k with s ± symmetry, if both electron-like and hole-like pockets are present 31,32 , the pairing arises from spin-fluctuations which are enhanced by nesting between the electron-like and hole-like pockets. But the presence of three hole pockets, combined with relatively weak spin fluctuations 33 , allow for several possible competing ground states in the presence of repulsive interactions. In ref. 34 , it was pointed out that, under these conditions, several s-wave channels are nearly degenerate. These channels include the s ± state where the signs on all hole pockets are the same 35,36 and opposite to the signs on the electron bands, so-called "orbital antiphase state" that occurs when the interaction is diagonal in orbital space 24 , and a distinct sign structure obtained when vertex corrections were included 36 . Reference 37 considered the question of whether these various proposed phases could be distinguished using HAEM based on Eq. (3) and concluded that it would be challenging.
Here we examine the relative signs of Δ α k in LiFeAs by using MAHAEM. Figure 3a shows the typical cleaved surface of LiFeAs. The scattering sites used in our analysis are Fe-atom vacancies which are non-magnetic ( Supplementary Fig. S3). The theoretical simulations for LiFeAs were performed from the experimentally fitted tight binding model 27 and anisotropic gap magnitude structure 18,30 . At wavevectors corresponding to electron-hole scattering in q-space, a "horn-shaped" feature in g(q, E) appears within which we focus on an exemplary scattering vector q eh indicated by a dashed arrow in Fig. 3b. Figure 3d then shows the theoretical, single-atom ρ −Th (q, E) integrated for the q in the brown oval in Fig. 3c for s ± and s ++ gaps, where sign of the gap was imposed by hand. The sign of ρ ÀTh s ± does not change for the energy values within the superconducting gap and its amplitude peaks at the energy E % Δ e1 Δ h1 , both characteristics of a sign changing gap 37 ; contrariwise ρ ÀTh sþþ changes sign indicative of same sign energy gaps throughout. For comparison, differential conductance g(r, E) imaging of LiFeAs is performed at T = 1.2 K. The typical g(E) spectrum consists of two gaps corresponding to Δ 1 = 5.3 meV and Δ 2 = 2.6 meV. The measured g(q, E) are shown in Fig. 4a and the feature at q eh expected from the theoretical model in Fig. 3c is indicated by a circle. We evaluate ρ À MA ðq; EÞ from Eq. (9) for N = 100 atomic scale Fe-atom vacancy sites (Supplementary Note 4). The resulting Fig. 4 LiFeAs superconducting energy-gap symmetry from MAHAEM. a The measured |g(q)| pattern recorded in the FOV with multiple atomic scattering sites. The hole-like to electron-like scattering as predicted in Fig. 3c is detected clearly and indicated by a brown circle at the same location as 3c. The features inside the circle do not appear identical because the intensity in real experimental data falls off at higher |q|, while the theoretical simulation is replicated about the q-space Brillouin zone boundary with the same intensity, creating a non-physical equal intensity reflected feature. A Gaussian mask of σ = 0.68 Å −1 = 0.8π/a is used to suppress the |g(q)| data in the |q| ≈ 0 core region, to allow clearer presentation of the high-|q|-under study. Much of this |q| ≈ 0 signal intensity is believed to emanate from long-range disorder and is, moreover, unrelated to the scientific objectives of this paper. b The measured ρ À MA ðq; E ¼ 3:33 meVÞ using Eq. (9); it is typical of all ρ À MA ðq; EÞ between 1 and 6 meV. The circle indicates the hole-like to electron-like scattering in a. We integrate the ρ À MA ðq; EÞ over the range of q within this region. The dashed lines are guide to eye to a feature which is consistent with intra h 3 scattering. Again, a Gaussian mask of σ ¼ 0:68 Å À1 ¼ 0:8π=a is used to suppress the intense core emanating from long range disorder, to allow clearer presentation of the ρ À MA q; E ð Þ information at high-|q|. This suppressed area corresponds to~18% of the total area of the first qspace Brillouin zone, and the unprocessed ρ À MA ðq; EÞ data for this figure are provided in its entirety at Supplementary Fig. S6a. c The resulting ρ À MA ðEÞ calculated by summing over the oval enclosed region in b (black dots), and the theory curves (solid) for s ± (black) and s ++ (pink) overlaid. Clearly, this demonstrates using MAHAEM that the superconducting energy gap symmetry of LiFeAs is s ± (black).
Of note in Fig. 4b is the variety of structures at q j j ( jq eh j, which are challenging to understand. The thin outer blue ring (indicated by dashed light blue curve as guide-to-eye) is located at a radius in q-space that corresponds well to the expected intraband scattering within pocket h 3 . Furthermore, much of the q-space within this ring is blue and of rather high intensity for 1 meV < |E| < 6 meV ( Supplementary Fig. S6a shows dashed contours for various possible inter-hole-band scatterings overlaid on the unprocessed ρ − (q, E)). The blue color, indicating signpreserving scattering, is consistent with the conventional s ± picture within a HAEM scenario, but the high intensity is not. As discussed in Supplementary Note 5 there are several possible explanations of these low |q| phenomena, including strong scattering, quasiparticle bound states, and antiphase holepocket gaps.
Nevertheless, when the high |q| scattering between hole-like and electron-like pockets (Fig. 3b, c) is integrated within the qspace region shown by a brown circle on the ρ À MA ðq; EÞ of Fig. 4a, it yields ρ À MA ðEÞ as plotted in Fig. 4c. The theoretically predicted ρ À ðEÞ curves are overlaid for comparison. It is clear that the experimental ρ À MA ðEÞ is consistent with the ρ ÀTh ± ðEÞ theory because it does not change sign and exhibits a peak at E % 3:7 meV % ffiffiffiffiffiffiffiffiffiffi Δ 1 Δ 2 p . In this way, the MAHAEM technique efficiently demonstrates that Δ α k changes sign between electronlike and hole-like bands of LiFeAs.

Conclusions
We report development and demonstration of an improved approach for signed Δ α k determination (Eq. (9)), but now for use with multiple impurity atoms or scattering centers. This MAHAEM technique for measuring ρ − (q, E) is based on a combination of the Fourier shift theorem and high-precision registry of scatterer locations. It extends the original HAEM approach 4 to more disordered superconductors (Figs. 2a, 3a), enables its application to far larger fields of view thereby enhancing q-space resolution (Fig. 4b), and greatly increases signal-to-noise ratios (Figs. 1d, 4b) by suppressing phase randomization in multi-atom scattering interference. Overall, MAHAEM now represents a powerful and general technique for Δ α k determination in complex superconductors.

Sample growth and preparation
FeSe samples with T c ≈ 8.7 K were prepared using KCl 3 /AlCl 3 chemicalvapor transport and LiFeAs samples with T c ≈ 15 K were grown using LiAs flux method. The highly reactive LiFeAs samples are prepared in a dry nitrogen atmosphere in a glove box.

SI-STM measurements and analysis
All samples are cleaved in situ in our ultra-high cryogenic vacuum STM at low temperature. The g(r, E) data were acquired with a 3 He-refrigeratorequipped STM. The picometer level atomic registration was performed before applying the HAEM technique as described in full detail in the Supplementary Note 2. Full details of the multi-atom HAEM analysis are presented in detail in Supplementary Note 3. Theoretical predictions for ρ − (E) curves were performed using the T-matrix formalism with energy gap on each band and normal state tight binding parameters fitted to experiments.

DATA AVAILABILITY
The datasets generated and/or analyzed during this study are available to qualified requestors from the corresponding author.