Direct visualization of the nematic superconductivity in Cu$_x$Bi$_2$Se$_3$

Cu$_x$Bi$_2$Se$_3$ hosts both topological surface states and bulk superconductivity. It has been identified recently as a topological superconductor (TSC) with an extraordinary nematic, i. e. C2-symmetric, superconducting state and odd-parity pairing. Here, using scanning tunneling microscopy (STM), we directly examine the response of the superconductivity of CuxBi2Se3 to magnetic field. Under out-of-plane fields, we discover elongated magnetic vortices hosting zero-bias conductance peaks consistent with the Majorana bound states expected in a TSC. Under in-plane fields, the average superconducting gap exhibits two-fold symmetry with field orientation; the long C2 symmetry axes are pinned to the dihedral mirror planes under B//=0.5 T but rotate slightly under B//=1.0 T. Moreover, a nodeless {\Delta}4x gap structure is semi-quantitatively determined for the first time. Our data paint a microscopic picture of the nematic superconductivity in CuxBi2Se3 and pose strong constraints on theory.


I. INTRODUCTION
Topological superconductors, as host of Majorana fermions and Majorana zero modes (MZM) [1], may facilitate topological quantum computing. However, in practice, TSCs are rare. Among various recipes for making a TSC [1,2], bulk superconductors that host topological surface states are the most natural candidates. However, unlike the theoretical triumph in predicting topological insulators [1,3], metals, and semi-metals, theory and experiment diverge when it comes to TSCs. Consequently, there is a large gap in our understanding of TSCs, and in particular, details of the microscopic behavior of TSCs are in urgent demand for improving theory.
CuxBi2Se3 is a prototypical example of the difficulties encountered. Since it hosts both topologically non-trivial surface states and bulk superconductivity [4,5], it has been proposed to be a TSC, likely with an odd-parity pairing symmetry [6]. However, experimental results on CuxBi2Se3 present challenges for this interpretation. The absence of Pauli limiting behavior in the upper critical field suggests spin-triplet superconductivity or an anisotropic spin-singlet state [7], and point-contact spectroscopy found a zero-bias conductance peak (ZBCP) on a cleaved surface [8], which was attributed to MZM. However, Andreev reflection spectroscopy on CuxBi2Se3 shows that the existence of the ZBCP depends on the barrier strength -its absence under finite barrier strength implies the absence of zero-energy Majorana fermions [9]. Furthermore, a low-temperature STM study observed a full superconducting gap without in-gap states, nor any ZBCP in the vortex core, discrediting CuxBi2Se3 as a TSC [10].
Recently, a nematic superconducting state was discovered in an NMR study of CuxBi2Se3 (Ref. [11]). The spin susceptibility in the superconducting state exhibits two-fold symmetry when rotating the in-plane magnetic field, which is surprising considering the three-fold symmetry of the lattice. This, together with the invariance of the Knight shift upon crossing Tc under an out-of-plane field, indicates the Cooper pairs are in a pseudospin-triplet state with a pinned d-vector direction [11]. Subsequently, measurements on Cu-, Sr-, and Nb-doped Bi2Se3 superconductors all revealed two-fold symmetry in the in-plane field angular dependence of the specific heat, upper critical field, critical current, and magnetoresistance below Tc (Refs. [12][13][14][15]). These remarkable observations in CuxBi2Se3 have established it as an odd-parity TSC state [16][17][18], with the d-vector pinned to a specific lattice axis by strong spinorbit coupling. However, contradictions still exist, especially in the gap symmetry and the pinning directions of the d-vectors among the reported bulk and/or macroscopic measurements [11][12][13][14][15]. The d-vector is reported to be pinned along the Se-Se bond direction in an NMR study by Matano et al. [11], but it lies along the dihedral mirror plane of the Se lattice in specific heat measurements by Yonezawa et al. [13] and magnetoresistance measurements by Pan et al. [12]. Thus, to reach a comprehensive microscopic understanding of this unique state, and to reconcile the previous contradicting reports, we investigate the superconducting properties of CuxBi2Se3 by STM at ultra-low temperature under magnetic fields with variable direction.

