RF谱PRA文献重复

Model 部分

Hamiltonian 为

\begin{align*} H =& \sum_{k,\sigma}\varepsilon_{k,\sigma} a_{k,\sigma}^{\dagger} a_{k,\sigma} + \sum_k (\varepsilon_{k,b}+\delta)b_k^{\dagger}b_k \\ & + \frac{\lambda}{\sqrt{\Omega}} \sum_{p,q}(b_{p+q}^{\dagger}a_{p,1}a_{q,2}+\mathrm{H.c.}) \end{align*}

前两项是原子和分子态的能量, 最后一项是态 1, 2 和分子态之间的耦合项. 其中

\begin{align*} \varepsilon_{k,\sigma} =& \frac{\hbar^2 k^2}{2m} - \mu_{\sigma} \\ \varepsilon_{k,b} =& \frac{\hbar^2 k^2}{4m} - \mu_1 -\mu_2 \\ \mu_1 =& \mu_3 = \mu ,\quad \mu_2 =0 \end{align*}

将耦合部分简记为

\begin{align*} \Lambda = \sum_{p,q}(b_{p+q}^{\dagger}a_{p,1}a_{q,2}+\mathrm{H.c.}) \end{align*}

微扰为

\begin{align*} V = \sum_k(a_{k,2}^{\dagger}a_{k,3} e^{-\mathrm{i} \omega't}+a_{k,3}^{\dagger}a_{k,2}e^{\mathrm{i}\omega't}) \end{align*}

Sum rules 部分

假设

假设系统的基态为 1 和 3 态的费米海

\begin{align*} |GS\rangle = |F \rangle = \prod_{k< k_F} a_{k,1}^{\dagger} a_{k,3}^{\dagger} |0\rangle \end{align*}

Transition rate form 3 to 2

\begin{align*} I(\omega) =& 2\pi \sum_f |\langle GS |V |f \rangle|^2 \delta(\omega+\mu -E_f +E_0) \\ =& 2 \sum_f |\langle GS |V |f \rangle|^2 \mathrm{Im}\frac{1}{\omega +\mu -E_f +E_0 -\mathrm{i}0^{ +}} \\ \propto &\mathrm{Im} \sum_f \langle GS |V \frac{1}{\omega +\mu -E_f +E_0 -\mathrm{i}0^{ +}} |f \rangle\langle f| V^{\dagger} |GS\rangle \\ =& \mathrm{Im} \langle GS |V \frac{1}{\omega +\mu -H +E_0 -\mathrm{i}0^{ +}} V^{\dagger} |GS\rangle \\ =& \mathrm{Im} \langle GS |V \frac{1}{\omega - \bar{H}} V^{\dagger} |GS\rangle \end{align*}

其中

\begin{align*} \bar{H} = H -E_0 -\mu \end{align*}

第二个等号用到了 Dirac Identity

\begin{align*} \frac{1}{x - x_0 \pm \mathrm{i}0^+} = \mathcal{P}\frac{1}{x-x_0} \mp \mathrm{i}\pi\delta(x-x_0) \end{align*}

零温公式推导

考虑中间态 \(|n\rangle\) 有两种

\begin{align*} |q \rangle =& a_{q,2}^{\dagger} a_{q,3}| F\rangle \quad(q< k_F)\\ |p,q \rangle =& b_{p+q}^{\dagger} a_{p,1} a_{q_3}|F \rangle \quad(p,q< k_F) \end{align*}

耦合作用到中间态上的结果为

\begin{align*} \Lambda |q \rangle =& \sum_p |p,q\rangle \\ \Lambda |p,q\rangle =& |q \rangle \end{align*}

作用是在两种中间态之间切换, 则 \(\bar{H}\) 在空间中的矩阵元为

\begin{align*} \langle q |\bar{H} |q \rangle =& \varepsilon_{q,2} -\varepsilon_{q,3} - \mu = 0 \\ \langle p,q|\bar{H} |p,q \rangle =& \delta + \varepsilon_{p+q,b} - \varepsilon_{p,1}-\varepsilon_{q,3}-\mu \\ \langle q| \bar{H}|p,q\rangle =& \langle p,q |\bar{H} | q\rangle = \lambda/\sqrt{\Omega} \end{align*}

Transition rate form 3 to 2 在上述空间中为

