以前、以下の記事でデバイスの電気伝導現象をランダウアー理論を使って計算する方法をご紹介したと思います。
今回はこれを使って普通の金属(=常伝導状態)と超伝導体を接合したデバイスに応用してみたいと思います。
今回使う超伝導のモデル
今回は１次元キタエフ模型と呼ばれる模型を使いたいと思います。
1次元キタエフ模型
タイトバインディング模型として次のハミルトニアンで表されます。
H = -\frac{t}{2} \sum_j (a^\dagger_j a_{j+1} + a^\dagger_{j+1} a_j ) - \mu \sum_j a^\dagger_j a_j + \frac{\Delta}{2} \sum_j ( a^\dagger_j a^\dagger_{j+1} + a_{j+1} a_j )
$\Delta$が超伝導ギャップとなっていて、$\Delta=0$ならば普通のs軌道が最外殻電子になっている金属のモデルです。ただ、普通は考えられるスピンのパラメーターがありません。これは例えば強磁性体で最外殻のすべての電子のスピンが揃っている状況などに該当します。実験では強磁性体のFe原子を超伝導体のPb上に1列に並べるとPbの超伝導状態がFeにトンネル効果で染み出して１次元キタエフモ模型特有の現象が観測されています：
周期的な場合
周期的な場合は周期を$L$として離散フーリエ変換して
H = \frac{1}{L}\sum_{k,k'} \left[ -\frac{t}{2} \sum_j (a^\dagger_k a_{k'} e^{i[jk-(j+1)k']} + a^\dagger_{k'} a_k e^{-i[jk-(j+1)k']} ) - \mu \sum_j a^\dagger_k a_{k'} e^{ij(k-k')} + \frac{\Delta}{2} \sum_j ( a^\dagger_k a^\dagger_{k'} e^{i[jk+(j+1)k']} + a_{k'} a_k e^{-i[jk+(j+1)k']} ) \right] \\
= \sum_{k} \left[ \frac{1}{L}\sum_{k',j} e^{ij(k-k')} \left\{ -\frac{t}{2} (a^\dagger_k a_{k'} e^{ -ik' } + a^\dagger_{k'} a_k e^{ik'} ) - \mu a^\dagger_k a_{k'} \right\} + \frac{1}{L}\sum_{k',j} e^{ij(k+k')} \frac{\Delta}{2} ( a^\dagger_k a^\dagger_{k'} e^{ik'} + a_{k'} a_k e^{-ik'} ) \right] \\
\sum_{k} \left[ -\frac{t}{2} (a^\dagger_k a_{k} e^{ -ik } + a^\dagger_{k} a_k e^{ik} ) - \mu a^\dagger_k a_{k} + \frac{\Delta}{2} ( a^\dagger_k a^\dagger_{-k} e^{-ik} + a_{-k} a_k e^{ik} ) \right]
これは周期性から$k\rightarrow -k$としても同じなので
H = \sum_{k} \left[ \frac{1}{2} \left\{ -\frac{t}{2} (a^\dagger_k a_{k} e^{ -ik } + a^\dagger_{k} a_k e^{ik} ) - \mu a^\dagger_k a_{k} + \frac{\Delta}{2} ( a^\dagger_k a^\dagger_{-k} e^{-ik} + a_{-k} a_k e^{ik} ) \right\} + \frac{1}{2} \left\{ -\frac{t}{2} (a^\dagger_{-k} a_{-k} e^{ ik } + a^\dagger_{-k} a_{-k} e^{-ik} ) - \mu a^\dagger_{-k} a_{-k} + \frac{\Delta}{2} ( a^\dagger_{-k} a^\dagger_{k} e^{ik} + a_{k} a_{-k} e^{-ik} ) \right\} \right] \\
= \frac{1}{2} \sum_{k} \left[ -(t \cos{k} + \mu ) ( a^\dagger_k a_{k} + a^\dagger_{-k} a_{-k} ) + i\Delta\sin{k} ( a^\dagger_k a^\dagger_{-k} - a_{-k} a_k ) \right] \\
= \frac{1}{2}\sum_k ( a^\dagger_\mathbf{k} \quad a_{-\mathbf{k}} )
\left( \begin{matrix}
-t\cos{k}-\mu & i\Delta\sin{k} \\
-i\Delta\sin{k} & t\cos{k}+\mu
\end{matrix} \right)
\begin{pmatrix} a_\mathbf{k} \\ a^\dagger_{-\mathbf{k}}
\end{pmatrix} -\frac{1}{2}\sum_k( t\cos{k} + \mu )
定数項の分のズレはエネルギー基準点を変えて0に取り直して
H = \frac{1}{2}\sum_k ( a^\dagger_k \quad a_{-k} )
\left( \begin{matrix}
-t\cos{k}-\mu & i\Delta\sin{k} \\
-i\Delta\sin{k} & t\cos{k}+\mu
\end{matrix} \right)
\begin{pmatrix} a_k \\ a^\dagger_{-k}
\end{pmatrix}
となります。以下$ \mathbf{\psi}_k = ( a_k \quad a^\dagger_{-k} )^T $とします。ここで
\left( \begin{matrix}
-t\cos{k}-\mu & i\Delta\sin{k} \\
-i\Delta\sin{k} & t\cos{k}+\mu
\end{matrix} \right) = -\frac{1}{2} \begin{pmatrix} 0 \\ \Delta\sin{k} \\ t\cos{k}+\mu
\end{pmatrix} \cdot \vec{\sigma} \\
= -\frac{1}{2} \vec{R}(k) \cdot \vec{\sigma}
と表されます。計算を簡単にするためにアダマール変換
\mathcal{H} = \frac{\sigma_1+\sigma_3}{\sqrt{2}} = \mathcal{H}^\dagger
を作用させると
\mathcal{H} \begin{pmatrix} 0 \\ \Delta\sin{k} \\ t\cos{k}+\mu
\end{pmatrix} \cdot \vec{\sigma} \mathcal{H} = \begin{pmatrix} t\cos{k}+\mu \\ \Delta\sin{k} \\ 0
\end{pmatrix} \cdot \vec{\sigma}
となるので、
H = -\frac{1}{2}\sum_k ( \mathcal{H} \mathbf{\psi}_k )^\dagger
\left( \begin{matrix}
0 & t\cos{k}+\mu - i\Delta\sin{k} \\
t\cos{k}+\mu + i\Delta\sin{k} & 0
\end{matrix} \right)
\mathcal{H} \mathbf{\psi}_k
と書けます。
これを対角化すると
q(k) = t\cos{k}+\mu + i\Delta\sin{k} \\
\theta(k) = \arg[q(k)]
として
E_\pm = \pm \frac{1}{2} |q(k)| = \pm \frac{1}{2} \sqrt{ (t\cos{k}+\mu)^2 + \Delta^2\sin^2{k} } \\
\gamma(k) = \begin{pmatrix} \gamma_-(k) \\ \gamma_+(k)
\end{pmatrix} = \frac{1}{\sqrt{2}} \left( \begin{matrix}
1 & -e^{i\theta(k)} \\
e^{-i\theta(k)} & 1
\end{matrix} \right)
\mathcal{H} \begin{pmatrix} a_k \\ a^\dagger_{-k}
\end{pmatrix} \\
となります。$\gamma(k)$は超伝導準粒子(ボゴリューボフ準粒子)であり、$\theta(k)=-\theta(-k)$なので
\gamma_+^\dagger(k) = \gamma_-(-k)
が成り立ちます。また、チャーン数を計算すると、$\hat{q} = q/|q|$として
C = \frac{1}{2\pi i} \int_{-\pi}^\pi dk \frac{d\log{\hat{q}}}{dk}= \frac{1}{2\pi i} \oint_{\partial D(t,\mu)} \frac{d\hat{q}}{\hat{q}}
と表されます。ここで領域$D(t,\mu)$は$\hat{q}(k)$の軌跡が囲む領域であり、$t < |\mu|$のときは $0 \notin D$、$t > |\mu|$のときは$0 \in D - \partial D$、$t=|\mu|$のときは$0 \in \partial D$となります。
したがって、
C = \begin{cases}
1 \quad if \quad t \leq |\mu| \\
0 \quad if \quad t > |\mu|
\end{cases}
となり、$t \leq \mu$でエッジモードが存在します。これはバルクエッジ対応と言って、周期系で計算したチャーン数が$C\neq 0$のとき、開放系で端のサイトから指数的に減衰するような表面に局在する$E=0$の固有状態が存在することが知られています。これはマヨラナゼロモードと呼ばれており、対称性から常にマヨラナゼロモードは$2C$重縮退しています。
バンド図
t=1 \\
\Delta = 0.1t
の条件のもと$\mu= 0, 1, 2$の３種類のバンド図を描いてみます。
ちなみに$\Delta=0,\mu=0$だと次のようになります。
これを$\Delta\neq 0$の場合と比べてみると、上下のバンドが交差していて、交差している点の縮退が$\Delta$によって溶けていると分かります。
この交差が$|\mu|=t$を境になくなり、それがチャーン数の変化に反映されていると分かります。
再帰グリーン関数法
のときとの違いとして、今回は右のリードにBdGハミルトニアンを使って同じ計算手続きを行います。
扱いとして、左の端子を常伝導、右の端子を超伝導、メインのシステムを超伝導と常伝導の界面とします。
超伝導体の端子
H_0 = (a^\dagger_i,a_i) \left(\begin{matrix}
-\mu & 0 \\
0 & \mu \\
\end{matrix}\right) \left(\begin{matrix}
a_i \\ a^\dagger_i
\end{matrix}\right) \\
V = (a^\dagger_i,a_i) \left(\begin{matrix}
-t & \Delta \\
-\Delta & t \\
\end{matrix}\right) \left(\begin{matrix}
a_{i+1} \\ a^\dagger_{i+1}
\end{matrix}\right)
とすれば、BdGハミルトニアンは
H_R = (a^\dagger_1,a_1,\cdots,a^\dagger_N,a_N)\left(\begin{matrix}
H_0 & V & 0 & & & & \\
V^\dagger & H_0 & V & & & \\
& V^\dagger & H_0 & V & & \\
& & & & \ddots & V \\
& & & & V^\dagger & H_0
\end{matrix}\right) \left(\begin{matrix}
a_1 \\ a^\dagger_1 \\ \vdots \\ a_N \\ a^\dagger_N
\end{matrix}\right)
となり、計算で導かれる$i,j$サイト間のグリーン関数はゴルコフのグリーン関数
\hat{G}_{ij}(i\omega_n) = \left(\begin{matrix}
G_{ij}(i\omega_n) & F_{ij}(i\omega_n) \\
\tilde{F}_{ij}(i\omega_n) & \tilde{G}_{ij}(i\omega_n)
\end{matrix}\right)
です。ここで、ゴルコフのグリーン関数の各成分は$T$を時間順序積の記号として
G_{ij}(\tau) = - \langle T a_i(\tau) a^\dagger_j \rangle \\
\tilde{G}_{ij}(\tau) = - \langle T a^\dagger_i(\tau) a_j \rangle \\
F_{ij}(\tau) = - \langle T a_i(\tau) a_j \rangle \\
\tilde{F}_{ij}(\tau) = - \langle T a^\dagger_i(\tau) a^\dagger_j \rangle
であり、これらを松原周波数でフーリエ変換し、さらに実周波数に解析接続したものを使います。
メインのシステム
これは界面の絶縁膜のことなので、左右のリードの端にポテンシャル障壁のサイトをつなげばよいことになります。
よってポテンシャル障壁のサイトは
H_B = (a^\dagger_b,a_b)\left(\begin{matrix}
U_b -\mu & 0 \\
0 & -U_b + \mu
\end{matrix}\right) \left(\begin{matrix}
a_b \\ a^\dagger_b
\end{matrix}\right)
となり、全体のハミルトニアンは
H = \left(\begin{matrix}
H_L & \hat{t} & 0 \\
\hat{t}^\dagger & H_B & V_{sc} \\
0 & V_{sc}^\dagger & H_R
\end{matrix}\right)
と書けます。ここで
\hat{t} = \left(\begin{matrix}
-t & 0 \\
0 & t \\
\end{matrix}\right)
です。
コンダクタンスの計算
BdGハミルトニアンに対してMeir-Wingreen公式をそのまま適用すれば超伝導の場合も計算できます。超伝導の場合の公式はBlonder-Tinkham-Klapwijk (BTK)公式と言われていて久保公式で直接計算して導かれます。ただ、今回は１サイトしかなくて左のリードとの間の伝導を考えるので、少し工夫が必要です。左の端子の右端のサイト$L$と絶縁膜のサイト$B$の間の伝導度を計算すると、
L = \frac{ G^{(A)} - G^{(R)} }{2i} = \left(\begin{matrix}
L_{LL} & L_{LB} \\
L_{BL} & L_{BB}
\end{matrix}\right) \\
L_{LL} = \frac{1}{2i}(G^{(A)}_{LL} - G^{(R)}_{LL}) = - \mathrm{Im}(G^{(R)}_{LL})\\
L_{BB} = \frac{1}{2i}(G^{(A)}_{BB} - G^{(R)}_{BB}) = - \mathrm{Im}(G^{(R)}_{BB})\\
L_{LB} = \frac{1}{2i}(G^{(A)}_{LB} - G^{(R)}_{LB}) = - \mathrm{Im}(G^{(R)}_{LB})\\
L_{BL} = \frac{1}{2i}(G^{(A)}_{BL} - G^{(R)}_{BL}) = - \mathrm{Im}(G^{(R)}_{BL})
として、電流演算子
J = \frac{et}{i\hbar} \left(\begin{matrix}
0 & -\mathbf{1} \\
\mathbf{1} & 0
\end{matrix}\right)
を用いると、微分コンダクタンスは
\sigma_S(E) = \frac{\hbar}{2\pi}\mathrm{Tr}[LJLJ] \\
= \frac{e^2 t^2}{h} \mathrm{Tr}[( - L_{LR}L_{LR} - L_{RL}L_{RL} + L_{LL}L_{RR} + L_{RR}L_{LL} )]
となります。
実装
上付きの$L/R$は左(右)接続グリーン関数を表し、下付きの$L$は左側の端子の右端のサイト、$R$は右の端子の左端のサイト、$B$は絶縁膜のサイトを表します。
端子の計算は
のコードを使います。これによって$G^L_L$と$G^R_R$が求められます。上の記事では通常のグリーン関数で表記していますが、そもそも端子単独ではこれらのサイトは片側しかつながっていないので、左(右)接続グリーン関数と同じものになります。
障壁の部分のグリーン関数は1サイトしかないので、今回は手計算で求めることにします。
G^{R}_B(E) = [ E - H_B - V_{sc}^\dagger G^{R}_{R} V_{sc} ]^{-1} \\
G_L(E) = [ (G^{L}_L)^{-1} - t^\dagger G^{R}_{B} t ]^{-1} \\
G_B(E) = [ E - H_B - V_{sc}^\dagger G^{R}_{R} V_{sc} - t (G^{L}_L) t^\dagger ]^{-1} = [ (G^{R}_{B} )^{-1} - t (G^{L}_L) t^\dagger ]^{-1} \\
G_{LB}(E) = G_L t^\dagger G^R_{B} \\
G_{BL}(E) = G_B t G^L_L \\
超伝導端子の解析
まず、超伝導鎖の状態密度を解析します。
import numpy as np
from numpy import linalg as LA
from matplotlib import pyplot as plt
def G_lead(E,H0,V,n_lead,delta):
dagger = lambda M: np.conj(M.T)
dim = len(H0)
I = np.identity(dim).astype("complex")
Zero = np.zeros((dim,dim)).astype("complex")
G0 = LA.inv( complex(E,delta)*I - H0 )
A = G0.dot(dagger(V)); B = G0.dot(V)
# multiplication of A[0]*...*A[n]
AA = I.copy(); BB = I.copy()
Asum = Zero.copy(); Bsum = Zero.copy()
for i in range(n_lead):
Asum += BB.dot(A); Bsum += AA.dot(B);
AA = AA.dot(A); BB = BB.dot(B)
W = LA.inv( I - A.dot(B) - B.dot(A) )
A = W.dot(A.dot(A)); B = np.dot( W, B.dot(B) )
GL_inv = complex(E,delta)*I - H0 - V.dot(Asum)
GR_inv = complex(E,delta)*I - H0 - dagger(V).dot(Bsum)
return LA.inv(GL_inv), LA.inv(GR_inv)
以下のように実行します。$E=0,\pm2\Delta$で急激に変化するのでこれらの周辺だけ点を細かく取るようにエネルギーメッシュをとっています。
E_mesh = 2**7
dEmesh = 0.1
E0_arr = np.append(np.linspace(-dEmesh,dEmesh,2*E_mesh),0)
E2_arr = np.linspace(2.-dEmesh,2.+dEmesh,2*E_mesh)
E_arr = np.sort(np.append(np.append(np.append(np.linspace(-3,3,E_mesh),E0_arr),E2_arr),-E2_arr))
t = 1.
sc =0.01
mu = 0.*t
H0 = np.array([[-mu,0],[0,mu]])
V0 = np.array([[-t,0],[0,t]])
Vsc = np.array([[-t,sc],[-sc,t]])
G_arr = np.array([ G_lead(sc*E,H0,Vsc,2**6,1e-2*sc)[1] for E in E_arr ])
G00 = np.array([ G[0,0] for G in G_arr ])
G01 = np.array([ G[0,1] for G in G_arr ])
G10 = np.array([ G[1,0] for G in G_arr ])
G11 = np.array([ G[1,1] for G in G_arr ])
plt.plot(E_arr,-G00.imag)
plt.hlines(y=0,xmin=-5,xmax=5,color="black",ls=":")
plt.xlim([-3,3])
plt.ylim([-0.1,4])
plt.xlabel("$E/\Delta$")
plt.ylabel(r"$\rho$")
plt.show()
結果は次図のようになります。フーリエ変換して波数空間になおすと超伝導ギャップの大きさは$2\Delta$なので、$-2\Delta \leq E \leq 2\Delta$の範囲では状態がほぼ０です。しかし、アンドレーエフ束縛状態(エッジモード)が存在するので、$E=0$では状態が存在します。
伝導度の計算
次に伝導度を計算します。
import numpy as np
from numpy import linalg as LA
from matplotlib import pyplot as plt
def G_lead(E,H0,V,n_lead,delta):
dagger = lambda M: np.conj(M.T)
dim = len(H0)
I = np.identity(dim).astype("complex")
Zero = np.zeros((dim,dim)).astype("complex")
G0 = LA.inv( complex(E,delta)*I - H0 )
A = G0.dot(dagger(V)); B = G0.dot(V)
# multiplication of A[0]*...*A[n]
AA = I.copy(); BB = I.copy()
Asum = Zero.copy(); Bsum = Zero.copy()
for i in range(n_lead):
Asum += BB.dot(A); Bsum += AA.dot(B);
AA = AA.dot(A); BB = BB.dot(B)
W = LA.inv( I - A.dot(B) - B.dot(A) )
A = W.dot(A.dot(A)); B = np.dot( W, B.dot(B) )
GL_inv = complex(E,delta)*I - H0 - V.dot(Asum)
GR_inv = complex(E,delta)*I - H0 - dagger(V).dot(Bsum)
return LA.inv(GL_inv), LA.inv(GR_inv)
def G_Barrier( E, HB, VB,Vsc, GL_L_r, GR_R_r, delta ):
dagger = lambda M: np.conj(M.T)
Z = complex(E,delta)*np.identity(len(HB)).astype("complex")
self_EL = VB.dot(GL_L_r).dot(dagger(VB))
self_ER = dagger(Vsc).dot(GR_R_r).dot(Vsc)
GR_B = LA.inv( Z - HB - self_ER )
G_L = LA.inv( LA.inv(GL_L_r) - VB.dot(GR_B).dot(dagger(VB)) )
G_B = LA.inv( Z - HB - self_EL - self_ER )
G_LB = G_L.dot(dagger(VB)).dot(GR_B)
G_BL = G_B.dot(VB).dot(GL_L_r)
return G_L, G_LB, G_BL, G_B
def BTK(GL_a,GLB_a,GBL_a,GB_a):
FLL = np.imag(GL_a); FLB = np.imag(GLB_a);
FBL = np.imag(GBL_a); FBB = np.imag(GB_a);
T = - FLB.dot(FLB) - FBL.dot(FBL) + FLL.dot(FLL) + FBB.dot(FBB)
return np.trace(T).real
def calc_conductance(E_per_sc,t,mu,sc,Ub,len_lead,delta):
E = 2*sc*E_per_sc
H0 = np.array([[-mu,0],[0,mu]])
V0 = np.array([[-t,0],[0,t]])
Vsc = np.array([[-t,sc],[-sc,t]])
HB = np.array([[Ub-mu,0],[0,-Ub+mu]])
GL_L_a = G_lead(E,H0,V0,len_lead,-delta)[0]
GR_R_a = G_lead(E,H0,Vsc,len_lead,-delta)[1]
GL_a, GLB_a, GBL_a, GB_a = G_Barrier( E, HB, V0, Vsc, GL_L_a, GR_R_a, -delta )
return BTK(GL_a,GLB_a,GBL_a,GB_a)/2
トポロジーが非自明な場合
以下の条件で計算すると、下図のような結果になります。
E_mesh = 2**7
dEmesh = 0.1
E0_arr = np.append(np.linspace(-dEmesh,dEmesh,4*E_mesh),0)
E2_arr = np.linspace(2.-dEmesh,2.+dEmesh,4*E_mesh)
E_arr = np.sort(np.append(np.append(np.append(np.linspace(-3,3,E_mesh),E0_arr),E2_arr),-E2_arr))
t = 1.
sc = 0.01*t
delta = 1e-4*sc
Ub = 3.*t
mu = 0.*t
len_lead = 2**7
conductance_arr = np.array([ calc_conductance(E,t,mu,sc,Ub,len_lead,delta) for E in E_arr ])
conductance_no_barrier_arr = np.array([ calc_conductance(E,t,mu,sc,0,len_lead,delta) for E in E_arr ])
from matplotlib import style
from matplotlib.pyplot import hlines
plt.plot(E_arr,conductance_arr,label="Ub="+str(Ub)+"t")
plt.plot(E_arr,conductance_no_barrier_arr,label="Ub=0")
plt.xlabel("$E/2\Delta$")
plt.ylabel("$\sigma h/2e^2$")
plt.legend()
plt.xlim([-1.5,1.5])
plt.ylim([0,1.2])
plt.hlines(y=1.,xmin=-10,xmax=10,colors="black",ls=":")
plt.vlines(x=[0.,-1,1],ymin=0,ymax=10,colors="black",ls=":")
plt.show()
この計算は$-2t \leq \mu \leq 2t$の場合になっています。これはトポロジーが非自明なエネルギー領域になっています。(超伝導ギャップと同様に端数空間に直すと、ほっっピングは$2t$)になっています。)
まず、$U_b = 0$の場合を見てみます。伝導度はマヨラナゼロモード１個に対して$e^2/h$だけあります。１次元キタエフモデルではマヨラナゼロモードは2個あるので、この領域で$\sigma= 2e^2/h$で量子化しています。
$U_b=3t$の場合ではポテンシャル障壁によって透過率が減衰しています。しかし$E=0$ではコンダクタンスの量子化が保たれています。このエネルギー帯では$E-0$での量子化は必ず起こります。
トポロジカル量子臨界点の近傍
t = 1.
sc = 0.01*t
delta = 1e-4*sc
Ub = 3.*t
mu = 1.9*t
len_lead = 2**7
$|\mu| \sim 2t$の領域はトポロジーが自明/非自明な領域の境目になっていて接合の境界で完全反射が起こります。よってポテンシャルバリアがあると、超伝導ギャップのエネルギー帯での伝導度はほとんど$0$になります。(ギリギリトポロジーが非自明な領域なので、$E=0$でのマヨラナゼロモードの伝導は生き残っています。)
トポロジーが自明な場合
逆に上記の領域から外れていると、次の条件のようにほとんど伝導度は0です。(上下のバンド交差が解けてしまった状態になっている。)
t = 1.
sc = 0.01*t
delta = 1e-4*sc
Ub = 3.*t
mu = 2.1*t
len_lead = 2**7
参考文献
- 田中由喜夫 「超伝導接合の物理」 名古屋大学出版会
- https://arxiv.org/abs/cond-mat/0010440
- https://qiita.com/hiro949/items/804fba02078f33cfdb8d
- https://iopscience.iop.org/article/10.1088/1361-648X/ab8bf9#:~:text=The%20Kitaev%20chain%20is%20a%20one%20dimensional%20model,of%20standard%20fermionic%20operators%2C%20is%20%5B%201%2C%202%5D
- https://arxiv.org/pdf/1008.2026.pdf
- https://journals.aps.org/rmp/abstract/10.1103/RevModPhys.88.035005
- https://journals.jps.jp/doi/10.1143/JPSJ.71.2005
- https://journals.aps.org/prb/abstract/10.1103/PhysRevB.101.024509
- https://arxiv.org/abs/1906.01682
- https://arxiv.org/abs/cond-mat/0505164
- https://arxiv.org/abs/0705.0829
- https://journals.aps.org/prl/pdf/10.1103/PhysRevLett.47.882
- https://arxiv.org/abs/0909.4654v3
- https://journals.aps.org/prb/abstract/10.1103/PhysRevB.63.052512
- https://journals.aps.org/prresearch/abstract/10.1103/PhysRevResearch.2.043430