Electronic structure of mesoscopic superconducting disk: Quasiparticle tunneling between the giant vortex core and disk edge

TL;DR

This paper investigates how quasiparticle tunneling influences the electronic structure of giant vortex states in mesoscopic superconducting disks, revealing a temperature-dependent crossover affecting tunneling conductance and potential cooling applications.

Contribution

It introduces a detailed analysis of quasiparticle tunneling effects in mesoscopic superconductors using the Usadel approach, highlighting a temperature-driven crossover in tunneling regimes.

Findings

Quasiparticle tunneling significantly alters local density of states.

A crossover from edge-dominated to core-dominated tunneling occurs with decreasing temperature.

Magnetic field influences quasiparticle cooling efficiency in mesoscopic devices.

Abstract

The electronic structure of the giant vortex states in a mesoscopic superconducting disk is studied in a dirty limit using the Usadel approach. The local density of states profiles are shown to be strongly affected by the effect of quasiparticle (QP) tunneling between the states localized in the vortex core and the ones bound to the sample edge. Decreasing temperature leads to a crossover between the edge-dominated and core-dominated regimes in the magnetic field dependence of the tunneling conductance. This crossover is discussed in the context of the efficiency of quasiparticle cooling by the magnetic field induced QP traps in various mesoscopic superconducting devices.

Figures

Equations

$$\Delta (\mathbf{r})={\Delta }_L(r){\mathrm{e}}^{iL\phi },$$

$$\mathcal{F}(\mathbf{r},{\omega }_n)={\mathcal{F}}_L(r,{\omega }_n){\mathrm{e}}^{iL\phi }.$$

$$\begin{gathered} -\frac{\hslash D}{2}\left[\frac{1}{r}\frac{d}{dr}\left(r\frac{d{\theta }_L}{dr}\right)-{\left(\frac{L-{\varphi }_r}{r}\right)}^2\sin {\theta }_L\cos {\theta }_L\right] \\ +{\omega }_n\sin {\theta }_L={\Delta }_L(r)\cos {\theta }_L. \end{gathered}$$

$$\frac{{\Delta }_L(r)}{g}-2\pi T\sum _{n\ge 0}\sin {\theta }_L=0.$$

$$\frac{1}{g}=\sum _{n=0}^{{\Omega }_{\mathrm{D}}/(2\pi T_{\mathrm{cs}})}\frac{1}{n+1/2}\simeq \ln [{\Omega }_{\mathrm{D}}/2\pi T_{\mathrm{cs}}]+2\ln 2+\gamma ,$$

$${\frac{d{\Delta }_L}{dr}|}_R={\frac{d{\theta }_L}{dr}|}_R=0.$$

$$\begin{gathered} -\frac{\hslash D}{2}\left[\frac{1}{r}\frac{d}{dr}\left(r\frac{d{\theta }_L}{dr}\right)-{\left(\frac{L-{\varphi }_r}{r}\right)}^2{\theta }_L\right] \\ +{\omega }_n{\theta }_L={\Delta }_L(r), \end{gathered}$$

$$\frac{{\Delta }_L(r)}{g}-2\pi T\sum _{n\ge 0}{\theta }_L=0.$$

$${\theta }_L(r,{\omega }_n)=\frac{{\Delta }_L(r)}{{\omega }_n+{\Omega }_L},$$

$${\Delta }_L(r)$$

$$f_L({\varphi }_r)$$

$$a_L=\frac{1}{2}\left(|L|-L+1-\frac{{\Phi }_0{\Omega }_L}{\pi \hslash DH}\right),b_L=|L|+1.$$

$$\begin{gathered} {\Gamma }_L(a_L,\varphi )=b_L\left(|L|-\varphi \right)K(a_L,b_L,\varphi ) \\ +2\varphi a_{L}K(a_L+1,b_L+1,\varphi )=0. \end{gathered}$$

$$\ln \frac{T_L}{T_{\mathrm{cs}}}=\Psi \left(\frac{1}{2}\right)-\Psi \left(\frac{1}{2}+\frac{{\Omega }_L}{2\pi T_L}\right),$$

$${\Omega }_{\mathrm{c}}=\underset{L,n}{\min }\left\{{\Omega }_L^{(n)}(\varphi ,R)\right\}$$

$${\Gamma }_L(a_L,{\varphi }_L)=0,{\Gamma }_{L+1}(a_{L+1},{\varphi }_L)=0,$$

$$H_{c2}=\underset{L}{\max }\left\{H_L(T)\right\},$$

