format: text

text:
# Theory of de Haas-van Alphen Effect in Type-II Superconductors

## I Introduction

The de Haas-van Alphen (dHvA) experiment on normal metals has been a unique and powerful tool to probe their Fermi surfaces.LK55 ; Shoenberg84  The main purpose of this paper is to establish theoretically that it can even be used to detect detailed gap structures of type-II superconductors.

Back in 1976, Graebner and Robbins discovered the dHvA oscillation in 2H-NbSe2 persisting down through the superconducting upper critical field $`H_{c2}`$.GR76  It was after 15 years later when Ōnuki et al. first reconfirmed it.Onuki92  Since then, however, a considerable number of materials have been found to display the dHvA oscillation in the vortex state. They include: $`A15`$ superconductors V3SiMueller92 ; Corcoran94  Nb3Sn,Harrison94  a borocarbide superconductor YNi2B2C,Heinecke95 ; Terashima95  heavy-fermion superconductors CeRu2,Hedo95  URu2Si2,Ohkuni97 ; Bergemann97  UPd2Al3,Inada99  CeCoIn5,Settai01  and an organic superconductor $`\kappa `$-(BEDT-TTF)2Cu(NCS)2;Wel94  see Refs. Janssen98,  and Inada99-2,  for a recent review. Basic features of the oscillation are be summarized as follows: (i) the dHvA frequencies remain unchanged through the transition; (ii) the oscillation amplitude experience an extra attenuation; (iii) the cyclotron mass does not change except for strongly correlated heavy fermion systems.

It is somewhat surprising that the dHvA oscillation is observable even in superconductors without a well-defined Fermi surface. Many theories have been presented to explain the persistent oscillation and the extra damping,Maki91 ; Miyake93 ; Gorkov98 ; Wasserman94 ; Stephen92 ; Curnoe00 ; Dukan95 ; Maniv92 ; Mineev00 ; Bruun97 ; Gvozdikov98 ; NMA95 ; Miller95  which may be classified into three categories.

The first approach applies a Bohr-Sommerfeld semiclassical quantization to either the Brandt-Pesch-TewordtBrandt67  (BPT) Green’s function near $`H_{c2}`$,Maki91  the electron number $`N`$ at $`H=0`$,Miyake93  or Dyson’s equation at $`H=0`$,Gorkov98  for obtaining the oscillatory behavior of the magnetization. As may be seen by this diversity of the applications, however, there is no unique semiclassical quantization scheme for quasiparticles superconductors; thus, the validity of the procedure is not clear. This category includes Maki’s theory,Maki91  which was later reproduced by Wasserman and SpringfordWasserman94  by treating the BTP self-energy as the extra broadening factor in the normal-state thermodynamic potential and then following Dingle’s procedure.Dingle52  However, the BTP self-energy itself is obtained within the quasiclassical approximation without the Landau-level structure so that this approximation may also be questionable.

The second approach relies on some approximate analytic solutions for the Bogoliubov-de Gennes (BdG) equations or the equivalent Gor’kov equations, such as the coherent-potential approximation,Stephen92 ; Curnoe00  the diagonal-pairing approximation,Dukan95  or the Ginzburg-Landau (GL) expansion for the free energy with respect to the order parameter.Maniv92 ; Mineev00 ; Bruun97  However, quantitative estimations of those approximations are yet to be carried out. It should also be noted that the GL expansion necessarily integrates out the quasiparticle degrees of freedom, so that the physical origin of the extra oscillation damping may be obscured in the GL approach.

The third approach solves the BdG equations numerically without approximations.NMA95 ; Miller95  Norman et alNMA95  thereby extracted an analytic formula for the dHvA oscillation damping through a fitting to their numerical data.NMA95  However, it is a two-dimensional calculation for the isotropic $`s`$-wave pairing where the number of Landau levels below the Fermi level is $`N_\mathrm{F}10`$ at $`H_{c2}`$. As may be realized from the appearance of the quantized Hall effect, two-dimensional systems in high magnetic fields may be qualitatively different from three-dimensional systems. Thus, the obtained formula may not be appropriate for describing real three-dimensional materials with $`N_\mathrm{F}1`$. On the other hand, another calculation by Miller and Györffy for a two-dimensional lattice model Miller95  corresponds to the low-field limit near $`H_{c1}`$comment1  and may not be suitable to explain the experiments. Moreover, calculations for lattice models have a flaw that they cannot yield continuous magnetic oscillation due to the commensurability condition between the underlying lattice and the vortex lattice.

Notice finally that most of the above theories consider only the isotropic $`s`$-wave pairing. Especially, no numerical studies have been performed yet for anisotropic pairings or three-dimensional systems.

With these observations, we will perform both numerical and analytic calculations for three-dimensional BdG equations with various gap structures. They can be solved efficiently by the Landau-level-expansion method (LLX), which was formulated for vortex-lattice states of arbitrary pairing symmetry Kita98-2  and used successfully to compare low-energy quasiparticle spectra between $`s`$\- and $`d`$-wave pairings.Yasui99  We will thereby clarify how the discrete Landau levels experience quantitative changes as the pair potential grows below $`H_{c2}`$. Another purpose is to find out the connection between the gap anisotropy and the extra dHvA amplitude attenuation. Terashima et alTerashima95  reported a dHvA experiment on YNi2B2C where an oscillation is observed to persist down to a field $`0.2H_{c2}`$. On the other hand, a specific-heat experiment at $`H=0`$ shows a power-law behavior $`T^3`$ at low temperatures,Movshovich94  indicating existence of gap anisotropy in this material. Indeed, Izawa et alIzawa02  recently reported presence of four point nodes in the gap based on a thermal-conductivity measurement. Miyake Miyake93  argued that point or line nodes along the extremal orbit may weaken the damping, and proposed to use the dHvA effect as a probe for gap anisotropy. We examine this possibility in full details and present a quantitative theory on the issue.

This paper is organized as follows. Section II provides a formulation to solve the BdG equations for vortex-lattice states. Section III presents calculated quasiparticle spectra and the dHvA oscillation for two-dimensional systems to clarify their basic features as well as the origin of the extra dHvA oscillation damping in the vortex state. In Sec. IV, it is demonstrated that the gap anisotropy at $`H=0`$ can be detectable via the dHvA oscillation in the vortex state based on both numerical and analytic calculations for three-dimensional systems with various gap structures. Section V presents estimations of energy gap for NbSe2, Nb3Sn, and YNi2B2C using the analytic formula obtained in Appendix C. Section VI concludes the paper with a brief summary. In Appendix A, we derive a convenient expression for the thermodynamic potential. Appendix B summarizes the expressions of basis functions and overlap integrals used in the numerical calculations. In Appendix C, we derive an analytic expression for the extra dHvA oscillation damping in the vortex state. A brief report of the contents was already presented in Ref. YK01, .

## II Formulation

### II.1 Bogoliubov-de Gennes equations

Throughout the paper we will rely on the mean-field BdG equations, which obtain the quasiparticle wavefunctions $`𝐮_s`$ and $`𝐯_s^{}`$ with a positive eigenvalue $`E_s>0`$ by

$`{\displaystyle 𝑑𝐫_2[\begin{array}{cc}\hfill \underset{¯}{}(𝐫_1,𝐫_2)& \underset{¯}{\mathrm{\Delta }}(𝐫_1,𝐫_2)\\ \hfill \underset{¯}{\mathrm{\Delta }}^{}(𝐫_1,𝐫_2)& \underset{¯}{}^{}(𝐫_1,𝐫_2)\end{array}][\begin{array}{c}𝐮_s(𝐫_2)\\ 𝐯_s^{}(𝐫_2)\end{array}]}`$ (5)
$`=E_s[\begin{array}{c}𝐮_s(𝐫_1)\\ 𝐯_s^{}(𝐫_1)\end{array}].`$ (8)
Here $`\underset{¯}{\mathrm{\Delta }}`$ is the pair potential and $`\underset{¯}{}`$ denotes the normal-state Hamiltonian in the magnetic field; both are $`2\times 2`$ matrices to describe the spin degrees of freedom. The symbol  denotes Hermitian conjugate in both the coordinate and spin variables as $`[\underset{¯}{\mathrm{\Delta }}^{}(𝐫_1,𝐫_2)]_{\sigma _1\sigma _2}=\mathrm{\Delta }_{\sigma _2\sigma _1}^{}(𝐫_2,𝐫_1)`$ with $`\sigma _j=,`$. With this definition, we can see immediately that the matrix in Eq. (8) is Hermitian.

We adopt the free-particle Hamiltonian for $`\underset{¯}{}`$:

$$\underset{¯}{}(𝐫_1,𝐫_2)=\delta (𝐫_1𝐫_2)\left\{\frac{\left[i\mathrm{}\mathbf{}_2+\frac{e}{c}𝐀(𝐫_2)\right]^2}{2m_\mathrm{e}}\epsilon _\mathrm{F}\right\}\underset{¯}{1},$$
(9)
where $`m_\mathrm{e}`$, $`e`$ $`(e>0)`$, $`c`$, and $`\epsilon _\mathrm{F}`$ are the electron mass, the electron charge, the light velocity, and the chemical potential, respectively. We will not consider the spin paramagnetism throughout. We also neglect the spatial variation of the magnetic field as appropriate for the relevant high-$`\kappa `$ materials. Then, the vector potential $`𝐀`$ can be expressed using the symmetric gauge as

$$𝐀(𝐫)=\frac{1}{2}B\widehat{𝐳}\times 𝐫,$$
(10)
where $`B`$ denotes the average flux density, and we have chosen the field along $`\widehat{𝐳}`$ for convenience.

The pair potential in turn is given with respect to the quasiparticle wavefunctions as

$$\underset{¯}{\mathrm{\Delta }}(𝐫_1,𝐫_2)=V(𝐫_1𝐫_2)\underset{¯}{\mathrm{\Phi }}(𝐫_1,𝐫_2),$$
(11)
where $`V`$ denotes the interaction and $`\underset{¯}{\mathrm{\Phi }}`$ is the order parameter defined by

$`\underset{¯}{\mathrm{\Phi }}(𝐫_1,𝐫_2)`$ $``$ $`{\displaystyle \underset{s}{}}\left[𝐮_s(𝐫_1)𝐯_s^\mathrm{T}(𝐫_2)𝐯_s(𝐫_1)𝐮_s^\mathrm{T}(𝐫_2)\right]`$ (12)
$`\times {\displaystyle \frac{1}{2}}\mathrm{tanh}{\displaystyle \frac{E_s}{2k_\mathrm{B}T}},`$
with $`T`$ the temperature and T denoting the transpose.

It is shown in Appendix A that the thermodynamic potential corresponding to Eqs. (8)-(11) is given by

$`\mathrm{\Omega }=k_\mathrm{B}T{\displaystyle \underset{s}{}}\mathrm{ln}(1+\mathrm{e}^{E_s/k_\mathrm{B}T}){\displaystyle \underset{s}{}}E_s{\displaystyle |𝐯_s(𝐫)|^2𝑑𝐫}`$
$`{\displaystyle \frac{1}{2}}{\displaystyle \mathrm{Tr}\underset{¯}{\mathrm{\Delta }}^{}(𝐫_1,𝐫_2)\underset{¯}{\mathrm{\Phi }}(𝐫_2,𝐫_1)d𝐫_1d𝐫_2},`$ (13)
where Tr denotes taking trace over spin variables. This expression will be useful to obtain an analytic expression for the extra dHvA amplitude attenuation in the vortex state. The magnetization is then calculated by

$$M=\frac{(\mathrm{\Omega }/𝒱)}{B},$$
(14)
where $`𝒱`$ is the volume of the system.

### II.2 Vortex lattices and basic vectors