II. EXPERIMENTAL METHODS
The CuxBi2Se3 (x=0.31) single crystals were prepared via an electrochemical intercalation method as described in Ref. [19]. Magnetic susceptibility measurement shows a superconducting transition temperature (Tc) of about 3 K and shielding fraction of about 17% at 1.8 K (See part 1A of Ref. [20]). STM measurements were conducted in a millikelvin STM system with vector magnetic field, at the base temperature of 20 mK. The effective electron temperature (Teff) of the system is calibrated to be 310 mK (see part 2 of Ref. [20]). The CuxBi2Se3 samples were cleaved in vacuum at 77 K and immediately transferred into the STM module. PtIr STM tips were used after being treated on an Au (111) surface. The dI/dV spectra are collected using a standard lock-in technique with modulation frequency f = 787 Hz and typical modulation amplitude ΔV = 30-50 µV.
CuxBi2Se3 crystals are thin-bar shaped with the sizes of about 1.7×1.3×0.30 mm 3 for sample 1 and 2.5×1.7×0.14 mm 3 for sample 2. The demagnetization effect due to the irregular sample shape is very weak as discussed in part 1B of Ref. [20], thus it is neglected in our experiments under magnetic fields. For the measurements under B//, the possible residual out-of-plane field component is carefully eliminated via vortex mapping (see part 3 of Ref. [20]). Moreover, the influence of in-plane magnetic vortex and spatially inhomogeneity of superconducting gap under B// were considered and discussed in detail in part 4 of Ref. [20].

