学术, 其他笔记

高陈数的蜂窝双层系统中的拓扑隐藏相变 Topological hidden phase transition in honeycomb bilayers with a high Chern number

发布时间：

最近修改时间：2026年4月18日

置顶说明

文章中涉及到实空间哈密顿量的结论（电导部分）可能存在部分问题，可以忽略；而文章的主题“拓扑隐藏相变现象”是对的，因为这是由倒空间的哈密顿量计算出的结论。另外，半反手性边缘态好像还是有的，这个物理图像是对的。

这个计算上的问题是由评论区指出，更正了蜂窝双层系统实空间的 modified_Haldane_model 部分，尤其是其 h01 的构建，涉及到向左向右跃迁方向的问题，不能只是简单修改相位，套用 Haldane_model 模板，详细请阅读下面的评论。另外，在论文中的附录 A 中，关于双层模型的倒空间哈密顿量，也特别提到了：“Note that to get the unit cell of the bilayer honeycomb lattice containing all possible layer hoppings [Fig.7(a)], here one layer (HM layer or mHM layer) needs to be applied with its counterpart of mirror symmetry along x direction.”。

评论区提到，按照更正后的计算结果，论文中 Fig3 的电导变化趋势是基本正确的，只是电导变为 2，即打开全局能隙的临界 tc 会减小，大致会在 3\sqrt{3}t_2 。

以下是由评论区给出的代码（推荐用这个版本作为参考和后续研究）：

def h00_of_modified_Haldane_model(N, M, t1, t2, phi, sign):
    h00 = np.zeros((4*N, 4*N), dtype=complex)
    for i0 in range(N):
        h00[i0*4+0, i0*4+0] = M
        h00[i0*4+1, i0*4+1] = -M
        h00[i0*4+2, i0*4+2] = M
        h00[i0*4+3, i0*4+3] = -M
        h00[i0*4+0, i0*4+1] = t1
        h00[i0*4+1, i0*4+0] = t1
        h00[i0*4+1, i0*4+2] = t1
        h00[i0*4+2, i0*4+1] = t1
        h00[i0*4+2, i0*4+3] = t1
        h00[i0*4+3, i0*4+2] = t1
        # h00[i0*4+0, i0*4+2] = t2*cmath.exp(-1j*phi*sign)
        h00[i0*4+0, i0*4+2] = t2*cmath.exp(1j*phi*sign)
        h00[i0*4+2, i0*4+0] = h00[i0*4+0, i0*4+2].conj()
        # h00[i0*4+1, i0*4+3] = t2*cmath.exp(1j*phi*sign)
        h00[i0*4+1, i0*4+3] = t2*cmath.exp(-1j*phi*sign)
        h00[i0*4+3, i0*4+1] = h00[i0*4+1, i0*4+3].conj()
    for i0 in range(N-1):
        h00[i0*4+3, (i0+1)*4+0] = t1
        h00[(i0+1)*4+0, i0*4+3] = t1
        # h00[i0*4+2, (i0+1)*4+0] = t2*cmath.exp(1j*phi*sign)
        h00[i0*4+2, (i0+1)*4+0] = t2*cmath.exp(-1j*phi*sign)
        h00[(i0+1)*4+0, i0*4+2] = h00[i0*4+2, (i0+1)*4+0].conj()
        # h00[i0*4+3, (i0+1)*4+1] = t2*cmath.exp(-1j*phi*sign)
        h00[i0*4+3, (i0+1)*4+1] = t2*cmath.exp(1j*phi*sign)
        h00[(i0+1)*4+1, i0*4+3] = h00[i0*4+3, (i0+1)*4+1].conj()
    return h00

def h01_of_modified_Haldane_model(N, M, t1, t2, phi, sign):
    h01 = np.zeros((4*N, 4*N), dtype=complex)
    for i0 in range(N):
        # h01[i0*4+1, i0*4+0] = t1
        # h01[i0*4+2, i0*4+3] = t1
        h01[i0*4+0, i0*4+1] = t1
        h01[i0*4+3, i0*4+2] = t1
        h01[i0*4+0, i0*4+0] = t2*cmath.exp(1j*phi*sign)
        h01[i0*4+1, i0*4+1] = t2*cmath.exp(1j*phi*sign)
        h01[i0*4+2, i0*4+2] = t2*cmath.exp(1j*phi*sign)
        h01[i0*4+3, i0*4+3] = t2*cmath.exp(1j*phi*sign)
        # h01[i0*4+1, i0*4+3] = t2*cmath.exp(-1j*phi*sign)
        # h01[i0*4+2, i0*4+0] = t2*cmath.exp(-1j*phi*sign)
        h01[i0*4+3, i0*4+1] = t2*cmath.exp(-1j*phi*sign)
        h01[i0*4+0, i0*4+2] = t2*cmath.exp(-1j*phi*sign)
        if i0 != 0:
            # h01[i0*4+1, (i0-1)*4+3] = t2*cmath.exp(-1j*phi*sign)
            h01[i0*4+0, (i0-1)*4+2] = t2*cmath.exp(-1j*phi*sign)
    for i0 in range(N-1):
        # h01[i0*4+2, (i0+1)*4+0] = t2*cmath.exp(-1j*phi*sign)
        h01[i0*4+3, (i0+1)*4+1] = t2*cmath.exp(-1j*phi*sign)
    return h01

原文

评价：一般来说，价带内部的拓扑不变量的变化很难体现在量子输运的测量中，但通过研究发现，在这个复合蜂窝双层结构中，半反手性边缘态 (half-antichiral edge states) 2e2/h 的电导平台的出现可以作为拓扑隐藏相变 (Topological hidden phase transition) 的一个重要指标或测量方法。

摘要