$$F_{L}=2\pi T\sum\limits_{\omega_{n}<\Omega_{\rm D}}\int\limits_{0}^{R}r\,dr\bigg{\{}\left(\frac{\partial\theta_{L}}{\partial r}\right)^{2}+\left(\frac{L-\phi_{r}}{r}\right)^{2}\,\sin^{2}\theta_{L}\bigg{.}\\ \bigg{.}-\frac{4}{\hbar D}\left(\omega_{n}\cos\theta_{L}+\Delta_{L}\sin\theta_{L}\right)\bigg{\}}+\frac{2}{\hbar D}\int\limits_{0}^{R}r\,dr\,\frac{\Delta_{L}^{2}}{g}\,.$$

$$G_L(\epsilon ,r,\varphi )$$

$$\frac{\hbar D}{2}\left(F_{L}-F_{\rm N}\right)=\int\limits_{0}^{R}r\,dr\,\bigg{\{}\frac{\Delta_{L}^{2}}{g}-2\pi T\sum\limits_{\omega_{n}<\Omega_{\rm D}}\Delta_{L}\theta_{L}\\ -2\pi T\sum\limits_{\omega_{n}<\Omega_{\rm D}}\bigg{[}\frac{\hbar D}{2}\left(\frac{L-\phi_{r}}{r}\right)^{2}\,\frac{\theta_{L}^{4}}{3}+\omega_{n}\frac{\theta_{L}^{4}}{12}-\frac{\Delta_{L}\theta_{L}^{3}}{3}\bigg{]}\bigg{\}}\,,$$

$$\frac{1}{g}=\sum _{{\omega }_{n\mathrm{c}}<{\Omega }_{\mathrm{D}}}\frac{2\pi T_{\mathrm{c}}(H)}{{\omega }_{n\mathrm{c}}+{\Omega }_L}$$

$$\frac{\hbar D}{4\pi}\left(F_{L}-F_{\rm N}\right)\equiv\frac{\hbar D}{4\pi}\left(-A\,C_{L}^{2}+B\,C_{L}^{4}\right)=\\ \int\limits_{0}^{R}rdr\bigg{\{}\Delta_{L}^{2}\left(\sum\limits_{\omega_{n{\rm c}}<\Omega_{\rm D}}\frac{T_{\rm c}(H)}{\omega_{n{\rm c}}+\Omega_{L}}-\sum\limits_{\omega_{n}<\Omega_{\rm D}}\frac{T}{\omega_{n}+\Omega_{L}}\right)\\ +\Delta_{L}^{4}T\sum\limits_{\omega_{n}<\Omega_{\rm D}}\left[\frac{1}{4(\omega_{n}+\Omega_{L})^{3}}+\frac{\Omega_{L}-2\hbar D\left(\frac{L-\phi_{r}}{r}\right)^{2}}{12(\omega_{n}+\Omega_{L})^{4}}\right]\bigg{\}}\,.$$

$$A=\frac{\Phi_{0}}{\pi\hbar DH}I_{2,0}\times\\ \Bigl{[}\Psi(\omega_{L,T_{\rm c}})-\Psi(\omega_{{\rm D},T_{\rm c}})-\Psi(\omega_{L,T})+\Psi(\omega_{{\rm D},T})\Bigr{]}$$

$$B=\frac{1}{6(2\pi T)^{3}}\Bigl{\{}\frac{\Phi_{0}I_{4,0}}{2\pi\hbar DH}\left[6\pi T\zeta_{3}(\omega_{L,T})+\Omega_{L}\zeta_{4}(\omega_{L,T})\right]\\ -\left(I_{4,1}-2LI_{4,0}+L^{2}I_{4,-1}\right)\zeta_{4}(\omega_{L,T})\Bigr{\}}\ ,$$

$$I_{n,k}={\int }_0^{\varphi }{f}_L^n({\varphi }_r){\varphi }_r^{k}d{\varphi }_r.$$

$$N_L(\epsilon ,r,\varphi )$$

$$N_{L+1}-N_{L}\bigg{|}_{\begin{subarray}{l}\epsilon=0\\ r=R\\ \phi=\phi_{L}\end{subarray}}\sim\frac{\Delta_{L}^{2}(R)-\Delta_{L+1}^{2}(R)}{2\Omega_{L}^{2}}\ .$$

$$\frac{d}{d{\varphi }_r}\left({\varphi }_r\frac{d{\theta }_L}{d{\varphi }_r}\right)-\frac{(L-{\varphi }_r{)}^2}{4{\varphi }_r}{\theta }_L+\frac{{\Phi }_0{\Omega }_L}{2\pi \hslash DH}{\theta }_L=0$$

$$\displaystyle\theta_{L}(r)=e^{-\phi_{r}/2}\phi_{r}^{|L|/2}W(\phi_{r})\$$

$${\varphi }_r\frac{d^{2}W}{d{\varphi }_r^2}+(b_L-{\varphi }_r)\frac{dW}{d{\varphi }_r}-a_{L}W=0,$$

