Orbitally resolved superconductivity in real space: FeSe

TL;DR

This paper introduces a novel low-temperature STM method to directly resolve orbital contributions to superconductivity in FeSe at atomic resolution, bypassing traditional scattering techniques and combining experimental data with theoretical calculations.

Contribution

It presents a new approach using STM to identify orbital-specific superconducting gaps in FeSe without relying on surface scatterers or QPI measurements.

Findings

Successfully resolved superconducting gaps within the unit cell.

Identified the orbital character of each superconducting gap.

Demonstrated the method's capability on defect-free crystals.

Abstract

Multi-orbital superconductors combine unconventional pairing with complex band structures, where different orbitals in the bands contribute to a multitude of superconducting gaps. We here demonstrate a fresh approach using low-temperature scanning tunneling microscopy (LT-STM) to resolve the contributions of different orbitals to superconductivity. This approach is based on STM's capability to resolve the local density of states (LDOS) with a combined high energy and sub unit-cell resolution. This technique directly determines the orbitals on defect free crystals without the need for scatters on the surface and sophisticated quasi-particle interference (QPI) measurements. Taking bulk FeSe as an example, we directly resolve the superconducting gaps within the units cell using a 30 mK STM. In combination with density functional theory calculations, we are able to identify the orbital character of each gap.

Equations

$$\rho \left(\mathbf{r},\omega \right)=\sum _{\alpha }{\rho }_{\alpha }\left(\mathbf{r},\omega \right).$$

$$\rho \left(\mathbf{r},\omega \right)$$

$$G_{\sigma {\sigma }^{\prime }}^R\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)=\sum _{\alpha \beta }\sum _{ij}{G}_{\sigma {\sigma }^{\prime },\alpha \beta }^R\left({\mathbf{R}}_i,{\mathbf{R}}_j,\omega \right){\varphi }_{\alpha }\left(\mathbf{r}-{\mathbf{R}}_i\right){\varphi }_{\beta }^{\ast }\left({\mathbf{r}}^{\prime }-{\mathbf{R}}_j\right)$$

$$G_{\sigma {\sigma }^{\prime },\alpha \ne \beta }^R\left({\mathbf{R}}_i,{\mathbf{R}}_i,\omega \right)=0,$$

$$H=\sum _{\mathbf{k}}\sum _{\alpha \beta }{ϵ}_{\alpha \beta }\left(\mathbf{k}\right)d_{\alpha ,\mathbf{k}}^{†}{d}_{\beta ,\mathbf{k}},$$

$$H=\sum _{\mathbf{k}}\sum _{\Gamma }{h}^{\left(\Gamma \right)}\left(\mathbf{k}\right)T^{\left(\Gamma \right)}\left(\mathbf{k}\right).$$

$$G_{\alpha\beta}\left(z,\mathbf{k}\right)=\big{[}z-\sum_{\Gamma}h^{\left(\Gamma\right)}\left(\mathbf{k}\right)\lambda_{\alpha\beta}^{\left(\Gamma\right)}\big{]}^{-1}.$$

$$\hat{G}\left(z,\mathbf{k}\right)$$

$$={\hat{G}}_0\left(z,\mathbf{k}\right)+{\hat{G}}_0\left(z,\mathbf{k}\right)\hat{V}\left(\mathbf{k}\right){\hat{G}}_0\left(z,\mathbf{k}\right)+\dots$$

Full text

