# Theory of vortex lattice effects on STM spectra in d-wave superconductors

## Abstract

Theory of scanning tunneling spectroscopy of low energy quasiparticle (QP) states in vortex lattices of d-wave superconductors is developed taking account of the effects caused by an extremely large extension of QP wavefunctions in the nodal directions and the band structure in the QP spectrum. The oscillatory structures in STM spectra, which correspond to van Hove singularities are analysed. Theoretical calculations carried out for finite temperatures and scattering rates are compared with recent experimental data for high-$`T_c`$ cuprates.

The electronic structure of the mixed state in d-wave superconductors reveals a number of fundamentally new features (see  and references therein) as compared to the case of s-wave compounds, where low lying quasiparticle (QP) states are bound to the vortex core and are weekly perturbed by the presence of neibouring vortices at magnetic fields $`HH_{c2}`$. The vanishing pair potential in the nodal directions results in the extremely large extension of QP wavefunctions which are sensitive to the superfluid velocity ($`𝐕_s`$) fields of all vortices and, thus, the electronic structure is influenced by the vortex lattice geometry. The resulting peculiarities of the local density of states (DOS) can be detected, e.g., by a scanning tunneling microscope (STM). In this paper we focus on the theory of scanning tunneling spectroscopy of low energy QP states in vortex lattices of d-wave superconductors and compare the theoretical calculations with recent experimental data  for high-$`T_c`$ cuprates, where the dominating order parameter is believed to be of d-wave symmetry. Hereafter we assume the Fermi surface (FS) to be two-dimensional (2D), take the gap function in the form $`\mathrm{\Delta }_𝐤=2\mathrm{\Delta }_0k_xk_y/k_F^2`$ (the $`x`$ axis makes an angle $`\pi /4`$ with the $`a`$ axis of the $`CuO_2`$ planes). Let us orient $`𝐇`$ along the $`c`$ axis ($`H_{c1}HH_{c2}`$) and consider two types of vortex lattices: (I) rectangular lattice with primitive translations $`𝐚_1=a𝐱_0`$, $`𝐚_2=\sigma a𝐲_0`$; (II) centered rectangular lattice with $`𝐚_1=a𝐱_0`$, $`𝐚_2=a(𝐱_0/2\sigma 𝐲_0)`$, where $`H\sigma a^2=\varphi _0`$ is the flux quantum, and $`𝐱_0`$, $`𝐲_0`$, $`𝐳_0`$ are the unit vectors of the coordinate system.

Van Hove singularities. Our consideration is based on the analysis of the Bogolubov-de Gennes (BdG) equations for low energy excitations with momenta close to a certain gap node direction (e.g., $`𝐤_1=k_F𝐱_0`$): $`(\widehat{H}_0+\widehat{H}^{})\widehat{g}=\epsilon \widehat{g}`$, where $`\widehat{g}=(u,v)`$ is the QP wavefunction, $`\widehat{H}_0=V_F\widehat{\sigma }_z\widehat{p}_x+V_\mathrm{\Delta }\widehat{\sigma }_x\widehat{p}_y`$, $`\widehat{\sigma }_x,\widehat{\sigma }_z`$ are the Pauli matrices, $`\widehat{H}^{}=MV_FV_{sx}(1+\widehat{\sigma }_z)+MV_\mathrm{\Delta }V_{sy}\widehat{\sigma }_x`$, $`M`$ is the electron effective mass, $`\widehat{𝐩}=i\mathrm{}e𝐀/c`$, $`V_F=\mathrm{}k_F/M`$, $`V_\mathrm{\Delta }=2\mathrm{\Delta }_0/(\mathrm{}k_F)`$, $`𝐇=H𝐳_0`$, $`𝐀=Hy𝐱_0`$, $`𝐕_s=(V_{sx},V_{sy})`$. The spectrum of the Dirac Hamiltonian $`\widehat{H}_0`$ can be obtained using the usual quantization rule for a cyclotron orbit (CO) area . The periodic potential $`\widehat{H}^{}`$ removes the degeneracy of the discrete energy levels with respect to the CO center and induces a band structure in the spectrum . The general solution can be written in the form of a magnetic Bloch wave:

$$\widehat{g}=\underset{n}{}e^{ix(q_x+2\pi n/a)+2in\sigma q_ya}\widehat{G}(y2n\sigma a,𝐪),$$
(1)   
where $`n`$ is an integer, and $`𝐪`$ is the quasimomentum lying within the first magnetic Brillouin zone (MBZ): $`\pi /(2a)<q_x<\pi /(2a)`$, $`\pi /(2\sigma a)<q_y<\pi /(2\sigma a)`$. The wavefunction $`\widehat{G}(y,𝐪)`$ is localized in the domain with the size $`L`$ determined by $`𝐪`$ and energy values. The potential $`\widehat{H}^{}`$ results in the splitting of the CO near MBZ boundaries (see Fig. 1) and the spectrum consists of branches which correspond to the splitted portions of the CO.

