摘要:本文用哈伯德模型研究各向异性ruby晶格中费米子行为, 在团簇动力学平均场理论框架内将格点模型映射为有效自洽场中的杂质模型后用连续-时间量子蒙特卡罗算法求解杂质模型. 基于自洽计算的结果, 用最大熵方法得到各向异性ruby晶格中具有相互作用的费米子体系的单粒子态密度和双占据数后讨论了温度(T )、相互作用(U )和各向异性参数(λ)对体系的金属-绝缘相变的影响. 最后给出各向异性ruby晶格中费米子体系的温度-相互作用相图, 研究结果表明, 低温和弱相互作用范围体系处在金属相, 而在高温和强相互作用下体系内出现Mott绝缘体相.
关键词: 各向异性ruby晶格/
哈伯德模型/
动力学平均场理论/
量子相变

English Abstract

全文HTML

--> --> -->

研究和发现二维强关联体系的拓扑绝缘体、量子霍尔效应、量子反常霍尔效应、玻色-爱因斯坦凝聚等新奇物相是凝聚态物理的重要内容 ^[1-7]. 拓扑绝缘体材料Bi₁₄Rh₃I₉^[8]的某一特定平面内存在的二维ruby晶格, 由于独特的晶格结构而引起研究人员兴趣. 拓扑绝缘体材料Bi₁₄Rh₃I₉的晶体结构中有周期性交错堆垛的铋-铑网格和绝缘层(图1(a)—(c)); 由共棱RhBi₈立方体覆盖六角晶格边所形成的金属间化合物的某一特定平面构成二维ruby晶格(图1(d)).
图 1 拓扑绝缘体材料Bi₁₄Rh₃I₉的晶格结构示意图　(a)?(c) Bi₁₄Rh₃I₉的晶体结构及其构成单元, 绝缘层[Bi₂I₈]^2–的zigzag链分离由共棱RhBi₈立方体构成的六角网格状金属间的 [(Rh₄Bi)₃I]²⁺层, 六角晶格的边由共棱RhBi₈立方体覆盖; (d) 二维ruby晶格与六角晶格结构俯视图^[8]
Figure1. Sketch of crystal structure of topological insulator Bi₁₄Rh₃I₉: (a)?(c) Triclinic crystal structure of Bi₁₄Rh₃I₉. Insulating layers of [Bi₂I₈]^2– zigzag chains separate the intermetallic [(RhBi₄)₃I]²⁺ layers that consist of hexagonal nets of edge-sharing RhBi₈ cubes; (d) honeycomb lattice of graphene with the structure of the intermetallic layer ^[8].

研究人员开展了诸多有关ruby晶格的研究工作并得到了丰富的成果. 文献[9, 10]介绍了用二维Ising模型描述ruby晶格和ruby晶格上二维冰剩余熵的研究. 近期研究人员在考虑自旋-轨道耦合效应的情况下用紧束缚模型或哈伯德模型研究二维ruby晶格后发现了拓扑绝缘体、分数量子霍尔效应和量子自旋液体等更为丰富的结果^[11-14]. 除此之外, 用Kitaev模型研究二维ruby晶格后发现了拓扑自旋液体相等的有趣结果^[15,16]. 然而, 目前为止还没有研究工作讨论各向异性ruby晶格中费米子间的在位相互作用和温度对该体系量子相变的影响.
本文以半满的单带哈伯德模型描述各向异性ruby晶格中具有在位相互作用的费米子行为, 结合团簇动力学平均场理论和连续-时间量子蒙特卡罗方法求解杂质模型.

2

-->

2.1.理论模型

二维ruby晶格是类六角晶格, 其结构相当于将六角晶格的顶点用三角晶格替换, 六角晶格的棱则用平方格子替换, 如图2(a)所示. 用哈伯德模型^[17-21]描述各向异性ruby晶格中具有在位相互作用的费米子行为的哈密顿量为
图 2 (a) Ruby晶格结构示意图; (b) ruby晶格第一布里渊区; (c) 各向异性ruby晶格($\lambda = 2.0$)中费米子体系$ T = 0.2$和$ U = 0 $时的态密度; (d) 各向异性ruby晶格($\lambda = 2.0$)中费米子体系无相互作用情况下色散关系
Figure2. (a) Sketch of ruby lattice; (b) first Brillouin zone of ruby lattice; (c) density of states of anisotropic ruby lattice ($\lambda = 2.0$) for $T = 0.2$, $U = 0$; (d) energy dispersion of anisotropic ruby lattice ($\lambda = 2.0$) for $T = 0.2$, $U = 0$.

$\begin{split} H = \;& - {t_1}\sum\limits_{ \langle i,j \rangle \in HL,\sigma } {(c_{i\sigma }^{\dagger } {c_{j\sigma }} + } c_{j\sigma }^{\dagger } {c_{i\sigma }}) \\& - {t_2}\sum\limits_{ \langle i,j \rangle \in HL \mapsto HL,\sigma } {(c_{i\sigma }^{\dagger } {c_{j\sigma }} + c_{j\sigma }^{\dagger } {c_{i\sigma }})} \\&+ U\sum\limits_I {n_{i \uparrow }^{\dagger } {n_{j \downarrow }}} - \mu \sum\limits_{i,\sigma } {{n_{i\sigma }}} , \end{split}$

