Inhomogeneous Knight shift in vortex cores of superconducting FeSe

The iron-based superconductors have attracted broad interest recently because they can host Majorana modes on the surface, at domain walls, and within vortex cores [1, 2, 3, 4]. Fe(Se,Te), and Li(Fe,Co)As contain bands with $p_{z}$ and $d_{xz}/d_{yz}$ character with non-trivial topologies, that give rise to both topological surface states as well as a bulk Dirac point near the $\Gamma$ point in k-space [5]. FeSe, although topologically trivial, is a particularly interesting case because the superconducting state emerges from a nematic phase that develops below $T_{nem}=91$K [6]. Moreover, the Fermi energy, $E_{F}$, in this system is usually small, such that this system lies close to the BCS-BEC crossover regime [7, 8]. Evidence has emerged that suggests FeSe exhibits a Fulde-Ferell-Larkin-Ovchinnikov (FFLO) phase at high magnetic fields [9]. The possibility of both FFLO and dispersive Majorana modes underlies the importance of a detailed understanding of the nature of the vortices in these materials.

To probe vortex matter it is important to first understand the underlying superconducting state. The spatial part of the superconducting wavefunction in FeSe is generally assumed to be either $s_{\pm}$ or $d$-wave. The nematic normal state gives rise to twin domains and in-plane anisotropy, and the Fermi surfaces contain different orbital characters in the two domains. The superconducting gap, $\Delta$, appears to correlate with the orbital content on the Fermi surface [10]. However, the presence of the domains may mask intrinsic properties about the density of states below $T_{c}$, and there are conflicting reports about the presence or lack of nodes and anisotropy of the superconducting gap function [11, 12, 13, 14].

Information about the spin component of the wavefunction can be gleaned from nuclear magnetic resonance (NMR) Knight shift measurements. The spin susceptibility of a condensate with singlet pairing vanishes, whereas that with triplet pairing can remain unchanged through $T_{c}$. Conventional magnetometry cannot discern these changes because the spin component is much smaller than the orbital component, however the Knight shift is usually dominated by the former and is thus one of the only experimental probes of the spin susceptibility of the condensate [7, 8]. To date, Knight shift measurements in the superconducting state have been inconclusive, revealing little or no change below $T_{c}$ [15, 16, 17, 18]. A recent study reported no change in the Knight shift along the $c$ axis in fields up to 16 T, which have been interpreted as evidence for highly spin-polarized Fermi liquid in the BCS-BEC regime [19]. A lack of suppression of the Knight shift may be evidence for spin-triplet pairing [20], but may also reflect thermal instability of the sample due to eddy-current heating from radiofrequency pulses [21]. In fact, spin-orbit coupling can give rise naturally to a spin-triplet component [22, 23]. To fully characterize the symmetry of the condensate, therefore, it is important to understand the full tensor nature of the Knight shift in the superconducting state.

Here we report ⁷⁷Se NMR on a high quality single crystal as a function of temperature and field. We find that between 3.6 and 11.7 T, the spin part of the planar Knight shift is reduced by $\sim 10-15\%$ from their normal state values below $T_{c}$, whereas the out-of-plane component shows no change within the experimental resolution. These results are consistent with spin singlet pairing in the presence of large orbital susceptibility and spin-orbit coupling. Surprisingly, the NMR linewidths broaden inhomogenously by more than a factor of two below $T_{c}$ for planar fields, but not for $\mathbf{H}_{0}\parallel\mathbf{c}$. Accompanying this broadening is a frequency-dependent spin-lattice relaxation rate, $T_{1}^{-1}$, that reveals electronic inhomogeneity in the superconducting state. This inhomogeneity cannot be explained by the presence of a conventional vortex lattice, but may reflect an enhanced local density of states within the vortex cores.