In this Letter, we show that by reducing the distance of the tip to the surface, orbitally resolved information on the gap can be obtained by lateral variations of the tunneling spectra within the unit cell. This technique may be very helpful to gain additional information on many multi-band superconductors and the direct observation of the orbitals may clarify some controversies regarding the nature of the involved bands. As a model system to demonstrate this approach, we chose bulk FeSe as it has the simplest crystalline structure among iron-based superconductors and does not require doping to become superconducting Hsu et al. (2008a); Hirschfeld et al. (2011); Böhmer & Kreisel (2017); Mizuguchi & Takano (2010); Coldea & Watson (2018). The unit cell of FeSe contains a layer of two Fe atoms and two layers Se atoms, one above and one below the Fe plane and the multi-band electronic structure near the Fermi level consists of three hole pockets at the $\Gamma$-point and two electron pockets at the M-point Watson et al. (2015a); McQueen et al. (2009); Watson et al. (2017a); Massat et al. (2016); Yamakawa et al. (2016); Ishizuka et al. (2018). It is established that the $d_{xz},d_{yz},d_{xy}$ orbitals of Fe dominate the Fermi surface. However, nematic order adds complexity to this material with a structural phase-transition at $T_{s}=$90 K  from a high-temperature tetragonal phase ($a=b\neq c$, space group: $P4/nmm$) to a low-temperature orthorhombic phase ($a\neq b\neq c$, space group: $Cmma$) with no magnetic ordering. Due to the 0.3% orthorhombic distortion, the unit cell is only of $C_{2}$ rotational symmetry with two inequivalent mirror planes (see Figure 1b). As a consequence, the degeneracy between the $d_{xz}$ and $d_{yz}$ orbitals is lifted. A splitting of 50 meV was observed with angle-resolved photoemission spectroscopy (ARPES) Zhang et al. (2015); Watson et al. (2015a, 2016); Nakayama et al. (2014); Maletz et al. (2014); Shimojima et al. (2014); Suzuki et al. (2015). The size of the splitting is larger than expected for an orthorhombic lattice distortion Suzuki et al. (2015); Watson et al. (2015b, 2017b, a) and was thus attributed to electronic nematicity Shimojima et al. (2014); Suzuki et al. (2015); Kostin et al. (2018).

Ultimately, a complex nature of superconductivity in FeSe arises, that occurs below a critical temperature of $\approx$8 K Hsu et al. (2008b). The superconducting gap symmetry of single-crystalline FeSe has been discussed widely in the literature. Some report a nodal gap Song et al. (2011); Kasahara et al. (2014) and more recent papers a nodeless pairing symmetry Jiao et al. (2016); Sprau et al. (2017); Kreisel et al. (2017); Benfatto et al. (2018). Even though there are some variations with respect to the size of the various gaps, a uniform observation is the multigap structure. Using QPI patterns, Sprau et al. could identify two anisotropically gapped bands without nodes of mainly $d_{yz}$ character Sprau et al. (2017). However, no consensus about the exact orbital contributions of each gap has been reached Kreisel et al. (2015); Mukherjee et al. (2015); Benfatto et al. (2018). To address these questions, we focus on STM images and spectra spatially resolved within the unit cell.

Figure 1a shows an STM topography of the sample cleaved in situ at low temperatures. The white dots reflect the upper (or lower) Se atoms of the surface layer, depending on the tunneling conditions Kreisel et al. (2016). The regular pattern illustrates the translational symmetry of the lattice of FeSe with practically no defects. To first order, the tunneling current in STM is given by the LDOS of the sample leaking out into the vacuum integrated over the bias energy window, as has been shown by Tersoff and Hamann Tersoff & Hamann (1985). Thus, STM topography in constant current mode reflects the iso-surface of the energy integrated LDOS. For low sample bias, the LDOS of a normal conductor does not vary significantly and can be considered a constant. Thus, STM topography at low bias directly represents the iso-surface of the electron density at the Fermi level $E_{F}$ decaying into the vacuum. In the superconducting state, only the quasi-particle states outside the gap $\Delta$ contribute to the LDOS. If the bias voltage $U$ is set significantly above $\Delta$, the integrated LDOS in the superconducting state is nearly that of the normal state. In that respect, the STM image reflects the LDOS near $E_{F}$ of the normal state. In theory, the LDOS of the normal state is given by the single particle excitations of the electrons, i.e. by the single-particle spectral function weighted by the form factor of Bloch states, at the respective energy. These states take into account the translational symmetry of the crystal as they are a product of a lattice periodic wave function $u_{\lambda}^{\boldsymbol{k}}(\boldsymbol{r})$ and a plane wave $e^{i\boldsymbol{kr}}$, where $\boldsymbol{k}$ is the wave vector, $\boldsymbol{r}$ is the position in real space and $\lambda$ denotes the quantum numbers related to the atomic states as the orbital degrees of freedom or the spin. The band structure describes the energy $E$ as function of $\boldsymbol{k}$ and $\lambda$. In order to study multi-band superconductors, either ARPES is used to directly measure $E(\boldsymbol{k})$. Alternatively, laterally resolved STM spectra of the sample including scatterers on the surface causing quasi-particle interference patterns can be imaged. By Fourier transforming the observed QPI patters into momentum space, information on the bands can be obtained. Since QPI shows the scattering intensity of the impurities as function of momentum transfer $\boldsymbol{q}$ between the initial and final wave vectors Sprau et al. (2017), one typically compares QPI patterns calculated from theoretical band structures with the experimental one to reveal the nature of the underlying bands. Thus, QPI mainly focusses on the band dispersion.