$$a_L=\frac{1}{2}\left(|L|-L+1-\frac{{\Phi }_0{\Omega }_L}{\pi \hslash DH}\right),b_L=|L|+1.$$

Full text

I Introduction

Vortex states in mesoscopic superconducting (SC) systems of the size comparable to the superconducting coherence length, have been well studied over the past few decades, mainly with the emphasis on the dependence of the vortex configuration on the size and geometry of the sample. In such nanoscale samples theory predicts that only a few vortices can be placed, and confinement effects result in different exotic vortex configurations unlike the triangular Abrikosov lattice fink ; Moshchalkov-Nature95 ; Peeters-Nature97 ; geim ; Buisson ; Peeters-PRL98 ; palacios ; bruy ; jadallah ; Moshchalkov-Nature00 ; Melnikov-Nature02 ; Melnikov-PRB02 ; Peeters-PRL06 . These exotic configurations are formed by the interplay between imposed boundary conditions and the repulsive interactions between vortices.

The most remarkable consequence of this interplay is the formation of the so-called giant vortex state or multiquantum vortex when all the vortices merge in the disk center predicted mostly within the Ginzburg-Landau formalism provided the disk size is of order of the coherence length. A variety of experimental methods have been used to verify these theoretical predictions: (i) Hall probe microscopy Peeters-Nature97 ; geim ; Moshchalkov-PRB04 ; Moshchalkov-PRB08 , (ii) Bitter decoration Grigorieva , (iii) scanning SQUID microscopy Kokubo-PRB10 , (iv) different tunneling experiments including scanning tunneling microscopy/spectroscopy studies Peeters-PRL04 ; Roditchev-PRL11 ; Roditchev-PRL09 ; Nishio-PRL08 ; Moshchalkov-PRB16 .

The latter experimental approach is known to be sensitive to the electronic structure of the sample, namely, to the local density of states of quasiparticle excitations and, thus, the phenomenological Ginzburg – Landau theory often appears to be insufficient for the interpretation of the experimental data. This clear demand to the microscopic theory has stimulated theoretical activity in the field concentrated mainly on the calculations based on the Bogolubov-de-Gennes theory Melnikov-Nature02 ; Melnikov-PRB02 ; tanaka1 ; tanaka2 ; rainer ; virt ; esch ; tanaka3 ; Melnikov-PRL ; mrs1 ; mrs2 ; Peeters-PRL12 , i.e., on the clean limit corresponding to the very large mean free path $\ell$ well exceeding both the coherence length $\xi_{0}$ and the sample size. Certainly, the predictions made within such approach may be difficult to use for most of the experimentally available samples for which the dirty limit conditions ($\ell\ll\xi_{0}$) are much more appropriate. In particular, it is natural to expect that all the density of states features associated, e.g., with the different anomalous spectral branches volovik in the giant vortex or with the mesoscopic oscillations of the Caroli-de Gennes-Matricon energy levels  Melnikov-PRL due to the finite sample size should be smeared by disorder. An adequate theoretical description of the sample electronic structure in this diffusive regime should be, of course, based on the Usadel-type theory. And indeed such calculations are known to provide an excellent tool for the analysis of the Abrikosov vortex lattices in unrestricted geometries (see, e.g., Golubov-PRL94 ). For multiquantum giant vortices these results have been generalized in Ref. Silaev-JPCM13 without accounting the effect of the sample boundary.

It is important to note that the demand in the theoretical explanation of the available data of scanning tunneling microscopy and spectroscopy (STM/STS) on the exotic vortex structures in mesoscopic samples (see, e.g., NS_vort_PRL ; NS_vort_JETP ) is not the only motivation for the continuing research work in the field. Nowadays, superconducting nanostructures have become an important element in designing devices for rapidly expanding fields of quantum computing, quantum memory, superconducting logic, and metrology and they are obviously the main building blocks for the superconducting electronics. However, superconductors are known to be easily poisoned by non-equilibrium quasiparticles and these extra excitations drastically affect the performance of above-mentioned quantum devices, e.g. via overheating or unwanted population in general. To suppress overheating in a superconductor different types of quasiparticle traps are used (see, e.g., Refs. Nguyen-NJP13 ; Ullom-APL98 ; Khaymovich-NatCom16 and references therein). One of the possible types of quasiparticle traps can be formed by regions with the reduced superconducting gap which appear in the Meissner and vortex states and can be successfully controlled by the external magnetic field (see Peltonen2011 ; Nsanzineza2014 ; Wang2014 ; Vool2014 ; Woerkom2015 ; Khaymovich-NatCom16 ; Nakamura17 ). Further progress in the field requires a quantitative theoretical description of both types of quasiparticle traps based on the vortex penetration as well as on Meissner currents flowing mostly at the sample edge. Thus, the main goal of our work is to analyze the behavior of the local density of states in giant vortices penetrating to a circular superconducting sample of order of the coherence length, with proper accounting of the sample edge effects. This analysis, to our mind, should provide an important step on the route to rather general model of quasiparticle traps in mesoscopic samples.