其中, $c_{i\sigma }^{\dagger } $和${c_{j\sigma }}$分别表示i格点上自旋为σ的粒子生成算符和j 格点上自旋为σ的费米子湮灭算符, $ {n_{i\sigma }} $是粒子数算符, σ是自旋指标, 其值为↑或↓. t₁和t₂ 分别代表同一个六角晶格(HL)内最近邻格点之间的跃迁和紧邻两个六角晶格的最近邻格点之间的跃迁. $\langle i, j \rangle \in HL$表示同一个六角晶格内所有最近邻格点, $\langle i, j \rangle \in HL \mapsto HL$表示紧邻两个六角晶格之间所有最近邻的格点. U是在位相互作用, μ是化学势. 计算中定义各向异性参数$ \lambda = {{{t_1}} \mathord{\left/ {\vphantom {{{t_1}} {{t_2}}}} \right. } {{t_2}}} $并选定${t_1} = 1$来讨论ruby晶格的各向异性程度对体系量子相变的影响.
2 -->

2.2.研究方法

在解析求解描述强关联体系的Ising模型、Heisenberg模型、t-J模型和哈伯德模型等理论模型时遇到了极大的挑战, 即便处理简单的哈伯德模型和Kondo晶格模型的时候也遇到了困难. 随着计算机数据处理能力的提高和储存空间的扩大, 精确对角化^[22,23]、蒙特卡罗方法^[24,25]、Lanczos方法^[26]及重整化群^[27,28]等数值计算方法被应用在强关联体系的研究中并给出了令人满意的结果, 但以上数值计算方法各有其不足之处. 相比上述数值计算方法, 动力学平均场理论在强关联体系的研究中因给出更加令人满意的结果而得到了研究人员的青睐. 动力学平均场理论的建立与发展为强关联电子体系的研究开辟了新的途径, 主要成果有反铁磁金属到反铁磁绝缘体相变的动力学平均场理论、铜氧化物中反铁磁和d波高温超导的团簇动力学理论研究、六角晶格上半满哈伯德模型Mott相变、一维扩展哈伯德模型量子相变的团簇动力学平均场理论研究等^[29-35].
动力学平均场理论是将格点问题映射到有效自洽场中杂质模型的粗粒化方法^[36]. 巨正则系综的配分函数为 ${\cal Z} = \displaystyle\int {\prod\limits_{j, \mu , \sigma } {{\cal D}c_{i\sigma }^\dagger {\cal D}{c_{i\sigma }}{{\rm{e}}^{ - S}}} }$, 体系的作用量为

$\begin{split} S =\;& \displaystyle\int_0^\beta {{\rm{d}}\tau \left( {\sum\limits_{i\mu ,j\upsilon } {c_{i\mu }^{\dagger } O_{\mu \upsilon }^{ij}{\partial _\tau }{c_{i\mu }} - H[c_{i\mu }^{\dagger } ,{c_{j\upsilon }}]} } \right)}\\\equiv\;& {S_c} + {S_{cb}} + {S_b}, \end{split}$

式中, i, j是团簇的序号(此处将所选团簇序号标为$j = 0$), $\mu , \upsilon $是每个团簇内格点的指标, ${S_c}$包含团簇内的所有项的作用量, ${S_b}$包括所选团簇之外所有项的作用量, ${S_{bc}}$则包含连接团簇和其所处自洽场的所有项的作用量. 做如下两个假设: 1) 所有超过团簇大小的相互作用不考虑; 2) 有效作用量中不考虑四次或更高阶的重整化. 通过计算在$j \ne 0$的所有变量${c_{j\mu }}$($c_{j\mu }^{\dagger } $)上的路径积分, 可以获得只包含所选团簇变量${c_{0\mu }}$($c_{0\mu }^{\dagger } $)的(团簇自由度的)有效作用量:

$ \frac{1}{{{{\cal{Z}}_{{\rm{eff}}}}}}{{\rm{e}}^{ - {S_{{\rm{eff}}}}\left[ {c_{0\mu \sigma }^{\dagger } ,{c_{0\mu \sigma }}} \right]}} \equiv \frac{1}{\cal{Z}}\int {\mathop \prod \limits_{j \ne 0,\mu ,\sigma } } Dc_{j\mu \sigma }^{\dagger } D{c_{j\mu \sigma }}{{\rm{e}}^{ - S}}, $

$\begin{split} {S}_{\rm{eff}} =\;& -{\displaystyle {\int }_{0}^{\beta }{\rm{d}}\tau {\rm{d}}{\tau }^{\prime }}{\displaystyle \sum _{\mu \upsilon \sigma }{c}_{\mu \sigma }^{?}(\tau )}{\cal{G}}_{0,\mu \upsilon \sigma }^{-1}(\tau -{\tau }^{\prime }){c}_{\mu \sigma }({\tau }^{\prime })\\&+{\displaystyle {\int }_{0}^{\beta }{\displaystyle \sum _{\mu =1}^{N{\rm{c}}{c}}U{n}_{\mu \uparrow }(\tau ){n}_{\mu \downarrow }}}(\tau ), \\[-25pt]\end{split}$