We here take the opposite approach to gain complementary information. In samples without defects, the Bloch states cause a LDOS that does not vary from unit cell to unit cell (no QPI patterns). Instead, the atomic part of the wave function causes LDOS variations within the unit cell directly reflecting the quantum numbers $\lambda$. Thus, the atomically resolved image of Fig. 1a basically shows the iso-surface of $\sum_{\lambda}|u_{\lambda}^{\boldsymbol{k}}(\boldsymbol{r})|^{2}$. Figure 1c shows a topographic map of the size of 2$\times$2 unit cells. The unit cell and its symmetry are indicated by solid back lines (translational symmetry) and dashed lines (mirror planes). The white dots represent the positions of the upper Se atoms. The image was created by averaging over 16 individual unit cells using the translational symmetry and by using the two mirror planes. This procedure significantly reduces the statistical noise in the STM data. Besides the Se atoms, a low intensity and finer structure becomes visible. Figure 1d displays an iso-LDOS map within the surface unit cell calculated from first-principles. It was obtained by determining the surface of constant LDOS in the vacuum in front of the surface corresponding to the setpoint conditions of the STM experiment (see Methods). The experimental and calculated maps agree reasonably well (see below for more details).

Figure 2a illustrates the energy dependence of the LDOS near $\Delta$ recorded at 30 mK Balashov et al. (2018). When taking relatively mild tunneling conditions for tip stabilization ($U$= 5 mV and $I$=320 pA, green curve), the $dI/dV$ curves, which are proportional to the LDOS, clearly shows a nodeless gap between $\pm\approx$150 $\mu$V, and two clear coherence peaks at $\approx$ 2.2 and 1.3 meV in agreement with previous results on the anisotropic gap of FeSe Sprau et al. (2017); Chen et al. (2017). Thus, we can validate the full gap with our energy resolution of $\approx$24 $\mu$eV. The two shoulders inside the lowest energy coherence peaks stem from the anisotropic gap. Under these mild tunneling conditions, we probe the evanescent Bloch states relatively far from the topmost atoms. It is well established, that far out the LDOS is dominated by s-electrons Hofer et al. (2003); Takahashi et al. (2016) as s-electrons decay the slowest into the vacuum. Note that due to the lifting of the continuous rotation symmetry by the crystal, the angular momentum is not a strict quantum number and all orbitals partially mix into s-states. When, however, increasing the tunneling conductance by approaching the tip further ($U$= 5 mV and $I$=1 nA, blue curve), the contributions of orbitals with larger orbital momentum to the LDOS will increase Takahashi et al. (2016). As a consequence, the gap spectrum gains more structure. The shoulders develop into coherence peaks and the wider coherence peaks split into several peaks. As will be shown later, this is due to states of different orbital character and different gaps. Most interestingly, when recording $dI/dV$ spectra at low distance at different lateral positions (see Fig. 2b) the intensities of the features but not their energy vary. Note that an increase of the intensity at positive bias goes hand in hand with an increase at the same negative bias. This hints for orbitally selective superconducting gaps.