根据体边对应关系，在量子反常霍尔系统中一维条带通常在两条边上表现出两个无能隙的边缘态，可以由陈数来描述。这里我们引入了一个模型，该模型展示了一种被称为“半反手性”边缘模式的现象。该模型为霍尔丹模型层和修改后的霍尔丹模型层组成的复合双层系统。“半反手性”边缘态在同一边缘上沿相同方向传播，并且通过体内的反向传播模式得到补偿。在强层间耦合的情况下，借助霍尔丹模型层中的手性边缘态，“半反手性”边缘态在整个系统中在无序下表现出比较强的稳定性，从而形成双电导平台。有趣的是，在这一过程中，我们观察到一种拓扑隐藏相变，其中，价带由陈数 C=±(1+1) 演变成一个具有高陈 数 C=±2 的价带和一个陈数为 C=0 的平凡价带，这和双电导平台密切相关。这些研究发现为发展层辅助的拓扑态以及在具有高陈数的系统中研究新的拓扑相变提供了更多的可能性。

目录 隐藏
置顶说明
原文
摘要
1. 主要研究背景
2. 蜂窝双层系统中的半反手性边缘态
2.1 哈密顿量
2.2 能带结构
2.3 波函数分布
2.4 量子输运性质
3. 蜂窝双层系统中的拓扑隐藏相变
4. 鲁棒的半反手性边缘态和拓扑隐藏相变的关系
5. 更多内容
6. 总结
7. 代码开源
代码一
代码二

1. 主要研究背景

Question & Motivation 1：

根据体边对应关系，拓扑手性边缘态的稳定性由价带的总陈数来决定。而价带内部的陈数变化，即拓扑隐藏相变 (topological hidden phase transition)，是否也可以反应在量子输运中？

Question & Motivation 2：

由于反手性边缘态 (antichiral edge states) 由手征对称性 (chiral symmetry) 保护，其抗无序的能力会比手性边缘态弱很多。那么，是否有办法提高反手性边缘态的量子输运性质？

手性边缘态和反手性边缘态的示意图：

手性边缘态和反手性边缘态的量子输运性质对比：

结论：虽然反手性边缘态已经在光学晶体、声子晶格、电路等系统中得到了观测，也具有一定的鲁棒性，但和手性边缘态对比来说，其抗无序的能力会弱很多。

2. 蜂窝双层系统中的半反手性边缘态

2.1 哈密顿量

蜂窝双层系统示意图：

实空间哈密顿量：

2.2 能带结构

Haldane model (HM) ribbon 的能带结构（可以看到手性边缘态）：

modified-Haldane model (mHM) ribbon 的能带结构（可以看到反手性边缘态）：

蜂窝双层系统的能带结构（在不同层间耦合强度 tc 下）：

2.3 波函数分布

可以看出：

  • E1, E2 的波函数分布只在条带的一个边缘上，且均匀分布在 HM, mHM 层上。
  • B1, B2 的波函数分布在条带的内部，但偏向条带的另一个边缘，它主要分布在 mHM 层上。

2.4 量子输运性质

无层间耦合（左图）和强层间耦合（右图）下的量子输运性质：

可以看出：

  • 在无层间耦合下，电导随无序的变化情况为两个单层的电导图的叠加，有一个 e2/h 的电导平台。
  • 有意思的是：在强层间耦合下，电导可以稳定到 2e2/h 的平台上。

不同层间耦合强度 tc 下的量子输运性质：

主要结论：随着层间耦合强度 tc 的增加，电导 G 先减小后增加，会从 e2/h 的平台（左图蓝线）到 2e2/h 的平台（左图红线）。这从右图中也能看出来。

这说明了存在两个散射过程的竞争：（1）HM 层的手性边缘态被体态散射后减弱量子输运；（2）mHM 层的反手性边缘态和 HM 层的手性边缘态互相散射增强输运。

不同条带宽度下的量子输运性质：

主要现象：

  • 现象1（左图）：随着条带宽度的增加，电导维持在 2e2/h 的平台能力越强。这说明了量子输运过程中受到体态散射的影响。
  • 现象2（右图）：可以看到当条带宽度 Ny0 足够大时，有一个临界值 tc，这个突变现象无法简单从散射过程的角度来解释。

3. 蜂窝双层系统中的拓扑隐藏相变

Haldane model 二维系统的能带结构（ky=0）：

modified-Haldane model 二维系统的能带结构（ky=0）：

二维蜂窝双层系统的能带结构（ky=0，在不同层间耦合强度 tc 下）：

对应的准一维条带系统的能带结构（条带宽度 Ny0 需要选取得足够大，这里选为 200）：

二维系统和准一维条带系统的能带结构对比：

可以看出：发生拓扑隐藏相变时，对应的条带能带的主要体带隙刚好为打开状态。这是因为在比较强的层间耦合下，如果要发生拓扑隐藏相变，那么原先的 mHM 半金属的能带结构（狄拉克锥）需要往导带或价带内部收，从而形成这个条带的间接带隙。

陈数 C、电导 G、带隙 ∆ 随层间耦合强度的变化情况：

结论：

  • 价带内部陈数（Cv1和Cv2）变化的 tc 临界值 和 价带内部的带隙 ∆v 变化的 tc 临界值是一致的，在 0.31 附近。
  • 电导 G 变化的 tc 临界值 和 条带系统能带结构的带隙 ∆cv_ribbon 变化的 tc 临界值是一致的，也在 0.31 附近。这个数值是偶然的吗？下面给出了具体的分析。

4. 鲁棒的半反手性边缘态和拓扑隐藏相变的关系

陈数 C、电导 G、带隙 ∆ 在不同次近邻跃迁强度 t2 下随层间耦合 tc 的变化情况：

可以发现，当次近邻跃迁强度 t2 越大，所需要的临界层间耦合强度 tc 也越大。