式中, $\hat{\cal{G}}_{0, \mu \upsilon \sigma }^{\rm{-1}}(\tau -{\tau }')$为Weiss自洽场, $\mu , \upsilon = 1, 2, \cdots , $$ N{\rm{c}}$是团簇内格点的指标. 通过有效作用量可计算一个团簇的格林函数${\rm{G}}_{\mu \upsilon }^c(\tau - \tau ') \equiv - \left\langle {{T_\tau }c_\mu ^{\dagger } (\tau ){c_\upsilon }(\tau ')} \right\rangle $和自能$\hat{\displaystyle\sum }_{\mu \upsilon }^{c}({\rm{i}}{{\omega }_{n}})\equiv {\cal{G}}_{\mu \upsilon }^{-1}({\rm{i}}{{\omega }_{n}})-G{{_{\mu \upsilon }^{c}}^{-1}}({\rm{i}}{{\omega }_{n}})$, 其中忽略自旋的自由度. 利用Dyson方程可以从新计算Weiss动力学自洽场:

$\begin{split} {\cal{G}}_{0}^{-1}({\rm{i}}{{\omega }_{n}})=\;&{{\left( \sum\limits_{K}{\frac{1}{({\rm{i}}{{\omega }_{n}}+\mu )\hat{1}-\hat{t}(K)-{{{\hat{\sum }}}^{c}}({\rm{i}}{{\omega }_{n}})}} \right)}^{-1}}\\&+{{\hat{\sum }}^{c}}({\rm{i}}{{\omega }_{n}}). \\[-10pt]\end{split}$

以方格子为例, 进一步说明动力学平均场理论将格点模型映射到有效自洽场杂质模型的基本思想. 对图3所示的$2 \times 2$的超级晶格, 其团簇自能是$4 \times 4$的矩阵, K是超级晶格的简约布里渊区, I表示$4 \times 4$的单位矩阵, μ为化学势. t(K)是超级晶格的跃迁矩阵, 保留团簇内部的指标, 其矩阵元为
图 3 方晶格模型映射到有效自洽场杂质模型的示意图　(a) 方晶格; (b) 具有四个格点的超级晶格; (c)有效自洽场中团簇杂质模型的示意图
Figure3. Sketch of mapping square lattice model to impurity model in self-consistent field: (a) Square lattice; (b) sketch of supper lattice consists of 4 lattice point; (c) sketch of cluster impurity in self-consistent field.

$ {t_{ij}}(K) = \mathop \sum \limits_{{\boldsymbol r}_1 - \boldsymbol r_2} t_{ij}^{\boldsymbol r_1\boldsymbol r _2} \exp ({\rm i} K \cdot ({\boldsymbol r _1} - \boldsymbol r _2)), $

${{\boldsymbol{r}}_1}$和${{\boldsymbol{r}}_2}$是超晶格基矢, $i, j$为团簇内部格点的指标.
将格点模型映射到有效自洽场中杂质模型后需要用杂质求解器解杂质模型. 本工作选用的连续-时间量子蒙特卡罗算法^[37-40]对具有非局域影响和相互作用依赖时间的体系是有效的和可行的杂质求解器. 连续-时间量子蒙特卡罗方法在随机变量及其分布的计算中所用到的主要数学基础包括Metropolis重要抽样、大数定理和中心极限定理. 以本论文所涉及的“用半满的单带哈伯德模型描述各向异性ruby晶格中费米子体系的哈密顿量”为例介绍连续-时间量子蒙特卡罗算法.

$\begin{split} & H = {H_0} + {H_1} \\=\;& - t\sum\limits_{ \langle i,j \rangle \sigma } {(C_{i\sigma }^{\dagger } {C_{j\sigma }}} + {\rm{H}}{\rm{.c}}) \\&+ \sum\limits_{ij\sigma {\sigma'}} {{u_{ij}}{n_{i\sigma }}{n_{j{\sigma'}}}} + {\rm{ }}\mu \sum\limits_i {C_{i\sigma }^{\dagger } {C_{i\sigma }}} , \end{split}$

其中${H_0}$和${H_1}$分别为相互作用绘景中不含时间的部分和微扰的部分. 相互作用$U = 0$时, ruby晶格的哈伯德模型转变为紧束缚模型的哈密顿量, 其哈密顿量在动量空间中的表达式则变为${H_0} = \displaystyle\sum\nolimits_0 {\varPsi _k^{\dagger } } {\cal{H}}_k^0{\Psi _k}$, 其中

$ {\cal{H}}_k^0 = - \left[ {\begin{array}{*{20}{c}} 0&{{t_1}}&{{t_2}{{\rm{e}}^{ - {{\rm{i}}}{{(a{k_x} + b{k_y})} \mathord{\left/ {\vphantom {{(a{k_x} + b{k_y})} 2}} \right. } 2}}}}&0&{{t_2}{{\rm{e}}^{ - {{\rm{i}}}2b{k_y}}}}&{{t_1}} \\ {{t_1}}&0&{{t_1}}&{{t_2}{{\rm{e}}^{ - {{\rm{i}}}2b{k_y}}}}&0&{{t_2}{{\rm{e}}^{{{\rm{i}}}{{(a{k_x} - b{k_y})} \mathord{\left/ {\vphantom {{(a{k_x} - b{k_y})} 2}} \right. } 2}}}} \\ {{t_2}{{\rm{e}}^{{{\rm{i}}}{{(a{k_x} + b{k_y})} \mathord{\left/ {\vphantom {{(a{k_x} + b{k_y})} 2}} \right. } 2}}}}&{{t_1}}&0&{{t_1}}&{{t_2}{{\rm{e}}^{{\rm{i}}{{(a{k_x} - b{k_y})} \mathord{\left/ {\vphantom {{(a{k_x} - b{k_y})} 2}} \right. } 2}}}}&0 \\ 0&{{t_2}{{\rm{e}}^{{{\rm{i}}}b{k_y}}}}&{{t_1}}&0&{{t_1}}&{{t_2}{{\rm{e}}^{{{\rm{i}}}{{(a{k_x} + b{k_y})} \mathord{\left/ {\vphantom {{(a{k_x} + b{k_y})} 2}} \right. } 2}}}} \\ {{t_2}{{\rm{e}}^{{{\rm{i}}}b{k_y}}}}&0&{{t_2}{{\rm{e}}^{ - {{\rm{i}}}{{(a{k_x} - b{k_y})} \mathord{\left/ {\vphantom {{(a{k_x} - b{k_y})} 2}} \right. } 2}}}}&{{t_1}}&0&{_{{t_1}}} \\ {{t_1}}&{{t_2}{{\rm{e}}^{ - {{\rm{i}}}{{(a{k_x} - b{k_y})} \mathord{\left/ {\vphantom {{(a{k_x} - b{k_y})} 2}} \right. } 2}}}}&0&{{t_2}{{\rm{e}}^{ - {{\rm{i}}}{{(a{k_x} + b{k_y})} \mathord{\left/ {\vphantom {{(a{k_x} + b{k_y})} 2}} \right. } 2}}}}&{{t_1}}&0 \end{array}} \right]. $