III. RESULTS AND DISCUSSION
The cleaved CuxBi2Se3 surface exhibits three types of regions ( Fig. 1), two superconducting (SC-I, SC-II) and one non-superconducting (NSC). The statistics on how often these regions are observed in different samples are summarized in Table I of Ref. [20]. The probability of SC regions being found is unexpectedly low, indicative of the intrinsically inhomogeneous superconductivity of CuxBi2Se3 [21].  Fig.  1(e). This height corresponds to the spacing between adjacent Bi-Se quintuple layers based on previous x-ray diffraction of CuxBi2Se3 (Ref. [4]). The terraces are atomically flat with two kinds of intrinsic defects ( Fig. 1(b)). Bright dots are most likely intercalated Cu atoms as seen in previous STM studies [10,22], while trefoil defects, commonly observed on Bi2Se3 surfaces, arise from the substitution of Se atoms by Bi (Ref. [23]). The inset of Fig. 1(b) shows the undistorted hexagonal Se atomic lattice with the lattice constant of about 0.41 nm. Here we define the x-axis as one of the Se-Se bond directions, and the y-axis as perpendicular to x within the same Se-plane, as indicated in the inset of Fig. 1(b). As for the electronic structure, a Dirac cone-like feature is suggested by linearly dispersing density of states (DOS)) observed in the dI/dV spectra shown in Fig. 1(d) [23,24], as well as the quasiparticle interference patterns along a step edge shown in part 5 of Ref. [20]. The extracted Dirac point is located at about 370 meV below the Fermi energy (EF). This is consistent with ARPES measurements on bulk CuxBi2Se3 crystals [5], indicating the existence of topologically non-trivial surface states in NSC regions. A nearly flat DOS is observed around EF ( Fig. 1(c)), precluding the existence of a superconducting gap.
Figures 1(f) and 1(g) display the surface morphologies of one kind of superconducting region (SC-I) observed in sample 1. Large terraces are observed here as well, but their heights are integer multiples of 1.56 nm ( Fig. 1(j)), much larger than in the NSC regions. Moreover, as shown in Fig. 1(g), the terraces are rough, with many Cu clusters distributed. It appears that the spacing between Bi2Se3 quintuple layers in this region is expanded considerably by Cu intercalation, which would lead to a more two-dimensional structure. The SC-I region is adjacent to the NSC region shown in Figs. 1(a) and 1(b). Although the surface lattice cannot be directly seen in SC-I region, the orientation of atomic step edges suggests it has the same lattice orientation with the nearby NSC region (see part 1D of Ref. [20] for details). The dI/dV spectrum obtained in SC-I regions ( Fig. 1(i)) exhibits a similar Dirac-cone like feature, but the Dirac point is now 650 meV below EF. This significant shift to lower energy compared to the NSC region indicates that the SC-I region is heavily doped with electrons. As shown in Fig. 1(h), a fully developed superconducting gap with pronounced coherence peaks is observed at EF, which is consistent with previous reports [10,19]. We note that the heavy electron doping effect may be critical for the origin of superconductivity here. Such a U-shaped gap is spatially homogeneous in SC-I region, and it is insensitive to Cu clusters and step edges, as shown in Figs. 1(k) and S6 in Ref. [20]. Usually for spin-triplet superconductors, such as Sr2RuO4, the superconductivity is very sensitive to impurities [25]. However, two theoretical studies suggest a different impurity scattering behavior for the unconventional superconductivity resided in a topological insulator, where the pair decoherence introduced by impurity scattering is strongly suppressed by spin-orbital coupling and spin helical nature in intercalated Bi2Se3 [26,27].
Besides SC-I regions, we observe another type of superconducting region (referred to here as SC-II) in sample 2, on islands with steep topographic variations as shown in Figs. 1(l) and 1(m). While the NSC and SC-I regions often exhibit large flat terraces with uniform terrace height, SC-II regions exhibit many irregular terraces of several to tens of nanometers width; the terrace height is non-uniform, even with fractional values observed, whose origin are discussed in part 1E of Ref. [20]. However, the terraces are atomically flat -with the atomic lattice shown in the inset of Fig. 1(m), which is the same as that of the NSC region. No clear Dirac-cone-like feature is observed in this region, while a sharp U-shaped superconducting gap appears at EF, as respectively shown in Figs. 1(o) and 1(n). Why superconductivity occurs in these islands is unclear, and we give some discussions in part 1E of Ref. [20].
The U-shaped spectrum suggests the superconducting gap is fully opened at the Fermi surface, however we find that it still has significant broadening beyond thermal effects, which can be reasonably accounted for by k-space anisotropy of the gap size [11][12][13][14][15]. We fitted the superconducting spectra by the Dynes formula [28] with a twofold anisotropic gap function: Δk=Δ0+Δ1|cosθ| and Teff = 310 mK (see part 6 of Ref. [20]). A small Dynes term Γ is used to account for finite quasi-particle lifetime broadening [28]. The fit to the gap in the SC-I region shown in Fig. 1(h) yields Δ0=0.256 meV, Δ1=0.206 meV and Γ=0.004 meV; the fit to a typical superconducting spectrum in the SC-II region ( Fig. 1(n)) yields Δ0=0.477 meV, Δ1=0.103 meV and Γ=0.001 meV. Evidently, these regions differ greatly in gap size and gap anisotropy ratio. Significant anisotropy of the superconducting gap is expected for the reported nematic superconducting state of CuxBi2Se3 (Refs. [11][12][13][14][15]).
To further examine the superconducting state, we applied an out-of-plane magnetic field and measured the SC-I regions. Figures 2(a) and S9 in Ref. [20] show zero-bias conductance (ZBC) mappings under various B⊥, where typical vortex lattices are revealed. It is notable that the vortices are significantly elongated along the y-axis and exhibit an elliptical shape. Meanwhile, the vortex lattices are also stretched along this direction. Exponential fits to the profiles along the long and short axes of a vortex under B⊥=0.2 T are shown in Fig. 2(b), which give Ginzburg-Landau coherence lengths (ξ) of 30.83 nm and 19.70 nm, respectively, a ratio of 1:0.64. Since the shape of the vortex is related to the superconducting gap structure [29], the elongated vortices observed here are consistent with the highly anisotropic superconducting gap (Δk) suggested by the fit in Fig. 1(h), where the ratio between the gap maximum and gap minimum in the SC-I region is estimated to be (Δ0+Δ1):Δ0=1:0.55. Moreover, upon increasing B⊥, as the vortex density increases, the vortex size and ξ decrease noticeably (Figs. 2(c) and S9(a-c) in Ref. [20]). Such a significant field dependence of the vortex size may suggest multiband effects [30][31][32][33], e.g., the ξ of the surface and bulk states of CuxBi2Se3 could be different. We note that a field-dependent coherence length was also observed in Bi2Te3 films grown on a NbSe2 substrate [34].