tc 临界值的总结：

结论：

  • 价带内部陈数（Cv1和Cv2）变化的 tc 临界值 和 价带内部的带隙 ∆v 变化的 tc 临界值是一致的。
  • 电导 G 变化的 tc 临界值 和 条带系统能带结构的带隙 ∆cv_ribbon 变化的 tc 临界值是一致的。
  • 随着 t2 的增加，电导达到 2e2/h 的平台所需要的 tc 值会略大于发生拓扑隐藏相变所需要的 tc 值。
  • 然而对于实际体系，次近邻跃迁强度 t2 通常是比较小的，因此它们保持了一致性。电导为 2e2/h 的平台可以作为拓扑隐藏相变的一个重要的判断指标。

一个较大的次近邻跃迁强度 t2 的例子（层间耦合强度 tc 选为拓扑隐藏相变点）：

可以发现，当次近邻跃迁强度 t2 比较大时，发生拓扑隐藏相变点时，条带的能带结构中体带隙还未完全打开，但这个偏差并不影响本篇文章的主要结论，这是因为：（1）当探测到 2e2/h 电导平台，系统必然已经发生拓扑隐藏相变。（2）次近邻跃迁强度 t2 的值一般比较小，因此拓扑隐藏相变和 2e2/h 电导平台的探测基本上保持一致性。

5. 更多内容

不同堆垛方式的影响：

通过量子点接触 (quantum point contact) 电调控半反手性边缘态：

6. 总结

  1. 随着层间耦合强度的增加，半反手性边缘态（half-antichiral edge states）变得鲁棒，电导达到了 2e2/h 的平台，准一维的条带系统发生了从半金属态向拓扑绝缘态的转变。
  2. 在这个过程中，二维系统价带的内部同时也发生了“拓扑隐藏相变”（topological hidden phase transition），价带陈数由 C=±(1+1) 转变成 C=±(2+0)。
  3. 在这个复合蜂窝双层结构中，半反手性边缘态 2e2/h 的电导平台的出现可以作为拓扑隐藏相变的一个重要指标或测量方法。

7. 代码开源

这里开源主要的实现代码。更多的代码未整理，暂不开源，如果需要更多代码，可联系我。

代码一

import numpy as np
import functools
from math import *
import cmath
import guan

# Installation: pip install --upgrade guan

def main():
    N = 15
    M = 0
    t1 = 1
    phi= pi/2
    t2= 0.03
    tc = 0.8
    layer_num = 2
    bands_of_bilayer_model(N, M, t1, t2, phi, tc, layer_num)
    conductance_of_bilayer_model_with_disorder_array(N, M, t1, t2, phi, tc, layer_num)

# Haldane model

def hamiltonian_of_Haldane_model(k, N, M, t1, t2, phi, sign):
    h00 = h00_of_Haldane_model(N, M, t1, t2, phi, sign)
    h01 = h01_of_Haldane_model(N, M, t1, t2, phi, sign)
    hamiltonian = h00 + h01*cmath.exp(1j*k) + h01.transpose().conj()*cmath.exp(-1j*k)
    return hamiltonian

def h00_of_Haldane_model(N, M, t1, t2, phi, sign):
    h00 = np.zeros((4*N, 4*N), dtype=complex)
    for i0 in range(N):
        h00[i0*4+0, i0*4+0] = M
        h00[i0*4+1, i0*4+1] = -M
        h00[i0*4+2, i0*4+2] = M
        h00[i0*4+3, i0*4+3] = -M
        h00[i0*4+0, i0*4+1] = t1
        h00[i0*4+1, i0*4+0] = t1
        h00[i0*4+1, i0*4+2] = t1
        h00[i0*4+2, i0*4+1] = t1
        h00[i0*4+2, i0*4+3] = t1
        h00[i0*4+3, i0*4+2] = t1
        h00[i0*4+0, i0*4+2] = t2*cmath.exp(1j*phi*sign)
        h00[i0*4+2, i0*4+0] = h00[i0*4+0, i0*4+2].conj()
        h00[i0*4+1, i0*4+3] = t2*cmath.exp(1j*phi*sign)
        h00[i0*4+3, i0*4+1] = h00[i0*4+1, i0*4+3].conj()
    for i0 in range(N-1):
        h00[i0*4+3, (i0+1)*4+0] = t1
        h00[(i0+1)*4+0, i0*4+3] = t1
        h00[i0*4+2, (i0+1)*4+0] = t2*cmath.exp(-1j*phi*sign)
        h00[(i0+1)*4+0, i0*4+2] = h00[i0*4+2, (i0+1)*4+0].conj()
        h00[i0*4+3, (i0+1)*4+1] = t2*cmath.exp(-1j*phi*sign)
        h00[(i0+1)*4+1, i0*4+3] = h00[i0*4+3, (i0+1)*4+1].conj()
    return h00

def h01_of_Haldane_model(N, M, t1, t2, phi, sign):
    h01 = np.zeros((4*N, 4*N), dtype=complex)
    for i0 in range(N):
        h01[i0*4+1, i0*4+0] = t1
        h01[i0*4+2, i0*4+3] = t1
        h01[i0*4+0, i0*4+0] = t2*cmath.exp(-1j*phi*sign)
        h01[i0*4+1, i0*4+1] = t2*cmath.exp(1j*phi*sign)
        h01[i0*4+2, i0*4+2] = t2*cmath.exp(-1j*phi*sign)
        h01[i0*4+3, i0*4+3] = t2*cmath.exp(1j*phi*sign)
        h01[i0*4+1, i0*4+3] = t2*cmath.exp(-1j*phi*sign)
        h01[i0*4+2, i0*4+0] = t2*cmath.exp(1j*phi*sign)
        if i0 != 0:
            h01[i0*4+1, (i0-1)*4+3] = t2*cmath.exp(-1j*phi*sign)
    for i0 in range(N-1):
        h01[i0*4+2, (i0+1)*4+0] = t2*cmath.exp(1j*phi*sign)
    return h01

