Possible superconductivity induced by a large spin-orbit coupling in carrier doped iridium oxide insulators: A weak coupling approach

TL;DR

This study investigates the potential for superconductivity in doped iridium oxide Sr$_{2}$IrO$_{4}$, revealing that large spin-orbit coupling favors $d_{x^{2}-y^{2}}$-wave singlet pairing in electron-doped cases and $s_{ ext{±}}$-wave pairing in hole-doped cases, with implications for strongly correlated materials.

Contribution

It provides a systematic analysis of pairing symmetries in doped iridium oxides considering strong spin-orbit coupling using a three-orbital Hubbard model and Eliashberg equations.

Findings

Large SOC favors $d_{x^{2}-y^{2}}$-wave pairing in electron doping.

Small SOC and large Hund's coupling favor $s_{ ext{±}}$-wave pairing in hole doping.

Pseudospin triplet pairing is unfavorable across studied parameters.

Abstract

We study possible superconductivity in a carrier-doped iridium oxide insulator Sr$_{2}$IrO$_{4}$ based on an effective $t_{2g}$ three-orbital Hubbard model on the square lattice with a large spin-orbit coupling (SOC). Numerically solving the linearized Eliashberg equation for the superconducting (SC) gap function with the random phase approximation, we systematically examine both singlet and triplet SC gap functions with possible pairing symmetries and the parameter dependence of the superconductivity. For the realistic SOC $\lambda$ and Hund's coupling $J/U$ relevant to Sr$_{2}$IrO$_{4}$, namely, for a large $\lambda$ and small $J/U$ region, we find that the intra-band antiferromagnetic (AF) pseudospin $\bm{j}_{\text{eff}} = -\bm{l}+\bm{s}$ fluctuations favor a $d_{x^{2}-y^{2}}$-wave pseudospin $j_{\text{eff}}=1/2$ singlet pairing in the electron-doping. We also find that the $d_{x^{2}-y^{2}}$-wave pairing is more stabilized with increasing the SOC and decreasing the Hund's coupling. Furthermore, we show for a small $\lambda$ and large $J/U$ region that an $s_{\pm}$-wave singlet pairing is favored in the hole-doped region. The origin of the $s_{\pm}$-wave pairing is due to the inter-band pair scattering arising from the intra-orbital AF spin $\bm s$ fluctuations. Although the possibility of a pseudospin triplet pairing is considered, we find it always unfavorable for all parameters studied here. The experimental consequences for other strongly correlated materials with a large SOC are also discussed.

Tables

Table 1. Table 1: Comparison between 3 ​ d 3 𝑑 3d transition metal oxide La 2 CuO 4 and 5 d 𝑑 d transition metal oxide Sr 2 IrO 4 .

	La2CuO4 [5, 4, 6, 7]	Sr2IrO4[54, 55, 9, 10, 13, 14, 53]
Crystal structure	K2NiF4 type	K2NiF4 type
Number of holes	one hole per $d_{x^2-y^2}$	one hole per $j_{\mathrm{eff}}=1/2$
Ground state	$s=1/2$ AF insulator	$j_{\mathrm{eff}}=1/2$ AF insulator
Neel temperature	$T_N=325$ K	$T_N=230$ K
Magnetic exchange	$J_{\mathrm{ex}}\approx 125$ meV	$J_{\mathrm{ex}}\approx 50-100$ meV
Carrier doping	high-$T_{\text{c}}$ superconductivity	??

Table 2. Table 2: Dominant components (left column) of the SC gap function Δ σ ​ σ ′ a ​ b ​ ( 𝒌 ) subscript superscript Δ 𝑎 𝑏 𝜎 superscript 𝜎 ′ 𝒌 \Delta^{ab}_{\sigma\sigma^{\prime}}(\bm{k}) in the t 2 ​ g subscript 𝑡 2 𝑔 t_{2g} basis found in the numerically obtained pseudospin singlet SC gap function D η ​ η ′ m ​ m ′ ​ ( 𝒌 ) subscript superscript 𝐷 𝑚 superscript 𝑚 ′ 𝜂 superscript 𝜂 ′ 𝒌 D^{mm^{\prime}}_{\eta\eta^{\prime}}(\bm{k}) and the sign change of Δ σ ​ σ ′ a ​ b ​ ( 𝒌 ) subscript superscript Δ 𝑎 𝑏 𝜎 superscript 𝜎 ′ 𝒌 \Delta^{ab}_{\sigma\sigma^{\prime}}(\bm{k}) under the parity transformation, orbital index exchange, and spin index exchange (three columns in the middle). The sign change of Δ σ ​ σ ′ a ​ b ​ ( 𝒌 ) subscript superscript Δ 𝑎 𝑏 𝜎 superscript 𝜎 ′ 𝒌 \Delta^{ab}_{\sigma\sigma^{\prime}}(\bm{k}) when it is transformed to Δ σ ′ ​ σ b ​ a ​ ( − 𝒌 ) subscript superscript Δ 𝑏 𝑎 superscript 𝜎 ′ 𝜎 𝒌 \Delta^{ba}_{\sigma^{\prime}\sigma}(-\bm{k}) is also shown in the right column.

Pseudospin singlet: $D_{\eta {\eta }^{\prime }}^{m{m}^{\prime }}\left(𝒌\right)=-D_{{\eta }^{\prime }\eta }^{m^{\prime }m}\left(-𝒌\right)$
${\Delta }_{\sigma {\sigma }^{\prime }}^{ab}\left(𝒌\right)$	Parity ($𝒌\leftrightarrow -𝒌$)	Orbital ($a\leftrightarrow b$)	Spin ($\sigma \leftrightarrow {\sigma }^{\prime }$)	${\Delta }_{{\sigma }^{\prime }\sigma }^{ba}\left(-𝒌\right)$
${\Delta }_{\sigma \overline{\sigma }}^{xyxy}\left(𝒌\right)$	+	+	-	-
${\Delta }_{\sigma \overline{\sigma }}^{yzyz}\left(𝒌\right)$	+	+	-	-
${\Delta }_{\sigma \overline{\sigma }}^{zxzx}\left(𝒌\right)$	+	+	-	-
${\Delta }_{\sigma \sigma }^{yzxy}\left(𝒌\right)$	+	-	+	-
${\Delta }_{\sigma \sigma }^{zxxy}\left(𝒌\right)$	+	-	+	-
${\Delta }_{\sigma \overline{\sigma }}^{yzzx}\left(𝒌\right)$	+	-	+	-

Equations

$$H=H_{\text{kin}}+H_{\text{SO}}+H_{\text{int}}.$$

$$H_{\text{kin}}=\sum _{\mathbf{k}}\sum _a\sum _{\sigma }({ϵ}_{\mathbf{k}}^a-\mu )c_{\mathbf{k}\sigma }^{a†}{c}_{\mathbf{k}\sigma }^a,$$

$$c_{\mathbf{k}\sigma }^{a†}=\frac{1}{\sqrt{N}}\sum _i{\mathrm{e}}^{\mathrm{i}\mathbf{k}\cdot {\mathbf{R}}_i}{c}_{i\sigma }^{a†}$$

$$c_{\mathbf{k}\sigma }^a=\frac{1}{\sqrt{N}}\sum _i{\mathrm{e}}^{-\mathrm{i}\mathbf{k}\cdot {\mathbf{R}}_i}{c}_{i\sigma }^a$$

$$\begin{split}H_{\text{SO}}&=\frac{\lambda}{2}\sum_{i}\sum_{a,b}\sum_{\sigma,\sigma^{\prime}}\langle a|\bm{l}|b\rangle\cdot\langle\sigma|\bm{s}|\sigma^{\prime}\rangle c^{a\dagger}_{i\sigma}c^{b}_{i\sigma^{\prime}}\\ &=\frac{\lambda}{2}\sum_{i}\sum_{\sigma}\left(\begin{array}[]{ccc}{c_{i\sigma}^{yz}}^{\dagger}\ \,{c_{i\sigma}^{zx}}^{\dagger}\ \,{c_{i\bar{\sigma}}^{xy}}^{\dagger}\end{array}\right)\\ &\quad\times\left(\begin{array}[]{ccc}0&\mathrm{i}s_{\sigma}&-s_{\sigma}\\ -\mathrm{i}s_{\sigma}&0&\mathrm{i}\\ -s_{\sigma}&-\mathrm{i}&0\end{array}\right)\left(\begin{array}[]{ccc}c_{i\sigma}^{yz}\\ c_{i\sigma}^{zx}\\ c_{i\bar{\sigma}}^{xy}\end{array}\right),\end{split}$$

$$\begin{array}{rl}{H}_{\text{int}}& =U\sum _i\sum _a{n}_{i\uparrow }^a{n}_{i\downarrow }^a+\frac{U^{\prime }}{2}\sum _i\sum _{a\ne b}{n}_i^a{n}_i^b\\ & -J_{\text{H}}\sum _i\sum _{a\ne b}\left({\mathbf{S}}_i^a\cdot {\mathbf{S}}_i^b+\frac{1}{4}{n}_i^a{n}_i^b\right)\\ & +\frac{J^{\prime }}{2}\sum _i\sum _{a\ne b}{c}_{i\uparrow }^{a†}{c}_{i\downarrow }^{a†}{c}_{i\downarrow }^b{c}_{i\uparrow }^b,\end{array}$$