To elucidate the key results of our study it is useful to note that both the giant vortex cores and the sample edge with the flowing Meissner screening currents can be clearly viewed as Andreev potential wells for quasiparticles in the clean limit Melnikov-Nature02 . On the other hand, the impurity scattering in the dirty limit surely modifies some spectral characteristics of these wells compared to the clean regime: (i) scattering broadens the discrete levels of the Caroli-de Gennes-Matricon energy branch CdGM , which crosses the Fermi level, suppressing the minigap in the spectrum Golubov-PRL94 , (ii) scattering can also result in the increase of minigap in the quasiparticle spectrum $E_{\rm g}$ at the sample edge because the changes in the quasiparticle momentum directions partially suppress the effect of the Doppler shift of the quasiparticle energy in the presence of the surface currents eg1 ; eg2 ; eg3 ; eg4 . The overall spectral characteristics and local density of states of the mesoscopic sample can be considered as an interplay of the subgap states, located in the vortices and in the regions with the reduced spectral gap $E_{\rm g}$ by Meissner currents, especially at the sample edge. To illustrate this interplay we consider for instance different contributions to the zero bias conductance (ZBC) at the sample edge (see Fig. 1), which can be experimentally accessed in tunneling transport measurements. The contribution of the giant vortex core states to this quantity can be estimated as follows: $\sim\exp(-R/d_{L})$, where $R$ is the distance from the vortex center to the boundary and $d_{L}$ is the effective decay length dependent on the vorticity $L$. For $L=1$ the latter length $d_{1}\simeq\xi_{0}$ is of order of the coherence length $\xi_{0}$. The contribution of the edge states should include the temperature activation exponent $\sim\exp(-E_{\rm g}/T)$ due to the finite spectral minigap $E_{\rm g}$. These two terms are comparable for a characteristic temperature $T^{*}(R)\simeq E_{\rm g}d_{L}/R$. Thus, we conclude that in a sample of certain size $R$ for the temperatures larger than $T^{*}(R)$ the core contribution is negligible at the sample edge and, consequently, the finite temperature masks the coupling of the Andreev wells in the vortex core and at the edge. In the opposite limit of small temperatures $T<T^{*}(R)$, quantum-mechanical tunneling of the subgap quasiparticles between the vortex and the edge traps becomes observable in the experimentally measurable quantities and dominating over thermally-activated processes. Here and further the Boltzmann’s constant is set to unity, $k_{\rm B}=1$.

The above estimates give us a simple criterion of the interplay of the core and edge state contributions which will be quantitatively confirmed by further calculations of the local density of states (LDOS) in a diffusive mesoscopic SC disk in a wide interval of magnetic fields, applied perpendicular to the sample plane. Note that these estimates can be of course applied not only for a vortex in a finite size sample but also for any experimental geometry with vortices positioned close to the superconductor edge (see, e.g., STM images in Ref. guill ).

The paper is organized as follows. In Sec. II we briefly discuss the basic equations. In Sec. III we calculate the superconducting critical temperature $T_{\rm c}$ and study the switching between the states with different vorticity $L$ while sweeping the magnetic field. In Sec. IV we find both analytically and numerically the spatially resolved LDOS and study the behavior of the jumps in ZBC which are attributed to the entrance of a vortex into the disk. We summarize our results in Sec. V.

II Model and basic equations

Hereafter we consider a thin superconducting disk of a finite radius $R$ of order of the coherence length at the temperature $T$, placed in external magnetic field $\mathbf{H}=H\mathbf{z}_{0}$ oriented perpendicular to the plane of the disk (Fig. 1). The disk thickness is assumed to be small compared to the London penetration depth, thus, the effective magnetic field penetration depth is large. This allows us to neglect the contributions to the magnetic field from supercurrents and, thus, $\mathrm{rot}\mathbf{A}=\mathbf{B}\equiv\mathbf{H}$. Using the notations $\tau^{-1}$ for the electron elastic scattering rate and $T_{\rm cs}$ for the bare superconductor transition temperature the dirty limit conditions can be written as $T_{\rm cs}\tau\ll 1$. In this regime the normal ($\mathcal{G}$) and anomalous ($\mathcal{F}$) quasiclassical Green’s functions are described by the Usadel equations Usadel-PRL70 , which are valid for the whole temperature and magnetic field range. Focusing on the axisymmetric multiquantum vortex states with the vortex core positioned in the center of the disk $r=0$

[TABLE]

we consider solutions, homogeneous along the $z$-axis and characterized by a certain angular momentum $L$, referred further as vorticity

[TABLE]