# modified_Haldane_model

def hamiltonian_of_modified_Haldane_model(k, N, M, t1, t2, phi, sign):
    h00 = h00_of_modified_Haldane_model(N, M, t1, t2, phi, sign)
    h01 = h01_of_modified_Haldane_model(N, M, t1, t2, phi, sign)
    hamiltonian = h00 + h01*cmath.exp(1j*k) + h01.transpose().conj()*cmath.exp(-1j*k)
    return hamiltonian

def h00_of_modified_Haldane_model(N, M, t1, t2, phi, sign):
    h00 = np.zeros((4*N, 4*N), dtype=complex)
    for i0 in range(N):
        h00[i0*4+0, i0*4+0] = M
        h00[i0*4+1, i0*4+1] = -M
        h00[i0*4+2, i0*4+2] = M
        h00[i0*4+3, i0*4+3] = -M
        h00[i0*4+0, i0*4+1] = t1
        h00[i0*4+1, i0*4+0] = t1
        h00[i0*4+1, i0*4+2] = t1
        h00[i0*4+2, i0*4+1] = t1
        h00[i0*4+2, i0*4+3] = t1
        h00[i0*4+3, i0*4+2] = t1
        h00[i0*4+0, i0*4+2] = t2*cmath.exp(-1j*phi*sign)
        h00[i0*4+2, i0*4+0] = h00[i0*4+0, i0*4+2].conj()
        h00[i0*4+1, i0*4+3] = t2*cmath.exp(1j*phi*sign)
        h00[i0*4+3, i0*4+1] = h00[i0*4+1, i0*4+3].conj()
    for i0 in range(N-1):
        h00[i0*4+3, (i0+1)*4+0] = t1
        h00[(i0+1)*4+0, i0*4+3] = t1
        h00[i0*4+2, (i0+1)*4+0] = t2*cmath.exp(1j*phi*sign)
        h00[(i0+1)*4+0, i0*4+2] = h00[i0*4+2, (i0+1)*4+0].conj()
        h00[i0*4+3, (i0+1)*4+1] = t2*cmath.exp(-1j*phi*sign)
        h00[(i0+1)*4+1, i0*4+3] = h00[i0*4+3, (i0+1)*4+1].conj()
    return h00

def h01_of_modified_Haldane_model(N, M, t1, t2, phi, sign):
    h01 = np.zeros((4*N, 4*N), dtype=complex)
    for i0 in range(N):
        h01[i0*4+1, i0*4+0] = t1
        h01[i0*4+2, i0*4+3] = t1
        h01[i0*4+0, i0*4+0] = t2*cmath.exp(1j*phi*sign)
        h01[i0*4+1, i0*4+1] = t2*cmath.exp(1j*phi*sign)
        h01[i0*4+2, i0*4+2] = t2*cmath.exp(1j*phi*sign)
        h01[i0*4+3, i0*4+3] = t2*cmath.exp(1j*phi*sign)
        h01[i0*4+1, i0*4+3] = t2*cmath.exp(-1j*phi*sign)
        h01[i0*4+2, i0*4+0] = t2*cmath.exp(-1j*phi*sign)
        if i0 != 0:
            h01[i0*4+1, (i0-1)*4+3] = t2*cmath.exp(-1j*phi*sign)
    for i0 in range(N-1):
        h01[i0*4+2, (i0+1)*4+0] = t2*cmath.exp(-1j*phi*sign)
    return h01

# bilayer model

def hamiltonian_of_bilayer_model(k, N, M, t1, t2, phi, tc, layer_num):
    modified_Haldane_model = hamiltonian_of_modified_Haldane_model(k=k, N=N, M=M, t1=t1, t2=t2, phi=phi, sign=1)
    Haldane_model = hamiltonian_of_Haldane_model(k=k, N=N, M=M, t1=t1, t2=t2, phi=phi, sign=1)
    hamiltonian = np.zeros((4*N*layer_num, 4*N*layer_num), dtype=complex)
    for layer in range(layer_num):
        if np.mod(layer,2) == 0:
            hamiltonian[layer*4*N+0:layer*4*N+4*N, layer*4*N+0:layer*4*N+4*N] = modified_Haldane_model
        if np.mod(layer,2) == 1:
            hamiltonian[layer*4*N+0:layer*4*N+4*N, layer*4*N+0:layer*4*N+4*N] = Haldane_model 
    for layer in range(layer_num-1):
        for i0 in range(N):
            # AB stacking
            hamiltonian[layer*4*N+i0*4+0, (layer+1)*4*N+i0*4+1] = tc
            hamiltonian[(layer+1)*4*N+i0*4+1, layer*4*N+i0*4+0] = tc
            hamiltonian[layer*4*N+i0*4+2, (layer+1)*4*N+i0*4+3] = tc
            hamiltonian[(layer+1)*4*N+i0*4+3, layer*4*N+i0*4+2] = tc
            
            # # # BA stacking
            # hamiltonian[layer*4*N+i0*4+1, (layer+1)*4*N+i0*4+0] = tc
            # hamiltonian[(layer+1)*4*N+i0*4+0, layer*4*N+i0*4+1] = tc
            # hamiltonian[layer*4*N+i0*4+3, (layer+1)*4*N+i0*4+2] = tc
            # hamiltonian[(layer+1)*4*N+i0*4+2, layer*4*N+i0*4+3] = tc

            # # AA stacking
            # hamiltonian[layer*4*N+i0*4+0, (layer+1)*4*N+i0*4+0] = tc
            # hamiltonian[(layer+1)*4*N+i0*4+0, layer*4*N+i0*4+0] = tc
            # hamiltonian[layer*4*N+i0*4+1, (layer+1)*4*N+i0*4+1] = tc
            # hamiltonian[(layer+1)*4*N+i0*4+1, layer*4*N+i0*4+1] = tc
            # hamiltonian[layer*4*N+i0*4+2, (layer+1)*4*N+i0*4+2] = tc
            # hamiltonian[(layer+1)*4*N+i0*4+2, layer*4*N+i0*4+2] = tc
            # hamiltonian[layer*4*N+i0*4+3, (layer+1)*4*N+i0*4+3] = tc
            # hamiltonian[(layer+1)*4*N+i0*4+3, layer*4*N+i0*4+3] = tc
    return hamiltonian