Single crystals of FeSe were grown by vapor transport with a tilted two-temperature zone tube furnace [24]. Several samples were characterized with magnetic susceptibility and resistivity measurements, with the best samples having $T_{c}=8.9$ K, and RRR defined as the resistance ratio between 250 and 10 K around 19, similar to reported high-quality samples [24]. A large crystal of dimensions 2.4 mm $\times$1.4 mm $\times$ 0.2 mm was selected and mounted in a custom-built NMR probe equipped with a dual-axis goniometer. The majority of the experiments were conducted within a variable-temperature cryostat in a high-homogeneity NMR magnet with a field of $H_{0}=11.7294$ T, and some experiments at lower fields were conducted in a PPMS system. In this field, $T_{c}$ is suppressed to $\sim 5.3$ K (measured by resistivity) [25]. Spectra (Fig. 1) were collected for field aligned along the tetragonal [110], [100] and [001] directions. The [110] ([1$\bar{1}$0]) direction corresponds to the Fe-Fe bond, and is the $\mathbf{a}$ ($\mathbf{b}$) direction in the nematic phase [26]. The spectra were measured at several temperatures down to 2.1 K using low-power rf pulses ($\pi/2$-pulse widths up to $80\mu$s), sweeping frequency and summing the Fourier transforms. Our results are consistent with previous reports [15, 17, 16, 27, 18], and reveal a splitting of the single ⁷⁷Se resonance below $T_{nem}=91$ K due to twinning. The resonance frequencies are given by $f=\gamma H_{0}(1+K)$, where $\gamma=8.118$ MHz/T is the gyromagnetic ratio and $K$ is the Knight shift. We fit each resonance to a Gaussian function, and Figs. 2(a,b) shows the temperature dependence of $K$ and the full-width half-maxima, FWHM, for several different field directions. Below $T_{c}$, the spectra exhibited a strong dependence on the pulse power, as illustrated in Fig. 1(c,d). The radiofrequency pulses induce eddy currents around the sample, which can lead to Joule heating. As a result, the temperature may temporarily exceed $T_{c}$ immediately after the pulse. Similar effects have been observed in other superconductors, leading to misinterpretations about the temperature dependence of the Knight shift [21]. The shifts reported in Fig. 2(b) were measured at 18.4 dBm, where there was no power-dependence to the spectra.

The Knight shift arises from the hyperfine interaction between the nuclear spin and the spin and orbital degrees of freedom of the electrons: $\mathcal{H}_{hf}=\mathbf{I}\cdot\mathbb{A}_{S}\cdot\mathbf{S}+\mathbf{I}\cdot\mathbb{A}_{L}\cdot\mathbf{L}$, where $\mathbb{A}_{S,L}$ are the hyperfine coupling tensors, and $\mathbf{S}$ and $\mathbf{L}$ are the spin and orbital angular momenta. The Knight shift is given by $K_{\alpha}=A^{L}_{\alpha\alpha}\chi_{\alpha\alpha}^{orb}+A^{S}_{\alpha\alpha}\chi_{\alpha\alpha}^{spin}+(A^{S}_{\alpha\alpha}+A^{L}_{\alpha\alpha})\chi_{\alpha\alpha}^{mixed}$, where $\chi_{\alpha\alpha}^{spin}$, $\chi_{\alpha\alpha}^{orb}$, and $\chi_{\alpha\alpha}^{mixed}$ are the static spin, orbital, and mixed susceptibilities at zero wavevector [28, 26]. In the absence of spin-orbit coupling, the mixed term vanishes and the Knight shift is usually decomposed as $K_{\alpha}=K_{\alpha 0}+A^{S}_{\alpha\alpha}\chi_{\alpha\alpha}^{spin}$. $K_{\alpha 0}$ is often considered to be a temperature-independent shift arising from a Van-Vleck orbital susceptibility, however, this decomposition breaks down in the presence of spin-orbit coupling [29]. Moreover, theoretical calculations have revealed that $\chi_{\alpha\alpha}^{orb}\gg\chi_{\alpha\alpha}^{spin},\chi_{\alpha\alpha}^{mixed}$ due to the multiorbital nature of the band structure and nematic instability [26]. As a result, the relationship between $K_{\alpha}$ and the bulk susceptibility, $\chi=\chi^{spin}+\chi^{orb}+2\chi^{mixed}$, is complicated. Nevertheless, we find that $K_{\alpha}$ varies linearly with $\chi$ above $T_{nem}$, as shown in the inset of Fig. 3. Linear fits to the data yield parameters close to previously reported values [18].