For large Dirac cone anisotropy $`\alpha =V_F/V_\mathrm{\Delta }1`$ ($`\alpha =k_F\xi _0/2`$) and $`\epsilon <0.5\epsilon ^{}`$ ($`\epsilon ^{}=\pi \mathrm{}V_F/a\mathrm{\Delta }_0\sqrt{H/H_{c2}}`$) the harmonics in Eq. (1) do not overlap ($`L<2\sigma a`$) and one can replace $`\widehat{H}^{}`$ by the potential $`\widehat{H}^{}_x`$ averaged in the $`x`$ direction (see ). Such a simplification is a natural consequence of a small size of the cyclotron orbit (CO1 in Fig. 1) in the nodal direction as compared to the size of the MBZ. The energy spectrum consists of branches $`\epsilon _n(q_x=\pi Q/a)=\epsilon ^{}E_n(Q,\sigma \alpha )`$, which are displayed in Fig. 2 in the first MBZ for $`\pi \sigma \alpha =50`$ and $`\pi \sigma \alpha =100`$.

The number of energy branches which cross the Fermi level can be determined as follows: $`N2q_{}^{}/\delta q_{}2\sqrt{\pi \sigma \alpha }`$, where $`q_{}^{}`$ is the minimum possible size of the CO in the $`q_{}`$ direction and $`\delta q_{}`$ is the distance between MBZ boundaries. Each energy branch has an extremum as a function of the momentum $`q_x`$ near the MBZ boundary at a certain $`\stackrel{~}{\epsilon }_n`$ (we neglect here additional extrema which appear due to the exponentially small splitting of energy levels near the points of intersection of the branches in the $`EQ`$ plane). Due to the one-dimensional (1D) nature of the low energy spectrum the divergent contributions to the DOS take the form: $`\delta N(\epsilon )|\epsilon \stackrel{~}{\epsilon }_n|^{1/2}`$ ($`\epsilon >\stackrel{~}{\epsilon }_n`$ for energy minima and $`\epsilon <\stackrel{~}{\epsilon }_n`$ for maxima). The distance between these peaks $`\delta \epsilon \epsilon ^{}/(2\sigma \alpha )`$ coincides with a characteristic energy scale corresponding to van Hove singularities which occur when the CO intersects MBZ boundaries in the $`q_y`$ direction (see Fig. 1). The crossover between 1D and 2D regimes in the band spectrum occurs at $`\epsilon _c0.5\epsilon ^{}`$, when the CO size in the $`q_{}`$ direction becomes larger than the size of the first MBZ (CO2 in Fig. 1). For $`\epsilon \stackrel{_>}{_{}}\epsilon _c`$ the $`q_y`$-dependence of energy becomes essential and results in the appearance of 2D critical points, i.e. 2D local maxima (or minima) and saddle-points. Thus, instead of square-root van Hove singularities we obtain a set of discontinuities and logarithmic peculiarities ($`\delta N(\epsilon )ln|1\epsilon /\stackrel{~}{\epsilon }_n|`$) in the DOS, respectively. Obviously, these 2D singularities are more sensitive to temperature and finite lifetime effects and, consequently, the suppression of the corresponding oscillatory structure in the DOS should be stronger in the high energy regime. The above analysis can be generalized for gap nodes at $`𝐤=\pm k_F𝐲_0`$: the corresponding energy scales take the form $`\delta \epsilon 0.5\epsilon ^{}/\alpha `$, $`\epsilon _c0.5\epsilon ^{}/\sigma `$.

Even in the low energy regime the DOS oscillations with the energy scale $`\delta \epsilon `$ are surely smeared due to a finite scattering rate $`\mathrm{\Gamma }`$ and temperature and can be observed only for a moderate Dirac cone anisotropy and rather large magnetic fields. Comparing our results with a numerical solution  of the BdG equations for $`\sigma =1`$ and $`\alpha =5/2`$ we find that the above mechanism gives a good estimate of the energy scale of the double-peak structure in the tunneling conductance at the core center at $`H/H_{c2}=0.3`$ ($`\delta \epsilon 0.1\mathrm{\Delta }_0`$) and can explain the absence of this structure at low fields $`H/H_{c2}=0.05`$ due to temperature broadening ($`T=0.1T_c>\delta \epsilon 0.05\mathrm{\Delta }_0`$). In principle, van Hove singularities may account for peaks with a large energy gap $`\mathrm{\Delta }_0/4`$ observed experimentally at the vortex centers in YBaCuO  at $`H6T`$ provided we assume $`\alpha 1`$. Unfortunately the latter assumption is not consistent with the results of thermal conductivity measurements  ($`\alpha 14`$), and, thus, the nature of the experimentally observed peaks is still unclear. It is also necessary to stress here that the critical points in the DOS are a direct consequence of perfect periodicity and the introduction of rather strong disorder surely remove these singularities.

