Drastic enhancement of the thermal Hall angle in a $d$-wave superconductor

TL;DR

This paper provides a microscopic calculation of the thermal Hall angle in $d$-wave superconductors, confirming its drastic enhancement at low temperatures and weak magnetic fields due to quasiparticle behavior near gap nodes.

Contribution

It derives the thermal Hall coefficient microscopically using augmented quasiclassical equations, explaining the enhancement without impurity scattering.

Findings

Thermal Hall angle is significantly enhanced at low temperatures.

Quasiparticle relaxation time increases near gap nodes.

Enhancement may occur in other nodal superconductors.

Abstract

A drastic enhancement of the thermal Hall angle in $d$-wave superconductors was observed experimentally in a cuprate superconductor and in CeCoIn$_5$ at low temperatures and very weak magnetic field [Phys. Rev. Lett. $\bf 86$, 890 (2001); Phys. Rev. B $\bf 72$, 214515 (2005)]. However, to the best of our knowledge, its microscopic calculation has not been performed yet. To study this microscopically, we derive the thermal Hall coefficient in extreme type-II superconductors with an isolated pinned vortex based on the augmented quasiclassical equations of superconductivity with the Lorentz force. Using it, we can confirm that the quasiparticle relaxation time and the thermal Hall angle are enhanced in $d$-wave superconductors without impurities of the resonant scattering because quasiparticles around the gap nodes which become dominant near zero temperature are restricted to the momentum in a specific orientation. This enhancement of the thermal Hall angle may also be observed in other nodal superconductors with large magnetic-penetration depth.

Equations

$${\left[\epsilon {\check{\tau }}_3-\check{\Delta }-\check{\sigma },\check{g}\right]}_{\circ }+i\hslash {\mathbf{v}}_{\mathrm{F}}\cdot \partial \check{g}$$

$$+\frac{i\hslash }{2}\left[e{\mathbf{v}}_{\mathrm{F}}\cdot \mathbf{E}\frac{\partial }{\partial \epsilon }+e\left({\mathbf{v}}_{\mathrm{F}}\times \mathbf{B}\right)\cdot \frac{\partial }{\partial {\mathbf{p}}_{\mathrm{F}}}\right]\left\{{\check{\tau }}_3,\check{g}\right\}=\check{0},$$

$$\check{g}=\left[\begin{array}{cc}{\hat{g}}^{\mathrm{R}}& {\hat{g}}^{\mathrm{K}}\\ \hat{0}& {\hat{g}}^{\mathrm{A}}\end{array}\right],\check{\Delta }=\left[\begin{array}{cc}\hat{\Delta }& \hat{0}\\ \hat{0}& \hat{\Delta }\end{array}\right],$$

$$\check{\sigma }=-\frac{i\hslash }{2\tau }⟨\check{g}{⟩}_{\mathrm{F}}=\left[\begin{array}{cc}{\hat{\sigma }}^{\mathrm{R}}& {\hat{\sigma }}^{\mathrm{K}}\\ \hat{0}& {\hat{\sigma }}^{\mathrm{A}}\end{array}\right],$$

$${\hat{g}}^{\mathrm{R},\mathrm{A}}=\left[\begin{array}{cc}{g}^{\mathrm{R},\mathrm{A}}& -i{f}^{\mathrm{R},\mathrm{A}}\\ i{\bar{f}}^{\mathrm{R},\mathrm{A}}& -{\bar{g}}^{\mathrm{R},\mathrm{A}}\end{array}\right],{\hat{g}}^{\mathrm{K}}=\left[\begin{array}{cc}{g}^{\mathrm{K}}& -i{f}^{\mathrm{K}}\\ -i{\bar{f}}^{\mathrm{K}}& {\bar{g}}^{\mathrm{K}}\end{array}\right],$$

$$\hat{\Delta }=\left[\begin{array}{cc}0& -\Delta \varphi \\ {\Delta }^{\ast }{\varphi }^{\ast }& 0\end{array}\right],$$

$${\check{\tau }}_3=\left[\begin{array}{cc}{\hat{\tau }}_3& \hat{0}\\ \hat{0}& {\hat{\tau }}_3\end{array}\right],{\hat{\tau }}_3=\left[\begin{array}{cc}1& 0\\ 0& -1\end{array}\right].$$

$${\hat{g}}^{\mathrm{A}}=-{\hat{\tau }}_3{\hat{g}}^{\mathrm{R}†}{\hat{\tau }}_3,{\hat{g}}^{\mathrm{K}}={\hat{\tau }}_3{\hat{g}}^{\mathrm{K}†}{\hat{\tau }}_3.$$

$$a\circ b$$

$$\displaystyle\equiv\exp\left[\frac{i\hbar}{2}\left(\frac{\partial}{\partial\varepsilon}\frac{\partial}{\partial t^{\prime}}-\frac{\partial}{\partial t}\frac{\partial}{\partial\varepsilon^{\prime}}\right)\right]a(\varepsilon,t)b(\varepsilon^{\prime},t^{\prime})\bigg{|}_{\varepsilon^{\prime}=\varepsilon,t^{\prime}=t}.$$

