Ayaya Notes

lecture 2 exact diagonalization

实际上做的是求解定态的 Schrodinger 方程

H∣Ψ⟩=E∣Ψ⟩H\ket{\Psi}=E\ket{\Psi}

Why ED? 这种算法对于选定的有限系统是数值精确的 (仍受浮点精度和迭代收敛误差限制), 可以作为其他算法的 benchmark. 如果得到完整的本征值和本征态, 原则上可以计算这个有限系统的全部静态/动力学性质; 只求低能谱时则不能直接获得全部信息。

但这种算法的根本困难来自于 Hilbert space 的 “exponential wall” , 计算复杂度随自由度指数增长, 对于复杂的 Hilbert space 计算是困难, 甚至是不可做到的 ED 通常只被限制到很小的系统

一种做法是做 finite size scaling, 通过逐渐计算自由度增大的相同系统, 来外推更大系统的结果

各种应用

计算尺度的例子: spin - 1/2 的某些模型可做到 N=48N=48 的低能谱计算, 这依赖于模型、对称性和计算资源, 并不是通用上限, 也不是对整个 Hilbert 空间做全谱对角化。常见多体模型通常只能处理几十个格点。

CODE

ED 的代码结构:

on the fly 即指比起考虑存下来, 不如重新计算一遍 一种以时间换空间的方式

我们以 spin - 1/2 作为例子 考虑 4 个 spin, 一共 16 个 basis

一种储存方法是把每个态对应到一个二进制数, 再映射到十进制数, 如

∣↑↑↑↑⟩↦∣1111⟩2↦∣15⟩\ket{\uparrow\uparrow\uparrow\uparrow}\mapsto\ket{1111}_2\mapsto\ket{15}

bit operations

对于逐 bit 二进制数操作, 如 AND , OR, XOR, NOR

也有 bit shift, 移动数位

numpy 提供了进制数的操作抽象

由于我们用二进制数表示了系统的基底, 也可以将算符抽象为对 bit 的操作

例如: 设置第 n 个位置的 spin 为 1 (bit 从 0 开始编号)

def SetBit(i,n):
    #输入态 i, 指定设置bit 位 n
    #输出设置后的态
    return i | (1 << n)

以及各种 bit 上的操作, 只要把基底的 spin 上下与 bit 位的 0/1 对应即可

现在我们考虑我们的系统: 取 J=1J=1, ℏ=1\hbar=1, Hamiltonian 为

H=∑i=0N−1Si⋅Si+1=∑i=0N−1(SixSi+1x+SiySi+1y+SizSi+1z)=∑i=0N−1[(Si+Si+1−+Si−Si+1+)/2+SizSi+1z]\begin{aligned}\mathrm{H}&=\sum_{i=0}^{N-1}\mathbf{S}_{i}\cdot\mathbf{S}_{i+1}=\sum_{i=0}^{N-1}\left(S_{i}^{x}S_{i+1}^{x}+S_{i}^{y}S_{i+1}^{y}+S_{i}^{z}S_{i+1}^{z}\right)\\&=\sum_{i=0}^{N-1}\left[\left(S_i^+S_{i+1}^-+S_i^-S_{i+1}^+\right)/2+S_i^zS_{i+1}^z\right]\end{aligned}

设置周期性边界条件

SN=S0\mathbf{S}_N=\mathbf{S}_0

Diagonal term

在上述 SzS^z 乘积基底中, 对应 bit 位, 可以定义为

⟨col∣SizSi+1z∣col⟩=⟨Siz⟩⟨Si+1z⟩⟨Siz⟩=ReadBit⁡(col,i)−0.5\begin{aligned} \langle\mathrm{col}|S_i^zS_{i+1}^z|\mathrm{col}\rangle&=\langle S_i^z\rangle\langle S_{i+1}^z\rangle\\ \langle S_i^z\rangle&=\operatorname{ReadBit}(\mathrm{col},i)-0.5 \end{aligned}

Off-diagonal term

容易得到