Zero-bias conductance. Hereafter we neglect the DOS oscillations, discussed above, and consider the peculiarities of the zero-bias tunneling conductance $`g(𝐫)`$ starting from a modified semiclassical model proposed in . According to this approach, the Doppler shift of the QP energy, which plays important role for $`\epsilon \stackrel{_<}{_{}}\mathrm{\Delta }_0\sqrt{H/H_{c2}}`$, appears to be averaged in the nodal direction due to an extremely large size of a semiclassical wave packet in this energy interval. Within such an approximation a diagonal (retarded) Green’s function can be written in the form:

$$G^R(𝐤,\epsilon ,𝐫)=\frac{\epsilon +i\mathrm{\Gamma }\mathrm{}𝐤_F𝐕_{av}+ϵ_𝐤}{(\epsilon +i\mathrm{\Gamma }\mathrm{}𝐤_F𝐕_{av})^2\mathrm{\Delta }_𝐤^2ϵ_𝐤^2},$$
(2)   
where $`ϵ_𝐤`$ is the normal state electron dispersion, $`𝐕_{av}=𝐕_s_x+𝐕_s_y`$. The scattering rate $`\mathrm{\Gamma }`$ should be determined self-consistently: $`\mathrm{\Gamma }=N(\mathrm{\Gamma },\epsilon )/(2N_F\tau )`$ (Born limit), $`\mathrm{\Gamma }=N_F\mathrm{\Gamma }_u/N(\mathrm{\Gamma },\epsilon )`$ (unitary limit), where $`2\tau `$ and $`N_F`$ are the relaxation time and DOS at the Fermi level in the normal state, $`\mathrm{\Gamma }_u=n_{imp}/(\pi N_F)`$, $`n_{imp}`$ is the concentration of of point potential scatterers, and $`N(\mathrm{\Gamma },\epsilon )=ImG^Rd^2k/(2\pi ^3)`$ is the local DOS.

Let us first consider the effect of finite temperature on the zero-bias conductance in the clean limit ($`\mathrm{\Gamma }0`$). The expression for the normalized conductance reads:

$$\stackrel{~}{g}(𝐫)=\frac{g(𝐫)}{g_N}=\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\frac{N(\mathrm{\Gamma }=0,\epsilon )d\epsilon }{4N_FTcosh^2\left(\frac{\epsilon }{2T}\right)}=$$

$$\frac{T\mathrm{𝑙𝑛}2}{\mathrm{\Delta }_0}+\frac{T}{2\mathrm{\Delta }_0}\underset{i=x,y}{}\mathrm{𝑙𝑛}cosh\frac{\epsilon _i^{}\mathrm{\Phi }_i}{4T}$$
(3)   
Here $`g_N`$ is the normal state conductance, $`\mathrm{\Phi }_x=\mathrm{\Phi }(x/R_x)`$, $`\mathrm{\Phi }_y=\mathrm{\Phi }(y/R_y)`$, $`\mathrm{\Phi }(z)=2z(2m+1)`$ for $`m<z<m+1`$ ($`m`$ is an integer), $`R_y`$ ($`R_x`$) is the distance between the lines parallel to the $`x`$ ($`y`$) axis and passing through the vortex centers, $`\epsilon _x^{}=\pi \mathrm{}V_FHR_x/\varphi _0`$, $`\epsilon _y^{}=\pi \mathrm{}V_F\sigma /R_y`$, For type I (II) lattices we have $`R_x=a`$, $`R_y=\sigma a`$ ($`R_x=a/2`$, $`R_y=\sigma a`$). One can separate two qualitatively different regimes in the behavior of the conductance:   
(i) superflow dominated regime $`T\epsilon _{x,y}^{}`$,

$$\stackrel{~}{g}\frac{1}{8}\sqrt{\frac{\pi \sigma H}{2H_{c2}}}F_1(x,y),$$
(4)   
(ii) temperature dominated regime $`T\epsilon _{x,y}^{}`$,