矩阵中$a = {{(3 + \sqrt 3 )} \mathord{\left/ {\vphantom {{(3 + \sqrt 3 )} 2}} \right. } 2}$, $b = (1 + \sqrt 3)/2,~~ \varPsi _k^{\dagger } = (c_{Ak}^{\dagger } , c_{Bk}^{\dagger } , c_{Ck}^{\dagger } , c_{Dk}^{\dagger } , c_{Ek}^{\dagger } , c_{Fk}^{\dagger } ),$, A, B, C, D, E和F点分别是团簇动力学平均场理论框架下所选定的团簇内格点的指标, 如图2(a)所示. 相互作用绘景中体系的配分函数为: $ {\cal{Z}} = {\rm{Tr}}\left[ {\exp ( - \beta H)} \right] = {\rm{Tr}}{T_\tau }\left[ {\exp ( - \beta {H_0}) \cdot S(\beta )} \right] , $ 其中

$\begin{split} S(\beta ) =\;& {T_\tau }\exp \left[ { - \int_0^\beta {{H_1}} (\tau ){\rm{d}}\tau } \right] , ~~{H_1}(\tau ) = {{\rm{e}}^{\tau {H_0}}}{H_1}{{\rm{e}}^{ - \tau {H_0}}} = \sum\limits_{ij\sigma {\sigma'}} {{u_{ij}}{{\rm{e}}^{\tau {H_0}}}{n_{i\sigma }}{{\rm{e}}^{ - \tau {H_0}}}{n_{j{\sigma'}}}} = \sum\limits_{ij\sigma {\sigma'}} {{u_{ij}}{n_{i\sigma }}(\tau ){n_{j{\sigma'}}}} (\tau ) . \end{split}$

对作用量级数展开后体系的配分函数形式如下:

$ \begin{split} {\cal{Z}} =\;& {\rm{Tr}}(\exp ( - \beta {H_0}) \cdot S(\beta )) = \rm{Tr}\Bigg\{ \exp ( - \beta {H_0}) \cdot {T_\tau }\exp \Bigg[ - \int_0^\beta {{\rm{d}}\tau } {H_1}(\tau )\Bigg]\Bigg\} \\=\;& {\rm{Tr\Bigg\{ }}\exp ( - \beta {H_0})\sum\limits_k {\frac{{{{( - 1)}^k}}}{{k!}}} \int {{\rm{d}}{\tau _1}} \int {{\rm{d}}{\tau _2}} \cdots \int {{\rm{d}}{\tau _k}} {T_\tau }\Bigg[{H_1}({\tau _1}){H_1}({\tau _2}) \cdots {H_1}({\tau _k})\Bigg]\Bigg\} \\=\;& {\rm{Tr\Bigg\{ }}\exp ( - \beta {H_0})\sum\limits_k {\frac{{{{( - U)}^k}}}{{k!}}} \int {{\rm{d}}{\tau _1}} \cdots \int {{\rm{d}}{\tau _k}} {T_\tau }\Bigg[\sum\limits_{{i_1}} {{n_{{i_1} \uparrow }}({\tau _1})} {n_{{i_1} \downarrow }}({\tau _1}) \cdots \sum\limits_{{i_k}} {{n_{{i_k} \uparrow }}({\tau _{kl}})} {n_{{i_k} \downarrow }}({\tau _k})\Bigg] \Bigg\} ,\end{split}$

其中Tr是求矩阵迹的符号, ${T_\tau }$为编时算符. 在不引起误解的情况下, 将(8)式中时间指标、空间指标和级数展开指标可以简化为$\kappa \equiv \{ k, i, \{ {\tau _i}\} \} $, 则体系配分函数可写为如下形式:

$\begin{split} {\cal{Z}}=\;&{\cal{Z}}_{0}{\displaystyle \sum _{\kappa }\frac{{(-U)}^{\kappa }}{\kappa !}}{\displaystyle \int {\rm{d}}I}\mathrm{\cdots}{\displaystyle \int {\rm{d}}\kappa }\langle {T}_{\tau }{n}_{\uparrow }(1)\cdots {n}_{\uparrow }(\kappa )\rangle _{0}{\langle {T}_{\tau }{n}_{\downarrow }(1)\cdots {n}_{\downarrow }(\kappa )\rangle }_{0}, \end{split}$

暂且假定(9)式中被积函数均为正, 其中${{\cal{Z}}_0} = $$ {\rm{Tr}}\exp ( - \beta {H_0})$. 使用维克定理可使配分函数中的被积函数简化, 且任意可观测量O的平均值为

$ {\left\langle O \right\rangle _0} = \frac{{{\rm{Tr}}\exp ( - \beta {H_0}) \cdot O}}{{{\rm{Tr}}\exp ( - \beta {H_0})}}, $

因此被积函数中自旋向上部分的期望值为(以$\kappa = 3$的随机行走为例, 如图4所示):
图 4 $\kappa = 3$时一个随机行走示意图
Figure4. Sketch of random walk for$\kappa = 3$.

