高次の電気伝導度の応答計算(理論パート2)

https://qiita.com/hiro949/items/c0124dc5ca1478d61b15

で紹介した計算の続きです。電場について高次の成分まで含んだ電気伝導度の一般的な表式を前回導出したので、今回は具体的な計算方法を紹介したいと思います。そして、1次、２次の項あたりを計算してみたいと思います。(長くなりすぎたので１次、２次項の計算はそれぞれ別の記事にします。)
やり方は

https://journals.aps.org/prb/abstract/10.1103/PhysRevB.99.045121

に倣っています。

# 電気伝導度の表式(復習)

電場の周波数成分$\mathbf{E}(\omega)$に対して電気伝導度$\sigma$を用いて電流の熱力学的平均値は

```math
\langle \hat{J}_\mathbf{E}^\mu(\omega) \rangle = \sum_{n=1}^{\infty} \int \prod_{i=1}^{n} d\omega_i \sigma^{\mu\alpha_1\cdots \alpha_n}(\omega;\omega_1,\cdots,\omega_n) \prod_{i=1}^{n} E_{\alpha_i}(\omega_i)
```

となります。ここで、(ギリシャ文字)$=x,y,z$であり、アインシュタインの縮約を使っています。このとき、電気伝導度は電場についての汎関数微分を使って

```math
\sigma^{\mu\alpha_1\cdots \alpha_n}(\omega;\omega_1,\cdots,\omega_n) = -i\omega \left( \frac{\delta^{n+1} \left\langle \hat{S}_c(E) \right\rangle_c }{\delta E_\mu(\omega) \delta E_{\alpha_1}(\omega_1) \cdots \delta E_{\alpha_n}(\omega_n) } \right)_{\mathbf{E} \rightarrow 0}
```

と表されます。式中の記号について、外場項は

```math
\hat{V}_E(t) = \sum_{n=1}^{\infty}\frac{1}{n!}\left[ \frac{e}{\hbar} \mathbf{A}(t) \cdot \hat{ \mathbf{D} } \right]^n \hat{H}_0
```

で、散乱行列

```math
\hat{S}_c(E) =  T\exp\left( -i \int dtV_E(t) \right)
```

で、$\langle * \rangle_c$はconnectedなファインマンダイヤグラムの総和です。

# ファインマンダイヤグラムと数式の対応

# 各部品で生じるファクター

まず、説明のために次の量を導入します。

```math
\hat{h}^{\alpha_1\cdots \alpha_n} = U ( \hat{D}^\alpha_1 \cdots \hat{D}^\alpha_n H_0 ) U^\dagger
```

### 外線と内線
今回は外場の項があるので、外のソースと電場を通じてエネルギー(=周波数)のやり取りをする波線(外線)とシステム内でエネルギーが保存されるの状態遷移を表す実線の矢印(内線)の２つで構成されます。

ちなみに内線は矢印が始点についている部分と終点の先っぽについている部分の２つがくっついてできています。前者が消滅演算子、後者が生成演算子に対応していて２つの始点と終点がくっついてグリーン関数

```math
G_{a}(\mathbf{k},\omega) = -i \int dt e^{-i\omega t}\langle T \gamma_a(\mathbf{k},t) \gamma_a^\dagger(\mathbf{k},0) \rangle
```

に対応しています。グリーン関数が生じる際に両社の添え字が一致しないと0になるので、両者の添え字が一致する成分のみが残ります。$\gamma$は１体ハミルトニアン$H_0$について対角化した消滅演算子です。

外線
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2445529/fecd6b3b-0055-bdc0-8487-f4916cd36e31.png)

内線
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2445529/ce555ef4-d9ba-e137-7c24-10519b331ddd.png)


### 電場とのカップリングの部分