i+1⟨0∣i⟨1∣(Si+Si+1−+Si−Si+1+)/2∣0⟩i∣1⟩i+1=0.5i+1⟨1∣i⟨0∣(Si+Si+1−+Si−Si+1+)/2∣1⟩i∣0⟩i+1=0.5\begin{aligned}&_{i+1}\langle0\mid_i\langle1\mid\left(S_i^+S_{i+1}^-+S_i^-S_{i+1}^+\right)/2\mid0\rangle_i\mid1\rangle_{i+1}=0.5\\&_{i+1}\langle1\mid_i\langle0\mid\left(S_i^+S_{i+1}^-+S_i^-S_{i+1}^+\right)/2\mid1\rangle_i\mid0\rangle_{i+1}=0.5\end{aligned}

采用稀疏矩阵的储存方式: row, column, value

我们的 sample code 为 (这里改用开边界, 所以不包含首尾相连的键)

from scipy import sparse


def ReadBit(i, n):
    return (i >> n) & 1


def FlipBit(i, n):
    return i ^ (1 << n)


def GetHopList(Ns): #这里是开边界版本
    HopList = []
    for i in range(Ns-1):
        HopList.append([i, i+1])
    return HopList


Ns = 4
NL = 1 << Ns
HopList = GetHopList(Ns)

HI = []
HJ = []
HV = []

for i0 in range(NL):
    for ih in range(len(HopList)):
        Pos0 = HopList[ih][0]
        Pos1 = HopList[ih][1]

        if ReadBit(i0, Pos0) != ReadBit(i0, Pos1):
            i1 = FlipBit(i0, Pos0)
            i1 = FlipBit(i1, Pos1)
            HI.append(i1)
            HJ.append(i0)
            HV.append(0.5)

        HI.append(i0)
        HJ.append(i0)
        HV.append((ReadBit(i0, Pos0)-0.5)*(ReadBit(i0, Pos1)-0.5))

Hamr = sparse.coo_matrix((HV, (HI, HJ)), shape=(NL, NL)).tocsc()

Symmetries

这个 Hilbert 空间如果有对称性, 与 Hamiltonian 对易的话, Hamiltonian 矩阵块在对应对称性的基底下会有更小的分块

在这个问题里, 具有 U(1)U(1) 对称性: 总 SzS^z 守恒, 等价于向上 spin 数 (bit 为 1 的个数) 守恒。

不断选择各种对称性, 就能利用 Hamiltonian 分块的性质压缩计算量

Implement U(1) symmetry: 即按照向上 spin 数来划分 basis, 则 Hamiltonian 自动按对应扇区分块对角。固定 N↑N_\uparrow 的子空间维数为 (NN↑)\binom{N}{N_\uparrow}。

Basis look up procedure: 把原始的基底编号映射到新的基底, 我们需要做到尽可能快

一种著名的方法为 Lin Tables 这是尝试解决在固定粒子数子空间的 lookup 方法: 考虑 6 个 spin - 1/2, 将构型分为两个部分 A,B 然后为 A,B 分别建立 lookup 表, 在选定的排序下, 将完整构型在固定粒子数子空间中的 index 写成两张表的查表值之和。

这里的查表值是专门构造的偏移量, 不能直接用两部分各自的普通编号相加; 得到的是子空间中的连续 index, 一般不等于原构型的二进制整数编号。

Fermionic systems

考虑 Hubbard model

H=−t∑⟨i,j⟩,σ(Ciσ†Cjσ+H.c.)+U∑ini↑ni↓ni,σ=Ci,σ†Ci,σH=-t\sum_{\braket{i,j},\sigma}\left(C_{i\sigma}^\dagger C_{j\sigma}+\mathrm{H.c.}\right)+ U\sum_{i} n_{i\uparrow}n_{i\downarrow}\quad n_{i,\sigma}=C_{i,\sigma}^\dagger C_{i,\sigma}

这里 tt 为实数, ⟨i,j⟩\braket{i,j} 对每条无向近邻键只计一次, 因而显式加上 Hermitian conjugate。

约定反对易关系