For spin-singlet pairing, $\chi^{spin}$ and $\chi^{mixed}$ should vanish in the superconducting state, giving rise to a suppression of $K$ below $T_{c}$, as observed in Fig. 3. For planar fields, $K_{a,b}$ is suppressed by about $100\pm 15$ ppm in both domains, as well as for the [100] direction oriented $45^{\circ}$ to the Fe-Fe bond direction. This magnitude of suppression does not change significantly at lower fields. For out of plane fields, any change in $K_{c}$ is within the noise, but is less than $\sim 10$ ppm. These results are also independent of applied field, and are consistent with previous reports [30, 19]. Note that $K_{\alpha}(T\rightarrow 0)\neq K_{\alpha 0}$, or in other words the low temperature limit of the shift does not equal the intercepts from the $K-\chi$ plot. In fact, $\chi^{orb}$ is strongly temperature dependent, so $K_{\alpha 0}$ does not represent a temperature independent Van-Vleck term. The low temperature shift reflects a finite $\chi^{orb}$, since the spin component vanishes for singlet pairing; however impurity states may play a role [31]. Determining how much $\chi^{spin}$, $\chi^{orb}$ and $\chi^{mixed}$ are suppressed below $T_{c}$ will likely require detailed theoretical calculations [26]. It is noteworthy that the difference $K_{a}-K_{b}$, shown in Fig. 2(c), exhibits a subtle enhancement below $T_{c}$. This observation suggests that the superconductivity is slightly anisotropic in the two domains, and may reflect an anisotropy in the coherence lengths, $\xi_{a,b}$.

Below $T_{c}$ the spectra for both domains broaden and become asymmetric with a high frequency tail, as observed in Fig. 1(c,d) and 2(a). At room temperature, the FWHM of the spectrum $(\sim 0.08$ kHz) is close to the second moment of the nuclear spin dipole moments of the lattice ($\sim 0.06$ kHz). The excess inhomogeneous broadening above $T_{c}$ may be due to either macroscopic or microscopic strain fields [32]. The crystals were initially secured to the goniometer with a light coat of superglue and the coil fit loosely around the sample. In this case, the observed linewidths were smaller (open points in Fig. 2(a)), reflecting the high quality of this crystal. In later measurements, the crystal was remounted and we observed the linewidth increase by a factor $\sim 2$ (solid points in Fig. 2(a)). It is possible that remounting the crystal introduced inhomogeneous macroscopic strain fields. Non-magnetic impurities such as Fe vacancies are known to exist in the lattice [10], which may also be a source of local strain and inhomogeneous broadening in the normal state.

Regardless of the linewidth in the normal state, an even larger increase of linewidth is observed below $T_{c}$, which is unexpected. A vortex lattice certainly gives rise to a distribution of local magnetic fields, $B(\mathbf{r})$, and in very low fields $B\ll B_{c2}$, the second moment of the field distribution can be estimated as $\Delta B^{2}\approx 0.00371\phi_{0}^{2}\lambda^{-4}$, where $\phi_{0}$ is the flux quantum, and the penetration depths are $\lambda_{a}=446$ nm, $\lambda_{c}=1320$ nm [33, 34]. There are important corrections to this expression in the higher fields where our measurements were conducted [35], however after accounting for these we estimate that the normal state spectra should broaden by only $\sim 8$ Hz, three orders of magnitude smaller than the enhancement observed in Fig. 2(a) [36].

Since the field distribution alone is unable to capture the asymmetric broadening, we hypothesize the presence of a spatially-varying Knight shift, $K_{\alpha}(\mathbf{r})$, that is equal to the normal state value within the vortex cores and decays to $K_{\alpha}(T\rightarrow 0)$ outside. The spectrum is given by the histogram of the local resonance frequency, $f(\mathbf{r})=\gamma B(\mathbf{r})(1+K_{\alpha}(\mathbf{r}))$. The exact shape of the spectrum depends on microscopic details, but if the spatial variation $\delta K$ is equal to the 100 ppm suppression observed in Fig. 3, the spectrum will broaden by $\sim 10$ kHz, which agrees well with the excess linewidth below $T_{c}$ in Fig. 2(a). These results suggest that the local spin susceptibility within the vortex cores is identical to that in the normal state.