また、以下のような生成、消滅演算子に対応した内線の間に外線がくっついているものも考えます。これは粒子が前の状態から遷移して電場とエネルギーのやり取りをして他の状態が生じる過程を表しています。生えている外線の本数は電場の何次の効果によるものかを表しています。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2445529/aa34e55a-fd00-a49d-62d0-8447433c40cb.png)

以上の図は内線の矢印から線が突き出ていますが、グリーン関数ではなく生成、消滅演算子の線なので注意してください。(この論文では線が矢尻から突き出たような描き方をしています。)

これらの項は電場が外線の頂点にくっついていない方の端が本来、電場と結合していて$\mathbf{E}\rightarrow 0$とすることから$E$で汎関数微分されたものだけが残ります。つまり外線１本につき汎関数微分を１回行っていることになります。

さらに、$\mathbf{E}(\omega)$は元々$\mathbf{A}(t)$で微分した部分だったので、微分した後に$-i\omega$を掛ける操作が加わっています。これを区別するために論文中ではノードを$\otimes$で表されています。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2445529/faf7a4c6-d963-64b5-da6f-48062761878b.png)

$i\omega$を掛けるという操作以外は他のノードと同じであり、外線が何本出ていても問題ありません。
ただし、そのうちの１本は周波数が$\omega$でなければいけません。

さらに、周波数$\omega$の外線は１本だけ特殊であり、注入する他と違ってエネルギーが出ていく出力部です。それゆえ、この出力部を$\otimes$ノードから生えている外線のうちのどれか一つを$E(\omega)$微分するときに区別が必要になります。つまり$\otimes$ノードから生える外線の本数分の自由度が生じます。

### 電場とのカップリングの部分で生じるファクター

前回の記事で$\langle \hat{S}_c \rangle_c$を$A^\mu(t)$で汎関数微分すると、$D^\mu$が一個出てきます。

```math
\frac{\delta \langle \hat{S}_c \rangle }{\delta A^\mu(t)} = \langle D^\mu \hat{W}_A H_0 \hat{S}_c \rangle \xrightarrow{E\rightarrow 0} \langle D^\mu H_0 \rangle = \langle h^\mu \rangle = \int d\mathbf{k} h^\mu_{ab}(\mathbf{k},t) \langle c_a^\dagger(\mathbf{k}) c_b(\mathbf{k}) \rangle
```

また、$E^\mu(\omega_i)$で微分すると、同様に$D^\mu$が出てきますが、それに加えて

```math
\frac{\delta A(t)}{\delta E(\omega_i)} = \frac{i e^{-i\omega_i t}}{2\pi\omega_i}
```

が出てきます。しかし、最後に$\int d\omega_i$をとるので、$i e^{-i\omega_i t}/2\pi$は最終的にデルタ関数になって消えて$i/\omega_i$が残ります。それと同時にこれによって生じたデルタ関数によって各ノードに入ってくる矢印の周波数と出ていく周波数の和が一致しなければならないという制約が発生します。(外線は$\omega$の線だけが出ていく線で、それ以外は入っていく線です。)また、$A(t)$で微分する部分も$E(\omega)$の形に直しましたが、$i/\omega$は式の先頭の$-i\omega$と互いに打ち消すので整合性が取れています。

最後に、$m$本の波線のギリシャ文字について和をとるときにラベルの入れ替えによる順列の重複が生じるので、$m!=\prod_{k=1}^m k$で割る必要があります。

最終的に出力部の選び方の自由度も合わせると、$\otimes$ノードから生える外線の本数が$l$本であるとき、各ダイヤグラムで重複を相殺するための重みは$\frac{l}{m!}$となります。



### ノードについて
ダイヤグラムには外線と内線が結合したノードがいくつかできますが、これらのポイントで内線の添え字や周波数について和(あるいは積分)をとる必要があります。一方、外線の周波数と添え字は入力/出力する量の成分や周波数を表しているので固定されています。

## Connectedなダイヤグラム