$$\begin{split}H_{0}&=\sum_{\bm{k},\sigma}\left(\begin{array}[]{ccc}c^{yz\dagger}_{\bm{k}\sigma}\ \,c^{zx\dagger}_{\bm{k}\sigma}\ \,c^{xy\dagger}_{\bm{k}\bar{\sigma}}\end{array}\right)\\ &\quad\times\left(\begin{array}[]{ccc}\epsilon^{yz}_{\bm{k}}-\mu&\mathrm{i}s_{\sigma}\lambda/2&-s_{\sigma}\lambda/2\\ -\mathrm{i}s_{\sigma}\lambda/2&\epsilon^{zx}_{\bm{k}}-\mu&\mathrm{i}\lambda/2\\ -s_{\sigma}\lambda/2&-\mathrm{i}\lambda/2&\epsilon^{xy}_{\bm{k}}-\mu\end{array}\right)\left(\begin{array}[]{ccc}c^{yz}_{\bm{k}\sigma}\\ c^{zx}_{\bm{k}\sigma}\\ c^{xy}_{\bm{k}\bar{\sigma}}\end{array}\right).\end{split}$$

$$\begin{array}{rl}{ϵ}_{\mathbf{k}}^{yz}& =-2t_5\cos k_x-2t_4\cos k_y,\\ {ϵ}_{\mathbf{k}}^{zx}& =-2t_4\cos k_x-2t_5\cos k_y,\end{array}$$

$$\begin{array}{rl}{ϵ}_{\mathbf{k}}^{xy}& =-2t_1(\cos k_x+\cos k_y)-4t_2\cos k_x\cos k_y\\ & -2t_3(\cos 2k_x+\cos 2k_y)+{\mu }_{xy}\end{array}$$

$$\begin{array}{rl}& (t_1,t_2,t_3,t_4,t_5,{\mu }_{xy},\lambda )\\ & =(0.36,0.18,0.09,0.37,0.06,-0.36,0.37)\mathrm{eV}.\end{array}$$

$$\begin{split}&\left(\begin{array}[]{ccc}d^{1\dagger}_{\bm{k}\Uparrow}\ \,d^{2\dagger}_{\bm{k}\Uparrow}\ \,d^{3\dagger}_{\bm{k}\Uparrow}\ \,\big{|}\ \,d^{1\dagger}_{\bm{k}\Downarrow}\ \,d^{2\dagger}_{\bm{k}\Downarrow}\ \,d^{3\dagger}_{\bm{k}\Downarrow}\end{array}\right)\\ &=\left(\begin{array}[]{ccc}c^{yz\dagger}_{\bm{k}\uparrow}\ \,c^{zx\dagger}_{\bm{k}\uparrow}\ \,c^{xy\dagger}_{\bm{k}\uparrow}\ \,\big{|}\ \,c^{yz\dagger}_{\bm{k}\downarrow}\ \,c^{zx\dagger}_{\bm{k}\downarrow}\ \,c^{xy\dagger}_{\bm{k}\downarrow}\end{array}\right)\mathcal{U}^{\dagger}_{\bm{k}},\end{split}$$

$$\begin{split}&\left(\begin{array}[]{ccc}a^{1\dagger}_{i\Uparrow}\ \,a^{2\dagger}_{i\Uparrow}\ \,a^{3\dagger}_{i\Uparrow}\ \,\big{|}\ \,a^{1\dagger}_{i\Downarrow}\ \,a^{2\dagger}_{i\Downarrow}\ \,a^{3\dagger}_{i\Downarrow}\end{array}\right)\\ &=\left(\begin{array}[]{ccc}c^{yz\dagger}_{i\uparrow}\ \,c^{zx\dagger}_{i\uparrow}\ \,c^{xy\dagger}_{i\uparrow}\ \,\big{|}\ \,c^{yz\dagger}_{i\downarrow}\ \,c^{zx\dagger}_{i\downarrow}\ \,c^{xy\dagger}_{i\downarrow}\end{array}\right)\mathcal{W}^{\dagger},\end{split}$$

$$\mathcal{W}^{\dagger}=\left(\begin{array}[]{ccc|ccc}0&0&0&\frac{1}{\sqrt{3}}&\frac{1}{\sqrt{2}}&-\frac{1}{\sqrt{6}}\\ 0&0&0&-\frac{\mathrm{i}}{\sqrt{3}}&\frac{\mathrm{i}}{\sqrt{2}}&\frac{i}{\sqrt{6}}\\ \frac{1}{\sqrt{3}}&0&-\frac{2}{\sqrt{6}}&0&0&0\\ \hline\cr\frac{1}{\sqrt{3}}&-\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{6}}&0&0&0\\ \frac{\mathrm{i}}{\sqrt{3}}&\frac{\mathrm{i}}{\sqrt{2}}&\frac{\mathrm{i}}{\sqrt{6}}&0&0&0\\ 0&0&0&-\frac{1}{\sqrt{3}}&0&-\frac{2}{\sqrt{6}}\end{array}\right).$$

$$\begin{array}{rl}{G}_{\alpha \beta }(k)& =-{\int }_0^{\beta }d\tau {\mathrm{e}}^{\mathrm{i}{\omega }_n\tau }⟨T_{\tau }{c}_{\mathbf{k}\alpha }(\tau )c_{\mathbf{k}\beta }^{†}(0)⟩,\\ {\bar{G}}_{\alpha \beta }(k)& =-{\int }_0^{\beta }d\tau {\mathrm{e}}^{\mathrm{i}{\omega }_n\tau }⟨T_{\tau }{c}_{-\mathbf{k}\alpha }^{†}(\tau )c_{-\mathbf{k}\beta }(0)⟩,\end{array}$$

$$\hat{G}_{0}(k)=\Big{[}\mathrm{i}\omega_{n}-\hat{\xi}_{\bm{k}}\Big{]}^{-1},$$

$$\begin{split}H_{0}&=\sum_{\bm{k}}\sum_{\alpha,\beta}c^{\dagger}_{\bm{k}\alpha}\xi_{\bm{k}\,\alpha\beta}c_{\bm{k}\beta}\\ &=\sum_{\bm{k}}\left(\begin{array}[]{ccc}c^{yz\dagger}_{\bm{k}\uparrow}\ \,c^{zx\dagger}_{\bm{k}\uparrow}\ \,c^{xy\dagger}_{\bm{k}\uparrow}\ \,\big{|}\ \,c^{yz\dagger}_{\bm{k}\downarrow}\ \,c^{zx\dagger}_{\bm{k}\downarrow}\ \,c^{xy\dagger}_{\bm{k}\downarrow}\end{array}\right)\\ &\times\left(\begin{array}[]{ccc|ccc}\xi^{yz}_{\bm{k}}&\mathrm{i}\frac{\lambda}{2}&0&0&0&-\frac{\lambda}{2}\\ -\mathrm{i}\frac{\lambda}{2}&\xi^{zx}_{\bm{k}}&0&0&0&\mathrm{i}\frac{\lambda}{2}\\ 0&0&\xi_{\bm{k}}^{xy}&\frac{\lambda}{2}&-\mathrm{i}\frac{\lambda}{2}&0\\ \hline\cr 0&0&\frac{\lambda}{2}&\xi_{\bm{k}}^{yz}&-\mathrm{i}\frac{\lambda}{2}&0\\ 0&0&\mathrm{i}\frac{\lambda}{2}&{\rm i}\frac{\lambda}{2}&\xi_{\bm{k}}^{zx}&0\\ -\frac{\lambda}{2}&-\mathrm{i}\frac{\lambda}{2}&0&0&0&\xi_{\bm{k}}^{xy}\end{array}\right)\left(\begin{array}[]{c}c^{yz}_{\bm{k}\uparrow}\\ c^{zx}_{\bm{k}\uparrow}\\ c^{xy}_{\bm{k}\uparrow}\\ \hline\cr c^{yz}_{\bm{k}\downarrow}\\ c^{zx}_{\bm{k}\downarrow}\\ c^{xy}_{\bm{k}\downarrow}\end{array}\right).\end{split}$$

$$\begin{array}{rl}{F}_{\alpha \beta }(k)& =-{\int }_0^{\beta }d\tau {\mathrm{e}}^{\mathrm{i}{\omega }_n\tau }⟨T_{\tau }{c}_{\mathbf{k}\alpha }(\tau )c_{-\mathbf{k}\beta }(0)⟩,\\ {\bar{F}}_{\alpha \beta }(k)& =-{\int }_0^{\beta }d\tau {\mathrm{e}}^{\mathrm{i}{\omega }_n\tau }⟨T_{\tau }{c}_{-\mathbf{k}\alpha }^{†}(\tau )c_{\mathbf{k}\beta }^{†}(0)⟩.\end{array}$$

$$F_{\alpha \beta }(k)=\sum _{{\alpha }^{\prime },{\beta }^{\prime }}{G}_{\alpha {\alpha }^{\prime }}(k){\Delta }_{{\alpha }^{\prime }{\beta }^{\prime }}(k){\bar{G}}_{{\beta }^{\prime }\beta }(k),$$

$${\Delta }_{\alpha \beta }(k)=\frac{1}{N\beta }\sum _{k^{\prime }}\sum _{{\alpha }^{\prime }{\beta }^{\prime }}{V}_{\text{P}\alpha {\alpha }^{\prime }{\beta }^{\prime }\beta }(k-k^{\prime })F_{{\alpha }^{\prime }{\beta }^{\prime }}(k^{\prime }),$$