Figures 2(d) and 2(e)
show the evolution of dI/dV spectra taken along the short and long axes of a vortex under B⊥=0.2 T, respectively. Far from the vortex, the superconducting gap is only slightly affected (blue curves in Figs. 2(d) and 2(e)), and it is gradually suppressed upon approaching the vortex core. Within a few nanometers of one vortex core center, weak ZBCPs are observed with a full-width at halfmaximum of about 0.29 meV (red curves in Figs. 2(d) and 2(e)). This ZBCP is observed for the first time in vortices of CuxBi2Se3 (Refs. [10,22]). It could be composed of both conventional Caroli-de Gennes-Matricon bound states and the long-predicted Majorana bound state expected at zero energy in 2D TSC. These are hard to distinguish because the energy separation of the bound states δE= Δ 2 /EF (Ref. [35]) is about 0.4 µeV here, assuming Δk=0.5 meV and EF=650 meV. Nevertheless, our observation of a ZBCP is consistent with the prediction of TSC in CuxBi2Se3. As shown in Figs. 2(f), S9(d), S9(e), and S10 in Ref. [20], ZBCPs have also been observed in some vortex cores under B⊥=0.3 T, but their intensities are slightly weaker than under B⊥=0.2 T. The spectrum near zero bias is distorted for B⊥=0.4 T, and the ZBCP is clearly absent for B⊥=0.6 T, as shown in Figs. 2(f) and S9 in Ref. [20]. This behavior resembles the evolution of the Majorana zero-energy mode observed in 5 QL Bi2Te3 films grown on NbSe2 (Ref. [36]). The weakening and disappearance of a ZBCP in the vortex core under high field is interpreted as a result of the coupling between states in adjacent vortices [36]. When the field is small, the distance between vortices is large compared to the vortex core size, so the interaction between the vortex bound states can be neglected. Upon increasing field, the intervortex distance shortens, enhancing these interactions and eventually destroying the Majorana zero mode [36].
To further illustrate the nematic superconductivity in CuxBi2Se3, we apply an inplane magnetic field and study the superconducting spectra within the same area of SC-I region as a function of θ, defined as the azimuthal angle of B// with respect to the x-axis of the Se lattice (see inset of Fig. 3(a)). Under various B//, no obvious in-plane vortices are observed on the surface, and the superconducting spectra are spatially homogeneous (see part 4 of Ref. [20] for more details). This is possibly because the superconductivity on the surface (likely resides in topological surface state) is more robust to the orbital depairing of an in-plane field, so that the bulk vortex lines are pushed away from the surface, as shown theoretically in Ref. [37]. Consequently, the in-plane field did not generate noticeable spatial inhomogeneity on the superconducting spectrum. However, we observed obvious variation of superconducting spectrum with the direction of B//. Figure 3(a) shows typical superconducting spectra when applying a B// of 0.5 T at θ = 90°, 40° and 10°. Clearly, the superconducting gap, as estimated from both the leading edge (Δ * LE) and coherence peaks (Δ * CP), varies with θ. We use asterisks (*) here to emphasize that these actually constitute weighted averages of Δk under B//. The angular dependence of Δ * CP and Δ * LE based on an extensive dataset is summarized in Fig. 3(b), and exhibits clear two-fold symmetry, with maxima around θ = 90° and 270° (coinciding with the y-axis) and minima near θ = 0° and 180° (coinciding with the x-axis). Figures  3(c) and 3(d) show similar two-fold symmetry for B//=1.0 T. In addition, Δ * LE along θ = 85°, 40° and 0° decrease monotonically with increasing B//, and no crossover is observed between them, as shown in Figs. 3(e) and S11 in Ref. [20]. Our detailed experiments indicate that the superconducting gap (Δ * CP and Δ * LE) of CuxBi2Se3 exhibits two-fold symmetry throughout the studied B// range from 0.25 to 1.5 T. The suppression of the superconducting gap by B// is much faster along the x-axis (θ = 0°) than nearly along the y-axis (θ = 85°), suggesting a smaller upper critical field #$ along the x-axis. Thus, a larger ξy is evaluated than ξx from the Ginzburg-Landau relations #$ % = ' ( $)* + * , , which is consistent with our observations that the vortex is elongated along the y-axis of the lattice shown in Fig. 2(a). Figure  3(f) displays the oscillation amplitudes of Δ * CP and Δ * LE under various B//. The ratios are different for Δ * CP and Δ * LE, but the trend is the same -both increase gradually as B// increases. For Δ * LE, the oscillation amplitude increases from 8.7% under B// = 0.25 T to 28.9% under B// = 1.5 T, while for Δ * CP it changes from 5.1% to 8.1%. The gradual smearing of the gap edge (Δ * LE) of the DOS under B// is consistent with the Volovik effect -the quasiparticle energy spectrum is Doppler shifted in the mixed state due to a supercurrent flow running around the vortices, which leads to a broadening of the gap edge even for a nodeless superconductor at finite temperature [38][39][40][41][42].
Our main results are summarized in Fig. 4. The data are plotted in polar coordinates and with respect to the underlying Se lattice. Figure 4(a) shows the profile of an elongated vortex core and the two-fold symmetry of Δ * LE measured under B//=0.5 T for the SC-I region. The long axis of Δ * LE is shown by the red dashed line, which actually nearly coincides with the long axis of the vortex, as well as the dihedral mirror plane (dash-dot lines) of the Se lattice. Figure 4(b) shows the angular distribution of Δ * LE under B//=1.0 T -its long axis is rotated by about 17.5° off the dihedral mirror plane. This can also be readily seen by comparing the data in Figs. 3(b) and 3(d) directly, as in Fig. S14(a) in Ref. [20].
One can retrieve information on the angular distribution of Δk based on the angular dependence of Δ * LE. Since the orbital depairing effect of B// (in-plane vortex) is not obvious here, the influence of B// on superconductivity mainly includes two effects: the Volovik [38][39][40][41][42] and Zeeman effects. The former has been extensively studied in nodal superconductors theoretically and experimentally [38][39][40], and leads to a DOS under B// that exhibits characteristic oscillatory behavior with field orientation. At low temperature and low field, as illustrated in Fig. 4(e), more quasiparticles are excited when the field is perpendicular to the gap-minimal direction of Δk than along other directions, because the field-induced supercurrents (which run perpendicular to the field) can lead to the strongest Doppler shift of the quasiparticle energy around the gap minimal region. Thus, the short axis of Δ * LE is actually perpendicular to the gap minimal direction, since this field orientation produces the most quasiparticle excitations. Meanwhile, the application of B// will induce Zeeman effects, whose energy scale is about 0.1 meV for B// = 1.0 T, comparable to the minimal gap magnitude (0.256 meV) in the SC-I region. For Eu-symmetry pairing, the d-vector is expected to be parallel to the maximal gap direction of Δk, and it is normal to the plane in which the spins of the equal-spin-paired electrons are aligned [16]. The Zeeman depairing effect is the strongest when B// is parallel to the d-vector, because the electron spins may deviate from the original plane under sufficient B// [25]. Therefore, both the Volovik effect and Zeeman effects will suppress Δ * LE much more strongly when B// lies along the maximal gap direction of Δk. Consequently, the distribution of Δk for the SC-I region can be qualitatively determined -as illustrated by the thick blue curve in Fig. 4(a), the gap maxima are along a Se-Se bond direction, which is rotated 90° from the long axis of Δ * LE. The tunneling spectra here and the specific heat measured by Yonezawa et al. [13] reflect the integrated DOS over the Fermi surface under the influence of B//, and they thus exhibit weaker anisotropy than Δk, whose anisotropy can be more precisely estimated by fitting the tunneling spectra or via the vortex anisotropy. The blue curve thus gives the first semi-quantitative sketch of an anisotropic nodeless gap structure in CuxBi2Se3. This angular distribution of the Δk is consistent with the orientation of the observed elongated vortices as, empirically, a vortex will extend out in the direction of a gap minimum (assuming only weak Fermi velocity anisotropy) [29].
For the SC-II region, qualitatively similar data were measured on two domains (see Figs. S12 and S13 in Ref. [20]). Δ * LE and Δ * CP in this region also vary noticeably with field orientation θ, exhibiting two-fold symmetry, although the Δk here are larger than those of the SC-I regions. Figure 4(c) gives the angular dependence of Δ * LE in domain A of the SC-II region under B//=0.5 T -the long axis of Δ * LE lies along the dihedral mirror plane of the Se lattice. According to above analysis, the angular distribution of Δk in domain A is illustrated by the blue curve, with its maxima again along a Se-Se bond direction. Similar to the case of SC-I region, in domain B of the SC-II region, which has the same lattice orientation as domain A, the long axis of Δ * LE under B//=1.0 T is rotated to an intermediate angle, about 20.7° off another dihedral mirror plane, as shown in Figs. 4(d) and S14(b) in Ref. [20]. The origin of such a rotation under B//=1.0 T is unclear and merits future investigation. Our data show that gap nematicity is ubiquitous in all the superconducting domains of CuxBi2Se3, the long C2 axis of Δk is pinned along one of the three equivalent Se-Se bond directions, which may vary for different SC domains (see more discussion in part 9 of Ref. [20]).
In the theory of Fu et al. describing the nematic superconductivity in CuxBi2Se3 (Ref. [16]), two types of odd-parity pairing symmetries were proposed, Δ4y and Δ4x. Δ4y is nodeless with its gap maxima along a dihedral mirror plane of the Se lattice (yaxis), while Δ4x has its gap maxima along a Se-Se bond direction (x-axis) with a pair of point nodes in the yz-plane that are protected by mirror symmetry. The Δ4y gap distribution is supported by specific heat [13] and magnetoresistance measurements [12]. However, our data suggest that the gap maxima of Δk are pinned along one of the three equivalent Se-Se bond directions, consistent with the NMR study [11], and supporting a Δ4x pairing symmetry without nodes. The differences among bulk measurements reminds us that the summed responses from all domains with different gap amplitudes and orientations may give an average gap maximum rotated from the actual ones, as the multi-domain effect is directly observed in our study on SC-II region (part 9 of Ref. [20]) and previously discussed by Yonezawa et al. in Fig. S8 of their paper [13]. Thus, our findings at the atomic scale facilitate an understanding of the complex phenomena reported in bulk measurements on doped Bi2Se3 materials [11][12][13][14][15], and provide more-direct constraints on theory.