考えるダイヤグラムはConnected、つまりすべての内線が矢印の向きを合わせて繋がっていなければいけません。(超伝導のときは向きが合わない矢印同士をくっつけますが、今回は考えません。)また、トポロジー的に同じ図形は考えません。

すると、１次の項(=外線が２本)は次の２つのダイヤグラム

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2445529/62d83c7f-a59f-84ed-a950-3ff515a7adea.png)

です。以下の記事で実際に計算してみました。

https://qiita.com/hiro949/items/4d1de786a0fc4e059b7d

2次(=外線が3本)は6つ(うち2つは以下の４つの図で2番目と4番目の$\omega_1$と$\omega_2$を入れ替えたもの)です。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2445529/92fa0850-0bd2-3ebc-a690-9d4fb3b28d34.png)

以下の記事で実際に計算してみました。

https://qiita.com/hiro949/items/381af92af56e72b74f03

３次(＝外線が４本)は以下のダイヤグラムになります。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2445529/52d1845c-4699-dfaf-5325-79a117e8154c.png)

# グリーン関数に関する積分値

全てのダイヤグラムから出てくるグリーン関数に関する積分は

```math
\begin{align}
I^{(n)}_{a_1 a_2 \cdots a_n}(\omega_2,\omega_3,\cdots,\omega_n; \mathbf{k}) &= \int d\omega' \prod_{i=1}^{n} G_{a_i} \left( \omega'+ \sum_{2\leq j \leq i} \omega_j ; \mathbf{k} \right) \\
&= \int d\omega' G_{a_1}(\omega'; \mathbf{k} ) G_{a_2}(\omega'+\omega_2; \mathbf{k} ) \cdots G_{a_n}(\omega'+\omega_2+ \cdots \omega_n; \mathbf{k} )
\end{align}
```

という形になります。

## ウィック回転

論文では実時間周波数で表示していますが、今回のグリーン関数は熱力学平均による温度グリーン関数であるため、虚時間と松原和を使って周波数について記述する必要があります。(ウィック回転)

```math
\begin{align}
\omega^\text{フェルミオン}_l &= \frac{(2l+1)\pi}{\beta} \\
\omega^\text{ボゾン}_l &= \frac{2l\pi}{\beta} \\
G(\tau) &= \frac{1}{\beta}\sum_{i\omega_l} G(i\omega_l) = -\int_{-i\infty}^{+i\infty} \frac{d(i\omega)}{2\pi} G(i\omega)
\end{align}
```

その後に$i\omega_l \rightarrow \omega + i \delta$と解析接続することで実周波数の解が得られます。$\delta$（+0と書くこともある）は微小量で、$2\delta$などもすべて$\delta$と置き換えてしまいます。

ここで内線は電子の状態遷移、外線は電場(=光子)に対応しているので、内線をウィック回転しするときはフェルミオンの松原周波数を、外線はボゾンの松原周波数を使って置き換える必要があります。



$I^{(n)}$の場合、積分をとる$\omega'$がフェルミオン、引数になっている$\omega_2 \cdots \omega_n$がボゾンの場の周波数になっているので、それぞれ対応する松原周波数に置き換えて、

```math
\begin{align}
\tilde{I}^{(n)}_{a_1 a_2 \cdots a_n}(i\omega_2,i\omega_3,\cdots,i\omega_n; \mathbf{k}) &= \frac{1}{\beta} \sum_{i\omega'} \prod_{i=1}^{n} G_{a_i} \left( i\omega'+ \sum_{2\leq j \leq i} i\omega_j ; \mathbf{k} \right) \\
&= \frac{1}{\beta} \sum_{i\omega'} \prod_{i=1}^{n} \frac{1}{ i\omega'+ \sum_{2\leq j \leq i} i\omega_j - E_{a_i}(\mathbf{k}) }
\end{align}
```

という量を計算することになります。これは留数積分を使って次図の経路$C_{Im}$上で

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2445529/b5e8ce77-07d7-66a6-63d5-0c4c14483925.png)