$$\stackrel{~}{g}\frac{T\mathrm{𝑙𝑛}2}{\mathrm{\Delta }_0}+\frac{\pi \mathrm{\Delta }_0\sigma H}{32TH_{c2}}F_2(x,y),$$
(5)   
where $`F_1(x,y)=|\mathrm{\Phi }_x|(R_x/R_y)+|\mathrm{\Phi }_y|`$ and $`F_2(x,y)=\mathrm{\Phi }_x^2(R_x/R_y)^2+\mathrm{\Phi }_y^2`$. In Fig. 3 we display the contour plots of the functions $`F_1(x,y)`$, $`F_2(x,y)`$ for a square lattice of type I (which is close to the one observed experimentally in YBaCuO ). There are two consequences of an increase in temperature: (i) first, the spatial dimensions of peaks in the local DOS become rather small comparing to the intervortex distance only for $`T>T^{}\mathrm{\Delta }_0\sqrt{H/H_{c2}}`$; (ii) second, the amplitude of the peaks appears to be essentially suppressed in the limit $`TT^{}`$. For magnetic fields $`H6T`$ (which is typically the field of STM experiment ) one obtains $`T^{}20K`$. Thus, we conclude that the finite temperature effects can not explain neither the narrow zero-bias conductance peaks observed in YBaCuO nor the absence of these peaks in BiSrCaCuO at $`T=4.2K`$. To explain these experimental facts it is necessary to take account of the finite lifetime effects which can stronly influence on the behavior of the DOS, as it follows from the results of Refs.  obtained on the basis of the usual semiclassical approach with a local Doppler shift. Starting from the modified semiclassical model (2) we obtain the following expression for the tunneling conductance at $`T=0`$:

$$\stackrel{~}{g}=\frac{N(\mathrm{\Gamma },\epsilon =0)}{N_F}=\frac{\mathrm{\Gamma }}{4\pi \mathrm{\Delta }_0}\left(4\mathrm{𝑙𝑛}\frac{\mathrm{\Delta }_0}{\mathrm{\Gamma }}+\underset{i=x,y}{}f_i\right)$$
(6)   
$$f_i=\frac{\epsilon _i^{}|\mathrm{\Phi }_i|}{\mathrm{\Gamma }}tan^1\frac{\epsilon _i^{}|\mathrm{\Phi }_i|}{2\mathrm{\Gamma }}\mathrm{𝑙𝑛}\left(1+\frac{(\epsilon _i^{}\mathrm{\Phi }_i)^2}{4\mathrm{\Gamma }^2}\right)$$

Obviously Born scatterers result only in a moderate change of the DOS (see ) since the corresponding $`\mathrm{\Gamma }`$ value for $`\mathrm{\Delta }_0\tau 1`$ is very small comparing to $`\epsilon _{x,y}^{}`$ and the conductance is given by Eq. (4). On the contrary, in the unitary limit the expression (4) is valid only in the clean case $`\mathrm{\Gamma }_u\mathrm{\Gamma }_{x,y}^{}0.1\epsilon _{x,y}^2/\mathrm{\Delta }_0`$ (for a square lattice $`\mathrm{\Gamma }_{x,y}^{}0.1\mathrm{\Delta }_0H/H_{c2}`$). In the dirty limit $`\mathrm{\Gamma }_u\mathrm{\Gamma }_{x,y}^{}`$ we obtain:

$$\stackrel{~}{g}\stackrel{~}{g}(H=0)\left(1+\frac{\mathrm{\Delta }_0H\sigma }{64\mathrm{\Gamma }_uH_{c2}}F_2(x,y)\right),$$
(7)   
where $`\stackrel{~}{g}(H=0)0.5\sqrt{\mathrm{\Gamma }_u/\mathrm{\Delta }_0}`$. In the vicinity of each vortex center the local DOS exhibits a fourfold symmetry with maxima along the nodal directions in a good agreement with numerical calculations based on the Eilenberger theory . For $`H=6T`$ finite lifetime effects become substantial if we assume $`\mathrm{\Gamma }_u\stackrel{_>}{_{}}10^2\mathrm{\Delta }_0`$. Thus, our approach allows to explain rather narrow conductance peaks (see Fig. 3b) observed near vortex centers in YBaCuO , even without taking account of the nontrivial structure of the tunneling matrix element, discussed in . With a further increase of the $`\mathrm{\Gamma }_u`$ value the amplitude of the peaks at the vortex centers vanishes: $`\delta \stackrel{~}{g}\sqrt{\mathrm{\Delta }_0/\mathrm{\Gamma }_u}(H/H_{c2})`$. Such a high sensitivity of the $`\delta \stackrel{~}{g}`$ value to finite lifetime effects can probably explain the difficulties in the observation of these peaks in the mixed state of BiSrCaCuO . Note in conclusion that according to Eq. (7) the spatially averaged DOS in the dirty limit varies as $`H`$ rather than $`HlnH`$ (the latter dependence has been predicted in  within the semiclassical approach taking account of the local $`𝐕_s`$ value).

I am pleased to acknowledge useful discussions with Dr.N.B.Kopnin, Dr.Yu.S.Barash, Dr.A.A.Andronov, Dr.I.D.Tokman, and Dr.D.A.Ryndyk.