\begin{align*} I(\omega) \propto \mathrm{Im}\sum_{n,n'} \langle GS| V | n \rangle \langle n| \frac{1}{\omega-\bar{H}} | n'\rangle \langle n' V^{\dagger}|GS\rangle \end{align*}

将耦合 \(\lambda/\sqrt{\Omega}\Lambda\) 看作小量,即 \(\bar{H}= \bar{H_0}+\lambda/\sqrt{\Omega}\Lambda\) , 可做如下展开

\begin{align*} \frac{1}{\omega-\bar{H}} &=\frac{1}{\omega - \bar{H_0}-\lambda/\sqrt{\Omega}\Lambda}\\ =& \frac{1}{\omega - \bar{H_0}}\cdot \left( \frac{1}{1- \frac{\lambda}{\sqrt{\Omega}}\frac{1}{\omega-\bar{H_0}}\Lambda} \right) \\ =& \frac{1}{\omega - \bar{H_0}}\cdot \sum_m \left( \frac{\lambda}{\sqrt{\Omega}}\frac{1}{\omega-\bar{H_0}}\Lambda \right)^m \\ \end{align*}

而且将微扰 \(V\) 的表达式代入有

\begin{align*} \langle GS| V |n\rangle = \sum_q e^{\mathrm{i}\omega't}\delta_{n,q}\\ \langle n'| V |GS\rangle = \sum_{q'} e^{-\mathrm{i}\omega't}\delta_{n',q'}\\ \end{align*}

所以

\begin{align*} I(\omega) \propto& \mathrm{Im}\sum_{q,q'} \langle q| \frac{1}{\omega-\bar{H}} | q'\rangle \\ =&\mathrm{Im}\sum_{q,q'} \langle q|\frac{1}{\omega - \bar{H_0}}\cdot \sum_m \left( \frac{\lambda}{\sqrt{\Omega}}\frac{1}{\omega-\bar{H_0}}\Lambda \right)^m | q'\rangle \\ \end{align*}

由前面得出的耦合 \(\Lambda\) 作用在中间态上的结果可知, 耦合作用奇数次会由于正交性使结果为 \(0\) ( \(\bar{H_0}\) 为本征态), 所以 \(m\) 只能取偶数. 而且不同的 \(|q\rangle\) 也是正交的, 所以

\begin{align*} I(\omega) \propto& \mathrm{Im}\sum_q \langle q|\frac{1}{\omega - \bar{H_0}}\cdot \sum_m \left( \frac{\lambda}{\sqrt{\Omega}}\frac{1}{\omega-\bar{H_0}}\Lambda \right)^{2m} | q\rangle \\ =& \mathrm{Im}\sum_q \langle q|\frac{1}{\omega - \bar{H_0}}\cdot \sum_m \left( \frac{\lambda^2}{\Omega}\frac{1}{\omega-\bar{H_0}}\Lambda\frac{1}{\omega-\bar{H_0}}\Lambda \right)^m | q\rangle \\ =& \mathrm{Im}\sum_q \langle q|\frac{1}{\omega - \bar{H_0}}\cdot \frac{1}{1- \frac{\lambda^2}{\Omega}\frac{1}{\omega-\bar{H_0}}\Lambda\frac{1}{\omega-\bar{H_0}}\Lambda } | q\rangle \\ =& \mathrm{Im}\sum_q \langle q|\frac{1}{\frac{1}{\omega - \bar{H_0}}- \frac{\lambda^2}{\Omega}\Lambda\frac{1}{\omega-\bar{H_0}}\Lambda } | q\rangle \\ =& \mathrm{Im}\sum_{q< q_F}\frac{1}{\omega+ \mu +\varepsilon_{q,3}-\varepsilon_{q,2}-\lambda^2 \theta(q,\omega)} \end{align*}

其中

\begin{align*} \theta(q,\omega) = \frac{1}{\Omega} \sum_{p< k_F} \frac{1}{\omega+ \mu + \varepsilon_{q,3}+\varepsilon_{p,1} - \varepsilon_{p+q,b}-\delta} \end{align*}

零温数值

化简公式

\begin{align*} \mu + \varepsilon_{q,3} -\varepsilon_{q,2} = 0 \end{align*}

所以

\begin{align*} I(\omega) \propto & \mathrm{Im} \sum_{q< q_{F}} \frac{1}{\omega -\lambda^2 \theta(q,\omega)} \\ =& \sum_{q< k_F} \frac{\lambda^2\mathrm{Im}\theta(q,\omega)}{\left[ \omega - \lambda^2 \mathrm{Re}\theta(q,\omega) \right]^2 + \left[ \lambda^2 \mathrm{Im}\theta(q,\omega) \right]^2} \end{align*}

\begin{align*} \theta(q,\omega) =& \frac{1}{\Omega}\sum_{p< k_F} \frac{1}{\omega + \frac{\hbar^2}{4m}\left( p-q\right)^2-\delta} \end{align*}

无量纲化

对以下量进行无量纲化:

\begin{align*} \tilde{\omega} = \frac{\omega}{\varepsilon_F} ;\quad \tilde{\delta} = \frac{\delta}{\varepsilon_F} ;\quad \tilde{p} = \frac{p}{k_F} ;\quad \\ \tilde{q} = \frac{q}{k_F}; \quad \tilde{\lambda}^2 = \frac{\lambda^2}{\Omega \varepsilon_F^2} \end{align*}

原式就变为

\begin{align*} \theta(\tilde{q},\tilde{\omega}) =& \frac{1}{\Omega \varepsilon_F} \sum_{\tilde{p}< 1} \frac{1}{\tilde{\omega} + \frac{1}{2}(\tilde{p}-\tilde{q})^2 -\tilde{\delta} } \\ =&\frac{1}{\Omega\varepsilon_F}(B+ \mathrm{i}A) \end{align*}

其中

\begin{align*} A =& \pi \sum_{\tilde{p}< 1} \delta\left( \tilde{\omega}-\tilde{\delta} + \frac{1}{2}(\tilde{p}-\tilde{q})^2 \right) \\ B =& \sum_{\tilde{p} < 1} \mathcal{P}\frac{1}{\tilde{\omega}-\tilde{\delta} + \frac{1}{2}(\tilde{p}-\tilde{q})^2 } \end{align*}

Transition rate 就记为:

\begin{align*} \varepsilon_F I(\tilde{\omega}) =\sum_{\tilde{q}< 1} \frac{ \tilde{\lambda}^2A}{\left( \tilde{\omega} -\tilde{\lambda}^2 B\right)^2+\tilde{\lambda}^4 A^2} \end{align*}

数值处理

Dirac Delta 函数和主值分别用其极限形式:

\begin{align*} \delta(x) =& \frac{1}{\pi} \lim_{b\to 0}\frac{b}{b^2+x^2} \\ \mathcal{P}\frac{1}{x} =& \lim_{b\to 0}\frac{x}{b^2+x^2} \end{align*}

计算中 \(b=0.0001\) .

求和近似为从 \(0\) 到 \(1\) 的积分.

数值取值为

\begin{align*} \tilde{\lambda} = 1 \end{align*}

结果

code


import numpy as np
from matplotlib import pyplot as plt
from scipy import integrate

# 定义变量
l = 1 # lambd_tilde

#画图的精度
N = 1000

# 定义delta函数
def delta(x):
b = .0001
f = b/(b**2 + x**2)/np.pi
return f
# 定义主值函数
def pv(x):
b = 0.0001
f = x/(b*2 + x**2)
return f


# 定义最终想要得到的函数
def I(o,d):

def I(q):#定义没有积分的结果

def A(q):#定义theta的虚部
def A(p):
x = o + (p-q)**2/2 - d
A =np.pi * delta(x)
return A
(fres,err) = integrate.quad(A,0,1)
return fres

def B(q):#定义theta的实部
def B(p):
x = o + (p-q)**2/2 - d
B = pv(x)
return B
(fres,err) = integrate.quad(B,0,1)
return fres

I = l**2*A(q) / ( (o-l**2*B(q))**2 + l**4*A(q)**2 )
return I

(fres,err) = integrate.quad(I,0,1)#对其积分,即得到最终结果
return(fres)



#画图的横坐标omega从-2取到2
omega = np.linspace(-2,2,N)

C = 33 # 画C根线
#求出想要的函数在横坐标取值区间内的结果
I_omega = np.linspace(0,0,N)

for j in range (C):
d = j-(C-1)/2
d = d/2 #间距是分母分之一
for i in range(N):
I_omega[i-1] = I(omega[i-1],d)
print(d)
print(i-1)
plt.plot(omega,I_omega+100*j, label='d= %.2f' %(d) , color='gray')

#plt.legend()

plt.show()

fig

结果如图

图中从上到下, \(\tilde{\delta}\) 在区间 \([8,-8]\) 每间隔 \(0.5\) 取一个值.

上图的峰一些尖. 如果 Dirac delta 函数中的参数 \(b\) 取得更大一些, 图就会平滑一些.

下图是 \(b=0.5\) 时的结果. 为了直观,线的间距也由 \(100\) 改为 \(10\) .

eps矢量图: fig4a.eps

零温解析

推导回顾

原公式为:

\begin{align*} I(\omega) \propto \mathrm{Im}\sum_{q< q_F}\frac{1}{\omega+ \mu +\varepsilon_{q,3}-\varepsilon_{q,2}-\lambda^2 \theta(q,\omega)} \end{align*}

其中

\begin{align*} \theta(q,\omega) = \frac{1}{\Omega} \sum_{p< k_F} \frac{1}{\omega+ \mu + \varepsilon_{q,3}+\varepsilon_{p,1} - \varepsilon_{p+q,b}-\delta} \end{align*}

进行化简, 首先在这个系统中有如下关系

\begin{align*} \mu + \varepsilon_{q,3} -\varepsilon_{q,2} = 0 \end{align*}

所以原式可以化为:

\begin{align*} I(\omega) \propto \sum_{q< k_F} \frac{\lambda^2\mathrm{Im}\theta(q,\omega)}{\left[ \omega - \lambda^2 \mathrm{Re}\theta(q,\omega) \right]^2 + \left[ \lambda^2 \mathrm{Im}\theta(q,\omega) \right]^2} \end{align*}

其中

\begin{align*} \theta(q,\omega) =& \frac{1}{\Omega}\sum_{p< k_F} \frac{1}{\omega + \frac{\hbar^2}{4m}\left( p-q\right)^2-\delta} \end{align*}

对以下量进行无量纲化:

\begin{align*} \tilde{\omega} = \frac{\omega}{\varepsilon_F} ;\quad \tilde{\delta} = \frac{\delta}{\varepsilon_F} ;\quad \tilde{p} = \frac{p}{k_F} ;\quad \\ \tilde{q} = \frac{q}{k_F}; \quad \tilde{\lambda}^2 = \frac{\lambda^2}{\Omega \varepsilon_F^2} \end{align*}

原式就变为

\begin{align} \varepsilon_F I(\tilde{\omega}) =\sum_{\tilde{q}< 1} \frac{ \tilde{\lambda}^2A}{\left( \tilde{\omega} -\tilde{\lambda}^2 B\right)^2+\tilde{\lambda}^4 A^2} \end{align}

其中

\begin{align} A =& \pi \sum_{\tilde{p}< 1} \delta\left( \tilde{\omega}-\tilde{\delta} + \frac{1}{2}(\tilde{p}-\tilde{q})^2 \right) \\ B =& \sum_{\tilde{p} < 1} \mathcal{P}\frac{1}{\tilde{\omega}-\tilde{\delta} + \frac{1}{2}(\tilde{p}-\tilde{q})^2 } \end{align}

解析计算

求和化积分

\begin{align*} A \approx& \pi \frac{Vk^3_F}{(2\pi)^3} \int_{|\tilde{p}| < 1}\mathrm{d}^3\tilde{p} \cdot \delta\left( \tilde{\omega}-\tilde{\delta} + \frac{1}{2}(\tilde{p}-\tilde{q})^2 \right) =\pi \frac{Vk^3_F}{(2\pi)^3} a\\ B \approx& \frac{Vk^3_F}{(2\pi)^3} \int_{|\tilde{p}| < 1}\mathrm{d}^3\tilde{p} \cdot\mathcal{P}\frac{1}{\tilde{\omega}-\tilde{\delta} + \frac{1}{2}(\tilde{p}-\tilde{q})^2 } = \frac{Vk^3_F}{(2\pi)^3} b \end{align*}

解析的计算积分 \(a\) 和 \(b\)

\begin{align*} a =&2\pi \int_{-1}^{1}\mathrm{d}x \int_0^1\mathrm{d}\tilde{p} \cdot \delta\left( \tilde{\omega}-\tilde{\delta} + \frac{\tilde{p}^2}{2} +\frac{\tilde{q}^2}{2} -\tilde{p}\tilde{q}x \right) \tilde{p}^2\\ b =&2\pi\int_{-1}^{1}\mathrm{d}x \int_0^1\mathrm{d}\tilde{p} \cdot\mathcal{P}\frac{1}{\tilde{\omega}-\tilde{\delta} + \frac{\tilde{p}^2}{2} +\frac{\tilde{q}^2}{2} -\tilde{p}\tilde{q}x }\tilde{p}^2 \end{align*}

积分 \(a\) 是一个对 \(\delta\) 函数在 \([0,1]\) 区间内的积分. 所以要对 \(\delta\) 函数内的 部分分情况讨论.

如果 \(\tilde{\omega}-\tilde{\delta} + \frac{\tilde{p}^2}{2} +\frac{\tilde{q}^2}{2} -\tilde{p}\tilde{q}x = 0\) 没有根, 那么积分 \(a=0\) . 如果有根, 记两个根为 \(r_1, r_2\) , 那么积分 \(a\) 化为

\begin{align*} a =&2\pi\int_{-1}^{1}\mathrm{d}x \int_0^1\mathrm{d}\tilde{p} \cdot \delta\left[ (\tilde{p}-r_1)(\tilde{p}-r_2) \right] \tilde{p}^2 \\ =&2\pi\frac{1}{|r_1-r_2|} \int_{-1}^{1}\mathrm{d}x \int_0^1\mathrm{d}\tilde{p} \cdot \left[ \delta(\tilde{p}-r_1) +\delta(\tilde{p}-r_2) \right] \tilde{p}^2 \end{align*}

积分 \(b\) 是一个对主值的积分,也需要分情况讨论.

如果 \(\tilde{\omega}-\tilde{\delta} +\frac{\tilde{p}^2}{2} +\frac{\tilde{q}^2}{2} -\tilde{p}\tilde{q}x =0\) 没有根, 那么它就可以当做普通的积分

\begin{align*} b =&\int_0^1\mathrm{d}\tilde{p} \cdot \frac{1}{\tilde{\omega}-\tilde{\delta} + \frac{1}{2}(\tilde{p}-\tilde{q})^2 } \\ =&\frac{1}{\tilde{\omega}-\tilde{\delta}} \int_0^1 \mathrm{d}\tilde{p} \cdot \frac{1}{1+\left( \frac{\tilde{p}- \tilde{q}}{ \sqrt{2(\tilde{\omega} -\tilde{\delta})} } \right)^2} \\ =&\frac{1}{\tilde{\omega}-\tilde{\delta}} \cdot \left[ \arctan(t) \right]_{t==\frac{-\tilde{q}}{ \sqrt{2(\tilde{\omega} -\tilde{\delta})} }} ^{t=\frac{1-\tilde{q}}{ \sqrt{2(\tilde{\omega} -\tilde{\delta})} }} \end{align*}

如有两个根 \(r_1, r_2\) 那么

\begin{align*} b =&\int_0^1\mathrm{d}\tilde{p} \cdot \mathcal{P} \frac{1}{(\tilde{p}-r_1)(\tilde{p}-r_2)} \\ =&\frac{1}{r_1-r_2}\int_0^1\mathrm{d}\tilde{p} \cdot\left[ \mathcal{P} \frac{1}{\tilde{p}-r_1}-\mathcal{P}\frac{1}{\tilde{p}-r_2}\right] \end{align*}

两个根 \(r_1, r_2\) 是否在积分区间 \([0,1]\) 之间也会对积分 \(a, b\) 的结果有影响.

根据根的情况讨论积分结果

\begin{align} \tilde{\omega}-\tilde{\delta} +\frac{1}{2}(\tilde{p}-\tilde{q})^2 = 0 \end{align}

如果有根, 记 \(r_1 = \tilde{q} + \sqrt{2(\tilde{\delta}-\tilde{\omega})}, r_1 = \tilde{q} + \sqrt{2(\tilde{\delta}-\tilde{\omega})}\) . 当 \(\tilde{p}, \tilde{\omega}\) 取不同的值时, 根的情况如图

no root 时

\begin{align*} a =&0 \\ b = &\frac{1}{\tilde{\omega}-\tilde{\delta}} \cdot \left[ \arctan(t) \right]_{t==\frac{-\tilde{q}}{ \sqrt{2(\tilde{\omega} -\tilde{\delta})a} }} ^{t=\frac{1-\tilde{q}}{ \sqrt{2(\tilde{\omega} -\tilde{\delta})} }} \end{align*}

\(r_1, r_2 \in [0,1]\) 时

\begin{align*} a =& \frac{1}{ \sqrt{2(\tilde{\delta}-\tilde{\omega})} } \\ b =& \frac{1}{r_1-r_2} \ln \frac{ (1-r_1)r_2 }{ (1-r_2)r_1 } \end{align*}

\(r_1\in [0,1]\) 和 \(r_2 \in [0,1]\) 时

\begin{align*} a =& \frac{1}{2 \sqrt{2(\tilde{\delta}-\tilde{\omega})} } \\ b =& \frac{1}{r_1-r_2} \ln \frac{ -(1-r_1)r_2 }{ (1-r_2)r_1 } \end{align*}

\(r_1, r_2 \notin [0,1]\) 时

\begin{align*} a =&0 \\ b =& \frac{1}{r_1-r_2} \ln \frac{ (1-r_1)r_2 }{ (1-r_2)r_1 } \end{align*}

代回原式

\begin{align*} \varepsilon_F I(\tilde{\omega}) =&\sum_{\tilde{q}< 1} \frac{ \tilde{\lambda}^2A}{\left( \tilde{\omega} -\tilde{\lambda}^2 B\right)^2+\tilde{\lambda}^4 A^2} \\ \approx& \int_0^1 \mathrm{d}\tilde{q}\cdot \frac{ \tilde{\lambda}^2 \pi A}{\left( \frac{(2\pi)^3}{Vk_F}\tilde{\omega} -\tilde{\lambda}^2 B\right)^2+\tilde{\lambda}^4 \pi^2 A^2} \end{align*}

积分和结果

code

import numpy as np
from matplotlib import pyplot as plt
from scipy import integrate

l = 1 # 定义变量

def delta(x):
b = .01
f = b/(b**2 + x**2)/np.pi
return f

def I(o,d):
def I(q):
r1 = q + np.sqrt( 2*(d-o) )
r2 = q - np.sqrt( 2*(d-o) )

a1 = 0 # 无根
b1 = 1/(o-d) * ( \
np.arctan( (1-q)/np.sqrt( 2*(o-d) ) ) - \
np.arctan( ( -q)/np.sqrt( 2*(o-d) ) ))

a2 = 1/np.sqrt( 2*(d-o) ) # 根都在积分区间内
b2 = 1/(r1-r2)*np.log( ( (1-r1)*r2 )/ \
( (1-r2)*r1 ) )

a3 = .5/np.sqrt( 2*(d-o) ) # 一内一外
b3 = 1/(r1-r2)*np.log(-( (1-r1)*r2 )/ \
( (1-r2)*r1 ) )

a4 = 0 # 都在外
b4 = b2

if o>d:
A, B = a1, b1
case = 1
elif o>(d-q**2/2) and o>(d-(1-q)**2/2):
A, B = a2, b2
case = 2
elif o>(d-q**2/2) or o>(d-(1-q)**2/2):
A, B = a3, b3
case = 3
else:
A, B = a4, b4
case = 4

if case == 2 or case == 3:
I = l**2*A / ( (o-l**2*B)**2 + l**4*A**2 )
else:
I = delta(o - l**2*B)
return I
(fres, err) = integrate.quad(I, 0, 1)
return fres




N = 1000
#画图的横坐标omega从-2取到2
omega = np.linspace(-2,2,N)

C = 33 # 画C根线
#求出想要的函数在横坐标取值区间内的结果
I_omega = np.linspace(0,0,N)

for j in range (C):
d = j-(C-1)/2
d = d/2 #间距是分母分之一
for i in range(N):
I_omega[i-1] = I(omega[i-1],d)
print(d)
print(i-1)
plt.plot(omega,I_omega+j*10, label='d= %.2f' %(d) , color='gray')


plt.show()

fig

有限温公式推导

2 粒子的松原函数(有自能修正)

\begin{align*} G_2(k, \mathrm{i} \omega_2) =& \frac{1}{G_2^0(k,\mathrm{i}\omega_2)^{-1} -\lambda^2 \Sigma (k, \mathrm{i}\omega_2)} \\ =& \frac{1}{\mathrm{i}\omega_2 - \epsilon_{k,2} -\lambda^2 \Sigma (k, \mathrm{i}\omega_2)} \end{align*}

3 粒子的松原函数(无相互作用)

\begin{align*} G_3^0(k,\mathrm{i}\omega_3) = \frac{1}{\mathrm{i}\omega_3 - \epsilon_{k,3}} \end{align*}

代入总的松原函数

\begin{align*} R(k , \mathrm{i}\omega_n) =& \frac{1}{\beta}\sum_{\omega_2} G_3(k,\mathrm{i}(\omega_2-\omega_n)) G_2(k,\mathrm{i}\omega_2)\\ =& \frac{1}{\beta} \sum_{\omega_2} \frac{1}{\mathrm{i}(\omega_2 -\omega_n) - \epsilon_{k,3}}\cdot \frac{1}{\mathrm{i}\omega_2 - \epsilon_{k,2} -\lambda^2 \Sigma (k, \mathrm{i}\omega_2)} \\ =&\frac{1}{\beta} \sum_{\omega_2}\left( \frac{1}{\mathrm{i}(\omega_2 -\omega_n) - \epsilon_{k,3}} + \frac{-1}{\mathrm{i}\omega_2 - \epsilon_{k,2} -\lambda^2 \Sigma (k, \mathrm{i}\omega_2)} \right)\cdot \frac{1}{\mathrm{i}\omega_n +\epsilon_{k,3}-\epsilon_{k,2}-\lambda^2 \Sigma (k, \mathrm{i}\omega_2)} \end{align*}

接下来单独计算第一项

\begin{align*} \frac{1}{\beta} \sum_{\omega_2} \frac{1}{\mathrm{i}(\omega_2 -\omega_n) - \epsilon_{k,3}}\cdot \frac{1}{\mathrm{i}\omega_n +\epsilon_{k,3}-\epsilon_{k,2}-\lambda^2 \Sigma (k, \mathrm{i}\omega_2)} \end{align*}

利用环路积分计算此项. 积分

\begin{align*} \oint_{R\to\infty}\frac{1}{\mathrm{i}(\omega_2 -\omega_n) - \epsilon_{k,3}}\cdot \frac{1}{\mathrm{i}\omega_n +\epsilon_{k,3}-\epsilon_{k,2}-\lambda^2 \Sigma (k, \mathrm{i}\omega_2)} \cdot \frac{1}{e^{\mathrm{i}\beta\omega_2} +1}\mathrm{d}\omega_2 = 0 \end{align*}

积分的留数为:

第三项对应的留数

\begin{align*} R_3 = \frac{1}{\beta} \sum_{\omega_2} \frac{1}{\mathrm{i}(\omega_2 -\omega_n) - \epsilon_{k,3}}\cdot \frac{1}{\mathrm{i}\omega_n +\epsilon_{k,3}-\epsilon_{k,2}-\lambda^2 \Sigma (k, \mathrm{i}\omega_2)} \end{align*}

其中 \(\omega_2 = (2m+1)\pi/\beta\) . 就是要求的总的松原函数 \(R(k,\mathrm{i}\omega)\) .

第一项对应的留数

\begin{align*} R_1 =& \frac{1}{\mathrm{i}\omega_n +\epsilon_{k,3}-\epsilon_{k,2}-\lambda^2 \Sigma (k, \mathrm{i}\omega_n +\epsilon_{k,3})} \cdot \frac{1}{e^{\beta(\mathrm{i}\omega_n + \epsilon_{k,3})} +1} \\ =& \frac{ f(\epsilon_{k,3}) }{\mathrm{i}\omega_n +\epsilon_{k,3}-\epsilon_{k,2}-\lambda^2 \Sigma (k, \mathrm{i}\omega_n +\epsilon_{k,3})} \end{align*}

有限温及零温的最终重复结果

code

T=0.py

T=1.py

T=1.py

import numpy as np
from matplotlib import pyplot as plt
from scipy import integrate
from scipy.special import roots_legendre as leg
def gauquad(f,a,b,n = 50):
'''
定义 Gaussian quadrature 积分
函数 f 的积分区间为 [a,b]
取 n 个 Legendre 的根
def Gaussian quadrature integration
integrate function f from a to b
take n Legendre roots
'''

ft = lambda t: f( (b-a)*t/2 +(a+b)/2 ) * (b-a)/2
x, w = leg(n)
I = 0
for i in range(n):
I = I + w[i]*ft(x[i])
err = 0
return I,err
def twovarplt(f,ax,bx,ay,by,Nx=100,Ny=10,gap = 2):
'''
两变量画图
f(x,y): 要绘制的函数
变量 x 的取值区间为 [ax,bx], 取点的个数为 Nx 个
变量 y 的取值区间为 [ay,by], 画 Ny 条线, 每条线对应一个 y 值
线与线之间的间隔为 gap
two variables plot
plot function f(x,y)
x takes Nx points from 'ax' to 'bx'
each line corresponds a fixd y, y takes Ny points from 'ay' to 'by'
the gap between lines is 'gap'
'''

x = np.linspace(ax,bx,Nx)
y = np.linspace(ay,by,Ny)
fx = np.linspace(0,0,Nx)
for j in range(Ny):
for i in range(Nx):
print('第',j+1,'条线,','第',i+1,'个点,','共',Ny,'条线,','每条线',Nx,'个点.')
fx[i] = f(x[i],y[j])
plt.plot(x,fx+j*gap )
# plt.plot(x,x*0+y[j],color=(j/Ny,1-j/Ny,j/Ny))
plt.xlabel(r'$\tilde{\omega}$')
plt.ylabel(r'$I(\tilde{\omega})$')
plt.yticks([])
plt.savefig('./T = 2.jpeg',writer='imagemagick')
plt.show()

def f(p,m,T): # Fermi distribution function
f = 1/(np.exp((p**2-m)/T)+1)
return f
def theta(x): # Heaviside theta function
if x>=0:
f = 1
else:
f = 0
return f
def delta(x,b = 1e-2): # Dirac delta function
b = 1e-2
f = b/(b**2 + x**2)/np.pi
return f

T = 2 # Temperature
m = 1 # Chemical potential mu
up = 10 # Integration cut off of

def A(o,d,k,p): # the unintegrated Imaginary part of self-energy
x1 = o-d+p**2/2+k**2/2-p*k
x2 = o-d+p**2/2+k**2/2+p*k
if x1>0 or x2<0:
A = 0
else:
A = np.pi*p/k*f(p,m,T)
return A
def AA(o,d,k): # the integrated Imaginary part of self-energy
AA, err = gauquad(lambda p:A(o,d,k,p),0,up)
return AA
def s(o,d,k,p): # the unintegrated Real part of self-energy
s = f(p,m,T)*p/k*( np.log( np.abs(o-d + (k+p)**2/2) ) - np.log( np.abs(o-d + (k-p)**2/2) ) )
return s
def S(o,d,k): # the integrated Real part of self-energy
S,err = gauquad(lambda p:s(o,d,k,p),0,up)
return S

def I(o,d,k): # the unintegrated spectral function
aa = AA(o,d,k)
bb = S(o,d,k)
if aa == 0:
I = k**2*f(k,m,T)*delta(o-bb,b=1)
else:
I = k**2*f(k,m,T)*aa/((o-bb)**2+aa**2)
return I
def II(o,d): # the integrated spectral function
if d>o:
II,err = gauquad(lambda k:I(o,d,k),0,up,n = 100)
else:
II, err = integrate.quad(lambda k:I(o,d,k),0,up,epsabs=1.49e-2)
return II

twovarplt(II,-2,2,-8,8,Nx = 1000,Ny=33,gap=3)

#twovarplt(lambda k ,o:S(o,0,k),0,10,0.1,2,Nx = 100,Ny = 25,gap=0)

'''
def A(o,d,k,x):
g = o-d+k**2/2
r1 = k*x + np.sqrt( k**2*x**2 -2*g )
r2 = k*x - np.sqrt( k**2*x**2 -2*g )
A = theta(-g)*theta(1+2*g-2*k*x)*r1**2*f(r1,m,T)
A = A+ theta(g)*theta(k*x-np.sqrt(2*g))*theta(x) * ( theta(1+2*g-2*k*x)*r1**2*f(r1,m,T) + r2**2*f(r2,m,T))
A = A/(r1-r2)
return A
def AA(o,d,k):
AA, err = gauquad(lambda x:A(o,d,k,x),-1,1)
AA = np.nan_to_num(AA)
AA = np.pi*AA
return AA
'''

'''
n = 500
k = np.linspace(1e-5,up,n)
for i in range(n):
dd = np.abs(o - S(o,d,k[i]))
# print(dd)
if dd<1e-1:
ko = k[i]
break
else:
ko = 0
II = ko**2*f(ko,m,T)
return II
'''

Supplementary


import numpy as np
from matplotlib import pyplot as plt

delta = .5
N = 100
o = np.linspace(-2, 2, N)
q = np.linspace(-2, 2, N)
fig, ax = plt.subplots()
#line = ax.plot(o,q)
#ax3 = fig.add_axes([0.1, 0.1, 0.2, 0.2])

#ax.plot(o,o*0+1)
#ax.plot(o,o*0)
ax.plot(o*0+delta,o, color = 'green')
ax.plot(o*0+delta-.5,o, color = 'gray', linestyle = '--')
ax.plot(delta - q**2/2, q, 'r', label = r'$\~\omega = \~\delta - \frac{\~q^2}{2}$')
ax.plot(delta - (q-1)**2/2, q, color = 'blue', label = r'$\~\omega = \~\delta - \frac{(\~q-1)^2}{2}$')

ax.text(-1.2, .5, r'$r_1,r_2\notin [0,1]$', {'fontsize':20})
ax.text(.05, .95, r'$r_2\in [0,1]$', {'fontsize':20})
ax.text(.07, .05, r'$r_1\in [0,1]$', {'fontsize':20})
ax.text(1, .5, r'no root', {'fontsize':20})

ax.arrow(1, .25, -.55, .25, width = .01)
ax.text(1, .25, r'$r_1,r_2\in [0,1]$', {'fontsize':20})

#ax.set_xlim(0,4*np.pi)
ax.set_xlim(-2,2)
ax.set_ylim(0,1)
ax.set_xticks([-2, delta-.5, delta, 2])
ax.set_yticks([0,1],)
ax.set_xticklabels(['$-2$', '$\~\delta-1/2$', '$\~\delta$', '$2$'],{'fontsize':20})
ax.set_yticklabels([0,1],{'fontsize':20})
ax.set_xlabel('$\~\omega$',{'fontsize':20})
ax.set_ylabel('$\~q$',{'fontsize':20})
ax.legend(fontsize = 20)
plt.show()

参考文献

Junjun Xu, Qiang Gu, and Erich J. Mueller Phys. Rev. A 88, 023604 (2013)