```math
\tilde{I}^{(n)}_{a_1 a_2 \cdots a_n}(i\omega_2,i\omega_3,\cdots,i\omega_n) = -\oint \frac{dz}{2\pi i} \frac{1}{e^{\beta z} + 1} \prod_{i=1}^{n} \frac{1}{ z + \sum_{2\leq j \leq i} i\omega_j - E_{a_i}(\mathbf{k}) }
```

$C_{Im}$は虚軸周りを無限遠まで反時計回りにぐるっと囲んだものになっていますが、コーシーの積分定理から経路を無限遠を通る円周に膨らませるように変形することができます。無限遠の円周の部分は0に収束しますが、もともと$C_{Im}$の外側にあった虚軸上以外の極（たとえば次図の$z=\epsilon_a$）を時計回りに囲んだ閉曲線たちが残ります。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2445529/ff9b3b9a-58ac-2f9c-1889-5f644692093f.png)

再び留数定理を使うと、時計回りの経路なので$-1$がついて

```math
\begin{align}
\tilde{I}^{(n)}_{a_1 a_2 \cdots a_n}(i\omega_2,i\omega_3,\cdots,i\omega_n) &= \sum_i \frac{1}{e^{\beta (E_{a_i} - i\Omega_i)} + 1} \prod_{j \neq i, 2\leq j \leq n} \frac{1}{  E_{a_i}(\mathbf{k}) - E_{a_j}(\mathbf{k}) + i\Omega_j - i\Omega_i } \\
\Omega_i &= \sum_{2\leq j \leq i} i\omega_j
\end{align}
```
と表すことができます。そして$\omega_i$はボゾンの松原周波数なので、その和である$\Omega_i$もボゾンの松原周波数であり、

```math
e^{\beta \Omega_i} = 1
```

です。よって、

```math
\begin{align}
\tilde{I}^{(n)}_{a_1 a_2 \cdots a_n}(i\omega_2,i\omega_3,\cdots,i\omega_n) &= \sum_i \frac{1}{e^{\beta E_{a_i} } + 1} \prod_{j \neq i, 1 \leq j \leq n} \frac{1}{  E_{a_i}(\mathbf{k}) - E_{a_j}(\mathbf{k}) + i\Omega_j - i\Omega_i } \\ 
&= \sum_i f_{a_i}(\mathbf{k}) \prod_{j \neq i, 1\leq j \leq n} \frac{1}{  E_{a_i}(\mathbf{k}) - E_{a_j}(\mathbf{k}) + i\Omega_j - i\Omega_i }
\end{align}
```

となります。$f_{a}(\mathbf{k}) = f(E_a(\mathbf{k}))$はフェルミ分布です。$i\omega_i \rightarrow \omega_i$と実周波数に置き直して（$\Omega_i$も同様にします。）

```math
I^{(n)}_{a_1 a_2 \cdots a_n}(\omega_2,\omega_3,\cdots,\omega_n)
= \sum_i f_{a_i}(\mathbf{k}) \prod_{j \neq i, 1\leq j \leq n} \frac{1}{  E_{a_i}(\mathbf{k}) - E_{a_j}(\mathbf{k}) + \Omega_j - \Omega_i + i\delta }
```

となります。$2\delta$など$\delta=+0$の係数は$\delta$が微小量であることから$k\delta \rightarrow \delta$と置き直してしまいます。


例えば、

## n=1

```math
I^{(1)}_{a} = f_{a}
```

## n=2

```math
\begin{align}
\tilde{I}^{(2)}_{ab}(i\omega) &= \frac{f_a}{  E_a - E_b + i\omega } + \frac{f_b}{  E_b - i\omega - E_a  } \\
&=  \frac{f_a-f_b}{  E_a - E_b + i\omega } \rightarrow \frac{f_a-f_b}{  E_a - E_b + \omega + i\delta }
\end{align}
```