{Ci,σ,Ci′,σ′†}=δii′δσσ′,{Ci,σ†,Ci′,σ′†}={Ci,σ,Ci′,σ′}=0\{C_{i,\sigma},C_{i^{\prime},\sigma^{\prime}}^\dagger\}=\delta_{ii^{\prime}}\delta_{\sigma\sigma^{\prime}},\{C_{i,\sigma}^\dagger,C_{i^{\prime},\sigma^{\prime}}^\dagger\}=\{C_{i,\sigma},C_{i^{\prime},\sigma^{\prime}}\}=0

并约定态的 convention , 产生算符的顺序:

C0,↑†C2,↑†C3,↑†C1,↓†C3,↓†∣0⟩C_{0,\uparrow}^\dagger C_{2,\uparrow}^\dagger C_{3,\uparrow}^\dagger C_{1,\downarrow}^\dagger C_{3,\downarrow}^\dagger|0\rangle

这里使用: spin 向上在前面, 同一方向 spin 内部从小到大排

在这里, 我们使用两组二进制数来存

按照我们的 convention, 有如下结果

Ci,↑†[⋅⋅⋅0i,↑⋅⋅⋅∣0⟩]=(−1)∑k=0i−1nk,↑[⋅⋅⋅Ci,↑†⋅⋅⋅∣0⟩]Ci,↑[⋅⋅⋅Ci,↑†⋅⋅⋅∣0⟩]=(−1)∑k=0i−1nk,↑[...0i,↑⋅⋅⋅∣0⟩]Ci,↓†[⋅⋅⋅0i,↓⋅⋅⋅∣0⟩]=(−1)∑k=0i−1nk,↓+∑k=0N−1nk,↑[⋅⋅⋅Ci,↓†⋅⋅⋅∣0⟩]Ci,↓[⋅⋅⋅Ci,↓†⋅⋅⋅∣0⟩]=(−1)∑k=0i−1nk,↓+∑k=0N−1nk,↑[⋅⋅⋅0i,↓⋅⋅⋅∣0⟩]\begin{gathered}C_{i,\uparrow}^\dagger\left[\cdot\cdot\cdot0_{i,\uparrow}\cdot\cdot\cdot|0\rangle\right]=(-1)^{\sum_{k=0}^{i-1}n_{k,\uparrow}}\left[\cdot\cdot\cdot C_{i,\uparrow}^\dagger\cdot\cdot\cdot|0\rangle\right]\\C_{i,\uparrow}\left[\cdotp\cdotp\cdotp C_{i,\uparrow}^{\dagger}\cdotp\cdotp\cdotp|0\rangle\right]=(-1)^{\sum_{k=0}^{i-1}n_{k,\uparrow}}\left[...0_{i,\uparrow}\cdotp\cdotp\cdotp|0\rangle\right]\\C_{i,\downarrow}^\dagger\left[\cdot\cdot\cdot0_{i,\downarrow}\cdot\cdot\cdot|0\rangle\right]=(-1)^{\sum_{k=0}^{i-1}n_{k,\downarrow}+\sum_{k=0}^{N-1}n_{k,\uparrow}}\left[\cdot\cdot\cdot C_{i,\downarrow}^\dagger\cdot\cdot\cdot|0\rangle\right]\\C_{i,\downarrow}\left[\cdot\cdot\cdot C_{i,\downarrow}^\dagger\cdot\cdot\cdot|0\rangle\right]=(-1)^{\sum_{k=0}^{i-1}n_{k,\downarrow}+\sum_{k=0}^{N-1}n_{k,\uparrow}}\left[\cdot\cdot\cdot0_{i,\downarrow}\cdot\cdot\cdot|0\rangle\right]\end{gathered}

这里 0i,σ0_{i,\sigma} 表示该轨道未占据, 占据轨道则显式写出对应产生算符。向已占据轨道产生粒子、或从空轨道湮灭粒子, 结果都为 0。其余符号来自将算符移到约定位置时的反对易关系。

HH 也有可能产生额外的符号

[C3↓†C0↓][C0↓†C1↓†∣0⟩]=C3↓†C1↓†∣0⟩=−C1↓†C3↓†∣0⟩\left[C_{3\downarrow}^\dagger C_{0\downarrow}\right]\left[C_{0\downarrow}^\dagger C_{1\downarrow}^\dagger|0\rangle\right]=C_{3\downarrow}^\dagger C_{1\downarrow}^\dagger|0\rangle=-C_{1\downarrow}^\dagger C_{3\downarrow}^\dagger|0\rangle