Solving the above equations for general non-uniform systems is a formidable task. For lattice states, however, it can be reduced into a numerically tractable problem using the symmetry that they are periodic with a single flux quantum $`\varphi _0hc/2e`$ per unit cell. We hence define a pair of basic vectors by

$$\{\begin{array}{c}𝐚_1(a_{1x},a_{1y},0)\hfill \\ 𝐚_2(0,a_2,0)\hfill \end{array},(𝐚_1\times 𝐚_2)\widehat{𝐳}=\frac{\varphi _0}{B}=\pi l_B^2,$$
(15)
where $`l_B\sqrt{\mathrm{}c/eB}`$ is the magnetic length. The basic vectors of the corresponding reciprocal lattice are then defined by

$$\{\begin{array}{c}𝐛_12(𝐚_2\times \widehat{𝐳})/l_B^2\hfill \\ 𝐛_22(\widehat{𝐳}\times 𝐚_1)/l_B^2\hfill \end{array}.$$
(16)
We now introduce magnetic Bloch vectors for quasiparticle eigenstates by Kita98-2

$$𝐤\frac{\mu _1}{𝒩_\mathrm{f}}𝐛_1+\frac{\mu _2}{𝒩_\mathrm{f}}𝐛_2\left(\frac{𝒩_\mathrm{f}}{4}<\mu _j\frac{𝒩_\mathrm{f}}{4}\right),$$
(17)
and those for the center-of-mass coordinate by

$$𝐪\frac{\mu _1}{𝒩_\mathrm{f}}𝐛_1+\frac{\mu _2}{𝒩_\mathrm{f}}𝐛_2\left(\frac{𝒩_\mathrm{f}}{2}<\mu _j\frac{𝒩_\mathrm{f}}{2}\right),$$
(18)
where $`𝒩_\mathrm{f}`$ is an even integer with $`𝒩_\mathrm{f}^2`$ denoting the number of flux quanta in the system. Notice that $`𝐪`$ covers an area four times as large as that of k.

### II.3 Landau-level-expansion method (LLX)

It has been shown Kita98-2  that the pair potential of the conventional Abrikosov lattice can be expanded in two ways with respect to $`(𝐫_1,𝐫_2)`$ and $`(𝐑,𝐫)(\frac{𝐫_1+𝐫_2}{2},𝐫_1𝐫_2)`$ as

$`\underset{¯}{\mathrm{\Delta }}(𝐫_1,𝐫_2)`$
$`={\displaystyle \underset{𝐤\alpha }{}}{\displaystyle \underset{N_1N_2}{}}\underset{¯}{\mathrm{\Delta }}_{N_1N_2}^{(𝐤p_z)}\psi _{N_1𝐤\alpha }(𝐫_1)\psi _{N_2𝐪𝐤\alpha }(𝐫_2){\displaystyle \frac{\mathrm{e}^{ip_z(z_1z_2)}}{L}}`$
$`={\displaystyle \frac{𝒩_\mathrm{f}}{\sqrt{2}}}{\displaystyle \underset{N_\mathrm{c}}{}}{\displaystyle \underset{N_\mathrm{r}mp_z}{}}(1)^{N_\mathrm{r}}\overline{\underset{¯}{\mathrm{\Delta }}}_{N_\mathrm{r}p_z}^{(N_\mathrm{c}m)}\psi _{N_\mathrm{c}𝐪}^{(\mathrm{c})}(𝐑)\psi _{N_\mathrm{r}m}^{(\mathrm{r})}(𝐫){\displaystyle \frac{\mathrm{e}^{ip_zz}}{L}}.`$
(19)
Here $`\psi _{N𝐤\alpha }`$ is a quasiparticle basis function with $`N`$ denoting the Landau level, $`𝐤`$ defined by Eq. (17), and $`\alpha `$ (=1,2) signifying signifying two-fold degeneracy of every orbital state. On the other hand, $`\psi _{N_\mathrm{c}𝐪}^{(\mathrm{c})}`$ and $`\psi _{N_\mathrm{r}m}^{(\mathrm{r})}`$ are basis functions for the center-of-mass and relative coordinates, respectively, with $`N_\mathrm{c}`$ and $`N_\mathrm{r}`$ denoting the corresponding Landau levels, $`𝐪`$ defined by Eq. (18), and $`m`$ an eigenvalue for the relative orbital angular momentum operator $`\widehat{l}_z`$. The quantities $`p_z`$ and $`L`$ are, respectively, the wavenumber and the system length along the $`z`$ direction parallel to the magnetic field; we adopt a notation of using $`𝐩`$ as a wavevector in zero field to distinguish it from the two-dimensional magnetic Bloch vector $`𝐤`$ perpendicular to the field. As noted in Ref. Kita98, , an arbitrary single $`𝐪`$ suffices to describe the conventional Abrikosov lattices due to the broken translational symmetry of the vortex lattice. Then the first expansion of Eq. (19) tells us that, by choosing $`𝐪=\mathrm{𝟎}`$, we get an complete analogy with the uniform system in that the pairing occurs between $`(𝐤,p_z)`$ and $`(𝐤,p_z)`$. Finally, the two expansion coefficients $`\underset{¯}{\mathrm{\Delta }}_{N_1N_2}^{(𝐤p_z)}`$ and $`\overline{\underset{¯}{\mathrm{\Delta }}}_{N_\mathrm{r}p_z}^{(N_\mathrm{c}m)}`$ are connected by

$`\underset{¯}{\mathrm{\Delta }}_{N_1N_2}^{(𝐤p_z)}={\displaystyle \frac{𝒩_\mathrm{f}}{2}}{\displaystyle \underset{N_\mathrm{c}N_\mathrm{r}}{}}N_1N_2|N_\mathrm{c}N_\mathrm{r}{\displaystyle \underset{m}{}}2𝐤𝐪|m+N_\mathrm{r}`$
$`\times (1)^{N_\mathrm{r}}\overline{\underset{¯}{\mathrm{\Delta }}}_{N_\mathrm{r}p_z}^{(N_\mathrm{c}m)},`$ (20a)
$`\overline{\underset{¯}{\mathrm{\Delta }}}_{N_\mathrm{r}p_z}^{(N_\mathrm{c}m)}={\displaystyle \frac{2}{𝒩_\mathrm{f}}}{\displaystyle \underset{N_1N_2}{}}N_\mathrm{c}N_\mathrm{r}|N_1N_2{\displaystyle \underset{𝐤}{}}m+N_\mathrm{r}|2𝐤𝐪`$
$`\times (1)^{N_\mathrm{r}}\underset{¯}{\mathrm{\Delta }}_{N_1N_2}^{(𝐤p_z)},`$ (20b)
where $`N_1N_2|N_\mathrm{c}N_\mathrm{r}`$ and $`2𝐤𝐪|m+N_\mathrm{r}`$ are elements of unitary matrices for the basis change, i.e. overlap integrals. Their explicit expressions together with those for $`\psi _{N𝐤\alpha }`$, $`\psi _{N_\mathrm{c}𝐪}^{(\mathrm{c})}`$, and $`\psi _{N_\mathrm{r}m}^{(\mathrm{r})}`$ are given in Appendix B.

A great advantage of using Eq. (19) is that it enables us to transform Eq. (8) into a numerically tractable problem, as mentioned before. Indeed, expanding the quasiparticle wavefunctions as

$`𝐮(𝐫)={\displaystyle \underset{N𝐤\alpha p_z}{}}𝐮_s(N)\psi _{N𝐤\alpha }(𝐫){\displaystyle \frac{\mathrm{e}^{ip_zz}}{\sqrt{L}}},`$ (21a)
$`𝐯(𝐫)={\displaystyle \underset{N𝐤\alpha p_z}{}}𝐯_s(N)\psi _{N𝐪𝐤\alpha }(𝐫){\displaystyle \frac{\mathrm{e}^{ip_zz}}{\sqrt{L}}},`$ (21b)
Eq. (8) is reduced to a separate matrix eigenvalue problem for each $`𝐤\alpha p_z`$, and the eigenstate is labelled by $`s=\nu 𝐤\alpha p_z\sigma `$ with $`\nu `$ and $`\sigma `$ denoting the quasiparticle band and its spin, respectively. Explicitly, Eq. (8) becomes

$$\underset{N_2}{}[\begin{array}{cc}\hfill \underset{¯}{}_{N_1N_2}^{(p_z)}& \underset{¯}{\mathrm{\Delta }}_{N_1N_2}^{(𝐤p_z)}\\ \hfill \underset{¯}{\mathrm{\Delta }}_{N_1N_2}^{(𝐤p_z)}& \underset{¯}{}_{N_1N_2}^{(p_z)}\end{array}][\begin{array}{c}𝐮_s(N_2)\\ 𝐯_s^{}(N_2)\end{array}]=E_s[\begin{array}{c}𝐮_s(N_1)\\ 𝐯_s^{}(N_1)\end{array}],$$
(22)
where $`\underset{¯}{\mathrm{\Delta }}_{N_1N_2}^{(𝐤p_z)}`$ is given by Eq. (20a), and $`\underset{¯}{}_{N_1N2}^{(p_z)}`$ is diagonal as

$$\underset{¯}{}_{N_1N2}^{(p_z)}=\delta _{N_1N_2}\left[(N_1+\frac{1}{2})\mathrm{}\omega _B+\frac{\mathrm{}^2p_z^2}{2m_\mathrm{e}}\epsilon _\mathrm{F}\right]\underset{¯}{1},$$
(23)
with $`\omega _BeB/Mc`$ the cyclotron frequency.

The self-consistency equation (11) are also simplified greatly. Let us define

$`V_{\mathrm{𝐩𝐩}^{}}{\displaystyle V(𝐫)\mathrm{e}^{\mathrm{i}(𝐩𝐩^{})𝐫}\mathrm{d}^3r}`$
$`=4\pi {\displaystyle \underset{\mathrm{}=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{m=\mathrm{}}{\overset{\mathrm{}}{}}}\overline{V}_{\mathrm{}}(p,p^{})Y_\mathrm{}m(\widehat{𝐩})Y_\mathrm{}m^{}(\widehat{𝐩}^{}),`$ (24)
where $`Y_\mathrm{}m(\widehat{𝐩})\mathrm{\Theta }_\mathrm{}m(\theta _𝐩)\frac{1}{\sqrt{2\pi }}\mathrm{exp}(im\phi _𝐩)`$ is the spherical harmonic.Schiff  We also expand both $`\underset{¯}{\mathrm{\Delta }}`$ and $`\underset{¯}{\mathrm{\Phi }}`$ in terms of the center-of-mass and relative coordinates as the last line of Eq. (19). Then Eq. (11) is transformed into an equation for the expansion coefficients of each $`(N_\mathrm{c},m)`$ as

$$\overline{\underset{¯}{\mathrm{\Delta }}}_{N_\mathrm{r}p_z}^{(N_\mathrm{c}m)}=\frac{1}{2\pi l_B^2L}\underset{\mathrm{}N_\mathrm{r}^{}p_z^{}}{}\overline{V}_{\mathrm{}}(p,p^{})\mathrm{\Theta }_\mathrm{}m(\theta _𝐩)\mathrm{\Theta }_\mathrm{}m(\theta _𝐩^{})\overline{\underset{¯}{\mathrm{\Phi }}}_{N_\mathrm{r}^{}p_z^{}}^{(N_\mathrm{c}m)}$$
(25)
with $`p=\sqrt{N_\mathrm{r}/l_B^2+p_z^2}`$ and $`\theta _𝐩=\mathrm{tan}^1\left(\frac{\sqrt{N_\mathrm{r}}/l_B}{p_z}\right)`$.

Let us further assume that a single $`\mathrm{}`$ is relevant in Eq. (24) and take the corresponding $`\overline{V}_{\mathrm{}}(p,p^{})`$ in a separable form as

$$\overline{V}_{\mathrm{}}(p,p^{})=V_{\mathrm{}}W_{\mathrm{}}(\xi )W_{\mathrm{}}(\xi ^{}),$$
(26)
where $`W_{\mathrm{}}(\xi )`$ is some cut-off function with respect to $`\xi \mathrm{}^2p^2/2m_\mathrm{e}\epsilon _\mathrm{F}`$ satisfying $`W_{\mathrm{}}(0)=1`$. Then can rewrite Eq. (25) as

$$\overline{\underset{¯}{\mathrm{\Delta }}}_{N_\mathrm{r}p_z}^{(N_\mathrm{c}m)}=\underset{¯}{\overset{~}{\mathrm{\Delta }}}^{(N_\mathrm{c}m)}W_{\mathrm{}}(\xi )\mathrm{\Theta }_\mathrm{}m(\theta _𝐩),$$
(27a)
with $`\xi =\mathrm{}^2(N_\mathrm{r}/l_B^2+p_z^2)/2m_\mathrm{e}\epsilon _\mathrm{F}`$ and
$$\underset{¯}{\overset{~}{\mathrm{\Delta }}}^{(N_\mathrm{c}m)}=\frac{V_{\mathrm{}}}{2\pi l_B^2L}\underset{N_\mathrm{r}^{}p_z^{}}{}W_{\mathrm{}}(\xi ^{})\mathrm{\Theta }_\mathrm{}m(\theta _𝐩^{})\overline{\underset{¯}{\mathrm{\Phi }}}_{N_\mathrm{r}^{}p_z^{}}^{(N_\mathrm{c}m)}.$$
(27b)
Thus, we only need a self-consistent solution for a set of discrete parameters $`\{\underset{¯}{\overset{~}{\mathrm{\Delta }}}^{(N_\mathrm{c}m)}(T,B)\}`$ through Eqs. (22) and (27) using Eq. (20).

It has been shown Ryan93 ; Kita98  that retaining a few $`N_\mathrm{c}`$’s, e.g., $`N_\mathrm{c}=0,6,12`$ for the hexagonal lattice, is sufficient to describe the vortex lattices of $`H0.05H_{c2}`$. Thus, the original problem of obtaining self-consistency for $`\underset{¯}{\mathrm{\Delta }}(𝐫_1,𝐫_2)`$ at all space points is now reduced to the one for a few expansion coefficients $`\{\underset{¯}{\overset{~}{\mathrm{\Delta }}}^{(N_\mathrm{c}m)}(T,B)\}`$. This situation is analogous to the zero-field case where a single parameter $`\mathrm{\Delta }_0(T)`$ specifies the pair potential.

The linearized self-consistency equation is obtained by substituting into Eq. (27b) the expression of $`\overline{\underset{¯}{\mathrm{\Phi }}}_{N_\mathrm{r}p_z}^{(N_\mathrm{c}m)}`$ linear in the pair potential:Kita98-2

$`\overline{\underset{¯}{\mathrm{\Phi }}}_{N_\mathrm{r}p_z}^{(N_\mathrm{c}m)}={\displaystyle \frac{1}{2}}{\displaystyle \underset{N_1N_2}{}}N_\mathrm{c}N_\mathrm{r}|N_1N_2{\displaystyle \frac{\mathrm{tanh}\frac{\xi _1}{2T}+\mathrm{tanh}\frac{\xi _2}{2T}}{\xi _1+\xi _2}}`$
$`\times {\displaystyle \underset{n}{}}(1)^nN_1N_2|N_\mathrm{c}+nN_\mathrm{r}n\overline{\underset{¯}{\mathrm{\Delta }}}_{N_\mathrm{r}np_z}^{(N_\mathrm{c}+nm+n)}.`$ (28)
Equation (27) with Eq. (28) determines the mean-field $`H_{c2}(T)`$, i.e. $`T_c(H)`$. If we use the asymptotic expression (78) for the overlap integral $`N_\mathrm{c}N_\mathrm{r}|N_1N_2`$ and replace the sum over $`N_1`$ by the integral over $`x(N_1N_2)/\sqrt{2(N_\mathrm{c}+N_\mathrm{r})}`$, we reproduce the smooth quasiclassical $`H_{c2}^{\mathrm{quasi}}(T)`$ obtained, for example, for the $`s`$-wave pairing by Helfand and Werthamer.HW66

Finally, Eq. (11) for two-dimensional systems can be transformed similarly. It is also obtained from Eqs. (24)-(27) by replacing $`V_{\mathrm{}}(p,p^{})V^{(m)}(p,p^{})`$, extending the summation over $`m`$ in Eq. (24) from $`\mathrm{}`$ to $`\mathrm{}`$, and finally restricting the summation over $`\mathrm{}`$ and $`p_z`$ only to $`\mathrm{}=0`$ and $`p_z=0`$, respectively.

## III Two-Dimensional calculations

We first consider a couple of two-dimensional models and perform fully self-consistent calculations. Our purposes in this section are summarized as follows: (i) to clarify the essential features of the results by self-consistent calculations; (ii) to see whether point nodes in the gap really enhances the dHvA signals as Miyake claims.Miyake93

### III.1 Models

The one-particle Hamiltonian (9) yields an isotropic Fermi surface specified by a unit vector $`\widehat{𝐩}=(\mathrm{cos}\phi _𝐩,\mathrm{sin}\phi _𝐩)`$. As for the pairing interaction (24), we adopt the following models:

$$V_{\mathrm{𝐩𝐩}^{}}=\{\begin{array}{c}V_0W(\xi )W(\xi ^{})\hfill \\ V_2W(\xi )(\widehat{𝐩}_x^2\widehat{𝐩}_y^2)W(\xi ^{})(\widehat{𝐩}_x^2\widehat{𝐩}_y^2)\hfill \end{array},$$
(29)
where

$$W(\xi )=\mathrm{exp}\left[\frac{1}{2}\left(\frac{\xi }{\mathrm{}\omega _\mathrm{D}}\right)^4\right]$$
(30)
is a smooth cut-off function with $`\omega _\mathrm{D}`$ a cut-off frequency. The second model of Eq. (29) is beyond the original isotropic interaction (24), but it is convenient for the above-mentioned purposes. In zero field, the two interactions yield the $`s`$\- and $`d_{x^2y^2}`$-wave gaps as

$$\underset{¯}{\mathrm{\Delta }}_𝐩=\{\begin{array}{cc}\mathrm{\Delta }_0W(\xi )i\underset{¯}{\sigma }_2\hfill & \text{}s\text{-wave}\hfill \\ \mathrm{\Delta }_0W(\xi )(\widehat{𝐩}_x^2\widehat{𝐩}_y^2)i\underset{¯}{\sigma }_2\hfill & \text{}d_{x^2y^2}\text{-wave}\hfill \end{array},$$
(31)
respectively. The corresponding $`\overline{\underset{¯}{\mathrm{\Delta }}}_{N_\mathrm{r}}^{(N_\mathrm{c}m)}`$ of Eq. (27) for $`𝐁\widehat{𝐳}`$ is given by

$$\overline{\underset{¯}{\mathrm{\Delta }}}_{N_\mathrm{r}}^{(N_\mathrm{c}m)}=\{\begin{array}{cc}\stackrel{~}{\mathrm{\Delta }}^{(N_\mathrm{c})}W(\xi )\delta _{m0}i\underset{¯}{\sigma }_2\hfill & \\ \stackrel{~}{\mathrm{\Delta }}^{(N_\mathrm{c})}W(\xi )\frac{\delta _{m2}+\delta _{m2}}{2}i\underset{¯}{\sigma }_2\hfill & \end{array},$$
(32a)
with $`\xi \mathrm{}^2N_\mathrm{r}/2m_\mathrm{e}l_B^2\epsilon _\mathrm{F}`$ and
$$\stackrel{~}{\mathrm{\Delta }}^{(N_\mathrm{c})}=\{\begin{array}{c}\frac{V_0}{4\pi l_B^2}\underset{N_\mathrm{r}^{}}{}W_{\mathrm{}}(\xi ^{})\overline{\mathrm{\Phi }}_{N_\mathrm{r}^{}}^{(N_\mathrm{c},0)}\hfill \\ \frac{V_2}{4\pi l_B^2}\underset{N_\mathrm{r}^{}}{}W_{\mathrm{}}(\xi ^{})\frac{\overline{\mathrm{\Phi }}_{N_\mathrm{r}^{}}^{(N_\mathrm{c},2)}+\overline{\mathrm{\Phi }}_{N_\mathrm{r}^{}}^{(N_\mathrm{c},2)}}{2}\hfill \end{array}.$$
(32b)
Here we have adopted a normalization for $`\stackrel{~}{\mathrm{\Delta }}^{(N_\mathrm{c})}`$ different from Eq. (27) so that $`\stackrel{~}{\mathrm{\Delta }}^{(N_\mathrm{c})}`$ acquires a direct correspondence to the maximum gap $`\mathrm{\Delta }_0`$ in Eq. (31); the factor $`\frac{1}{2}`$ in the second case stems from $`\mathrm{cos}2\phi _𝐩`$ in Eq. (31).

We have chosen $`V_{\mathrm{}}`$ in Eq. (29) as

$$g_{\mathrm{}}N(0)V_{\mathrm{}}=0.5,$$
(33)
where $`N(0)=m_e/2\pi \mathrm{}^2`$ is the density of states per spin at the Fermi level. Another important parameter is the zero-temperature coherence length defined by

$`\xi _0\mathrm{}v_\mathrm{F}/\mathrm{\Delta }_0,`$ (34)
with $`v_\mathrm{F}`$ the Fermi velocity. We have adopted $`p_\mathrm{F}\xi _0=5`$ for our calculations. The above two quantities fix our models completely; the cut-off $`\mathrm{}\omega _\mathrm{D}`$ in Eq. (30) and $`T_c(B=0)`$ are calculated using the gap equation.

It should be noted that choosing $`p_\mathrm{F}\xi _0`$ also determines the following quantities: (i) the ratio $`\mathrm{}\omega _{H_{c2}^{\mathrm{quasi}}}/k_\mathrm{B}T_c`$, where $`\mathrm{}\omega _{H_{c2}^{\mathrm{quasi}}}`$ is the zero-temperature cyclotron energy at the quasiclassical upper critical field $`H_{c2}^{\mathrm{quasi}}`$; (ii) the number $`N_\mathrm{F}`$ of the Landau levels below the Fermi level at $`H_{c2}^{\mathrm{quasi}}`$. Indeed, using the usual cut-off model $`W(\xi )=\theta (\mathrm{}\omega _\mathrm{D}|\xi _𝐩|)`$ with $`\theta `$ the step function, we obtain the following results for the $`s`$-wave pairing:

$`\{\begin{array}{c}\mathrm{}\omega _{H_{c2}^{\mathrm{quasi}}}/k_\mathrm{B}T_c=6.28/p_\mathrm{F}\xi _0,\hfill \\ N_\mathrm{F}\epsilon _\mathrm{F}/\mathrm{}\omega _{H_{c2}^{\mathrm{quasi}}}=0.140(p_\mathrm{F}\xi _0)^2.\hfill \end{array}`$ (37)
To reproduce the experimental situation $`\mathrm{}\omega _{H_{c2}^{\mathrm{quasi}}}/k_\mathrm{B}T_c=13`$ within the present model, we should go into the quantum limit $`p_\mathrm{F}\xi _0=62`$, but we then have $`N_\mathrm{F}=51`$. In real materials, however, $`\mathrm{}\omega _{H_{c2}^{\mathrm{quasi}}}/k_\mathrm{B}T_c`$ and $`p_\mathrm{F}\xi _0`$ are apparently independent parameters due to effects not covered by the free-particle model such as the energy band structure. Indeed, $`p_\mathrm{F}\xi _0`$ is of the order of $`30`$ (NbSe2) or even larger, whereas $`\mathrm{}\omega _{H_{c2}^{\mathrm{quasi}}}/k_\mathrm{B}T_c=13`$; see Table 1 below. Also, $`N_\mathrm{F}1`$ for those materials. The failures to describe these situations are among the main difficulties of the free-particle model of Eq. (9).

Motivated by these observations, we also perform another calculations with much more Landau levels below the Fermi level. This is achieved by including the effect of the band dispersion. Specifically, we apply the Onsager-Lifshiz (OL) quantization scheme to $`\underset{¯}{}`$ of Eq. (8), i.e. the procedure which has been very successful for describing the dHvA oscillations in the normal state.LK55 ; Shoenberg84  Given the density of states per spin $`N(\epsilon )`$ and the average flux density $`B`$, the $`N`$th Landau level $`\epsilon _N`$ ($`N=0,1,2,\mathrm{}`$) is determined by

$`2(N+{\displaystyle \frac{1}{2}}){\displaystyle \frac{\mathrm{}^2}{l_B^2}}=4\pi \mathrm{}^2{\displaystyle _0^{\epsilon _N}}N(\epsilon ^{})𝑑\epsilon ^{}.`$ (38)
We fix $`\underset{¯}{}_{N_1N_2}^{(p_z)}`$ of Eq. (22) in this way assuming it is diagonal. In contrast, we use the same expression for $`\underset{¯}{\mathrm{\Delta }}_{N_1N_2}^{(𝐤p_z)}`$ of Eq. (22) as the free-particle case. Finally, we choose $`\xi =\epsilon _{N_\mathrm{r}/2}\epsilon _\mathrm{F}`$ in Eq. (32) based on the consideration of the free-particle model.Kita98-2  With these prescriptions together with the transformation (20), the coupled equations (22) and (32) are defined unambiguously. We adopt the model density of states:

$`N(\epsilon )={\displaystyle \frac{m_\mathrm{e}}{2\pi \mathrm{}^2}}\left(1+{\displaystyle \frac{\alpha \mathrm{\Gamma }}{\epsilon ^2+\mathrm{\Gamma }^2}}\right),`$ (39)
and choose the numerical constants $`(\alpha ,\mathrm{\Gamma })=(2.1,2.7)`$ and $`(5.0,1.0)`$ for the $`s`$\- and $`d`$-wave models of Eq. (31), respectively. We also use Eq. (33) for the pairing interaction and fix $`\mathrm{}\omega _\mathrm{D}=0.5\epsilon _\mathrm{F}`$ in Eq. (30). We thereby obtain $`\mathrm{}\omega _{H_{c2}^{\mathrm{quasi}}}k_\mathrm{B}T_c`$ at $`T=0`$ and $`N_\mathrm{F}30`$ at $`H=H_{c2}^{\mathrm{quasi}}`$, which describe the experimental situation much better than the free-particle model.

### III.2 Numerical procedures

Coupled equations (22) and (32) are solved iteratively with the help of the transformation (20) to obtain self-consistent $`\{\stackrel{~}{\mathrm{\Delta }}^{(N_\mathrm{c})}\}`$’s and the corresponding quasiparticle eigenstates. The hexagonal (square) vortex lattice is assumed for the $`s`$-wave ($`d_{x^2y^2}`$-wave) model, as expected theoretically in high magnetic fieldsWon96 ; Ichioka99  and observed recently in La1.83Sr0.17CuO4+δ.Gilardi02  It should be noted, however, that the precise lattice structure is not important for the theory of the dHvA oscillation in superconductors. We set $`𝐪=\frac{1}{2}(𝐛_1+𝐛_2)`$ in the relevant equations so that a core of the pair potential is located at the origin $`𝐑=\mathrm{𝟎}`$.Kita98  An advantage of this choice is that the corresponding quasiparticle energies have the rotational symmetry of the hexagonal (square) lattice around $`𝐤=\mathrm{𝟎}`$; other choices would shift the rotation axis from the origin. We then perform calculations of Eqs. (22) and (32) for a set of discrete $`𝐤`$’s defined by Eq. (17), where $`𝒩_\mathrm{f}`$ is chosen as a multiple of $`12`$ to include all the high-symmetry points $`\mathrm{\Gamma }`$, $`M`$, and $`K`$ ($`\mathrm{\Gamma }`$, $`X`$, and $`M`$) of the hexagonal (square) lattice. Three different values $`𝒩_\mathrm{f}=12,24,36`$ are used to see the size dependence, and it has been checked that the results do not differ for the three cases. The hexagonal (square) symmetry enables us to restrict the summation over $`𝐤`$ into approximately $`1/12`$ ($`1/8`$) area of the Brillouin zone. Thus, the calculations can be reduced greatly with due care on the degeneracy of high-symmetry points. Finally, the obtained eigenvalues and eigenstates are substituted into Eq. (13) to calculate the magnetization by Eq. (14). All the calculations are performed at $`T=0.1T_c`$.

### III.3 Results

The above self-consistent procedure is known to give rise to oscillatory singular behaviors in both $`H_{c2}`$ and $`\stackrel{~}{\mathrm{\Delta }}^{(N_\mathrm{c})}`$ in the field range where the dHvA oscillation persists.Rasolt87 ; MacDonald92 ; Rasolt92  Figure 1 displays $`H_{c2}(T)`$ calculated self-consistently for the $`s`$\- and $`d`$-wave models; it is normalized by the quasiclassical upper critical field $`H_{c2}^{\mathrm{quasi}}(T=0)`$. An oscillatory behavior sets in around $`T0.2T_c`$, and $`H_{c2}`$ deviates substantially from the smooth Helfand-Werthamer behaviorHW66  predicted by the quasiclassical theory. The number of the Landau levels below the Fermi level is $`N_\mathrm{F}10`$ around $`H_{c2}^{\mathrm{quasi}}(0)`$, which is considerably smaller than those for the real materials. Figure 2 shows $`\stackrel{~}{\mathrm{\Delta }}^{(N_\mathrm{c})}`$ as a function of $`B`$ at $`T=0.1T_c`$ for the $`s`$-wave hexagonal lattice (first column) and the $`d`$-wave square lattice (second column); they are real and finite only for $`N_\mathrm{c}=0,6,12,\mathrm{}`$ $`(0,4,8,\mathrm{})`$ for the hexagonal (square) lattice,Ryan93 ; Kita98  as already mentioned. We observe that $`\stackrel{~}{\mathrm{\Delta }}^{(N_\mathrm{c})}`$’s are also singular, and the dominant $`\stackrel{~}{\mathrm{\Delta }}^{(0)}`$ component cannot be described by the square-root behavior near $`H_{c2}^{\mathrm{quasi}}(0)`$ expected from the quasiclassical theory. However, those singular behaviors disappear gradually as $`B`$ decreases.

Figure 3 displays the quasiparticle energies in the magnetic Brillouin zone for the $`s`$-wave hexagonal lattice (first column) and the $`d`$-wave square lattice (second column) at $`B/H_{c2}^{\mathrm{quasi}}(0)=0.8`$, $`0.5`$, and $`0.1`$. At $`B/H_{c2}^{\mathrm{quasi}}(0)=0.8`$, we already observe large dispersion for $`E2\mathrm{\Delta }_0`$ where the pair potential is effective. In contrast, the flat Landau-level structure remains for $`E2\mathrm{\Delta }_0`$ where the pair potential vanishes in the present cut-off model of Eq. (30). Thus, the dispersion is caused clearly by the scattering from the growing pair potential, and as will be discussed below, it is the origin of the extra dHvA oscillation damping in the vortex state. We also notice that, for $`B/H_{c2}^{\mathrm{quasi}}(0)0.5`$, almost no qualitative difference can be seen between the $`s`$\- and $`d`$-wave cases. At a lower field of $`B/H_{c2}^{\mathrm{quasi}}(0)=0.1`$, however, a marked difference grows around $`E\mathrm{\Delta }_0`$. The $`s`$-wave energy bands of $`E0.7\mathrm{\Delta }_0`$ are flat and occur in pairs with the level spacing of the order of $`\mathrm{\Delta }_0^2/\epsilon _\mathrm{F}`$. As already pointed out by Norman et al.,NMA95  these corresponds to the bound core states of an isolated vortex with little tunneling probability between adjacent cores. In contrast, the $`d`$-wave bands in the same region are densely packed with large dispersion, indicating the extended nature of the corresponding quasiparticle wavefunctions. From this comparison, we conclude that no bound states exist for the $`d`$-wave model even in the zero-field limit of an isolated vortex, in agreement with the result of Franz and Tešanović.Franz98  This difference in the low-energy dispersion at low fields was already reported in Ref. Yasui99, .

Figure 4 shows oscillatory part of the magnetization $`M_{\mathrm{osc}}`$ calculated numerically by Eq. (14), where curves of the corresponding normal state are also plotted for comparison. The damping starts from above $`H_{c2}^{\mathrm{quasi}}(0)`$ where $`\stackrel{~}{\mathrm{\Delta }}^{(0)}`$ is already finite as in Fig. 2, and develops rapidly as $`\stackrel{~}{\mathrm{\Delta }}^{(0)}`$ grows in decreasing $`B`$. Thus, the mean-field theory predicts that the dHvA oscillation comes together with the oscillatory singular behaviors in $`H_{c2}`$ and $`\stackrel{~}{\mathrm{\Delta }}^{(N_\mathrm{c})}`$. Combined with the energy dispersion given in Fig. 3, we are now able to attribute the origin of the extra damping unambiguously to the Landau-level broadening due to the pair potential. The oscillations are rather irregular in both the $`s`$-wave and $`d`$-wave cases, in accordance with the singular behaviors of $`\mathrm{\Delta }^{(N_\mathrm{c})}`$ in Fig. 2. We also see no qualitative difference between the two cases.

However, the free-particle model has several inappropriate points as discussed already around Eq. (37). For example, the number of Landau levels below the Fermi level $`N_\mathrm{F}`$ is necessarily $`N_\mathrm{F}10`$ at $`H_{c2}^{\mathrm{quasi}}`$ for $`p_\mathrm{F}\xi _0=5`$, which is much smaller than the values of the materials displaying the dHvA oscillation. Hence the above numerical results may not be sufficient to say anything quantitative about the dHvA attenuation or the differences between the $`s`$\- and $`d`$-wave cases. We have thus performed another calculations for the model described around Eqs. (38) and (39) where $`\mathrm{}\omega _{H_{c2}^{\mathrm{quasi}}}/k_\mathrm{B}T_c1`$ and $`N_\mathrm{F}30`$ at $`B=H_{c2}^{\mathrm{quasi}}(0)`$.

Figure 5 shows the field-dependence of the expansion coefficients $`\stackrel{~}{\mathrm{\Delta }}^{(N_\mathrm{c})}`$ calculated self-consistently for the $`s`$-wave hexagonal lattice (first column) and the $`d`$-wave square lattice (second column). Singular oscillatory behaviors are manifest in both cases as in the case of the quadratic dispersion, which originate from the singular density of states of Landau levels.Rasolt92  For example, the dominant $`\stackrel{~}{\mathrm{\Delta }}^{(0)}`$ component have a nonzero value from above $`H_{c2}^{\mathrm{quasi}}`$ and deviates substantially from the quasiclassical square-root behavior (dotted lines).

Figure 6 displays the corresponding oscillatory part of the magnetization $`M_{\mathrm{osc}}`$ calculated numerically by Eq. (14), where normal-state results (dotted lines) are also plotted for comparison. The main features are summarized as follows: (i) The oscillations are seen to decrease from above the quasiclassical $`H_{c2}^{\mathrm{quasi}}`$, due to the reentrant behavior of $`\stackrel{~}{\mathrm{\Delta }}^{(N_\mathrm{c})}`$, to be reduced considerably around $`B0.8H_{c2}^{\mathrm{quasi}}`$, i.e. $`H_{c2}^{\mathrm{quasi}}/B1.25`$. However, they do not disappear completely in lower fields. (ii) This extra attenuation is due to the broadening of the Landau levels caused by the pair potential, as in the case of the quadratic dispersion. Indeed, we have obtained quasiparticle spectra similar to those of Fig. 3. (iii) The period of the oscillation remains unchanged above $`0.8H_{c2}^{\mathrm{quasi}}`$, but some irregularity appears in lower fields. These features are in agreement with the results by Norman et al.NMA95  (iv) Little difference can be seen between the $`s`$\- and $`d`$-wave attenuations.

### III.4 Summary of Two-Dimensional Calculations

Let us summarize results and conclusions from our two-dimensional calculations. (i) Combining Figs. 3 and 4, we are now able to attribute the origin of the extra dHvA oscillation damping unambiguously to the Landau-level broadening due to the scattering by the pair potential. (ii) As may be realized from Fig. 6, presence of point nodes along the extremal orbit does not weaken the attenuation, contrary to the statement by Miyake.Miyake93  This fact suggests that the attenuation is determined by the average gap along the extremal orbit. (iii) The mean-field theory predicts that the dHvA oscillation comes together with the oscillatory behaviors in $`H_{c2}`$ and $`\stackrel{~}{\mathrm{\Delta }}^{(N_\mathrm{c})}`$. This will be so in three dimensional models where $`H_{c2}(T)`$ also shows an oscillatory behavior.Rasolt92  However, such singular behaviors of $`H_{c2}`$ have never been identified definitely in any materials displaying the dHvA oscillation, and reported $`H_{c2}`$ curves show more or less the smooth quasiclassical behavior. This discrepancy between the mean-field theory and the dHvA experiments remains a puzzle to be resolved in the future. (iv) The oscillation attenuates considerably around $`B0.8H_{c2}^{\mathrm{quasi}}`$, i.e. $`H_{c2}^{\mathrm{quasi}}/B1.25`$, although we have set $`\mathrm{}\omega _{H_{c2}^{\mathrm{quasi}}}/k_\mathrm{B}T_c1`$ and $`N_\mathrm{F}1`$ at $`H_{c2}^{\mathrm{quasi}}`$. Thus, the two dimensional models fail to explain the experiment by Terashima et alTerashima95  which shows a persistent oscillation down to $`0.2H_{c2}`$. In addition, the models cannot say anything about whether presence of a line node along the extremal orbit weakens the attenuation. (v) The approximation of retaining only $`\stackrel{~}{\mathrm{\Delta }}^{(0)}`$ works excellently for calculating $`M_{\mathrm{osc}}`$. Indeed, we have checked that the curves of $`M_{\mathrm{osc}}`$ thereby obtained are almost indistinguishable from those of Fig. 6. (vi) The discrepancy mentioned in (iii) above suggests that we should rather use $`\stackrel{~}{\mathrm{\Delta }}^{(0)}`$ obtained quasiclassically to reproduce the smooth behaviors of $`H_{c2}`$ in real materials. Figure 7 plots curves of $`M_{\mathrm{osc}}`$ calculated using quasiclassical $`\stackrel{~}{\mathrm{\Delta }}^{(0)}`$, i.e. the dotted lines of Fig. 5. The oscillations are seen more regular than those of Fig. 6, but the amplitudes attenuate almost similarly and are reduced considerably around $`H_{c2}^{\mathrm{quasi}}/B1.25`$. We hence realize that using the quasiclassical $`\stackrel{~}{\mathrm{\Delta }}^{(0)}`$ suffices for the theory of the oscillation damping. This statement is especially true in the low-field region $`0.2H_{c2}^{\mathrm{quasi}}`$ where $`\stackrel{~}{\mathrm{\Delta }}^{(0)}`$ approaches to the quasiclassical behavior, as may be realized from Fig. 5.

## IV Three-Dimensional calculations

Having clarified basic features of the dHvA oscillation for two-dimensional models as well as the mechanism of the extra oscillation damping, we proceed to consider three-dimensional models with various gap structures which are more relevant to real materials. Our purposes in this section are summarized as follows: (i) to clarify the connection between the extra dHvA oscillation damping and the gap anisotropy by numerical calculations; (ii) to obtain an analytic formula for the extra oscillation damping; (iii) to estimate the gap magnitudes of various materials using the obtained analytic formula.

### IV.1 Model

The one-particle Hamiltonian (9) yields a spherical Fermi surface in the normal state. As for the pairing interaction, we consider three different models:

$`V_{\mathrm{𝐩𝐩}^{}}=\{\begin{array}{c}V_0W(\xi )W(\xi ^{})\hfill \\ V_2W(\xi )(\widehat{𝐩}_x^2\widehat{𝐩}_y^2)W(\xi ^{})(\widehat{𝐩}_x^2\widehat{𝐩}_y^2)\hfill \\ V_1W(\xi )\widehat{𝐩}\widehat{𝐜}W(\xi ^{})\widehat{𝐩}^{}\widehat{𝐜}\hfill \end{array}.`$ (43)
Here $`W(\xi )`$ is a cut-off function given by Eq. (30), $`\widehat{𝐩}(\mathrm{sin}\theta _𝐩\mathrm{cos}\phi _𝐩,\mathrm{sin}\theta _𝐩\mathrm{sin}\phi _𝐩,\mathrm{cos}\theta _𝐩)`$ specifies a point on the Fermi surface, and $`\widehat{𝐜}(\mathrm{sin}\theta _𝐜,0,\mathrm{cos}\theta _𝐜)`$ denotes the direction of the crystal $`c`$-axis, in the coordinate frame where $`𝐁\widehat{𝐳}`$. Again the latter two models of Eq. (43) are beyond the original spherical interaction (24), but they are convenient for the above-mentioned purposes. In zero field, Eq. (43) yield

$`\underset{¯}{\mathrm{\Delta }}_𝐩=\{\begin{array}{cc}\mathrm{\Delta }_0W(\xi )i\underset{¯}{\sigma }_2\hfill & \text{}s\text{-wave}\hfill \\ \mathrm{\Delta }_0W(\xi )(\widehat{𝐩}_x^2\widehat{𝐩}_y^2)i\underset{¯}{\sigma }_2\hfill & \text{}d_{x^2y^2}\text{-wave}\hfill \\ \mathrm{\Delta }_0W(\xi )\widehat{𝐩}\widehat{𝐜}i\underset{¯}{\sigma }_3\underset{¯}{\sigma }_2\hfill & \text{}p_z\text{-wave}\hfill \end{array},`$ (47)
which denote the isotropic $`s`$-wave state, a three-dimensional $`d_{x^2y^2}`$-wave state with four point nodes in the extremal orbit perpendicular to $`𝐁`$, and the $`p`$-wave polar state with a line node perpendicular to $`\widehat{𝐜}`$, respectively. The corresponding $`\overline{\underset{¯}{\mathrm{\Delta }}}_{N_\mathrm{r}p_z}^{(N_\mathrm{c}m)}`$ in Eq. (25) can be written as

$$\overline{\underset{¯}{\mathrm{\Delta }}}_{N_\mathrm{r}p_z}^{(N_\mathrm{c}m)}=\{\begin{array}{cc}\stackrel{~}{\mathrm{\Delta }}^{(N_\mathrm{c})}W(\xi )\delta _{m0}i\underset{¯}{\sigma }_2\hfill & \\ \stackrel{~}{\mathrm{\Delta }}^{(N_\mathrm{c})}W(\xi )\mathrm{sin}^2\theta _𝐩\frac{\delta _{m2}+\delta _{m2}}{2}i\underset{¯}{\sigma }_2\hfill & \\ \stackrel{~}{\mathrm{\Delta }}^{(N_\mathrm{c})}W(\xi )[\mathrm{cos}\theta _𝐩\mathrm{cos}\theta _𝐜\delta _{m0}\hfill & \\ +\mathrm{sin}\theta _𝐩\mathrm{sin}\theta _𝐜\frac{\delta _{m1}+\delta _{m1}}{2}]i\underset{¯}{\sigma }_3\underset{¯}{\sigma }_2\hfill & \end{array},$$
(48a)
where $`\xi \mathrm{}^2(N_\mathrm{r}/l_B^2+p_z^2)/2m_\mathrm{e}\epsilon _\mathrm{F}`$, $`\theta _𝐩\mathrm{tan}^1\left(\frac{\sqrt{N_\mathrm{r}}/l_B}{p_z}\right)`$, and $`\stackrel{~}{\mathrm{\Delta }}^{(N_\mathrm{c})}`$ is defined by
$$\stackrel{~}{\mathrm{\Delta }}^{(N_\mathrm{c})}=\{\begin{array}{c}\frac{V_0}{4\pi l_B^2}\underset{N_\mathrm{r}^{}p_z^{}}{}W_{\mathrm{}}(\xi ^{})\overline{\mathrm{\Phi }}_{N_\mathrm{r}^{}p_z^{}}^{(N_\mathrm{c},0)}\hfill \\ \frac{V_2}{4\pi l_B^2}\underset{N_\mathrm{r}^{}p_z^{}}{}W_{\mathrm{}}(\xi ^{})\mathrm{sin}^2\theta _𝐩^{}\frac{\overline{\mathrm{\Phi }}_{N_\mathrm{r}^{}p_z^{}}^{(N_\mathrm{c},2)}+\overline{\mathrm{\Phi }}_{N_\mathrm{r}^{}p_z^{}}^{(N_\mathrm{c},2)}}{2}\hfill \\ \frac{V_1}{4\pi l_B^2}\underset{N_\mathrm{r}^{}p_z^{}}{}W_{\mathrm{}}(\xi ^{})[\mathrm{cos}\theta _𝐩^{}\mathrm{cos}\theta _𝐜\overline{\mathrm{\Phi }}_{N_\mathrm{r}^{}p_z^{}}^{(N_\mathrm{c},0)}\hfill \\ +\mathrm{sin}\theta _𝐩^{}\mathrm{sin}\theta _𝐜\frac{\overline{\mathrm{\Phi }}_{N_\mathrm{r}^{}p_z^{}}^{(N_\mathrm{c},1)}+\overline{\mathrm{\Phi }}_{N_\mathrm{r}^{}p_z^{}}^{(N_\mathrm{c},1)}}{2}]\hfill \end{array}.$$
(48b)
Here we have adopted a normalization for $`\stackrel{~}{\mathrm{\Delta }}^{(N_\mathrm{c})}`$ different from Eq. (27) so that this quantity acquires a direct correspondence to the maximum gap $`\mathrm{\Delta }_0`$ in Eq. (47). The factors $`\frac{1}{2}`$ in the second and third cases stem from $`\mathrm{cos}2\phi _𝐩`$ and $`\mathrm{cos}\phi _𝐩`$ in Eq. (47), respectively.

The coefficients $`\mathrm{\Delta }^{(N_\mathrm{c})}=\mathrm{\Delta }^{(N_\mathrm{c})}(B,T)`$ in Eq. (48) completely specify the pair potential, as already mentioned. Based on the reasoning given in Sec. III.4(vi), we here adopt a quasiclassical $`\stackrel{~}{\mathrm{\Delta }}^{(N_\mathrm{c})}`$ rather than the fully self-consistent one. Then the dominant $`\stackrel{~}{\mathrm{\Delta }}^{(0)}`$ near $`H_{c2}`$ follows the mean-field square-root behavior to an excellent approximation:

$`\stackrel{~}{\mathrm{\Delta }}^{(0)}=a(1B/H_{c2})^{1/2}.`$ (49)
See the dotted lines in Figs. 2 and 5, for example. In addition, other components $`\stackrel{~}{\mathrm{\Delta }}^{(N_\mathrm{c}>0)}`$ can be neglected for the relevant region $`B0.1H_{c2}`$, as pointed out in Sec. III.4(v). We hence use the lowest-Landau-level approximation of retaining only $`\stackrel{~}{\mathrm{\Delta }}^{(0)}`$. The coefficient $`a=a(T)`$ in Eq. (49) is determined by requiring that the maximum of

$`{\displaystyle \frac{1}{𝒱}}{\displaystyle 𝑑𝐑\left|𝑑𝐫\mathrm{\Delta }(𝐫_1,𝐫_2)\mathrm{e}^{i𝐩𝐫/\mathrm{}}\right|^2}`$ (50)
be equal to $`\mathrm{\Delta }_0^2(1B/H_{c2})`$, where $`\mathrm{\Delta }_0(T)`$ denotes the maximum gap obtained from the weak-coupling theory. This procedure yields

$`a\sqrt{0.5},`$ (51)
for all the three cases of Eq. (48). Substituting Eq. (49) into Eq. (20a) with the choice $`\mathrm{}\omega _\mathrm{D}10\mathrm{\Delta }_0(T=0)`$, the off-diagonal elements of Eq. (22) are fixed completely.

The above non-self-consistent procedure has another advantage that we can choose $`\mathrm{}\omega _{B=H_{c2}}`$ and $`\epsilon _\mathrm{F}`$ in Eq. (23) independently. We have set

$$\mathrm{}\omega _{H_{\mathrm{c2}}}=k_\mathrm{B}T_c\text{at }T=0,$$
(52)
in accordance with $`\mathrm{}\omega _{H_{\mathrm{c2}}}/k_\mathrm{B}T_c=13`$ and $`N_\mathrm{F}1`$ for relevant materials (see Table 1 below). Also, we have chosen $`\epsilon _\mathrm{F}`$ in such a way that there are about $`50`$ Landau levels below $`\epsilon _\mathrm{F}`$ for the extremal orbit at $`H_{c2}`$. Now, the matrix elements of Eq. (22) are specified completely. Hexagonal, square, and hexagonal lattices are assumed for the three cases of Eq. (48), respectively.

### IV.2 Numerical Procedures

The wavevector $`p_z`$ of $`0p_z1.2p_\mathrm{F}`$ is discretized into $`1000`$ points with an equal interval. For each of them, we have diagonalized Eq. (22) with the same procedure as described in Sec. III.2. The obtained results are substituted into Eq. (13) to calculate the magnetization by Eq. (14). All the calculations are performed at $`T=0`$.

### IV.3 Numerical Results

We first focus on the dHvA oscillation of the $`s`$-wave model in Eq. (47) to clarify its basic features. Using our numerical data, we also test the applicability of various theoretical formulas presented so far.

Figure 8(a) presents oscillation of the $`s`$-wave magnetization (blue line) as compared with the normal-state one (sky-blue line). With $`\mathrm{}\omega _{H_{c2}}=k_\mathrm{B}T_c`$, the oscillation is observed to persist down to a fairly low field of $`H_{c2}/B1.8`$, i.e., $`B0.55H_{c2}`$, which is smaller than $`0.8H_{c2}`$ around which $`\mathrm{}\omega _B`$ becomes equal to the spatial average of the energy gap, Eq. (50). This is partly because the gap is smaller within the extremal orbit, as shown quasiclassically by Brandt et al.Brandt67  Indeed, Fig. 9 calculated at $`B=0.968H_{c2}`$ demonstrates that the dispersion for $`p_z=0`$ is smaller than that for $`p_z=0.9p_\mathrm{F}`$. This tendency remains in the high-field region of $`B0.5H_{c2}`$.

It has become conventional to express this extra attenuation in the vortex state by introducing an additional factor $`R_\mathrm{s}`$ for the dHvA oscillation amplitude:

$$R_\mathrm{s}=\mathrm{exp}\left(\frac{\pi }{\omega _B\tau _\mathrm{s}}\right)=\mathrm{exp}\left(\frac{2\pi ^2k_\mathrm{B}T_\mathrm{\Delta }}{\mathrm{}\omega _B}\right),$$
(53)
where the parameters $`\tau _\mathrm{s}`$ and $`T_\mathrm{\Delta }`$ are directly connected with the extra Landau-level broadening $`\mathrm{\Gamma }_\mathrm{s}`$ in the vortex state as $`\mathrm{\Gamma }_\mathrm{s}=\mathrm{}/2\tau _\mathrm{s}=\pi k_\mathrm{B}T_\mathrm{\Delta }`$. The points with error bars in Fig. 8(b) shows $`\mathrm{ln}R_\mathrm{s}`$ as a function of $`1/B`$, i.e. the Dingle plot, obtained by numerical differentiation. This extra damping at high fields shows the behavior $`1B/H_{c2}`$ in the logarithmic scale, but irregularity sets in around $`0.55H_{c2}`$ where the oscillation disappears. We attribute this irregularity to the effect of the bound-state formation in the core region.

The lines in Fig. 8(b) are the predictions from various theoretical formulas. Maki’s formulaMaki91  reproduces the correct functional behavior $`1B/H_{c2}`$ at high fields, but the prefactor is seen too large. The NMA formula,NMA95  deduced from the two-dimensional self-consistent numerical results with $`N_\mathrm{F}10`$ at $`H_{c2}`$, predicts a more rapid attenuation incompatible with our numerical data. One reason for this discrepancy may originate from the fact that their numerical data with $`N_\mathrm{F}10`$ at $`H_{c2}`$ are not appropriate for obtaining an analytic formula by fitting. Another may be attributed to the difference in dimensions. Indeed, the dHvA oscillation in three dimensions differs from that in two dimensions on the point that some finite region $`\delta p_z`$ around the extremal orbit is relevant. Most of the Landau levels in the region do not satisfy the particle-hole symmetry with respect to $`\epsilon _\mathrm{F}`$, so that the effect of the pair potential becomes smaller than that in two dimensions. Another theory by Dukan and Tešanović,Dukan95  which would predict $`R_\mathrm{s}=0`$ in the clean limit of $`T=0`$, is also inconsistent with the data.

The red line in Fig. 8(b) is due to our formula for the extra Dingle temperature:

$`k_\mathrm{B}T_\mathrm{\Delta }=0.5\stackrel{~}{\mathrm{\Gamma }}|\mathrm{\Delta }_𝐩|^2_{\mathrm{eo}}{\displaystyle \frac{m_\mathrm{b}c}{\pi e\mathrm{}}}{\displaystyle \frac{1B/H_{c2}}{B}},`$ (54)
which is derived in Appendix C based on the second-order perturbation with respect to the pair potential. Here $`|\mathrm{\Delta }_𝐩|^2_{\mathrm{eo}}`$ denotes the average gap along the extremal orbit at $`B=0`$, and $`m_\mathrm{b}`$ is the band mass. The numerical constant $`0.5`$ stems from Eq. (51), and $`\stackrel{~}{\mathrm{\Gamma }}`$ is a dimensionless quantity characterizing the Landau-level broadening due to the pair potential. This unknown parameter $`\stackrel{~}{\mathrm{\Gamma }}`$ is determined by a best fit to the $`s`$-wave numerical data, i.e. the points with error bars in Fig. 8(b). This procedure yields

$`\stackrel{~}{\mathrm{\Gamma }}=0.125.`$
We observe in Fig. 8(b) that Eq. (54), which predicts the dependence $`1B/H_{c2}`$ for $`\mathrm{ln}R_\mathrm{s}`$, agrees with the numerical results. This formula will be seen below to reproduce other numerical data excellently without any adjustable parameters.

A difference of Eq. (54) from Maki’s formulaMaki91  lies in the prefactor where the Fermi velocity $`v_\mathrm{F}`$ is absent. Indeed, a dimensional analysis on the second-order perturbation tells us that the Landau-level broadening in the vortex state should be of order $`|\stackrel{~}{\mathrm{\Delta }}^{(0)}(B)|^2/\mathrm{}\omega _B`$, where $`\stackrel{~}{\mathrm{\Delta }}^{(0)}(B)\sqrt{|\mathrm{\Delta }_𝐩|^2_{\mathrm{eo}}(1B/H_{c2})}`$ is essentially the average gap along the extremal orbit. This leads to Eq. (54) except for the numerical constant.

We now turn our attention to see how the presence of point nodes affect the dHvA oscillation. Figure 10(a) shows the oscillation of the $`d`$-wave magnetization (blue line) as compared with the normal-state one (sky-blue line). Although the $`d`$-wave gap in Eq. (47) has four point nodes on the Fermi surface along the extremal orbit, the damping is seen strong and not much different from the $`s`$-wave case. From this fact, we may conclude that it is the average gap along the extremal orbit which is relevant for the extra dHvA oscillation damping. Figure 10(b) presents the corresponding Dingle plot (points with error bars), which is compared with the prediction of Eq. (54). The formula with the average gap $`|\mathrm{\Delta }_𝐩|^2_{\mathrm{eo}}`$ reproduces the numerical result for $`H_{c2}/B1.8`$ excellently without adjustable parameters, thereby providing a strong support for the above statement.

This $`d`$-wave result is in disagreement with Miyake’s theory that point nodes in the extremal orbit should weaken the attenuation.Miyake93  Indeed, Miyake’s theory is based on a semiclassical quantization for the expression of the electron number $`N_e`$ at $`B=0`$. Neither his starting point $`N_e(B=0)`$ nor the use of the semiclassical quantization may be justified for describing the dHvA oscillation observed mainly near $`H_{c2}`$.

We finally consider the $`p`$-wave model with a line node in Eq. (47) to double-check the applicability of Eq. (54). The left figures in Fig. 11 display the dHvA oscillation for the line-node model, where the crystal $`c`$-axis is tilted from the magnetic-field direction by $`\theta _𝐜=0`$ (top), $`\theta _𝐜=\pi /6`$ (second), and $`\theta _𝐜=\pi /4`$ (bottom). The damping is seen weakest in the top figure where the gap vanishes along the extremal orbit, but increases gradually as finite gap opens along the orbit for $`\theta _𝐜=0\pi /4`$. These results indicate conclusively that the average gap along the extremal orbit is relevant for the extra dHvA attenuation. However, the non-zero extra damping in the top figure implies that not only the extremal orbit alone but some finite region around it contributes to the extra damping. Theoretically, this corresponds to the fact that we have to perform the Fresnel integral $`_{\mathrm{}}^{\mathrm{}}\mathrm{exp}[i(\sqrt{2\pi N_\mathrm{F}}p_z/p_\mathrm{F})^2]𝑑p_z`$ for obtaining the LK formula in the normal state. Our data show that this off-extremal-orbit contribution cannot be neglected in the case where the gap vanishes exactly and completely at the extremal orbit. However, this off-extremal contribution is expected to become less important where finite gap is present along the extremal orbit. The right figures in Fig. 11 show the corresponding Dingle plot (points with error bars), which is compared with the prediction of Eq. (54). Except for the weak damping of $`\theta _𝐜=0`$ due to the off-extremal-orbit contribution, the formula is observed to reproduce the numerical results excellently.

## V  Estimation of Energy Gap

Our calculations in Sec. IV.3 have clarified that (i) the gap anisotropy can be detected by measuring the extra dHvA oscillation damping in the vortex state, and (ii) Eq. (54) is particularly useful for this purpose. Using the formula, we finally provide quantitative estimations of the average gap along the extremal orbit for several materials displaying the dHvA oscillation in the vortex state. Table 1 summarizes parameters describing three relevant materials. These materials commonly have fairly high $`T_c`$’s, and the ratio $`\mathrm{}\omega _{H_{c2}}/k_\mathrm{B}T_c`$ ranges from $`1`$ to $`3`$. These features seem to be basic conditions for observing the dHvA oscillation in the vortex state. The values for $`\sqrt{|\mathrm{\Delta }_\mathrm{p}|^2_{\mathrm{eo}}}`$ are obtained by applying Eq. (54) to the observed dHvA attenuation in the vortex state. In doing so, we have adopted as $`m_b`$ in Eq. (54) the values from dHvA experiments rather than those from band calculations, as indicated by the theory of Luttinger.Luttinger61  For comparison, we have also listed the values $`\mathrm{\Delta }(0)`$ estimated by a specific-heat experiment,Movshovich94  a far-infrared measurement,Clayman71  a tunneling experiment,Lee75  or a Raman-scattering experiment.Hackl83  Thus, $`\mathrm{\Delta }(0)`$ is expected to represent band- and/or angle-averaged energy gap. As seen in Table 1, the two quantities coincide excellently for NbSe2 and Nb3Sn, indicating uniformly opened gap in these materials. On the other hand, $`\sqrt{|\mathrm{\Delta }_\mathrm{p}|^2_{\mathrm{eo}}}=1.5`$ for the $`\alpha `$ band of YNi2B2C is considerably smaller than $`\mathrm{\Delta }(0)=2.5`$ from a specific-heat experiment.Movshovich94  This fact implies that YNi2B2C have large band- and/or angle-dependent gap anisotropy. Indeed, Bintley et alMeeson02  have recently carried out a detailed dHvA experiment on this material, rotating the field direction and observing the extra attenuation. They have reported a large angle dependence of the attenuation magnitude. They have also pointed out that their result is in agreement with the model with point nodes presented by Izawa et alIzawa02  based on a thermal-conductivity measurement.

## VI  Summary

We have carried out the first three-dimensional numerical calculations on the dHvA oscillation in the vortex state for various gap structures. We have thereby clarified the relation between gap anisotropy and persistence of the oscillation. We have also derived an analytic formula for the extra dHvA attenuation in the vortex state.

Our main results are given by Figs. 8-11 and Eq. (54). Those figures indicate clearly that the extra dHvA attenuation in the vortex state is directly connected with the average gap along the extremal orbit at $`B=0`$. The derived formula (54) have been shown to reproduce the numerical results excellently. Our theory attributes the origin of the extra dHvA damping to the Landau-level broadening caused by the pair potential. Hence the periodicity of the vortex lattice assumed here is almost irrelevant, and the theory is applicable also to the cases with irregularity such as a random array of vortices. Using Eq. (54), we have estimated average gap amplitudes along the extremal orbit for NbSe2, Nb3Sn, and YNi2B2C. The results indicate presence of large gap anisotropy in YNi2B2C.

Thus, we have shown explicitly that the dHvA effect in the vortex state can be a powerful tool to probe the average gap along the extremal orbit. Our results imply that, by rotating the field direction and observing the attenuation amplitude, we can obtain unique information on the band- and/or angle-dependent gap structure. Such an experiment has recently been performed on UPd2Al3 by Inada et alInada99  and on YNi2B2C by Bintley et al., Meeson02  and the latter group indeed has detected large gap anisotropy in the $`ab`$ plane. Equation (54) will be useful in similar experiments for estimating band- and/or angle-dependent gap amplitudes.

###### Acknowledgements.

We are grateful to F. J. Ohkawa for enlightening discussions. This research is supported by Grant-in-Aid for Scientific Research from the Ministry of Education, Culture, Sports, Science, and Technology of Japan.

## Appendix A Thermodynamic Potential

The Luttinger-Ward thermodynamic potential corresponding to Eq. (8) is given by Kita96

$`\mathrm{\Omega }={\displaystyle \frac{k_\mathrm{B}T}{2}}{\displaystyle \underset{n}{}}\mathrm{Tr}\mathrm{ln}\left[\begin{array}{cc}i\epsilon _n& \mathrm{\Delta }\\ \mathrm{\Delta }^{}& ^{}i\epsilon _n\end{array}\right]`$ (57)
$`\times \left[\begin{array}{cc}\mathrm{e}^{i\epsilon _n0_+}& 0\\ 0& \mathrm{e}^{i\epsilon _n0_+}\end{array}\right]{\displaystyle \frac{1}{2}}\mathrm{Tr}\mathrm{\Delta }^{}\mathrm{\Phi },`$ (60)
where we have adopted a compact notation of using $`x𝐫\sigma `$ with $`\mathrm{\Delta }_{\sigma _1\sigma _2}(𝐫_1,𝐫_2)\mathrm{\Delta }(x_1,x_2)`$, etc., and Tr also implies both integration and summation over $`𝐫`$ and $`\sigma `$, respectively. The quantity $`\epsilon _n/\mathrm{}`$ denotes the Matsubara frequency, and $`0_+`$ is an infinitesimal positive constant. Now, Eq. (8) tells us that the first matrix in Eq. (60) can be diagonalized as Kita96

$`\left[\begin{array}{cc}(x,x^{})& \mathrm{\Delta }(x,x^{})\\ \mathrm{\Delta }^{}(x,x^{})& ^{}(x,x^{})\end{array}\right]`$ (63)
$`={\displaystyle \underset{s}{}}[\begin{array}{cc}u_s^{}(x)& v_s(x)\\ v_s^{}(x)& u_s(x)\end{array}][\begin{array}{cc}E_s& 0\\ 0& E_s\end{array}]\left[\begin{array}{cc}u_s(x^{})& v_s(x^{})\\ v_s^{}(x^{})& u_s^{}(x^{})\end{array}\right].`$ (70)
(71)
Substituting Eq. (71) into Eq. (60), the first term on the right-hand side becomes

$`{\displaystyle \frac{k_\mathrm{B}T}{2}}{\displaystyle \underset{ns}{}}{\displaystyle }\mathrm{d}x\{[|u_s(x)|^2\mathrm{e}^{z_n0_+}+|v_s(x)|^2\mathrm{e}^{z_n0_+}]`$
$`\times \mathrm{ln}(E_sz_n)`$
$`+[|v_s(x)|^2\mathrm{e}^{z_n0_+}+|u_s(x)|^2\mathrm{e}^{z_n0_+}]\mathrm{ln}(E_sz_n)\},`$
(72)
with $`z_ni\epsilon _n`$. The summation over $`n`$ are then transformed with a standard techniqueLW60  into a contour integral just above and below the real axis, using $`f(z)(\mathrm{e}^{z/k_\mathrm{B}T}+1)^1`$ and $`f(z)`$ for the terms with $`\mathrm{e}^{z_n0_+}`$ and $`\mathrm{e}^{z_n0_+}`$, respectively. Considering the poles inside the two contours and using $`[|u_s(x)|^2+|v_s(x)|^2]dx=1`$, we obtain Eq. (13).

## Appendix B Basis Functions and Overlap Integrals

We here present explicit expressions for the quantities appearing in Eqs. (19) and (20); see Ref. Kita98-2,  for their detailed derivations. It should be noted that we here adopt the symmetric gauge (10) which is more convenient than the Landau gauge used in Ref. Kita98-2, . Hence there is an extra factor due to the gauge transformation in every expression of the basis functions, such as $`\mathrm{e}^{\mathrm{i}xy/2l_B^2}`$ in Eq. (73) below.

The basis function $`\psi _{N𝐤\alpha }`$ $`(N=0,1,2,\mathrm{};\alpha =1,2)`$ in Eq. (19) is defined by

$`\psi _{N𝐤\alpha }(𝐫)=`$ $`{\displaystyle \underset{n=𝒩_\mathrm{f}/2+1}{\overset{𝒩_\mathrm{f}/2}{}}}\mathrm{e}^{i[k_y(y+l_B^2k_x/2)+na_{1x}(y+l_B^2k_xna_{1y}/2)/l_B^2]}`$ (73)
$`\times \mathrm{e}^{ixy/2l_B^2(xl_B^2k_yna_{1x})^2/2l_B^2+i(\alpha 1)n\pi }`$
$`\times \sqrt{{\displaystyle \frac{a_{1x}/l_B}{2^NN!\sqrt{\pi }𝒮}}}H_N({\displaystyle \frac{xl_B^2k_yna_{1x}}{l_B}}),`$
where $`𝒮\pi l_B^2𝒩_\mathrm{f}^2`$, and $`H_N(x)\mathrm{e}^{x^2}\left(\frac{\mathrm{d}}{\mathrm{d}x}\right)^N\mathrm{e}^{x^2}`$ is the Hermite polynomial.Abramowitz  The basis function $`\psi _{N𝐪}^{(\mathrm{c})}`$ for the center-of-mass coordinates is obtained from Eq. (73) by putting $`𝐤𝐪`$, $`\alpha =1`$, and $`l_Bl_\mathrm{c}l_B/\sqrt{2}`$ as

$`\psi _{N𝐪}^{(\mathrm{c})}(𝐫)=`$ $`{\displaystyle \underset{n=𝒩_\mathrm{f}/2+1}{\overset{𝒩_\mathrm{f}/2}{}}}\mathrm{e}^{i[q_y(y+l_\mathrm{c}^2q_x/2)+na_{1x}(y+l_\mathrm{c}^2q_xna_{1y}/2)/l_\mathrm{c}^2]}`$ (74)
$`\times \mathrm{e}^{ixy/2l_\mathrm{c}^2(xl_\mathrm{c}^2q_yna_{1x})^2/2l_\mathrm{c}^2}`$
$`\times \sqrt{{\displaystyle \frac{a_{1x}/l_\mathrm{c}}{2^NN!\sqrt{\pi }𝒮}}}H_N({\displaystyle \frac{xl_\mathrm{c}^2q_yna_{1x}}{l_\mathrm{c}}}).`$
Finally, the basis function $`\psi _{Nm}`$ for the relative coordinates, which is conveniently chosen as an eigenstate of the orbital angular momentum operator $`\widehat{l}_z`$, is given by

$$\psi _{Nm}^{(\mathrm{r})}(𝐫)=\frac{(1)^N}{\sqrt{2\pi }l_\mathrm{r}}\sqrt{\frac{N!}{(N+m)!}}\zeta ^m\mathrm{e}^{|\zeta |^2/2}L_N^{(m)}(|\zeta |^2),$$
(75)
where $`\zeta (x+\mathrm{i}y)/\sqrt{2}l_\mathrm{r}`$ with $`l_\mathrm{r}\sqrt{2}l_B`$, and $`L_N^{(m)}(x)\frac{1}{N!}\mathrm{e}^xx^m\left(\frac{\mathrm{d}}{\mathrm{d}x}\right)^N\mathrm{e}^xx^{N+m}`$ is the generalized Laguerre polynomial satisfying $`N+m0`$.Abramowitz

We next provide expressions of the overlap integrals in Eq. (20). The first one, which was obtaind by Rajagopal and Ryan,Rajagopal91  is given by

$`N_\mathrm{c}N_\mathrm{r}|N_1N_2`$
$`\delta _{N_1+N_2,N_\mathrm{c}+N_\mathrm{r}}\sqrt{{\displaystyle \frac{N_1!N_2!N_\mathrm{c}!N_\mathrm{r}!}{2^{N_1+N_2}}}}`$
$`\times {\displaystyle \underset{n=\mathrm{max}(0,N_\mathrm{c}N_2)}{\overset{\mathrm{min}(N_1,N_\mathrm{c})}{}}}{\displaystyle \frac{(1)^{N_2N_\mathrm{c}+n}}{n!(N_1n)!(N_\mathrm{c}n)!(N_2N_\mathrm{c}+n)!}},`$
(76)
which satisfies

$$N_\mathrm{c}N_\mathrm{r}|N_1N_2=(1)^{N_\mathrm{r}}N_\mathrm{c}N_\mathrm{r}|N_2N_1,$$
(77a)
$$\underset{N_1N_2}{}N_\mathrm{c}^{}N_\mathrm{r}^{}|N_1N_2N_1N_2|N_\mathrm{c}N_\mathrm{r}=\delta _{N_\mathrm{c}^{}N_\mathrm{c}}\delta _{N_\mathrm{r}^{}N_\mathrm{r}}.$$
(77b)
The asymptotic expression of Eq. (76) for $`N_\mathrm{c}N_1,N_2`$ is given byKita98-2

$`N_\mathrm{c}N_\mathrm{r}|N_1N_2\delta _{N_1+N_2,N_\mathrm{c}+N_\mathrm{r}}(1)^{N_2}\left[{\displaystyle \frac{2}{\pi (N_\mathrm{c}+N_\mathrm{r})}}\right]^{1/4}`$
$`\times \mathrm{e}^{x^2/2}{\displaystyle \frac{H_{N_\mathrm{c}}(x)}{\sqrt{2^{N_\mathrm{c}}N_\mathrm{c}!}}},`$ (78)
with $`x(N_1N_2)/\sqrt{2(N_\mathrm{c}+N_\mathrm{r})}`$. This expression forms the basis for quasiclassical approximations.

The second overlap integral is given by

$$𝐪|m+N_\mathrm{r}=\sqrt{2\pi }l_\mathrm{r}(1)^{m+N_\mathrm{r}}\left[\psi _{m+N_\mathrm{r}𝐪}^{(\mathrm{r})}(\mathrm{𝟎})\right]^{},$$
(79)
where $`\psi _{m+N_\mathrm{r}𝐪}^{(\mathrm{r})}`$ is another basis function of the relative coordinates defined by

$`\psi _{N𝐪}^{(\mathrm{r})}(𝐫)=`$ $`{\displaystyle \underset{n=𝒩_\mathrm{f}/4+1}{\overset{𝒩_\mathrm{f}/4}{}}}\mathrm{e}^{i[q_y(y+l_\mathrm{r}^2q_x/4)/2+2na_{1x}(y+l_\mathrm{r}^2q_x/2na_{1y})/l_\mathrm{r}^2]}`$
$`\times \mathrm{e}^{ixy/2l_\mathrm{r}^2(xl_\mathrm{r}^2q_y/22na_{1x})^2/2l_\mathrm{r}^2}`$
$`\times \sqrt{{\displaystyle \frac{2a_{1x}/l_\mathrm{r}}{2^NN!\sqrt{\pi }𝒮}}}H_N({\displaystyle \frac{xl_\mathrm{r}^2q_y/22na_{1x}}{l_\mathrm{r}}}).`$
## Appendix C Analytic Formula for dHvA oscillation damping

In order to derive an analytic expression for the dHvA oscillation damping in the vortex states, we start from the thermodynamic potential of Eq. (13). The last term $`\frac{1}{2}\text{Tr}\underset{¯}{\mathrm{\Delta }}\underset{¯}{\mathrm{\Phi }}`$ may be expressed solely with respect to the pair potential, so that they can be neglected in the present model to consider the oscillatory part. Since we are interested in the extra damping in the vortex state, we adopt as $`E_s`$ and $`𝐯_s`$ the expressions from the second-order perturbation with respect to $`\underset{¯}{\mathrm{\Delta }}`$. They are obtained as

$$E_{N𝐤\alpha p_z\sigma }=|\xi _{Np_z\sigma }|+\eta _{N𝐤p_z}^{(1)}\mathrm{sign}(\xi _{Np_z\sigma }),$$
(81)
$$|𝐯_{N𝐤\alpha p_z\sigma }(𝐫)|^2𝑑𝐫=\theta (\xi _{Np_z\sigma })+\eta _{N𝐤p_z}^{(2)}\mathrm{sign}(\xi _{Np_z\sigma }),$$
(82)
where $`\theta `$ is the step function, and $`\eta _{N𝐤p_z}^{(n)}`$ is defined by using Eq. (20a) as

$`\eta _{N𝐤p_z}^{(n)}{\displaystyle \underset{N^{}}{}}{\displaystyle \frac{|\mathrm{\Delta }_{NN^{}}^{(𝐤p_z)}|^2}{(\xi _{Np_z}+\xi _{N^{}p_z})^n}}.`$ (83)
The first terms on the right-hand side of Eqs. (81) and (82) are just the normal-state results. The second terms, on the other hand, denote the finite quasiparticle dispersion in the magnetic Brillouin zone and the smearing of the Fermi surface, respectively, due to the scattering by the growing pair potential. It is useful to express $`\eta _{N𝐤p_z}^{(n)}`$ in terms of $`\stackrel{~}{\mathrm{\Delta }}^{(0)}(B)`$ and the cyclotron energy $`\mathrm{}\omega _B`$ of the extremal orbit as

$`\eta _{N𝐤p_z}^{(n)}={\displaystyle \frac{|\stackrel{~}{\mathrm{\Delta }}^{(0)}(B)|^2}{(\mathrm{}\omega _B)^n}}\stackrel{~}{\eta }_{N𝐤p_z}^{(n)}.`$ (84)
The quantity $`\stackrel{~}{\eta }_{N𝐤p_z}^{(n)}`$ thus defined is dimensionless, and we realize that the main $`B`$ dependence in Eq. (84) lies in the prefactor $`|\stackrel{~}{\mathrm{\Delta }}^{(0)}(B)|^2/(\mathrm{}\omega _B)^n`$. The explicit expression of $`\stackrel{~}{\eta }_{N𝐤p_z}^{(n)}`$ is obtained by using Eqs. (20a) and (48). Considering the case $`\theta _𝐜=0`$, for simplicity, it is given by

$`\stackrel{~}{\eta }_{N𝐤p_z}^{(n)}=`$ $`{\displaystyle \frac{𝒩_\mathrm{f}^2}{4}}{\displaystyle \underset{N^{}mm^{}}{}}{\displaystyle \frac{|NN^{}|0N+N^{}|^2N+N^{}+m|2𝐤𝐪}{[N+N^{}2(N_\mathrm{F}+\delta )]^n}}`$ (88)
$`\times 2𝐤𝐪|N+N^{}+m^{}\times \{\begin{array}{c}\delta _{m0}\delta _{m^{}0}\hfill \\ \delta _{m,\pm 2}\delta _{m^{},\pm 2}\mathrm{sin}^4\theta _𝐩\hfill \\ \delta _{m0}\delta _{m^{}0}\mathrm{cos}^2\theta _𝐩\hfill \end{array},`$
where the quantity $`\delta =\delta (B,p_z)`$ ($`|\delta |<1/2`$) specifies the location of $`\epsilon _\mathrm{F}`$ between the two closest Landau levels, and the overlap integrals are defined by Eqs. (76) and (79). The corresponding normalized density of states:

$`D_{Np_z}^{(n)}(\stackrel{~}{\eta }){\displaystyle \frac{2}{𝒩_\mathrm{f}^2}}{\displaystyle \underset{𝐤\alpha }{}}\delta (\stackrel{~}{\eta }\stackrel{~}{\eta }_{N𝐤p_z}^{(n)}),`$ (90)
will play a central role in the following.

Substituting Eqs. (81) and (82) into Eq. (13), we find that the terms containing $`\eta _{N𝐤p_z}^{(2)}`$ may be neglected due to the cancellation between the particle and hole contributions. The remaining term can be transformed with the standard procedure.LK55  We thereby obtain, for the first harmonic of $`\mathrm{\Omega }/V`$, the expression:

$`{\displaystyle \frac{\mathrm{\Omega }_1}{V}}`$ $`={\displaystyle \frac{k_\mathrm{B}T}{2\pi ^2l_B^2}}{\displaystyle \underset{\sigma }{}}{\displaystyle _{1/2}^{\mathrm{}}}𝑑N\mathrm{cos}(2\pi N){\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑p_z{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\stackrel{~}{\eta }`$ (91)
$`\times D_{Np_z}^{(1)}(\stackrel{~}{\eta })\mathrm{ln}\left[1+\mathrm{e}^{(\xi _{Np_z\sigma }+\stackrel{~}{\eta }|\stackrel{~}{\mathrm{\Delta }}^{(0)}(B)|^2/\mathrm{}\omega _B)/k_\mathrm{B}T}\right].`$
The function $`D_{Np_z}^{(1)}(\stackrel{~}{\eta })`$ depends on $`(N,p_z)`$, but may be replaced by a representative one $`\overline{D}_{\mathrm{}}^{(1)}(\stackrel{~}{\eta })`$ to be placed outside the $`N`$ and $`p_z`$ integrals,comment2  where the recovered index $`\mathrm{}`$ specifies the $`s`$-, $`d`$-, or $`p`$-wave case of Eq. (LABEL:etaT). It may also be acceptable to use a Lorenzian for it: $`\overline{D}_{\mathrm{}}^{(1)}(\stackrel{~}{\eta })=\stackrel{~}{\mathrm{\Gamma }}_{\mathrm{}}/\pi (\stackrel{~}{\eta }^2+\stackrel{~}{\mathrm{\Gamma }}_{\mathrm{}}^2)`$.comment3  We thereby obtain an expression for the magnetization which carries an extra damping factor:

$`R_\mathrm{s}(B)`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}\overline{D}_{\mathrm{}}^{(1)}(\stackrel{~}{\eta })\mathrm{exp}\left[2\pi i\stackrel{~}{\eta }|\stackrel{~}{\mathrm{\Delta }}^{(0)}(B)|^2/(\mathrm{}\omega _B)^2\right]𝑑\stackrel{~}{\eta }`$ (92)
$`=\mathrm{exp}\left[2\pi \stackrel{~}{\mathrm{\Gamma }}_{\mathrm{}}|\stackrel{~}{\mathrm{\Delta }}^{(0)}(B)|^2/(\mathrm{}\omega _B)^2\right].`$
Thus, the superconductivity gives rise to an extra Dingle temperature of $`k_\mathrm{B}T_\mathrm{\Delta }\stackrel{~}{\mathrm{\Gamma }}_{\mathrm{}}|\stackrel{~}{\mathrm{\Delta }}^{(0)}(B)|^2/\pi \mathrm{}\omega _B`$, or equivalently, the extra scattering rate of $`\tau _\mathrm{s}^12\pi k_\mathrm{B}T_\mathrm{\Delta }/\mathrm{}`$.

Equation (92) has an advantage that one can trace the origin of the extra dHvA damping definitely to the growing pair potential, which brings finite quasiparticle dispersion in the magnetic Brillouin zone as Eq. (LABEL:etaT), and the corresponding Landau-level broadening as Eq. (90). Moreover, Eq. (LABEL:etaT) reveals that this broadening near $`H_{c2}`$ is closely connected with the zero-field gap structure given by Eq. (47).

There seems to be no analytic way to estimate $`\stackrel{~}{\mathrm{\Gamma }}_s`$, so we fix it through the best fit to the numerical data of Fig. 8(b). Using Eq. (49) with $`a^2=0.5\mathrm{\Delta }_0^2`$, the procedure yields $`\stackrel{~}{\mathrm{\Gamma }}_s=0.125`$, as mentioned before. It is also clear both from Figs. 10 and 11 and from Eq. (LABEL:etaT) that the average gap around the extremal orbit is relevant for the extra attenuation. We hence put $`a^2\stackrel{~}{\mathrm{\Gamma }}_{\mathrm{}}=0.5|\mathrm{\Delta }_𝐩|^2_{\mathrm{eo}}\stackrel{~}{\mathrm{\Gamma }}_s`$. We thereby obtain Eq. (54), which yields excellent fits to the $`d`$\- and $`p`$-wave numerical data without any adjustable parameters, as seen in Figs. 10 and 11.