## n=3

$\omega_{12} = \omega_1 + \omega_2$として

```math
\begin{align}
\tilde{I}^{(3)}_{abc}(i\omega_1,i\omega_2) &= \frac{f_a}{ E_a + i\omega_1 - E_b} \frac{1}{ E_a + i\omega_{12} - E_c} + \frac{f_b}{ E_b - i\omega_1 - E_a} \frac{1}{ E_b - i\omega_1 + i\omega_{12} - E_c} +
\frac{f_c}{ E_c - i\omega_{12} - E_a} \frac{1}{ E_c - i\omega_{12} + i\omega_1 - E_b} \\
&= \frac{ f_a(E_b - E_c + i\omega_2) - f_b( E_a - E_c + i\omega_{12} ) + f_c( E_a - E_b + i\omega_1 ) }{ ( E_a - E_b  + i\omega_1 )(E_b - E_c + i\omega_2)( E_a - E_c + i\omega_{12} ) } \\
&= \frac{ (f_c-f_b)E_a + (f_a-f_c)E_b + (f_b-f_a)E_c + (f_c-f_b)i\omega_1 + (f_a-f_b)i\omega_2 }{ ( E_a - E_b  + i\omega_1 )(E_b - E_c + i\omega_2)( E_a - E_c + i\omega_{12} ) } \\
&= \frac{ (f_c-f_b)(E_a + i\omega_1) + (f_a-f_c)E_b + (f_b-f_a)(E_c-i\omega_2) }{ ( E_a - E_b  + i\omega_1 )(E_b - E_c + i\omega_2)( E_a - E_c + i\omega_{12} ) } \\
&= \frac{ (f_c-f_b)(E_a + i\omega_1) + (f_a-f_b+f_b-f_c)E_b + (f_b-f_a)(E_c-i\omega_2) }{ ( E_a - E_b  + i\omega_1 )(E_b - E_c + i\omega_2)( E_a - E_c + i\omega_{12} ) } \\
&= \frac{ (f_c-f_b)(E_a-E_b + i\omega_1) + (f_a-f_b)(E_b-E_c+i\omega_2) }{ ( E_a - E_b  + i\omega_1 )(E_b - E_c + i\omega_2)( E_a - E_c + i\omega_{12} ) } \\

&= \frac{ f_c - f_b  }{ ( E_b - E_c + i\omega_2 )( E_a - E_c + i\omega_1 + i\omega_2 ) } \\
&+  \frac{ f_a - f_b }{ ( E_a - E_b + i\omega_1 )( E_a - E_c + i\omega_1 + i\omega_2 ) } \\
&= \frac{ \tilde{I}^{(2)}_{ab}(i\omega_1;\mathbf{k}) - \tilde{I}^{(2)}_{bc}(i\omega_2;\mathbf{k}) }{ E_a - E_c + i\omega_1 + i\omega_2 } 
\end{align}
```

よって

```math
I^{(3)}_{abc}(\omega_1,\omega_2) = \frac{ I^{(2)}_{ab}(\omega_1;\mathbf{k}) - I^{(2)}_{bc}(\omega_2;\mathbf{k}) }{ E_a - E_c + \omega_1 + \omega_2 + i\delta } 
```
となります。

## 一般のn

$G_{a_1}(\omega)G_{a_n}(\omega + \sum_i\omega_i)$の項について部分分数分解すると、

```math
\frac{1}{ (i\omega' - E_{a_1} )( i\omega' + \sum_{i=2}^n i\omega_i - E_{a_n} ) }  = \left( \frac{1}{ i\omega' - E_{a_1} } - \frac{1}{ i\omega' + \sum_{i=2}^n i\omega_i - E_{a_n} } \right)　\frac{1}{ ( E_{a_1} - E_{a_n}) + \sum_{i=2}^n i\omega_i  }
```