This interpretation is supported by $T_{1}^{-1}$ measurements. Fig. 4(a) shows $(T_{1}T)^{-1}$ versus temperature. The data in the normal state agree well with published results [15, 37, 18]. This quantity drops due to the superconductivity, and becomes inhomogeneous in the mixed phase. Fig. 4(b) shows that $T_{1}^{-1}$ increases by nearly a factor of two in the high frequency tails of the spectra in the superconducting state, which correspond with the vortex cores. Localized Caroli-deGennes-Matricon (CdGM) electronic states normally exist within isolated cores [38]. At higher fields quasiparticles from different cores can propagate coherently across multiple vortices, and the energy spectrum becomes dispersive, with gapless excitations remaining within the vortex cores that give rise to a finite local density of states (LDOS) which should be manifest in any technique sensitive to low energy excitations [39, 40, 41, 42]. Indeed NMR studies have identified enhanced spin-lattice-relaxation rate within the vortex cores of both conventional [43] and unconventional superconductors [44, 45].

There are, however, important differences between FeSe and previous observations on other superconductors. In the cuprates, the excess relaxation rate has been attributed to antiferromagnetic fluctuations from a competing ground state to superconductivity [44, 45], as well as from Doppler-shifted quasiparticles associated with d-wave nodes [46]. In such cases $(T_{1}T)^{-1}$ exhibits a strong Curie-Weiss divergence within the cores, whereas outside the cores $(T_{1}T)^{-1}$ remains temperature independent. In the s-wave superconductor LaRu₄P₁₂, $(T_{1}T)^{-1}$ in the cores is also strongly temperature dependent, and even exceeds the value in the normal state [43]. In the case of FeSe, $(T_{1}T)^{-1}$ exhibits Curie-Weiss behavior in the normal state (dotted lines in Fig. 4a), but drops below $T_{c}$. This behavior has been attributed to antiferromagnetic spin fluctuations that are gapped by the superconductivity [47, 48]. The open circles (squares) in Fig. 4(a,c) show the temperature dependence at the upper end of the spectra in the superconducting state. $(T_{1}T)^{-1}$ in the vortex cores changes only by a factor of two from the background rate, remaining well below the normal state value, and exhibits the same trend with temperature as the background. These results suggest the absence of any spin fluctuations within the normal cores of FeSe.

FeSe appears unique in that there is a $\mathbf{q}=0$ spin response in the vortex cores. It is unclear whether this behavior could be related to either a proximity to the BCS-BEC crossover, or either a Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) phase [8, 9] or a field-induced spin density wave [49] for parallel fields above $H^{*}=24$ T. A true FFLO phase should exhibit both segmented vortex lines and normal planes where the LDOS reaches the normal state values, giving rise to inhomogeneously broadened NMR spectra. Although $H_{0}\sim 0.5H^{*}$ in our experiments, the inhomogeneity we observe already indicates the presence of large spin polarization in spatial regions where the superconducting order vanishes. It is noteworthy that the $\mathbf{q}=0$ susceptibility in FeSe is dominated by orbital contributions, whereas the finite $\mathbf{q}$ response is dominated by spin fluctuations [26]. Condensation of singlet pairs enables us to probe the small spin response at $\mathbf{q}=0$. In Ba(Fe,Co)₂As₂, finite $\mathbf{q}$ spin fluctuations can freeze and exhibit long range antiferromagnetism in vortex cores [50]. In FeSe we see no evidence for such behavior, which may be due to the presence of nematic order and the different contribution of orbital versus spin susceptibility. The absence of such fluctuations suggests that the high field phase is unrelated to SDW order [49, 51].

In summary, we find that the Knight shift is suppressed below $T_{c}$ for in-plane fields, but see little to no suppression for field along the $c$-axis. The spectra are inhomogeneously broadened below $T_{c}$, and $T_{1}^{-1}$ becomes frequency-dependent. These observations are consistent with a finite LDOS within the vortex cores. We find no evidence of competing antiferromagnetic fluctuations in the vortex cores. Further studies at higher fields or with Te doping should shed light on the unusual nature of the vortex states in this system.

Acknowledgment. We acknowledge helpful discussions with B. Andersen, E. da Silva Neto, M. Walker, and R. Fernandes, and thank P. Klavins for assistance in the lab. Work at UC Davis was supported by the NSF under Grants No. DMR-1807889, and synthesis of single crystals was supported by the UC-Lab fees program.