在这里先整理一个最简答的问题，再知道了体系的哈密顿量后，怎么求解系统的格林函数并求得态密度。

数值角度考虑

在这里以一个陈绝缘体的哈密顿量为例来演示具体每一步是如何进行的。 \(H(k_x,k_y)=t(\sin(k_x)s_x+\sin(k_y)s_y)+(1-\cos(k_x)-\cos(k_y))s_z\) 这里的$s_i$代表Pauli矩阵，这是一个2*2的Hamiltonian，通常格林函数为$G(k_x,k_y,\omega)=\frac{1}{i\omega-\hat{H}}$,虽然这里在写法上是个分数的形式，但注意这里分母上的\hat{H}其实是个算符，这里请理解这一点，先不要将它就当作是这里的哈密顿量。接下来我们要计算格林函数，这里是要求倒数，那么对应到矩阵形式就是要求解$i\omega-H$这个矩阵的逆，这个时候就已经把算符$\hat{H}$写成具体的矩阵形式了，这个关系要捋清楚。

那么矩阵求逆就没什么好讲的了，对于这个2*2的矩阵解析和数值都是可以很轻松的求解出来。态密度的话就是把$k_x$和$k_y$对应下的格林函数积分即可 \(A(\omega)=\int_0^{2\pi} G(k_x,k_y,\omega)d\mathbf{k}\) 下面附上Mathematica的代码作为示例，可以进行参考

t=1.0;
i=PauliMatrix[0];x=PauliMatrix[1];y=PauliMatrix[2];z=PauliMatrix[3];
H[kx_,ky_]:=t(Sin[kx]x+Sin[ky]y)+(1-Cos[kx]-Cos[ky])z;
G2[kx_,ky_,\[Omega]_]:=Inverse[\[Omega] i-H[kx,ky]]
spec2[kx_,ky_,\[Omega]_]:=-Im@Tr@G2[kx,ky,\[Omega]]
dos2[\[Omega]_]:=1/(2\[Pi]) Sum[spec2[kx,ky,\[Omega]+I*0.01]*0.1^2,{kx,-Pi,Pi,0.1},{ky,-Pi,Pi,0.1}]

解析角度考虑

通过上面的内容已经可以知道$G=\frac{1}{i\omega-\hat{H}}$，在实际的计算中H一般会是个矩阵形式，当它在分母上的时候求格林函数就要对矩阵求逆，淡然了求逆是很不快捷的，如果能把整个矩阵放到分子，那么所有问题就变成了矩阵相乘的计算。将分母当作代数表达式，那么利用平方和公式把在分母的问题转换到分子，变成矩阵相乘的问题 \(G=\frac{1}{i\omega-H}=\frac{i\omega+H}{(i\omega)^2-H^2}\) 这时候分子上的H仍然是个矩阵形式，但是分母上的H需要用矩阵的本征值来代替，所以这个方式只是将矩阵求逆变成了矩阵本征值求解。通常哈密顿量是用Pauli矩阵写出来的，那么它的本征值是比较容易求解的。这里的H是个2*2的矩阵，本征值有两个，一正一负，平方之后是相同的。在这里只是考虑了大小为2的矩阵问题，对于更高维度的矩阵，两种方法是等价的，但是涉及到求态密度的时候，可能在这里要考虑到底是哪个本征值来贡献态密度了，这个问题在之后会单独讨论。下面上代码演示

t=1.0;
i=PauliMatrix[0];x=PauliMatrix[1];y=PauliMatrix[2];z=PauliMatrix[3];
H[kx_,ky_]:=t(Sin[kx]x+Sin[ky]y)+(1-Cos[kx]-Cos[ky])z;
val=Eigenvalues[H[kx,ky]];
G[kx_,ky_,\[Omega]_]:=Evaluate[(\[Omega]+H[kx,ky])/(\[Omega]^2-val[[1]]^2)]
spec[kx_,ky_,\[Omega]_]:=-Im@Tr@G[kx,ky,\[Omega]]
dos[\[Omega]_]:=1/(2\[Pi]) Sum[spec[kx,ky,\[Omega]+I*0.01]*0.1^2,{kx,-Pi,Pi,0.1},{ky,-Pi,Pi,0.1}]

通过这两种方式计算得到的态密度结果是相同的。解析角度的考虑更加适合在解析推导时可以更加方便的计算，当然了数值的考虑也同可可以解析的计算，只不过解析求解矩阵的逆矩阵时可能过程相对于求矩阵本征值要稍微繁琐一些，这里所说的矩阵是哈密顿量对应的矩阵，并不是一般任意的矩阵。

公众号

相关内容均会在公众号进行同步，若对该Blog感兴趣，欢迎关注微信公众号。