def h00_of_bilayer_model(N, M, t1, t2, phi, tc, layer_num):
    h00_modified = h00_of_modified_Haldane_model(N, M, t1, t2, phi, sign=1)
    h00_Haldane = h00_of_Haldane_model(N, M, t1, t2, phi, sign=1)
    h00  = np.zeros((4*N*layer_num, 4*N*layer_num), dtype=complex) 
    for layer in range(layer_num):
        if np.mod(layer,2) == 0:
            h00[layer*4*N+0:layer*4*N+4*N, layer*4*N+0:layer*4*N+4*N] = h00_modified
        if np.mod(layer,2) == 1:
            h00[layer*4*N+0:layer*4*N+4*N, layer*4*N+0:layer*4*N+4*N] = h00_Haldane
    for layer in range(layer_num-1):
        for i0 in range(N):
            # AB stacking
            h00[layer*4*N+i0*4+0, (layer+1)*4*N+i0*4+1] = tc
            h00[(layer+1)*4*N+i0*4+1, layer*4*N+i0*4+0] = tc
            h00[layer*4*N+i0*4+2, (layer+1)*4*N+i0*4+3] = tc
            h00[(layer+1)*4*N+i0*4+3, layer*4*N+i0*4+2] = tc

            # # BA stacking
            # h00[layer*4*N+i0*4+1, (layer+1)*4*N+i0*4+0] = tc
            # h00[(layer+1)*4*N+i0*4+0, layer*4*N+i0*4+1] = tc
            # h00[layer*4*N+i0*4+3, (layer+1)*4*N+i0*4+2] = tc
            # h00[(layer+1)*4*N+i0*4+2, layer*4*N+i0*4+3] = tc

            # # AA stacking
            # h00[layer*4*N+i0*4+0, (layer+1)*4*N+i0*4+0] = tc
            # h00[(layer+1)*4*N+i0*4+0, layer*4*N+i0*4+0] = tc
            # h00[layer*4*N+i0*4+1, (layer+1)*4*N+i0*4+1] = tc
            # h00[(layer+1)*4*N+i0*4+1, layer*4*N+i0*4+1] = tc
            # h00[layer*4*N+i0*4+2, (layer+1)*4*N+i0*4+2] = tc
            # h00[(layer+1)*4*N+i0*4+2, layer*4*N+i0*4+2] = tc
            # h00[layer*4*N+i0*4+3, (layer+1)*4*N+i0*4+3] = tc
            # h00[(layer+1)*4*N+i0*4+3, layer*4*N+i0*4+3] = tc
    return h00

def h01_of_bilayer_model(N, M, t1, t2, phi, tc, layer_num):
    h01_modified = h01_of_modified_Haldane_model(N, M, t1, t2, phi, sign=1)
    h01_Haldane = h01_of_Haldane_model(N, M, t1, t2, phi, sign=1)
    h01  = np.zeros((4*N*layer_num, 4*N*layer_num), dtype=complex)
    for layer in range(layer_num):
        if np.mod(layer,2) == 0:
            h01[layer*4*N+0:layer*4*N+4*N, layer*4*N+0:layer*4*N+4*N] = h01_modified
        if np.mod(layer,2) == 1:
            h01[layer*4*N+0:layer*4*N+4*N, layer*4*N+0:layer*4*N+4*N] = h01_Haldane
    return h01

def bands_of_bilayer_model(N, M, t1, t2, phi, tc, layer_num):
    k_array = np.linspace(0, 2*pi, 300)
    hamiltonian = functools.partial(hamiltonian_of_bilayer_model, N=N, M=M, t1=t1, t2=t2, phi=phi, tc=tc, layer_num=layer_num)
    eigenvalue_array = guan.calculate_eigenvalue_with_one_parameter(k_array, hamiltonian, print_show=1)
    plt, fig, ax = guan.import_plt_and_start_fig_ax(labelsize=25)
    guan.plot_without_starting_fig_ax(plt, fig, ax, k_array, eigenvalue_array, xlabel='\mathrm{k_x}', ylabel='\mathrm{E}', style='k', fontsize=30,  y_max=0.2, y_min=-0.2,linewidth=None, markersize=None, color=None, fontfamily='Times New Roman')
    plt.show()

def conductance_of_bilayer_model_with_disorder_array(N, M, t1, t2, phi, tc, layer_num):
    h00 = h00_of_bilayer_model(N, M, t1, t2, phi, tc, layer_num)
    h01 = h01_of_bilayer_model(N, M, t1, t2, phi, tc, layer_num)
    fermi_energy = 0
    disorder_intensity_array = np.arange(0, 2.5, .05)
    conductance_array = guan.calculate_conductance_with_disorder_intensity_array(fermi_energy, h00, h01, disorder_intensity_array, length=100, calculation_times=3, print_show=1) # length=2000, calculation_times=20
    guan.plot(disorder_intensity_array, conductance_array, xlabel='\mathrm{W_d}', ylabel='\mathrm{G(e^2/h)}', style='o-')

if __name__ == '__main__':
    main()