In order to disentangle the orbital contributions, we carried out the following experiments and calculations. First, we recorded $dI/dV$ spectra within 12 unit cells with a resolution of 11x11 points in the individual unit cell. Similar to the topography, we use translational, rotational and mirror symmetries to average the data. Figure 3a shows the averaged and symmetrized extracted LDOS. Note, however, that the individual spectra vary within the unit cell, as will be discussed below, but the peak positions observed in the low current spectra can all be found in the averaged high-current spectrum. Next, we decompose the LDOS near $E_{F}$ from the DFT calculations into contributions of the different orbitals on the afore-mentioned iso-plane of the total LDOS of Fig. 1d. The patterns within the unit cell are displayed in Figure 3c. The calculations include all s, p and d- orbitals. It can be seen, that the composition of the LDOS varies dramatically within the unit cell for the different orbitals. Note that due to the different quantum numbers of the orbitals, no interference terms for the electrons tunneling into the different orbitals (final states) is expected. Thus, the LDOS can be written as a linear combination of the partial LDOS of the different orbitals (see Supplementary). This gives the opportunity to decompose the experimentally observed LDOS patters into their constituents regarding the orbital quantum numbers. Figure 3d displays pairs of the observed pattern and the simulated pattern as well as the weights of the composition of the different orbitals (bar graphs) at energies as indicated. While the s-state still dominates the current, clearly variations of the weights of the other orbitals are observed. These are responsible for the variation of the patterns. The experimental and simulated patterns agree well within the capabilities of STM and our first-principles calculations. The p-states are in general of higher intensity than the d-states. Moreover, the d-states that extend more into the vacuum ($d_{z^{2}}$ more than $d_{xz}$ and $d_{yz}$) have a higher intensity. The in-plane $d_{xy}$ and $d_{x^{2}-y^{2}}$ states have undetectable weight. This agrees with the expectation of the spatial distribution of the orbitals. Figure 3b plots the contributions of all detected orbitals as function of energy. The s-spectrum agrees nicely with the low-current spectrum of Fig. 2a. Most interestingly, the individual coherence peaks in the data of Fig. 3b show clear selectivity to specific orbitals. For example, the lowest coherence peak at 0.28 meV is mostly of $p_{x}$ character, and the highest peak at 1.08 meV is mostly of $p_{z}$ character and the shoulder at 1.12 meV has a large $p_{y}$ character. These states stem from the p-orbitals of Se. Similarly, the contribution of the Fe d-states to the different coherence peaks varies largely. The coherence peak at 0.82 meV is composed equally of $d_{xz}$ and $d_{yz}$. In the nematic phase, this is characteristic for the hole pocket at the $\Gamma$-point Watson et al. (2015a, 2017b); Kreisel et al. (2015); Mukherjee et al. (2015); Benfatto et al. (2018); Guterding et al. (2017); Yamakawa et al. (2016); Ishizuka et al. (2018). Thus we can attribute this gap to the hole pocket. The coherence peak at 1.12 meV shows a sizeable contribution of $d_{yz}$ but not of $d_{xz}$. Thus, the symmetry is broken which is characteristic for the elliptical electron pockets at the M-point and we can associate the gap to the electron pocket.

We hope that our analysis regarding the p-orbitals of Se and their contributions to the bands will stimulate further calculations on the hybridization of the p- and d-states. This would allow to identify more of the coherence peaks. Finally, when the contribution of specific states peaks at a certain bias voltage, QPI patterns recorded at that voltage may also provide momentum information. As such, we see this approach as a method to extend the use of STM in superconductivity research, especially in combination with QPI.

Vanishing of orbital interference terms in the local density of states

In the following, we provide evidence for the assumption that no interference terms for tunneling electrons tunneling into orbitals exist. It implies that the local density of states is given by the sum of partial density of states of different orbitals $\alpha$,

[TABLE]

Our analysis assumes tetragonal crystal systems and is therefore applicable to bulk or surface electron systems hosted in FeSe as discussed in the main text. The result is furthermore based on the assumption that there is no (weak) spin orbit coupling. Finite spin orbit coupling may potentially alter the result such that interference terms are present. This aspect will be discussed in more detail at the end of this section.

The following analysis relies on group theoretical arguments and is therefore independent on the exact microscopic model. We start by expressing the local density of states probed in STM measurements in terms of the retarded single-particle correlator

[TABLE]

with spin index $\sigma$. Here, averaging over spin degrees of freedom assumes implicitly spin-rotational symmetry of the tip/probe system. Expressed in localized Wannier orbitals $\phi_{\alpha}\left(\mathbf{r}\right)$ forming a complete set of single-particle states, the single-particle propagator reads

[TABLE]

with orbital indices $\alpha,\beta$ and the location of Bravais lattice site $\mathbf{R}_{i}$ (the extension to a $n$-dimensional basis is straightforward). In the following, we assume that the Wannier orbitals decay exponentially on length scales of the lattice spacing. It is therefore sufficient to consider the overlap of wave functions at the same lattice site $i=j$ only. In what follows, we will show that orbital off-diagonal terms, $\alpha\neq\beta$, of the single-particle correlator vanish,

[TABLE]

which leads to the vanishing of orbital interference terms in Eq. (1).

We consider the single-particle Hamilton operator represented in the Bloch basis by

[TABLE]