Here we choose the cylindrical coordinate system $(r,\varphi,z)$ and the gauge $\mathbf{A}=(0,\,A_{\varphi},\,0)$, $A_{\varphi}=rH/2$. Due to the symmetry of Usadel equations $\mathcal{F}$ is an even function of $\omega_{n}$, $\mathcal{F}(r,-\omega_{n})=\mathcal{F}(r,\omega_{n})$, so that it is enough to treat only positive $\omega_{n}$ values. In the standard trigonometrical parametrization $\mathcal{G}=\cos\theta_{L}$, $\mathcal{F}=\sin\theta_{L}\,\mathrm{e}^{iL\varphi}$, $\mathcal{F}^{\dagger}=\sin\theta_{L}\,\mathrm{e}^{-iL\varphi}$ the Usadel equations take the form

[TABLE]

The self-consistency equation for the singlet superconducting order parameter function reads

[TABLE]

Here $D=v_{\rm F}l/3$ is the diffusion coefficient, $\Phi_{0}=\pi\hbar c/e$ is the flux quantum, $\omega_{n}=\pi T(2n+1)$ is the Matsubara frequency at the temperature $T$, $\phi_{r}=\pi r^{2}H/\Phi_{0}$ is a dimensionless flux of the external magnetic field $\mathbf{H}$ threading the circle of certain radius $r$, and the pairing parameter $g$ determines the bare critical temperature $T_{\rm cs}$ as

[TABLE]

with the Debye frequency $\Omega_{\rm D}$ and the Euler – Mascheroni constant $\gamma\simeq 0.5772$. The coherence length $\xi_{0}=\sqrt{\hbar D/2\Delta_{0}}$ plays the role of a typical lengthscale in the Usadel equations.

The equations (3, 4) should be supplemented with the boundary conditions at the disk edge $r=R$:

[TABLE]

III Critical temperature of superconducting transitions with different vorticities

For the temperatures close to the critical temperature of the superconducting transition $T\lesssim T_{\rm c}(H)$, we can restrict ourselves by the solution of the Usadel equations Eqs. (3, 4) linearized in the anomalous Green function ($\sin\theta_{L}\simeq\theta_{L}$):

[TABLE]

[TABLE]

In these linearized equations the relation between the anomalous Green function $\theta_{L}(r)$ and the order parameter $\Delta_{L}(r)$ can be written in the standard form

[TABLE]

where $\Omega_{L}$ is the depairing parameter depending on the disk radius $R$ and the external magnetic field ${\bf H}$. Thus, the solution of Eq. (7) in the region $r\leq R$ can be expressed via the confluent hypergeometric function of the first kind (Kummer’s function $K(a,b,z)$ Abramowitz-Handbook )

[TABLE]

Here $C_{L}$ is a constant, and the parameters $a_{L}$ and $b_{L}$ depend on the vorticity $L$ as follows (see Appendix A for details):

[TABLE]

The boundary condition (6) for the orbital mode $L$ written in the form (10) results in the following algebraic equation

[TABLE]

The equation (11) determines the implicit dependence of the parameter $a_{L}$ on the flux $\phi(H)=\pi R^{2}H/\Phi_{0}\equiv H/H_{0}$ through the disk for a fixed value of vorticity $L$ normalized to the flux quantum. Here we introduced a characteristic field $H_{0}=\Phi_{0}/\pi R^{2}$.

The solutions $a_{L}^{(n)}$ of the equation (11) give a set of values $\Omega_{L}^{(n)}$ which depend on the normalized flux threading the whole disk $\phi$ and the disk radius $R$: $\Omega_{L}=\Omega_{L}(\phi,\,R)$. Finally, substituting the expression (9) into the self-consistency condition (8) one obtains the following equation for the critical temperature $T_{L}$ of the state with the vorticity $L$:

[TABLE]

where $\Psi$ is the digamma function.

In accordance with the self-consistency equation (12), the minimal value of the depairing parameter

[TABLE]

determines the vorticity $L_{\rm c}$ and the critical temperature $T_{\rm c}=T_{L_{\rm c}}$ of the orbital mode, which nucleates in the disk of the radius $R$ placed in the external magnetic $\mathbf{H}$.

Figure 2 shows typical dependencies of the critical temperature $T_{\rm c}$ and the depairing parameter $\Omega_{\rm c}$ on the external magnetic flux $\phi$ across the disk for a fixed value of the disk radius $R$. The phase boundary $T_{\rm c}(\phi)$ exhibits an oscillatory behavior similar to the well-known Little-Parks oscillations Little-PRL62 ; Parks-PRA64 , caused by the transitions between the states with different angular momenta $L$. The values of the normalized flux through the disk $\phi_{L}$, where the switching of the orbital modes $L\rightleftarrows L+1$ takes place, obey the equations