$$\displaystyle{\bm{\partial}}\equiv\left\{\begin{array}[]{ll}{\bm{\nabla}}&{\rm on}\ g^{\rm R,A,K},\bar{g}^{\rm R,A,K}\\ \displaystyle{\bm{\nabla}}-i\frac{2e{\bm{A}}}{\hbar}&{\rm on}\ f^{\rm R,A,K},\Delta\\ \displaystyle{\bm{\nabla}}+i\frac{2e{\bm{A}}}{\hbar}&{\rm on}\ \bar{f}^{\rm R,A,K},\bar{\Delta}\end{array}\right..$$

$$\Delta =\frac{{\mathrm{g}}_0}{4i}{\int }_{-{\epsilon }_{\mathrm{c}}}^{{\epsilon }_{\mathrm{c}}}⟨f^{\mathrm{K}}{\varphi }^{\ast }{⟩}_{\mathrm{F}}d\epsilon ,$$

$${\mathbf{j}}_{\mathrm{Q}}=-\frac{N(0)}{2}{\int }_{-\infty }^{\infty }\epsilon ⟨{\mathbf{v}}_{\mathrm{F}}{g}^{\mathrm{K}}{⟩}_{\mathrm{F}}d\epsilon .$$

$$i\hslash {\mathbf{v}}_{\mathrm{F}}\cdot \partial \check{g}=i\hslash {\mathbf{v}}_{\mathrm{F}}\cdot \nabla \check{g}+\left[e{\mathbf{v}}_{\mathrm{F}}\cdot \mathbf{A}{\check{\tau }}_3,\check{g}\right],$$

$$\left[e{\mathbf{v}}_{\mathrm{F}}\cdot \mathbf{A}{\check{\tau }}_3,\check{g}\right]+\frac{i\hslash }{2}e{\mathbf{v}}_{\mathrm{F}}\cdot \mathbf{E}\frac{\partial }{\partial \epsilon }\left\{{\check{\tau }}_3,\check{g}\right\}={\left[e{\mathbf{v}}_{\mathrm{F}}\cdot \mathbf{A}{\check{\tau }}_3,\check{g}\right]}_{\circ },$$

$${\left[(\epsilon +e{\mathbf{v}}_{\mathrm{F}}\cdot \mathbf{A}){\check{\tau }}_3-\check{\Delta }-\check{\sigma },\check{g}\right]}_{\circ }+i\hslash {\mathbf{v}}_{\mathrm{F}}\cdot \nabla \check{g}$$

$$+\frac{i\hslash }{2}e\left({\mathbf{v}}_{\mathrm{F}}\times \mathbf{B}\right)\cdot \frac{\partial }{\partial {\mathbf{p}}_{\mathrm{F}}}\left\{{\check{\tau }}_3,\check{g}\right\}=\check{0}.$$

$$\check{h}=\check{\Delta }+\check{\sigma }-e{\mathbf{v}}_{\mathrm{F}}\cdot \mathbf{A}{\check{\tau }}_3=\left[\begin{array}{cc}{\hat{h}}^{\mathrm{R}}& {\hat{h}}^{\mathrm{K}}\\ \hat{0}& {\hat{h}}^{\mathrm{A}}\end{array}\right],$$

$${\left[\epsilon {\hat{\tau }}_3-{\hat{h}}^{\mathrm{R},\mathrm{A}},{\hat{g}}^{\mathrm{R},\mathrm{A}}\right]}_{\circ }+i\hslash {\mathbf{v}}_{\mathrm{F}}\cdot \nabla {\hat{g}}^{\mathrm{R},\mathrm{A}}$$

$$+\frac{i\hslash }{2}e\left({\mathbf{v}}_{\mathrm{F}}\times \mathbf{B}\right)\cdot \frac{\partial }{\partial {\mathbf{p}}_{\mathrm{F}}}\left\{{\hat{\tau }}_3,{\hat{g}}^{\mathrm{R},\mathrm{A}}\right\}=\hat{0},$$

$${\left[\epsilon {\hat{\tau }}_3,{\hat{g}}^{\mathrm{K}}\right]}_{\circ }-{\hat{h}}^{\mathrm{R}}\circ {\hat{g}}^{\mathrm{K}}-{\hat{h}}^{\mathrm{K}}\circ {\hat{g}}^{\mathrm{A}}+{\hat{g}}^{\mathrm{R}}\circ {\hat{h}}^{\mathrm{K}}+{\hat{g}}^{\mathrm{K}}\circ {\hat{h}}^{\mathrm{A}}$$

$$+i\hslash {\mathbf{v}}_{\mathrm{F}}\cdot \nabla {\hat{g}}^{\mathrm{K}}+\frac{i\hslash }{2}e\left({\mathbf{v}}_{\mathrm{F}}\times \mathbf{B}\right)\cdot \frac{\partial }{\partial {\mathbf{p}}_{\mathrm{F}}}\left\{{\hat{\tau }}_3,{\hat{g}}^{\mathrm{K}}\right\}=\hat{0}.$$

$$\left[\epsilon {\hat{\tau }}_3-\hat{\Delta }-{\hat{\sigma }}^{\mathrm{R},\mathrm{A}},{\hat{g}}^{\mathrm{R},\mathrm{A}}\right]+i\hslash {\mathbf{v}}_{\mathrm{F}}\cdot \nabla {\hat{g}}^{\mathrm{R},\mathrm{A}}$$

$$+\frac{i\hslash }{2}e\left({\mathbf{v}}_{\mathrm{F}}\times \mathbf{B}\right)\cdot \frac{\partial }{\partial {\mathbf{p}}_{\mathrm{F}}}\left\{{\hat{\tau }}_3,{\hat{g}}^{\mathrm{R},\mathrm{A}}\right\}=\hat{0},$$

$$g_0^{\mathrm{R},\mathrm{A}}$$

$${\tilde{f}}_0^{\mathrm{R},\mathrm{A}}$$

$${\epsilon }^{\mathrm{R}}$$

$${\epsilon }^{\mathrm{A}}$$

$${\Delta }^{\mathrm{R},\mathrm{A}}$$

$$\tilde{\epsilon }\equiv \epsilon -m{\mathbf{v}}_{\mathrm{F}}\cdot {\mathbf{v}}_{\mathrm{s}},$$

$${\mathbf{v}}_{\mathrm{s}}\equiv -\frac{\hslash }{2m}\nabla \phi .$$

$${\mathbf{v}}_{\mathrm{F}}\cdot \nabla g_1^{\mathrm{R},\mathrm{A}}$$

$$g_1^{\mathrm{R},\mathrm{A}}$$

$$f_1^{\mathrm{R},\mathrm{A}}$$

$${\left[\epsilon {\hat{\tau }}_3-{\hat{h}}^{\mathrm{R},\mathrm{A}},\delta {\hat{g}}^{\mathrm{R},\mathrm{A}}\right]}_{\circ }-{\left[\delta {\hat{h}}^{\mathrm{R},\mathrm{A}},{\hat{g}}^{\mathrm{R},\mathrm{A}}\right]}_{\circ }+i\hslash {\mathbf{v}}_{\mathrm{F}}\cdot \nabla \delta {\hat{g}}^{\mathrm{R},\mathrm{A}}$$

$$+\frac{i\hslash }{2}e\left({\mathbf{v}}_{\mathrm{F}}\times \mathbf{B}\right)\cdot \frac{\partial }{\partial {\mathbf{p}}_{\mathrm{F}}}\left\{{\hat{\tau }}_3,\delta {\hat{g}}^{\mathrm{R},\mathrm{A}}\right\}$$

$$+\frac{i\hslash }{2}e\left({\mathbf{v}}_{\mathrm{F}}\times \delta \mathbf{B}\right)\cdot \frac{\partial }{\partial {\mathbf{p}}_{\mathrm{F}}}\left\{{\hat{\tau }}_3,{\hat{g}}^{\mathrm{R},\mathrm{A}}\right\}=\hat{0},$$

$${\left[\epsilon {\hat{\tau }}_3,\delta {\hat{g}}^{\mathrm{K}}\right]}_{\circ }+i\hslash {\mathbf{v}}_{\mathrm{F}}\cdot \nabla {\hat{g}}^{\mathrm{K}}+i\hslash {\mathbf{v}}_{\mathrm{F}}\cdot \nabla \delta {\hat{g}}^{\mathrm{K}}$$

$$-{\hat{h}}^{\mathrm{R}}\circ \delta {\hat{g}}^{\mathrm{K}}-\delta {\hat{h}}^{\mathrm{R}}\circ {\hat{g}}^{\mathrm{K}}-{\hat{h}}^{\mathrm{K}}\circ \delta {\hat{g}}^{\mathrm{A}}-\delta {\hat{h}}^{\mathrm{K}}\circ {\hat{g}}^{\mathrm{A}}$$

Full text

I Introduction

It was observed experimentally that the thermal conductivity and the thermal Hall angle are greatly enhanced in the superconducting states of cuprates Cohn ; Krishana95 ; Zeini99 ; Krishana99 ; Zhang00 ; Ocana ; Zhang01 ; Zeini01 and CeCoIn5 Kasahara . These materials are regarded as $d$-wave superconductors, and their thermal Hall conductivity is also consistent with the scaling relation for a $d$-wave pairing based on the Bogoliubov–de Gennes equation Simon . Hirschfeld et al. also calculated the longitudinal thermal conductivity in a $d$-wave superconductor, including the effect of antiferromagnetic spin fluctuations within the random phase approximation (RPA), and studied the large peak structure of the longitudinal thermal conductivity in YBa2Cu3O7-δ (YBCO) Hirschfeld96 . The enhancement of the thermal Hall angle has been explained to originate from the enhancement of the quasiparticle mean free path Zhang01 ; Kasahara . However, to the best of our knowledge, its microscopic calculation has not been carried out yet. Therefore, the origin of the enhancement of the quasiparticle mean free path needs to be clarified using a fully microscopic treatment. The purpose of this present paper is to develop a theoretical formalism for investigating the thermal Hall conductivity microscopically within the quasiclassical theory and to study this drastic enhancement of the thermal Hall angle in the superconducting states of YBCO and CeCoIn5.

Numerous theoretical studies have been carried out, such as on the longitudinal component of the thermal conductivities in both conventional ($s$-wave) Bardeen ; Graf ; Adachi and unconventional superconductors Hirschfeld96 ; Schmitt ; Hirschfeld88 ; Graf ; Kubert ; Franz ; Vekhter ; Vorontsov ; Adachi ; Hara ; Kobayashi and the spontaneous thermal Hall conductivity due to the intrinsic Read ; Nomura ; Sumiyoshi ; Imai16 ; Yoshioka and extrinsic Yip ; Yilmaz ; Ngampruetikorn ; Ngampruetikorn20 mechanisms in chiral superconductors. On the other hand, several studies have been carried out on the thermal Hall conductivity in $d$-wave superconductors, in terms of the cross section of quasiparticle scattering from a single vortex Durst , due to the Berry phase acquired by a quasiparticle in the vortex lattice state Vafek01 ; Cvetkovic ; Murray ; Vafek15 , based on the Kubo formula approach and taking into account the Doppler shift effect Shahzamanian . The thermal Hall conductivity due to the Lorentz force in type-II superconductors is not fully understood, this may be because the Lorentz force is missing from the quasiclassical Eilenberger equations Eilenberger , which is a powerful tool for investigating inhomogeneous and nonequilibrium superconductors microscopically Eschrig . More precisely, the component of the magnetic Lorentz force balanced with the Hall electric field and/or force induced by a transverse temperature gradient may be missing from the standard Eilenberger equations, since there is the component of the magnetic Lorentz force balanced with the hydrodynamic force in the Ginzburg–Landau equations Kato16 .

Here we derive the thermal Hall coefficient in extreme type-II superconductors with an isolated pinned vortex, based on the augmented quasiclassical equations of superconductivity with the Lorentz force in the Keldysh formalism Kita01 by incorporating a next-to-leading-order contribution in the expansion of the Gor’kov equations Gor'kov ; Gor'kov2 in terms of the quasiclassical parameter $\delta\equiv 1/k_{\rm F}\xi_{0}(0)$, where $k_{\rm F}$ and $\xi_{0}(0)$ denote the Fermi wave number and the coherence length at zero temperature, respectively. The calculations of linear responses based on the augmented equations obtained from this derivation Kita01 have not been performed yet, as far as we know, except for the flux-flow Hall effect Arahata . In extreme type-II superconductors with an isolated pinned vortex, the thermal conductivity is dominated by the contribution from quasiparticles outside the core. Thus, we consider the contribution of the Doppler shifted quasiparticles due to the circulating supercurrent around the core (i.e., the Volovik effect) Kubert ; Volovik ; Kopnin and neglect that of the Andreev-reflected quasiparticles in the core Dahm ; Nagai . The newly derived expression for the thermal conductivity is an extension of the expression proposed by Graf et al. Graf and can also describe the thermal Hall effect in both fully gapped and nodal superconductors. Based on this formalism, we calculate the thermal conductivity and the thermal Hall angle in $d$-wave superconductors and address the enhancement of $\lim_{H\to 0}\tan\theta_{\rm H}/H$ at low temperatures in YBCO and CeCoIn5 measured by Zhang et al. in Ref. Zhang01, and by Kasahara et al. in Ref. Kasahara, , where $\tan\theta_{\rm H}$ and $H$ denote the thermal Hall angle and the external magnetic field, respectively. It is also imperative to consider impurity scattering close to the unitarity limit to explain the experimental results on the longitudinal thermal conductivity in heavy-fermion superconductors and Zn-doped YBCO Coffey ; Hirschfeld96 ; Schmitt ; Hirschfeld88 . Thus, we also study the impurity effect on the thermal Hall conductivity and the thermal Hall angle in a $d$-wave superconductor, using the impurity self-energy by $t$-matrix approximation Hirschfeld96 ; Schmitt ; Hirschfeld88 ; Graf ; Vorontsov ; Kato00 .

This paper is organized as follows. In Sec. II, we present the augmented quasiclassical equations of superconductivity with the Lorentz force in the Keldysh formalism. In Sec. III, we derive the thermal Hall coefficient in extreme type-II superconductors with an isolated pinned vortex based on the augmented quasiclassical equations of superconductivity with the Lorentz force and the linear response theory. In Sec. IV, we present numerical results for the thermal Hall effect in $d$-wave superconductors. In Sec. V, we provide concluding remarks.

II Augmented Quasiclassical Equations

We consider type-II superconductors with pinned vortices and neglect the pair-potential-gradient force Ohuchi ; Masaki17 ; Masaki18 and the pressure difference arising from the slope in the density of states (DOS) Ueki18 , which only contribute to the vortex-core charging in this case, since the charge in a pinned vortex does not contribute to thermal conductivity. For simplicity, we also restrict ourselves to the spin-singlet pairing without spin paramagnetism. Then the quasiclassical equations of superconductivity with the Lorentz force are given in the Keldysh formalism by Kita01

[TABLE]

where $e<0$ is the electron charge, ${\bm{E}}$ is the electric field, ${\bm{B}}$ is the magnetic field, $\varepsilon$ is the excitation energy, ${\bm{v}}_{\rm F}$ is the Fermi velocity, and ${\bm{p}}_{\rm F}$ is the Fermi momentum. The Green’s functions $\check{g}$, the pair potential $\check{\Delta}$, and the Born-type impurity self-energy $\check{\sigma}$ can be written as

[TABLE]

where $\tau$ is the relaxation time and $\langle\cdots\rangle_{\rm F}$ denotes the Fermi surface average with $\langle 1\rangle_{\rm F}=1$. This Born-type impurity self-energy will later be rewritten to the impurity self-energy by the $t$-matrix approximation Graf ; Vorontsov ; Kato00 . The $2\times 2$ retarded and Keldysh Green’s functions and the $2\times 2$ pair potential can be written as

[TABLE]

where the barred function in the Keldysh formalism is defined generally by $\bar{X}(\varepsilon,{\bm{p}}_{\rm F},{\bm{r}},t)\equiv X^{*}(-\varepsilon,-{\bm{p}}_{\rm F},{\bm{r}},t),$ $\Delta=\Delta({\bm{r}},t)$ denotes the amplitude of the energy gap, and $\phi=\phi({\bm{p}}_{\rm F})$ is the basis function on the Fermi surface normalized as $\langle|\phi|^{2}\rangle_{\rm F}=1$. Matrix $\check{\tau}_{3}$ is given by

[TABLE]

The Green’s functions satisfy the following symmetry relations:

[TABLE]

We choose the gauge ${\bm{E}}=-\partial{\bm{A}}/\partial t$ and ${\bm{B}}={\bm{\nabla}}\times{\bm{A}}$ with $\Phi=0$, where ${\bm{A}}$ and $\Phi$ denote the vector and scalar potentials. Notations $[a,b]_{\circ}$ and $\{a,b\}$ are then given by $[a,b]_{\circ}\equiv a\circ b-b\circ a$ and $\{a,b\}\equiv ab+ba$ with

[TABLE]

The gauge invariant derivative ${\bm{\partial}}$ is defined by

[TABLE]

The equation of the gap amplitude for the weak-coupling limit is given by

[TABLE]

where $\varepsilon_{\rm c}$ and ${\rm g}_{0}$ denote the cutoff energy and the coupling constant, respectively, defined by ${\rm g}_{0}\equiv-N(0)V_{l}^{({\rm eff})}$ with $V_{l}^{({\rm eff})}$ and $N(0)$ denoting the constant effective potential and the normal-state DOS per spin and unit volume at the Fermi level. Using the Green’s function, we can also express the heat flux density ${\bm{j}}_{\rm Q}$ as

[TABLE]

Now, we use the following relations:

[TABLE]

to rewrite Eq. (1) as

[TABLE]

Introducing matrix $\check{h}$ as Eschrig

[TABLE]

we obtain the equations for $\hat{g}^{\rm R,A,K}$ as

[TABLE]

III Thermal Hall Conductivity

We consider the linear response $\hat{X}=\hat{X}^{\rm le}+\delta\hat{X}$ to the thermal gradient ${\bm{\nabla}}T$, where $\hat{X}^{\rm le}$ is the solution in local equilibrium and $\delta\hat{X}$ is a term of the first order in ${\bm{\nabla}}T$ Vorontsov ; Graf . We also neglect the electric field ${\bm{E}}$, since the quasiparticle current and supercurrent counterflow in order to keep the charge current zero Ginzbrug ; Izawa07 , except the circulating supercurrent of a vortex, and the charge in a pinned vortex due to the Lorentz force Ueki16 does not contribute to thermal conductivity. Assuming extreme type-II superconductors as $\lambda_{0}\gg\xi_{0}(0)$ with an isolated pinned vortex, we furthermore use the Doppler shift method for the quasiclassical equations Dahm , where $\lambda_{0}$ denotes magnetic penetration depth defined by $\lambda_{0}\equiv[\mu_{0}N(0)e^{2}\langle v_{\rm F}^{2}\rangle_{\rm F}]^{-1/2}$, the coherence length $\xi_{0}(0)$ is defined by $\xi_{0}(0)\equiv\hbar\langle v_{\rm F}\rangle_{\rm F}/\Delta_{0}(0)$, $\Delta_{0}(0)$ denotes the gap amplitude in clean superconductors at zero temperature and zero magnetic field, and $\mu_{0}$ is the vacuum permeability. Then we can neglect the spatial derivatives of the phase-transformed Green’s functions, the vector potential terms. Hereafter, we remove the superscript “le” from these equations.

III.1 Local equilibrium

Equations for the retarded and advanced Green’s functions $\hat{g}^{\rm R,A}$ in local equilibrium are written as

[TABLE]

where notation $[a,b]$ is defined by $[a,b]\equiv ab-ba$. We first expand formally in the quasiclassical parameter $\delta\equiv\hbar/\langle p_{\rm F}\rangle_{\rm F}\xi_{0}(0)$ as $g^{\rm R,A}=g_{0}^{\rm R,A}+g_{1}^{\rm R,A}+\cdots$, $f^{\rm R,A}=f_{0}^{\rm R,A}+f_{1}^{\rm R,A}+\cdots$, and $\Delta=\Delta_{0}+\Delta_{1}+\cdots$ Kita09 . We next take the direction of the magnetic field as ${\bm{B}}=B\hat{\bm{z}}$, and transform $f_{0}^{\rm R,A}$ and $\Delta_{0}$ as $f_{0}^{\rm R,A}=\tilde{f}_{0}^{\rm R,A}{\rm e}^{-i\varphi}$ and $\Delta_{0}=\tilde{\Delta}_{0}{\rm e}^{-i\varphi}$ with the azimuth angle $\varphi$ in the cylindrical coordinate system $(\rho\cos\varphi,\rho\sin\varphi,z)$. Neglecting the spatial dependence of $\tilde{\Delta}_{0}$ and the imaginary part in $g_{0}^{\rm R,A}$ and $\tilde{f}_{0}^{\rm R,A}$, then we obtain

[TABLE]

$\varepsilon^{\rm R,A}$ and $\Delta^{\rm R,A}$ are defined by

[TABLE]

where $\eta$ denotes an infinitesimal positive constant and $\tilde{\varepsilon}$ is defined by

[TABLE]

with $m$ denoting the electron mass and ${\bm{v}}_{\rm s}$ denoting the superfluid velocity defined by

[TABLE]

Solving the equations for $g_{1}^{\rm R,A}$ and $f_{1}^{\rm R,A}$ in the almost same manner as Ref. Kita09, , we also obtain

[TABLE]

III.2 First-order response

Equations for the retarded and advanced Green’s functions $\delta\hat{g}^{\rm R,A}$ in the first-order response are written as

[TABLE]

and equations for the Keldysh Green’s functions $\delta\hat{g}^{\rm K}$ in the first-order response are given by

[TABLE]

Note that operator ${\bm{\nabla}}$ to $\hat{g}^{\rm K}$ and $\delta\hat{g}^{\rm K}$ only affects temperature in the distribution functions and the Green’s functions except the distribution functions, respectively Vorontsov . Let us introduce $\delta\hat{g}^{\rm a}$ as Graf ; Eschrig

[TABLE]

with

[TABLE]

The Green’s functions $\delta g^{\rm R,A}$ are used to calculate the supercurrent, and $\delta g^{\rm a}$ is related to the quasiparticle current. Thus, we derive equation for $\delta g^{\rm a}$ to obtain the heat current. Substituting Eq. (25) into Eq. (24), and using $\hat{g}^{\rm K}=(\hat{g}^{\rm R}-\hat{g}^{\rm A})\tanh\varepsilon/2k_{\rm B}T$ and Eq. (23), the equation for $\delta\hat{g}^{\rm a}$ can be expressed as

[TABLE]

We also neglect the time derivative terms and the self-energy correction terms as $\langle\delta g^{\rm a}\rangle_{\rm F}$ and $\langle\delta f^{\rm a}\rangle_{\rm F}$ Graf . Then the equations for $\delta g^{\rm a}$ and $\delta f^{\rm a}$ are obtained from the $(1,1)$ and $(1,2)$ components in Eq. (27) as

[TABLE]

We first calculate the zeroth-order quantity $\delta g_{0}^{\rm a}+\delta\bar{g}_{0}^{\rm a}$ in $\delta$, since the heat current [Eq. (12)] can be rewritten as

[TABLE]

Transforming the zeroth-order quantity $\delta f_{0}^{\rm a}$ in $\delta$ as $\delta f_{0}^{\rm a}=\delta\tilde{f}_{0}^{\rm a}{\rm e}^{-i\varphi}$, and neglecting the spatial derivatives of $\delta g_{0}^{\rm a}$, the equation for the zeroth-order quantity $\delta g_{0}^{\rm a}$ in $\delta$ is obtained from Eq. (28a) as

[TABLE]

We can obtain the equation for $\delta\bar{g}_{0}^{\rm a}$ from the barred Eq. (30a) and $\delta\bar{\tilde{f}}_{0}^{\rm a}$ from the barred Eq. (30b). We solve the equations for $\delta g_{0}^{\rm a}$, $\delta\bar{g}_{0}^{\rm a}$, $\delta\tilde{f}_{0}^{\rm a}$, and $\delta\bar{\tilde{f}}_{0}^{\rm a}$ choosing $\phi=\phi^{*}$ and using Eqs. (5) and (18). Then we obtain $\delta g^{\rm a}+\delta\bar{g}^{\rm a}$ as

[TABLE]

We next calculate the first-order quantity $\delta g_{1}^{\rm a}+\delta\bar{g}_{1}^{\rm a}$ in $\delta$. Transforming the first-order quantity $\delta f_{1}^{\rm a}$ in $\delta$ as $\delta f_{1}^{\rm a}=\delta\tilde{f}_{1}^{\rm a}{\rm e}^{-i\varphi}$, neglecting the spatial derivatives of $\delta g_{1}^{\rm a}$ and $\delta\tilde{f}_{1}^{\rm a}$, and using Eq. (22c), the equations for the first-order quantities $\delta g_{1}^{\rm a}$ and $\delta f_{1}^{\rm a}$ in $\delta$ are obtained from Eq. (28) as

[TABLE]

We solve the equations for $\delta g_{1}^{\rm a}$, $\delta\bar{g}_{1}^{\rm a}$, $\delta\tilde{f}_{1}^{\rm a}$, and $\delta\bar{\tilde{f}}_{1}^{\rm a}$ choosing $\phi=\phi^{*}$ and using Eqs. (5), (18), and (22b). Then we also obtain $\delta g_{1}^{\rm a}+\delta\bar{g}_{1}^{\rm a}$ as

[TABLE]

III.3 Coefficient of thermal conductivity

We finally calculate the thermal conductivity in extreme type-II superconductors with an isolated pinned vortex. The thermal conductivity $\underline{\kappa}$ is given by taking the spatial average of the local thermal conductivity $\underline{\kappa}^{\rm loc}({\bm{r}})$ defined by ${\bm{j}}_{\rm Q}=-\underline{\kappa}^{\rm loc}{\bm{\nabla}}T$. We thereby obtain the coefficient of thermal conductivity as

[TABLE]

where $\langle\cdots\rangle$ denotes the Fermi surface and spatial average. The first line in Eq. (34) for $B=0$ (i.e., ${\bm{v}}_{\rm s}={\bm{0}}$) is the same as the longitudinal thermal conductivity proposed in Ref. Graf, , and the second and third lines are the newly derived thermal Hall conductivity. Neglecting the Doppler shift effect as ${\bm{v}}_{\rm s}={\bm{0}}$, the expression for the thermal conductivity in $d$-wave superconductors can be rewritten as the following simpler form:

[TABLE]

We can use the impurity self-energy by the $t$-matrix approximation, changing $\tau$ in Eqs. (19) and (34) as Graf ; Vorontsov ; Kato00

[TABLE]

$\tau_{0}^{\rm R}$ is defined by

[TABLE]

where $\tau_{0}$ and $\delta_{0}$ are the relaxation time in the normal state and the scattering phase shift given by

[TABLE]

with $n_{a}$ and $U^{\rm imp}$ denoting the density of impurities and the impurity potential. Using $g_{0}^{\rm R}\to 1$ for the normal-state limit $\Delta\to 0$, we obtain the expression for the thermal conductivity in normal metal with a spherical Fermi surface as

[TABLE]

where $\sigma_{xx}^{(n)}=\tau_{0}e^{2}n/m$ denotes the dc conductivity in normal state obtained from the Drude model, $n=(2/3)mN(0)v_{\rm F}^{2}$ is the electron density.

IV Numerical Results

Here we present numerical results for the thermal Hall conductivity in $d$-wave superconductors with a cylindrical Fermi surface. Then $\langle\cdots\rangle$ is given by

[TABLE]

where $\varphi_{\bm{p}}$ denotes the angle of the ${\bm{p}}$ vector in the two-dimensional polar coordinate system as ${\bm{p}}=(p\cos\varphi_{\bm{p}},p\sin\varphi_{\bm{p}})$, and $R$ denotes the radius of the vortex lattice unit cell area given by $R\approx\xi_{0}(0)\sqrt{H_{\rm c2}/H}$ for $R\gg\xi_{0}(0)$ with $H_{c2}$ denoting the upper critical magnetic-field Kubert . We consider extreme type-II superconductors with an isolated vortex, and assume that the magnetic field $B$ is a very small constant within $r\leq R$.

To discuss the Volovik and impurity effects qualitatively, we first solve Eqs. (11) and (18) self-consistently to obtain $(g_{0}^{\rm R},f_{0}^{\rm R},\Delta_{0})$, and substituting the resulting solutions into Eq. (34) or (35), we can obtain the thermal conductivity and the thermal Hall angle. Hereafter, we choose the cutoff energy in the gap equation [Eq. (11)] as $\varepsilon_{\rm c}=40k_{\rm B}T_{\rm c}$, and fix the parameter to $\eta=0.00001\Delta_{0}(0)$, where $T_{\rm c}$ denotes the superconducting transition temperature without impurities and external fields. We also adopt a model $d$-wave pairing as $\phi=\sqrt{2}\cos 2\varphi_{\bm{p}}$, and the direction of applied thermal gradient is the same as the measurement in Ref. Zhang01, .

IV.1 Volovik effect

We first discuss the Doppler shift effect on the thermal conductivity and the thermal Hall angle in $d$-wave superconductors. It is shown that the Doppler shift effect can be neglected at very weak magnetic fields. In this subsection, we fix the parameters as $\sin^{2}\delta_{0}=0$ (the Born-type impurity) and $\tau_{0}=10\hbar/\Delta_{0}(0)$.

Figure 1 plots the longitudinal thermal conductivity $\kappa_{xx}$ normalized by the value at the transition temperature without impurities and external fields $T_{\rm c}$, for different vortex radii $R=7\xi_{0}(0),10\xi_{0}(0),\dots,\infty$ as a function of temperature. The thermal conductivity for $R\to\infty$ are the same as that for ${\bm{v}}_{\rm s}={\bm{0}}$. We observe that the normalized longitudinal thermal conductivity decreases as temperature decreases from $T=T_{\rm c}$, and $T$-linear behavior for $d$-wave superconductors at low temperatures. These behaviors at low temperatures are well-known results Coffey ; Graf .

Figure 2 plots the thermal Hall conductivity $\kappa_{xy}$ normalized by the value at the transition temperature without impurities and external fields $T_{\rm c}$, for the different vortex radii $R=7\xi_{0}(0),10\xi_{0}(0),\dots,\infty$ as a function of temperature. We observe that the normalized thermal Hall conductivity varies more as a function of the vortex radius compared with the longitudinal component. The normalized thermal Hall conductivity has a finite value even near $T=0$, and decreases significantly around $T=0.05T_{\rm c}$ as the vortex radius decreases from $R\to\infty$. We here note that the thermal Hall conductivity (without the normalization) is proportional to the magnetic field in the temperature range where the normalized thermal Hall conductivity does not change as a function of the vortex radius, since the thermal Hall conductivity at $T=T_{\rm c}$ used to normalize is also proportional to the magnetic field as seen in Eq. (39b). It means that the Doppler shift effect is very small in the temperature range where the normalized thermal Hall conductivity does not almost change.

Figure 3 plots the thermal Hall angle $\tan\theta_{\rm H}\equiv\kappa_{xy}/\kappa_{xx}$ normalized by the value at the transition temperature without impurities and external fields $T_{\rm c}$, for the different vortex radii $R=7\xi_{0}(0),10\xi_{0}(0),\dots,\infty$ as a function of temperature. We observe that the normalized thermal Hall angle greatly enhanced in $d$-wave superconductors.

We confirmed that the normalized thermal conductivity and the normalized thermal Hall angle approach that for ${\bm{v}}_{\rm s}={\bm{0}}$ as the vortex radius $R$ increases, and the normalized thermal conductivity and the normalized thermal Hall angle for ${\bm{v}}_{\rm s}={\bm{0}}$ are almost equal to that at $R=1000\xi_{0}(0)$ within $T\geq 0.01T_{\rm c}$.

IV.2 Impurity effect

We finally discuss the impurity effect on the thermal conductivity and the thermal Hall angle normalized by the value at the transition temperature without impurities and external fields $T_{\rm c}$ in $d$-wave superconductors for ${\bm{v}}_{\rm s}={\bm{0}}$, since the experimental results on the longitudinal thermal conductivity in heavy-fermion superconductors and Zn-doped YBCO were explained by taking into account impurity scattering close to the unitarity limit Coffey ; Hirschfeld96 ; Schmitt ; Hirschfeld88 , and $\tan\theta_{\rm H}/H$ and $\kappa_{xy}/H$ in the zero magnetic-field limit $H\to 0$ were estimated experimentally by the initial slope in the $H$ dependence of $\kappa_{xy}$ in Refs. Zhang01, and Kasahara, . The thermal Hall conductivity and the thermal Hall angle at $T=T_{\rm c}$ used to normalize are proportional to the “internal” magnetic field $B$ [see Eq. (39b)], and the internal magnetic field $B$ in extreme type-II superconductors also approaches zero as the external magnetic field $H$ approaches zero. Thus, the normalized thermal Hall conductivity and the normalized thermal Hall angle in extreme type-II superconductors for ${\bm{v}}_{\rm s}={\bm{0}}$ are the same as evaluating $\lim_{H\to 0}\tan\theta_{\rm H}/H$ and $\lim_{H\to 0}\kappa_{xy}/H$. Here we discuss the contribution from the scattering of quasiparticles on vortices Durst . The thermal Hall conductivity due to the vortex scattering at weak magnetic fields is given by $\kappa_{xy}^{\rm vortex}\sim T\sqrt{H}$. Since $\kappa_{xy}^{\rm vortex}/H$ diverges for $H\to 0$ regardless of temperature, and the longitudinal thermal conductivity is given by $\kappa_{xx}\sim T$, $\tan\theta_{\rm H}/H$, and $\kappa_{xy}/H$ for $H\to 0$ are dominated by the Lorentz force and the impurities. Therefore, we can study the impurity effect on thermally excited quasiparticles outside the core in extreme type-II superconductors, estimating $\lim_{H\to 0}\tan\theta_{\rm H}/H$ and $\lim_{H\to 0}\kappa_{xy}/H$ experimentally.

Figures 4, 5, 6, and 7 plot the longitudinal thermal conductivity $\kappa_{xx}$ and $\kappa_{xx}/T$, the thermal Hall conductivity $\kappa_{xy}$, and the thermal Hall angle $\tan\theta_{\rm H}$ normalized by the value at the transition temperature without impurities and external fields $T_{\rm c}$, for the different scattering phase shifts $\sin^{2}\delta_{0}=0,0.2,\dots,1$ at the normal relaxation time $\tau_{0}=10\hbar/\Delta_{0}(0)$ as a function of temperature. We reproduce the longitudinal thermal conductivity $\kappa_{xx}$ proposed by Graf et al. Graf , and show that the behavior of the normalized thermal Hall conductivity is similar to that of the normalized longitudinal thermal conductivity. We also observe that the enhancement of the normalized thermal Hall angle near zero temperature decreases as $\sin^{2}\delta_{0}$ increases from $\sin^{2}\delta_{0}=0$, and does not occur at $\sin^{2}\delta_{0}=1$. Figures 8 and 9 plot the thermal Hall angle normalized by the value at the transition temperature without impurities and external fields $T_{\rm c}$, for the different normal relaxation times $\tau_{0}=1.5\hbar/\Delta_{0}(0),2\hbar/\Delta_{0}(0),\dots,100\hbar/\Delta_{0}(0)$ in the Born and unitarity limit as $\sin^{2}\delta_{0}=0$ and $1$, respectively, as a function of temperature. It is shown that the normalized thermal Hall angle is greatly enhanced in the Born limit as $\sin^{2}\delta_{0}=0$ near zero temperature but is suppressed in the unitarity limit as $\sin^{2}\delta_{0}=1$ near zero temperature even at the large normal relaxation time $\tau_{0}=100\hbar/\Delta_{0}(0)$. We can explain the normal relaxation time and temperature dependence of the normalized thermal Hall angle in the Born and unitarity limit as follows, using the corresponding quasiparticle DOS in Figs. 10 and 11.

We first consider the Born-type impurity as $\sin^{2}\delta_{0}=0$. Then the quasiparticle relaxation time $\tau_{\rm QP}$ can be written from Eq. (35) as $\tau_{\rm QP}=\tau_{0}/\langle{\rm Re}g_{0}^{\rm R}\rangle_{\rm F}$ near zero temperature, since quasiparticles around the gap nodes at $\phi=0$ become dominant. Furthermore, in the $d$-wave pairing case, since the quasiparticle DOS has a finite small value near $\varepsilon=0$, the impurity scattering for quasiparticles near $\varepsilon=0$ is very small. We can also explain that the quasiparticle impurity scattering is very small near zero temperature, since quasiparticles around the gap nodes are restricted to the momentum in a specific orientation. In the Born limit, the normalized thermal Hall angle near zero temperature can be roughly written as

[TABLE]

where the $\varepsilon$ integration signifies an integration over $-\infty\leq\varepsilon\leq\infty$, except $\langle{\rm Re}g_{0}^{\rm R}\rangle_{\rm F}=0$. Thus, the normalized thermal Hall angle is more enhanced in the Born limit at the large normal relaxation time near zero temperature, since the quasiparticle DOS near $\varepsilon=0$ becomes small as the normal relaxation time increases, as seen in Fig. 10. On the other hand, when taking the unitarity limit as $\sin^{2}\delta_{0}=1$, the quasiparticle relaxation time $\tau_{\rm QP}$ near zero temperature can be given by $\tau_{\rm QP}=\tau_{0}\langle{\rm Re}g_{0}^{\rm R}\rangle_{\rm F}$, and the restriction on the direction of the quasiparticle momentum is relaxed due to the strong scattering effect. In the unitarity limit, the normalized thermal Hall angle near zero temperature can be also roughly written as

[TABLE]

Thus, the normalized thermal Hall angle in the the unitarity limit is more suppressed at the large normal-relaxation-time near zero temperature inversely to the Born-type impurity, since the quasiparticle DOS near $\varepsilon=0$ is smaller as $\tau_{0}$ increases, as seen Fig. 11. However, since the thermal conductivity in $d$-wave superconductors with scattering close to the unitarity limit is also almost zero at low temperature, it may be difficult to observe the suppression of the Hall angle experimentally.

In our results, the normalized thermal Hall angle near zero temperature was up to about $100$ times larger than that at the transition temperature. From the figures in Refs. Zhang01, and Kasahara, estimating roughly the ratio of the thermal Hall angle at zero temperature to that at the transition temperature, that in Ref. Zhang01, is about 100, and that in Ref. Kasahara, is about 15. Here we note that these thermal-Hall-angle ratios were estimated in the limit of zero magnetic field. The thermal-Hall-angle ratio in Ref. Zhang01, is approximately consistent with our result for the large normal relaxation time in the Born limit, and that in Ref. Kasahara, is almost the same as the result for $\tau_{0}=100\hbar/\Delta_{0}(0)$ and $\sin^{2}\delta_{0}=0.8$, or $\tau_{0}=10\hbar/\Delta_{0}(0)$ and $\sin^{2}\delta_{0}=0.6$, or …, or $\tau_{0}=2\hbar/\Delta_{0}(0)$ and $\sin^{2}\delta_{0}=0$. Thus, we find from calculations of the thermal Hall angle that YBCO in Ref. Zhang01, is very clean, and materials in Ref. Kasahara, is relatively dirty within this our present theory.

V Concluding remarks

We derived the thermal conductivity in extreme type-II superconductors with an isolated pinned vortex based on the augmented quasiclassical equations of superconductivity with the Lorentz force. Using it, we calculated the thermal conductivity and the thermal Hall angle in $d$-wave superconductors, and addressed the enhancement of $\lim_{H\to 0}\tan\theta_{\rm H}/H$ at low temperatures in YBCO and CeCoIn5 measured by Zhang et al. Zhang01 and by Kasahara et al. Kasahara . They estimated $\lim_{H\to 0}\tan\theta_{\rm H}/H$ using the initial slope in the $H$ dependence of $\kappa_{xy}$ measured experimentally. We also should point out that since the superfluid velocity ${\bm{v}}_{\rm s}$ is proportional to $1/\rho$, the Doppler shift effect has the potential to contribute to the spatial average of the thermal conductivity and the thermal Hall angle. Thus, we calculated the magnetic-field dependence of $\tan\theta_{\rm H}$ and $\kappa_{xy}$ due to the Doppler shift effect, and confirmed that the Doppler shift effect can be neglected in the limit $H\to 0$.

The thermal Hall conductivity due to the vortex scattering at weak magnetic fields proposed by Durst et al. in Ref. Durst, is given by $\kappa_{xy}^{\rm vortex}=C_{0}T\sqrt{H}$, where $C_{0}$ is a constant introduced in Ref. Durst, , and it has been used successfully to explain the experimental results on the magnetic field dependence of the thermal conductivity in YBa2Cu3O6.99 in the weak magnetic-field region. However, they did not address the enhancement of $\lim_{H\to 0}\tan\theta_{\rm H}/H$ at low temperatures even in Refs. Ganeshan, and Kulkarni, , which are their more recent papers, and we cannot explain this enhancement using their theories. Unlike our theory, since their theory takes into account only the contribution of quasiparticles around the gap nodes, it is impossible to investigate the wide temperature dependence from the transition temperature to zero temperature. Moreover, $\kappa_{xy}^{\rm vortex}$ is proportional to $\sqrt{H}$, and $\kappa_{xy}^{\rm vortex}/H$ diverges for $H\to 0$, i.e., the quasiparticle relaxation time due to the vortex scattering diverges for $H\to 0$ regardless of temperature, since the longitudinal thermal conductivity becomes $\kappa_{xx}\sim T$ for $H\to 0$. On the other hand, the thermal Hall conductivity due to the Lorentz force $\kappa_{xy}^{\rm Lorentz}$ is proportional to $H$ [see Eq. (34)], and $\kappa_{xy}^{\rm Lorentz}/H$ has a finite value for $H\to 0$. Therefore, $\lim_{H\to 0}\tan\theta_{\rm H}/H$ and $\lim_{H\to 0}\kappa_{xy}/H$ are dominated by the Lorentz force and the impurities. We observed that the thermal Hall angle is greatly enhanced in $d$-wave superconductors without impurities of the resonant scattering at low temperatures and very weak magnetic field. This great enhancement of the thermal Hall angle has been observed experimentally in YBCO Zhang01 and CeCoIn5 Kasahara , and may also be observed experimentally in other nodal superconductors with large magnetic-penetration-depth.

On the other hand, our absolute values of the thermal conductivity are different from the experimental values. To discuss the experimental values more quantitatively, we may need to consider inelastic scattering as considered in the calculation of the longitudinal component, such as scattering by the antiferromagnetic spin fluctuations within the RPA Hirschfeld96 or the fluctuation exchange approximation Hara in YBCO. The relaxation times for inelastic scattering in the superconducting phase of UPt3 and CeCoIn5 are also given in Refs. Fledderjohann, and Graf96, , and Truncik, . Here we emphasize that the great enhancement of the electrical and longitudinal thermal conductivity in the superconducting state of YBCO is dominated by the contribution from inelastic scattering Hirschfeld96 , but thermal Hall angle is greatly enhanced in nodal superconductors even without inelastic scattering.

We also discuss the contribution from quasiparticles in the vortex core. For clean superconductors in the vortex lattice state under larger magnetic field, we should consider the contribution of quasiparticles in the vortex core. The thermal Hall angle in the vortex lattice state may decrease due to the increase in quasiparticle scattering in the vortex core as the external magnetic field increases Kasahara , but we also need to calculate it microscopically and quantitatively to clarify the suppression of the thermal Hall angle due to external magnetic field. We may confirm this by solving the augmented quasiclassical equations directly in the vortex lattice system Adachi ; Ichioka99 (see Ref. Ueki18, for the pair-potential-gradient force and the pressure difference arising from the slope in the DOS). This method can include not only the contribution of quasiparticles outside the core but also that of the Andreev-reflected quasiparticles Dahm ; Nagai in the thermal conductivity. We may also calculate the thermal Hall conductivity in high magnetic fields more easily, using the Brandt–Pesch–Tewordt approximation for the augmented quasiclassical equations Houghton . Moreover, we may also need to confirm the contributions from skew scattering due to a vortex Durst ; Kulkarni and the Berry phase acquired by a quasiparticle on traversing around a vortex Vafek01 ; Cvetkovic ; Murray ; Vafek15 ; Ganeshan in the calculations under larger magnetic field.

Our approach for the study of thermal Hall conductivity in type-II superconductors can be applied to the microscopic study of quasiparticle transport in a variety of superconductors under magnetic field.

Acknowledgements.

H. U. thanks K. Izawa and J. Goryo for extensive discussions on the thermal conductivity in superconductors, Y. Masaki and W. Kohno for sound advice about the impurity and Volovik effects, T. Matsushita and A. Kirikoshi for helpful discussions on the formulation in this work, and E. S. Joshua for many discussions and comments on our present paper. This work was partially supported by JSPS KAKENHI Grant No. 15H05885 (J-Physics). The computation in this work was partially carried out using the facilities of the Supercomputer Center, the Institute for Solid State Physics, the University of Tokyo.