$$\begin{array}{rl}{\lambda }_{\text{e}}{\Delta }_{\alpha \beta }(\mathbf{k})& =-\frac{1}{N}\sum _{{\mathbf{k}}^{\prime }}\sum _{{\alpha }^{\prime }{\beta }^{\prime }}\sum _{\mu \nu }{V}_{\text{P}\alpha {\alpha }^{\prime }{\beta }^{\prime }\beta }(\mathbf{k}-{\mathbf{k}}^{\prime })\\ & \times {\varphi }_{{\alpha }^{\prime }{\beta }^{\prime }\mu \nu }({\mathbf{k}}^{\prime }){\Delta }_{\mu \nu }({\mathbf{k}}^{\prime }),\end{array}$$

$${\varphi }_{{\alpha }^{\prime }{\beta }^{\prime }\mu \nu }({\mathbf{k}}^{\prime })=-\frac{1}{\beta }\sum _{\mathrm{i}{\omega }_n}{G}_{0{\alpha }^{\prime }\mu }({\mathbf{k}}^{\prime },\mathrm{i}{\omega }_n){\bar{G}}_{0\nu {\beta }^{\prime }}({\mathbf{k}}^{\prime },\mathrm{i}{\omega }_n).$$

$${\Delta }_{\alpha \beta }^{ab}(\mathbf{k})=-{\Delta }_{\beta \alpha }^{ba}(-\mathbf{k})$$

$$\begin{array}{rl}{\chi }_{\alpha \beta \gamma \delta }(q)& =\frac{1}{N}{\int }_0^{\beta }d\tau {\mathrm{e}}^{\mathrm{i}{ϵ}_m\tau }\\ & \times \sum _{\mathbf{k},{\mathbf{k}}^{\prime }}⟨T_{\tau }{c}_{\mathbf{k}\beta }^{†}(\tau )c_{\mathbf{k}+\mathbf{q}\alpha }(\tau )c_{{\mathbf{k}}^{\prime }\gamma }^{†}(0)c_{{\mathbf{k}}^{\prime }-\mathbf{q}\delta }(0){⟩}_{\text{c}},\end{array}$$

$${\chi }_{0\alpha \beta \gamma \delta }(q)=-\frac{1}{N\beta }\sum _k{G}_{0\alpha \gamma }(q+k)G_{0\delta \beta }(k),$$

$${\hat{\chi }}_0(q)=\left(\begin{array}{cccc}{\hat{\chi }}_{0\uparrow \uparrow \uparrow \uparrow }& {\hat{\chi }}_{0\uparrow \uparrow \uparrow \downarrow }& {\hat{\chi }}_{0\uparrow \uparrow \downarrow \uparrow }& {\hat{\chi }}_{0\uparrow \uparrow \downarrow \downarrow }\\ {\hat{\chi }}_{0\uparrow \downarrow \uparrow \uparrow }& {\hat{\chi }}_{0\uparrow \downarrow \uparrow \downarrow }& {\hat{\chi }}_{0\uparrow \downarrow \downarrow \uparrow }& {\hat{\chi }}_{0\uparrow \downarrow \downarrow \downarrow }\\ {\hat{\chi }}_{0\downarrow \uparrow \uparrow \uparrow }& {\hat{\chi }}_{0\downarrow \uparrow \uparrow \downarrow }& {\hat{\chi }}_{0\downarrow \uparrow \downarrow \uparrow }& {\hat{\chi }}_{0\downarrow \uparrow \downarrow \downarrow }\\ {\hat{\chi }}_{0\downarrow \downarrow \uparrow \uparrow }& {\hat{\chi }}_{0\downarrow \downarrow \uparrow \downarrow }& {\hat{\chi }}_{0\downarrow \downarrow \downarrow \uparrow }& {\hat{\chi }}_{0\downarrow \downarrow \downarrow \downarrow }\end{array}\right)(q),$$

$$\hat{V}=\left(\begin{array}{cccc}{\hat{V}}_{\uparrow \uparrow \uparrow \uparrow }& 0& 0& {\hat{V}}_{\uparrow \uparrow \downarrow \downarrow }\\ 0& {\hat{V}}_{\uparrow \downarrow \uparrow \downarrow }& 0& 0\\ 0& 0& {\hat{V}}_{\downarrow \uparrow \downarrow \uparrow }& 0\\ {\hat{V}}_{\downarrow \downarrow \uparrow \uparrow }& 0& 0& {\hat{V}}_{\downarrow \downarrow \downarrow \downarrow }\end{array}\right),$$