$ \begin{split} &{\left\langle {{T_\tau }{n_ \uparrow }(1){n_ \uparrow }(2){n_ \uparrow }(3)} \right\rangle _0} \\ =\;& {\left\langle {{T_\tau }c_ \uparrow ^{\dagger } (1){c_ \uparrow }(1)c_ \uparrow ^{\dagger } (2){c_ \uparrow }(2)c_ \uparrow ^{\dagger } (3){c_ \uparrow }(3)} \right\rangle _0} \\ =\;& \left| {\begin{array}{*{20}{c}} {{\cal{G}}_ \uparrow ^0(1,1)}&{{\cal{G}}_ \uparrow ^0(1,2)}&{{\cal{G}}_ \uparrow ^0(1,3)} \\ {{\cal{G}}_ \uparrow ^0(2,1)}&{{\cal{G}}_ \uparrow ^0(2,2)}&{{\cal{G}}_ \uparrow ^0(2,3)} \\ {{\cal{G}}_ \uparrow ^0(3,1)}&{{\cal{G}}_ \uparrow ^0(3,2)}&{{\cal{G}}_ \uparrow ^0(3,3)} \end{array}} \right|, \\[-6pt] \end{split} $

$ {{\cal{G}}^0}(i, j) $为相互作用大小等于0的情况下体系的格林函数. 自旋向下部分的期望值同理可得.
对于κ阶的情况,

$ \begin{split}&{\langle {T}_{\tau }{n}_{\uparrow }(1){n}_{\uparrow }(2)\cdots {n}_{\uparrow }(\kappa )\rangle }_{0}\\=\;&\left|\begin{array}{cccc}{\cal{G}}_{\uparrow }^{0}(1,1)& {\cal{G}}_{\uparrow }^{0}(1,2)& \cdots & {\cal{G}}_{\uparrow }^{0}(1,\kappa )\\ {\cal{G}}_{\uparrow }^{0}(2,1)& {\cal{G}}_{\uparrow }^{0}(2,2)& \cdots & {\cal{G}}_{\uparrow }^{0}(\kappa ,\kappa )\\ \vdots& \vdots& \vdots& \vdots\\ {\cal{G}}_{\uparrow }^{0}(\kappa ,1)& {\cal{G}}_{\uparrow }^{0}(\kappa ,2)& \cdots & {\cal{G}}_{\uparrow }^{0}(\kappa ,\kappa )\end{array}\right|\\ =\;&\mathrm{det}{D}_{\uparrow }(\kappa ), \end{split} $

自旋向下部分的期望值同理可得.
随机行走的过程中需要考虑减顶点和加顶点两种情况. Metropolis重要抽样算法的细致平衡条件为:

$ \frac{{{P_{k \to k + 1}}}}{{{P_{k + 1 \to k}}}} = - \frac{{U\beta N}}{{\kappa + 1}}\frac{{\displaystyle\prod\nolimits_\sigma {\det D(\kappa + 1)} }}{{\displaystyle\prod\nolimits_\sigma {\det D(\kappa )} }} . $

增加一个顶点的Metropolis重要抽样的接收率为

$ R = \min \left[ {1,\frac{{{P_{k \to k + 1}}}}{{{P_{k + 1 \to k}}}}} \right], $

具体计算中在0到1的区间取均匀分布的随机数, 直到随机数小于R方可结束随机行走. 减顶点和加顶点概率的计算中对行列式比值$ \dfrac{{\det D(\kappa + 1)}}{{\det D(\kappa )}} $的计算非常重要. 在已有的κ顶点增加一个顶点时, 快速计算$ \dfrac{{\det D(\kappa + 1)}}{{\det D(\kappa )}} $的方法如下:

$ \begin{split} \;& \dfrac{{\det D(\kappa + 1)}}{{\det D(\kappa )}} = \det (D(\kappa + 1)M(\kappa )) = \det (I + (D(\kappa + 1) - D(\kappa ))M(\kappa )) \\ =\;& \det \left[ {\begin{array}{*{20}{c}} 1&0& \cdots &{{{\cal{G}}^0}(1,\kappa + 1)} \\ 0&1& \cdots &{{{\cal{G}}^0}(2,\kappa + 1)} \\ 0&0& \cdots &{{{\cal{G}}^0}(3,\kappa + 1)} \\ \vdots & \vdots & \vdots & \vdots \\ {\displaystyle\sum\limits_i {{{\cal{G}}^0}(\kappa + 1,i)M{{(\kappa )}_{i,1}}} }&{\displaystyle\sum\limits_i {{{\cal{G}}^0}(\kappa + 1,i)M{{(\kappa )}_{i,2}}} }& \cdots &{{{\cal{G}}^0}(\kappa + 1,\kappa + 1)} \end{array}} \right] \end{split} , $

$ \dfrac{{\det D(\kappa + 1)}}{{\det D(\kappa )}} = \det \left[ {\begin{array}{*{20}{c}} {\boldsymbol{A}}&{\boldsymbol{D}} \\ {\boldsymbol{C}}&{\boldsymbol{B}} \end{array}} \right] $, 其中$ \left[ {\begin{array}{*{20}{c}} A&D \\ C&B \end{array}} \right] $由四块儿矩阵构成, 此块矩阵的元素分别是: A为左上角的单位矩阵, B为右下角的$ (\kappa {\kern 1 pt} {\kern 1 pt} {\kern 1 pt} + 1, \kappa {\kern 1 pt} {\kern 1 pt} {\kern 1 pt} + 1) $矩阵元, C为左下角的$({\rm{1}} \times \kappa )$的行矩阵, D是右上角的$(\kappa \times 1)$列矩阵. 因此,