代码二

import numpy as np
from math import *
import cmath
import functools
import guan 

# Installation: pip install --upgrade guan

def main():
    M = 0
    t1 = 1
    phi= pi/2
    t2 = 0.03
    tc = 0.31
    bilayer_bands(M, t1, t2, phi, tc)
    chern_number(M, t1, t2, phi)

def hamiltonian_of_Haldane(k1, k2, M, t1, t2, phi, a=1/sqrt(3)): 
    h0 = np.zeros((2, 2), dtype=complex)
    h1 = np.zeros((2, 2), dtype=complex)
    h2 = np.zeros((2, 2), dtype=complex)

    h1[1, 0] = t1*(cmath.exp(1j*k2*a)+cmath.exp(1j*sqrt(3)/2*k1*a-1j/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j/2*k2*a))
    h1[0, 1] = h1[0, 1].conj()

    h2[0, 0] = t2*cmath.exp(1j*phi)*(cmath.exp(1j*sqrt(3)*k1*a)+cmath.exp(-1j*sqrt(3)/2*k1*a+1j*3/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j*3/2*k2*a))
    h2[1, 1] = t2*cmath.exp(-1j*phi)*(cmath.exp(1j*sqrt(3)*k1*a)+cmath.exp(-1j*sqrt(3)/2*k1*a+1j*3/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j*3/2*k2*a))
    matrix = h0 + h1 + h2 + h2.transpose().conj()
    return matrix

def hamiltonian_of_modified_Haldane(k1, k2, M, t1, t2, phi, a=1/sqrt(3)): 
    h0 = np.zeros((2, 2), dtype=complex)
    h1 = np.zeros((2, 2), dtype=complex)
    h2 = np.zeros((2, 2), dtype=complex)

    k1 = -k1  # kx  # Note that to get the unit cell the bilayer honeycomb lattice containing all possible layer hoppings, here one layer (HM layer or mHM layer） needs to be applied with its counterpart of mirror symmetry along x direction.

    h1[1, 0] = t1*(cmath.exp(1j*k2*a)+cmath.exp(1j*sqrt(3)/2*k1*a-1j/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j/2*k2*a))
    h1[0, 1] = h1[1, 0].conj()

    h2[0, 0] = t2*cmath.exp(1j*phi)*(cmath.exp(1j*sqrt(3)*k1*a)+cmath.exp(-1j*sqrt(3)/2*k1*a+1j*3/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j*3/2*k2*a))
    h2[1, 1] = t2*cmath.exp(1j*phi)*(cmath.exp(1j*sqrt(3)*k1*a)+cmath.exp(-1j*sqrt(3)/2*k1*a+1j*3/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j*3/2*k2*a))

    matrix = h0 + h1 + h2 + h2.transpose().conj()
    return matrix

def hamiltonian_of_bilayer(k1, k2, M, t1, t2, phi, tc):
    hamiltonian1 = hamiltonian_of_Haldane(k1, k2, M, t1, t2, phi, a=1/sqrt(3))
    hamiltonian2 = hamiltonian_of_modified_Haldane(k1, k2, M, t1, t2, phi, a=1/sqrt(3))
    hamiltonian = np.zeros((4, 4), dtype=complex)
    hamiltonian[0:2, 0:2] = hamiltonian1
    hamiltonian[2:4, 2:4] = hamiltonian2

    # AB stacking
    hamiltonian[1, 2] = tc   # B on HM， A on mHM
    hamiltonian[2, 1] = tc

    # BA stacking
    # hamiltonian[0, 3] = tc  # A on HM， B on mHM
    # hamiltonian[3, 0] = tc

    # AA stacking
    # hamiltonian[0, 2] = tc
    # hamiltonian[2, 0] = tc
    # hamiltonian[1, 3] = tc
    # hamiltonian[3, 1] = tc
    return hamiltonian

def bilayer_bands(M, t1, t2, phi, tc):
    k1_array = np.linspace(0, 4*pi, 3000)
    k2 = 0
    hamiltonian = functools.partial(hamiltonian_of_bilayer, k2=k2, M=M, t1=t1, t2=t2, phi=phi, tc=tc)
    eigenvalue_array = guan.calculate_eigenvalue_with_one_parameter(k1_array, hamiltonian)
    plt, fig, ax = guan.import_plt_and_start_fig_ax(labelsize=25)
    guan.plot_without_starting_fig_ax(plt, fig, ax, k1_array, eigenvalue_array[:, 0],  linewidth=2)
    guan.plot_without_starting_fig_ax(plt, fig, ax, k1_array, eigenvalue_array[:, 1],  linewidth=2)
    guan.plot_without_starting_fig_ax(plt, fig, ax, k1_array, eigenvalue_array[:, 2],  linewidth=2, color='k')
    guan.plot_without_starting_fig_ax(plt, fig, ax, k1_array, eigenvalue_array[:, 3],  linewidth=2, color='k')
    guan.plot_without_starting_fig_ax(plt, fig, ax, [], [], fontsize=30, y_max=1.5, y_min=-1.5, xlabel='\mathrm{k_x}', ylabel='\mathrm{E}')
    plt.show()

def chern_number(M, t1, t2, phi):
    tc_array = np.arange(0.05, 0.8, 0.05)
    chern_number_array = []
    for tc in tc_array:
        print(tc)
        hamiltonian = functools.partial(hamiltonian_of_bilayer, M=M, t1=t1, t2=t2, phi=phi, tc=tc)
        chern_number = guan.calculate_chern_number_for_honeycomb_lattice(hamiltonian, a=1/sqrt(3), precision=500, print_show=0)
        print(chern_number, '\n')
        chern_number_array.append(chern_number[0:2])
    guan.plot(tc_array, np.real(chern_number_array), xlabel='\mathrm{t_c}', ylabel='\mathrm{C}', style='o-')