となります。これを$\tilde{I}^{(n)}_{a_1 a_2 \cdots a_n}(i\omega_2,i\omega_3,\cdots,i\omega_n; \mathbf{k})$に代入すると、

```math
\begin{align}
\tilde{I}^{(n)}_{a_1 a_2 \cdots a_n}(i\omega_2,i\omega_3,\cdots,i\omega_n; \mathbf{k}) &= \frac{1}{\beta} \sum_{i\omega'} \prod_{i=1}^{n} \frac{1}{ i\omega'+ \sum_{2\leq j \leq i} i\omega_j - E_{a_i} } \\
&= \left( \frac{1}{\beta} \sum_{i\omega'} \prod_{i=1}^{n-1} \frac{1}{ i\omega'+ \sum_{2\leq j \leq i} i\omega_j - E_{a_i} } - \frac{1}{\beta} \sum_{i\omega'} \prod_{i=2}^{n} \frac{1}{ i\omega'+ \sum_{2\leq j \leq i} i\omega_j - E_{a_i} } \right)　\frac{1}{ ( E_{a_1} - E_{a_n}) + \sum_{i=2}^n i\omega_i  }
\end{align}
```

となります。$()$の中の第２項について$i\omega'' = i\omega' + i\omega_2$とするとこれは依然としてフェルミオンの松原周波数であり、

```math
\begin{align}
\tilde{I}^{(n)}_{a_1 a_2 \cdots a_n}(i\omega_2,i\omega_3,\cdots,i\omega_n; \mathbf{k}) &= \left( \frac{1}{\beta} \sum_{i\omega'} \prod_{i=1}^{n-1} \frac{1}{ i\omega'+ \sum_{2\leq j \leq i} i\omega_j - E_{a_i} } - \frac{1}{\beta} \sum_{i\omega''} \prod_{i=2}^{n} \frac{1}{ i\omega''+ \sum_{3\leq j \leq i} i\omega_j - E_{a_i} } \right)　\frac{1}{ ( E_{a_1} - E_{a_n}) + \sum_{i=2}^n i\omega_i  } \\
&= \frac{ \tilde{I}^{(n-1)}_{a_1 a_2 \cdots a_{n-1}}(i\omega_2,i\omega_3,\cdots,i\omega_{n-1} ) - \tilde{I}^{(n-1)}_{a_2 a_3 \cdots a_n}(i\omega_3,i\omega_4,\cdots,i\omega_n ) }{ E_{a_1} - E_{a_n} + \sum_{i=2}^n i\omega_i  } 
\end{align}
```

実数の周波数に解析接続して最終的に以下の漸化式を得ることができます。

```math
\begin{align}
I^{(1)}_{a_1} &= f_{a_1} \\
I^{(n)}_{a_1 a_2 \cdots a_n}(\omega_2,\omega_3,\cdots,\omega_n) &= \frac{ I^{(n-1)}_{a_1 a_2 \cdots a_{n-1}}(\omega_2,\omega_3,\cdots,\omega_{n-1} ) - I^{(n-1)}_{a_2 a_3 \cdots a_n}(\omega_3,\omega_4,\cdots,\omega_n ) }{ E_{a_1} - E_{a_n} + \sum_{i=2}^n \omega_i+ i\delta }
\end{align}
```
ここで$\delta=+0$は微小量なので$n\delta \rightarrow \delta$と$\delta$の係数は全部無視することにしています(この部分を実際の数値計算でどう取り扱うべきかは僕自身ちょっとよく分かっていません。)。
先に示した$n=2,3$の場合の結果も確かにこの漸化式に従っています。($\omega$の添え字は変更されていますが)


# 参考文献
- https://journals.aps.org/prb/abstract/10.1103/PhysRevB.99.045121
- http://www.kookai.pc.uec.ac.jp/kotaibutsuri/practical_guide_vol_2.pdf