$ \begin{split}&\dfrac{\mathrm{det}D(\kappa +1)}{\mathrm{det}D(\kappa )} = \mathrm{det}\left[\begin{array}{cc}A& D\\ C& B\end{array}\right] = \mathrm{det}A\cdot {\rm det}(B - C{A}^{-1}D ) = {\cal{G}}^{0}(\kappa +1,\kappa +1) - \displaystyle \sum _{ij}{\cal{G}}_{\sigma }^{0}(\kappa +1,i) {M}_{\sigma }{(\kappa )}_{i,j}{\cal{G}}_{\sigma }^{0}(j,\kappa +1), \end{split} $

从已有的$\kappa + 1$个顶点减去一个顶点的概率为

$\begin{split} \frac{{{P_{k \to k + 1}}}}{{{P_{k + 1 \to k}}}} = - \frac{{U\beta N}}{{\kappa + 1}}\prod\limits_\sigma \Big[ {{\cal{G}}^0}(\kappa + 1, \kappa + 1) - \sum\limits_{ij} {\cal{G}}_\sigma ^0(\kappa + 1,i) {M_\sigma }{{(\kappa )}_{i,j}}{\cal{G}}_\sigma ^0(j,\kappa + 1) \Big] . \end{split}$

至此未知的是$ {M_\sigma }(\kappa ) $矩阵, 而$ {M_\sigma }(\kappa ) $和${M_\sigma }\left(\right.\kappa + 1\left.\right)$的关系可由(15)式得到:

$\begin{split}&{M_\sigma }(\kappa + 1) = \det \left[ {\begin{array}{*{20}{c}} \cdot & \cdot & \cdot &{ - {\lambda ^{ - 1}}{L_{1,\kappa + 1}}} \\ \cdot &{{{M'}_\sigma }}& \cdot &{ - {\lambda ^{ - 1}}{L_{2,\kappa + 1}}} \\ \cdot & \cdot & \cdot &{ - {\lambda ^{ - 1}}{L_{3,\kappa + 1}}} \\ \vdots & \vdots & \vdots & \vdots \\ { - {\lambda ^{ - 1}}{R_{\kappa + 1,1}}}&{ - {\lambda ^{ - 1}}{R_{\kappa + 1,2}}}& \cdots &{ - {\lambda ^{ - 1}}} \end{array}} \right] , \end{split}$

矩阵元素为$ {M'_{\sigma ij}} = M{(\kappa )_{ij}} + {L_{i, \kappa + 1}}{\lambda ^{ - 1}}{R_{\kappa + 1, j}} , $ 其中

$\begin{split} {R_{ij}} =\;& \sum\limits_n {{\cal{G}}_\sigma ^0(i, n)M{{(\kappa )}_{nj}}} ,~~ {L_{ij}} = \sum\limits_n {M{{(\kappa )}_{in}}{\cal{G}}_\sigma ^0(n, j)}, ~~\lambda = \frac{{\det D(\kappa + 1)}}{{\det D(\kappa )}} .\end{split}$

完成抽样后, 对可观测的量进行计算. 以自旋向上的情况为例, 体系单粒子格林函数如下:

$\begin{split} G_{\uparrow }({\tau })_{{i}{j}}=\;&-\langle {T}_{\tau }{c}_{\uparrow }(i){c}_{\uparrow }^{?}(j)\rangle =-\frac{1}{\cal{Z}}{\rm{Tr}}\left({T}_{\tau }{c}_{\uparrow }(i){c}_{\uparrow }^{?}(j){\rm{e}}^{-\beta H}\right) \\=\;&-\frac{1}{\cal{Z}}{\rm{Tr}}\Bigg({T}_{\tau }{c}_{\uparrow }(i){c}_{\uparrow }^{?}(j)\mathrm{exp}(-\beta {H}_{0}){T}_{\tau }\mathrm{exp}\Bigg({\displaystyle {\int }_{\beta }^{0}U{\displaystyle \sum _{i}{n}_{i\uparrow }(\tau ){n}_{i\downarrow }(\tau )}}{\rm{d}}\tau \Bigg)\Bigg) \\=\;&-\frac{1}{\cal{Z}}{\rm{Tr}}[{T}_{\tau }{c}_{\uparrow }(i){c}_{\uparrow }^{?}(j)\mathrm{exp}(-\beta {H}_{0}){\displaystyle \sum _{\kappa }\frac{{(-U)}^{\kappa }}{\kappa !}}{\displaystyle \int {\rm{d}}l\cdots }{\displaystyle \int {\rm{d}}\kappa }{T}_{\tau }{n}_{i\uparrow }(1)\cdots {n}_{i\uparrow }(\kappa ){T}_{\tau }{n}_{i\downarrow }(1)\cdots {n}_{i\downarrow }(\kappa )] \\=\;&-\frac{{\cal{Z}}_{0}}{\cal{Z}}{\displaystyle \sum _{\kappa }\frac{{(-U)}^{\kappa }}{\kappa !}}{\displaystyle \int {\rm{d}}l\cdots }{\displaystyle \int {\rm{d}}\kappa }{\langle {T}_{\tau }{c}_{\uparrow }(i){c}_{\uparrow }^{?}(j){n}_{i\uparrow }(1)\cdots {n}_{i\uparrow }(\kappa )\rangle }_{0}{\langle {T}_{\tau }{n}_{i\downarrow }(1)\cdots {n}_{i\downarrow }(\kappa )\rangle }_{0}, \end{split}$

将维克定理应用于$ {\left\langle {{T_\tau }{c_ \uparrow }(i)c_ \uparrow ^{\dagger } (j){n_{i \uparrow }}(1) \cdots {n_{i \uparrow }}(\kappa )} \right\rangle _0} $ 中, 可得其值为$ \det D(\kappa + 1) $.

$ \begin{split} {G_ \uparrow }{(\tau )_{ij}} =\;& \dfrac{{\displaystyle\sum\limits_\kappa {\dfrac{{{{( - U)}^\kappa }}}{{\kappa !}}} \int {{\rm{d}}l \cdots } \int {{\rm{d}}\kappa } \det [{D_ \uparrow }(\kappa + 1){D_ \downarrow }(\kappa )]}}{{\displaystyle\sum\limits_\kappa {\frac{{{{( - U)}^\kappa }}}{{\kappa !}}} \int {{\rm{d}}l \cdots } \int {{\rm{d}}\kappa } \det [{D_ \uparrow }(\kappa ){D_ \downarrow }(\kappa )]}} \\ =\;& \dfrac{{\displaystyle\sum\limits_\kappa {\dfrac{{{{( - U)}^\kappa }}}{{\kappa !}}} \int {{\rm{d}}l \cdots } \int {{\rm{d}}\kappa } \det [{D_ \uparrow }(\kappa ){D_ \downarrow }(\kappa )]\frac{{\det [{D_ \uparrow }(\kappa + 1){D_ \downarrow }(\kappa )]}}{{\det [{D_ \uparrow }(\kappa ){D_ \downarrow }(\kappa )]}}}}{{\displaystyle\sum\limits_\kappa {\frac{{{{( - U)}^\kappa }}}{{\kappa !}}} \int {{\rm{d}}l \cdots } \int {{\rm{d}}\kappa } \det [{D_ \uparrow }(\kappa ){D_ \downarrow }(\kappa )]}}, \end{split} $

其中$\dfrac{{\det [{D_ \uparrow }(\kappa + 1){D_ \downarrow }(\kappa )]}}{{\det [{D_ \uparrow }(\kappa ){D_ \downarrow }(\kappa )]}}$为测量函数,

$\begin{split} &{D_ \uparrow }(\kappa + 1) \\=\;& \left[ {\begin{array}{*{20}{c}} {{\cal{G}}_ \uparrow ^0(1,1)}&{{\cal{G}}_ \uparrow ^0(1,2)}& \cdots &{{\cal{G}}_ \uparrow ^0(1,\kappa )}&{{\cal{G}}_ \uparrow ^0(1,j)} \\ {{\cal{G}}_ \uparrow ^0(2,1)}&{{\cal{G}}_ \uparrow ^0(2,2)}& \cdots &{{\cal{G}}_ \uparrow ^0(2,\kappa )}&{{\cal{G}}_ \uparrow ^0(2,j)} \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ {{\cal{G}}_ \uparrow ^0(\kappa ,1)}&{{\cal{G}}_ \uparrow ^0(\kappa ,2)}& \cdots &{{\cal{G}}_ \uparrow ^0(\kappa ,\kappa )}&{{\cal{G}}_ \uparrow ^0(\kappa ,j)} \\ {{\cal{G}}_ \uparrow ^0(i,1)}&{{\cal{G}}_ \uparrow ^0(i,2)}& \cdots &{{\cal{G}}_ \uparrow ^0(i,\kappa )}&{{\cal{G}}_ \uparrow ^0(i,j)} \end{array}} \right] . \end{split}$

因此

$\begin{split}\;&{G_ \uparrow }{(\tau )_{ij}} = \frac{{\det [{D_ \uparrow }(\kappa + 1){D_ \downarrow }(\kappa )]}}{{\det [{D_ \uparrow }(\kappa ){D_ \downarrow }(\kappa )]}} \\= \;&{\cal{G}}_ \uparrow ^0(i,j) - \sum\limits_{pq} {{\cal{G}}_ \uparrow ^0(i,p)} M{(\kappa )_{p,q}}{\cal{G}}_ \uparrow ^0(q,j) . \end{split}$

结合团簇动力学平均场理论和连续-时间量子蒙特卡罗算法进行数值计算的流程如下.
1) 利用微扰论给出一个小的初始团簇自由能.
2) 利用杂质求解器, 如: 连续-时间量子蒙特卡罗算法进行求解, 得到团簇的格林函数.
3) 对团簇格林函数做一次傅里叶变换.
4) 利用Dyson 方程

$\begin{split} {\cal{G}}_0^{ - 1}({{\rm{i}}}{\omega _n}) =\;& {\left( {\sum\limits_K {\frac{1}{{({{\rm{i}}}{\omega _n} + \mu )\hat{1} - \hat t(K) - {\displaystyle{\hat \sum }^c}({{\rm{i}}}{\omega _n})}}} } \right)^{ - 1}} \\&+ {\hat \sum ^c}({{\rm{i}}}{\omega _n}) \end{split}$

得到Wiess自洽场.
5) 从方程$\hat{\displaystyle\sum }_{\mu \upsilon }^{c}({\rm{i}}{{\omega }_{n}})\equiv {\cal{G}}_{\mu \upsilon }^{-1}({\rm{i}}{{\omega }_{n}})-G{{_{\mu \upsilon }^{c}}^{-1}}({\rm{i}}{{\omega }_{n}})$可以得到了新的自能.
循环运行上述计算步骤直到前后两个自能的差值达到所要精度求方可结束自洽计算, 自洽计算的具体流程可见图5.
图 5 自洽计算的流程图
Figure5. Flow chart of self-consistent calculation.

态密度(density of state)和双占据数(double occupancy)是二维强关联体系Mott转变研究中的两个重要参数. 本文用连续-时间量子蒙特卡罗方法求解自洽场中的杂质模型, 从虚时格林函数$G(\tau )$出发, 用最大熵方法^[39]计算得到各向异性ruby晶格中费米子体系的态密度

$\begin{split} D(\omega ) = - \dfrac{1}{{\rm{π }}} \displaystyle\sum\limits_{i = 1}^6 {\left[ {{\rm Im} {G_{ii}}(\omega - {{\rm{i}}}\delta )} \right]} , \end{split}$

其中i是用团簇动力学平均场理论将格点模型映射到自洽场中杂质模型后的团簇内格点的序号.
首先给出各向异性ruby晶格($\lambda = 2.0$)中费米子体系的固定温度($T = 0.2$)情况下对应不同相互作用的态密度(图6(a))和固定相互作用($U = $$ 8.0$)情况下对应不同温度的态密度(图6(b)). 之后比较了温度($T = 0.2$)和相互作用($U = 8.0$)都固定的情况下对应不同各向异性参数λ的态密度的演化(图7).
图 6 各向异性ruby晶格($\lambda =2.0 $)中费米子体系的态密度　(a) T = 0.2时不同相互作用对应的态密度; (b) U = 8时不同温度对应的态密度
Figure6. Density of states of anisotropic ruby lattice ($\lambda =2.0 $) with fermions: (a) Density of states for different interaction at $T = 0.2$; (b) density of states for different temperature at $U = 8.0 $.

图 7 固定温度($T = 0.2$)和固定相互作用($U = 8.0 $)情况下, 各向异性参数对ruby晶格中费米子体系态密度的影响
Figure7. Comparison of the effect of anisotropic parameter λ on the density of states of fermions in ruby lattice on for $T = 2 $ and $ U = 0 $.

如图7所示, 温度$T = 0.2$和相互作用$U = 8.0 $的情况下, 随着各向异性参数λ的增大费米面两侧态密度谱峰逐渐演化到最后在$\lambda \approx 1$时费米面处出现能隙. $\lambda = 2$和$\lambda = 0.67$时态密度的演化形式类似, 不出现近藤峰. $\lambda = 1$态密度演化过程中出现准粒子峰, 即近藤峰, 其特征是松原频率$\omega = 0$处两侧出现具有类似肩膀的准粒子峰. 由态密度演化形式可以推测, 在$\lambda = 2$到$\lambda = 1$的过程中, 体系中会出现近藤金属. 各向异性参数对态密度的演化发现和赝能隙的形成具有显著的影响. 在固定的排斥相互作用能情况下, 通过比较不同各向异性参数所对应的态密度发现, 松原频率$\omega = 0$处态密度随着各向异性参数的增大而减小. 严格来说, 在有限温度下由$\omega = 0$处态密度的赝能隙决定的相变, 实际上是一个转变(crossover).
双占据数是用半满的哈伯德模型描述强关联体系金属-绝缘相变的另一个重要参数. 双占据数定义为$\mathrm{D}\mathrm{o}\mathrm{c}\mathrm{c}=\partial F/\partial U=\dfrac{1}{6}{\displaystyle\sum }_{i}\left\langle{{n}_{i\uparrow }{n}_{i\downarrow }}\right\rangle$, 其中F是自由能, U是相互作用. 图8(a), (b)分别给出不同相互作用情况下各向异性ruby晶格($\lambda = 2$)的双占据数随温度的变化和不同温度情况下各向异性ruby晶格($\lambda = 2$)的双占据数随相互作用的变化. 如图8(a)所示, 随着温度的降低, 温度对于双占据数的影响趋于不明显. 由图8(b)可知, 随着相互作用的增加体系的双占据数趋于0, 意味着费米子的局域化程度增强.
图 8 各向异性ruby晶格中费米子体系的双占据数　(a) 不同相互作用下体系双占据数随温度的变化; (b) 不同温度下体系双占据数随相互作用的变化
Figure8. Double occupancy of anisotropic ruby lattice with fermions: (a) Comparison between double occupancy for different temperature with the change of interaction; (b) comparison between double occupancy for different interaction with the change of temperature.

图9所示为各向异性参数对体系双占据数的影响. 在固定温度($T = 0.2$)和固定相互作用($T = $$ 8.0$)情况下, ruby晶格中费米子体系的双占据数随着各向异性参数的增大趋于0, 这一趋势说明费米子局域化加强, 即体系趋于Mott绝缘体相, 但趋势放缓.
图 9 温度和相互作用固定的情况下, ruby晶格中费米子体系的双占据数随各向异性参数的变化
Figure9. Trend of double occupancy with the change of anisotropic parameter for fixed T and U.

最后, 基于体系的态密度和双占据数, 给出了体现温度和相互作用对各向异性ruby晶格($\lambda = $$ 2.0 $)中相互作用费米子体系金属-绝缘相变的影响, 即温度-相互作用相图. 如图10所示, 相图被二阶Mott相变线划分为两个区域, 即低温和弱相互作用区的金属相、高温和强相互作用区域的Mott 绝缘体相. 随着温度的降低, 体系金属-绝缘相变所对应的相互作用也减小.
图 10 各向异性ruby晶格中费米子体系的金属-绝缘相图
Figure10. Metal-Insulator phase diagram of fermions in anisotropic ruby lattice.

本文中用半满的哈伯德模型描述各向异性ruby晶格中费米子. 结合团簇动力学平均场理论和连续时间量子蒙特卡罗方法求解杂质模型. 在自洽计算的基础上, 通过进一步系统的运算后给出体系的单粒子态密度和双占据数. 结合态密度和双占据数的演变趋势, 最终得到了各向异性ruby晶格中费米子体系的温度-相互作用相图. 研究结果发现了体系低温弱相互作用区的金属相和高温强相互作用区的Mott绝缘相.