if __name__ == '__main__':
    main()

1,829 次浏览

【说明：本站主要是个人的一些笔记和代码分享，内容可能会不定期修改。为了使全网显示的始终是最新版本，这里的文章未经同意请勿转载。引用请注明出处：https://www.guanjihuan.com】

头像

Published by guanjihuan

14 thoughts on “高陈数的蜂窝双层系统中的拓扑隐藏相变 Topological hidden phase transition in honeycomb bilayers with a high Chern number”

  1. 关老师您好！感谢您开源这份代码。我最近拜读了您的这篇论文，您开源的代码对我复现和学习这篇文章很有帮助。我在阅读您的代码时，发现了一个细节，想在这里留言和你探讨一下。
    我注意到，在代码一中，您在构建mHM哈密顿量的代码时，直接复用了HM的代码，只是翻转了一类子格点间次近邻跃迁的相位，这意味着您的mHM的原胞和HM的原胞是同样的结构；同时，您在构建双层系统的哈密顿量中，将全部层间耦合包含在了胞内耦合中，我认为这可能是矛盾的。如若mHM和HM的原胞是同一结构，那么mHM的原子3和同一原胞内的HM的原子2存在耦合，但是mHM的原子1会和相邻原胞的原子0存在耦合，这部分应该包含在h01内；如若要在一个原胞内包含所有层间耦合，那么我们不得不翻转两层中的一个原胞，此时两层的哈密顿量将不是简单的改变一类子格点间的跃迁相位。这可能是论文中准一维条带的能带和代码二绘制的ky=0能带存在不一致的原因。
    为验证这个猜想，举个最明显的例子，您可以用我后文附的代码替换掉您原代码中的h00_of_modified_Haldane_model和h01_of_modified_Haldane_model. 设置参数与论文中的Fig6(f)(g)(h)相同（t2=0.07, tc=0.73），您应该能看到准一维条带的能带和Fig6(f)能很好的对应，在N较大时，也能清晰地看到两条价带的带隙闭合时相交于一点。
    冒昧提出这个发现，希望您能抽出时间验证。再次感谢您的优秀工作和无私开源。附修改的代码：

    def h00_of_modified_Haldane_model(N, M, t1, t2, phi, sign):
    h00 = np.zeros((4*N, 4*N), dtype=complex)
    for i0 in range(N):
        h00[i0*4+0, i0*4+0] = M
        h00[i0*4+1, i0*4+1] = -M
        h00[i0*4+2, i0*4+2] = M
        h00[i0*4+3, i0*4+3] = -M
        h00[i0*4+0, i0*4+1] = t1
        h00[i0*4+1, i0*4+0] = t1
        h00[i0*4+1, i0*4+2] = t1
        h00[i0*4+2, i0*4+1] = t1
        h00[i0*4+2, i0*4+3] = t1
        h00[i0*4+3, i0*4+2] = t1
        # h00[i0*4+0, i0*4+2] = t2*cmath.exp(-1j*phi*sign)
        h00[i0*4+0, i0*4+2] = t2*cmath.exp(1j*phi*sign)
        h00[i0*4+2, i0*4+0] = h00[i0*4+0, i0*4+2].conj()
        # h00[i0*4+1, i0*4+3] = t2*cmath.exp(1j*phi*sign)
        h00[i0*4+1, i0*4+3] = t2*cmath.exp(-1j*phi*sign)
        h00[i0*4+3, i0*4+1] = h00[i0*4+1, i0*4+3].conj()
    for i0 in range(N-1):
        h00[i0*4+3, (i0+1)*4+0] = t1
        h00[(i0+1)*4+0, i0*4+3] = t1
        # h00[i0*4+2, (i0+1)*4+0] = t2*cmath.exp(1j*phi*sign)
        h00[i0*4+2, (i0+1)*4+0] = t2*cmath.exp(-1j*phi*sign)
        h00[(i0+1)*4+0, i0*4+2] = h00[i0*4+2, (i0+1)*4+0].conj()
        # h00[i0*4+3, (i0+1)*4+1] = t2*cmath.exp(-1j*phi*sign)
        h00[i0*4+3, (i0+1)*4+1] = t2*cmath.exp(1j*phi*sign)
        h00[(i0+1)*4+1, i0*4+3] = h00[i0*4+3, (i0+1)*4+1].conj()
    return h00