Finite Temperature

于是可以计算期望值, 其中 β=1/(kBT)\beta=1/(k_BT), On=⟨ψn∣O∣ψn⟩O_n=\langle\psi_n|O|\psi_n\rangle, nn 遍历所选系综内的完整能量本征基。若分对称性扇区计算, 需要合并该系综允许的所有扇区; 只保留低能态时是低温截断近似。

⟨O⟩=Tr(e−βHO)Tr(e−βH)=∑nOne−βEn∑ne−βEn\langle O\rangle=\frac{Tr\left(e^{-\beta H}O\right)}{Tr\left(e^{-\beta H}\right)}=\frac{\sum_{n}O_{n}e^{-\beta E_{n}}}{\sum_{n}e^{-\beta E_{n}}}

电荷密度

⟨n⟩=⟨1N∑i=1N(ni↑+ni↓)⟩=1N∑i=1N(⟨ni↑⟩+⟨ni↓⟩)\langle n\rangle=\left\langle\frac{1}{N}\sum_{i=1}^{N}\left(n_{i\uparrow}+n_{i\downarrow}\right)\right\rangle=\frac{1}{N}\sum_{i=1}^{N}\left(\langle n_{i\uparrow}\rangle+\langle n_{i\downarrow}\rangle\right)

计算磁化, 设格点分成 a,b 两个子晶格

Ma(b)=∑i∈a(b)Siz=∑i∈a(b)(ni↑−ni↓)/2M_{a(b)}=\sum_{i\in a(b)}S_i^z=\sum_{i\in a(b)}\left(n_{i\uparrow}-n_{i\downarrow}\right)/2

均匀磁化与交替磁化

mu=1N(⟨Ma⟩+⟨Mb⟩)=12N∑i=1N(⟨ni↑⟩−⟨ni↓⟩)ms=1N(⟨Ma⟩−⟨Mb⟩)=12N[∑i∈a(⟨ni↑⟩−⟨ni↓⟩)−∑i∈b(⟨ni↑⟩−⟨ni↓⟩)]\begin{aligned}&m_{u}=\frac{1}{N}\left(\langle M_{a}\rangle+\langle M_{b}\rangle\right)=\frac{1}{2N}\sum_{i=1}^{N}\left(\langle n_{i\uparrow}\rangle-\langle n_{i\downarrow}\rangle\right)\\&m_{s}=\frac{1}{N}\left(\langle M_{a}\rangle-\langle M_{b}\rangle\right)=\frac{1}{2N}\left[\sum_{i\in a}\left(\langle n_{i\uparrow}\rangle-\langle n_{i\downarrow}\rangle\right)-\sum_{i\in b}\left(\langle n_{i\uparrow}\rangle-\langle n_{i\downarrow}\rangle\right)\right]\end{aligned}

磁化率

令 Mu=Ma+MbM_u=M_a+M_b, Ms=Ma−MbM_s=M_a-M_b, 扰动为 H(hu,hs)=H−huMu−hsMsH(h_u,h_s)=H-h_uM_u-h_sM_s, 则均匀磁化率与交错磁化率分别对对应的场求导:

χu=∂mu∂hu∣hu=hs=0,χs=∂ms∂hs∣hu=hs=0\chi_u=\left.\frac{\partial m_u}{\partial h_u}\right|_{h_u=h_s=0},\qquad \chi_s=\left.\frac{\partial m_s}{\partial h_s}\right|_{h_u=h_s=0}

这里 [H,Mu]=0[H,M_u]=0, 所以均匀磁化率可以写成等时涨落:

χu=βN[⟨Mu2⟩−⟨Mu⟩2]\chi_u=\frac{\beta}{N}\left[\langle M_u^2\rangle-\langle M_u\rangle^2\right]

但一般 [H,Ms]≠0[H,M_s]\neq0, 交错磁化率不能直接写成 βVar⁡(Ms)/N\beta\operatorname{Var}(M_s)/N, 而应使用虚时关联函数:

χs=1N∫0βdτ [⟨Ms(τ)Ms(0)⟩−⟨Ms⟩2],Ms(τ)=eτHMse−τH\chi_s=\frac{1}{N}\int_0^\beta d\tau\,\left[\langle M_s(\tau)M_s(0)\rangle-\langle M_s\rangle^2\right],\qquad M_s(\tau)=e^{\tau H}M_se^{-\tau H}

只有在 MsM_s 也与 HH 对易时, 才能化为同样的等时方差公式。热平均需要允许相应磁化涨落; 只取单一固定 SzS^z 扇区会使 MuM_u 的方差为零, 不能代替未限制磁化扇区的热平衡响应。

python package

see https://scipy.org

Real time evolution

量子系统的动力学演化

iℏ∂∂t∣ψ(t)⟩=H∣ψ(t)⟩i\hbar\frac{\partial}{\partial t}|\psi(t)\rangle=H|\psi(t)\rangle

对于不含时的 HH, 我们可以使用 ED 求解定态, 给出形式解

∣ψ(t)⟩=e−iHt/ℏ∣ψ(0)⟩=∑ne−iEnt/ℏ∣ψn⟩⟨ψn∣ψ(0)⟩=∑nCne−iEnt/ℏ∣ψn⟩\left|\psi(t)\right\rangle=e^{-iHt/\hbar}\left|\psi(0)\right\rangle=\sum_n e^{-iE_nt/\hbar}\left|\psi_n\right\rangle\left\langle\psi_n|\psi(0)\right\rangle=\sum_n C_ne^{-iE_nt/\hbar}\left|\psi_n\right\rangle

Hubbard model 的无相互作用极限 (U=0U=0), 这里固定一个自旋分量 σ\sigma, 取 hopping 振幅为 1, 考虑单粒子波包:

H=−∑j(aj,σ†aj+1,σ+H.c.),∣ψ(t=0)⟩=1Ω∑je−(α2/2)(j−N0)2eik0jaj,σ†∣0⟩H=-\sum_j\left(a_{j,\sigma}^\dagger a_{j+1,\sigma}+\mathrm{H.c.}\right),\quad|\psi(t=0)\rangle=\frac{1}{\sqrt{\Omega}}\sum_je^{-\left(\alpha^2/2\right)\left(j-N_0\right)^2}e^{ik_0j}a_{j,\sigma}^\dagger|0\rangle

Non-Hermitian Hamiltonian

若 Hamiltonian 不是 Hermitian

H†≠HH^\dagger\neq H

若 HH 可对角化, 可以选取完备的 bi - orthogonal 左右本征基 (简并子空间内也需选择互为对偶的基)。若遇到不可对角化的 exceptional point, 下面的展开不再直接适用, 需要广义本征矢和 Jordan 分解。

H∣Fn⟩=ωn∣Fn⟩⟨Gn∣H=⟨Gn∣ωnorH†∣Gn⟩=ωn∗∣Gn⟩\begin{aligned}&H\left|F_n\right\rangle=\omega_n\left|F_n\right\rangle\\&\left\langle G_n\right|H=\left\langle G_n\right|\omega_n\quad\mathrm{or}\quad H^\dagger\left|G_n\right\rangle=\omega_n^*\left|G_n\right\rangle\end{aligned}

单位算子

I=∑n∣Fn⟩⟨Gn∣⟨Gn∣Fn⟩\mathbb{I}=\sum_n\frac{\left|F_n\right\rangle\left\langle G_n\right|}{\left\langle G_n|F_n\right\rangle}

时间演化算子

e−iHt/ℏ=∑n∣Fn⟩e−iωnt/ℏ⟨Gn∣⟨Gn∣Fn⟩e^{-iHt/\hbar}=\sum_n\frac{\left|F_n\right\rangle e^{-i\omega_nt/\hbar}\left\langle G_n\right|}{\left\langle G_n|F_n\right\rangle}

在上述条件下可构造线性时间演化。非厄米演化通常不保持态的范数, 不能直接沿用幺正演化的概率解释。