IV. CONCLUSION
In summary, our data have revealed a unique nematic superconducting state in CuxBi2Se3, and give a semi-quantitative description of its anisotropic gap distribution, which is likely a nodeless Δ4x pairing symmetry. On the one hand, we have observed a Dirac-cone-like surface DOS in the SC-I regions, a ZBCP in vortex cores under 0.2~0.3 T out-of-plane fields, and a nematic superconducting order parameter with an anisotropic gap and a preferential symmetry axis. A nematic superconducting state with an anisotropic gap has only been predicted for odd-parity superconductivity [16][17][18], and the preferential direction of the C2 symmetry axis of the gap and its highly nontrivial response to field strength indeed suggest a pseudo-spin triplet pairing [11]. The ZBCP might be the long-sought Majorana bound state, so our findings further support topologically non-trivial superconductivity in this material. On the other hand, many aspects of our data are not readily explained by current theories for TSC. For example, we neither observed in-gap states (Figs. 1(k) and S15 in Ref. [20]) at the step edge as would be expected for a two-dimensional TSC [8,16,43], nor did we observe in-gap states corresponding to Majorana Fermions on the surface as expected for a three-dimensional TSC [16,18]. It is still mysterious why the predicted point nodes appear to be gapped out, and why the long axis of Δ * LE rotates under high B//. In short, our findings reveal the microscopic behavior of the nematic superconductivity in CuxBi2Se3, which will facilitate an understanding of its extraordinary properties and microscopic pairing mechanism, and more importantly, our data impose strong constraints for further improving theories of topological superconductivity. Red squares are the experimental data, while the magenta symbols in (d) denote the symmetric data obtained by shifting the experimental data by 180°. Red dashed lines show the long axis of Δ * LE through its two gap maxima. Polar coordinates are plotted out by gray dashed circles, on a scale that differs among the panels; the sinusoidal fit functions describing the two-fold symmetry of Δ * LE are also listed to show the oscillation amplitudes and the directions of the C2 axes. The rotation of the long C2 axis of Δ * LE under B//=1.0 T is marked in (b) and (d). The angular distribution of Δk in various superconducting regions is sketched by thick blue curves, whose anisotropies are extracted from the fitted gap anisotropy in Figs. 1(h), 1(n) and S8. The fitting functions for Δk are also listed. (e) Schematics of the in-plane field angular dependent quasiparticle excitation (red dots) in the low-temperature and low-field condition [39].