def h01_of_modified_Haldane_model(N, M, t1, t2, phi, sign):
    h01 = np.zeros((4*N, 4*N), dtype=complex)
    for i0 in range(N):
        # h01[i0*4+1, i0*4+0] = t1
        # h01[i0*4+2, i0*4+3] = t1
        h01[i0*4+0, i0*4+1] = t1
        h01[i0*4+3, i0*4+2] = t1
        h01[i0*4+0, i0*4+0] = t2*cmath.exp(1j*phi*sign)
        h01[i0*4+1, i0*4+1] = t2*cmath.exp(1j*phi*sign)
        h01[i0*4+2, i0*4+2] = t2*cmath.exp(1j*phi*sign)
        h01[i0*4+3, i0*4+3] = t2*cmath.exp(1j*phi*sign)
        # h01[i0*4+1, i0*4+3] = t2*cmath.exp(-1j*phi*sign)
        # h01[i0*4+2, i0*4+0] = t2*cmath.exp(-1j*phi*sign)
        h01[i0*4+3, i0*4+1] = t2*cmath.exp(-1j*phi*sign)
        h01[i0*4+0, i0*4+2] = t2*cmath.exp(-1j*phi*sign)
        if i0 != 0:
            # h01[i0*4+1, (i0-1)*4+3] = t2*cmath.exp(-1j*phi*sign)
            h01[i0*4+0, (i0-1)*4+2] = t2*cmath.exp(-1j*phi*sign)
    for i0 in range(N-1):
        # h01[i0*4+2, (i0+1)*4+0] = t2*cmath.exp(-1j*phi*sign)
        h01[i0*4+3, (i0+1)*4+1] = t2*cmath.exp(-1j*phi*sign)
    return h01

    回复

    1. 感谢您的阅读，观察很细致。我在论文中的附录 A 中也提到了“Note that to get the unit cell of the bilayer honeycomb lattice containing all possible layer hoppings [Fig.7(a)], here one layer (HM layer or mHM layer) needs to be applied with its counterpart of mirror symmetry along x direction.”。在双层的倒空间哈密顿量的构建上，我使用的是相近的哈密顿量，在相位变化的基础上，另外在不同层上做了个x方向的镜像翻转。

      在双层的实空间哈密顿量代码实现方面，我提供的代码 mHM 层和 HM 层的原子编号顺序可能会不大一样，一层是顺时针的1, 2, 3, 4 编号，一层是逆时针的 1, 2, 3, 4 编号，使得双层的耦合能在同一个元胞之内，所以能直接给出正确的结果。哈密顿量如果有不同形式，这应该跟原子编号的顺序有关，只要逻辑整体是自洽的，结果就能对上，在代码上主要体现为几个赋值位置多个负号。

      如果还有疑问或者有其他问题，可以随时留言或邮箱/微信联系我。

      回复

      1. 我能理解您说的一层顺时针编号，另一层是逆时针编号的含义，这和我在留言中说的将一个原胞翻转是相同的。但是这样对哈密顿量的影响不只是简单的几个相位添加负号，例如您在HM和mHM的h01中都写道：

        h01[i0*4+1, i0*4+0] = t1
        h01[i0*4+2, i0*4+3] = t1

        在顺时针编号的情况下，一列的0原子和左边的1原子连接，2原子和右边的3原子连接。但在逆时针编号的情况下，一列的0原子和右边的1原子连接，2原子和左边的3原子连接。两层的连接正好相反，h01也应该是不同的，我在代码中修改了这些内容，计算出的结果与您是有显著不同的。这一结果和我在保持均为顺时针编号，考虑h01内也存在耦合的情况是一直的。

        希望您能拨冗运行一下我修改的代码，例如在t_2=0.07, t_c=0.73, N=30的例子，与您Fig 6 (f)(g)(h)应该是相同的参数，您将看到以下结果：
        倒空间哈密顿量方法计算的能隙\Delta_{cv}和实空间哈密顿量计算的能隙\Delta_{cv_ribbon}是一致的。
        实空间哈密顿量计算的能带中也能看到两条价带间带隙的闭合。

        回复

        1. 嗯，我运行了您修改后的代码，以及认真思考了您的观察和分析，目前我感觉您的思路是对的。如果按这个纠正过后的双层模型，文章中涉及到实空间哈密顿量的结论（电导部分）就存在问题，您可以直接忽略；而文章的主题“拓扑隐藏相变现象”是对的，因为这是由倒空间的哈密顿量计算出的结论。另外，半反手性边缘态好像还是有的，这个物理图像是对的，只是电导不一定能作为一种探测方式，可能一开始就进入了强拓扑，因为价带总陈数始终是 2。

          我会更新博文，在置顶的位置特别强调您说的这个问题，感谢指出！很抱歉给您带来误导，您可以用自己搭建的这个模型开展相关的工作。

          如果还有其他的问题，可以随时保持联系，谢谢。

          回复

          1. 感谢您的亲自验证和回复。按照我的计算结果，论文中Fig3的电导变化趋势是基本正确的，只是电导变为2，即打开全局能隙的临界tc会减小，大致会在 3\sqrt{3}t_2 。
            祝您科研生活一切顺利。

            回复

        2. 非常感谢你发现的这个问题。另外，你可以考虑基于这个改正后的代码，研究多层的半反手性边缘态的演变和调控规律，或者其他的一些现象。我目前没有继续往这方面做，你有感兴趣的话可以深入研究。

          回复

  2. 关老师，请问这个双层AA堆叠的霍尔丹模型，怎么判断它是不是二阶拓扑绝缘体
    我想的可以通过构建一个有限大小的双层AA堆叠六边形霍尔丹模型，然后算他的能级图，和态密度分布，看有没有角态，从而判断它是不是二阶拓扑绝缘体吗？
    感觉这种是从现象来判断的，有没有从原理判断的

    回复

      1. 好的，感谢关老师的解答！我还有一个疑问，我在仿真有限大小的双层霍尔丹模型的时候，感觉不太会设置哈密顿量，不知道关老师有没有相关代码，我想参考一下,单层双层有限大小的霍尔丹模型都行

        回复

  3. 关老师，我想仿真一下有限大小的单层和双层霍尔丹模型的能带图，可以仿真吗，如果可以，应该以两个格点为基础，还是以四个格点为基础进行设置呀，我想象不出有限大小的蜂窝晶格的形状，不知道怎么设置哈密顿量

    回复

    1. 如果是有限大小的单层和双层，那么是没法画能带图的，只能画能级或态密度分布等。这是因为如果为有限大小，那么是没法傅里叶变换到 k 空间，不存在动量 k。另外有限大小的蜂窝晶格也有边界形状的选择问题，例如有 zigzag 边，armchair 边等。

      如果要算能带图，可以考虑二维的情况。选取最小的单元（元胞），单层应该是两个格点，双层的应该是四个格点。如果有其他掺杂，那么需要扩胞，元胞会比较大，在选取最小单元后写出该元胞的哈密顿量以及相邻元胞之间的跃迁矩阵。

      回复

发表评论

您的邮箱地址不会被公开。 必填项已用 * 标注

Captcha Code