$$\begin{array}{rl}& V_{\uparrow \uparrow \uparrow \uparrow }^{abcd}=V_{\downarrow \downarrow \downarrow \downarrow }^{abcd}=\{\begin{array}{ll}{U}^{\prime }-J& (a=c\ne b=d)\\ -(U^{\prime }-J)& (a=b\ne c=d)\\ 0& (\text{otherwise})\end{array},\\ & V_{\uparrow \uparrow \downarrow \downarrow }^{abcd}=V_{\downarrow \downarrow \uparrow \uparrow }^{abcd}=\{\begin{array}{ll}-U& (a=b=c=d)\\ -J& (a=c\ne b=d)\\ -U^{\prime }& (a=b\ne c=d)\\ -J& (a=d\ne b=c)\\ 0& (\text{otherwise})\end{array},\\ & V_{\uparrow \downarrow \uparrow \downarrow }^{abcd}=V_{\downarrow \uparrow \downarrow \uparrow }^{abcd}=\{\begin{array}{ll}U& (a=b=c=d)\\ U^{\prime }& (a=c\ne b=d)\\ J& (a=b\ne c=d)\\ J& (a=d\ne b=c)\\ 0& (\text{otherwise})\end{array},\end{array}$$

$$[\hat{A}\hat{B}{]}_{{\sigma }_1{\sigma }_2{\sigma }_3{\sigma }_4}^{abcd}=\sum _{m,n}\sum _{\sigma ,{\sigma }^{\prime }}{A}_{{\sigma }_1{\sigma }_2\sigma {\sigma }^{\prime }}^{abmn}{B}_{\sigma {\sigma }^{\prime }{\sigma }_3{\sigma }_4}^{mncd}.$$

$$\hat{\chi }(q)={\left(1-{\hat{\chi }}_0(q)\hat{V}\right)}^{-1}{\hat{\chi }}_0(q).$$

$$\hat{V}_{\text{P}}(q)=-\Bigg{[}\frac{\hat{V}}{2}+\hat{V}\hat{\chi}(q)\hat{V}\Bigg{]},$$

Full text

1 Introduction

Superconductivity is one of the most fundamental and intriguing quantum mechanical phenomena in condensed matter physics, and the discovery of new exotic superconductors such as high-$T_{\text{c}}$ cuprate superconductors, [1] organic superconductors, [2] and iron-based superconductors [3] has stimulated extensive experimental as well as theoretical studies. [4, 5, 6] The origins of these exotic superconductivity have been still under active debate. However, the current consensus is that the high-$T_{\text{c}}$ cuprate superconductor is an example of strongly correlated systems, and the antiferromagnetic (AF) fluctuations favor the unconventional $d$-wave superconducting (SC) state once mobile carriers are introduced into the AF Mott insulating phase. The highest critical temperature ($T_{\text{c}}$) at ambient pressure to date is achieved by hole doping in high-$T_{\text{c}}$ cuprates. [7]

Recently, 5$d$ transition metal oxide Sr2IrO4 and other related iridium oxides have attracted much interests both experimentally [8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25] and theoretically [26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52] because of the unexpected properties and novel phenomena arising from highly entangled spin and orbital degrees of freedom due to the large relativistic spin-orbit coupling (SOC) [53]. Sr2IrO4 is in a layered perovskite structure forming a two-dimensional square lattice composed of Ir atoms [54, 55], and thus it is similar to a parent compound of high-$T_{\text{c}}$ cuprate superconductors such as La2CuO4 [56]. The five 5$d$ electrons of Ir ion occupy the $t_{2g}$ orbitals, which are energetically lower than the $e_{g}$ orbitals due to the crystal field. The large SOC lifts the sixfold degenerate $t_{2g}$ orbitals (including spin degree of freedom) into the fourfold degenerate effective total angular momentum $j_{\text{eff}}=|-\bm{l}+\bm{s}|=3/2$ states and the doubly degenerate $j_{\text{eff}}=1/2$ states; Four electrons fill the $j_{\text{eff}}=3/2$ states and the remaining electron occupies one of the $j_{\text{eff}}=1/2$ states.

Several experiments have revealed that the states around the Fermi level is contributed mostly by these $j_{\text{eff}}=1/2$ states [10, 18, 21, 24], which are antiferromagnetically ordered in the two-dimensional plane at low temperatures with weak ferromagnetic moment [57, 58, 59, 60]. Furthermore, the antiferromagnetic exchange has been estimated to be as large as 50–100 meV [9, 10, 13, 14]. Therefore, there are many similarities between La2CuO4 and Sr2IrO4, as summarized in Table 1. Because of these similarities, the possibility of superconductivity in the Ir oxides has been proposed theoretically once mobile carriers are introduced into the $j_{\text{eff}}=1/2$ antiferromagnetic insulator [61, 62, 63, 64, 65]. Although the superconductivity has not been experimentally observed yet in these Ir oxides [66, 67, 68, 69, 70, 71, 72, 73], anomalous normal state properties similar to those observed often for under doped cuprates have been reported [74, 75, 76, 77, 78, 79, 80, 81, 82, 83].

In the previous theoretical studies, various numerical techniques such as variational Monte Carlo method, functional renormalization group method, and dynamical mean field theory with continuous time quantum Monte Carlo method are employed to show that a $d_{x^{2}-y^{2}}$-wave pairing is favored in the electron-doped region[61, 62] and a $s_{\pm}$-wave[64] or $p$-wave[63] pairing in the hole-doped region. However, in these studies, the parameter dependence of these SC states has not been thoroughly explored. This is precisely one of our purposes in this paper.

To investigate possible superconductivity in carrier doped 5$d$ transition metal oxide Sr2IrO4, we consider a $t_{2g}$ three-orbital Hubbard model with the SOC term on the square lattice as an effective model for Sr2IrO4. For this end, the linearized Eliashberg equation for the SC gap function is solved with the random phase approximation (RPA). To take account of the SOC in the RPA framework, we formulate the spin-dependent and multi-orbital RPA which can apply to a multi-orbital Hubbard model with the SOC. With this formulation, we discuss the possible pairing symmetries of the SC gap functions and their pairing mechanisms arising from the fluctuations based on the weak-coupling theory. The stability of the possible SC states is systematically studied by varying the model parameters such as the SOC $\lambda$, carrier concentration, and Hund’s coupling $J/U$.

Our weak-coupling analysis reveal that that a $d_{x^{2}-y^{2}}$-wave pseudospin $j_{\text{eff}}=1/2$ singlet pairing is favored for the parameters relevant to Sr2IrO4, namely, for large $\lambda$ and small $J/U$. We find that the $d_{x^{2}-y^{2}}$-wave pseudospin $j_{\text{eff}}=1/2$ singlet pairing is induced by the intra-band AF pseudospin $\bm{j}_{\text{eff}}$ fluctuations, and that the $d_{x^{2}-y^{2}}$-wave pairing is stabilized by the large SOC especially in the electron-doped region. In contrast, we find that a $s_{\pm}$-wave pseudospin singlet pairing is favored in the hole-doped region when the SOC is small and $J/U$ is large. The origin of the $s_{\pm}$-wave pairing is attributed to the inter-band pair scattering arising from the intra-orbital AF spin $\bm{s}$ fluctuations. We also find that the Hund’s coupling suppresses (enhances) the $d_{x^{2}-y^{2}}$-wave ($s_{\pm}$-wave pairing) pseudospin singlet pairing. Although the possibility of a pseudospin triplet pairing is considered, we find that it cannot be the most stable solution of the linearized Eliashberg equation with the maximum eigenvalue. We also argue that the enhancement of the $d_{x^{2}-y^{2}}$ pseudospin $j_{\text{eff}}=1/2$ singlet pairing with increasing the SOC $\lambda$ would be related to the purification of the orbital discussed in the context of cuprate superconductors [84, 85, 86]. We expect that the superconductivity with $s_{\pm}$-wave pairing symmetry found here is more relevant to hole doped Rh oxide Sr2RhO4.

The rest of this paper is organized as follows. In Sec. 2, we introduce a $t_{2g}$ three-orbital Hubbard model with the SOC as an effective model for Sr2IrO4, and derive the RPA scheme and the linearized Eliashberg equation for this multi-band model with the SOC. We also introduce the pseudospin picture which is a better description in the large SOC limit. In Sec. 3, we show the numerical results to examine the possible superconductivity, including the pairing symmetry, and the mechanism mediated by different fluctuations. We also explore the parameter dependence of the superconductivity on the SOC and Hund’s coupling. In Sec. 4, we conclude the paper with summary and discussion.

2 Model and Formulation

2.1 Effective model for Sr2IrO4

We consider an effective $t_{2g}$ three-orbital Hubbard model on the square lattice described by the following Hamiltonian:

[TABLE]

The first term $H_{\text{kin}}$ is the kinetic energy term,

[TABLE]

where

[TABLE]

and

[TABLE]

are the Fourier transformation of the creation and annihilation operators $c^{a\dagger}_{i\sigma}$ and $c^{a}_{i\sigma}$, respectively, of an electron in $t_{2g}$ orbital $a$ ($=yz,zx,xy$) with spin $\sigma$ ($=\uparrow,\downarrow$) at position $\bm{R}_{i}$ of site $i$ on the square lattice composed of $N$ sites. Here $\epsilon^{a}_{\bm{k}}$ represents the energy dispersion for $t_{2g}$ orbital $a$ and $\mu$ is the chemical potential. In this paper, we consider the simplest case in the absence of structural distortion, i.e., rotation of IrO6 octahedra, [54, 55] as in Ba2IrO4. [87, 88]

The second term $H_{\text{SO}}$ is the SOC term with the coupling constant $\lambda$ described as

[TABLE]

where $\bm{l}$ ($\bm{s}$) is the orbital (spin) angular momentum operator, $\bar{\sigma}$ the opposite spin of $\sigma$, and $s_{\sigma}=1\,(-1)$ for $\sigma=\uparrow$ ($\downarrow$).

The third term $H_{\text{int}}$ is the Coulomb interaction term,

[TABLE]

which contains the intra- and inter-orbital on-site Hubbard interactions ($U$ and $U^{\prime}$, respectively), Hund’s coupling ($J_{\rm H}$), and pair hopping ($J^{\prime}$). Here $n^{a}_{i\sigma}=c^{a\dagger}_{i\sigma}c^{a}_{i\sigma}$, $n^{a}_{i}=n^{a}_{i\uparrow}+n^{a}_{i\downarrow}$, and $\bm{S}^{a}_{i}=\sum_{\sigma,\sigma^{\prime}}\langle\sigma|\bm{s}|\sigma^{\prime}\rangle c^{a\dagger}_{i\sigma}c^{b}_{i\sigma^{\prime}}$ is the spin operator at site $i$. We impose the constraint for $d$-orbital electrons in the atomic limit by setting $U=U^{\prime}+2J$ and $J=J_{\rm H}=J^{\prime}$. [89, 90] Although the realistic $U$ value for Ir oxides is estimated as large as 2–3 eV, [30] we set smaller $U$ because the RPA calculation usually overestimates the instability toward ordered states.

2.2 Noninteracting limit

The one-body part of the Hamiltonian, $H_{0}=H_{\text{kin}}+H_{\text{SO}}$, is written in the matrix form as

[TABLE]

According to the previous studies [28, 62], the energy band structure of Sr2IrO4 obtained by the first-principles calculation [9] can be well reproduced using the following tight-binding energy dispersion for the three $t_{2g}$ orbitals:

[TABLE]

and

[TABLE]

with the parameters

[TABLE]

Diagonalizing the ($3\times 3$) Hamiltonian matrix in Eq. (7), we can obtain the three doubly degenerate energy band dispersions $E^{m}_{\bm{k}}$ ($m=1,2,3$). Note that each band is doubly degenerate because of the time-reversal symmetry and the inversion symmetry.

In this study, we consider the SOC $\lambda$, electron concentration $n$, and Hund’s coupling $J$ as parameters to systematically study the possibility of superconductivity in Sr2IrO4. Figure 1 shows the energy band dispersion, the density of states (DOS), and Fermi surfaces (FSs) for $H_{0}$ with a parameter set given in Eq. (10) but varying $\lambda$. The Fermi level $E_{\text{F}}$ for different $n$ is also indicated in Fig. 1. The energy band basis that diagonalizes the Hamiltonian matrix in Eq. (7) is obtained with the $\bm{k}$-dependent ($6\times 6$) unitary matrix $\mathcal{U}^{\dagger}_{\bm{k}}$ as

[TABLE]

where $d^{m\dagger}_{\bm{k}\eta}$ and $d^{m}_{\bm{k}\eta}$ denote the creation and annihilation operator of an electron in the $m$-th band with $\eta\,(=\Uparrow,\Downarrow)$ indexing the doubly degenerate states.

It is instructive to consider the basis set in the atomic limit, i.e., $\epsilon^{a}\to 0$ or $\lambda\rightarrow\infty$, which is achieved by diagonalizing $H_{\rm SO}$. This basis set is also the eigenstates of the pseudospin $\bm{j}_{\text{eff}}=-\bm{l}+\bm{s}$ and its $z$ component with the eigenvalues $j_{\text{eff}}$ and $j^{z}_{\text{eff}}$, respectively, thus often called the effective total angular momentum $j_{\text{eff}}$ basis. This basis set is obtained by the site independent unitary transformation,

[TABLE]

where the $6\times 6$ unitary matrix $\mathcal{W}^{\dagger}$ is given as

[TABLE]

The creation and annihilation operators $a^{m\dagger}_{i\eta}$ and $a^{m}_{i\eta}$ in the $j_{\text{eff}}$ basis with $(m,\eta)=(1,\Uparrow\Downarrow)$, $(2,\Uparrow\Downarrow)$, and $(3,\Uparrow\Downarrow)$ correspond to $(j_{\text{eff}},j^{z}_{\text{eff}})=(1/2,\pm 1/2)$, $(3/2,\mp 3/2)$, and $(3/2,\pm 1/2)$, respectively. In the atomic limit, these two unitary matrices becomes the same, i.e., $\mathcal{U}^{\dagger}_{\bm{k}}\to\mathcal{W}^{\dagger}$, and the energy band eigenstates described by $d^{m\dagger}_{\bm{k}\eta}$ and $d^{m}_{\bm{k}\eta}$ with $(m,\eta)=(1,\Uparrow\Downarrow)$, $(2,\Uparrow\Downarrow)$, and $(3,\Uparrow\Downarrow)$ approach to the $j_{\text{eff}}$ bases with $(j_{\rm eff},j_{\rm eff}^{z})=(1/2,\pm 1/2)$, $(3/2,\mp 3/2)$, and $(3/2,\pm 1/2)$.

In fact, we can observe in Fig. 2 that the energy band $E^{m}_{\bm{k}}$ with $m=1$, which crosses the Fermi level at $n=5$, forming electron-like FS, and is the most important for the later discussion, is approximately described by the $j_{\text{eff}}=1/2$ basis when $\lambda=0.37\ \text{eV}\sim 1.03t_{1}$, relevant to Sr2IrO4, because almost all the spectral weight projected onto the $j_{\text{eff}}=1/2$ basis (i.e., partial spectral function) is concentrated on the band $E^{m}_{\bm{k}}$ with $m=1$. On the other hand, the spectral weight projected onto the $j_{\text{eff}}=3/2$ bases are concentrated mostly on the bands $E^{m}_{\bm{k}}$ with $m=2$ and 3. We can also notice that the band $E^{m}_{\bm{k}}$ with $m=2$ crosses the Fermi level at $n=5$ with this $\lambda$ value, forming the hole-like FS around $\bm{k}=(\pi,\pi)$ [see Fig. 1(b)], which has the strong character of $(j_{\text{eff}},j^{z}_{\text{eff}})=(3/2,\pm 3/2)$.

The number of FSs crossing the Fermi level depends on the SOC $\lambda$ as well as the electron concentration $n$. Since the SOC tends to increase the separation between the $j_{\rm eff}=1/2$ based band and the $j_{\rm eff}=3/2$ based bands, there is only one FS crossing the Fermi level for the larger SOC when $n$ is around 5, as shown in Fig. 1. It is also apparent that when the number of electron increases, only the highest band, i.e., the $j_{\rm eff}=1/2$ based band, is relevant and thus there is only a single FS crossing the Fermi level (see Fig. 1). This feature is summarized in Fig. 3 for the case of the tight-binding parameters given in Eq. (10) except that here the SOC $\lambda$ is varied. In Fig. 3, we also indicate $(\lambda,n)$ at which the van Hove singularity (vHS) is located exactly at the Fermi level. The vHS is originated from the inflection point of the energy band $E^{m}_{\bm{k}}$ with $m=1$ at $\bm{k}=(\pm\pi,0)$ and $(0,\pm\pi)$, which is always located in the electron-doped side for the parameters considered here.

2.3 Linearized Eliashberg equation in RPA

Here, we show the procedure of the RPA calculations to treat the electron correlations. The RPA incorporates the spin and charge fluctuations by including the infinite summation of bubble and ladder diagrams [91, 92, 93, 94, 95]. Since our system $H$ is a multi-band system with the SOC, it requires a spin-dependent multi-band formalism, [96] in which the spin SU(2) symmetry is broken due to the presence of the SOC.

Let us first derive the linearized Eliashberg equation for a SC gap function. The SC properties such as the pairing symmetry and the transition temperature can be determined by solving the linearized Eliashberg equation. For this purpose, we first define the normal Green’s functions $\hat{G}(k)$ and $\hat{\bar{G}}(k)$ in the multi-band system as

[TABLE]

where $k=(\bm{k},i\omega_{n})$ denotes a set of wave vector $\bm{k}$ and Matsubara frequency $\omega_{n}=(2n+1)\pi/\beta$ for fermions ($n$: integer), and $T$ ($\beta=1/k_{\text{B}}T$) denotes the (inverse) temperature. Note that we temporarily treat the Greek indices as a set of the orbital and spin degrees of freedom: $\alpha=(a,\sigma)$ contains the orbital index $a=yz,zx,xy$ and the spin index $\sigma=\uparrow,\downarrow$, and thus $c_{\bm{k}\alpha}=c^{a}_{\bm{k}\sigma}$ and $c^{\dagger}_{\bm{k}\alpha}=c^{a\dagger}_{\bm{k}\sigma}$. The time dependence of the field operators in the imaginary-time $\tau$ formalism is defined as $c_{\bm{k}\alpha}(\tau)=\mathrm{e}^{H\tau}c_{\bm{k}\alpha}\mathrm{e}^{-H\tau}$ and $c^{\dagger}_{\bm{k}\alpha}(\tau)=\mathrm{e}^{H\tau}c^{\dagger}_{\bm{k}\alpha}\mathrm{e}^{-H\tau}$. $\langle\cdots\rangle=\mathrm{Tr}[\mathrm{e}^{-\beta H}\cdots]/Z$ is the statistical average of the system $H$ and $Z=\mathrm{Tr}\ \mathrm{e}^{-\beta H}$ is the partition function. From the definition of $\hat{G}(k)$ and $\hat{\bar{G}}(k)$, we can easily show that $\bar{G}_{\alpha\beta}(k)=-G_{\beta\alpha}(-k)$.

The noninteracting Green’s function is obtained from the one-body dispersion matrix $\hat{\xi}_{\bm{k}}$ as

[TABLE]

where the matrix elements of $\hat{\xi}_{\bm{k}}$ can be given by rewriting the one-body Hamiltonian $H_{0}$ in this representation as

[TABLE]

The noninteracting Green’s functions also satisfy that $\bar{G}_{0\,\alpha\beta}(k)=-G_{0\,\beta\alpha}(-k)$.

We also define the anomalous Green’s functions $\hat{F}(k)$ and $\hat{\bar{F}}(k)$ in the multi-band system as

[TABLE]

The SC order parameter, i.e., the SC gap function $\hat{\Delta}(k)$ can be interpreted as the self-energy of the anomalous Green’s function $\hat{F}(k)$. From the Dyson’s equation for $\hat{G}(k)$ and $\hat{F}(k)$, we can obtain the following relation:

[TABLE]

where we assume that the temperature is around the critical temperature $T_{\text{c}}$ and thus $\hat{F}(k)$ and $\hat{\Delta}(k)$ are infinitesimally small.

The SC gap function $\hat{\Delta}(k)$ can be obtained from the functional derivative of the Luttinger–Ward functional $\Phi$ with respect to the anomalous Green’s function $\hat{\bar{F}}(k)$ [97, 98, 99]. Here the SC gap function can be given as

[TABLE]

where $\hat{V}_{\text{P}}(q)$ represents the SC pairing interaction, i.e., the effective interaction that induces the self-energy of the anomalous Green’s function $\hat{F}(k)$ (see Fig. 4), and $q=(\bm{q},i\varepsilon_{m})$ with $\varepsilon_{m}=2m\pi/\beta$ being a Matsubara frequency for bosons ($m$: integer). $\hat{V}_{\text{P}}(q)$ can be uniquely determined once the Luttinger–Ward functional $\Phi$ is chosen. In the RPA, the infinite summation of bubble and ladder diagrams connected to each other via the Coulomb interactions are taken into account in $\Phi$.

By combining Eqs. (18) and (19), we can obtain the linearized Eliashberg equation. To solve this equation, first we neglect the frequency dependence of the gap function, i.e., $\hat{\Delta}(k)\rightarrow\hat{\Delta}(\bm{k})$, and the pairing interaction, i.e., $\hat{V}_{\text{P}}(q)\rightarrow\hat{V}_{\text{P}}(\bm{q})\equiv\hat{V}_{\text{P}}(\bm{q},i\varepsilon_{m}=0)$, for simplicity. In addition to this, in the RPA calculation, the Green’s functions $\hat{G}(k)$ and $\hat{\bar{G}}(k)$ are replaced with the noninteracting Green’s functions $\hat{G}_{0}(k)$ and $\hat{\bar{G}}_{0}(k)$, respectively. Thus, the linearized Eliashberg equation can be given as

[TABLE]

where the summation over the Matsubara frequencies $\omega_{n}$ is gathered in $\hat{\phi}$ as

[TABLE]

Note that we have to impose the constraint

[TABLE]

on the solution $\hat{\Delta}(\bm{k})$ of the linearized Eliashberg equation in Eq. (20) because of the anticommutation relation for fermions. Note also that the mixture of singlet and triplet pairing is permitted as a consequence of the spin SU(2) symmetry breaking due to a finite SOC. Moreover, it is convenient to consider the eigenvalue $\lambda_{\text{e}}$ of the linearized Eliashberg equation in Eq. (20). The SC transition occurs when $\lambda_{\text{e}}=1$ at the SC critical temperature $T_{\text{c}}$. The value of $\lambda_{\text{e}}$ itself gives us a criterion for the strength of the superconductivity at temperature $T$ around $T_{\text{c}}$.

To obtain the SC paring interaction $\hat{V}_{\text{P}}(q)$, we next consider the dynamical correlation function $\hat{\chi}(q)$ defined as

[TABLE]

where the subscript “c” indicates that only connected diagrams are included. In the weak-coupling theory, $\hat{\chi}(q)$ can be evaluated from the infinite summation of the bubble and ladder diagrams connected to each other via the Coulomb interactions, which is a solution of the Bethe–Salpeter equation for the particle-hole correlation function. Furthermore, in the RPA calculations, the bubble and ladder diagrams are replaced with those composed of the noninteracting Green’s functions $\hat{G}_{0}(k)$. One single bubble/ladder diagram in the RPA thus corresponds to the noninteracting dynamical correlation function,

[TABLE]

and this bubble/ladder diagram is connected to others repeatedly via the Coulomb interaction $\hat{V}$ in $\hat{\chi}(q)$.

Showing explicitly the orbital and spin indices, we can now express the noninteracting dynamical correlation function $\chi_{0\,\alpha\beta\gamma\delta}(q)=\chi^{abcd}_{0\,\sigma_{1}\sigma_{2}\sigma_{3}\sigma_{4}}(q)$ in a $2^{2}\times 2^{2}$ block matrix form with respect to the spin indices (the number of spin indices being $2$), each element of the block matrix being composed of $\hat{\chi}_{0\,\sigma_{1}\sigma_{2}\sigma_{3}\sigma_{4}}(q)$, i.e.,

[TABLE]

where the four rows (columns) correspond to $\sigma_{1}\sigma_{2}\,\,(\sigma_{3}\sigma_{4})=\uparrow\uparrow,\,\uparrow\downarrow,\,\downarrow\uparrow,\,\downarrow\downarrow$. Each element of the block matrix, $\hat{\chi}_{0\,\sigma_{1}\sigma_{2}\sigma_{3}\sigma_{4}}(q)$, is a $3^{2}\times 3^{2}$ matrix $\chi^{abcd}_{0\,\sigma_{1}\sigma_{2}\sigma_{3}\sigma_{4}}(q)$ with respect to the orbital indices (the number of orbital indices being $3$ in the three-orbital case), where the nine rows (columns) correspond to $ab\,\,(cd)=yz\ yz,\,yz\ zx,\,yz\ xy,\,zx\ yz,\,zx\ zx,\,zx\ xy,\,xy\ yz,\,xy\ zx,\,xy\ xy$.

In the same way, the Coulomb interaction $\hat{V}$ can also be expressed in a $2^{2}\times 2^{2}$ block matrix form, composed of $\hat{V}_{\sigma_{1}\sigma_{2}\sigma_{3}\sigma_{4}}$, with respect to the spin indices, i.e.,

[TABLE]

where each element of the block matrix, $\hat{V}_{\sigma_{1}\sigma_{2}\sigma_{3}\sigma_{4}}$, can be given explicitly as

[TABLE]

and other elements of the block matrix are zero.

Therefore, these matrices $\hat{\chi}_{0}(q)$ and $\hat{V}$ can be multiplied simply by the matrix (or more precisely tensor) product with respect to the spin and orbital indices, i.e.,

[TABLE]

Thus, the dynamical correlation function $\hat{\chi}(q)$ in the RPA, which is the infinite series $\hat{\chi}(q)=\hat{\chi}_{0}(q)+\hat{\chi}_{0}(q)\hat{V}\hat{\chi}_{0}(q)+\cdots=\hat{\chi}_{0}(q)+\hat{\chi}_{0}(q)\hat{V}\hat{\chi}(q)$, can be given as

[TABLE]

In the RPA, the superconductivity is mediated by spin and orbital (charge) fluctuations that are included in the dynamical correlation function $\hat{\chi}(q)$. Therefore, the SC pairing interaction $\hat{V}_{\text{P}}(q)$ is finally given as

[TABLE]

where the factor $1/2$ in the first term is required in order to avoid the double counting in the first order diagrams.

3 Numerical Results

To investigate the possible superconductivity in the $t_{2g}$ three-orbital Hubbard model with the SOC, an effective model for Sr2IrO4, we solve the linearized Eliashberg equation numerically with the RPA. In the numerical calculations, we set the number of momentum mesh points in the first Brillouin zone as large as $64\times 64$ points and the temperature $T=0.03$ eV $\sim 300$ K.

3.1 Possible Superconductivity

For the analysis of possible superconductivity, it is useful to represent the SC gap function $\hat{\Delta}(\bm{k})$, represented in terms of the $t_{2g}$ basis in Eq. (20), by using the energy band basis $m$ that diagonalizes the noninteracting Hamiltonian $H_{0}$, i.e.,

[TABLE]

where the unitary matrix $\hat{\mathcal{U}}_{\bm{k}}$ is introduced in Eq. (11). Recall that the pseudospin ($=\Uparrow,\Downarrow$), specifying the doubly degenerated energy band, is a good quantum number, not the electron spin $\sigma\,(\uparrow,\downarrow)$, because of the presence of the SOC. Since there exists the spatial inversion symmetry in the Hamiltonian $H$, the SC gap function $\hat{D}(\bm{k})$ represented in terms of the energy band basis $m$ can be characterized as the “pseudospin” singlet or triplet pairing. This is an analogy to the intrinsic spin singlet or triplet pairing when the SOC is absent. In the numerical calculations, we can distinguish the pseudospin singlet and triplet pairings by imposing the even or odd parity on the SC gap function, i.e.,

[TABLE]

Note that this is equivalent to $\hat{\Delta}(-\bm{k})=\pm\left[\hat{\Delta}(\bm{k})\right]^{T}$.

Figure 5 and Figure 6 show the results of the largest eigenvalue $\lambda_{\text{e}}$ of the linearized Eliashberg equation with the pseudospin singlet and triplet pairings for the representative two sets of the Coulomb interactions $U$ and $J$. We find that the pseudospin singlet pairing is always favored as compared to the pseudospin triplet pairing for all the parameter sets studied here because the maximum eigenvalue $\lambda_{\text{e}}$ for the pseudospin singlet pairing is always larger than that for the pseudospin triplet pairing. It is also interesting to notice that the $\lambda$ dependence of the eigenvalue $\lambda_{\text{e}}$ is opposite in Fig. 5 and Fig. 6. The eigenvalue $\lambda_{\text{e}}$ increases (decreases) with $\lambda$ for a given electron concentration $n$ in Fig. 5 (Fig. 6) where a relatively small (large) $J$ is set.

Furthermore, as shown in Fig. 7, we find that the dominant SC gap function has $d_{x^{2}-y^{2}}$-wave symmetry for the pseudospin singlet pairing and $p$-wave symmetry for the pseudospin triplet pairing in the parameter region shown in Fig. 5. These pairings are formed mainly in the $m=1$ energy band that has a strong character of $j_{\rm eff}=1/2$ (see Fig. 2). These pairings are much favored for larger $\lambda$ and particularly in the electron doped region. This is because the single band nature is enhanced with increasing $\lambda$ in the electron doped region and also because there is the vHS in the electron-doped side (see Fig. 3).

On the other hand, as shown in Fig. 8, we find that the dominant SC gap functions have $s_{\pm}$-wave symmetry for the pseudospin singlet pairing and $p_{\pm}$-wave symmetry for the pseudospin triplet pairing in the parameter region shown in Fig. 6. These pairing are formed in the $m=1$ and $m=2$ energy bands that have the strong character of $j_{\rm eff}=1/2$ and $j_{\rm eff}=3/2$, respectively (see Fig. 2). Notice also in Fig. 8 that the SC gap functions alter the sign between $m=1$ and $m=2$. Therefore, the symmetry of the SC gap functions is $s_{\pm}$-wave or $p_{\pm}$-wave, instead of simply $s$-wave or $p$-wave. These pairings are favored for small $\lambda$ where the two band character is enhanced at the Fermi level by forming the large electron FS and the relatively small hole FS in the noninteracting limit (see Fig. 1).

3.2 Pseudospin Pairing

To better understand the internal structure of the pseudospin pairing, let us decompose the pseudospin singlet pairing found in Figs. 5(a) and 6(a). The dominant components of the numerically obtained SC gap function $\hat{D}(\bm{k})$ can be decomposed into a mixture of the spin singlet and triplet pairing when it is represented in the $t_{2g}$ orbital basis with $\Delta_{\alpha\beta}(\bm{k})=\Delta^{ab}_{\sigma\sigma^{\prime}}(\bm{k})$. As shown in Table 2, the intra-orbital elements $\Delta^{aa}_{\sigma\bar{\sigma}}(\bm{k})$ $(a=yz,zx,xy)$ are spin singlet, i.e., $\Delta^{aa}_{\sigma\bar{\sigma}}(\bm{k})=-\Delta^{aa}_{\bar{\sigma}\sigma}(\bm{k})$, while the inter-orbital ones are spin triplet, i.e, $\Delta^{ab}_{\sigma\sigma^{\prime}}(\bm{k})=\Delta^{ab}_{\sigma^{\prime}\sigma}(\bm{k})$ for $a\neq b$, but the SC gap function changes the sign when the $t_{2g}$ orbital indices are exchanged, i.e., $\Delta^{yz\,zx}_{\sigma\bar{\sigma}}(\bm{k})=-\Delta^{zx\,yz}_{\sigma\bar{\sigma}}(\bm{k})$, $\Delta^{yz\,xy}_{\sigma\sigma}(\bm{k})=-\Delta^{xy\,yz}_{\sigma\sigma}(\bm{k})$, and $\Delta^{zx\,xy}_{\sigma\sigma}(\bm{k})=-\Delta^{xy\,zx}_{\sigma\sigma}(\bm{k})$ for the dominant components. Since the spin SU(2) symmetry of the Hamiltonian is broken due to the presence of the SOC, the minus sign of the SC gap function $\hat{D}(\bm{k})$ for the psuedospin singlet pairing in Eq. (32) can be imposed on either orbital or spin degree of freedom in $\Delta_{\alpha\beta}(\bm{k})=\Delta^{ab}_{\sigma\sigma^{\prime}}(\bm{k})$ represented in terms of the $t_{2g}$ basis to satisfy $\hat{\Delta}(-\bm{k})=-\left[\hat{\Delta}(\bm{k})\right]^{\rm T}$.

The dominant components of the SC gap function $\Delta^{ab}_{\sigma\sigma^{\prime}}(\bm{k})$ found in the numerically obtained pseudospin singlet SC gap function $D^{mm^{\prime}}_{\eta\eta^{\prime}}(\bm{k})$ is listed in Table 2. Here, $m$ ($m^{\prime}$) and $\eta$ ($\eta^{\prime}$) are the band and pseudospin indices, respectively, for the energy band basis that diagonalizes the noninteracting Hamiltonian $H_{0}$ in Eq. (7), as defined in Eq. (11). To understand these dominant components, let us consider the $d_{x^{2}-y^{2}}$-wave pseudospin singlet pairing formed by the energy band $m=1$ in the large SOC limit, as a simple example. As described in Sec. 2.2, the energy band basis is smoothly connected to the local atomic $j_{\rm eff}$ basis in the large SOC limit. Therefore, in this limit, the $d_{x^{2}-y^{2}}$-wave pseudospin singlet pairing is formed by $a^{1}_{\bm{k}\Uparrow}$ and $a^{1}_{-\bm{k}\Downarrow}$, where the index “1” corresponds to $j_{\rm eff}=1/2$ [see Eq. (12)]. Since the local atomic $j_{\rm eff}$ basis and the $t_{2g}$ orbital basis are transformed via the unitary matrix $\mathcal{W}$ in Eq. (13), the $d_{x^{2}-y^{2}}$-wave pseudospin singlet pairing in the large SOC limit is represented as

[TABLE]

We can now observe that the $d_{x^{2}-y^{2}}$-wave pseudospin singlet pairing can be expressed as a linear combination of the spin singlet and triplet pairings formed by the $t_{2g}$ orbitals and these components correspond exactly to those of $\Delta^{ab}_{\sigma\sigma^{\prime}}(\bm{k})$ listed in Table 2. In other words, this is a consequence of the fact that the pseudospin singlet superconductivity found in our calculations is formed essentially by the pseudospin $j_{\rm eff}=1/2$ quasiparticles.

3.3 $d_{x^{2}-y^{2}}$-wave pseudospin singlet pairing

In order to understand the pairing mechanism of the pseudospin singlet pairing, here we calculate, in analogy to the intrinsic spin susceptibility, the the longitudinal (out-of-plane) pseudospin susceptibility

[TABLE]

and the transverse (in-plane) pseudospin susceptibility

[TABLE]

where the $z$-component of the pseudospin operator $J^{z\,m}_{\text{eff}\,\bm{q}}$ is defined as

[TABLE]

and the correspond transverse components $J^{\pm\,m}_{\text{eff}\,\bm{q}}$ are defined as

[TABLE]

and

[TABLE]

with $\left[J^{-\,m}_{\text{eff}\,\bm{q}}\right]^{\dagger}=J^{+\,m}_{\text{eff}\,-\bm{q}}$. Here, $a^{m\dagger}_{\bm{k}\eta}$ ($a^{m}_{\bm{k}\eta}$) is the Fourier transform of the $j_{\rm eff}$ operator $a^{m\dagger}_{i\eta}$ ($a^{m}_{i\eta}$) defined in Eq. (12).

Figure 9 shows the pseudospin susceptibilities for $\lambda=0.5$ eV, $U=0.75$ eV, $J/U=0.05$, and $n$=5. Note that the susceptibilities are multiplied by $t_{1}$ to make them dimensionless. Although the transverse and longitudinal pseudospin susceptibilities are quantitatively different due to a finite Hund’s coupling, [26] they both exhibit peak structures around momentum $\bm{Q}=(\pi,\pi)$ and other symmetrically equivalent momenta, and thus the AF pseudospin instability is enhanced for this set of parameters. This instability can be understood by the FS topology shown in Fig. 1(a): The shape of the FS satisfies nearly the perfect nesting condition with the nesting vector $\bm{Q}=(\pi,\pi)$ and other symmetrically equivalent momenta. In addition, we notice in Fig. 9 that the transverse pseudospin susceptibility shows a larger peak value than the longitudinal one. This implies that the in-plane AF order is more favored than the out-of-plane one, which is consistent with the experimental observation for Sr2IrO4 [9, 10].

Figure 10 displays the $n$ dependence of the transverse pseudospin susceptibility. As $n$ departs away from $n=5$, the peak values monotonically decrease and thus the pseudospin AF fluctuations are suppressed. Accordingly, the pairing interaction forming the Cooper pairs becomes weaker. We can also notice in Fig. 10 that the doping dependence of the pseudospin susceptibility from $n=5$ is qualitatively different for the electron-doped case ($n>5$) and the hole-doped case ($n<5$). The single peak appearing at $\bm{Q}=(\pi,\pi)$ for $n=5$ moves slightly along $x$ and $y$ directions to $(\pi,\pi\pm\delta)$ and $(\pi\pm\delta,\pi)$ for the electron doping $n>5$, while the single peak structure remains at $\bm{Q}=(\pi,\pi)$ for the hole doping $n<5$. Similar behavior of the spin susceptibility is found in the single-band Hubbard model for cuprate superconductors [99]. Note also that the hole pockets around $\bm{Q}=(\pi,\pi)$ and other symmetrically equivalent momenta appear for $n<4.9$ in the noninteracting limit (see Fig. 3), which also reduces significantly $\chi^{11}_{j^{\pm}_{\text{eff}}}(\bm{q},i\epsilon_{m}=0)$ and $\chi^{11}_{j^{zz}_{\text{eff}}}(\bm{q},i\epsilon_{m}=0)$ since other scattering processes arise due to the appearance of the additional Fermi surface.

From these results of the pseudospin susceptibilities, we can now understand the doping and $\lambda$ dependence of $\lambda_{\text{e}}$ in Fig. 5(a). For the electron doping with $n>5$, the DOS at the Fermi level increases toward the vHS around $n=5.2$–$5.3$ depending on $\lambda$ (see Fig. 3), while the pseudospin AF instability decreases with departing away from $n=5$. As a result of this competition, $\lambda_{\text{e}}$ exhibits a maximum around $n\sim 5.1$, as shown in Fig. 5(a). Furthermore, the peak value of $\lambda_{\text{e}}$ is enhanced with increasing $\lambda$, especially for $\lambda>0.4$ eV. This is because the hole pockets around $\bm{k}=(\pi,\pi)$ and symmetry equivalent momenta disappear for $\lambda>0.4$ eV and only the single band crosses the Fermi level in the noninteracting limit (Fig. 3). Therefore, it is understood that the $d_{x^{2}-y^{2}}$-wave pseudospin singlet pairing found here in our calculations is mediated by the AF pseudospin fluctuations. This is an analogy to cuprate superconductors where the $d_{x^{2}-y^{2}}$-wave spin $s=1/2$ singlet pairing is mediated by the AF spin fluctuations.

One of the important differences from the spin $s=1/2$ singlet pairing in cuprate superconductors is that the pseudospin singlet pairing found here is rather sensitive to the Hund’s coupling $J$. Figure 11 shows the $J/U$ dependence of $\lambda_{\text{e}}$ for the electron doping at $n=5.2$ with several $\lambda$ values. As we have discussed above, the $d_{x^{2}-y^{2}}$-wave pseudospin $j_{\text{eff}}=1/2$ singlet pairing becomes weaker with decreasing $\lambda$. Moreover, it is also suppressed with increasing $J/U$. This is because of the suppression of the pseudospin AF fluctuations due to the decrease of the inter-orbital on-site interaction $U^{\prime}$ through $U^{\prime}=U-2J$.

Figure 12 shows the longitudinal and transverse components,

[TABLE]

and

[TABLE]

respectively, of the total pseudospin susceptibility for small and large values of $J$, i.e., $J/U=0.05$ and $J/U=0.15$. It is clearly observed that the pseudospin susceptibilities are significantly suppressed with increasing $J/U$, and accordingly the $d_{x^{2}-y^{2}}$-wave pseudospin singlet pairing is also suppressed.

Finally, we systematically explore the stability of the $d_{x^{2}-y^{2}}$-wave pseudospin singlet pairing in the parameter space of the electron concentration $n$ and the SOC $\lambda$. The obtained phase diagram in Fig. 13 shows the contour plots of the largest eigenvalue of the linearized Eliashberg equation for the $d_{x^{2}-y^{2}}$-wave pseudospin singlet pairing in the parameter set of $U=0.8$ eV and $J/U=0.05$, where the $d_{x^{2}-y^{2}}$-wave pseudospin singlet pairing is favored over the $s_{\pm}$-wave pseudospin singlet pairing. The red colored region represents the $d_{x^{2}-y^{2}}$-wave pseudospin singlet SC phase and the blue one is the $x$-AF phase where the AF spin moment is aligned within the $xy$-plane. The SC and AF phases here are identified when $\lambda_{\text{e}}\geq 1$ and $\max[\hat{V}\hat{\chi}_{0}(q)]\geq 1$, respectively.

3.4 $s_{\pm}$-wave pseudospin singlet pairing

As we discussed in Sec. 3.1, the $s_{\pm}$-wave pseudospin singlet pairing is more favored than the $d_{x^{2}-y^{2}}$-wave pseudospin singlet pairing when the single band character is less pronounced and the Hund’s coupling $J/U$ is large. To understand the pairing mechanism of the $s_{\pm}$-wave pseudospin singlet pairing, let us consider the longitudinal and transverse spin susceptibilities,

[TABLE]

and

[TABLE]

where the longitudinal and transverse components of the spin operator, $S^{z\,a}_{\bm{q}}$ and $S^{\pm\,a}_{\bm{q}}$, for $t_{2g}$ orbital $a\,=(yz,\,zx,\,xy)$ are given as

[TABLE]

[TABLE]

and

[TABLE]

with $\left[S^{-\,a}_{\bm{q}}\right]^{\dagger}=S^{+\,a}_{-\bm{q}}$.

Figure 14(a) shows the transverse spin susceptibility in the $yz$ orbital $\chi^{yz\,yz}_{S^{\pm}}(\bm{q},\mathrm{i}\epsilon_{m}=0)$ for $n=4.7$ with $\lambda=0.3$ eV, $U=0.9$ eV, and $J/U=0.3$. Note that the spin susceptibility in the $zx$ orbital $\chi^{zx\,zx}_{S^{\pm}}(\bm{q},\mathrm{i}\epsilon_{m}=0)$ is obtained from $\chi^{yz\,yz}_{S^{\pm}}(\bm{q},\mathrm{i}\epsilon_{m}=0)$ simply by rotating the $yz$-orbital by 90 degrees. We can clearly observe in Fig. 14(a) that $\chi^{yz\,yz}_{S^{\pm}}(\bm{q},\mathrm{i}\epsilon_{m}=0)$ exhibits peaks at $\bm{Q}^{\prime}\sim(5\pi/8,5\pi/8)$ and symmetrically equivalent momenta. Note however that the spin susceptibility in the $xy$-orbital $\chi^{xy\,xy}_{S^{\pm}}(\bm{q},\mathrm{i}\epsilon_{m}=0)$ is found small compared to those in the $yz$ and $yz$ orbitals simply because the $xy$-orbital energy level is located lower than the $yz$- and $yz$-orbital energy levels and thus it is almost complete occupied by electrons. As shown in Fig. 14(b), this momentum $\bm{Q}^{\prime}$ corresponds to the nesting vector connecting the electron-like and hole-like FSs with the large component of the $yz$ orbital character. This nesting feature is exactly the same for the $zx$ orbital [see Fig. 14(b)], which thus causes the same four peak structures in $\chi^{zx\,zx}_{S^{\pm}}(\bm{q},\mathrm{i}\epsilon_{m}=0)$. Similarly, the longitudinal spin susceptibilities in the $yz$ and $zx$ orbitals exhibits four peaks at $\bm{Q}^{\prime}$ and symmetry equivalent momenta. Furthermore, we find that the nesting vector $\bm{Q}^{\prime}$ tends to change the sigh of the SC gap functions because the inter-band scattering through the nesting vector $\bm{Q}^{\prime}$ is repulsive and thus the sign of the SC gap functions in the energy bands $m=1$ and $2$ is opposite. Therefore, the spin fluctuations with the nesting vector $\bm{Q}^{\prime}$ induce the $s_{\pm}$-wave pseudospin singlet pairing. This is analogous to the pairing mechanism for the iron-based superconductors: they also possess multi-band FSs and the inter-band repulsive scattering is expected to be dominant, which favors the $s_{\pm}$-wave intrinsic spin singlet pairing. [100, 101]

It is now apparent that the two band feature is essential for stabilizing the $s_{\pm}$-wave pseudospin singlet pairing, which explains why the small $\lambda$ value is favored for this singlet pairing. In addition, we find that the $s_{\pm}$-wave pseudospin singlet pairing is more favored when the Hund’s coupling $J$ is large. Figure 15 shows the $J/U$ dependence of the eigenvalue $\lambda_{\text{e}}$ of the linearized Eliashberg equation in the hole doping at $n=4.7$ for several $\lambda$ values. We can clearly observe that the $s_{\pm}$-wave pseudospin singlet pairing is enhanced with increasing $J/U$ as well as with decreasing $\lambda$. This is because the spin fluctuations with the nesting vector $\bm{Q}^{\prime}$ is enhanced as $J/U$ is increased.

Figure 16 shows the longitudinal and transverse components of the total spin susceptibility,

[TABLE]

and

[TABLE]

respectively, for small and large values of $J$, i.e., $J/U=0.15$ and $J/U=0.3$. It is indeed clearly observed that the total spin susceptibilities are significantly enhanced with increasing $J/U$. Accordingly, the $s_{\pm}$-wave pseudospin singlet pairing is also enhanced.

Finally, we systematically explore the stability of the $s_{\pm}$-wave pseudospin singlet pairing in the parameter space of the electron concentration $n$ and the SOC $\lambda$. The obtained phase diagram in Fig. 17 shows the contour plots of the largest eigenvalue of the linearized Eliashberg equation for the $s_{\pm}$-wave pseudospin singlet pairing in the parameter set of $U=0.9$ eV and $J/U=0.3$, where the $s_{\pm}$-wave pseudospin singlet pairing is favored over the $d_{x^{2}-y^{2}}$-wave pseudospin singlet pairing. The red colored region represents the $s_{\pm}$-wave pseudospin singlet SC phase where $\lambda_{\text{e}}\geq 1$.

4 Summary and Discussions

To summarize, we have studied possible superconductivity in carrier doped 5$d$ transition metal oxide Sr2IrO4 based on an effective $t_{2g}$ three-orbital Hubbard model on the square lattice with a large SOC. The linearized Eliashberg equation for the SC gap function was solved numerically with the RPA. We revealed that the $d_{x^{2}-y^{2}}$-wave pseudospin singlet pairing is expected to appear in the realistic parameter region, i.e., large SOC $\lambda$ and small Hund’s coupling $J/U$. The $d_{x^{2}-y^{2}}$-wave pseudospin singlet pairing is induced by the AF pseudospin fluctuations and is favored for large SOC especially in the electron doped region. Furthermore, we have found that the $s_{\pm}$-wave pseudospin singlet pairing is stabilized for small $\lambda$ and large $J/U$, and is more favored in the hole doped region. The $s_{\pm}$-wave pseudospin singlet pairing is induced by the inter band scattering arising from the intra-orbital AF spin fluctuations.

The $d_{x^{2}-y^{2}}$-wave pseudospin singlet pairing found in large $\lambda$ and small $J/U$ agrees with the previous studies [62, 64, 63]. In this region, the single-band nature with strong character of pseudospin $j_{\text{eff}}=1/2$ is enhanced especially in the electron doped region. Therefore, the $d_{x^{2}-y^{2}}$-wave pseudospin singlet pairing found here is an analogy to the $d_{x^{2}-y^{2}}$-wave spin $s=1/2$ singlet pairing in cuprate superconductors. The realistic parameter for Sr2IrO4 is in this parameter region because it is estimated that $\lambda\sim 0.37\ \text{eV}$ and $J/U\sim 0.07$ [28, 30, 62]. Therefore, the $d_{x^{2}-y^{2}}$-wave pseudospin singlet pairing in the electron doped region is the most promising candidate for possible superconductivity in Sr2IrO4.

The feature found here is similar to the case of cuprate superconductors from the viewpoint of purification of the $d_{x^{2}-y^{2}}$ orbital. H. Sakakibara et al. [84, 85, 86] have suggested that SC critical temperature $T_{\text{c}}$ of cuprate superconductors tends to increase with reducing the hybridization of the $d_{3z^{2}-r^{2}}$ orbital from the $d_{x^{2}-y^{2}}$ one, thus purifying the $d_{x^{2}-y^{2}}$ orbital, which mainly contributes to the $d_{x^{2}-y^{2}}$-wave spin singlet pairing. In the same way, we found here that the largest eigenvalue of the linearized Eliashberg equation increases with increasing the SOC, purifying the band with the strong character of $j_{\text{eff}}=1/2$, which contributes most to the $d_{x^{2}-y^{2}}$-wave pseudospin singlet pairing.

We found the $s_{\pm}$-wave pseudospin singlet pairing in the small $\lambda$ and large $J/U$ region, which is different from the results in the previous studies [64, 63]. Z. Y. Meng et al. [63] have suggested a $p$-wave pairing in a hole doped region. However, we found that a pseudospin singlet pairing is always favored over a pseudospin triplet pairing. On the other hand, Y. Yang et al. [64] have proposed a $s_{\pm}$-wave pairing using the functional renormalization group method, and they have also suggested that ferromagnetic (FM) fluctuations are responsible for the pairing mechanism. FM fluctuations may be enhanced when the Fermi pocket is small. However, the enhancement cannot be found in our RPA calculations. Instead, our calculations show that the $s_{\pm}$-wave pairing is induced by the inter-band scattering arising from the intra-orbital spin fluctuations, without strong FM fluctuations.

Finally, we briefly comment on $4d$ transition metal oxide Sr2RhO4. When the tight-binding parameters are estimated form the first-principles band structure calculations [9], it is found that the ratios of the tight-binding parameters of the three $t_{2g}$ orbitals, $(t_{1},t_{2},t_{3},t_{4},t_{5},\mu_{xy})/t_{1}$, for Sr2RhO4 are almost the same as those for Sr2IrO4, except that the ratio of the SOC to $t_{1}$, $\lambda/t_{1}\sim 0.5$, is smaller than that for Sr2IrO4 ($\lambda/t_{1}\sim 1.0$). Therefore, based on our calculations here, we expect that the superconductivity with $s_{\pm}$-wave pairing symmetry would be favorable for Sr2RhO4 especially in the hole doped region, provided that the Hund’s coupling is sufficiently large.

Acknowledgments

We would like to thank Toshihiro Sato, Kazuhiro Seki, Toshikaze Kariyado, and Hirofumi Sakakibara for useful discussions. The calculations were carried out using the supercomputer in Advanced Center for Computing and Communication centers of RIKEN and the HOKUSAI supercomputer at RIKEN. This work was supported in part by Grants-in-Aid for Scientific Research from JSPS (Projects No. 24740251 and No. 25287096) of Japan. T.S. acknowledges Simons Foundation for financial support (award no. 534160).