[TABLE]

and do not depend on the disk radius $R$: $\phi_{L}\simeq 1.92;\,3.40;\,4.74;\,6.04;\,7.30;\ldots$ for $L=0\div 5\ldots$. The magnetic field of the switching between modes $L$ and $L+1$ is determined by the expression $H_{\rm s}=H_{0}\phi_{L}$. The values of the dimensionless fluxes corresponding to the vorticity switching coincide with the ones found in Ref. jadallah for a superconducting disk within the Ginzburg–Landau theory. This coincidence comes from the obvious fact that the linearized Usadel equation (7) after the substitution of the expression (9) becomes similar to the linearized Ginzburg–Landau equation. Surely, this similarity does not extend to the full behavior of the $T_{L}(H)$ curve determined by Eq. (12). Note also that both the depairing factor and, thus, the critical temperature depend strongly on the disk radius $R$: $\Omega_{L}\sim R^{-2}\tilde{\Omega}_{L}(\phi)$, where $\tilde{\Omega}_{L}(\phi)$ is a certain function of the dimensionless flux $\phi$ only. One can see that the decrease in the $R$ value results in the decrease in the number of observable different vortex states.

IV Density of States

For a fixed temperature $T$ the orbital mode $L$ exists in the interval of the magnetic field values $0\leq H_{L1}\leq H\leq H_{L2}$ which satisfy the condition $T_{L}(\phi(H))\geq T$ (see the inset in Fig. 3). Figure 3 shows a typical temperature dependence of the upper critical field for the disk,

[TABLE]

affected by the transitions between different orbital states.

In order to analyze the characteristics of the sample far from the phase transition line we return back to the nonlinear Usadel theory and consider full free energy functional:

[TABLE]

Focusing now on the effect of the switching between different vortex states on the density of states we should note that experimentally this quantity can be most directly probed by the measurements of the local differential conductance:

[TABLE]

where $V$ is the applied bias voltage, $(dI/dV)_{N}$ is a conductance of the normal metal junction, and $f(\varepsilon)=1/\left(1+\exp(\varepsilon/T)\right)$ is the Fermi function.

IV.1 High magnetic field: $H\lesssim H_{\rm c2}$

We start our analysis from the limit of high magnetic fields close to the phase transition line $H_{\rm c2}(T)$ shown in Fig. 3 when the solution of the Usadel equations can be significantly simplified due to the smallness of the functions $\Delta_{L}$ and $\theta_{L}$. In this case one can use the solution of the linearized theory (9, 10). The constant $C_{L}$ in Eq. (10a) should be found from the nonlinear Usadel theory (3). For this purpose we write the corresponding free energy up to the fourth power of $\Delta_{L}$ and $\theta_{L}$:

[TABLE]

where $F_{\rm N}$ is the free energy of the normal state. Using the above self-consistency equation for $T_{\rm c}(H)$

[TABLE]

with $\omega_{n{\rm c}}=\pi T_{\rm c}(H)(2n+1)$, and the relation (9), we obtain

[TABLE]

Here the second (third) line corresponds to the quadratic (quartic) terms in $\Delta_{L}=C_{L}f_{L}(\phi_{r})$.

Finally the amplitude $C_{L}$ which minimizes the above functional $F_{L}=F_{\rm N}-A\,C_{L}^{2}+B\,C_{L}^{4}$ takes the form $C_{L}^{2}=A/(2B)$, with

[TABLE]

[TABLE]

where $\zeta_{k}(a)=\sum_{n\geq 0}1/(n+a)^{k}$ is the zeta function, $\omega_{L,T}=\Omega_{L}/(2\pi T)+1/2$, $\omega_{{\rm D},T}=(\Omega_{\rm D}+\Omega_{L})/(2\pi T)+3/2$, and

[TABLE]

Substituting now $\omega_{n}=-i\varepsilon$ in the relation (9), one obtains the following expressions for the LDOS valid in the first order in $\Delta_{L}^{2}(r)$ :

[TABLE]

The transitions between different vortex states while sweeping the magnetic field up, are visualized by abrupt changes (or jumps) in ZBC Nishio-PRL08 ; Roditchev-PRL09 ; Moshchalkov-PRB16 , which are determined by the LDOS $N(\varepsilon,r,\phi)$ at the Fermi level $\varepsilon=0$. Let us consider a certain point ($T=T_{\rm s}\,,H=H_{\rm s}$) at the phase diagram Fig. 3, where switching of the orbital modes $L\rightleftarrows L+1$ takes place, $H_{\rm s}=H_{0}\phi_{L}$. Since the depairing parameters of the orbital modes $L$ and $L+1$ coincide $\Omega_{L}=\Omega_{L+1}$ the corresponding jump in the LDOS $N_{L+1}-N_{L}$ at the disk edge $r=R$ can be estimated as follows:

[TABLE]

Figure 4 shows the magnetic field dependence of the normalized LDOS at the disk edge. The transitions between different vortex states ($L\to L+1$) are accompanied by the abrupt reduction in LDOS at the disk edge while sweeping the magnetic field up. Similar jumps of the LDOS, which are attributed to the entrance of a vortex inside the disk, have been observed in measurements of the normalized ZBC on $\mathrm{Pb}$ nano-island Roditchev-PRL09 and $\mathrm{MoGe}$ nanostructures Moshchalkov-PRB16 .

IV.2 An arbitrary magnetic field: $0<H<H_{\rm c2}$

As a next step, we analyze the conductance behavior as a function of magnetic field and temperature at arbitrary magnetic fields, $0<H<H_{\rm c2}$. The Usadel equations (3 - 6) have been solved numerically for different vorticities which allowed us to calculate and compare the values of the free energy.

Figure 5 shows the magnetic field dependence of the free energy $F$ (15) and the zero bias conductance $G_{L}(0\,,R\,,\phi)$ (IV) at the Fermi level for a small disk radius $R=4\xi_{0}$ and two temperatures $T=0.1T_{\rm cs}$ and $T=0.2T_{\rm cs}$. All three curves illustrate the switching between the states with different vorticities $L=0\div 4$, which are similar Little-Parks-like switching of the critical temperature $T_{\rm c}(H)$, Fig. 2. Sequential entries of vortices produce a set of branches $F_{L}$ with different vorticity $L$ on the $F(H)$ and $dI/dV(H)$ curves. The transitions between different vortex states are accompanied by an abrupt change in the ZBC, which is attributed to the entry/exit of a vortex inside the disk while sweeping the magnetic field. We observe the Meissner state when the total vorticity $L=0$ for $H\lesssim H_{{\rm s}0}=2.24H_{0}$, and a single-vortex state $L=1$ in the field range $H_{{\rm s}0}\lesssim H\lesssim H_{{\rm s}1}=3.84H_{0}$. In the Meissner state the ZBC is suppressed and spatially homogeneous: the ZBC value at the disk edge is slightly higher then ZBC value in the center. In the increasing magnetic field the gap in the tunneling spectra gradually fills with the quasiparticle states. This effect is more pronounced near the disk edge where the screening superconducting currents have higher density. The smooth evolution of ZBC continues till $H/H_{0}\simeq 2.24$ where it is interrupted by a vortex entry. At higher fields $H>H_{{\rm s}1}$ the multivortex states $L=2\div 4$ become energetically favorable.

Note, that the field values $H_{{\rm s}L}$ at which the jumps in vorticity ($L\to L+1$) occur are always larger than the values $H_{0}\,\phi_{L}$ found from the calculations of the critical temperature behavior.

Figure 5 illustrates an important point noted in Introduction, i.e., the temperature crossover between different regimes in the behavior of the conductance vs magnetic field. Indeed, one can clearly see that the change in temperature from $0.1T_{\rm cs}$ to $0.2T_{\rm cs}$ is accompanied by the change of the direction of jumps in the dependence of zero bias conductance vs magnetic field. The upward jumps in conductance for the lower temperature, $T<T^{*}(R)$, can be associated with the core dominated regime $e^{-E_{\rm g}/T}<e^{-R/d_{L}}$, see Fig. 1, when the conductance increases with the increase in the number of vortices trapped in the center of the sample and therefore in the parameter $d_{L}$. The downward jumps in conductance at higher temperatures, $T>T^{*}(R)$, are caused by the increase in the (soft) spectral gap value $E_{\rm g}$ at the sample edge as the vortex enters, which should result in the suppression of the subgap conductance $G_{L}\sim e^{-E_{\rm g}/T}$. The change in vorticity $L\to L+1$ in this case results in the decrease of the screening current density and the corresponding enhancement of superconductivity at the edge of the disk. Assuming the crossover temperature $T^{*}(R)$ to be in the interval $0.1T_{\rm cs}<T^{*}(R)<0.2T_{\rm cs}$ and taking $R=4\xi_{0}$ one can estimate the value $E_{\rm g}\sim 0.8\Delta_{0}$ which is in good agreement with the behavior of the energy dependence of the local density of states in Fig. 6. Indeed, the position of the maximal slope of the energy dependence of the density of states roughly gives the value of the minigap at the edge: $E_{\rm g}\lesssim 0.8\Delta_{0}$.