with electronic creation-/annihilation operator $d_{\alpha,\mathbf{k}}^{\left(\dagger\right)}$ and orbital dependent dispersion relations $\epsilon_{\alpha\beta}\left(\mathbf{k}\right)$. Electronic single-particle states are labeled by the crystal momentum $\mathbf{k}$ and orbital indices $\alpha,\beta$. We rewrite the Hamilton operator using operator bilinear forms $T^{\left(\Gamma\right)}\left(\mathbf{k}\right)=\sum_{\alpha\beta}d_{\alpha,\mathbf{k}}^{\dagger}\lambda_{\alpha\beta}^{\left(\Gamma\right)}d_{\beta,\mathbf{k}}$ with $(\lambda_{\alpha\beta}^{\left(\Gamma\right)})$ being matrices in orbital space, which transform under the point symmetry operations of the underlying lattice according to the one-dimensional irreducible representation $\Gamma\in\left\{A_{1},A_{2},B_{1},B_{2}\right\}$,

[TABLE]

Dispersion relations are contained in $h^{\left(\Gamma\right)}\left(\mathbf{k}\right)$ which have well defined transformation behaviors depending on $\Gamma$. The single-particle correlator is consequently given by

[TABLE]

We now distinguish between contributions of trivial ($\Gamma=A_{1}$) and non-trivial ($\Gamma\neq A_{1}$) transformation behavior and rewrite the the previous equation by introducing the orbital matrix $\hat{G}=(G_{\alpha\beta})$ as

[TABLE]

with $\hat{G}_{0}\left(z,\mathbf{k}\right)=\left[z-h^{(A_{1})}\left(\mathbf{k}\right)\lambda^{(A_{1})}\right]^{-1}$ and $\hat{V}=\sum_{\Gamma\neq A_{1}}h^{\left(\Gamma\right)}\left(\mathbf{k}\right)\lambda^{\left(\Gamma\right)}$.

We now consider the sum over crystal momenta of the correlator, $\sum_{\mathbf{k}}G_{\alpha\beta}(\mathbf{k},z)=G_{\alpha\beta}(\mathbf{R}_{i},\mathbf{R}_{i},z)$, which is identical to its previously introduced real space version. By inspecting Eq. (9), we find that $\sum_{\mathbf{k}}G_{\alpha\beta}(z,\mathbf{k})=0$ for $\alpha\neq\beta$, where each term vanishes individually. This finding is traced back to the fact that either (I) the resulting matrix in orbital space is diagonal, or (II) the sum over crystal momenta of an non-trivially transforming function vanishes, \sum_{\mathbf{k}}G^{\left(\Gamma\right)}\left(\mathbf{k}\right)\big{|}_{\Gamma\neq A_{1}}=0. In particular, the first term vanishes because of (I); the second because of (II), $\sum_{\mathbf{k}}[h^{(A_{1})}(\mathbf{k})]^{2}h^{(\Gamma\neq A_{1})}\left(\mathbf{k}\right)=0$; the third because of (I) for $\Gamma=\Gamma^{\prime}$ and because of (II) for $\Gamma\neq\Gamma^{\prime}$, \sum_{\gamma\mathbf{k}}h^{(\Gamma)}\left(\mathbf{k}\right)h^{(\Gamma^{\prime})}\left(\mathbf{k}\right)\lambda_{\alpha\gamma}^{(\Gamma)}\lambda_{\gamma\beta}^{(\Gamma^{\prime})}\big{|}_{\Gamma\neq\Gamma^{\prime}}=0. The identical reasoning can be applied to any higher order term by noticing that if the orbital matrix is non-diagonal the product of $h$’s contains an odd number of non-trivially transforming functions. This eventually implies Eq. (4).

In the case of finite spin orbit coupling, when the spin degrees of freedom cease to be good quantum numbers, the result is potentially altered. In the presence of inversion symmetry, which holds true for electronic bulk states in the tetragonal crystal system, all bands are still doubly degenerate and the Hamilton operator can be rewritten in terms of bilinear form transforming according to one-dimensional representations. However, this potentially changes when the inversion symmetry is broken, which is especially the case when electronic surface states being probed in STM measurements. Here, the surface state band’s degeneracy may be lifted. Furthermore, the orbital character of electronic states of the STM tip has to be taken into account. The average over spin states in Eq. (2) is replace by a weighted average depending on the STM tip’s orbital character and a cancelation of contributions is not guaranteed. For this more complex scenario, the previously drawn conclusion manifesting in Eq. (1) loses its rigidity and orbital interference terms may appear as function of the strength of spin-orbit interaction.