Figure 6 also illustrates the switching between the states with hard and soft gaps with the increase in the magnetic field. In the Meissner state ($H<H_{{\rm s}0}$) the hard minigap $\Delta_{\rm m}$ in the spectrum exists ($N(\varepsilon<\Delta_{\rm m},r,\phi)=0$) till the first vortex entry. The density of states in the center of the disk $N(\varepsilon,\mathrm{0},\phi)$ is equal to the electronic density of states at the Fermi level $N_{0}$ for any vortex state $L\geq 1$, indicating a full suppression of the spectral gap in the disk center due to the vortex entry. At the same time, at the edge of the disk the superconductivity survives though the gap becomes soft, $0<N(0,R,\phi)<N_{0}$.

Figure 7 presents the radial distributions of the SC order parameter $\Delta_{L}(r)$ and the ZBC $dI/dV$ at $T=0.1\,T_{\rm c}$ for different values of the magnetic field $H$ corresponding to the switching between the states with different vorticity $L$. The profiles of ZBC in multiquantum vortices $L>1$ reveal a plateau near the vortex center, which can be considered as a hallmark of the multiquantum vortex formation in dirty mesoscopic superconductors Roditchev-PRL11 ; Roditchev-PRL09 ; Silaev-JPCM13 .

The electronic properties of the vortex states look to be rather different if the radius of the disk $R$ is much larger than the coherence length $\xi_{0}$. In this case the core of a multiquantum vortex does not extend to the edge of the disk, and quasiparticles in the vortex core remain well localized near the disk center. Clearly, in this case the temperature crossover between the core-dominated and edge-dominated regimes accompanied by the change in the direction of the jumps in the local ZBC at the sample edge becomes much more difficult to observe due to the exponentially small values of the factors $\exp(-R/d_{L})$ and $\exp(-E_{\rm g}/T)$ near the crossover.

V Conclusions

To sum up, we have analyzed the behavior of the LDOS $N(\varepsilon\,,r\,,H)$ and conductance on an external magnetic field $H$ in a mesoscopic superconducting disk on the basis of Usadel equations. We have demonstrated that transitions between the superconducting states with different vorticities provoke abrupt changes (jumps) in the local zero bias conductance $dI/dV$ at the edge of the disk. These jumps of the ZBC are attributed to the entry/exit of vortices while sweeping the magnetic field. The transitions between different vortex states can be accompanied both by the decrease and increase in the ZBC while sweeping the magnetic field up. The direction of jumps in ZBC attributed to the vortex entry depends on the disk radius $R$ and the temperature $T$ and is determined by two opposite in sign contributions to conductance: (i) the entrance of a vortex into the disk is accompanied by the reduction of the supercurrents flowing along the sample edge and, thus, improves superconductivity at the edge; (ii) the entrance of a vortex increases the number of subgap quasiparticle states in the multiquantum vortex core which provide an additional contribution to the conductance because of the quasiparticle tunneling between the vortex core and the sample edge. To the best of our knowledge, the systematic experimental analysis of the direction of the ZBC jumps has not been done yet. However, these measurements can provide an additional information about the soft gap value governing not only the contribution to the tunneling transport, but also the one to the thermal relaxation mechanisms (see, e.g., Khaymovich-NatCom16 ; Nakamura17 ) and also about the classical-to-quantum interplay in quasiparticle tunneling in mesoscopic superconducting samples. These results are directly related to the quantitative characterization of the quasiparticle traps appearing in the Meissner and vortex states of superconductors (see, e.g., Peltonen2011 ; Nsanzineza2014 ; Wang2014 ; Vool2014 ; Woerkom2015 ; Khaymovich-NatCom16 ; Nakamura17 ), especially in the different types of single-electron sources based on hybrid superconducting junctions and working far from equilibrium Khaymovich-NatCom16 ; Nakamura17 ; parity_NISIN2015 ; Khaymovich_Basko2016 ; vanZanten2016 .

Acknowledgements.

This work was supported, in part, by the Russian Foundation for Basic Research under Grant No. 17-52-12044. In the part concerning the LDOS numerical calculations, the work was supported by Russian Science Foundation (Grant No. 17-12-01383). IMK acknowledges the support of the German Research Foundation (DFG) Grant No. KH 425/1-1.

Appendix A Solution of Eq. (7) via the Kummer’s functions

Substituting the expression of the order parameter $\Delta_{L}(r)=\theta_{L}(r)\left(\omega_{n}+\Omega_{L}\right)$, (9), into Eq. (7) and rewriting the latter in terms of the renormalized flux $\phi_{r}=\pi r^{2}H/\Phi_{0}$,

[TABLE]

one can easily obtain the equation for the function $W(\phi_{r})$ defined as

[TABLE]

[TABLE]

with the parameters $a_{L}$ and $b_{L}$ given by

[TABLE]

The solution of Eq. (27) in the region $r\leq R$ is confluent hypergeometric function of the first kind (Kummer’s function), $W=K(a_{L},\,b_{L},\,\phi_{r})$, which after substitution to the expression (26) for $\theta_{L}(r)$ gives the result (10) from the main text.